Customizing the thermodynamic database

Types of data

• Columns 1-8 have character data.
• Columns 9-21 have numeric data.
  • Columns 9-13 have standard-state thermodynamic properties at 25 °C and 1 bar.
  • Columns 14-20 have parameters for calculating thermodynamic properties at other temperatures and pressures. There are two main models:
    • Revised Helgeson-Kirkham-Flowers “equation-of-state” parameters for aqueous species (HKF)
    • Heat capacity coefficients for crystalline, gaseous, and liquid species (CGL).
    • The columns are named by combining the the names of the HKF and CGL coefficients, separated by a dot.
  • Column 21 has the charge used in the HKF EOS or the maximum temperature for CGL species. The HKF EOS parameter for charge (Z), which appears in the g function for the temperature derivatives of the ω parameter (Shock et al., 1992), is taken from this column and not from the chemical formula of the species.

Ranges of HKF and CGL models

To a first approximation, the revised HKF equations of state are applicable within the stability region of liquid water or the supercritical fluid with a density greater than 0.35 g/cm3, and not exceeding the ranges of 0 to 1000 °C and 1 to 5000 bar (see Shock et al., 1992 for details). There are two ways in which these limits are enforced in CHNOSZ:

• The default source of water properties (H2O92D Fortran subroutine modified from SUPCRT92) yields NA values beyonds these T and P ranges. Note that the Deep Earth Water (DEW) model is available to extend the applicable range to pressures of up to 60 kbar (6 GPa) (Sverjensky et al., 2014).
• subcrt() generates NA values beyond the density limit; the exceed.rhomin argument can be used to enable calculations at lower density in the supercritical region.

The first three terms in the CGL heat capacity equation correspond to the Maier-Kelley equation for heat capacity (Maier and Kelley, 1932); the additional terms are useful for representing heat capacities of minerals (Robie and Hemingway, 1995) and organic gases and liquids (Helgeson et al., 1998). The upper temperature limit for applicability of the CGL heat-capacity equation corresponds to the temperature of phase transition or to a temperature beyond which extrapolations from experimental data becomes highly uncertain. A maximum temperature can be stored in the final column of OBIGT, and is used by subcrt() to replace computed values with NA above the maximum.

Required and optional data

REQUIRED:

• All species need a name and a state. The state can be one of aq, gas, or cr. For minerals with higher-temperature polymorphs, they are named cr2, cr3, etc. NOTE: cr stands for “crystalline”; this naming convention (which was inherited from SUPCRT92 data files) refers to any solid phases including amorphous SiO₂ and other minerals.
• A chemical formula is required to do almost anything useful in CHNOSZ (e.g. check reaction balancing with subcrt() and add species with basis() or species()). NOTE: The convention has been adopted that the name of aqueous species in OBIGT is the same as the chemical formula. Most minerals, gases, and liquids have a name that is a common name.
• Likewise, E_units needs to be specified to perform any calculations of thermodynamic properties. The value can be J for Joules or cal for calories.

OPTIONAL: Everything else. Really, it depends on what you need. For instance, if you just want to use subcrt() to calculate log K of a reaction from ΔG° of species at 25 °C, then G is the only parameter that is needed.

OPTIONAL but useful:

• abbrv can be an abbreviation (e.g. Qtz for quartz). It is used by info() (together with name and formula) to look up species in the database.
• date is a timestamp for the data entry (YYYY-MM-DD format in the default OBIGT database).
• ref1 and ref2 are bibliographic reference keys. They have matching entries in extdata/OBIGT/refs.csv, which is used by thermo.refs() to display references, and in vignettes/OBIGT.bib, which is used in the OBIGT.Rmd vignette to produce a reference list.

NOTE: Other functions in CHNOSZ do not depend on date, ref1, and ref2, so you can put anything there that is convenient for you.

NA or 0?

If a character value or thermodynamic parameter (Columns 1-13) is unknown, use NA. Note that a missing (blank) value in the file is treated as NA.

• Unknown values for character values (Columns 1-8) should be NA.
• If you have only two of G, H, and S, then the missing one should be NA (definitely not 0). NOTE: info() and dependent functions (particularly subcrt()) “know” the equation DG = DGH - TDS; together with the entropies of the elements, this equation is used to compute a single missing one of G, H, or S from the other two.
• If you don’t have Cp or V, then set it to NA.
  • If HKF or CGL parameters are present: they will be used to calculate Cp, so thermodynamic properties can be calculated at T > 25 °C.
  • If HKF or CGL parameters aren’t present: thermodynamic properties can’t be calculated at T > 25 °C (NAs will propagate to higher T).

If an “equation-of-state” or heat capacity coefficient (Columns 14-20) is unknown, use 0.

• Furthermore, if you would like to assume that Cp or V is 0, then set it to 0.
  • Thermodynamic properties will be extrapolated to T > 25 °C and P > 1 bar assuming that Cp and V are 0.

More detail on the inner working of the functions: For both HKF and CGL, if at least one parameter for a species is provided, any NA values of the other parameters are taken to be zero. If all EOS parameters are NA, but values of Cp and/or V are present, they are assumed to be constants for extrapolating thermodynamic properties (e.g. ΔG°) as a function of temperature and pressure.

OOM scaling and info()

To make data entry easier, some values in the the CSV files and OBIGT data frame are scaled by order-of-magnitude (OOM) factors. This is particularly applicable to HKF parameters, for which OOM scaling is nearly always used in published data tables. OOM scaling is also applied to the CGL parameters. See ?thermo for details of the OOM scaling.

info() provides a simple user interface to the OBIGT database and is called by other functions in CHNOSZ to retrieve unscaled values from the database. This is a summary of its main features:

• Remove OOM scaling. This is used primarily by other functions in CHNOSZ to get a set of unscaled model parameters for calculating thermodynamic properties as a function of T and P.
• Extract the HKF or CGL parts of column names (only if all matching species have the same model).
• Calculate a missing one of G, H, or S if two of them are present.
• Cross-check G, H, and S if all of them are present, and print a message if the difference is above a threshold (see ?checkGHS).
• Calculate a missing Cp or V from the model parameters, if possible.
• Cross-check Cp or V (if present) against the model parameters, if possible, and print a message if the difference is above a threshold (see ?checkEOS).

NOTE: info() does NOT change the units; the values it displays (including possibly calculated ones) correspond to the E.units() for that species. On the other hand, subcrt() outputs values in the units selected with E.units().

add.OBIGT() with optional data files

TODO

add.OBIGT() with other CSV files

add.OBIGT() can also be used to add data from a user-specified file to the OBIGT database. The file must be a CSV (comma separated value) file with column headers that match those in the main database. As an example, here are the contents of BZA10.csv, which has parameters taken from Bazarkina et al. (2010). Missing values are indicated by NA:

file <- system.file("extdata/adds/BZA10.csv", package = "CHNOSZ")
read.csv(file, as.is = TRUE)

##      name abbrv formula state  ref1 ref2       date       G  H     S     Cp
## 1   CdCl+    NA   CdCl+    aq BZA10   NA 2010-07-03  -52629 NA  7.06  11.12
## 2   CdCl2    NA   CdCl2    aq BZA10   NA 2010-07-03  -84883 NA 25.72 116.01
## 3  CdCl3-    NA  CdCl3-    aq BZA10   NA 2010-07-03 -115399 NA 45.15  97.78
## 4 CdCl4-2    NA CdCl4-2    aq BZA10   NA 2010-07-03 -145583 NA 50.61  42.52
##       V    a1.a    a2.b    a3.c    a4.d    c1.e    c2.f omega.lambda z.T
## 1  2.20  2.2303 -2.3357  6.6681 -2.6824 16.6723 -0.7693       0.4372   1
## 2 42.21  7.5221 10.5852  1.5895 -3.2166 73.7023 20.5956      -0.0495   0
## 3 63.47 10.8045 18.5994 -1.5605 -3.5479 72.0244 16.8832       0.9378  -1
## 4 81.35 13.8329 25.9938 -4.4669 -3.8536 53.6766  5.6267       2.4766  -2

Loading the data with add.OBIGT() produces a message that the new data replace existing species. We can then use subcrt() to calculate the equilibrium constant for a reaction involving the new species. Note the decrease in the stepwise stability constant for the second cadmium chloride complex with increasing pressure (Bazarkina et al., 2010, Fig. 4).

iCd <- add.OBIGT(file)

## add.OBIGT: read 4 rows; made 4 replacements, 0 additions [energy units: cal]

subcrt(c("CdCl+", "Cl-", "CdCl2"), c(-1, -1, 1), T = 25, P = c(1, 2000))

## subcrt: 3 species at 2 values of T (ºC) and P (bar) (wet) [energy units: J]

## $reaction
##     coeff  name formula state ispecies
## 392    -1 CdCl+   CdCl+    aq      392
## 29     -1   Cl-     Cl-    aq       29
## 393     1 CdCl2   CdCl2    aq      393
## 
## $out
##    T    P      rho      logK         G       H       S       V      Cp
## 1 25    1 0.997061 0.6413808 -3661.000 2669.61 21.3384 22.6816 561.306
## 2 25 2000 1.071783 0.0618531  -353.058 2732.39 10.4541 12.2632 539.264

After running reset() we can look up the source of data in the default OBIGT database (Sverjensky et al., 1997). Running the reaction now with the default values, we see that the equilibrium constant is not as sensitive to pressure:

reset()

## reset: resetting "thermo" object

## OBIGT: loading default database with 1904 aqueous, 3445 total species

thermo.refs(iCd)[, 1:3]

##      key                                           author year
## 86 SSH97 D. A. Sverjensky, E. L. Shock and H. C. Helgeson 1997

subcrt(c("CdCl+", "Cl-", "CdCl2"), c(-1, -1, 1), T = 25, P = c(1, 2000))

## subcrt: 3 species at 2 values of T (ºC) and P (bar) (wet) [energy units: J]

## $reaction
##     coeff  name formula state ispecies
## 392    -1 CdCl+   CdCl+    aq      392
## 29     -1   Cl-     Cl-    aq       29
## 393     1 CdCl2   CdCl2    aq      393
## 
## $out
##    T    P      rho     logK        G        H        S        V      Cp
## 1 25    1 0.997061 0.619391 -3535.48 -8522.81 -16.7360 10.04043 214.488
## 2 25 2000 1.071783 0.425897 -2431.02 -9603.90 -24.0664  2.34196 191.416

mod.OBIGT() for aqueous species

Let’s add data for CoCl₄^-2 from Liu et al. (2011). The values are taken from Table 5 of that paper; as is common for parameters in the HKF model, they are reported in caloric units, and some of the values have multipliers, which are kept when entering the data. The entry includes the date in ISO 8601 extended format (e.g. 2020-08-16); Sys.Date() is used in this example to get the current date.

mod.OBIGT("CoCl4-2", formula = "CoCl4-2", state = "aq", ref1 = "LBT+11", E_units = "cal",
  date = as.character(Sys.Date()), G = -134150, H = -171558, S = 19.55, Cp = 72.09, V = 27.74)

## mod.OBIGT: updated CoCl4-2(aq)

## [1] 845

The function prints a message saying that the species was added, and returns the species index of the new species. Now let’s modify the new species by adding the HKF coefficients. The z at the end refers to the charge of the species, and is used only for calculating the “g function” in the revised HKF model, not for balancing reactions.

mod.OBIGT("CoCl4-2", a1 = 6.5467, a2 = 8.2069, a3 = 2.0130, a4 = -3.1183,
  c1 = 76.3357, c2 = 11.6389, omega = 2.9159, z = -2)

## mod.OBIGT: no change for CoCl4-2(aq)

## [1] 845

Let us now calculate the equilibrium constant for the formation of CoCl₄^-2 from Co⁺² and Cl^-.

T <- c(25, seq(50, 350, 50))
sres <- subcrt(c("Co+2", "Cl-", "CoCl4-2"), c(-1, -4, 1), T = T)
round(sres$out$logK, 2)

## [1] -3.20 -2.96 -2.02 -0.74  0.77  2.50  4.57  7.29

The calculated values of logK are identical to those in Table 9 of Liu et al. (2011), which provides a good indication that the thermodynamic parameters were entered correctly. Nevertheless, this isn’t a guarantee that the thermodynamic parameters are consistent with the provided values of C_P° and V°. We can see this by running info() to cross-check the parameters for the new CoCl₄^-2 species:

inew <- info("CoCl4-2")
info(inew)

## checkEOS: Cp of CoCl4-2(aq) differs by 1.33 cal K-1 mol-1 from tabulated value

## checkEOS: V of CoCl4-2(aq) differs by -1.19 cm3 mol-1 from tabulated value

The messages indicate that the given values of C_P° and V° differ slightly from those calculated using the HKF parameters.

mod.OBIGT() for minerals

Let’s add data for magnesiochromite from Klemme et al. (2000). The parameters in this paper are reported in Joules, so we set the E.units() to J. The value for volume, in cm³ mol^-1, is from Robie and Hemingway (1995).

H <- -1762000
S <- 119.6
V <- 43.56
mod.OBIGT("magnesiochromite", formula = "MgCr2O4", state = "cr", ref1 = "KOSG00",
          date = as.character(Sys.Date()), E_units = "J", H = H, S = S, V = V)

## mod.OBIGT: added magnesiochromite(cr) with energy units of J

## [1] 3446

Here are the heat capacity parameters for the “Haas-Fisher” polynomial equation (Cp = a + bT + cT^− 2 + dT^− 0.5 + eT²). Order-of-magnitude multipliers are required for the values of b and c (and also e if it is present; see the description for columns 14-21 of the OBIGT data frame in ?thermo). Note that an additional f term is available, which can have any exponent given by lambda, but it is not needed here. 1500 K is a generic value for the high-temperature limit; experimental heat capacities were only reported up to 340 K (Klemme et al., 2000).

a <- 221.4
b <- -0.00102030 * 1e3
c <- -1757210 * 1e-5
d <- -1247.9
mod.OBIGT("magnesiochromite", E_units = "J", a = a, b = b, c = c, d = d,
          e = 0, f = 0, lambda = 0, T = 1500)

## mod.OBIGT: updated magnesiochromite(cr)

## [1] 3446

Now we can use subcrt() to calculate the heat capacity of magnesiochromite. For this calculation, we set the temperature units to Kelvin. We also specify a pressure of 1 bar to prevent an error in the calculation of P_sat below the freezing temperature of water.

T.units("K")

## changed temperature units to K

subcrt("magnesiochromite", property = "Cp", T = c(250, 300, 340), P = 1)

## subcrt: 1 species at 3 values of T (K) and P (bar) [energy units: J]

## $species
##                  name formula state ispecies
## 3446 magnesiochromite MgCr2O4    cr     3446
## 
## $out
## $out$magnesiochromite
##     T P      Cp
## 1 250 1 114.105
## 2 300 1 129.522
## 3 340 1 138.175

The calculated values are with 0.1 J K^-1 mol^-1 of those shown in Fig. 1 of Klemme et al. (2000).

Finally, let’s restore the units setting for later calculations with subcrt(). (Another way would be to run reset(), which also resets the OBIGT database.)

T.units("C")

## changed temperature units to C

Akinfiev-Diamond model

The Akinfiev-Diamond model for aqueous species (Akinfiev and Diamond, 2003) is activated by setting in thermo() for a given aqueous species. The database must also include a corresponding gasesous species with the same name or chemical formula.

TODO

TODO: Add van’t Hoff example?

Fitting formation constants

logB_to_OBIGT() requires three things:

• Experimental decimal logarithms of formation constants (log β) as a function of temperature (at least three different temperatures);
• The stoichiometry of the formation reaction in terms of known species (the new species must be last);
• The experimental temperature and pressure.

logB_to_OBIGT() does three things:

• Combines the formation constants with standard Gibbs energies (ΔG°) of the known species to calculate ΔG° of the new species;
• Fits ΔG° of the new species using a three-term equation with constant heat capacity (i.e., 25 °C G, S, and Cp parameters in OBIGT);
• Adds the fitted parameters to OBIGT for use by other functions in CHNOSZ.

First we set the pressure for all log β data.

P <- "Psat"

Add first species: HWO₄^- (Wang et al., 2019).

T <- c(250, 300, 350)
logB <- c(5.58, 6.51, 7.99)
species <- c("WO4-2", "H+", "HWO4-")
coeff <- c(-1, -1, 1)
logB_to_OBIGT(logB, species, coeff, T, P)

## mod.OBIGT: updated HWO4-(aq)

## logB_to_OBIGT: mean absolute difference between logB (experimental) and logK (calculated) is 0.0199

## [1] 442

Add second species: H₃WO₄F₂^- (Wang et al., 2021).

T <- seq(100, 250, 25)
logB <- c(17.00, 17.11, 17.46, 17.75, 18.17, 18.71, 19.23)
# Species and coefficients in the formation reaction
species <- c("H+", "WO4-2", "F-", "H3WO4F2-")
coeff <- c(-3, -1, -2, 1)
logB_to_OBIGT(logB, species, coeff, T, P)

## mod.OBIGT: updated H3WO4F2-(aq)

## logB_to_OBIGT: mean absolute difference between logB (experimental) and logK (calculated) is 0.0307

## [1] 3447

Add third species: H₂WO₄ (Wang et al., 2021). We have to increase the tolerance because there is considerable scatter in the experimental values.

logB <- c(7.12, 7.82, 7.07, 7.76, 7.59, 7.98, 8.28)
species <- c("H+", "WO4-2", "H2WO4")
coeff <- c(-2, -1, 1)
logB_to_OBIGT(logB, species, coeff, T, P, tolerance = 0.3)

## mod.OBIGT: updated H2WO4(aq)

## logB_to_OBIGT: mean absolute difference between logB (experimental) and logK (calculated) is 0.2035

## [1] 884

We can tell from the low species indices of HWO₄^- (442) and H₂WO₄ (884) that logB_to_OBIGT() replaced parameters for these species that were already present in OBIGT.

Diagram 1: Constant molality of F^-

Now we’re ready to make a speciation diagram. Our target is to reproduce Fig. 7b of Wang et al. (2021), which is made for 300 °C. A constant molality of F^- is based on the assumption of complete dissociation of 0.1 m NaF (we’ll change this later). An ionic strength of 1.2 mol/kg is estimated for a solution with 1.8 m NaCl (use NaCl(300, m_tot = 1.8)). NOTE: because the ionic strength is non-zero, the calculations here refer to molality instead of activity of species (see An Introduction to CHNOSZ).

basis(c("H+", "WO4-2", "F-", "H2O", "O2"))
basis("F-", log10(0.1))
iaq <- retrieve("W", c("O", "H", "F"), "aq")
species(iaq)
a <- affinity(pH = c(2, 7), T = 300, IS = 1.2)
e <- equilibrate(a)
col <- c(1, 4, 5, 2)
diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")

This isn’t quite the diagram we were looking for. The published diagram shows a broad region of coexistence of H₃WO₄F₂^- and HWO₄^- at pH < 5 and increasing abundance of H₂WO₄ at lower pH.

Diagram 2: Variable molality of F^-

In reality, the molality of F^- depends strongly on pH according to the reaction H⁺ + F^- = HF. With a little algebra, we can calculate the molality of F^- (a_F in the code below) from the equilbrium constant of this reaction for a given total F concentration (F_tot). NOTE: It is important to call subcrt() with a non-zero IS so that it returns effective equilibrium constants corrected for ionic strength (see what hapens to the diagram if you set IS = 0).

T <- 300
pH <- seq(2, 7, 0.1)
logK_HF <- subcrt(c("H+", "F-", "HF"), c(-1, -1, 1), T = T, IS = 1.2)$out$logK

## info.character: found HF(aq); also available in gas

## subcrt: 3 species at 300 ºC and 85.84 bar (wet) [energy units: J]

## nonideal: calculations for F- (Bdot equation)

## nonideal: calculations for HF (Setchenow equation)

F_tot <- 0.1
a_F <- F_tot / (1 + 10^(logK_HF - pH))

Now that we have the molality of F^- as a function of pH, we can provide it in the call to affinity().

basis(c("H+", "WO4-2", "F-", "H2O", "O2"))
iaq <- retrieve("W", c("O", "H", "F"), "aq")
species(iaq)
a <- affinity(pH = pH, "F-" = log10(a_F), T = T, IS = 1.2)
e <- equilibrate(a)
diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")

That’s more like it. We have captured the basic geometry of Fig. 7b in Wang et al. (2021). For instance, in accord with the published diagram, HWO₄^- plateaus at around 40% of total W, and H₂WO₄ and H₃WO₄F₂^- are nearly equally abundant at pH = 2.

The highest experimental temperature for the formation constants of H₂WO₄ and H₃WO₄F₂^- is 250 °C, but this diagram is drawn for 300 °C. Wang et al. (2021) used the modified Ryzhenko-Bryzgalin (MRB) model to extrapolate to 300 °C. In contrast, we used a different model but obtained quite similar results.

NOTE: The coefficients in the model used by logB_to_OBIGT() are 25 °C values of G, S, and Cp. In a conservative interpretation, these should be regarded only as fitting parameters, and should not be used to compute thermodynamic properties close to 25 °C unless they were fit to experimental data in that temperature range.