demo/DEW.R

In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry

# CHNOSZ/demo/DEW.R
# Demo for the Deep Earth Water (DEW) model in CHNOSZ 20170927
library(CHNOSZ)

# Set up subplots
opar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2), mar = c(3.0, 3.5, 2.5, 1.0), mgp = c(1.7, 0.3, 0), las = 1, tcl = 0.3, xaxs = "i", yaxs = "i")

# Activate DEW model
oldwat <- water("DEW")

###########
#### Plot 1: quartz solubility at high pressure
## after Figure 7D of Sverjensky et al., 2014a [SHA14]
## (Geochim. Cosmochim. Acta, https://doi.org/10.1016/j.gca.2013.12.019)
###########

# Load SiO2 and Si2O4 data taken from DEW spreadsheet
iSi <- add.OBIGT("DEW", c("SiO2", "Si2O4"))
# Print the data references to confirm we got the right ones
thermo.refs(iSi)
# Set temperature ranges for different pressures
# data.frame is used to make P and T the same length
PT0.5 <- data.frame(P = 500, T = seq(200, 550, 10))
PT1.0 <- data.frame(P = 1000, T = seq(200, 700, 10))
PT2.0 <- data.frame(P = 2000, T = seq(200, 700, 10))
PT5.0 <- data.frame(P = 5000, T = seq(200, 850, 10))
PT10.0 <- data.frame(P = 10000, T = seq(200, 825, 10))
PT20.0 <- data.frame(P = 20000, T = seq(200, 800, 10))
PT <- rbind(PT0.5, PT1.0, PT2.0, PT5.0, PT10.0, PT20.0)
# Reaction 1: quartz = SiO2(aq) [equivalent to quartz + 3 H2O = Si(OH)4]
SiO2_logK <- suppressWarnings(subcrt(c("quartz", "SiO2"), c("cr", "aq"), c(-1, 1), P = PT$P, T = PT$T))$out$logK
# Reaction 2: 2 quartz = Si2O4(aq) [equivalent to 2 quartz + 3 H2O = Si2O(OH)6]
Si2O4_logK <- suppressWarnings(subcrt(c("quartz", "Si2O4"), c("cr", "aq"), c(-2, 1), P = PT$P, T = PT$T))$out$logK
# Plot the sum of molalities (== activities) for each pressure
plot(c(200, 1000), c(-2.5, 0.5), type = "n", xlab = axis.label("T"), ylab = "log molality")
for(P in unique(PT$P)) {
  icond <- PT$P == P
  SiO2_logm <- SiO2_logK[icond]
  Si2O4_logm <- Si2O4_logK[icond]
  logm <- log10(10^SiO2_logm + 10^Si2O4_logm)
  lines(PT$T[icond], logm)
  # add text label
  lastT <- tail(PT$T[icond], 1)
  Pkb <- paste(format(P/1000, nsmall = 1), "kb")
  text(lastT+25, tail(logm, 1), Pkb, adj = 0)
}
t1 <- quote("Solubility of"~alpha*"-quartz")
t2 <- "after Sverjensky et al., 2014a"
mtitle(as.expression(c(t1, t2)))
# TODO: lines are a little low at highest P and T ...

###########
#### Plot 2: correlations between non-solvation volume and HKF a1 parameter
## after Figures 12B and 12C of Sverjensky et al., 2014a [SHA14]
###########

# We can't use the DEW water model at 25 degC
reset()
# Load the fitted parameters for species as used by SHA14
# TODO: also use their Ca+2??
# NOTE: don't load NaCl, NH4+, or HS- here because the DEW spreadsheet lists a1 from the correlation
iDEW <- add.OBIGT("DEW", c("CO3-2", "BO2-", "MgCl+", "SiO2", "HCO3-", "Si2O4"))
# Override the check in subcrt() for the DEW water model 20220929
lapply(iDEW, mod.OBIGT, model = "HKF")
# Set up the plot
V0nlab <- expression(Delta * italic(V) * degree[n]~~(cm^3~mol^-1))
a1lab <- expression(italic(a)[1]%*%10~~(cal~mol~bar^-1))
plot(c(-25, 50), c(-4, 12), type = "n", xlab = V0nlab, ylab = a1lab)
# A function to get the HKF parameters, calculate nonsolvation volume, plot points, labels, error bars, and correlation lines
plotfun <- function(species, col, pch, cex, dy, error, xlim, corrfun) {
  # Get HKF a1 parameter
  a1 <- info(info(species), check.it = FALSE)$a1 * 10
  # Set omega to zero to calculate non-solvation volume 20220326
  mod.OBIGT(species, omega = 0)
  Vn <- do.call(rbind, subcrt(species, T = 25)$out)$V
  points(Vn, a1, col = col, pch = pch, cex = cex)
  for(i in 1:length(species)) text(Vn[i], a1[i] + dy, expr.species(species[i]))
  arrows(Vn, a1 - error, Vn, a1 + error, length = 0.03, angle = 90, code = 3, col = col)
  lines(xlim, corrfun(xlim), col = col)
}
# Monovalent ions: Na+, K+, Cl-, Br-
monofun <- function(Vn) 2.0754 + 0.10871 * Vn
# For easier reading, set y-offset to NA so the labels aren't plotted
plotfun(c("Na+", "K+", "Cl-", "Br-"), "red", 19, 1, NA, 0.5, c(-7, 35), monofun)
# Divalent ions: Mg+2, Ca+2, CO3-2, SO4-2
difun <- function(Vn) 3.5321 + 0.23911 * Vn
plotfun(c("Mg+2", "Ca+2", "CO3-2", "SO4-2"), "black", 15, 1, 1.2, 0.7, c(-20, 25), difun)
# Complexes and neutral molecules: BO2-, MgCl+, SiO2, NaCl, HCO3-, Si2O4, NH4+, HS-
compfun <- function(Vn) 1.5204 + 0.19421 * Vn
plotfun(c("MgCl+", "SiO2", "NaCl", "HCO3-", "Si2O4"), "blue1", 18, 1.5, 1, 0.5, c(-20, 50), compfun)
# For easier reading, put some labels below the points
plotfun(c("BO2-", "NH4+", "HS-"), "blue1", 18, 1.5, -1.2, 0.5, c(-20, 50), compfun)
# Include an empty subscript for better spacing between the lines
t1 <- quote("Correlations between non-solvation"[])
t2 <- quote("volume and HKF "*italic(a)[1]*" parameter")
mtitle(as.expression(c(t1, t2)))

# Clean up for next plot 
OBIGT()
water("DEW")

###########
#### Plot 3: logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure
## after Figure 1b of Sverjensky et al., 2014b [SSH14]
## (Nature Geoscience, https://doi.org/10.1038/NGEO2291)
###########

# Define system
basis("CHNOS+")
inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4")
organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
species(c(inorganic, organic), 0)
fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic)))

# A function to make the diagrams
dfun <- function(T = 600, P = 50000, res = 300) {
  a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P)
  # Set total C concentration to 0.03 molal
  # (from EQ3NR model for eclogite [Supporting Information of SSH14])
  e <- equilibrate(a, loga.balance = log10(0.03))
  diagram(e, fill = fill)
  dp <- describe.property(c("     T", "     P"), c(T, P), digits = 0)
  legend("bottomleft", legend = dp, bty = "n")
}

water("DEW")
## Not run: make diagram using CHNOSZ default database
## (not recommended for high P)
#dfun()
# Make diagram using CO2, HCO3-, CO3-2, CH4, and acetic acid data from DEW spreadsheet
# (the acetate field disappears if we also use the data for acetate from the spreadsheet 20200629)
#add.OBIGT("DEW", c("CO2", "HCO3-", "CO3-2", "CH4", "acetic acid"))
add.OBIGT("DEW")
dfun()
mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal"))

###########
#### Plot 4: speciation of carbon as a function T, logfO2 and pH (added 20171008)
## after SSH14 Fig. 3
###########

# Conditions:
# T = 600, 700, 800, 900, 1000 degC
# P = 5.0GPa (50000 bar)
# fO2 = QFM - 2
# pH set by jadeite + kyanite + coesite
# output from EQ3NR calculations (SSH14 Supporting Information)
# dissolved carbon: 0.03, 0.2, 1, 4, 20 molal
# true ionic strength: 0.39, 0.57, 0.88, 1.45, 2.49
# pH: 3.80, 3.99, 4.14, 4.25, 4.33
## Activate DEW model
reset()
water("DEW")
# Add species data for DEW
inorganics <- c("CH4", "CO2", "HCO3-", "CO3-2")
organics <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
add.OBIGT("DEW")
## Set basis species
basis(c("Fe", "SiO2", "CO3-2", "H2O", "oxygen", "H+"))
## Calculate logfO2 in QFM buffer
basis("O2", "QFM")
T <- seq(600, 1000, 100)
buf <- affinity(T = T, P = 50000, return.buffer = TRUE)
## Add species
species(c(inorganics, organics))
## Generate spline functions from IS, pH, and molC values at every 100 degC
IS <- c(0.39, 0.57, 0.88, 1.45, 2.49)
pH <- c(3.80, 3.99, 4.14, 4.25, 4.33)
molC <- c(0.03, 0.2, 1, 4, 20)
## Use extended Debye-Huckel equation with b_gamma set to zero
nonideal("bgamma0")
## Calculate affinities on the T-logfO2-pH-IS transect
a <- affinity(T = T, O2 = buf$O2 - 2, IS = IS, pH = pH, P = 50000)
## Calculate metastable equilibrium activities using the total
## Carbon molality as an approximation of total activity
e <- equilibrate(a, loga.balance = log10(molC))
## Make the diagram; don't plot names of low-abundance species
names <- c(inorganics, organics)
names[c(4, 5, 7, 9)] <- ""
col <- rep("black", length(names))
col[c(1, 3, 6, 8, 10)] <- c("red", "darkgreen", "purple", "orange", "navyblue")
if(packageVersion("CHNOSZ") > "1.1.3") {
  diagram(e, alpha = "balance", names = names, col = col, ylim = c(0, 0.8), ylab = "carbon fraction", spline.method = "natural")
} else {
  diagram(e, alpha = "balance", names = names, col = col, ylim = c(0, 0.8), ylab = "carbon fraction")
}

## Add legend and title
ltxt1 <- "P = 50000 bar"
ltxt2 <- substitute(logfO2=="QFM-2", list(logfO2 = axis.label("O2")))
pH <- seq(3.8, 4.3, length.out = length(T))
legend("left", legend = as.expression(c(ltxt1, ltxt2)), bty = "n")
t1 <- "Aqueous carbon speciation"
t2 <- "after Sverjensky et al., 2014b"
mtitle(c(t1, t2))

### Additional checks

# The maximum absolute pairwise difference between x and y
maxdiff <- function(x, y) max(abs(y - x))
## Check that we're within 0.1 of the QFM-2 values used by SSH14
stopifnot(maxdiff((buf$O2-2), c(-17.0, -14.5, -12.5, -10.8, -9.4)) < 0.1)

# Here are the logKs of aqueous species dissociation reactions at 600 degC and 50000 bar,
# values from EQ3NR output in Supporting Information of the paper (p. 103-109):
inorganic.logK <- c(24.4765, -9.0784, -5.3468, 0)
organic.logK <- c(1.7878, 2.5648, 15.3182, 16.9743, 30.4088, 28.9185)
# Calculate equilibrium constants of the reactions in CHNOSZ; use a negative sign to change from formation to dissociation
logK.calc <- -unlist(affinity(T = 600, P = 50000, property = "logK")$values)
logK.calc - c(inorganic.logK, organic.logK)
## Except for acetate, we're within 0.021 of the logK values used by SSH14
stopifnot(maxdiff(logK.calc[-8], c(inorganic.logK, organic.logK)[-8]) < 0.021)

## Check that we get similar activity coefficients
# Activity coefficients for monovalent species from EQ3NR output
loggamma <- c(-0.15, -0.18, -0.22, -0.26, -0.31)
# Activity coefficients calculated in CHNOSZ
sres <- subcrt("propanoate", T = seq(600, 1000, 100), P = 50000, IS = c(0.39, 0.57, 0.88, 1.45, 2.49))
stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.023)
# If m_star in nonideal() was zero, we could decrease the tolerance here
#stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.004)

# Reset OBIGT database
reset()

par(opar)