Choosing the water model (DEW)

This vignette shows the conversion of a GWB data file to use the Deep Earth Water (DEW) model.

This vignette was compiled on 2023-09-15 with logKcalc 0.1.1-2 and CHNOSZ 2.0.0-25.

The Deep Earth Water (DEW) model includes correlations for H₂O and revised thermodynamic parameters of aqueous species for estimating thermodynamic properties to pressures of up to 6 GPa (Sverjensky et al., 2014a). This model was used to predict a significant stability region for organic ions, such as acetate and propionate, at high pressures corresponding to upper mantle conditions (Sverjensky et al., 2014b).

Let us use the DEW model in CHNOSZ in order to calculate equilibrium constants for dissociation reactions in a GWB thermodynamic data file. Separate commands are needed to change the water model and load the parameters of the revised Helgeson-Kirkham-Flowers (HKF) equations for aqueous species that are taken from the DEW spreadsheet.

reset()
water("DEW")
if(packageVersion("CHNOSZ") >= "1.4.0") add.OBIGT("DEW") else add.obigt("DEW")
# modOBIGT(c("addSUPCRT", "steam"))

Note that other changes to the database to add the properties of steam and minerals from SUPCRT92 (Johnson et al., 1992) are not needed for this calculation.

The thermo_12elements.tdat data file includes aqueous CH₄, CO₂, HCO₃^-, CO₃^‑‑, CH₃COO^- (acetate) and HCH₃COO (acetic acid), which will be updated during the conversion. We can use the ispecies argument to add other species from OBIGT, in particular formic and propanoic acid and their anions, which were predicted to be major species by Sverjensky et al. (2014b). For completeness, we include other aqueous organic species, although they are not expected to appear in the diagrams. However, we do not include glycinate (added in the DEW 2019 dataset) to maintain backwards compatibility with CHNOSZ versions before 1.4.0.

major <- c("formic acid", "formate", "propanoic acid", "propanoate")
o1 <- c("benzene", "diglycine", "diketopiperazine")
o2 <- c("ethane", "ethanol", "ethylene")
o3 <- c("glutamate", "glutamic acid", "H-glutarate", "glutaric acid")
#o4 <- c("glutamine", "glycinate", "glycine", "glycolate", "glycolic acid")
o4 <- c("glutamine", "glycine", "glycolate", "glycolic acid")
o5 <- c("H-succinate", "hexane", "isobutane")
o6 <- c("lactate", "lactic acid", "leucine", "methanol")
o7 <- c("propane", "propanol", "toluene", "urea")
ispecies <- info(c(major, o1, o2, o3, o4, o5, o6, o7))

In the GWB thermodynamic data files, pressure is not an independent variable, but is expanded in terms of the principal temperatures (Bethke and Farrell, 2020). In order to set up calculations to compare different pressures, we choose a series of closely spaced temperatures that approximate a constant temperature of 600 °C. The pressures are set to 0.5, 1, 1.5, 2, 2.5, 3, 4, 5 GPa, converted to units of bar.

T <- 593:600
P <- c(5000, 10000, 15000, 20000, 25000, 30000, 40000, 50000)

Now we perform the conversion. To calculate the extended term parameter of the Debye-Hückel equation at high pressure (see Manning et al., 2013; Dick, 2019), we use the bgamma method. This does not affect the activity diagrams shown below, but is important for speciation calculations that may be possible with GWB.

infile <- system.file("extdata/thermo_12elements.tdat", package = "logKcalc")
outfile <- file.path(tempdir(), "logKcalc_vig4.tdat")
logKcalc(infile, outfile, T = T, P = P, ispecies = ispecies, DH.method = "bgamma", maxprint = 10)

## The temperatures in the input file are:

## [1] "0, 25, 60, 100, 150, 200, 250, 300"

## The pressures in the input file are:

## [1] "1.0134, 1.0134, 1.0134, 1.0134, 4.76, 15.549, 39.776, 85.927"

## The temperatures for the logK calculation are:

## [1] "593, 594, 595, 596, 597, 598, 599, 600"

## The pressures for the logK calculation are:

## [1] "5000, 10000, 15000, 20000, 25000, 30000, 40000, 50000"

## The 13 available basis species have 13 elements:

## [1] "Al, Au, C, Cl, Fe, H, K, N, Na, O, S, Si, Z"

## Can't map these redox species to the OBIGT database:

## [1] "(O-phth)--"

## The basis or redox species for these aqueous species are not available:

## [1] "Al(O-phth)+, Al(O-phth)2-, H(O-phth)-, H2(O-phth), Na(O-phth)-"

## Can't map 32 aqueous species to the OBIGT database.

## The basis or redox species for these mineral species are not available:

## [1] "O-phth acid(c)"

## Can't map 53 mineral species to the OBIGT database.

## Reaction among pyrrhotite,H+,Fe+2,HS-,Fe+3 was unbalanced, missing Fe0.125

## Removing species with unavailable logKs:

## [1] "Alunite, Goethite, Pyrrhotite, Sulfur-Rhmb"

## Can't map these gas species to the OBIGT database:

## [1] "Steam"

## Adding aqueous species formic acid with a0 = 0

## Adding aqueous species formate with a0 = 3.72

## Adding aqueous species propanoic acid with a0 = 0

## Adding aqueous species propanoate with a0 = 3.72

## Adding aqueous species benzene with a0 = 0

## info.numeric: H of diglycine(aq) is NA; set to -171414.91 cal mol-1

## Adding aqueous species diglycine with a0 = 0

## info.numeric: H of diglycine(aq) is NA; set to -171414.91 cal mol-1
## info.numeric: H of diglycine(aq) is NA; set to -171414.91 cal mol-1

## info.numeric: H of diketopiperazine(aq) is NA; set to -101562.28 cal mol-1

## Adding aqueous species diketopiperazine with a0 = 0

## info.numeric: H of diketopiperazine(aq) is NA; set to -101562.28 cal mol-1
## info.numeric: H of diketopiperazine(aq) is NA; set to -101562.28 cal mol-1

## Adding aqueous species ethane with a0 = 0

## Adding aqueous species ethanol with a0 = 0

## Adding aqueous species ethylene with a0 = 0

## info.numeric: H of glutamate(aq) is NA; set to -234035.26 cal mol-1

## Adding aqueous species glutamate with a0 = 3.72

## info.numeric: H of glutamate(aq) is NA; set to -234035.26 cal mol-1
## info.numeric: H of glutamate(aq) is NA; set to -234035.26 cal mol-1

## Adding aqueous species glutamic acid with a0 = 0

## Adding aqueous species H-glutarate with a0 = 3.72

## Adding aqueous species glutaric acid with a0 = 0

## Adding aqueous species glutamine with a0 = 0

## Adding aqueous species glycine with a0 = 0

## Adding aqueous species glycolate with a0 = 3.72

## Adding aqueous species glycolic acid with a0 = 0

## Adding aqueous species H-succinate with a0 = 3.72

## Adding aqueous species hexane with a0 = 0

## Adding aqueous species isobutane with a0 = 0

## Adding aqueous species lactate with a0 = 3.72

## Adding aqueous species lactic acid with a0 = 0

## Adding aqueous species leucine with a0 = 0

## Adding aqueous species methanol with a0 = 0

## Adding aqueous species propane with a0 = 0

## Adding aqueous species propanol with a0 = 0

## Adding aqueous species toluene with a0 = 0

## Adding aqueous species urea with a0 = 0

## Writing output to /tmp/RtmpeEM4Fc/logKcalc_vig4.tdat ...

## [1] "13 basis species"
## [1] "10 redox species"
## [1] "69 aqueous species"
## [1] "1 electron species"
## [1] "32 mineral species"
## [1] "7 gas species"
## [1] "8 oxide species"

Scripts and output of GWB runs.

This script for the Act2 program plots the speciation of carbon as a function of oxygen fugacity (f_O2) and pH. The activity of HCO₃^- is set to 0.03, which is close to the concentration of carbon in the fluid calculated for a model eclogite at 600 °C and 5 GPa (Supplementary Information of Sverjensky et al., 2014b). So as not to exceed the maximum pressure accepted by the Act2 program, the script sets a pressure of 5000 bar. This pressure affects the water limits and stability fields of gases but not the log K values for reactions between aqueous species, which were calculated for higher pressures as described above. To make these diagrams, we turn off the water limits and gas and mineral fields in Act2.

The first diagram shows a small stability field for formate at 0.5 GPa (5 kbar). This field is present both for the DEW model and for the water model from SUPCRT92 (Johnson et al., 1992) with the default OBIGT database in CHNOSZ, which has parameters of the organic acids from Shock (1995). (The GWB data file was created as shown above except the water and add.OBIGT commands were omitted and all pressures were set to 5000 bar, the maximum for the revised HKF equations without the DEW model.)

The next pair of diagrams is for a pressure of 5 GPa (50 kbar) using the DEW model. In the first, propanoate and formate are predicted to predominate at some redox and pH conditions. In the second, the formation of both of these species is suppressed; the diagram shows that no other organic species can predominate at equilibrium.

The presence of predominance fields for propanoate and formate is consistent with the model shown in Figure 3 of Sverjensky et al. (2014b), where these species are the most abundant C-bearing species in the fluid at 600 °C. It is unclear why these species are not present in the logf_O2–pH diagram in Figure 1 of that paper, which instead shows a large predominance field for acetate. However, it may not be surprising that we have trouble reproducing some results of Sverjensky et al. (2014b); according to a comment in the DEW spreadsheet the acetate parameters were revised after that publication, on January 26th, 2016.

Comparison of equilibrium constants

To check the accuracy of the data in OBIGT and equations in CHNOSZ, let’s compare the calculations with the DEW spreadsheet. The dissociation reactions of the species in the GWB file are written in terms of the basis species HCO₃^-, H₂O, O₂(aq), and H⁺. We can use logKcomp to retrieve the equilibrium constants calculated for the reactions at 600 °C and 5 GPa, which is the 8th T,P pair. For comparison, the 2019 version of the DEW spreadsheet was used to calculate log K for the same reactions and conditions. HCO₃^- is not included; because it is a basis species, its dissociation reaction has a log K of zero by definition.

species <- c("CH4", "CO2", "CO3--", "formate", "CH3COO-", "propanoate")
# logK values calculated with CHNOSZ
calc <- logKcomp(outfile, iTP = 8)
cvals <- calc[species]
# logK values calculated with the DEW spreadsheet
dvals <- c(48.4582, -3.7302, 5.3468, 12.5707, 48.6838, 77.5725)

logKs <- t(data.frame(dvals, cvals, dvals - cvals))
colnames(logKs) <- paste0("`", species, "`")
rownames(logKs) <- c("spreadsheet", "CHNOSZ", "difference")
knitr::kable(logKs)

ratios <- t(data.frame(dvals / cvals))
colnames(ratios) <- paste0("`", species, "`")
rownames(ratios) <- "ratio"
knitr::kable(ratios)

The differences between the CHNOSZ calculation and DEW spreadsheet are roughly proportional to the magnitude of log K. The ratios are close to 1.000666, the ratio between the values of the gas constant (R) used in CHNOSZ (roughly 1.9872 cal mol^-1 K^-1, converted from 8.314463 J mol^-1 K^-1) and in the DEW spreadsheet (1.9858775 cal mol^-1 K^-1). These small differences have a negligible effect on the diagrams shown above.