Processing a GWB thermodynamic data file

This vignette shows how to update a GWB thermodynamic data file using log K values computed with the OBIGT database in CHNOSZ.

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

The logKcalc function updates the equilibrium constants for species in a GWB data file. In this example, we convert thermo_24elements.tdat, which is a modfied version of GWB’s thermo.tdat file reduced to 24 elements (Ag, Al, As, Au, C, Ca, Cl, Cu, F, Fe, H, K, Mg, Mn, N, Na, O, P, Pb, S, Si, Sn, U, Zn, the same as in the previous vignette). Some species in the GWB data file may not be present in the default OBIGT database, but optional databases in CHNOSZ can also be used. Here, we use modOBIGT() to add some minerals from SUPCRT92, an entry for steam that has been removed from OBIGT, and antigorite with a chemical formula to match the one used in thermo.tdat.

reset()
modOBIGT(c("addSUPCRT", "steam", "antigorite/2"))

Some basis species used in the GWB file are not in the default OBIGT database. They can be added to OBIGT by fitting their thermodynamic parameters (ΔG°, S° and C_p°) to the log K values for their dissociation reactions from a different GWB file where they have been swapped out of the basis. This depends on the parameters for As(OH)₃ (Perfetti et al., 2008), which were added to OBIGT in CHNOSZ 1.4.0, so they are defined here if needed for an earlier version.

if(packageVersion("CHNOSZ") < "1.4.0") mod.obigt("As(OH)3", E_units = "J",
  G = -639800, H = -737300, S = 212.4, Cp = 85, V = 47.5,
  a1 = 45.2, a2 = 16.1, a3 = 43.7, a4 = -18, c1 = 186.2, c2 = -47.7, omega = 0.52, z = 0)
addOBIGT("As(OH)4-")
addOBIGT("Sn++++")

Here we set the file paths and perform the log K calculation. The maxprint argument is used to prevent the printing of long lists of species that can’t be found in OBIGT.

infile <- system.file("extdata/thermo_24elements.tdat", package = "logKcalc")
outfile <- file.path(tempdir(), "thermo_OBIGT.tdat")
logKcalc(infile, outfile, 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] "0, 25, 60, 100, 150, 200, 250, 300"

## The pressures for the logK calculation are:

## [1] "1, 1, 1, 1.0133, 4.7573, 15.5366, 39.7366, 85.8379"

## The 25 available basis species have 25 elements:

## [1] "Ag, Al, As, Au, C, Ca, Cl, Cu, F, Fe, H, K, Mg, Mn, N, Na, O, P, Pb, S, Si, Sn, U, Z, Zn"

## 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-, Ca(O-phth), H(O-phth)-, H2(O-phth), Na(O-phth)-"

## Can't map 162 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 211 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] "Amesite-14A, Daphnite-7A, Heulandite, Minnesotaite, Pyrrhotite"

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

## [1] "25 basis species"
## [1] "19 redox species"
## [1] "106 aqueous species"
## [1] "1 electron species"
## [1] "131 mineral species"
## [1] "8 gas species"
## [1] "21 oxide species"

Note that species are removed if their log K values are all NA. This includes a few minerals taken from SUPCRT92, as well as pyrrhotite, whose chemical formula in the GWB file (Fe.875S) is different from that in OBIGT (FeS).

Now we can make some plots to compare the log K values in the input and output files.

Only species with the same name and dissociation reaction in both files are considered for the plot. The plot titles indicate how many of these species have log K values that are not available (NA) at this particular temperature and pressure (indicated by a “500” in the GWB file). Labels are added to points with a difference of greater than 1 log K unit.

lab1 <- "thermo.tdat"
lab2 <- "OBIGT"
plot1 <- logKcomp(infile, outfile, "aqueous", 4, lab1, lab2)
plot2 <- logKcomp(infile, outfile, "mineral", 4, lab1, lab2, c(-50, 100), c(NA, 5))
gridExtra::grid.arrange(plot1, plot2, ncol = 2)

OBIGT has generally lower log K values for dissociation reactions of minerals, which is similar to the case for the UNITHERM database. We also find a larger stability field for kaolinite in the solubility diagram made using thermo_OBIGT.tdat compared to the default thermo.tdat in GWB.

See the next vignette for an example of adding new species to the GWB thermodynamic data file.