[CHNOSZ-commits] r647 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 3 10:57:49 CET 2021
Author: jedick
Date: 2021-03-03 10:57:48 +0100 (Wed, 03 Mar 2021)
New Revision: 647
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/demo/Pourbaix.R
pkg/CHNOSZ/demo/zinc.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/TODO
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/solubility.Rd
Log:
Add solubilities(), used in demo/Pourbaix.R and demo/zinc.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/DESCRIPTION 2021-03-03 09:57:48 UTC (rev 647)
@@ -1,6 +1,6 @@
-Date: 2021-03-02
+Date: 2021-03-03
Package: CHNOSZ
-Version: 1.4.0-16
+Version: 1.4.0-17
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/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/NAMESPACE 2021-03-03 09:57:48 UTC (rev 647)
@@ -58,7 +58,9 @@
"CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "AkDi", "moles",
"lNaCl", "lS", "lT", "lP", "lTP", "lex",
# added 20200716 or later
- "mash", "mix", "rebalance"
+ "mash", "mix", "rebalance",
+# added 20210303 or later
+ "solubilities"
)
# Load shared objects
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/R/diagram.R 2021-03-03 09:57:48 UTC (rev 647)
@@ -41,10 +41,13 @@
### argument handling ###
- ## check that eout is a valid object
+ ## Check that eout is a valid object
efun <- eout$fun
if(length(efun)==0) efun <- ""
- if(!(efun %in% c("affinity", "equilibrate") | grepl("solubility", efun))) stop("'eout' is not the output from affinity(), equilibrate(), or solubility()")
+ if(!(efun %in% c("affinity", "equilibrate") | grepl("solubilit", efun)))
+ stop("'eout' is not the output from affinity(), equilibrate(), solubility(), or solubilities()")
+ # For solubilities(), type is always loga.balance 20210303
+ if(identical(efun, "solubilities")) type <- "loga.balance"
## 'type' can be:
# 'auto' - property from affinity() (1D) or maximum affinity (affinity 2D) (aout) or loga.equil (eout)
@@ -51,7 +54,7 @@
# 'loga.equil' - equilibrium activities of species of interest (eout)
# name of basis species - equilibrium activity of a basis species (aout)
# 'saturation' - affinity=0 line for each species (2D)
- # 'loga.balance' - activity of balanced basis species (eout from solubility())
+ # 'loga.balance' - activity of balanced basis species (eout from solubility() or solubilities())
eout.is.aout <- FALSE
plot.loga.basis <- FALSE
if(type %in% c("auto", "saturation")) {
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/R/examples.R 2021-03-03 09:57:48 UTC (rev 647)
@@ -33,7 +33,7 @@
demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "findit", "ionize", "buffer", "protbuff", "glycinate",
"mosaic", "copper", "arsenic", "solubility", "gold", "contour", "sphalerite", "zinc",
- "Shh", "saturation", "adenine", "DEW", "lambda", "potsassium", "TCA", "aluminum",
+ "Shh", "saturation", "adenine", "DEW", "lambda", "potassium", "TCA", "aluminum",
"AkDi", "comproportionation", "Pourbaix"), save.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
Modified: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/R/solubility.R 2021-03-03 09:57:48 UTC (rev 647)
@@ -181,3 +181,50 @@
c(aout, list(balance=bout$balance, m.balance=bout$n.balance, n.balance=bout$n.balance, in.terms.of=in.terms.of,
loga.balance=loga.balance, Astar=loga.equil, loga.equil=loga.equil))
}
+
+# Calculate solubilities of multiple minerals 20210303
+# a_cr: affinities for minerals (all bearing the same metal)
+# i_aq: aqueous species that can be produced by dissolution of the minerals
+# FIXME: what to do about 'dissociation' argument?
+solubilities <- function(a_cr, i_aq, in.terms.of = NULL, dissociation = NULL) {
+ # If a_cr is the output from mosaic(), just get the species' affinities
+ is.mosaic <- FALSE
+ m_cr <- a_cr
+ if(identical(a_cr$fun, "mosaic")) {
+ a_cr <- a_cr$A.species
+ is.mosaic <- TRUE
+ }
+ # Find all stable minerals across diagram
+ d_cr <- diagram(a_cr, plot.it = FALSE)
+ d_cr.stable <- d_cr$species$ispecies[unique(as.vector(d_cr$predominant))]
+ # Use basis species in a_cr as a template for the solubility calculations
+ ispecies <- a_cr$basis$ispecies
+ logact <- a_cr$basis$logact
+
+ # Make a list to store the calculated solubilities for each mineral
+ slist <- list()
+ # Loop over stable minerals
+ for(i in seq_along(d_cr.stable)) {
+ # Define basis species with the mineral first (so it will be dissolved)
+ ispecies[1] <- d_cr.stable[i]
+ basis(ispecies, logact)
+ # Add aqueous species (no need to define activities here - they will be calculated)
+ species(i_aq)
+ # Calculate affinities of formation reactions of this mineral at same conditions as a_cr (argument recall)
+ if(is.mosaic) a <- mosaic(m_cr)$A.species else a <- affinity(a_cr)
+ # Calculate solubility of this mineral
+ s <- solubility(a, in.terms.of = in.terms.of, dissociation = dissociation)
+ # Store the solubilities in the list
+ slist[[i]] <- s$loga.balance
+ }
+
+ # The overall solubility is the *minimum* among all the minerals
+ smin <- do.call(pmin, slist)
+ # Put this into the last-computed 'solubility' object
+ s$loga.balance <- smin
+ # Change the function name stored in the object so diagram() plots loga.balance automatically
+ s$fun <- "solubilities"
+ # Return the object
+ s
+}
+
Modified: pkg/CHNOSZ/demo/Pourbaix.R
===================================================================
--- pkg/CHNOSZ/demo/Pourbaix.R 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/demo/Pourbaix.R 2021-03-03 09:57:48 UTC (rev 647)
@@ -83,35 +83,11 @@
# Calculate affinities for minerals
species(i_cr)
a_cr <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P, IS = IS)
+# Calculate overall solubility (i.e. minimum solubility given all candidate minerals)
+s <- solubilities(a_cr, i_aq, in.terms.of = element)
-# Find all stable minerals across diagram
-d_cr <- diagram(a_cr, plot.it = FALSE)
-d_cr.stable <- d_cr$species$name[unique(as.vector(d_cr$predominant))]
-
-# Make a list to store the calculated solubilities for each mineral
-slist <- list()
-# Loop over stable minerals
-for(i in seq_along(d_cr.stable)) {
- # Define basis species with mineral to dissolve
- basis(c(d_cr.stable[i], "H2O", "H+", "e-"))
- # Add aqueous species (no need to define activities here - they will be calculated)
- species(i_aq)
- # Calculate affinities of formation reactions
- a <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P, IS = IS)
- # Calculate solubility of this mineral
- # FIXME: what to do about 'dissociation' argument?
- s <- solubility(a, in.terms.of = element, dissociation = FALSE)
- # Store the solubilities in the list
- slist[[i]] <- s$loga.balance
-}
-
-# The overall solubility is the *minimum* among all the minerals
-smin <- do.call(pmin, slist)
-# Put this into the last-computed 'solubility' object
-s$loga.balance <- smin
-
# Plot diagram (LAYER 2: equisolubility lines)
-diagram(s, type = "loga.balance", levels = levels, contour.method = "flattest", add = TRUE, lwd = 1.7)
+diagram(s, levels = levels, contour.method = "flattest", add = TRUE, lwd = 1.5)
# Calculate affinities for aqueous species
# FIXME: should be able to remove cr species from previous affinity object
@@ -167,3 +143,4 @@
Pexpr <- lP(P)
main <- bquote(.(element)*"-O-H at "*.(Texpr)*" and "*.(Pexpr))
title(main = main)
+
Modified: pkg/CHNOSZ/demo/zinc.R
===================================================================
--- pkg/CHNOSZ/demo/zinc.R 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/demo/zinc.R 2021-03-03 09:57:48 UTC (rev 647)
@@ -1,5 +1,5 @@
# CHNOSZ/demo/zinc.R
-# make Zn solubility diagram with multiple minerals
+# Make Zn solubility diagram with multiple minerals
# 20190530 jmd first version (plot_Zn.R)
# 20201008 combine solubility contours for different minerals
# 20201014 added to CHNOSZ
@@ -8,7 +8,7 @@
opar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
-# system variables
+# System variables
res <- 300
T <- 100
P <- "Psat"
@@ -15,15 +15,15 @@
Stot <- 1e-3
pH <- c(0, 14, res)
O2 <- c(-62, -40, res)
-# mass percent NaCl in saturated solution at 100 degC, from CRC handbook
+# Mass percent NaCl in saturated solution at 100 degC, from CRC handbook
w2 <- 0.2805
-# set up basis species
+# Set up basis species
basis(c("ZnO", "H2S", "Cl-", "oxygen", "H2O", "H+"))
basis("H2S", log10(Stot))
-# molality of NaCl in saturated solution
+# Molality of NaCl in saturated solution
m2 <- 1000 * w2 / (mass("NaCl") * (1 - w2))
-# estimate ionic strength and molality of Cl-
+# Estimate ionic strength and molality of Cl-
sat <- NaCl(T = 100, m_tot = m2)
basis("Cl-", log10(sat$m_Cl))
@@ -66,34 +66,15 @@
diagram(mcr$A.species, col = 2)
label.figure("B")
-## Solubility plot 20201008
-
-# Make a list to store the calculated solubilities for each mineral
-slist <- list()
-# Loop over minerals
-minerals <- c("zincite", "sphalerite")
-for(i in seq_along(minerals)) {
- # Define basis species with mineral to dissolve
- basis(c(minerals[i], "H2S", "Cl-", "oxygen", "H2O", "H+"))
- basis("H2S", log10(Stot))
- basis("Cl-", log10(sat$m_Cl))
- # Add aqueous species (no need to define activities here - they will be calculated)
- species(iaq)
- # Calculate affinities of formation reactions, using mosaic() to speciate the S-bearing basis species
- m <- mosaic(bases, pH = pH, O2 = O2, T = T, IS = sat$IS)
- # Calculate solubility of this mineral
- s <- solubility(m$A.species, in.terms.of = "Zn", dissociation = FALSE)
- # Store the solubilities in the list
- slist[[i]] <- s$loga.balance
-}
-
-# The overall solubility is the *minimum* among all the minerals
-smin <- do.call(pmin, slist)
-# Put this into the last-computed 'solubility' object
-s$loga.balance <- smin
+# Calculate *minimum* solubility among all the minerals 20201008
+# (i.e. saturation condition for the solution)
+# Use solubilities() 20210303
+# FIXME: why do we need to set dissociation = FALSE here?
+s <- solubilities(mcr, iaq, in.terms.of = "Zn", dissociation = FALSE)
# Specify contour levels
levels <- seq(-12, 9, 3)
diagram(s, type = "loga.balance", levels = levels, contour.method = "flattest")
+
# Show the mineral stability boundaries
diagram(mcr$A.species, names = NA, add = TRUE, lty = 2, col = 2)
title("Solubilities of 2 minerals", font.main = 1, line = 1.5)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2021-03-03 09:57:48 UTC (rev 647)
@@ -10,26 +10,40 @@
\newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
\newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
-\section{Changes in CHNOSZ version 1.4.0-16 (2021-03-02)}{
+\section{Changes in CHNOSZ version 1.4.0-17 (2021-03-03)}{
- \subsection{CHANGES}{
+ \subsection{NEW FEATURES}{
\itemize{
- \item Remove limSolve from Suggests in DESCRIPTION.
+ \item Add \samp{demo/Pourbaix.R} (Eh-pH diagrams with isosolubility
+ lines, after \href{https://www.worldcat.org/oclc/563921897}{Pourbaix,
+ 1974}). This demo depends on the following two changes.
+ \item Add \code{solubilities()} for calculating overall (i.e. minimum)
+ solubility for multiple minerals. This function is also now used in
+ \samp{demo/zinc.R}.
+
+ \item Restore \samp{lty.aq} and \samp{lty.cr} arguments to
+ \code{diagram()} to control plotting of aq-aq and cr-cr field boundaries.
+
+ }
+ }
+
+ \subsection{OTHER CHANGES}{
+ \itemize{
+
\item In the return value of \code{diagram()}, \code{predominant.values}
previously only contained affinities (extracted from the \code{values}
element of the \code{eout} argument). Now it has the values for the
- predominant species extracted from the actual values used to make the plot
- (the \code{plotvals} list in the output), i.e. affinities divided by the
- balancing coefficients if \code{eout} is the output of \code{affinity}, or
- activities if \code{eout} is the output of \code{equilibrate}. Therefore it
- can now be used to draw contours or a color image showing the activities of
- the predominant species.
+ predominant species extracted from the actual values used to make the
+ plot (the \code{plotvals} list in the output), i.e. affinities divided
+ by the balancing coefficients if \code{eout} is the output of
+ \code{affinity}, or activities if \code{eout} is the output of
+ \code{equilibrate}. Therefore it can now be used to draw contours or a
+ color image showing the activities of the predominant species. This is
+ used for a diagram in a new preprint
+ (\href{https://doi.org/10.1101/2021.01.29.428804}{Dick, 2021}).
- \item Use standard abbreviations for hematite (Hem) and magnetite (Mag) in
- OBIGT/Berman_cr.csv.
-
\item Revise multi-metal.Rmd: Improve mineral abbreviations and placement
of labels; use updated DFT energies from Materials Project; add Mosaic
Stacking 2 (minerals and aqueous species); add \emph{K}\S{eff}
@@ -36,15 +50,13 @@
calculation (\href{https://doi.org/10.1016/j.gca.2021.01.038}{Robinson et
al., 2021}); add Δ\emph{G}\s{pbx} color scale.
- \item Restore \samp{lty.aq} and \samp{lty.cr} arguments to
- \code{diagram()} to control plotting of aq-aq and cr-cr field boundaries.
+ \item Use standard abbreviations for hematite (Hem) and magnetite (Mag) in
+ OBIGT/Berman_cr.csv.
- \item Add \samp{demo/Pourbaix.R} (Eh-pH diagrams with equisolubility
- lines, after \href{https://www.worldcat.org/oclc/563921897}{Pourbaix,
- 1974}).
-
\item Rename \samp{demo/berman.R} to \samp{potassium.R}.
+ \item Remove unused limSolve package from Suggests in DESCRIPTION.
+
}
}
Modified: pkg/CHNOSZ/inst/TODO
===================================================================
--- pkg/CHNOSZ/inst/TODO 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/inst/TODO 2021-03-03 09:57:48 UTC (rev 647)
@@ -83,10 +83,13 @@
- Fix legend in demo/mosaic.R (-4 and -6 contours are swapped)
-[20210227]
+[20210301]
-- Add demo: equisolubility lines (p. 312 of Pourbaix, 1974)
+- anintro.Rmd: change red to blue for sulfur basis species
-[20210301]
+[20210303]
-- anintro.Rmd: change red to blue for sulfur basis species
+- Add demo for effective equilibrium constant for diglycine on T-pH plot at
+ high pressure (DEW model; Robinson et al., 2021)
+
+- Combine solubility() and solubilities() into single function (for CHNOSZ 1.5.0).
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/man/diagram.Rd 2021-03-03 09:57:48 UTC (rev 647)
@@ -93,7 +93,7 @@
\details{
This function displays diagrams representing either chemical affinities, or equilibrium chemical activities of species.
-The first argument is the output from \code{\link{affinity}}, \code{\link{equilibrate}}, or \code{\link{solubility}}.
+The first argument is the output from \code{\link{affinity}}, \code{\link{equilibrate}}, \code{\link{solubility}}, or \code{\link{solubilities}}.
0-D diagrams, at a single point, are shown as \code{\link{barplot}}s.
1-D diagrams, for a single variable on the x-axis, are plotted as lines.
2-D diagrams, for two variables, are plotted as predominance fields.
@@ -128,7 +128,10 @@
This paragraph describes the effect of the \code{type} argument when the output from \code{solubility} is being used.
For \code{type} set to \samp{auto}, the equilibrium activities (i.e. computed solubilities) of each species are plotted.
If \code{type} is \samp{loga.balance}, the activity of the balancing basis species (i.e. total solubility) is plotted; this is represented by contours on a 2-D diagram.
-For examples that use these features, see \code{\link{solubility}} and various \code{\link{demos}}: \samp{DEW}, \samp{contour}, \samp{gold}, \samp{solubility}, \samp{sphalerite}, \samp{zinc}.
+For examples that use these features, see \code{\link{solubility}} and various \code{\link{demos}}: \samp{DEW}, \samp{contour}, \samp{gold}, \samp{solubility}, \samp{sphalerite}.
+
+When the output from \code{\link{solubilities}} is provided, the only supported option is to plot \samp{loga.balance}.
+See \code{demo("Pourbaix")} and \code{demo("zinc")}.
}
\section{1-D diagrams}{
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/man/examples.Rd 2021-03-03 09:57:48 UTC (rev 647)
@@ -17,7 +17,7 @@
"density", "ORP", "findit", "ionize", "buffer", "protbuff",
"glycinate", "mosaic", "copper", "arsenic", "solubility", "gold",
"contour", "sphalerite", "zinc", "Shh", "saturation",
- "adenine", "DEW", "lambda", "potsassium", "TCA", "aluminum", "AkDi",
+ "adenine", "DEW", "lambda", "potassium", "TCA", "aluminum", "AkDi",
"comproportionation", "Pourbaix"),
save.png=FALSE)
}
Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd 2021-03-02 02:36:42 UTC (rev 646)
+++ pkg/CHNOSZ/man/solubility.Rd 2021-03-03 09:57:48 UTC (rev 647)
@@ -1,14 +1,16 @@
\encoding{UTF-8}
\name{solubility}
\alias{solubility}
+\alias{solubilities}
\title{Equilibrium Chemical Activities of Species}
\description{
-Calculate chemical activities of species in equilibrium with a soluble basis species.
+Calculate chemical activities of aqueous species in equilibrium with a mineral or gas.
}
\usage{
solubility(aout, dissociation = NULL, find.IS = FALSE, in.terms.of = NULL,
codeanal = FALSE)
+ solubilities(a_cr, i_aq, in.terms.of = NULL, dissociation = NULL)
}
\arguments{
@@ -17,18 +19,46 @@
\item{find.IS}{logical, find the equilibrium ionic strength by iteration?}
\item{in.terms.of}{character, express the total solubility in terms of moles of this species}
\item{codeanal}{logical, print diagnostic messages and values of internal variables for code analysis}
+ \item{a_cr}{list, output from \code{\link{affinity}} or \code{\link{mosaic}} for one or more minerals}
+ \item{i_aq}{character or numeric, aqueous species that can be formed by dissolving the minerals}
}
\details{
-This function performs a simple task: from the values of \code{\link{affinity}} of formation reactions of species at given activity, it works backward to find the activities of species that make the affinities zero.
+\code{solubility} performs a simple task: given the values of \code{\link{affinity}} of formation reactions of species at specified activity, it works backward to find the activities of species that make the affinities zero.
This corresponds to complete equilibrium with all of the basis species.
Usually, the basis species should be set up so that the first basis species represents the substance being dissolved (a mineral such as CaCO\s3 or gas such as CO\s2).
Internally, this is treated as the conserved basis species, so it must be present in all of the formation reactions of the species.
-It is also possible to set the conserved basis species as other than the first one (see \code{demo(gold)}), but this implies that dissociation reactions are not occurring (see below).
-The \code{\link{species}} should be defined to represent one set of ions (anions or cations or their complexes) formed in solution, all involving the conserved basis species.
+To generate \code{aout}, the \code{\link{species}} should be defined to represent one set of aqueous species (including ions and/or neutral complexes) formed in solution, all involving the conserved basis species.
For a dissociation reaction, the second basis species should be used to represent the counterion (cation or anion).
+Other variables (pH, ionic strength, activities of other basis species) should be defined in the call to \code{\link{affinity}} to make \code{aout}.
+A more complex calculation involves setting \code{find.IS} to TRUE in order to determine the final ionic strength of dissolving a substance in pure water.
+This works by calculating the ionic strength from the equilibrium solubility calculation, then re-running \code{affinity} with those values.
+Note that for dissociation reactions, the ionic strength is calculated from both the ions present in the species definition and the counter ion, which should be the second basis species.
+The calculation is iterated until the ionic strength deviation at every point is lower than a preset tolerance (1e-4).
+Alternatively, speciation of counterions (e.g. ionized forms of carbonate or sulfate) can also be accomplished by using the \code{\link{mosaic}} function instead of \code{affinity}; this is used in the second example below.
+
+The output of \code{solubility} has the same format as that of \code{equilibrate}, and can be used by \code{\link{diagram}} with \code{type = "loga.balance"} to plot the solubilities, or with \code{type = NULL} to plot the activities of species.
+The value of \code{loga.balance} reflects the activity (or molality) of the conserved basis species, i.e. the thing being dissolved.
+Use \code{in.terms.of} to express this value in terms of another species.
+For example, the solubility of corundum (Al\s{2}O\s{3}) can be expressed in terms of the moles of Al\S{+3} in solution (see the vignette \viglink{anintro}).
+
+Whereas \code{solubility} is limited to calculations for a single mineral at a time, \code{solubilities} combines the results of solubility calculations for multiple minerals.
+\code{a_cr} should contain the affinities of only minerals with any set of basis species, as long as the first one contains the main element.
+\code{i_aq} gives the names or species indices of all the aqueous species that all have that element in combination with the other basis species (possibly generated using \code{\link{retrieve}}).
+The function calculates the solubilities for each of the minerals present in code{a_cr}.
+Then, it takes the \emph{minimum} solubility among all the minerals at each grid location.
+This corresponds to the maximum activity of dissolved species.
+At this point the solution is saturated (not oversaturated) with respect to the most stable mineral.
+
+The result of \code{solubilities} is stored in the \code{loga.balance} component of the returned list object.
+Only this component is used by \code{\link{diagram}} for plotting; the other components, while present, refer to the solubility calculations for only the last mineral.
+This function is used for generating isosolubility (aka. equisolubility) lines in \code{demo("Pourbaix")} and \code{demo{"zinc"}}.
+The latter demo combines the calculation of solubilities with a \code{\link{mosaic}} calculation to account for the speciation of aqueous sulfur.
+}
+
+\section{Dissociation Reactions}{
The function perfoms some additional steps to calculate the solubility of something that dissociates (not just dissolves).
For example, the dissolution of calcite (CaCO\s3), involves the release of both calcium ions and different forms of carbonate in solution, depending on the pH.
The equilibrium calculation must take account of the \emph{total} activity of the shared ion (Ca\S{+2}), which is unknown at the start of the calculation.
@@ -40,23 +70,10 @@
An example reaction of this type can be found in demo(solubility): CaCO3 (first basis species) = Ca+2 (second basis species) + CO3-2 (first species).
Note that if the conserved basis species is not the first basis species, then the automatic detection of \code{dissociation} will always return FALSE.
Therefore, a reaction corresponding to Au (fourth basis species) + ... = ... gives \code{dissociation} = FALSE (see \code{demo(gold)}).
-This algorithm for determining whether dissociation occurs is prone to error, so \code{dissociation} can be explicitly set in the arguments.
+This algorithm for determining whether dissociation occurs is not perfect, so \code{dissociation} can be explicitly set in the arguments.
A \emph{not recommended} alternative is to set \code{dissociation} to a numeric value corresponding to the stoichiometry of released species (i.e. 2 for a 1:1 electrolyte).
This setting indicates to calculate activities on a per-reaction basis, where each reaction has its own (independent) activity of Ca\S{+2}.
That does not give a complete equilibrium in the system, but may be required to reproduce some published diagrams.
-
-Note that other variables (pH, ionic strength, activities of other basis species) should be defined in the preceding call to \code{\link{affinity}}.
-However, for dissolving a substance in pure water, \code{find.IS} can be set to TRUE to determine the final ionic strength.
-This works by calculating the ionic strength from the equilibrium solubility calculation, then re-running \code{affinity} with those values.
-Note that for dissociation reactions, the ionic strength is calculated from both the ions present in the species definition and the counter ion, which should be the second basis species.
-The calculation is iterated until the ionic strength deviation at every point is lower than a preset tolerance (1e-4).
-Alternatively, speciation of counterions (e.g. ionized forms of carbonate or sulfate) can also be accomplished by using the \code{\link{mosaic}} function instead of \code{affinity}.
-See the second example for this method.
-
-The output of \code{solubility} has the same format as that of \code{equilibrate}, and can be used by \code{\link{diagram}} with \code{type = "loga.balance"} to plot the solubilities, or with \code{type = NULL} to plot the activities of species.
-The value of \code{loga.balance} reflects the activity (or molality) of the conserved basis species, i.e. the thing being dissolved.
-Use \code{in.terms.of} to express this value in terms of another species.
-For example, the solubility of corundum (Al\s{2}O\s{3}) can be expressed in terms of the moles of Al\S{+3} in solution (see the vignette \viglink{anintro}).
}
\section{Warning}{
@@ -65,31 +82,34 @@
}
\seealso{
-\code{demo("solubility")} adds \T-pH diagrams to the CO\s{2} and calcite example here.
+\code{demo("solubility")} shows \T-pH diagrams related to the CO\s{2} and calcite example here.
\code{demo("gold")} shows solubility calculations for Au in aqueous solutions with hydroxide, chloride, and hydrosulfide complexes.
-\code{\link{equilibrate}} calculates equilibrium chemical activities of species given a constant value of \code{loga.balance} (the logarithm of total activity of the conserved basis species).
+In the latter demo, the conserved basis species is not the first one (this might have been done to prevent the incorrect automatic detection of dissociation reactions).
+
+\code{\link{retrieve}} provides a way to list all of the aqueous species in the database that have the specified elements.
+In contrast to \code{solubility}, \code{\link{equilibrate}} calculates equilibrium chemical activities of species given a constant value of \code{loga.balance} (the logarithm of total activity of the conserved basis species).
}
\examples{\dontshow{reset()}
-## solubility of CO2 and calcite as a function of pH
+## Solubility of CO2 and calcite as a function of pH
opar <- par(mfrow = c(1, 2))
-## set pH range and resolution, constant temperature and ionic strength
+## Set pH range and resolution, constant temperature and ionic strength
pH <- c(0, 14)
res <- 100
T <- 25
IS <- 0
-## start with CO2 (not a dissociation reaction)
+## Start with CO2 (not a dissociation reaction)
basis(c("carbon dioxide", "H2O", "O2", "H+"))
-# ca. atmospheric PCO2
+# This is close to atmospheric PCO2
basis("CO2", -3.5)
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T, IS = IS)
s <- solubility(a)
-# first plot total activity line
+# First plot total activity line
diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
-# add activities of species
+# Add activities of species
diagram(s, ylim=c(-10, 4), add = TRUE, dy = 1)
-# add legend
+# Add legend
lexpr <- as.expression(c("total", expr.species("CO2", state = "aq"),
expr.species("HCO3-"), expr.species("CO3-2")))
legend("topleft", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
@@ -97,7 +117,7 @@
title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
list(what = expr.species("CO2"), T = T)), line = 1.5)
mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
-## now do calcite (a dissociation reaction)
+## Now do calcite (a dissociation reaction)
basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T, IS = IS)
@@ -111,33 +131,33 @@
mtext("cf. Fig. 4A of Manning et al., 2013")
par(opar)
-## two ways to calculate pH-dependent solubility of calcite
+## Two ways to calculate pH-dependent solubility of calcite
## with ionic strength determination
-## method 1: CO2 and carbonate species as formed species
+## Method 1: CO2 and carbonate species as formed species
basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
species(c("CO2", "HCO3-", "CO3-2"))
-# ionic strength calculations don't converge below around pH=3
+# Ionic strength calculations don't converge below around pH = 3
a <- affinity(pH = c(4, 14))
sa0 <- solubility(a)
saI <- solubility(a, find.IS = TRUE)
-## method 2: CO2 and carbonate species as basis species
+## Method 2: CO2 and carbonate species as basis species
basis(c("calcite", "CO2", "H2O", "O2", "H+"))
species(c("Ca+2"))
m <- mosaic(c("CO2", "HCO3-", "CO3-2"), pH = c(4, 14))
sm0 <- solubility(m)
smI <- solubility(m, find.IS = TRUE)
-## plot the results
+## Plot the results
plot(0, 0, xlab="pH", ylab="solubility, log mol", xlim = c(4, 14), ylim = c(-5, 2))
-# method 1 with/without ionic strength
+# Method 1 with/without ionic strength
lines(a$vals[[1]], saI$loga.balance, lwd=5, col="lightblue")
lines(a$vals[[1]], sa0$loga.balance, lwd=5, col="pink")
-# method 2 with/without ionic strength
+# Method 2 with/without ionic strength
lines(a$vals[[1]], smI$loga.balance, lty=2)
lines(a$vals[[1]], sm0$loga.balance, lty=2)
legend("topright", c("I = 0", "I = calculated", "mosaic method"),
col = c("pink", "lightblue", "black"), lwd = c(5, 5, 1), lty = c(1, 1, 2))
title(main = "Solubility of calcite: Ionic strength and mosaic method")
-# the two methods give nearly equivalent results
+# The two methods give nearly equivalent results
stopifnot(all.equal(sa0$loga.balance, sm0$loga.balance))
stopifnot(all.equal(saI$loga.balance, smI$loga.balance, tolerance = 0.003))
## NOTE: the second method (using mosaic) takes longer, but is
More information about the CHNOSZ-commits
mailing list