Overview

In today’s lab we will work with water sample data two different aqueous chemsitry models. First, we will model mineral stability in R using the CHNOSZ package. Eh-pH diagrams show which species are thermodynamically favored at a range of pH and Eh conditions. These diagrams assume a fixed temperature, pressure, and treat dissolved species as ideal solutions where activitiy is equal to concentration. These diagrams are great tools for determining the thermodynamically stable phase under specific conditions.Once you have become familar with the functionality of the program we will then compare these saturation calculations to the the Gibbs free energy calculation.

Lab Submission

At the end of lab today you will upload your assignment as both a R Notebook file (.Rmd) and html file to our Canvas site. There is a template on the Canvas site to get you started with some tips on formatting you document. Your lab must be nicely organized and have a title and a header. - name your file using the format LabName_Lastname - for example for todays lab my file would be names “Lab5_Olson”

A very helpful guide on formatting R Markdown documents for assignments is available at http://www.stat.cmu.edu/~cshalizi/rmarkdown/.

Gettting Started

Go to https://docs.posit.co/cloud/get_started/. Log into your account and create a new R Studio project.

Part 1 Working with CHNOSZ

First we will install and load the package.

install.packages("CHNOSZ", repos = "https://cloud.r-project.org")
## 
## The downloaded binary packages are in
##  /var/folders/q8/hm5d074s0y91sv06nfv6qd1w0000gn/T//Rtmp6MYfbA/downloaded_packages
library(CHNOSZ)

Lets start with a simple Eh-pH plot

# Use CHNOSe basis (includes e-)
basis("CHNOSe")
##     C H N O S  Z ispecies logact state
## CO2 1 0 0 2 0  0      430     -3    aq
## H2O 0 2 0 1 0  0        1      0   liq
## NH3 0 3 1 0 0  0       64     -4    aq
## H2S 0 2 0 0 1  0       65     -7    aq
## e-  0 0 0 0 0 -1        2     -7    aq
## H+  0 1 0 0 0  1        3     -7    aq
#Minimal species set (just to generate an Eh–pH plot context)
species("H2O")   # harmless placeholder; ensures a valid diagram context
##   CO2 H2O NH3 H2S e- H+ ispecies logact state  name
## 1   0   1   0   0  0  0        1      0   liq water
#Compute affinities and plot
a <- affinity(pH = c(0, 14), Eh = c(-0.8, 1.2))
diagram(a, main = "Eh–pH framework (25 °C, 1 bar)")

#Add water stability lines
water.lines(a)

### Question

  1. What does the region between the water lines represent?

Now lets look at some species in this space. You can look up the species in the thermodynamic database using the info() function. You can use chemical formula or name to do this. For example info (“CH4”) and info(“methane”) will both pull up the number of aqueous species. You can use the info() function recursively to find the thermodynamic parameters info(info(“CH4”)). Or access fuzzy search using partial names for example info(“ribose+”). The subcrt() function will show calculate the standard thermodynamic properties from 0-350°C. Lets look at ferrous iron.

#First find the species using fuzzy search 
info("Fe+")
## [1] NA
#Now calculate the thermodynamic properties. 
subcrt("Fe+3")
## $species
##     name formula state ispecies model
## 255 Fe+3    Fe+3    aq      255   HKF
## 
## $out
## $out$`Fe+3`
##         T          P       rho       logK          G          H         S
## 1    0.01   1.000000 0.9998289  4.6019402 -24066.225  -46807.88 -267.6304
## 2   25.00   1.000000 0.9970614  3.0199808 -17238.080  -49580.40 -277.3992
## 3   50.00   1.000000 0.9880295  1.6540877 -10233.227  -51227.31 -282.7118
## 4   75.00   1.000000 0.9748643  0.4661679  -3107.124  -52802.67 -287.4043
## 5  100.00   1.013220 0.9583926 -0.5797922   4141.958  -54708.96 -292.6844
## 6  125.00   2.320144 0.9390726 -1.5127476  11530.897  -57206.84 -299.1385
## 7  150.00   4.757169 0.9170577 -2.3568854  19093.375  -60484.45 -307.0828
## 8  175.00   8.918049 0.8923427 -3.1308773  26862.062  -64603.81 -316.4791
## 9  200.00  15.536499 0.8647434 -3.8493110  34868.377  -69884.02 -327.8382
## 10 225.00  25.478603 0.8338733 -4.5260162  43164.442  -77343.63 -343.0138
## 11 250.00  39.736493 0.7990719 -5.1752426  51833.058  -88049.19 -363.6557
## 12 275.00  59.431251 0.7592362 -5.8078188  60948.409 -101454.73 -388.1554
## 13 300.00  85.837843 0.7124075 -6.4205474  70451.501 -111930.89 -406.0833
## 14 325.00 120.457572 0.6545772 -6.9791770  79921.608 -100135.64 -385.2236
## 15 350.00 165.211289 0.5746875 -7.3805360  88050.210   71534.85 -107.3251
##             V          Cp
## 1   -36.62951  -169.53528
## 2   -37.79023   -76.71380
## 3   -39.24367   -60.81813
## 4   -41.14509   -67.65675
## 5   -43.69977   -86.56386
## 6   -47.15404  -114.71067
## 7   -51.81179  -148.44772
## 8   -57.93142  -184.34609
## 9   -66.14503  -244.82388
## 10  -79.22608  -348.19110
## 11 -103.63536  -479.96106
## 12 -149.93737  -535.17079
## 13 -230.86491  -143.79280
## 14 -352.77114  2495.46959
## 15 -410.69579 40229.06809

Question

  1. What is the logK term in the table above? What does a positive logK term indicate?
# Now lets see when dissolves into water forming aqueous Fe+3

subcrt(c("pyrite", "Fe+3"), c("cr", "aq"), c(-1, 1), T = 25)
## $reaction
##      coeff   name formula state ispecies   model
## 1976    -1 pyrite    FeS2    cr     1976 CGL_Ttr
## 255      1   Fe+3    Fe+3    aq      255     HKF
## 65       2    H2S     H2S    aq       65     HKF
## 2       -1     e-   (Z-1)    aq        2     HKF
## 3       -4     H+      H+    aq        3     HKF
## 
## $out
##    T P       rho     logK        G        H         S        V       Cp
## 1 25 1 0.9970614 -15.2663 87140.17 46643.23 -144.6267 8.172212 219.5268

Question

  1. Based on the results of the calculation above. Is the reaction likely to occur at the conditions set (i.e t = 25 and activities -1 and 1 )? Note reaction coefficients are negative for reactants and positive for products. Why or why not is the reaction favored?
# Now lets change the activities and see what happens

subcrt(c("pyrite", "Fe+3"), c("cr", "aq"), c(1, -1), T = 25)
## $reaction
##      coeff   name formula state ispecies   model
## 1976     1 pyrite    FeS2    cr     1976 CGL_Ttr
## 255     -1   Fe+3    Fe+3    aq      255     HKF
## 65      -2    H2S     H2S    aq       65     HKF
## 2        1     e-   (Z-1)    aq        2     HKF
## 3        4     H+      H+    aq        3     HKF
## 
## $out
##    T P       rho    logK         G         H        S         V        Cp
## 1 25 1 0.9970614 15.2663 -87140.17 -46643.23 144.6267 -8.172212 -219.5268

Questions

  1. Based on the results of the calculation above. Is the reaction likely to occur at the conditions set (i.e t = 25 and activities 1 and -1 )? Note reaction coefficients are negative for reactants and positive for products. Why or why not is the reaction favored?

  2. Now consider the dissolution of carbon dioxide “CO2” in water from the “gas” to “aq” phases using the same coefficient and temperatures from the example above in the code block below.

#your code here for question 5

Lets look at an Eh-pH diagram for carbonate.

thermo(reset = TRUE)
## $reset
## NULL
thermo(t = 25, p = 1)
## $t
## NULL
## 
## $p
## NULL
# Basis for C-H-O with charge (Eh uses e-)
basis(c("CO2", "H2O", "H+", "e-"))
##     C H O  Z ispecies logact state
## CO2 1 0 2  0      430      0    aq
## H2O 0 2 1  0        1      0   liq
## H+  0 1 0  1        3      0    aq
## e-  0 0 0 -1        2      0    aq
# Carbonate aqueous species
species(c("CO2", "HCO3-", "CO3-2"))
##   CO2 H2O H+ e- ispecies logact state  name
## 1   1   0  0  0      430     -3    aq   CO2
## 2   1   1 -1  0       13     -3    aq HCO3-
## 3   1   1 -2  0       14     -3    aq CO3-2
# Eh–pH grid
a <- affinity(pH = c(0, 14, 200), Eh = c(-0.8, 1.2, 200))

# Predominance diagram 
diagram(a,
        xlab = "pH", ylab = "Eh (V)",
        main = "Carbonate system (CO2 / HCO3- / CO3^2-)")

# Add water stability bounds
water.lines(a)

### Question

  1. Based on the plot above, is carbonate speciation in natural waters controlled primarily by pH or by Eh? Explain your reasoning.

Now lets look at the speciation of carboante in waters.

# Carbonate speciation vs pH
basis(c("CO2", "H2O", "H+", "e-", "Ca+2"))
##      C Ca H O  Z ispecies logact state
## CO2  1  0 0 2  0      430      0    aq
## H2O  0  0 2 1  0        1      0   liq
## H+   0  0 1 0  1        3      0    aq
## e-   0  0 0 0 -1        2      0    aq
## Ca+2 0  1 0 0  2       10      0    aq
species(c("CO2", "HCO3-", "CO3-2","CaCO3"))
##   CO2 H2O H+ e- Ca+2 ispecies logact state  name
## 1   1   0  0  0    0      430     -3    aq   CO2
## 2   1   1 -1  0    0       13     -3    aq HCO3-
## 3   1   1 -2  0    0       14     -3    aq CO3-2
## 4   1   1 -2  0    1      108     -3    aq CaCO3
# 25 degrees C
a <- affinity(pH = c(0, 14))
e <- equilibrate(a)
diagram(e, alpha = TRUE)
# 0 degrees C
a <- affinity(pH = c(4, 12), T = 0)
e <- equilibrate(a)
diagram(e, alpha = TRUE, add = TRUE, col = 2, names = NA)
# Add legend
legend <- as.expression(list(lT(25), lT(0)))
legend("left", legend = legend, lty = 1, col = c(1, 2))

Question

  1. In the plot above alpha is the abundance of species. Based on the plot above is calcite (CaCO3) precipitation in an aqueous solution more or less kinetically favored in cold climates?
# Set conditions
thermo(reset = TRUE)
## $reset
## [1] TRUE
thermo(t = 10, p = 1)  # 10 °C, 1 bar
## $t
## [1] 25
## 
## $p
## [1] 1
# Define dissolution reactions (open CO2 form) 
# Calcite: CaCO3 + CO2 + H2O = Ca2+ + 2 HCO3-
sp_calc  <- c("calcite", "CO2", "H2O", "Ca+2", "HCO3-")
cf_calc  <- c(   -1,      -1,    -1,     +1,     +2)

# Dolomite: CaMg(CO3)2 + 2 CO2 + 2 H2O = Ca2+ + Mg2+ + 4 HCO3-
sp_dol   <- c("dolomite", "CO2", "H2O", "Ca+2", "Mg+2", "HCO3-")
cf_dol   <- c(   -1,      -2,    -2,     +1,     +1,     +4)

#Get logK from CHNOSZ thermodynamic database 
logK_calc <- subcrt(sp_calc, coeff = cf_calc, T = 10)$out$logK
logK_dol  <- subcrt(sp_dol,  coeff = cf_dol,  T = 10)$out$logK

# Convert logK -> ΔG° (kJ/mol reaction)
R <- 8.314462618  # J/mol/K
TK <- 273.15 + 10
dG_calc_kJ <- -2.303 * R * TK * logK_calc / 1000
dG_dol_kJ  <- -2.303 * R * TK * logK_dol  / 1000

# plot data 
par(mfrow = c(1, 2), mar = c(5, 4.5, 3, 1))

barplot(c(logK_calc, logK_dol),
        names.arg = c("Calcite", "Dolomite"),
        ylab = expression(log[10]*K~"(dissolution, open "~CO[2]*")"),
        main = "Equilibrium constant at 10 °C")

barplot(c(dG_calc_kJ, dG_dol_kJ),
        names.arg = c("Calcite", "Dolomite"),
        ylab = expression(Delta*G^degree~"(kJ/mol reaction)"),
        main = expression("Standard driving force at 10"~degree*C))

mtext("More positive logK and more negative ΔG° = more thermodynamically favored dissolution", side = 1, line = 4, cex = 0.9)

Question

  1. Using the information from the code above, which mineral is more thermodynamically favored to dissolve ?
# Set thermodyanic conditions
thermo(reset = TRUE)
## $reset
## [1] TRUE
thermo(t = 10, p = 1)   # 10 °C, 1 bar
## $t
## [1] 10
## 
## $p
## [1] 1
# Define acid-neutralization reactions
# Calcite + acid
# CaCO3 + H+ = Ca2+ + HCO3-
sp_cal <- c("calcite", "H+", "Ca+2", "HCO3-")
cf_cal <- c(   -1,     -1,    +1,     +1)

# Dolomite + acid
# CaMg(CO3)2 + 2 H+ = Ca2+ + Mg2+ + 2 HCO3-
sp_dol <- c("dolomite", "H+", "Ca+2", "Mg+2", "HCO3-")
cf_dol <- c(   -1,      -2,    +1,     +1,     +2)


# Pull equilibrium constants from CHNOSZ

logK_cal <- subcrt(sp_cal, coeff = cf_cal, T = 10)$out$logK
logK_dol <- subcrt(sp_dol, coeff = cf_dol, T = 10)$out$logK

# Convert logK -> ΔG° (kJ/mol reaction)
R  <- 8.314462618
TK <- 273.15 + 10
dG_cal_kJ <- -2.303 * R * TK * logK_cal / 1000
dG_dol_kJ <- -2.303 * R * TK * logK_dol / 1000


# Plot

par(mfrow = c(1, 2), mar = c(6, 4.5, 3, 1))

# logK comparison
barplot(c(logK_cal, logK_dol),
        names.arg = c("Calcite\n+ sulfuric acid",
                      "Dolomite\n+ sulfuric acid"),
        las = 2,
        ylab = expression(log[10]*K),
        main = expression("Carbonate dissolution at 10"~degree*C))

# ΔG° comparison
barplot(c(dG_cal_kJ, dG_dol_kJ),
        names.arg = c("Calcite\n+ sulfuric acid",
                      "Dolomite\n+ sulfuric acid"),
        las = 2,
        ylab = expression(Delta*G^degree~"(kJ/mol reaction)"),
        main = expression("Thermodynamic driving force"))

mtext("More positive logK and more negative ΔG° = more favorable dissolution",
      side = 1, line = 4, cex = 0.9)

# Print values for discussion
cat("\n--- 10 °C, sulfuric-acid conditions ---\n")
## 
## --- 10 °C, sulfuric-acid conditions ---
cat("Calcite:  logK =", logK_cal, "  ΔG° =", dG_cal_kJ, "kJ/mol rxn\n")
## Calcite:  logK = 2.233834   ΔG° = -12.11143 kJ/mol rxn
cat("Dolomite: logK =", logK_dol, "  ΔG° =", dG_dol_kJ, "kJ/mol rxn\n")
## Dolomite: logK = 3.774217   ΔG° = -20.4631 kJ/mol rxn

Question

  1. In the above example we added pyrite oxidation to the chemical reactions. Given the results what can you say about the results? How do they compare the prior scenario? What are the environmental implications for this?

Graduate level Question (undergraduate extra credit)

  1. Below create an Eh-pH diagram for another element of interest. Provide one environmental scenario where these species play a role.
# Your code here for questiono 10