Description Usage Arguments Details Value References Examples

Calculate activity coefficients and adjusted (non-ideal) molal properties of aqueous species.

1 2 3 |

`species` |
names or indices of species for which to calculate nonideal properties |

`speciesprops` |
list of dataframes of species properties |

`IS` |
numeric, ionic strength(s) used in nonideal calculations, mol kg |

`T` |
numeric, temperature (K) |

`P` |
numeric, pressure (bar); required for Helgeson method |

`A_DH` |
numeric, A Debye-Huckel coefficient; required for Helgeson method |

`B_DH` |
numeric, B Debye-Huckel coefficient; required for Helgeson method |

`method` |
character, Alberty, Helgeson, or Helgeson0 |

`TC` |
numeric, temperature (\degC) |

`showsplines` |
character, show isobaric (T) or isothermal (P) splines |

`nonideal`

takes a list of dataframes (in `speciesprops`

) containing the standard molal properties of the identified `species`

.
The function calculates the *adjusted* properties for given ionic strength (`IS`

); they are equal to the *standard* values at IS=0.
The function bypasses (leaves unchanged) properties of all species whose charge (determined by the number of Z in their `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 `thermo$opt$ideal.H`

or `ideal.e`

to FALSE.
The lengths of `IS`

and `T`

supplied in the arguments should be equal to the number of rows of each dataframe in `speciesprops`

, or length one to use single values throughout.

If `method`

is Alberty, the values of `IS`

are combined with Alberty's (2003) equation 3.6-1 (extended Debye-Hückel equation) 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 G, H, S and Cp; any columns in the dataframes of `speciesprops`

with other names are left untouched.

If `method`

is Helgeson, the “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 “ion size parameter”) is set to 3.72 Angstrom, which is appropriate for NaCl-dominated solutions (Helgeson et al., 1981 Table 2).
In addition to `IS`

and `T`

, this method depends on values of `P`

, `A_DH`

, and `B_DH`

given in the arguments.
The calculation of “B-dot”, also used in the equations, is made within `nonideal`

by calling the `Bdot`

function.
For some uses, it is desirable to set the “B-dot” parameter to zero; this can be done by setting the method to `Helgeson0`

.
Currently, G is the only adjusted molal property that is calculated (but this can be used by `subcrt`

to calculate adjusted equilibrium constants).

`Bdot`

calculates the “B-dot” deviation function (Helgeson, 1969) a.k.a. 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 `water.DEW`

calculations.
If all `P`

correspond to one of the isobaric conditions, the values of `Bdot`

at `T`

are calculated by spline fits to the isobaric data.
Otherwise, particular (dependent on the `T`

) isobaric spline fits are themselves used to construct isothermal splines for the given values of `T`

; the isothermal splines are then used to generate the values of `Bdot`

for the given `P`

.
To see the splines, set `showsplines`

to T to make the first set (isobaric splines) along with the data points, or 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.

One (G) or more (H, S, Cp; currently only with the Alberty method) standard thermodynamic properties (at IS=0) in `speciesprops`

are replaced by the corresponding adjusted thermodynamic properties (at higher IS).
For all affected species, a column named `loggam`

(common (base-10) logarithm of gamma, the activity coefficient) is appended to the output dataframe of species properties.

Alberty, R. A. (2003) *Thermodynamics of Biochemical Reactions*, John Wiley & Sons, Hoboken, New Jersey, 397 p. http://www.worldcat.org/oclc/51242181

Helgeson, H. C. (1969) Thermodynamics of hydrothermal systems at elevated temperatures and pressures. *Am. J. Sci.* **267**, 729–804. 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. *Am. J. Sci.* **281**, 1249–1516. https://doi.org/10.2475/ajs.281.10.1249

Hückel, E. (1925). The theory of concentrated, aqueous solutions of strong electrolytes. *Physikalische Zeitschrift* **26**, 93–147.

Manning, C. E. (2013) Thermodynamic modeling of fluid-rock interaction at mid-crustal to upper-mantle conditions. *Rev. Mineral. Geochem.* **76**, 135–164. 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. *Rev. Mineral. Geochem.* **75**, 109–148. https://doi.org/10.2138/rmg.2013.75.5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 | ```
### Examples following Alberty, 2003
### (page numbers given below)
## the default method setting is Helgeson;
## change it to Alberty
oldnon <- nonideal("Alberty")
## using nonideal() directly
# p. 273-276: activity coefficient (gamma)
# as a function of ionic strength and temperature
IS <- seq(0, 0.25, 0.005)
T <- c(0, 25, 40)
lty <- 1:3
species <- c("H2PO4-", "HADP-2", "HATP-3", "ATP-4")
col <- rainbow(4)
thermo.plot.new(xlim = range(IS), ylim = c(0, 1),
xlab = axis.label("IS"), ylab = "gamma")
for(j in 1:3) {
# use subcrt to generate speciesprops
speciesprops <- subcrt(species, T = rep(T[j], length(IS)))$out
# use nonideal to calculate loggamma; this also adjusts G, H, S, Cp,
# but we don't use them here
nonidealprops <- nonideal(species, speciesprops, IS = IS, T = convert(T[j], "K"))
for(i in 1:4) lines(IS, 10^(nonidealprops[[i]]$loggam), lty=lty[j], col=col[i])
}
t1 <- "Activity coefficient (gamma) of -1,-2,-3,-4 charged species"
t2 <- quote("at 0, 25, and 40 "*degree*"C, after Alberty, 2003")
mtitle(as.expression(c(t1, t2)))
legend("topright", lty=c(NA, 1:3), bty="n",
legend=c(as.expression(axis.label("T")), 0, 25, 40))
legend("top", lty=1, col=col, bty="n",
legend = as.expression(lapply(species, expr.species)))
## more often, the 'IS' argument of subcrt() is used to compute
## adjusted properties at given ionic strength
# p. 16 Table 1.3: adjusted pKa of acetic acid
# set ideal.H to FALSE to calculate activity coefficients for the proton
# (makes for better replication of the values in Alberty's book)
thermo$opt$ideal.H <<- FALSE
(sres <- subcrt(c("acetate", "H+", "acetic acid"), c(-1, -1, 1),
IS=c(0, 0.1, 0.25), T=25, property="logK"))
# we're within 0.01 log of Alberty's pK values
Alberty_logK <- c(4.75, 4.54, 4.47)
stopifnot(maxdiff(sres$out$logK, Alberty_logK) < 0.01)
# reset option to default
thermo$opt$ideal.H <<- TRUE
### An example using IS with affinity():
## speciation of phosphate as a function of ionic strength
opar <- par(mfrow=c(2, 1))
basis("CHNOPS+")
Ts <- c(25, 100)
species(c("PO4-3", "HPO4-2", "H2PO4-"))
for(T in Ts) {
a <- affinity(IS=c(0, 0.14), T=T)
e <- equilibrate(a)
if(T==25) diagram(e, ylim=c(-3.0, -2.6), legend.x=NULL)
else diagram(e, ylim=c(-3.0, -2.6), add=TRUE, col="red")
}
title(main="Non-ideality model for phosphate species")
dp <- describe.property(c("pH", "T", "T"), c(7, Ts))
legend("topright", lty=c(NA, 1, 1), col=c(NA, "black", "red"), legend=dp)
text(0.07, -2.76, expr.species("HPO4-2"))
text(0.07, -2.90, expr.species("H2PO4-"))
## phosphate predominance f(IS,pH)
a <- affinity(IS=c(0, 0.14), pH=c(6, 13), T=Ts[1])
d <- diagram(a, fill=NULL)
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)
### finished with Alberty equation, let's look at Helgeson
# this is the default setting, but is needed here because
# we set "Alberty" above
nonideal(oldnon) # same as nonideal("Helgeson")
## activity coefficients for monovalent ions at 700 degC, 10 kbar
# after Manning, 2010, 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")
wprop <- water(c("A_DH", "B_DH"), T=convert(700, "K"), P=10000)
water(oldwat)
# just use an empty table for a single species
speciesprops <- list(data.frame(G=rep(0, length(IS))))
# choose any monovalent species
(nonidealprops <- nonideal("Na+", speciesprops, IS=IS,
T=convert(700, "K"), P=10000, A_DH=wprop$A_DH, B_DH=wprop$B_DH))
# we get the nonideal Gibbs energy contribution and
# 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
stopifnot(maxdiff(gamma[1], Manning_gamma[1]) < 0.01)
stopifnot(maxdiff(gamma, Manning_gamma) < 0.23)
## data and splines used for calculating B-dot
## (extended term parameter)
Bdot(showsplines = "T")
Bdot(showsplines = "P")
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.