library(logKcalc)
library(knitr)
## colorize messages 20171031
## adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
  function(x, options) sprintf('<pre style="color:%s">%s</pre>', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
## use pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL 
## logK with a thin space 20200627
logK <- "log&thinsp;<i>K</i>"

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

This vignette was compiled on r Sys.Date() with logKcalc r packageDescription("logKcalc")$Version and CHNOSZ r packageDescription("CHNOSZ")$Version.

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 r logK 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)~3~ [@PPB+08], which were added to OBIGT in CHNOSZ 1.4.0, so they are defined here if needed for an earlier version.

addOBIGT("As(OH)4-")
addOBIGT("Sn++++")

Here we set the file paths and perform the r logK 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)

Note that species are removed if their r logK 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 r logK 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 r logK 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 r logK 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 r logK 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.

include_graphics("Act2_vig2.png")

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

References



jedick/logKcalc documentation built on Jan. 19, 2025, 10:41 a.m.