[CHNOSZ-commits] r252 - in pkg/CHNOSZ: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 12 10:21:01 CEST 2017
Author: jedick
Date: 2017-10-12 10:21:01 +0200 (Thu, 12 Oct 2017)
New Revision: 252
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/man/water.Rd
Log:
add Bdot() for calculating Debye-Huckel/Helgeson extended term parameter
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-11 15:40:05 UTC (rev 251)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-12 08:21:01 UTC (rev 252)
@@ -1,6 +1,6 @@
-Date: 2017-10-10
+Date: 2017-10-12
Package: CHNOSZ
-Version: 1.1.0-50
+Version: 1.1.0-51
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2017-10-11 15:40:05 UTC (rev 251)
+++ pkg/CHNOSZ/NAMESPACE 2017-10-12 08:21:01 UTC (rev 252)
@@ -59,7 +59,7 @@
# added 20170301 or later
"GHS_Tr", "calculateDensity", "calculateGibbsOfWater",
"calculateEpsilon", "calculateQ", "water.DEW", "berman",
- "maxdiff", "expect_maxdiff"
+ "maxdiff", "expect_maxdiff", "Bdot"
)
# Load shared objects
@@ -67,7 +67,7 @@
# Imports from default packages
importFrom("grDevices", "dev.cur", "dev.off", "extendrange",
- "heat.colors", "png", "rainbow")
+ "heat.colors", "png", "rainbow", "topo.colors")
importFrom("graphics", "abline", "axTicks", "axis", "barplot", "box",
"contour", "image", "legend", "lines", "mtext", "par", "plot",
"plot.new", "plot.window", "points", "rect", "text", "title")
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2017-10-11 15:40:05 UTC (rev 251)
+++ pkg/CHNOSZ/R/nonideal.R 2017-10-12 08:21:01 UTC (rev 252)
@@ -33,6 +33,12 @@
else if(prop=='S') return(- R * T * loggamma(eval(DD(A,'T',1)),Z,I,B))
else if(prop=='Cp') return(R * T^2 *loggamma(eval(DD(A,'T',2)),Z,I,B))
}
+ Helgeson <- function() {
+ # "distance of closest approach" of ions in NaCl solutions
+ # HKF81 Table 2
+ acirc <- 3.72 # Angstrom
+ }
+
if(!is.numeric(species[[1]])) species <- info(species,'aq')
proptable <- as.list(proptable)
# which gamma function to use
@@ -67,3 +73,247 @@
return(proptable)
}
+Bdot <- function(TC=25, P=1, showsplines="") {
+ # 20171012 calculate B-dot (bgamma) using P, T, points from:
+ # Helgeson, 1969 (doi:10.2475/ajs.267.7.729)
+ # Helgeson et al., 1981 (doi:10.2475/ajs.281.10.1249)
+ # Manning et al., 2013 (doi:10.2138/rmg.2013.75.5)
+ # T in degrees C
+ T <- TC
+ # are we at a pre-fitted constant pressure?
+ uP <- unique(P)
+ is1 <- identical(uP, 1) & all(T==25)
+ is500 <- identical(uP, 500)
+ is1000 <- identical(uP, 1000)
+ is2000 <- identical(uP, 2000)
+ is3000 <- identical(uP, 3000)
+ is4000 <- identical(uP, 4000)
+ is5000 <- identical(uP, 5000)
+ is10000 <- identical(uP, 10000)
+ is20000 <- identical(uP, 20000)
+ is30000 <- identical(uP, 30000)
+ is40000 <- identical(uP, 40000)
+ is50000 <- identical(uP, 50000)
+ is60000 <- identical(uP, 60000)
+ isoP <- is1 | is500 | is1000 | is2000 | is3000 | is4000 | is5000 | is10000 | is20000 | is30000 | is40000 | is50000 | is60000
+ # values for Bdot x 100 from Helgeson (1969), Figure (P = Psat)
+ if(!isoP | showsplines != "") {
+ T0 <- c(23.8, 49.4, 98.9, 147.6, 172.6, 197.1, 222.7, 248.1, 268.7)
+ B0 <- c(4.07, 4.27, 4.30, 4.62, 4.86, 4.73, 4.09, 3.61, 1.56)
+ # we could use the values from Hel69 Table 2 but extrapolation of the
+ # their fitted spline function turns sharply upward above 300 degC
+ #T0a <- c(25, 50, 100, 150, 200, 250, 270, 300)
+ #B0a <- c(4.1, 4.35, 4.6, 4.75, 4.7, 3.4, 1.5, 0)
+ S0 <- splinefun(T0, B0)
+ }
+ # values for bgamma x 100 from Helgeson et al., 1981 Table 27
+ if(is500 | !isoP | showsplines != "") {
+ T0.5 <- seq(0, 400, 25)
+ B0.5 <- c(5.6, 7.1, 7.8, 8.0, 7.8, 7.5, 7.0, 6.4, 5.7, 4.8, 3.8, 2.6, 1.0, -1.2, -4.1, -8.4, -15.2)
+ S0.5 <- splinefun(T0.5, B0.5)
+ if(is500) return(S0.5(T))
+ }
+ if(is1000 | !isoP | showsplines != "") {
+ T1 <- seq(0, 500, 25)
+ B1 <- c(6.6, 7.7, 8.7, 8.3, 8.2, 7.9, 7.5, 7.0, 6.5, 5.9, 5.2, 4.4, 3.5, 2.5, 1.1, -0.6, -2.8, -5.7, -9.3, -13.7, -19.2)
+ S1 <- splinefun(T1, B1)
+ if(is1000) return(S1(T))
+ }
+ if(is2000 | !isoP | showsplines != "") {
+ # 550 and 600 degC points from Manning et al., 2013 Fig. 11
+ T2 <- c(seq(0, 500, 25), 550, 600)
+ B2 <- c(7.4, 8.3, 8.8, 8.9, 8.9, 8.7, 8.5, 8.1, 7.8, 7.4, 7.0, 6.6, 6.2, 5.8, 5.2, 4.6, 3.8, 2.9, 1.8, 0.5, -1.0, -3.93, -4.87)
+ S2 <- splinefun(T2, B2)
+ if(is2000) return(S2(T))
+ }
+ if(is3000 | !isoP | showsplines != "") {
+ T3 <- seq(0, 500, 25)
+ B3 <- c(6.5, 8.3, 9.2, 9.6, 9.7, 9.6, 9.4, 9.3, 9.2, 9.0, 8.8, 8.6, 8.3, 8.1, 7.8, 7.5, 7.1, 6.6, 6.0, 5.4, 4.8)
+ S3 <- splinefun(T3, B3)
+ if(is3000) return(S3(T))
+ }
+ if(is4000 | !isoP | showsplines != "") {
+ T4 <- seq(0, 500, 25)
+ B4 <- c(4.0, 7.7, 9.5, 10.3, 10.7, 10.8, 10.8, 10.8, 10.7, 10.6, 10.5, 10.4, 10.3, 10.2, 10.0, 9.8, 9.6, 9.3, 8.9, 8.5, 8.2)
+ S4 <- splinefun(T4, B4)
+ if(is4000) return(S4(T))
+ }
+ if(is5000 | !isoP | showsplines != "") {
+ # 550 and 600 degC points from Manning et al., 2013 Fig. 11
+ T5 <- c(seq(0, 500, 25), 550, 600)
+ B5 <- c(0.1, 6.7, 9.6, 11.1, 11.8, 12.2, 12.4, 12.4, 12.4, 12.4, 12.4, 12.3, 12.3, 12.2, 12.1, 11.9, 11.8, 11.5, 11.3, 11.0, 10.8, 11.2, 12.52)
+ S5 <- splinefun(T5, B5)
+ if(is5000) return(S5(T))
+ }
+ # 10, 20, and 30 kb points from Manning et al., 2013 Fig. 11
+ # here, one control point at 10 degC is added to make the splines curve down at low T
+ if(is10000 | !isoP | showsplines != "") {
+ T10 <- c(25, seq(300, 1000, 50))
+ B10 <- c(12, 17.6, 17.8, 18, 18.2, 18.9, 21, 23.3, 26.5, 28.8, 31.4, 34.1, 36.5, 39.2, 41.6, 44.1)
+ S10 <- splinefun(T10, B10)
+ if(is10000) return(S10(T))
+ }
+ if(is20000 | !isoP | showsplines != "") {
+ T20 <- c(25, seq(300, 1000, 50))
+ B20 <- c(16, 21.2, 21.4, 22, 22.4, 23.5, 26.5, 29.2, 32.6, 35.2, 38.2, 41.4, 44.7, 47.7, 50.5, 53.7)
+ S20 <- splinefun(T20, B20)
+ if(is20000) return(S20(T))
+ }
+ if(is30000 | !isoP | showsplines != "") {
+ T30 <- c(25, seq(300, 1000, 50))
+ B30 <- c(19, 23.9, 24.1, 24.6, 25.2, 26.7, 30.3, 32.9, 36.5, 39.9, 43, 46.4, 49.8, 53.2, 56.8, 60)
+ S30 <- splinefun(T30, B30)
+ if(is30000) return(S30(T))
+ }
+ # 40-60 kb points extrapolated from 10-30 kb points of Manning et al., 2013
+ if(is40000 | !isoP | showsplines != "") {
+ T40 <- c(seq(300, 1000, 50))
+ B40 <- c(25.8, 26, 26.4, 27.2, 28.9, 33, 35.5, 39.2, 43.2, 46.4, 49.9, 53.4, 57.1, 61.2, 64.4)
+ S40 <- splinefun(T40, B40)
+ if(is40000) return(S40(T))
+ }
+ if(is50000 | !isoP | showsplines != "") {
+ T50 <- c(seq(300, 1000, 50))
+ B50 <- c(27.1, 27.3, 27.7, 28.5, 30.5, 34.8, 37.3, 41.1, 45.5, 48.7, 52.4, 55.9, 59.8, 64.3, 67.5)
+ S50 <- splinefun(T50, B50)
+ if(is50000) return(S50(T))
+ }
+ if(is60000 | !isoP | showsplines != "") {
+ T60 <- c(seq(300, 1000, 50))
+ B60 <- c(28, 28.2, 28.6, 29.5, 31.6, 36.1, 38.6, 42.5, 47.1, 50.4, 54.1, 57.6, 61.6, 66.5, 69.7)
+ S60 <- splinefun(T60, B60)
+ if(is60000) return(S60(T))
+ }
+ # show points and spline(T) curves
+ if(showsplines == "T") {
+ thermo.plot.new(c(0, 1000), c(-20, 70), xlab=axis.label("T"), ylab=expression(b[gamma]%*%100))
+ points(T0, B0, pch=0)
+ points(T0.5, B0.5, pch=1)
+ points(T1, B1, pch=1)
+ points(T2[-c(22:23)], B2[-c(22:23)], pch=1)
+ points(T2[c(22:23)], B2[c(22:23)], pch=2)
+ points(T3, B3, pch=1)
+ points(T4, B4, pch=1)
+ points(T5[-c(22:23)], B5[-c(22:23)], pch=1)
+ points(T5[c(22:23)], B5[c(22:23)], pch=2)
+ points(T10[-1], B10[-1], pch=2)
+ points(T20[-1], B20[-1], pch=2)
+ points(T30[-1], B30[-1], pch=2)
+ points(T10[1], B10[1], pch=5)
+ points(T20[1], B20[1], pch=5)
+ points(T30[1], B30[1], pch=5)
+ points(T40, B40, pch=6)
+ points(T50, B50, pch=6)
+ points(T60, B60, pch=6)
+ col <- rev(topo.colors(13))
+ T0 <- seq(0, 350, 5); lines(T0, S0(T0), col=col[1])
+ T0.5 <- seq(0, 500, 5); lines(T0.5, S0.5(T0.5), col=col[2])
+ T1 <- seq(0, 500, 5); lines(T1, S1(T1), col=col[3])
+ T2 <- seq(0, 600, 5); lines(T2, S2(T2), col=col[4])
+ T3 <- seq(0, 600, 5); lines(T3, S3(T3), col=col[5])
+ T4 <- seq(0, 600, 5); lines(T4, S4(T4), col=col[6])
+ T5 <- seq(0, 600, 5); lines(T5, S5(T5), col=col[7])
+ T10 <- c(25, seq(100, 1000, 5)); lines(T10, S10(T10), col=col[8])
+ T20 <- c(80, seq(100, 1000, 5)); lines(T20, S20(T20), col=col[9])
+ T30 <- c(125, seq(200, 1000, 5)); lines(T30, S30(T30), col=col[10])
+ T40 <- c(175, seq(300, 1000, 5)); lines(T40, S40(T40), col=col[11])
+ T50 <- c(225, seq(300, 1000, 5)); lines(T50, S50(T50), col=col[12])
+ T60 <- c(250, seq(300, 1000, 5)); lines(T60, S60(T60), col=col[13])
+ legend("topleft", pch=c(0, 1, 2, 5, 6),
+ legend=c("Helgeson, 1969", "Helgeson et al., 1981", "Manning et al., 2013", "spline control point", "high-P extrapolation"))
+ legend("bottomright", col=c(NA, rev(col)), lty=1,
+ legend=c("kbar", "60", "50", "40", "30", "20", "10", "5", "4", "3", "2", "1", "0.5", "Psat"))
+ } else if(showsplines=="P") {
+ #thermo.plot.new(c(0, 60000), c(-20, 70), xlab=axis.label("P"), ylab=expression(b[gamma]%*%100))
+ thermo.plot.new(c(0, 5), c(-20, 70), xlab="log P", ylab=expression(b[gamma]%*%100))
+ # pressures that are used to make the isothermal splines (see below)
+ P25 <- c(1, 500, 1000, 2000, 3000, 4000, 5000)
+ P100 <- c(1, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000)
+ P200 <- c(15, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000)
+ P300 <- c(86, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ P400 <- c(500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ P500 <- c(1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ P600 <- c(2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ P700 <- c(10000, 20000, 30000, 40000, 50000, 60000)
+ P800 <- c(10000, 20000, 30000, 40000, 50000, 60000)
+ P900 <- c(10000, 20000, 30000, 40000, 50000, 60000)
+ P1000 <- c(10000, 20000, 30000, 40000, 50000, 60000)
+ # plot the pressure and B-dot values used to make the isothermal splines
+ points(log10(P25), Bdot(25, P25) * 100)
+ points(log10(P100), Bdot(100, P100) * 100)
+ points(log10(P200), Bdot(200, P200) * 100)
+ points(log10(P300), Bdot(300, P300) * 100)
+ points(log10(P400), Bdot(400, P400) * 100)
+ points(log10(P500), Bdot(500, P500) * 100)
+ points(log10(P600), Bdot(600, P600) * 100)
+ points(log10(P700), Bdot(700, P700) * 100)
+ points(log10(P800), Bdot(800, P800) * 100)
+ points(log10(P900), Bdot(900, P900) * 100)
+ points(log10(P1000), Bdot(1000, P1000) * 100)
+ # plot the isothermal spline functions
+ col <- tail(rev(rainbow(12)), -1)
+ P <- c(1, seq(50, 5000, 50)); lines(log10(P), Bdot(25, P) * 100, col=col[1])
+ P <- c(1, seq(50, 20000, 50)); lines(log10(P), Bdot(100, P) * 100, col=col[2])
+ P <- c(1, seq(50, 40000, 50)); lines(log10(P), Bdot(200, P) * 100, col=col[3])
+ P <- c(1, seq(50, 60000, 50)); lines(log10(P), Bdot(300, P) * 100, col=col[4])
+ P <- seq(500, 60000, 50); lines(log10(P), Bdot(400, P) * 100, col=col[5])
+ P <- seq(1000, 60000, 50); lines(log10(P), Bdot(500, P) * 100, col=col[6])
+ P <- seq(2000, 60000, 50); lines(log10(P), Bdot(600, P) * 100, col=col[7])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(700, P) * 100, col=col[8])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(800, P) * 100, col=col[9])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(900, P) * 100, col=col[10])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(1000, P) * 100, col=col[11])
+ legend("topleft", col=c(NA, col), lty=1, legend=c("degrees C", 25, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000))
+ legend("bottomright", pch=1, legend="points from iso-P splines")
+ } else {
+ # make T and P the same length
+ ncond <- max(length(T), length(P))
+ T <- rep(T, length.out=ncond)
+ P <- rep(P, length.out=ncond)
+ # loop over P, T conditions
+ Bdot <- numeric()
+ lastT <- NULL
+ for(i in 1:length(T)) {
+ # make it fast: skip splines at 25 degC and 1 bar
+ if(T[i]==25 & P[i]==1) Bdot <- c(Bdot, 0.041)
+ else {
+ if(!identical(T[i], lastT)) {
+ # get the spline fits from particular pressures for each T
+ if(T[i] >= 700) {
+ PT <- c(10000, 20000, 30000, 40000, 50000, 60000)
+ B <- c(S10(T[i]), S20(T[i]), S30(T[i]), S40(T[i]), S50(T[i]), S60(T[i]))
+ } else if(T[i] >= 600) {
+ PT <- c(2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ B <- c(S2(T[i]), S3(T[i]), S4(T[i]), S5(T[i]), S10(T[i]), S20(T[i]), S30(T[i]), S40(T[i]), S50(T[i]), S60(T[i]))
+ } else if(T[i] >= 500) {
+ PT <- c(1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ B <- c(S1(T[i]), S2(T[i]), S3(T[i]), S4(T[i]), S5(T[i]), S10(T[i]), S20(T[i]), S30(T[i]), S40(T[i]), S50(T[i]), S60(T[i]))
+ } else if(T[i] >= 400) {
+ PT <- c(500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ B <- c(S0.5(T[i]), S1(T[i]), S2(T[i]), S3(T[i]), S4(T[i]), S5(T[i]), S10(T[i]), S20(T[i]), S30(T[i]), S40(T[i]), S50(T[i]), S60(T[i]))
+ } else if(T[i] >= 300) {
+ # here the lowest P is Psat
+ PT <- c(86, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
+ B <- c(S0(T[i]), S0.5(T[i]), S1(T[i]), S2(T[i]), S3(T[i]), S4(T[i]), S5(T[i]), S10(T[i]), S20(T[i]), S30(T[i]), S40(T[i]), S50(T[i]), S60(T[i]))
+ } else if(T[i] >= 200) {
+ # drop highest pressures because we get into ice
+ PT <- c(16, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000)
+ B <- c(S0(T[i]), S0.5(T[i]), S1(T[i]), S2(T[i]), S3(T[i]), S4(T[i]), S5(T[i]), S10(T[i]), S20(T[i]), S30(T[i]), S40(T[i]))
+ } else if(T[i] >= 100) {
+ PT <- c(1, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000)
+ B <- c(S0(T[i]), S0.5(T[i]), S1(T[i]), S2(T[i]), S3(T[i]), S4(T[i]), S5(T[i]), S10(T[i]), S20(T[i]))
+ } else if(T[i] >= 0) {
+ PT <- c(1, 500, 1000, 2000, 3000, 4000, 5000)
+ B <- c(S0(T[i]), S0.5(T[i]), S1(T[i]), S2(T[i]), S3(T[i]), S4(T[i]), S5(T[i]))
+ }
+ # make a new spline as a function of pressure at this T
+ ST <- splinefun(PT, B)
+ # remember this T; if it's the same as the next one, we won't re-make the spline
+ lastT <- T[i]
+ }
+ Bdot <- c(Bdot, ST(P[i]) / 100)
+ }
+ }
+ return(Bdot)
+ }
+}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-11 15:40:05 UTC (rev 251)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-12 08:21:01 UTC (rev 252)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-50 (2017-10-11)
+CHANGES IN CHNOSZ 1.1.0-51 (2017-10-12)
---------------------------------------
MAJOR CHANGES:
@@ -36,9 +36,15 @@
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.
+ calculate 'A_DH' and 'B_DH', i.e. A and B coefficients in the extended
+ Debye-Huckel equation of Helgeson, 1969.
+- Add Bdot() to calculate the "B-dot" (or b_gamma) extended term
+ parameter for activity coefficients in NaCl solutions at high
+ temperature and pressure (Helgeson et al., 1981) including
+ high-pressure extrapolations > 10 kbar based on data from Manning et
+ al., 2013.
+
- 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/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2017-10-11 15:40:05 UTC (rev 251)
+++ pkg/CHNOSZ/man/eos.Rd 2017-10-12 08:21:01 UTC (rev 252)
@@ -85,7 +85,7 @@
\references{
- 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\eqn{^{\circ}}{°}C and 5 Kb. \emph{Am. J. Sci.} \bold{281}, 1249--1516. \url{http://www.ajsonline.org/cgi/content/abstract/281/10/1249}
+ 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}
Helgeson, H. C., Owens, C. E., Knox, A. M. and Richard, L. (1998) Calculation of the standard molal thermodynamic properties of crystalline, liquid, and gas organic molecules at high temperatures and pressures. \emph{Geochim. Cosmochim. Acta} \bold{62}, 985--1081. \url{https://doi.org/10.1016/S0016-7037(97)00219-6}
@@ -97,7 +97,7 @@
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{https://doi.org/10.1039/FT9928800803}
- Tanger, J. C. IV and Helgeson, H. C. (1988) Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: Revised equations of state for the standard partial molal properties of ions and electrolytes. \emph{Am. J. Sci.} \bold{288}, 19--98. \url{http://www.ajsonline.org/cgi/content/abstract/288/1/19}
+ Tanger, J. C. IV and Helgeson, H. C. (1988) Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: Revised equations of state for the standard partial molal properties of ions and electrolytes. \emph{Am. J. Sci.} \bold{288}, 19--98. \url{https://doi.org/10.2475/ajs.288.1.19}
}
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2017-10-11 15:40:05 UTC (rev 251)
+++ pkg/CHNOSZ/man/nonideal.Rd 2017-10-12 08:21:01 UTC (rev 252)
@@ -1,6 +1,7 @@
\encoding{UTF-8}
\name{nonideal}
\alias{nonideal}
+\alias{Bdot}
\title{Activity coefficients of aqueous species}
\description{
Calculate activity coefficients and non-ideal contributions to apparent standard molal properties of aqueous species.
@@ -8,13 +9,17 @@
\usage{
nonideal(species, proptable, IS, T)
+ Bdot(TC, P, showsplines = "")
}
\arguments{
\item{species}{names or indices of species for which to calculate nonideal properties}
\item{proptable}{list of dataframes of species properties}
\item{IS}{numeric, ionic strength(s) used in nonideal calculations, mol kg\eqn{^{-1}}{^-1}}
- \item{T}{numeric, temperature (K) (\code{lines.water}, \code{nonideal})}
+ \item{T}{numeric, temperature (K)}
+ \item{TC}{numeric, temperature (\degC)}
+ \item{P}{numeric, pressure (bar)}
+ \item{showsplines}{character, show isobaric (\samp{T}) or isothermal (\samp{P}) splines}
}
\details{
@@ -25,6 +30,14 @@
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{proptable}, or one to use single values throughout.
The apparent molal properties that can be calculated include \code{G}, \code{H}, \code{S} and \code{Cp}; any columns in the dataframes of \code{proptable} with other names are left untouched.
A column named \code{loggam} (logarithm of gamma, the activity coefficient) is appended to the output dataframe of species properties.
+
+\code{Bdot} calculates the \dQuote{B-dot} deviation function (Helgeson, 1969) or extended term parameter (written as b_gamma; Helgeson et al., 1981) for activity coefficients in NaCl 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.
+If all \code{P} correspond to one of the isobaric conditions, the values of \code{Bdot} at \code{T} are calculated by spline fits to the isobaric data.
+Otherwise, particular (dependent on the \code{T}) isobaric spline fits are themselves used to construct isothermal splines for the given values of \code{T}; the isothermal splines are then used to generate the values of \code{Bdot} for the given \code{P}.
+To see the splines, set \code{showsplines} to \samp{T} to make the first set (isobaric splines) along with the data points, or \samp{P} for examples of isothermal splines at even temperature intervals (here, the symbols are not data, but values generated from the isobaric splines).
+This is a crude method of kriging the data, but produces fairly smooth interpolations without adding any external dependencies.
}
\section{Warning}{
@@ -100,10 +113,21 @@
a <- affinity(IS=c(0, 0.14), pH=c(6, 13), T=Ts[2])
d <- diagram(a, add=TRUE, names=NULL, col="red")
par(opar)
+
+## data and splines used for calculating B-dot
+## (extended term parameter)
+Bdot(showsplines = "T")
+Bdot(showsplines = "P")
}
\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}
+
+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}
+
+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}
}
\concept{Secondary thermodynamic modeling}
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2017-10-11 15:40:05 UTC (rev 251)
+++ pkg/CHNOSZ/man/water.Rd 2017-10-12 08:21:01 UTC (rev 252)
@@ -160,11 +160,11 @@
Haar, L., Gallagher, J. S. and Kell, G. S. (1984) \emph{NBS/NRC Steam Tables}. Hemisphere, Washington, D. C., 320 p. \url{http://www.worldcat.org/oclc/301304139}
-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. 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{https://doi.org/10.2475/ajs.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. 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{https://doi.org/10.2475/ajs.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}
More information about the CHNOSZ-commits
mailing list