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>"
## run reset() here so we don't use the DEW model (from vig4.Rmd)
reset()

This vignette shows how to use data for species from a GWB file in a diagram made with CHNOSZ.

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

One of the example scripts for the Act2 program in The Geochemist's Workbench® is a mosaic diagram for the As-S-O-H system [@GWB20a p.98]. A similar diagram can be produced using the OBIGT database and CHNOSZ, but some differences are apparent because of different thermodynamic data compared to the thermo.tdat data file in GWB. In particular, the default OBIGT database has parameters for AsH~3~ from @SSW01, aqueous As(OH)~3~ and AsO(OH)~3~ from @PPB+08, and ionized As-bearing species from @NA03, while thermo.tdat is based on the LLNL (Lawrence Livermore National Laboratory) compilation from various earlier sources.

This diagram is made using the default OBIGT database. However, the parameters for As(OH)~3~ are defined here if needed (they were added to OBIGT in CHNOSZ 1.4.0).


The major differences from the diagram made using thermo.tdat in GWB are a larger As(OH)~3~ field, the absence of orpiment and presence of native arsenic, and different aqueous speciation under relatively reducing conditions at high pH.

In order to more closely reproduce the example from GWB, we can use addOBIGT() to fit thermodynamic parameters to the r logK values of dissociation reactions of species from a GWB data file. This function uses a modified version of GWB's thermo.tdat in which the basis species are all available in OBIGT (thermo_24swapped.tdat). The fitted parameters (ΔG°, S° and C~p~°) of the formed species are added to the OBIGT database, so they can be used in diagrams made in CHNOSZ.

Here, we fit the parameters for minerals and aqueous species other than As(OH)~3~, which is a basis species in thermo_24swapped.tdat. Although it has no stability field in the diagram, AsH~3~ is also updated because it appears in the dissociation reaction for realgar. Because it uses just three parameters, the fit may introduce some error, so we increase the tolerance (in r logK units) for H~2~AsO~4~^-^ to allow the checks for this species to pass. We run reset() first to make sure we are starting with the default OBIGT database and add As(OH)~3~ again if needed.

reset()
addOBIGT("AsH3(aq)")
addOBIGT("H2AsO4-", tolerance = 0.15)
addOBIGT("HAsO4--")
addOBIGT("As(OH)4-")
addOBIGT("AsO2OH--")
addOBIGT("AsO4---")
addOBIGT("Orpiment", "As2S3")
addOBIGT("Realgar", "AsS")

The warning about an unbalanced reaction occurs because of the inexact representation of reaction coefficients in the GWB file due to rounding, and can be ignored.

For completeness we can even use fitted thermodynamic parameters for the S-bearing basis species that are speciated across the diagram. This does not include SO~4~^-2^, since it is a basis species in thermo_24swapped.tdat. For this example only, we include S^-2^, which has been recommended to be removed from thermodynamic databases [@MBH+18] and is not present in the default OBIGT database.

addOBIGT("HSO4-")
addOBIGT("H2S(aq)")
addOBIGT("HS-")
addOBIGT("S--")

Note that repeated '-' or '+' signs at the end of formulas are converted into e.g. '-2' for compatibility with the OBIGT database.

Next we identify As-bearing species in the As-S-O-H system. We remove species with data from @NA03, as they may conflict with the GWB dataset. However, we leave As(OH)~3~ in the system; although its thermodynamic parameters in OBIGT are still from @PPB+08, the parameters of the other species were made to be consistent with the r logK values for the dissociation reactions in thermo_24swapped.tdat that involve As(OH)~3~ as a basis species.

iaq <- retrieve("As", c("S", "O", "H"), "aq")
iaq <- iaq[!grepl("NA03", info(iaq)$ref1)]
iaq <- c(info("As(OH)3"), iaq)
icr <- retrieve("As", c("S", "O", "H"), "cr")
icr <- icr[!grepl("NA03", info(icr)$ref1)]

Now we write some code to set up the system in CHNOSZ and make a loga~O2~--pH mosaic diagram at 100 °C. The activity of As-bearing aqueous species is set to 10^-3^ and that of SO~4~^-2^ to 10^-4^ (this value applies to the S-bearing basis species as a whole).

# Set up basis species and system
basis(c("As(OH)3", "SO4-2", "H2O", "O2", "H+"))
basis("SO4-2", -4)
logact <- c(rep(-3, length(iaq)), rep(0, length(icr)))
species(c(iaq, icr), logact)
# Set pH and loga(O2) range and grid resolution
pH <- c(0, 14, 500)
O2 <- c(-70, -40, 500)
# Calculate reaction affinities for changing basis species
bases <- c("SO4-2", "HSO4-", "H2S", "HS-", "S-2")
m <- mosaic(bases, pH = pH, O2 = O2, T = 100, blend = TRUE)
# Make the diagram
par(cex = 2.5, mar = c(2.8, 3, 1.4, 1))
# Web colors lightcyan and lavenderblush with some transparency
fill <- c(rep("#E0FFFF88", length(iaq)), rep("#FFF0F588", length(icr)))
d <- diagram(m$A.species, fill = fill, lwd = 2)
diagram(m$A.bases, add = TRUE, col = 4, lty = 3, lwd = 2,
        col.names = 4, cex.names = 0.7, italic = TRUE)
water.lines(d)
legend("bottomleft", describe.property("T", 100), bty = "n")
title(main = "blend = TRUE", font.main = 1)

The diagram is made with the default setting of blend = TRUE to represent the continuous change of speciation between basis species. Using blend = FALSE would instead portray abrupt transitions between basis species, similar to the diagram from the Mosaic.ac2 example in GWB.


References



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