[CHNOSZ-commits] r251 - in pkg/CHNOSZ: . R inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 11 17:40:05 CEST 2017
Author: jedick
Date: 2017-10-11 17:40:05 +0200 (Wed, 11 Oct 2017)
New Revision: 251
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/tests/testthat/test-nonideal.R
pkg/CHNOSZ/vignettes/obigt.Rmd
pkg/CHNOSZ/vignettes/obigt.bib
Log:
water(): add A_DH and B_DH Debye-Huckel coefficients
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-10 15:22:41 UTC (rev 250)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-11 15:40:05 UTC (rev 251)
@@ -1,6 +1,6 @@
Date: 2017-10-10
Package: CHNOSZ
-Version: 1.1.0-49
+Version: 1.1.0-50
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2017-10-10 15:22:41 UTC (rev 250)
+++ pkg/CHNOSZ/R/water.R 2017-10-11 15:40:05 UTC (rev 251)
@@ -20,10 +20,10 @@
if(length(P) < length(T)) P <- rep(P, length.out = length(T))
else if(length(T) < length(P)) T <- rep(T, length.out = length(P))
}
- # turn 273.15 K to 273.16 K (needed for water.SUPCRT92 at Psat)
- T[T == 273.15] <- 273.16
wopt <- get("thermo")$opt$water
if(grepl("SUPCRT", wopt)) {
+ # change 273.15 K to 273.16 K (needed for water.SUPCRT92 at Psat)
+ if(identical(P, "Psat")) T[T == 273.15] <- 273.16
# get properties using SUPCRT92
w.out <- water.SUPCRT92(property, T, P)
}
@@ -49,7 +49,8 @@
"Speed", "alpha", "beta", "epsilon", "visc",
"tcond", "surten", "tdiff", "Prndtl", "visck", "albe",
"ZBorn", "YBorn", "QBorn", "daldT", "XBorn",
- "V", "rho", "Psat", "E", "kT")
+ "V", "rho", "Psat", "E", "kT",
+ "A_DH", "B_DH")
if(is.null(property)) return(available_properties)
# check for availability of properties
iprop <- match(property, available_properties)
@@ -121,7 +122,12 @@
Psat=P.out
E <- V*w.out$alpha
kT <- V*w.out$beta
- w.out <- cbind(w.out, data.frame(V=V, rho=rho, Psat=Psat, E=E, kT=kT))
+ # A and B parameters in Debye-Huckel equation:
+ # Helgeson (1969) doi:10.2475/ajs.267.7.729
+ # Manning (2013) doi:10.2138/rmg.2013.76.5
+ A_DH <- 1.8246e6 * rho.out^0.5 / (w.out$epsilon * T)^1.5
+ B_DH <- 50.29e8 * rho.out^0.5 / (w.out$epsilon * T)^0.5
+ w.out <- cbind(w.out, data.frame(V=V, rho=rho, Psat=Psat, E=E, kT=kT, A_DH=A_DH, B_DH=B_DH))
}
# tell the user about any problems
if(any(err.out==1)) {
@@ -143,7 +149,8 @@
available_properties <- c("A", "G", "S", "U", "H", "Cv", "Cp",
"Speed", "epsilon",
"YBorn", "QBorn", "XBorn", "NBorn", "UBorn",
- "V", "rho", "Psat", "de.dT", "de.dP", "pressure")
+ "V", "rho", "Psat", "de.dT", "de.dP", "pressure",
+ "A_DH", "B_DH")
if(is.null(property)) return(available_properties)
# to get the properties of water via IAPWS-95
message(paste("water.IAPWS95: calculating", length(T), "values for"), appendLF=FALSE)
@@ -292,6 +299,13 @@
my.rho <- rho.IAPWS95(T, P, get("thermo")$opt$IAPWS.sat)
rho <- function() my.rho
}
+ # get epsilon and A_DH, B_DH (so we calculate epsilon only once)
+ if(any(property %in% c("epsilon", "A_DH", "B_DH"))) {
+ my.epsilon <- epsilon()
+ epsilon <- function() my.epsilon
+ A_DH <- function() 1.8246e6 * (my.rho/1000)^0.5 / (my.epsilon * T)^1.5
+ B_DH <- function() 50.29e8 * (my.rho/1000)^0.5 / (my.epsilon * T)^0.5
+ }
for(i in 1:length(property)) {
if(property[i] %in% c("E", "kT", "alpha", "daldT", "beta")) {
# E and kT aren't here yet... set them to NA
@@ -305,7 +319,7 @@
w.out[, i] <- wnew
}
message("")
- # use uppercase property names (including properties available in SUPCRT that might be NA here)
+ # include properties available in SUPCRT that might be NA here
wprop <- unique(c(water.SUPCRT92(), available_properties))
iprop <- match(property, wprop)
property[!is.na(iprop)] <- wprop[na.omit(iprop)]
@@ -315,7 +329,7 @@
# get water properties from DEW model for use by subcrt() 20170925
water.DEW <- function(property = NULL, T = 373.15, P = 1000) {
- available_properties <- c("G", "epsilon", "QBorn", "V", "rho", "beta")
+ available_properties <- c("G", "epsilon", "QBorn", "V", "rho", "beta", "A_DH", "B_DH")
if(is.null(property)) return(available_properties)
# we can't use Psat here
if(identical(P, "Psat")) stop("Psat isn't supported in this implementation of the DEW model. Try setting P to at least 1000 bar.")
@@ -331,16 +345,19 @@
out <- matrix(NA, ncol=length(property), nrow=ncond)
out <- as.data.frame(out)
colnames(out) <- property
- # calculate rho if it's needed for any other properties
- if(any(c("rho", "V", "QBorn", "epsilon", "beta") %in% property)) rho <- calculateDensity(pressure, temperature)
+ # calculate rho and epsilon if they're needed for any other properties
+ if(any(c("rho", "V", "QBorn", "epsilon", "beta", "A_DH", "B_DH") %in% property)) rho <- calculateDensity(pressure, temperature)
+ if(any(c("epsilon", "A_DH", "B_DH") %in% property)) epsilon <- calculateEpsilon(rho, temperature)
# fill in columns with values
if("rho" %in% property) out$rho <- rho*1000 # use kg/m^3 (like SUPCRT)
if("V" %in% property) out$V <- 18.01528/rho
if("G" %in% property) out$G <- calculateGibbsOfWater(pressure, temperature)
if("QBorn" %in% property) out$QBorn <- calculateQ(rho, temperature)
- if("epsilon" %in% property) out$epsilon <- calculateEpsilon(rho, temperature)
+ if("epsilon" %in% property) out$epsilon <- epsilon
# divide drhodP by rho to get units of bar^-1
if("beta" %in% property) out$beta <- calculate_drhodP(rho, temperature) / rho
+ if("A_DH" %in% property) out$A_DH <- 1.8246e6 * rho^0.5 / (epsilon * T)^1.5
+ if("B_DH" %in% property) out$B_DH <- 50.29e8 * rho^0.5 / (epsilon * T)^0.5
# use SUPCRT-calculated values below 100 degC and/or below 1000 bar
ilow <- T < 373.15 | P < 1000
if(any(ilow)) {
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-10 15:22:41 UTC (rev 250)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-11 15:40:05 UTC (rev 251)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-49 (2017-10-10)
+CHANGES IN CHNOSZ 1.1.0-50 (2017-10-11)
---------------------------------------
MAJOR CHANGES:
@@ -35,6 +35,10 @@
"percent carbon" plots for systems where the species have different
carbon numbers.
+- All three water options (SUPCRT92, IAPWS95, DEW) can be used to
+ calculate 'A_DH' and 'B_DH', i.e. A and B coefficients in the B-dot
+ (extended Debye-Huckel) equation of Helgeson, 1969.
+
- For minerals with phase transitions (states 'cr2' 'cr3' etc.) in
thermo$obigt (i.e. the Helgeson minerals), it is now possible to use
the minerals in basis(), species(), affinity() with proper accounting
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2017-10-10 15:22:41 UTC (rev 250)
+++ pkg/CHNOSZ/man/water.Rd 2017-10-11 15:40:05 UTC (rev 251)
@@ -83,6 +83,8 @@
\code{de.dP} \tab Pressure derivative \tab bar\eqn{^{-1}}{^-1} \tab * \tab NA \tab NA \cr
\tab of dielectric constant \tab \tab \tab \tab \cr
\code{P} \tab Pressure \tab bar \tab * \tab NA \tab NA \cr
+ \code{A_DH} \tab A Debye-Huckel parameter \tab kg\eqn{^{0.5}}{0.5} mol\eqn{^{-0.5}}{-0.5} \tab * \tab * \tab * \cr
+ \code{B_DH} \tab B Debye-Huckel parameter \tab kg\eqn{^{0.5}}{0.5} mol\eqn{^{-0.5}}{-0.5} cm\eqn{^{-1}}{-1} \tab * \tab * \tab * \cr
}
Call \code{water.SUPCRT92}, \code{water.IAPWS95}, or \code{water.DEW} with no arguments to list the available properties.
@@ -104,6 +106,8 @@
The function also contains routines for calculating the Born functions as numerical derivatives of the static dielectric constant (from \code{water.AW90}).
For compatibility with geochemical modeling conventions, the values of Gibbs energy, enthalpy and entropy output by \code{IAPWS95} are converted by \code{water.IAPWS95} to the triple point reference state adopted in \code{SUPCRT92} (Johnson and Norton, 1991; Helgeson and Kirkham, 1974).
\code{water.IAPWS95} also accepts setting \code{P} to \samp{Psat}, with the saturation pressure calculated from \code{\link{WP02.auxiliary}}; by default the returned properties are for the liquid, but this can be changed to the vapor in \code{\link{thermo}$opt$IAPWS.sat}.
+
+\code{A_DH} and \code{B_DH} are solvent parameters in the \dQuote{B-dot} (extended Debye-Huckel) equation (Helgeson, 1969; Manning, 2013).
}
@@ -112,7 +116,7 @@
}
\seealso{
-Equations for thermodynamic properties of species other than water are coded in \code{\link{hkf}} and \code{\link{cgl}}.
+For calculating properties of reactions, \code{\link{subcrt}} coordinates the calculation of properties among \code{water} and \code{\link{hkf}} and \code{\link{cgl}} for other species.
}
\examples{\dontshow{data(thermo)}
@@ -158,12 +162,16 @@
Helgeson, H. C. and Kirkham, D. H. (1974) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. I. Summary of the thermodynamic/electrostatic properties of the solvent. \emph{Am. J. Sci.} \bold{274}, 1089--1098. \url{http://www.ajsonline.org/cgi/content/abstract/274/10/1089}
+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}
+
Johnson, J. W. and Norton, D. (1991) Critical phenomena in hydrothermal systems: state, thermodynamic, electrostatic, and transport properties of H\eqn{_2}{2}O in the critical region. \emph{Am. J. Sci.} \bold{291}, 541--648. \url{http://www.ajsonline.org/cgi/content/abstract/291/6/541}
Johnson, J. W., Oelkers, E. H. and Helgeson, H. C. (1992) SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000\eqn{^{\circ}}{°}C. \emph{Comp. Geosci.} \bold{18}, 899--947. \url{https://doi.org/10.1016/0098-3004(92)90029-Q}
Levelt-Sengers, J. M. H., Kamgarparsi, B., Balfour, F. W. and Sengers, J. V. (1983) Thermodynamic properties of steam in the critical region. \emph{J. Phys. Chem. Ref. Data} \bold{12}, 1--28. \url{https://doi.org/10.1063/1.555676}
+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}
+
Sverjensky, D. A., Harrison, B. and Azzolini, D. (2014) Water in the deep Earth: The dielectric constant and the solubilities of quartz and corundum to 60 kb and 1,200 \degC. \emph{Geochim. Cosmochim. Acta} \bold{129}, 125--145. \url{https://doi.org/10.1016/j.gca.2013.12.019}
Wagner, W. and Pruss, A. (2002) The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. \emph{J. Phys. Chem. Ref. Data} \bold{31}, 387--535. \url{https://doi.org/10.1063/1.1461829}
Modified: pkg/CHNOSZ/tests/testthat/test-nonideal.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-nonideal.R 2017-10-10 15:22:41 UTC (rev 250)
+++ pkg/CHNOSZ/tests/testthat/test-nonideal.R 2017-10-11 15:40:05 UTC (rev 251)
@@ -1,7 +1,49 @@
context("nonideal")
+# 20161219
test_that("loggam and logK values are consistent", {
rxn1 <- subcrt(c("anhydrite", "Ca+2", "SO4-2"), c(-1, 1, 1), P=1, T=25, I=0)
rxn2 <- subcrt(c("anhydrite", "Ca+2", "SO4-2"), c(-1, 1, 1), P=1, T=25, I=0.7)
expect_equal(rxn1$out$logK, rxn2$out$loggam + rxn2$out$logK)
})
+
+# 20171011
+test_that("A and B parameters are calculated correctly", {
+ ## Psat
+ # values from Helgeson, 1967 "Solution chemistry and metamorphism"
+ # (chapter in http://www.worldcat.org/oclc/152725534)
+ T <- c(25, 50, 100, 200, 300, 350)
+ A <- c(0.5095, 0.5354, 0.6019, 0.8127, 1.2979, 1.9936)
+ B <- c(0.3284, 0.3329, 0.3425, 0.3659, 0.4010, 0.4300)
+ SUP <- water.SUPCRT92(c("A_DH", "B_DH"), T=convert(T, "K"), P="Psat")
+ IAP <- water.IAPWS95(c("A_DH", "B_DH"), T=convert(T, "K"), P="Psat")
+ expect_maxdiff(SUP$A_DH, A, 0.18)
+ expect_maxdiff(IAP$A_DH, A, 0.11)
+ expect_maxdiff(SUP$B_DH / 1e8, B, 0.012)
+ expect_maxdiff(IAP$B_DH / 1e8, B, 0.008)
+
+ # values digitized from Fig. 10 of Manning et al., 2013
+ # doi: 10.2138/rmg.2013.75.5
+ ## 5 kbar
+ T5 <- c(25, seq(100, 1000, 100))
+ A5 <- c(0.441, 0.49, 0.583, 0.685, 0.817, 0.983, 1.164, 1.409, 1.673, 1.938, 2.187)
+ B5 <- c(0.328, 0.336, 0.350, 0.363, 0.377, 0.391, 0.405, 0.421, 0.434, 0.445, 0.453)
+ SUP5 <- water.SUPCRT92(c("A_DH", "B_DH"), T=convert(T5, "K"), P=rep(5000, 11))
+ IAP5 <- water.IAPWS95(c("A_DH", "B_DH"), T=convert(T5, "K"), P=rep(5000, 11))
+ DEW5 <- water.DEW(c("A_DH", "B_DH"), T=convert(T5, "K"), P=rep(5000, 11))
+ # DEW is the winner here
+ expect_maxdiff(SUP5$A_DH, A5, 0.47)
+ expect_maxdiff(IAP5$A_DH, A5, 0.26)
+ expect_maxdiff(DEW5$A_DH, A5, 0.14)
+ expect_maxdiff(SUP5$B_DH / 1e8, B5, 0.036)
+ expect_maxdiff(IAP5$B_DH / 1e8, B5, 0.021)
+ expect_maxdiff(DEW5$B_DH / 1e8, B5, 0.013)
+
+ ## 30 kbar
+ T30 <- seq(700, 1000, 100)
+ A30 <- c(0.625, 0.703, 0.766, 0.815)
+ B30 <- c(0.386, 0.400, 0.409, 0.416)
+ DEW30 <- water.DEW(c("A_DH", "B_DH"), T=convert(T30, "K"), P=rep(30000, 4))
+ expect_maxdiff(DEW30$A_DH, A30, 0.06)
+ expect_maxdiff(DEW30$B_DH / 1e8, B30, 0.024)
+})
Modified: pkg/CHNOSZ/vignettes/obigt.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/obigt.Rmd 2017-10-10 15:22:41 UTC (rev 250)
+++ pkg/CHNOSZ/vignettes/obigt.Rmd 2017-10-11 15:40:05 UTC (rev 251)
@@ -125,7 +125,7 @@
# Recent additions (late 2017)
-* Mineral data using the [Berman (1988)]() equations are listed under **Solids** / **Berman**.
+* Mineral data using the [Berman (1988)](https://doi.org/10.1093/petrology/29.2.445) equations are listed under **Solids** / **Berman**.
* Aqueous species data intended for use with the [Deep Earth Water](http://www.dewcommunity.org/) model are listed under **Optional Data** / **DEW**.
Modified: pkg/CHNOSZ/vignettes/obigt.bib
===================================================================
--- pkg/CHNOSZ/vignettes/obigt.bib 2017-10-10 15:22:41 UTC (rev 250)
+++ pkg/CHNOSZ/vignettes/obigt.bib 2017-10-11 15:40:05 UTC (rev 251)
@@ -36,15 +36,15 @@
}
@Article{AS01,
- author = {Amend, Jan P. and Shock, Everett L.},
- journal = {FEMS Microbiology Reviews},
- title = {{E}nergetics of overall metabolic reactions of thermophilic and hyperthermophilic {A}rchaea and {B}acteria},
- year = {2001},
- volume = {25},
- number = {2},
- pages = {175--243},
- doi = {10.1111/j.1574-6976.2001.tb00576.x},
- issn = {0168-6445},
+ author = {Amend, Jan P. and Shock, Everett L.},
+ journal = {FEMS Microbiology Reviews},
+ title = {{E}nergetics of overall metabolic reactions of thermophilic and hyperthermophilic {A}rchaea and {B}acteria},
+ year = {2001},
+ volume = {25},
+ number = {2},
+ pages = {175--243},
+ doi = {10.1111/j.1574-6976.2001.tb00576.x},
+ issn = {0168-6445},
}
@Article{BH83,
@@ -164,25 +164,25 @@
}
@Article{HDNB78,
- author = {Helgeson, Harold C. and Delany, Joan M. and Nesbitt, H. Wayne and Bird, Dennis K.},
- journal = {American Journal of Science},
- title = {{S}ummary and critique of the thermodynamic properties of rock-forming minerals},
- year = {1978},
- volume = {278A},
- pages = {1--229},
- url = {http://www.worldcat.org/oclc/13594862},
+ author = {Helgeson, Harold C. and Delany, Joan M. and Nesbitt, H. Wayne and Bird, Dennis K.},
+ journal = {American Journal of Science},
+ title = {{S}ummary and critique of the thermodynamic properties of rock-forming minerals},
+ year = {1978},
+ volume = {278A},
+ pages = {1--229},
+ url = {http://www.worldcat.org/oclc/13594862},
}
@Article{HOKR98,
- author = {Helgeson, Harold C. and Owens, Christine E. and Knox, Annette M. and Richard, Laurent},
- journal = {Geochimica et Cosmochimica Acta},
- title = {{C}alculation of the standard molal thermodynamic properties of crystalline, liquid, and gas organic molecules at high temperatures and pressures},
- year = {1998},
- volume = {62},
- number = {6},
- pages = {985--1081},
- doi = {10.1016/S0016-7037(97)00219-6},
- size = {97 p.},
+ author = {Helgeson, Harold C. and Owens, Christine E. and Knox, Annette M. and Richard, Laurent},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {{C}alculation of the standard molal thermodynamic properties of crystalline, liquid, and gas organic molecules at high temperatures and pressures},
+ year = {1998},
+ volume = {62},
+ number = {6},
+ pages = {985--1081},
+ doi = {10.1016/S0016-7037(97)00219-6},
+ size = {97 p.},
}
@Article{HRMNS09,
@@ -239,14 +239,14 @@
}
@Article{LH06a,
- author = {LaRowe, Douglas E. and Helgeson, Harold C.},
- journal = {Geochimica et Cosmochimica Acta},
- title = {{B}iomolecules in hydrothermal systems: {C}alculation of the standard molal thermodynamic properties of nucleic-acid bases, nucleosides, and nucleotides at elevated temperatures and pressures},
- year = {2006},
- volume = {70},
- number = {18},
- pages = {4680--4724},
- doi = {10.1016/j.gca.2006.04.010},
+ author = {LaRowe, Douglas E. and Helgeson, Harold C.},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {{B}iomolecules in hydrothermal systems: {C}alculation of the standard molal thermodynamic properties of nucleic-acid bases, nucleosides, and nucleotides at elevated temperatures and pressures},
+ year = {2006},
+ volume = {70},
+ number = {18},
+ pages = {4680--4724},
+ doi = {10.1016/j.gca.2006.04.010},
}
@Article{LH06b,
@@ -726,7 +726,7 @@
title = {{W}ater in the deep {E}arth: {T}he dielectric constant and the solubilities of quartz and corundum to 60 kb and 1,200 $^\circ${C}},
year = {2014},
volume = {129},
- pages = {125-145},
+ pages = {125--145},
doi = {10.1016/j.gca.2013.12.019},
issn = {0016-7037},
}
@@ -782,7 +782,7 @@
year = {2002},
volume = {31},
number = {2},
- pages = {387 -- 535},
+ pages = {387--535},
doi = {10.1063/1.1461829},
}
@@ -828,7 +828,7 @@
year = {1988},
volume = {29},
number = {2},
- pages = {445},
+ pages = {445--522},
doi = {10.1093/petrology/29.2.445},
}
@@ -928,7 +928,7 @@
year = {2014},
volume = {132},
number = {Supplement C},
- pages = {375 - 390},
+ pages = {375--390},
doi = {10.1016/j.gca.2014.01.030},
issn = {0016-7037},
}
@@ -952,7 +952,7 @@
year = {1992},
volume = {56},
number = {9},
- pages = {3435 - 3467},
+ pages = {3435--3467},
doi = {10.1016/0016-7037(92)90390-5},
issn = {0016-7037},
}
More information about the CHNOSZ-commits
mailing list