[CHNOSZ-commits] r79 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 6 12:12:22 CET 2015
Author: jedick
Date: 2015-03-06 12:12:22 +0100 (Fri, 06 Mar 2015)
New Revision: 79
Added:
pkg/CHNOSZ/demo/solubility.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/examples.Rd
Log:
new demo: solubility.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-03-06 06:52:51 UTC (rev 78)
+++ pkg/CHNOSZ/DESCRIPTION 2015-03-06 11:12:22 UTC (rev 79)
@@ -1,6 +1,6 @@
Date: 2015-03-06
Package: CHNOSZ
-Version: 1.0.3-16
+Version: 1.0.3-17
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2015-03-06 06:52:51 UTC (rev 78)
+++ pkg/CHNOSZ/R/examples.R 2015-03-06 11:12:22 UTC (rev 79)
@@ -31,7 +31,7 @@
demos <- function(which=c("sources", "NaCl", "cordierite",
"phosphate", "nucleobase", "orp",
- "findit", "CO2Ac", "nonideal", "TPX", "mosaic")) {
+ "findit", "CO2Ac", "nonideal", "TPX", "mosaic", "solubility")) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
# say something so the user sees where we are
Modified: pkg/CHNOSZ/R/util.expression.R
===================================================================
--- pkg/CHNOSZ/R/util.expression.R 2015-03-06 06:52:51 UTC (rev 78)
+++ pkg/CHNOSZ/R/util.expression.R 2015-03-06 11:12:22 UTC (rev 79)
@@ -36,7 +36,7 @@
# write a designation of physical state
# use the state given in log if it's a gas or neutral aqueous species
if(log %in% c("g", "gas")) state <- "g"
- else if(!"Z" %in% names(elements)) state <- log
+ else if(!"Z" %in% names(elements) & !missing(log)) state <- log
if(state != "") {
# subscript it if we're not in a log expression
if(log != "") expr <- substitute(a*group('(',italic(b),')'),list(a=expr, b=state))
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2015-03-06 06:52:51 UTC (rev 78)
+++ pkg/CHNOSZ/demo/00Index 2015-03-06 11:12:22 UTC (rev 79)
@@ -9,3 +9,4 @@
nonideal activity coefficient of charged species, using the IS argument of subcrt
TPX metastablilities of selected ionized bacterial thiol peroxidases as a chemical activity buffer
mosaic Eh-pH diagram for iron oxides, sulfides and carbonate with two sets of changing basis species
+solubility solubility of calicite or CO2(gas) as a function of pH
Added: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R (rev 0)
+++ pkg/CHNOSZ/demo/solubility.R 2015-03-06 11:12:22 UTC (rev 79)
@@ -0,0 +1,83 @@
+# demo showing how to calculate CO2(gas) or calcite solubility and aqueous carbonate speciation
+# the affinity() ... equilibrate() sequence in CHNOSZ gives *metastable equilibrium activities*
+# (activities for a total activity of the balanced component given in loga.balance) ...
+# here we are interested in finding the value of loga.balance itself
+# (the total activity of [CO2, HCO3-, CO3-2] species in the aqueous phase).
+# this total activity is the solubility of CO2(gas) or calcite if the affinities of the
+# aqueous species as formed from CO2(gas) or calcite are all equal to zero.
+# note that the affinities for species in metastable equilibrium are all equal.
+# Afun() calculates the metastable equilibrium affinities for a given loga.balance
+# and uniroot() finds the loga.balance where they are zero
+# additionally, if we are reacting calcite, the activity of Ca+2 should be set equal to loga.balance
+
+# for comparison with published calcite solubility plot, see Fig. 4A in
+# Manning et al., 2013, Reviews in Mineralogy & Geochemistry, v. 75, pp. 109-148
+# (doi: 10.2138/rmg.2013.75.5)
+
+# for comparison with published CO2 solubility plot, see Fig. 4.5 in
+# Stumm and Morgan, 1996, Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters
+# (New York: John Wiley & Sons), 3rd edition
+
+# set this to CO2 or calcite
+#what <- "CO2"
+what <- "calcite"
+
+# function to return the affinity of the metastable equilibrium species
+Afun <- function(loga.balance=-3, T=25) {
+ if(what=="calcite") basis("Ca+2", loga.balance)
+ a <- affinity(T=T)
+ e <- equilibrate(a, loga.balance=loga.balance)
+ # set metastable activities and re-calculate the affinity
+ species(1:3, unlist(e$loga.equil))
+ a <- affinity(T=T)
+ # check they're actually equal
+ stopifnot(all(abs(unlist(a$values) - a$values[[1]]) < 1e-10))
+ return(a$values[[1]])
+}
+
+# set up system
+if(what=="CO2") {
+ basis("CHNOS+")
+ basis("CO2", "gas")
+ # ca. atmospheric PCO2
+ basis("CO2", -3.5)
+} else if(what=="calcite") {
+ basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
+}
+species(c("CO2", "HCO3-", "CO3-2"))
+T <- 25
+# decrease this for higher resolution
+pHstep <- 1
+
+# where we'll store the results
+loga.tot <- numeric()
+loga.CO2 <- loga.HCO3 <- loga.CO3 <- numeric()
+
+# loop over pH range
+pHs <- seq(0, 14, pHstep)
+for(pH in pHs) {
+ print(paste("pH =", pH))
+ basis("pH", pH)
+ # this is for the solubility
+ loga.balance <- suppressMessages(uniroot(Afun, c(-10, 10), T=T)$root)
+ loga.tot <- c(loga.tot, loga.balance)
+ # this is for the speciation
+ if(what=="calcite") basis("Ca+2", loga.balance)
+ a <- affinity(T=T)
+ e <- equilibrate(a, loga.balance=loga.balance)
+ loga <- unlist(e$loga.equil)
+ loga.CO2 <- c(loga.CO2, loga[1])
+ loga.HCO3 <- c(loga.HCO3, loga[2])
+ loga.CO3 <- c(loga.CO3, loga[3])
+}
+
+# make plot
+ylim <- c(-10, 4)
+thermo.plot.new(xlim=range(pHs), ylim=ylim, xlab="pH", ylab="log a")
+lines(pHs, loga.tot, lwd=3)
+lines(pHs, loga.CO2, lwd=2)
+lines(pHs, loga.HCO3, lty=2, lwd=2)
+lines(pHs, loga.CO3, lty=3, lwd=2)
+legend(ifelse(what=="calcite", "topright", "topleft"), lty=c(1, 1:3), lwd=c(3, 2, 2, 2),
+ legend=as.expression(c("total", expr.species("CO2", state="aq"), expr.species("HCO3-"), expr.species("CO3-2"))))
+title(main=substitute(what~"solubility at"~T~degree*"C", list(what=what, T=T)))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-03-06 06:52:51 UTC (rev 78)
+++ pkg/CHNOSZ/inst/NEWS 2015-03-06 11:12:22 UTC (rev 79)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.3-16 (2015-03-06)
+CHANGES IN CHNOSZ 1.0.3-17 (2015-03-06)
---------------------------------------
- Add files with average amino acid compositions of proteins from Bison
@@ -50,6 +50,8 @@
- Add thermo$opt$varP option as a flag for subcrt() to calculate Gibbs
energies of gases using a variable-pressure standard state.
+- Add 'solubility.R' demo.
+
CHANGES IN CHNOSZ 1.0.3 (2014-01-12)
------------------------------------
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2015-03-06 06:52:51 UTC (rev 78)
+++ pkg/CHNOSZ/man/examples.Rd 2015-03-06 11:12:22 UTC (rev 79)
@@ -15,7 +15,7 @@
examples(do.png = FALSE)
demos(which = c("sources", "NaCl", "cordierite",
"phosphate", "nucleobase", "orp", "findit",
- "CO2Ac", "nonideal", "TPX", "mosaic"))
+ "CO2Ac", "nonideal", "TPX", "mosaic", "solubility"))
}
\details{
@@ -41,6 +41,7 @@
\code{nonideal} \tab activity coefficient of charged species (Alberty, 2003), using the \code{IS} argument of \code{subcrt} \cr
\code{TPX} \tab metastablilities of selected ionized bacterial thiol peroxidases as a chemical activity buffer \cr
\code{mosaic} \tab Eh-pH diagram for iron oxides, sulfides and carbonate with two sets of changing basis species (Garrels and Christ, 1965) \cr
+ \code{solubility} \tab solubility of calcite (cf. Manning et al., 2013 Fig. 4A) or CO2(gas) (cf. Stumm and Morgan, 1996 Fig 4.5) as a function of pH \cr
}
}
@@ -50,7 +51,11 @@
Garrels, R. M. and Christ, C. L. (1965) \emph{Solutions, Minerals, and Equilibria}, Harper & Row, New York, 450 p. \url{http://www.worldcat.org/oclc/517586}
+ Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. \emph{Rev. Mineral. Geochem.} \bold{75}, 109--148. \url{http://dx.doi.org/10.2138/rmg.2013.75.5}
+
Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 \eqn{^{\circ}}{°}C and 5 kbar. \emph{J. Chem. Soc. Faraday Trans.} \bold{88}, 803--826. \url{http://dx.doi.org/10.1039/FT9928800803}
+
+ Stumm, W. and Morgan, J. J. (1996) \emph{Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters}, John Wiley & Sons, New York, 1040 p. \url{http://www.worldcat.org/oclc/31754493}
}
More information about the CHNOSZ-commits
mailing list