[CHNOSZ-commits] r653 - in pkg/CHNOSZ: . R demo vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 19 09:35:15 CET 2021
Author: jedick
Date: 2021-03-19 09:35:15 +0100 (Fri, 19 Mar 2021)
New Revision: 653
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
Update anintro.Rmd and demo/gold.R for new solubility() arguments
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2021-03-19 07:50:28 UTC (rev 652)
+++ pkg/CHNOSZ/DESCRIPTION 2021-03-19 08:35:15 UTC (rev 653)
@@ -1,6 +1,6 @@
Date: 2021-03-19
Package: CHNOSZ
-Version: 1.4.0-22
+Version: 1.4.0-23
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2021-03-19 07:50:28 UTC (rev 652)
+++ pkg/CHNOSZ/R/solubility.R 2021-03-19 08:35:15 UTC (rev 653)
@@ -30,6 +30,7 @@
logact <- basis()$logact
# The current formed species are the minerals to be dissolved
mineral <- species()
+ if(is.null(mineral)) stop("please load minerals or gases with species()")
# Make a list to store the calculated solubilities for each mineral
slist <- list()
@@ -40,7 +41,11 @@
# Define basis species with the mineral first (so it will be dissolved)
ispecies[1] <- mineral$ispecies[i]
logact[1] <- mineral$logact[i]
- basis(ispecies, logact)
+ # Use numeric values first and put in buffer names second (needed for demo/gold.R)
+ loga.numeric <- suppressWarnings(as.numeric(logact))
+ basis(ispecies, loga.numeric)
+ is.na <- is.na(loga.numeric)
+ if(any(is.na)) basis(rownames(basis())[is.na], logact[is.na])
# Add aqueous species (no need to define activities here - they will be calculated by solubility_calc)
species(iaq)
if(is.mosaic) a <- suppressMessages(mosaic(...)$A.species) else a <- suppressMessages(affinity(...))
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2021-03-19 07:50:28 UTC (rev 652)
+++ pkg/CHNOSZ/demo/gold.R 2021-03-19 08:35:15 UTC (rev 653)
@@ -4,23 +4,23 @@
# 20190127 update Au species in OBIGT, not here
library(CHNOSZ)
-# set up system
-# use H2S here: it's the predominant species at the pH of the QMK buffer -- see sulfur()
-basis(c("Al2O3", "quartz", "Fe", "Au", "K+", "Cl-", "H2S", "H2O", "oxygen", "H+"))
+# Set up system
+# Use H2S here: it's the predominant species at the pH of the QMK buffer -- see sulfur()
+basis(c("Au", "Al2O3", "quartz", "Fe", "K+", "Cl-", "H2S", "H2O", "oxygen", "H+"))
# set molality of K+ in completely dissociated 0.5 molal KCl
# NOTE: This value is used only for making the legend;
# activities corrected for ionic strength are computed below
basis("K+", log10(0.5))
-# create a pH buffer
+# Create a pH buffer
mod.buffer("QMK", c("quartz", "muscovite", "K-feldspar"), "cr", 0)
-# define colors for Au(HS)2-, AuHS, AuOH, AuCl2-
+# Define colors for Au(HS)2-, AuHS, AuOH, AuCl2-
# after Williams-Jones et al., 2009
# (doi:10.2113/gselements.5.5.281)
col <- c("#ED4037", "#F58645", "#0F9DE2", "#22CC88")
-# sulfur logfO2-pH diagrams showing redox and pH buffers at four temperatures 20181031
+# Sulfur logfO2-pH diagrams showing redox and pH buffers at four temperatures 20181031
sulfur <- function() {
species(c("H2S", "HS-", "HSO4-", "SO4-2"))
T <- c(200, 300, 400, 500)
@@ -51,21 +51,21 @@
# log(m_Au)-pH diagram like Fig. 7 of Akinfiev and Zotov, 2001
# (http://pleiades.online/cgi-perl/search.pl/?type=abstract&name=geochem&number=10&year=1&page=990)
Au_pH1 <- function() {
- species(c("Au(HS)2-", "AuHS", "AuOH"))
- # apply PPM buffer for fO2 and aH2S
+ # Apply PPM buffer for fO2 and aH2S
basis("O2", "PPM")
basis("H2S", "PPM")
- # calculate affinity and solubility
+ # Calculate solubility of gold
+ species("Au")
+ iaq <- c("Au(HS)2-", "AuHS", "AuOH")
# (set IS = 0 for diagram to show "log m" instead of "log a")
- a <- affinity(pH = c(3, 8), T = 300, P = 1000, IS = 0)
- s <- solubility(a)
- # make diagram and show total log molality
+ s <- solubility(iaq, pH = c(3, 8), T = 300, P = 1000, IS = 0)
+ # Make diagram and show total log molality
diagram(s, ylim = c(-10, -5), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
- # add neutral pH line
+ # Add neutral pH line
pH <- -subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = 300, P = 1000)$out$logK/2
abline(v = pH, lty = 3)
- # make legend and title
+ # Make legend and title
dprop <- describe.property(c("T", "P", "IS"), c(300, 1000, 0))
legend("topleft", dprop, bty = "n")
dbasis <- describe.basis(ibasis = c(9, 7))
@@ -77,20 +77,20 @@
# (doi:10.1016/j.gca.2004.04.006)
Au_pH2 <- function() {
species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
- # apply PPM buffer for fO2 and aH2S
+ # Apply PPM buffer for fO2 and aH2S
basis("O2", "PPM")
basis("H2S", "PPM")
- # calculate affinity and solubility
+ # Calculate affinity and solubility
# (set IS = 0 for diagram to show "log m" instead of "log a")
a <- affinity(pH = c(3, 8), T = 450, P = 1000, IS = 0)
s <- solubility(a)
- # make diagram and show total log molality
+ # Make diagram and show total log molality
diagram(s, ylim = c(-8, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
- # add neutral pH line
+ # Add neutral pH line
pH <- -subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = 450, P = 1000)$out$logK/2
abline(v = pH, lty = 3)
- # make legend and title
+ # Make legend and title
dprop <- describe.property(c("T", "P", "IS"), c(450, 1000, 0))
legend("topleft", dprop, bty = "n")
dbasis <- describe.basis(ibasis = c(6, 9, 7))
@@ -98,14 +98,14 @@
title(main=("After Stef\u00e1nsson and Seward, 2004, Fig. 12b"), font.main = 1, cex.main = 1.1)
}
-# estimate the Cl- molality and ionic strength for a hypothetical
+# Dstimate the Cl- molality and ionic strength for a hypothetical
# NaCl solution with total chloride equal to specified NaCl + KCl solution,
# then estimate the molality of K+ in that solution 20181109
chloride <- function(T, P, m_NaCl, m_KCl) {
NaCl <- NaCl(T = T, P = P, m_tot = m_NaCl + m_KCl)
- # calculate logK of K+ + Cl- = KCl, adjusted for ionic strength
+ # Calculate logK of K+ + Cl- = KCl, adjusted for ionic strength
logKadj <- subcrt(c("K+", "Cl-", "KCl"), c(-1, -1, 1), T = T, P = P, IS = NaCl$IS)$out$logK
- # what is the molality of K+ from 0.5 mol KCl in solution with 2 mol total Cl
+ # What is the molality of K+ from 0.5 mol KCl in solution with 2 mol total Cl
m_K <- m_KCl / (10^logKadj * NaCl$m_Cl + 1)
list(IS = NaCl$IS, m_Cl = NaCl$m_Cl, m_K = m_K)
}
@@ -114,20 +114,20 @@
# (doi:10.2113/gselements.5.5.281)
Au_T1 <- function() {
species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
- # apply PPM buffer for fO2 and aH2S
+ # Apply PPM buffer for fO2 and aH2S
basis("O2", "PPM")
basis("H2S", "PPM")
- # apply QMK buffer for pH
+ # Apply QMK buffer for pH
basis("H+", "QMK")
- # estimate solution composition for 1.5 m NaCl and 0.5 m KCl
+ # Estimate solution composition for 1.5 m NaCl and 0.5 m KCl
chl <- chloride(T = seq(150, 550, 10), P = 1000, m_NaCl = 1.5, m_KCl = 0.5)
- # calculate affinity and solubility
+ # Calculate affinity and solubility
a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
s <- solubility(a)
- # make diagram and show total log molality
+ # Make diagram and show total log molality
diagram(s, ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
- # make legend and title
+ # Make legend and title
dP <- describe.property("P", 1000)
dNaCl <- expression(italic(m)[NaCl] == 1.5)
dKCl <- expression(italic(m)[KCl] == 0.5)
@@ -144,24 +144,24 @@
# (doi:10.1144/SP402.4)
Au_T2 <- function() {
species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
- # total S = 0.01 m
+ # Total S = 0.01 m
basis("H2S", -2)
- # apply HM buffer for fO2
+ # Apply HM buffer for fO2
basis("O2", "HM")
- # apply QMK buffer for pH
+ # Apply QMK buffer for pH
basis("H+", "QMK")
- # estimate solution composition for 1.5 m NaCl and 0.5 m KCl
+ # Estimate solution composition for 1.5 m NaCl and 0.5 m KCl
chl <- chloride(T = seq(150, 550, 10), P = 1000, m_NaCl = 1.5, m_KCl = 0.5)
-# # calculate affinity and solubility, considering speciation of sulfur
+# # Calculate affinity and solubility, considering speciation of sulfur
# bases <- c("H2S", "HS-", "SO4-2", "HSO4-")
# m <- mosaic(bases, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
# s <- solubility(m$A.species)
a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
s <- solubility(a)
- # make diagram and show total log molality
+ # Make diagram and show total log molality
diagram(s, ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
- # make legend and title
+ # Make legend and title
dP <- describe.property("P", 1000)
dNaCl <- expression(italic(m)[NaCl] == 1.5)
dKCl <- expression(italic(m)[KCl] == 0.5)
@@ -173,7 +173,7 @@
title(main=("After Williams-Jones et al., 2009, Fig. 2A"), font.main = 1)
}
-# make plots
+# Make plots
opar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
Au_pH1()
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2021-03-19 07:50:28 UTC (rev 652)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2021-03-19 08:35:15 UTC (rev 653)
@@ -906,26 +906,30 @@
Although this assumption of metastable equilibrium is useful for making many types of diagrams for aqueous species, the aqueous solutions would not be in equilibrium with other phases, including gases or solids such as carbon dioxide or calcite.
Moving away from metastable equilibrium to complete equilibrium actually involves a simplification of the computations.
-Instead of finding the activities of aqueous species where the affinities of formation reactions are equal to each other, equilibrium corresponds to the configuration where the affinities are equal to zero.
-However, this simplification comes with a different constraint: the reactions should be balanced on something that anchors them to a real quantity.
-It makes sense to think of this quantity as a conserved basis species (one that is present in the formation reactions of all species considered), which corresponds to a pure substance that is being dissolved.
+Instead of finding the activities of aqueous species where the affinities of formation reactions are equal to each other, equilibrium occurs when the affinities of all formation reactions are equal to zero.
+However, the solubility calculation comes with a different constraint: all reactions between the pure substance that is being dissolved and any of the possible formed aqueous species must be able to be balanced on the same element (usually the main metal in the system).
+This balancing constraint can be expressed as a conserved basis species, i.e. one that is present in non-zero quantity in the formation reactions of all species considered.
-The <span style="color:green">`solubility()`</span> function provides a way to compute activities of aqueous species in equilibrium with a solid or gas, which is usually defined as the first basis species.
+The <span style="color:green">`solubility()`</span> function provides a way to compute activities of aqueous species in equilibrium with one or more minerals or gases.
The following example for corundum (Al<sub>2</sub>O<sub>3</sub>) is based on Figure 15 of @Man13.
-To reproduce the diagram, we use superseded thermodynamic data that are kept in the `SLOP98` optional data file.
-Corundum is set as the first basis species, and the formed species all contain Al.
-The affinities of the formation reactions, for unit activity of the solid basis species (corundum) and any given activities of the formed aqueous species, are input to <span style="color:green">`solubility()`</span>.
+To more closely reproduce the diagram, we use superseded thermodynamic data that are kept in the `SLOP98` optional data file.
+
+The basis species are defined with the main element (Al) first.
+The mineral we want to dissolve, corundum, is loaded as a formed species.
+The aqueous species that can form by dissolution of corundum are listed in the `iaq` argument of <span style="color:green">`solubility()`</span>.
+The next arguments describe the variables for the <span style="color:green">`affinity()`</span> calculations.
+Note that setting `IS` to 0 has no effect on the calculations, but signals <span style="color:green">`diagram()`</span> to label the *y* axis with logarithm of molality instead of logarithm of activity.
An additional argument, `in.terms.of`, is used to compute the total molality of Al in solution, that is, twice the number of moles of Al<sub>2</sub>O<sub>3</sub> that are dissolved.
+
<span style="color:green">`diagram()`</span> is used twice, first to plot the total molality of Al, then the concentrations of the individual species, using `adj` and `dy` to adjust the positions of labels in the *x*- and *y*-directions.
-Note that setting `IS` to 0 in <span style="color:green">`affinity()`</span> has no effect on the calculations, but signals <span style="color:green">`diagram()`</span> to label the *y* axis with logarithm of molality instead of logarithm of activity.
At the end of the calculation, we use <span style="color:red">`reset()`</span> to restore the default thermodynamic database.
```{r corundum, fig.margin=TRUE, fig.width=4, fig.height=4, dpi=dpi, out.width="100%", results="hide", message=FALSE, cache=TRUE, fig.cap="Solubility of corundum (green line) and equilibrium concentrations of aqueous species (black lines).", pngquant=pngquant, timeit=timeit}
add.OBIGT("SLOP98")
-basis(c("corundum", "H2O", "H+", "O2"))
-species(c("Al+3", "AlO2-", "AlOH+2", "AlO+", "HAlO2"))
-a <- affinity(pH = c(0, 10), IS = 0)
-s <- solubility(a, in.terms.of = "Al+3")
+basis(c("Al+3", "H2O", "H+", "O2"))
+species("corundum")
+iaq <- c("Al+3", "AlO2-", "AlOH+2", "AlO+", "HAlO2")
+s <- solubility(iaq, pH = c(0, 10), IS = 0, in.terms.of = "Al+3")
diagram(s, type = "loga.balance", ylim = c(-10, 0), lwd = 4, col = "green3")
diagram(s, add = TRUE, adj = c(0, 1, 2.1, -0.2, -1.5), dy = c(0, 0, 4, -0.3, 0.1))
legend("topright", c("25 °C", "1 bar"), text.font = 2, bty = "n")
@@ -934,8 +938,8 @@
Other examples of using <span style="color:green">`solubility()`</span> are available in CHNOSZ.
See [<span style="color:blue">`demo(solubility)`</span>](../demo) for calculations of the solubility of CO<sub>2(<i>gas</i>)</sub> and calcite as a function of pH and temperature.
-The calculation assumes the stoichiometric dissolution of calcite, in which CaCO<sub>3</sub> dissociates to form equal quantities of Ca<sup>+2</sup> and CO<sub>3</sub><sup>-2</sup> ions.
-Adding in activity coefficients, a different example in <span style="color:blue">`?solubility`</span> uses the `find.IS` option to find the final ionic strength for dissolving calcite into pure water.
+The calculation considers the stoichiometric dissolution of calcite, in which CaCO<sub>3</sub> dissociates to form equal quantities of Ca<sup>+2</sup> and CO<sub>3</sub><sup>-2</sup> ions.
+Adding in activity coefficients, an example in <span style="color:blue">`?solubility`</span> uses the `find.IS` option to find the final ionic strength for dissolving calcite into pure water.
[<span style="color:blue">`demo(gold)`</span>](../demo) shows calculations of the solubility of gold as a function of pH and *T* as well as oxygen fugacity set by diferent mineral buffers, and considers ionic strength effects on activity coefficients, so that activities are transformed to molalities ([see below](#transformation-of-variables)).
# Activity coefficients
More information about the CHNOSZ-commits
mailing list