[CHNOSZ-commits] r343 - in pkg/CHNOSZ: . R data demo inst man man/macros tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 6 06:24:28 CET 2018
Author: jedick
Date: 2018-11-06 06:24:26 +0100 (Tue, 06 Nov 2018)
New Revision: 343
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/NaCl.R
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/data/opt.csv
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/NaCl.Rd
pkg/CHNOSZ/man/data.Rd
pkg/CHNOSZ/man/macros/macros.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/tests/testthat/test-logmolality.R
pkg/CHNOSZ/tests/testthat/test-nonideal.R
Log:
nonideal(): add calculations for neutral species (Setch?\195?\169now equation)
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-06 05:24:26 UTC (rev 343)
@@ -1,6 +1,6 @@
-Date: 2018-11-05
+Date: 2018-11-06
Package: CHNOSZ
-Version: 1.1.3-50
+Version: 1.1.3-51
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/NaCl.R
===================================================================
--- pkg/CHNOSZ/R/NaCl.R 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/R/NaCl.R 2018-11-06 05:24:26 UTC (rev 343)
@@ -3,7 +3,8 @@
# given a total molality of NaCl
# taking account of ion association: Na+ + Cl- = NaCl(aq)
# 20181102 jmd first version
-# 20181105 add activity coefficients of Na+
+# 20181105 use activity coefficient of Na+
+# 20181106 use activity coefficient of NaCl
NaCl <- function(T=seq(100, 500, 100), P=1000, m_tot=2) {
# define a function for the reaction quotient
@@ -24,22 +25,27 @@
ISout <- a_Cl <- numeric(N)
# initial guess for m_Cl and ionic strength assuming complete dissociation of NaCl
IS <- m_Cl <- rep(m_tot, N)
- # the species indices for Cl- and Na+
- ispecies <- info(c("Na+", "Cl-"))
+ # the corresponding total molality of dissolved species (NaCl + Cl- + Na+)
+ m_star <- (m_tot - m_Cl) + 2*m_Cl
+ # the species indices for Na+, Cl-, and NaCl(aq)
+ ispecies <- info(c("Na+", "Cl-", "NaCl"))
# we start by doing calculations for all temperatures
doit <- !logical(N)
while(any(doit)) {
# calculate activity coefficient at given ionic strength
speciesprops <- rep(list(data.frame(G=numeric(N))), length(ispecies))
- gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH))
+ gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH, m_star=m_star))
gam_Na <- 10^gammas[[1]]$loggam
gam_Cl <- 10^gammas[[2]]$loggam
+ gam_NaCl <- 10^gammas[[3]]$loggam
# solve for m_Cl
- for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gam_NaCl=1, gam_Na=gam_Na[i], gam_Cl=gam_Cl[i], logK=logK[i])$root
+ for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gam_NaCl=gam_NaCl[i], gam_Na=gam_Na[i], gam_Cl=gam_Cl[i], logK=logK[i])$root
+ # calculate new total molality
+ m_star <- (m_tot - m_Cl) + 2*m_Cl
# calculate new ionic strength and deviation
ISnew <- m_Cl
dIS <- ISnew - IS
- # set net ionic strength
+ # set new ionic strength
IS <- ISnew
# keep going until the deviation in ionic strength at any temperature is less than 0.01
doit <- abs(dIS) > 0.01
@@ -49,5 +55,5 @@
gam_Na <- 10^gammas[[1]]$loggam
gam_Cl <- 10^gammas[[2]]$loggam
# return the calculated values
- list(IS=IS, m_Cl=m_Cl, gam_Na=gam_Na, gam_Cl=gam_Cl)
+ list(IS=IS, m_Cl=m_Cl, gam_Na=gam_Na, gam_Cl=gam_Cl, gam_NaCl=gam_NaCl)
}
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/R/nonideal.R 2018-11-06 05:24:26 UTC (rev 343)
@@ -62,7 +62,7 @@
else if(prop=="Cp") return(R * T^2 *loggamma(eval(DD(A, "T", 2)), Z, I, B))
}
- # function for Debye-Huckel equation with B-dot extended term parameter (Helgeson, 1969)
+ # function for Debye-Huckel equation with b_gamma or B-dot extended term parameter (Helgeson, 1969)
Helgeson <- function(prop = "loggamma", Z, I, T, A_DH, B_DH, acirc, m_star, bgamma) {
loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) - log10(1 + 0.0180153 * m_star) + bgamma * I
R <- 1.9872 # gas constant, cal K^-1 mol^-1
@@ -71,6 +71,14 @@
# note the log(10) (=2.303) ... use natural logarithm to calculate G
}
+ # function for Setchenow equation with b_gamma or B-dot extended term parameter (Shvarov and Bastrakov, 1999) 20181106
+ Setchenow <- function(prop = "loggamma", I, T, m_star, bgamma) {
+ loggamma <- - log10(1 + 0.0180153 * m_star) + bgamma * I
+ R <- 1.9872 # gas constant, cal K^-1 mol^-1
+ if(prop=="loggamma") return(loggamma)
+ else if(prop=="G") return(R * T * log(10) * loggamma)
+ }
+
# get species indices
if(!is.numeric(species[[1]])) species <- info(species, "aq")
# loop over species #1: get the charge
@@ -79,7 +87,7 @@
# 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)
+ # no charge if Z is absent from the formula or equal to zero
if(is.na(thisZ)) next
if(thisZ==0) next
Z[i] <- thisZ
@@ -111,7 +119,7 @@
# "distance of closest approach" of ions in NaCl solutions (HKF81 Table 2)
acirc <- rep(3.72e-8, length(species))
}
- # get b_gamma (or Bdot)
+ # get b_gamma or B-dot
if(method=="bgamma") bgamma <- bgamma(convert(T, "C"), P)
else if(method=="Bdot") bgamma <- Bdot(convert(T, "C"))
else if(method %in% c("Bdot0", "bgamma0")) bgamma <- 0
@@ -120,35 +128,55 @@
iH <- info("H+")
ie <- info("e-")
speciesprops <- as.list(speciesprops)
- ndid <- 0
+ ncharged <- nneutral <- 0
for(i in 1:length(species)) {
myprops <- speciesprops[[i]]
# 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(!pname %in% c("G", "H", "S", "Cp")) next
- if(method=="Alberty") {
- myprops[, j] <- myprops[, j] + Alberty(pname, Z[i], IS, T)
- didit <- TRUE
- } else {
- myprops[, j] <- myprops[, j] + Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, bgamma)
- didit <- TRUE
+ didcharged <- didneutral <- FALSE
+ # logic for neutral and charged species 20181106
+ if(Z[i]==0) {
+ for(j in 1:ncol(myprops)) {
+ pname <- colnames(myprops)[j]
+ if(!pname %in% c("G", "H", "S", "Cp")) next
+ if(get("thermo")$opt$Setchenow == "bgamma") {
+ myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma)
+ didneutral <- TRUE
+ } else if(get("thermo")$opt$Setchenow == "bgamma0") {
+ myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma = 0)
+ didneutral <- TRUE
+ }
}
+ } else {
+ for(j in 1:ncol(myprops)) {
+ pname <- colnames(myprops)[j]
+ if(!pname %in% c("G", "H", "S", "Cp")) next
+ if(method=="Alberty") {
+ myprops[, j] <- myprops[, j] + Alberty(pname, Z[i], IS, T)
+ didcharged <- TRUE
+ } else {
+ myprops[, j] <- myprops[, j] + Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, bgamma)
+ didcharged <- TRUE
+ }
+ }
}
- # append a loggam column if we did any nonideal calculations of thermodynamic properties
- if(didit) {
+ # append a loggam column if we did any calculations of adjusted thermodynamic properties
+ if(didcharged) {
if(method=="Alberty") myprops <- cbind(myprops, loggam = Alberty("loggamma", Z[i], IS, T))
else myprops <- cbind(myprops, loggam = Helgeson("loggamma", Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, bgamma))
}
+ if(didneutral) {
+ if(get("thermo")$opt$Setchenow == "bgamma") myprops <- cbind(myprops, loggam = Setchenow("loggamma", IS, T, m_star, bgamma))
+ else if(get("thermo")$opt$Setchenow == "bgamma0") myprops <- cbind(myprops, loggam = Setchenow("loggamma", IS, T, m_star, bgamma = 0))
+ }
+ # save the calculated properties and increment progress counters
speciesprops[[i]] <- myprops
- if(didit) ndid <- ndid + 1
+ ncharged <- ncharged + sum(didcharged)
+ nneutral <- nneutral + sum(didneutral)
}
- if(ndid > 0) message("nonideal: calculated activity coefficients for ", ndid, " species (", mettext(method), ")")
+ if(ncharged > 0) message("nonideal: calculations for ", ncharged, " charged species (", mettext(method), ")")
+ if(nneutral > 0) message("nonideal: calculations for ", nneutral, " neutral species (Setchenow equation)")
return(speciesprops)
}
Modified: pkg/CHNOSZ/data/opt.csv
===================================================================
--- pkg/CHNOSZ/data/opt.csv 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/data/opt.csv 2018-11-06 05:24:26 UTC (rev 343)
@@ -1,2 +1,2 @@
-cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS.sat,paramin,ideal.H,ideal.e,nonideal,Berman,maxcores
-1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE,Bdot,NA,2
+cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS.sat,paramin,ideal.H,ideal.e,nonideal,Setchenow,Berman,maxcores
+1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE,Bdot,bgamma0,NA,2
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/demo/gold.R 2018-11-06 05:24:26 UTC (rev 343)
@@ -30,8 +30,10 @@
# set up system
# use H2S here: it's the predominant species at the pH of the QMK buffer -- see sulfur()
-basis(c("Al2O3", "SiO2", "Fe", "Au", "K+", "Cl-", "H2S", "H2O", "oxygen", "H+"))
-# set activity of K+ for 0.5 molal KCl assuming complete dissociation
+basis(c("Al2O3", "quartz", "Fe", "Au", "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
@@ -135,8 +137,12 @@
# calculate solution composition for 2 mol/kg NaCl
NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
a_Cl <- NaCl$m_Cl * NaCl$gam_Cl
+ # using this ionic strength, calculate the activity of K+
+ # assuming complete dissociation of 0.5 mol/kg KCl
+ gam_K <- 10^subcrt("K+", T = seq(150, 550, 10), P = 1000, IS=NaCl$IS)$out$`K+`$loggam
+ a_K <- 0.5 * gam_K
# calculate affinity, equilibrate, solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), P = 1000, IS = NaCl$IS)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), `K+` = log10(a_K), P = 1000, IS = NaCl$IS)
e <- equilibrate(a)
s <- solubility(e)
# make diagram and show total log molality
@@ -145,7 +151,7 @@
# make legend and title
dP <- describe.property("P", 1000)
dNaCl <- expression(NaCl == 2~mol~kg^-1)
- dK <- describe.basis(ibasis=5)
+ dK <- describe.basis(ibasis=5, use.molality=TRUE)
legend("topleft", c(dP, dNaCl, dK), bty = "n")
dbasis <- describe.basis(ibasis = c(9, 7, 10))
legend("topright", dbasis, bty = "n")
@@ -166,8 +172,12 @@
# calculate solution composition for 2 mol/kg NaCl
NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
a_Cl <- NaCl$m_Cl * NaCl$gam_Cl
+ # using this ionic strength, calculate the activity of K+
+ # assuming complete dissociation of 0.5 mol/kg KCl
+ gam_K <- 10^subcrt("K+", T = seq(150, 550, 10), P = 1000, IS=NaCl$IS)$out$`K+`$loggam
+ a_K <- 0.5 * gam_K
# calculate affinity, equilibrate, solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), P = 1000, IS = NaCl$IS)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), `K+` = log10(a_K), P = 1000, IS = NaCl$IS)
e <- equilibrate(a)
s <- solubility(e)
# make diagram and show total log molality
@@ -176,7 +186,7 @@
# make legend and title
dP <- describe.property("P", 1000)
dNaCl <- expression(NaCl == 2~mol~kg^-1)
- dK <- describe.basis(ibasis=5)
+ dK <- describe.basis(ibasis=5, use.molality=TRUE)
legend("topleft", c(dP, dNaCl, dK), bty = "n")
dbasis <- describe.basis(ibasis = c(9, 7, 10))
legend("topright", dbasis, bty = "n")
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-06 05:24:26 UTC (rev 343)
@@ -1,8 +1,20 @@
-CHANGES IN CHNOSZ 1.1.3-49 (2018-11-05)
+CHANGES IN CHNOSZ 1.1.3-51 (2018-11-06)
---------------------------------------
NEW FEATURES
+- Add solubility(). Run this after equilibrate() to calculate the
+ solubility (loga.balance) of the balanced basis species.
+
+- Revise demo/solubility.R to show solubility calculations for CO2(gas)
+ and calcite as a function of T and pH.
+
+- 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.,
+ 2009). This depends on the revised nonideal() and new NaCl() functions
+ described next.
+
- Reorganize and expand options for activity coefficient calculations
(set in thermo$opt$nonideal: Bdot, Bdot0, bgamma, bgamma0, or Alberty).
The previous default, which corresponds to 'bgamma' (T- and
@@ -15,21 +27,16 @@
(Shvarov and Bastrakov, 1999). Parameters are from Table 2.7 of
Garrels and Christ, 1965.
+- nonideal() now calulates activity coefficients of neutral species,
+ using the Setchénow equation. Whether the extended-term parameter in
+ this equation is taken to be zero or is taken from the value for
+ charged species (see above) is controlled by setting
+ 'thermo$opt$Setchenow' to bgamma0 (default) or bgamma.
+
- Add NaCl(), implementing a first-order calculation of the speciation
of NaCl in water, taking account of activity coefficients and the
reaction Na+ + Cl- = NaCl(aq).
-- Add solubility(). Run this after equilibrate() to calculate the
- solubility (loga.balance) of the balanced basis species.
-
-- Revise demo/solubility.R to show solubility calculations for CO2(gas)
- and calcite as a function of T and pH.
-
-- 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.,
- 2009).
-
THERMODYNAMIC DATA
- The Berman data (Berman, 1988 and later additions) have replaced the
Modified: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/man/NaCl.Rd 2018-11-06 05:24:26 UTC (rev 343)
@@ -23,7 +23,7 @@
The algorithm starts by calculating the equilibrium constant (\emph{K}) of the reaction and assuming complete dissociation of NaCl(aq).
This also permits calculating the ionic strength from the molalities of Na\S{+} and Cl\S{-}.
Then, \code{\link{uniroot}} is used to find the equilibrium molality of Cl\S{-}; that is, where the affinity of the reaction (log(\emph{K}/\emph{Q})) becomes zero.
-The activity quotient (\emph{Q}) is evaluated taking account of activity coefficients of Na\S{+} and Cl\S{-} calculated for the nominal ionic strength (see \code{\link{nonideal}}).
+The activity quotient (\emph{Q}) is evaluated taking account of activity coefficients of Na\S{+}, Cl\S{-}, and NaCl(aq) calculated for the nominal ionic strength (see \code{\link{nonideal}}).
The calculated molality of Cl\S{-} yields a new estimate of the ionic strength of the system.
The calculations are iterated until the deviation in ionic strength at all temperatures is less than 0.01.
}
@@ -43,7 +43,7 @@
\examples{\dontshow{data(thermo)}
# ionic strength of solution and activity coefficient of Cl-
-# from HCh (Shvarov and Bastrakov, 1999) at 1000 bar,
+# from HCh version 3.7 (Shvarov and Bastrakov, 1999) at 1000 bar,
# 100, 200, and 300 degress C, and 1 to 6 molal NaCl
m.HCh <- 1:6
IS.HCh <- list(`100`=c(0.992, 1.969, 2.926, 3.858, 4.758, 5.619),
@@ -67,8 +67,11 @@
# plot ionic strength from HCh and NaCl() as points and lines
par(mfrow=c(2, 1))
col <- c("black", "red", "orange")
-plot(c(1,6), c(0,6), xlab="NaCl (mol/kg)", ylab="I (mol/kg)", type="n")
+plot(c(1,6), c(0,6), xlab="NaCl (mol/kg)", ylab=axis.label("IS"), type="n")
for(i in 1:3) {
+ # NOTE: the differences are probably mostly due to different models
+ # for the properties of NaCl(aq) (HCh: B.Ryhzenko model;
+ # CHONSZ: revised HKF with parameters from Shock et al., 1997)
points(m.HCh, IS.HCh[[i]], col=col[i])
lines(m_tot, IS.calc[, i], col=col[i])
}
@@ -77,14 +80,14 @@
dprop <- describe.property(rep("T", 3), c(100, 300, 500))
legend("topleft", dprop, lty=1, pch=1, col=col)
title(main="H2O + NaCl; HCh (points) and 'NaCl()' (lines)")
+plot(c(1,6), c(0,0.8), xlab="NaCl (mol/kg)", ylab=expression(gamma[Cl^"-"]), type="n")
# plot activity coefficient (gamma)
-plot(c(1,6), c(0,1), xlab="NaCl (mol/kg)", ylab=expression(gamma[Cl^"-"]), type="n")
for(i in 1:3) {
points(m.HCh, gam_Cl.HCh[[i]], col=col[i])
lines(m_tot, gam_Cl.calc[, i], col=col[i])
}
# we should be fairly close
-stopifnot(maxdiff(unlist(gam_Cl.calc[seq(1,11,2), ]), unlist(gam_Cl.HCh)) < 0.034)
+stopifnot(maxdiff(unlist(gam_Cl.calc[seq(1,11,2), ]), unlist(gam_Cl.HCh)) < 0.033)
}
\references{
Modified: pkg/CHNOSZ/man/data.Rd
===================================================================
--- pkg/CHNOSZ/man/data.Rd 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/man/data.Rd 2018-11-06 05:24:26 UTC (rev 343)
@@ -58,7 +58,8 @@
\code{paramin} \tab integer \tab Minimum number of calculations to launch parallel processes [1000] (see \code{\link{palply}}) \cr
\code{ideal.H} \tab logical \tab Should \code{\link{nonideal}} ignore the proton? [\code{TRUE}] \cr
\code{ideal.e} \tab logical \tab Should \code{\link{nonideal}} ignore the electron? [\code{TRUE}] \cr
- \code{nonideal} \tab character \tab Method for \code{\link{nonideal}} [\code{Helgeson}] \cr
+ \code{nonideal} \tab character \tab Option for charged species in \code{\link{nonideal}} [\code{Bdot}] \cr
+ \code{Setchenow} \tab character \tab Option for neutral species in \code{\link{nonideal}} [\code{bgamma0}] \cr
\code{Berman} \tab character \tab User data file for mineral parameters in the Berman equations [\code{NA}] \cr
\code{maxcores} \tab numeric \tab Maximum number of cores for parallel calculations with \code{\link{palply}} [\code{2}]
}
Modified: pkg/CHNOSZ/man/macros/macros.Rd
===================================================================
--- pkg/CHNOSZ/man/macros/macros.Rd 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/man/macros/macros.Rd 2018-11-06 05:24:26 UTC (rev 343)
@@ -46,3 +46,4 @@
\newcommand{\Psat}{\ifelse{latex}{\eqn{P_{\mathrm{SAT}}}}{\ifelse{html}{\out{<i>P</i><sub>SAT</sub>}}{Psat}}}
\newcommand{\Delta}{\ifelse{latex}{\eqn{\Delta}}{\ifelse{html}{\out{Δ}}{Δ}}}
\newcommand{\aacute}{\ifelse{latex}{\out{\'{a}}}{\ifelse{html}{\out{á}}{á}}}
+\newcommand{\eacute}{\ifelse{latex}{\out{\'{e}}}{\ifelse{html}{\out{é}}{é}}}
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-11-06 05:24:26 UTC (rev 343)
@@ -4,7 +4,7 @@
\alias{bgamma}
\title{Activity coefficients of aqueous species}
\description{
-Calculate activity coefficients and adjusted (non-ideal) molal properties of aqueous species.
+Calculate activity coefficients and adjusted (transformed) molal properties of aqueous species.
}
\usage{
@@ -28,13 +28,17 @@
}
\details{
-\code{nonideal} takes a list of dataframes (in \code{speciesprops}) containing the standard molal properties of the identified \code{species}.
+This function calculates activity coefficients and adjusted thermodynamic proeprties (i.e. transformed standard Gibbs energy) for charged and neutral aqueous species.
+The species are identified by name or species index in \code{species}.
+\code{speciesprops} is a list of dataframes containing the input standard molal properties; normally, at least one column is \samp{G} for Gibbs energy.
The function calculates the *adjusted* properties for given ionic strength (\code{IS}); they are equal to the *standard* values only at IS=0.
-The function bypasses (leaves unchanged) properties of all species whose charge (determined by the number of Z in their \code{\link{makeup}}) is equal to zero.
-The proton (\Hplus) and electron (\eminus) are also bypassed by default; this makes sense if you are setting the pH, i.e. activity of \Hplus, to some value.
-To apply the calculations to H+ and/or e-, change \code{\link{thermo}$opt$ideal.H} or \code{ideal.e} to FALSE.
The lengths of \code{IS} and \code{T} supplied in the arguments should be equal to the number of rows of each dataframe in \code{speciesprops}, or length one to use single values throughout.
+}
+\section{Charged Species}{
+Calculations for the proton (\Hplus) and electron (\eminus) are skipped by default; this makes sense if you are setting the pH, i.e. activity of \Hplus, to some value.
+To apply the calculations to H+ and/or e-, change \code{\link{thermo}$opt$ideal.H} or \code{ideal.e} to FALSE.
+
If \code{method} is \samp{Alberty}, the values of \code{IS} are combined with Alberty's (2003) equation 3.6-1 (extended Debye-Hückel equation; Hückel, 1925) and its derivatives, to calculate adjusted molal properties at the specified ionic strength(s) and temperature(s).
The adjusted molal properties that can be calculated include \samp{G}, \samp{H}, \samp{S} and \samp{Cp}; any columns in the dataframes of \code{speciesprops} with other names are left untouched.
@@ -52,7 +56,16 @@
The distance of closest approach is set to a constant (3.72e-8 cm) (e.g., Manning et al., 2013).
The extended term parameter is calculated by calling the \code{bgamma} function.
Alternatively, set the extended term parameter to zero with \samp{bgamma0}.
+}
+\section{Neutral Species}{
+For neutral species, the Setch{\eacute}now equation is used, as described in Shvarov and Bastrakov, 1999.
+If \code{\link{thermo}$opt$Setchenow} is \samp{bgamma0} (the default), the extended term parameter is set to zero and the only non-zero term is the mole fraction to molality conversion factor (using the value of \code{m_star}).
+If \code{thermo$opt$Setchenow} is \samp{bgamma}, the extended term paramter is taken from the setting for the charged species (which can be either \samp{Bdot} or \samp{bgamma}).
+Set \code{thermo$opt$Setchenow} to any other value to disable the calculations for neutral species.
+}
+
+\section{b_gamma}{
\code{bgamma} calculates the extended term parameter (written as b_gamma; Helgeson et al., 1981) for activity coefficients in NaCl-dominated solutions at high temperature and pressure.
Data at \Psat and 0.5 to 5 kb are taken from Helgeson (1969, Table 2 and Figure 3) and Helgeson et al. (1981, Table 27) and extrapolated values at 10 to 30 kb from Manning et al. (2013, Figure 11).
Furthermore, the 10 to 30 kb data were used to generate super-extrapolated values at 40, 50, and 60 kb, which may be encountered using the \code{\link{water.DEW}} calculations.
Modified: pkg/CHNOSZ/tests/testthat/test-logmolality.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-06 05:24:26 UTC (rev 343)
@@ -3,39 +3,40 @@
test_that("non-zero ionic strength transforms variables from activity to molality", {
# what happens with activity coefficients when using subcrt() to calculate affinity,
# and in the rest of the main workflow of CHNOSZ?
- # 20171025
+ # 20171025 first version
+ # 20181106 include non-zero activity coefficient of CO2(aq)
- # first get the activity coefficients of H+ and HCO3-
- # the long way...
+ ### first get the activity coefficients of H+ and HCO3-
+ ## the long way...
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)
+ speciesprops <- subcrt(c("H+", "HCO3-", "CO2"), T=25)$out
+ nonid <- nonideal(c("H+", "HCO3-", "CO2"), speciesprops, IS=1, T=298.15, P=1, A_DH=wprop$A_DH, B_DH=wprop$B_DH)
# compare with a precalculated value:
expect_maxdiff(nonid[[2]]$loggam, -0.1868168, 1e-7)
- # the short way...
- loggam <- subcrt(c("H+", "HCO3-"), T=25, IS=1)$out[[2]]$loggam
- expect_equal(nonid[[2]]$loggam, loggam)
-
+ ## the short way...
+ out1 <- subcrt(c("H+", "HCO3-", "CO2"), T=25, IS=1)$out
+ loggam_HCO3 <- out1[[2]]$loggam
+ loggam_CO2 <- out1[[3]]$loggam
+ expect_equal(nonid[[2]]$loggam, loggam_HCO3)
+ expect_equal(nonid[[3]]$loggam, loggam_CO2)
## take-home message -1: with default settings, the activity coefficient of H+ is always 1
- # how do activity coefficient affect the value of G?
+ ### how do activity coefficient affect the value of G?
# let's step back and look at the *standard Gibbs energy* at IS = 0
- out0 <- subcrt(c("H+", "HCO3-"), T=25)$out
- # and at IS = 1
- out1 <- subcrt(c("H+", "HCO3-"), T=25, IS=1)$out
+ out0 <- subcrt(c("H+", "HCO3-", "CO2"), T=25)$out
# the adjusted standard Gibbs energy is less than the standard Gibbs energy
# by an amount determined by the activity coefficient
- expect_equal(out1[[2]]$G - out0[[2]]$G, -convert(loggam, "G"))
-
+ expect_equal(out1[[2]]$G - out0[[2]]$G, -convert(loggam_HCO3, "G"))
+ expect_equal(out1[[3]]$G - out0[[3]]$G, -convert(loggam_CO2, "G"))
## take-home message 0: setting IS in subcrt() gives adjusted standard Gibbs energy
- # now, what is the equilibrium constant for the reaction CO2 + H2O = H+ + HCO3-?
+ # what is the equilibrium constant for the reaction CO2 + H2O = H+ + HCO3-?
# (this is the standard state property at IS=0)
logK <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T=25)$out$logK
# we get logK = -6.344694 (rounded)
expect_maxdiff(logK, -6.344694, 1e-6)
- # now, what is the affinity of the reaction at pH=7 and molalities of HCO3- and CO2 = 10^-3?
-
+ ### what is the affinity of the reaction at pH=7 and molalities of HCO3- and CO2 = 10^-3?
## case 1: ionic strength = 0, so gamma = 0 and activity = molality
# first calculate it by hand from 2.303RTlog(K/Q)
# logQ = (logaH+ + logaHCO3-) - (logaH2O + logaCO2)
@@ -47,37 +48,33 @@
A0subcrt <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T=25, logact=c(-3, 0, -7, -3))$out$A
# we get the same affinity!
expect_equal(A0subcrt, A0manual)
-
- ## case 1: ionic strength = 0, so activity = molality * gamma
- logaHCO3 = -3 + loggam
- logQ1 <- (-7 + logaHCO3) - (0 + -3)
+ ## case 2: ionic strength = 1, so activity = molality * gamma
+ logaHCO3 <- -3 + loggam_HCO3
+ logaCO2 <- -3 + loggam_CO2
+ logQ1 <- (-7 + logaHCO3) - (0 + logaCO2)
A1manual <- -convert(logK - logQ1, "G")
A1subcrt <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T=25, logact=c(-3, 0, -7, -3), IS=1)$out$A
expect_equal(A1subcrt, A1manual)
-
## take-home message 1: using subcrt with IS not equal to zero, the "logact"
## argument is logmolal in affinity calculations for charged aqueous species
- # now, calculate the affinities using affinity()
+ ### now calculate the affinities using affinity()
basis("CHNOS+") # pH=7, logaCO2 = -3
species(c("CO2", "HCO3-")) # logactivities = -3
-
## case 1: IS = 0
a0 <- affinity()
# that gives us values in log units; convert to energy
# (HCO3- is species #2)
A0affinity <- -convert(a0$values[[2]], "G")
expect_equal(A0affinity[[1]], A0subcrt)
-
## case 2: IS = 1
a1 <- affinity(IS=1)
A1affinity <- -convert(a1$values[[2]], "G")
expect_equal(A1affinity[[1]], A1subcrt)
-
## take-home message 2: using affinity() with IS not equal to zero, the "logact"
## set by species() is logmolal in affinity calculations for charged aqueous species
- # now, swap HCO3- for CO2 in the basis
+ ### now, swap HCO3- for CO2 in the basis
swap.basis("CO2", "HCO3-")
basis("HCO3-", -3)
a0 <- affinity()
@@ -99,24 +96,26 @@
ACO2_1manual <- -convert(logKrev - logQrev1, "G")
expect_equal(ACO2_0manual, ACO2_0affinity[[1]])
expect_equal(ACO2_1manual, ACO2_1affinity[[1]])
-
## take-home message 3: using affinity() with IS not equal to zero, the "logact"
## set by basis() is logmolal in affinity calculations for charged aqueous species
- # now look at equilibrate()
+ ### now look at equilibrate()
e0 <- equilibrate(a0)
e1 <- equilibrate(a1)
# using the equilibrated values, calculate affinity of the reaction CO2 + H2O = H+ + HCO3-
# case 1: IS = 0
- logQeq0 <- (-7 + e0$loga.equil[[2]]) - (e0$loga.equil[[1]] + 0)
+ logact_HCO3 <- e0$loga.equil[[2]]
+ logact_CO2 <- e0$loga.equil[[1]]
+ logQeq0 <- (-7 + logact_HCO3) - (logact_CO2 + 0)
Aeq0 <- -convert(logK - logQeq0, "G") # zero!
expect_equal(Aeq0[[1]], 0)
# case 2: IS = 1
+ logact_HCO3 <- e1$loga.equil[[2]]
+ logact_CO2 <- e1$loga.equil[[1]]
# here, loga.equil is the *molality*, so we must multiply by loggam
- logQeq1 <- (-7 + e1$loga.equil[[2]] + loggam) - (e1$loga.equil[[1]] + 0)
+ logQeq1 <- (-7 + logact_HCO3 + loggam_HCO3) - (logact_CO2 + loggam_CO2 + 0)
Aeq1 <- -convert(logK - logQeq1, "G") # zero!
expect_equal(Aeq1[[1]], 0)
-
## take-home message 4: using affinity() with IS not equal to zero, the "loga.equil"
## returned by equilibrate() is logmolal for speciation calculations with charged aqueous species
@@ -124,7 +123,6 @@
a.balance <- 10^e1$loga.balance
m.total <- sum(10^unlist(e1$loga.equil))
expect_equal(a.balance, m.total)
-
## take-home message 5: using affinity() with IS not equal to zero, the "loga.balance"
## used by equilibrate() is the logarithm of total molality of the balancing basis species
})
Modified: pkg/CHNOSZ/tests/testthat/test-nonideal.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-nonideal.R 2018-11-05 11:32:04 UTC (rev 342)
+++ pkg/CHNOSZ/tests/testthat/test-nonideal.R 2018-11-06 05:24:26 UTC (rev 343)
@@ -79,9 +79,11 @@
# 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
+ # calculated with HCh version 3.7 (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)
+ # and the default setting for the Setchenow equation,
+ # for which the only non-zero term is the mole fraction to molality conversion factor
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))
@@ -91,19 +93,28 @@
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))
+ gamNaCl.HCh <- list(`100`=c(0.965, 0.933, 0.904, 0.876, 0.850, 0.827),
+ `300`=c(0.968, 0.941, 0.915, 0.892, 0.870, 0.849),
+ `500`=c(0.977, 0.955, 0.935, 0.915, 0.897, 0.879))
# 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.07)
expect_maxdiff(gamCl.300, gamCl.HCh$`300`, 0.03)
expect_maxdiff(gamCl.500, gamCl.HCh$`500`, 0.009)
- # calculate activity coefficent of Cl- at each temperature
+ # calculate activity coefficent of Na+ 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.08)
expect_maxdiff(gamNa.300, gamNa.HCh$`300`, 0.03)
expect_maxdiff(gamNa.500, gamNa.HCh$`500`, 0.013)
+ # calculate activity coefficent of NaCl at each temperature
+ gamNaCl.100 <- 10^subcrt("NaCl", T=100, P=1000, IS=IS.HCh$`100`)$out$`NaCl`$loggam
+ gamNaCl.300 <- 10^subcrt("NaCl", T=300, P=1000, IS=IS.HCh$`300`)$out$`NaCl`$loggam
+ gamNaCl.500 <- 10^subcrt("NaCl", T=500, P=1000, IS=IS.HCh$`500`)$out$`NaCl`$loggam
+ expect_maxdiff(gamNaCl.100, gamNaCl.HCh$`100`, 0.09)
+ expect_maxdiff(gamNaCl.300, gamNaCl.HCh$`300`, 0.09)
+ expect_maxdiff(gamNaCl.500, gamNaCl.HCh$`500`, 0.10)
})
More information about the CHNOSZ-commits
mailing list