[CHNOSZ-commits] r340 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 5 08:08:13 CET 2018
Author: jedick
Date: 2018-11-05 08:08:12 +0100 (Mon, 05 Nov 2018)
New Revision: 340
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/NaCl.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/tests/testthat/test-logmolality.R
pkg/CHNOSZ/tests/testthat/test-nonideal.R
Log:
nonideal(): use ion-size parameters for different species
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-04 14:44:13 UTC (rev 339)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-05 07:08:12 UTC (rev 340)
@@ -1,6 +1,6 @@
-Date: 2018-11-04
+Date: 2018-11-05
Package: CHNOSZ
-Version: 1.1.3-47
+Version: 1.1.3-48
Title: Thermodynamic Calculations and Diagrams for Geo(bio)chemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2018-11-04 14:44:13 UTC (rev 339)
+++ pkg/CHNOSZ/R/nonideal.R 2018-11-05 07:08:12 UTC (rev 340)
@@ -54,9 +54,9 @@
}
# function for Debye-Huckel equation with B-dot extended term parameter (Helgeson, 1969)
- Helgeson <- function(Z, I, T, P, A_DH, B_DH, prop = "loggamma") {
- # "distance of closest approach" of ions in NaCl solutions (HKF81 Table 2)
- acirc <- 3.72e-8 # cm
+ Helgeson <- function(Z, I, T, P, A_DH, B_DH, prop = "loggamma", acirc) {
+ ## "distance of closest approach" of ions in NaCl solutions (HKF81 Table 2)
+ #acirc <- 3.72e-8 # cm
if(method=="Helgeson") loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) + Bdot * I
else if(method=="Helgeson0") loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5)
R <- 1.9872 # gas constant, cal K^-1 mol^-1
@@ -70,39 +70,69 @@
# get species indices
if(!is.numeric(species[[1]])) species <- info(species, "aq")
+ # loop over species #1: get the charge
+ Z <- numeric(length(species))
+ for(i in 1:length(species)) {
+ # force a charge count even if it's zero
+ mkp <- makeup(c("Z0", species[i]), sum=TRUE)
+ thisZ <- mkp[match("Z", names(mkp))]
+ # don't do anything for neutral species (Z absent from formula or equal to zero)
+ if(is.na(thisZ)) next
+ if(thisZ==0) next
+ Z[i] <- thisZ
+ }
+ # get species formulas to assign acirc 20181105
+ if(grepl("Helgeson", method)) {
+ formula <- get("thermo")$obigt$formula[species]
+ # "ion size paramter" taken from UT_SIZES.REF of HCh package (Shvarov and Bastrakov, 1999),
+ # based on Table 2.7 of Garrels and Christ, 1965
+ acircdat <- c("Rb+"=2.5, "Cs+"=2.5, "NH4+"=2.5, "Tl+"=2.5, "Ag+"=2.5,
+ "K+"=3, "Cl-"=3, "Br-"=3, "I-"=3, "NO3-"=3,
+ "OH-"=3.5, "F-"=3.5, "HS-"=3.5, "BrO3-"=3.5, "IO3-"=3.5, "MnO4-"=3.5,
+ "Na+"=4, "HCO3-"=4, "H2PO4-"=4, "HSO3-"=4, "Hg2+2"=4, "SO4-2"=4, "SeO4-2"=4, "CrO4-2"=4, "HPO4-2"=4, "PO4-3"=4,
+ "Pb+2"=4.5, "CO3-2"=4.5, "SO4-2"=4.5, "MoO4-2"=4.5,
+ "Sr+2"=5, "Ba+2"=5, "Ra+2"=5, "Cd+2"=5, "Hg+2"=5, "S-2"=5, "WO4-2"=5,
+ "Li+"=6, "Ca+2"=6, "Cu+2"=6, "Zn+2"=6, "Sn+2"=6, "Mn+2"=6, "Fe+2"=6, "Ni+2"=6, "Co+2"=6,
+ "Mg+2"=8, "Be+2"=8,
+ "H+"=9, "Al+3"=9, "Cr+3"=9, "La+3"=9, "Ce+3"=9, "Y+3"=9, "Eu+3"=9,
+ "Th+4"=11, "Zr+4"=11, "Ce+4"=11, "Sn+4"=11)
+ acirc <- as.numeric(acircdat[formula])
+ acirc[is.na(acirc)] <- 4.5
+ # make a message
+ nZ <- sum(Z!=0)
+ if(nZ > 1) message("nonideal: using ", paste(acirc[Z!=0], collapse=" "), " for ion size parameters of ", paste(formula[Z!=0], collapse=" "))
+ else if(nZ==1) message("nonideal: using ", acirc[Z!=0], " for ion size parameter of ", formula[Z!=0])
+ # use correct units for ion size paramter
+ acirc <- acirc * 10^-8
+ }
+ # loop over species #2: activity coefficient calculations
iH <- info("H+")
ie <- info("e-")
speciesprops <- as.list(speciesprops)
ndid <- 0
- # loop over species
for(i in 1:length(species)) {
myprops <- speciesprops[[i]]
- # get the charge from the chemical formula
- # force a charge count even if it's zero
- mkp <- makeup(c("Z0", species[i]), sum=TRUE)
- Z <- mkp[match("Z", names(mkp))]
- # don't do anything for neutral species (Z absent from formula or equal to zero)
- if(is.na(Z)) next
- if(Z==0) next
# to keep unit activity coefficients of the proton and electron
if(species[i] == iH & get("thermo")$opt$ideal.H) next
if(species[i] == ie & get("thermo")$opt$ideal.e) next
+ # skip neutral species
+ if(Z[i]==0) next
didit <- FALSE
for(j in 1:ncol(myprops)) {
pname <- colnames(myprops)[j]
if(method=="Alberty" & pname %in% c("G", "H", "S", "Cp")) {
- myprops[, j] <- myprops[, j] + Alberty(Z, IS, T, pname)
+ myprops[, j] <- myprops[, j] + Alberty(Z[i], IS, T, pname)
didit <- TRUE
} else if(grepl("Helgeson", method) & pname %in% c("G", "H", "S", "Cp")) {
- myprops[, j] <- myprops[, j] + Helgeson(Z, IS, T, P, A_DH, B_DH, pname)
+ myprops[, j] <- myprops[, j] + Helgeson(Z[i], IS, T, P, A_DH, B_DH, pname, acirc[i])
didit <- TRUE
}
}
# append a loggam column if we did any nonideal calculations of thermodynamic properties
if(didit) {
- if(method=="Alberty") myprops <- cbind(myprops, loggam = Alberty(Z, IS, T))
+ if(method=="Alberty") myprops <- cbind(myprops, loggam = Alberty(Z[i], IS, T))
else if(grepl("Helgeson", method)) {
- myprops <- cbind(myprops, loggam = Helgeson(Z, IS, T, P, A_DH, B_DH))
+ myprops <- cbind(myprops, loggam = Helgeson(Z[i], IS, T, P, A_DH, B_DH, "loggamma", acirc[i]))
}
}
speciesprops[[i]] <- myprops
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-04 14:44:13 UTC (rev 339)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-05 07:08:12 UTC (rev 340)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.3-47 (2018-11-01)
+CHANGES IN CHNOSZ 1.1.3-48 (2018-11-05)
---------------------------------------
NEW FEATURES
@@ -10,6 +10,10 @@
of NaCl in water, taking account of activity coefficients and the
reaction Na+ + Cl- = NaCl(aq).
+- nonideal() now uses specific ion-size parameters for different ions,
+ in accord with the HCh package (Shvarov and Bastrakov, 1999).
+ Parameters are from Table 2.7 of Garrels and Christ, 1965.
+
- Add demo/gold.R for calculations of Au solubility in hydrothermal
chloride and sulfide solutions (based on diagrams from Akinfiev and
Zotov, 2001, Stefánsson and Seward, 2004, and Williams-Jones et al.,
Modified: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd 2018-11-04 14:44:13 UTC (rev 339)
+++ pkg/CHNOSZ/man/NaCl.Rd 2018-11-05 07:08:12 UTC (rev 340)
@@ -42,7 +42,7 @@
}
\examples{\dontshow{data(thermo)}
-# ionic strength and activity coefficient of Cl-
+# ionic strength of solution and activity coefficient of Cl-
# from HCh (Shvarov and Bastrakov, 1999) at 1000 bar,
# 100, 200, and 300 degress C, and 1 to 6 molal NaCl
m.HCh <- 1:6
@@ -60,7 +60,7 @@
gam.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
# NaCl() is *not* vectorized over m_tot, so we use a loop here
for(i in 1:length(m_tot)) {
- NaCl.out <- NaCl(c(100, 300, 500), m_tot=m_tot[i])
+ NaCl.out <- NaCl(c(100, 300, 500), P=1000, m_tot=m_tot[i])
IS.calc[i, ] <- NaCl.out$IS
gam.calc[i, ] <- NaCl.out$gamma
}
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-11-04 14:44:13 UTC (rev 339)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-11-05 07:08:12 UTC (rev 340)
@@ -39,7 +39,7 @@
If \code{method} is \samp{Helgeson}, the \dQuote{B-dot} equation is used.
This equation seems to have been originally proposed by Huckel, 1925; parameters were derived for use at high temperature and pressure by Helgeson, 1969; Helgeson et al., 1981; Manning, 2013.
-The distance of closest approach (the \dQuote{ion size parameter}) is set to 3.72 Angstrom, which is appropriate for NaCl-dominated solutions (Helgeson et al., 1981 Table 2).
+The distance of closest approach (the \dQuote{ion size parameter}) is taken from the UT_SIZES.REF file of the HCh package (Shvarov and Bastrakov, 1992), which is based on Table 2.7 of Garrels and Christ, 1965.
In addition to \code{IS} and \code{T}, this method depends on values of \code{P}, \code{A_DH}, and \code{B_DH} given in the arguments.
The calculation of \dQuote{B-dot}, also used in the equations, is made within \code{nonideal} by calling the \code{Bdot} function.
For some uses, it is desirable to set the \dQuote{B-dot} parameter to zero; this can be done by setting the method to \code{Helgeson0}.
@@ -138,7 +138,7 @@
nonideal(oldnon) # same as nonideal("Helgeson")
## activity coefficients for monovalent ions at 700 degC, 10 kbar
-# after Manning, 2010, Fig. 7
+# after Manning, 2013, Fig. 7
IS <- c(0.001, 0.01, 0.1, 1, 2, 2.79)
# we're above 5000 bar, so need to use IAPWS-95 or DEW
oldwat <- water("DEW")
@@ -153,10 +153,9 @@
# the activity coefficient; check values of the latter
Manning_gamma <- c(0.93, 0.82, 0.65, 0.76, 1.28, 2)
gamma <- 10^nonidealprops[[1]]$loggam
-# we're getting progressively further from his values with
-# higher IS; not sure why
+# the error is larger at higher IS
stopifnot(maxdiff(gamma[1], Manning_gamma[1]) < 0.01)
-stopifnot(maxdiff(gamma, Manning_gamma) < 0.23)
+stopifnot(maxdiff(gamma, Manning_gamma) < 0.36)
## data and splines used for calculating B-dot
## (extended term parameter)
@@ -167,6 +166,8 @@
\references{
Alberty, R. A. (2003) \emph{Thermodynamics of Biochemical Reactions}, John Wiley & Sons, Hoboken, New Jersey, 397 p. \url{http://www.worldcat.org/oclc/51242181}
+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}
+
Helgeson, H. C. (1969) Thermodynamics of hydrothermal systems at elevated temperatures and pressures. \emph{Am. J. Sci.} \bold{267}, 729--804. \url{https://doi.org/10.2475/ajs.267.7.729}
Helgeson, H. C., Kirkham, D. H. and Flowers, G. C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600\degC and 5 Kb. \emph{Am. J. Sci.} \bold{281}, 1249--1516. \url{https://doi.org/10.2475/ajs.281.10.1249}
@@ -176,6 +177,8 @@
Manning, C. E. (2013) Thermodynamic modeling of fluid-rock interaction at mid-crustal to upper-mantle conditions. \emph{Rev. Mineral. Geochem.} \bold{76}, 135--164. \url{https://doi.org/10.2138/rmg.2013.76.5}
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{https://doi.org/10.2138/rmg.2013.75.5}
+
+Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. \emph{Australian Geological Survey Organisation} \bold{1999/25}.
}
\concept{Thermodynamic calculations}
Modified: pkg/CHNOSZ/tests/testthat/test-logmolality.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-04 14:44:13 UTC (rev 339)
+++ pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-05 07:08:12 UTC (rev 340)
@@ -10,10 +10,10 @@
wprop <- water(c("A_DH", "B_DH"), P=1)
nonid <- nonideal(c("H+", "HCO3-"), subcrt(c("H+", "HCO3-"), T=25)$out, IS=1, T=298.15, P=1, A_DH=wprop$A_DH, B_DH=wprop$B_DH)
# here we have a small difference due to rounding of the expected value:
- expect_maxdiff(nonid[[2]]$loggam, -0.1890084, 1e-7)
+ expect_maxdiff(nonid[[2]]$loggam, -0.1798625, 1e-7)
# the short way...
loggam <- subcrt(c("H+", "HCO3-"), T=25, IS=1)$out[[2]]$loggam
- # we get that loggam(H+)=0 and loggam(HCO3-)=-0.189
+ # we get that loggam(H+)=0 and loggam(HCO3-)=-0.180
expect_equal(nonid[[2]]$loggam, loggam)
## take-home message -1: with default settings, the activity coefficient of H+ is always 1
Modified: pkg/CHNOSZ/tests/testthat/test-nonideal.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-nonideal.R 2018-11-04 14:44:13 UTC (rev 339)
+++ pkg/CHNOSZ/tests/testthat/test-nonideal.R 2018-11-05 07:08:12 UTC (rev 340)
@@ -75,3 +75,35 @@
expect_type(subcrt(c("Zn+2", "Cl-", "ZnCl+"), c(-1, -1, 1), T=200, P=16, IS=0.05), "list")
})
+
+# 20181105
+test_that("activity coefficients are similar to those from HCh", {
+ # ionic strength of solution and activity coefficients of Na+ and Cl-
+ # from HCh (Shvarov and Bastrakov, 1999) at 1000 bar,
+ # 100, 200, and 300 degress C, and 1 to 6 molal NaCl
+ # using the default "B-dot" activity coefficient model (Helgeson, 1969)
+ IS.HCh <- list(`100`=c(0.992, 1.969, 2.926, 3.858, 4.758, 5.619),
+ `300`=c(0.807, 1.499, 2.136, 2.739, 3.317, 3.875),
+ `500`=c(0.311, 0.590, 0.861, 1.125, 1.385, 1.642))
+ gamCl.HCh <- list(`100`=c(0.565, 0.545, 0.551, 0.567, 0.589, 0.615),
+ `300`=c(0.366, 0.307, 0.275, 0.254, 0.238, 0.224),
+ `500`=c(0.19, 0.137, 0.111, 0.096, 0.085, 0.077))
+ gamNa.HCh <- list(`100`=c(0.620, 0.616, 0.635, 0.662, 0.695, 0.730),
+ `300`=c(0.421, 0.368, 0.339, 0.318, 0.302, 0.288),
+ `500`=c(0.233, 0.180, 0.155, 0.138, 0.126, 0.117))
+ # calculate activity coefficent of Cl- at each temperature
+ gamCl.100 <- 10^subcrt("Cl-", T=100, P=1000, IS=IS.HCh$`100`)$out$`Cl-`$loggam
+ gamCl.300 <- 10^subcrt("Cl-", T=300, P=1000, IS=IS.HCh$`300`)$out$`Cl-`$loggam
+ gamCl.500 <- 10^subcrt("Cl-", T=500, P=1000, IS=IS.HCh$`500`)$out$`Cl-`$loggam
+ # TODO: get lower differences by adjusting the activity coefficient model in CHNOSZ
+ expect_maxdiff(gamCl.100, gamCl.HCh$`100`, 0.73)
+ expect_maxdiff(gamCl.300, gamCl.HCh$`300`, 0.22)
+ expect_maxdiff(gamCl.500, gamCl.HCh$`500`, 0.04)
+ # calculate activity coefficent of Cl- at each temperature
+ gamNa.100 <- 10^subcrt("Na+", T=100, P=1000, IS=IS.HCh$`100`)$out$`Na+`$loggam
+ gamNa.300 <- 10^subcrt("Na+", T=300, P=1000, IS=IS.HCh$`300`)$out$`Na+`$loggam
+ gamNa.500 <- 10^subcrt("Na+", T=500, P=1000, IS=IS.HCh$`500`)$out$`Na+`$loggam
+ expect_maxdiff(gamNa.100, gamNa.HCh$`100`, 0.67)
+ expect_maxdiff(gamNa.300, gamNa.HCh$`300`, 0.18)
+ expect_maxdiff(gamNa.500, gamNa.HCh$`500`, 0.06)
+})
More information about the CHNOSZ-commits
mailing list