| nonideal {CHNOSZ} | R Documentation |
Calculate activity coefficients and adjusted molal properties of aqueous species.
nonideal(species, speciesprops, IS, T, P, A_DH, B_DH,
m_star = NULL, method = thermo()$opt$nonideal)
bgamma(TC, P, showsplines = "")
species |
name of method to use, or 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 B-dot or b_gamma equation |
A_DH |
numeric, A Debye-Huckel coefficient; required for B-dot or b_gamma equation |
B_DH |
numeric, B Debye-Huckel coefficient; required for B-dot or b_gamma equation |
m_star |
numeric, total molality of all dissolved species |
method |
character, ‘Bdot’, ‘Bdot0’, ‘bgamma’, ‘bgamma0’, or ‘Alberty’ |
TC |
numeric, temperature (°C) |
showsplines |
character, show isobaric (‘T’) or isothermal (‘P’) splines |
nonideal calculates activity coefficients and adjusted thermodynamic properties for charged and neutral aqueous species.
At the user level, the main use of this function is to set the method for activity coefficient calculations that gets used by other functions in CHNOSZ.
See the “Charged Species” section for a description of the available methods.
Activity coefficient calculations are activated by setting the IS argument of subcrt or affinity.
Those functions then call nonideal with all the arguments needed to perform the calculations.
The default is to not apply calculations for the proton (H+) and electron (e-); this makes sense if you are setting the pH, i.e. activity of H+, to some value.
To apply the calculations to H+ and/or e-, change thermo()$opt$ideal.H or ideal.e to FALSE (see examples).
For the ‘Alberty’ method, the values of IS are combined with Alberty's (2003) equation 3.6-1 (extended Debye-Hückel equation with an empirical term valid up to 0.25 M ionic strength) and its derivatives (Alberty, 2001), to calculate adjusted molal properties at the specified ionic strength(s) and temperature(s).
The calculations use the equation for the Debye-Hückel constant given by Clarke and Glew, 1980, which is valid between 0 and 150 °C at saturated water vapor pressure (PSAT).
For the ‘Bdot’ method (the default), the “B-dot” form of the extended Debye-Hückel equation is used.
This equation is valid at ionic strengths up to approximately 3 mol / kg (Hörbrand et al., 2018).
The distance of closest approach for different ions (the “ion size parameter”) is taken from thermo()$Bdot_acirc; any species not listed in this file is assigned a value of 4.5 Å.
The extended term parameter for NaCl-dominated solutions, known as “B-dot”, is calculated as a function only of temperature (Helgeson, 1969).
To set the extended term parameter to zero, use the ‘Bdot0’ method.
For the ‘bgamma’ method, the “b_gamma” equation is used.
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 bgamma function.
Alternatively, set the extended term parameter to zero with ‘bgamma0’.
For neutral species, the Setchénow equation is used, as described in Shvarov and Bastrakov, 1999.
If thermo()$opt$Setchenow is ‘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 m_star).
If thermo()$opt$Setchenow is ‘bgamma’, the extended term paramter is taken from the setting for the charged species (which can be either ‘Bdot’ or ‘bgamma’).
Set thermo()$opt$Setchenow to any other value to disable the calculations for neutral species.
This information, about the arguments and return values used to perform the calculations, is not normally needed by the user (but the usage is shown in the last example).
For nonideal, the species can be identified by name or species index in species.
speciesprops is a list of dataframes containing the input standard molal properties; normally, at least one column is ‘G’ for Gibbs energy.
The function calculates the *adjusted* properties for given ionic strength (IS); they are equal to the *standard* values only at IS=0.
The adjusted molal properties that can be calculated include ‘G’, and (currently only with the Alberty method) ‘H’, ‘S’ and ‘Cp’; values of any columns with other names are left untouched.
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.
In addition to IS and T, the ‘Bdot’ and ‘bgamma’ methods depend on values of P, A_DH, B_DH, and m_star given in the arguments.
m_star, the total molality of all dissolved species, is used to compute the mole fraction to molality conversion factor.
If m_star is NULL, it is taken to be equal to IS, which is an underestimate.
For these methods, ‘G’ is currently the only adjusted molal property that is calculated (but this can be used by subcrt to calculate adjusted equilibrium constants).
The return value is the same as the input in speciesprops, except the input standard thermodynamic properties (at IS=0) are replaced by adjusted 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.
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 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 basic method of interpolating the data without adding any external dependencies.
Alberty, R. A. (2001) Effect of temperature on standard transformed Gibbs energies of formation of reactants at specified pH and ionic strength and apparent equilibrium constants of biochemical reactions. J. Phys. Chem. B 105, 7865–7870. doi:10.1021/jp011308v
Alberty, R. A. (2003) Thermodynamics of Biochemical Reactions, John Wiley & Sons, Hoboken, New Jersey, 397 p. https://www.worldcat.org/oclc/51242181
Clarke, E. C. W. and Glew, D. N. (1980) Evaluation of Debye-Hückel limiting slopes for water between 0 and 150 °C. J. Chem. Soc. Faraday Trans. 76, 1911–1916. doi:10.1039/f19807601911
Helgeson, H. C. (1969) Thermodynamics of hydrothermal systems at elevated temperatures and pressures. Am. J. Sci. 267, 729–804. doi: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°C and 5 Kb. Am. J. Sci. 281, 1249–1516. doi:10.2475/ajs.281.10.1249
Hörbrand, T., Baumann, T. and Moog, H. C. (2018) Validation of hydrogeochemical databases for problems in deep geothermal energy. Geotherm. Energy 6, 20. doi:10.1186/s40517-018-0106-3
Manning, C. E. (2013) Thermodynamic modeling of fluid-rock interaction at mid-crustal to upper-mantle conditions. Rev. Mineral. Geochem. 76, 135–164. doi: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. doi:10.2138/rmg.2013.75.5
Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. Australian Geological Survey Organisation 1999/25. https://pid.geoscience.gov.au/dataset/ga/25473
## Each of the available methods
nonideal("Alberty")
nonideal("bgamma0")
nonideal("bgamma")
nonideal("Bdot0")
nonideal("Bdot") # the default
## What's the activity coefficient of Na+ at
## 25 degC and 1 bar and an ionic strength of 0.7?
sres <- subcrt("Na+", T = 25, IS = 0.7)
# Exponentiate to convert log10(gamma) to gamma
10^sres$out[[1]]$loggam
# Now use a different method
nonideal("bgamma")
sres <- subcrt("Na+", T = 25, IS = 0.7)
10^sres$out[[1]]$loggam
## What are activity coefficients of -3, -2, -1, 0, +1, +2, +3 charged species
## as a function of ionic strength and temperature?
# First choose the method
nonideal("Bdot")
# Define the ionic strength and temperature increments
IS <- c(0.001, 0.01, 0.1, 0.7)
T <- seq(0, 100, 25)
# Use species with charges -3, -2, -1, 0, +1, +2, +3
species <- c("PO4-3", "HPO4-2", "H2PO4-", "H3PO4", "Na+", "Ca+2", "Al+3")
# Initialize empty output table for T (rows) and charge (columns)
gamtab <- matrix(nrow = length(T), ncol = length(species))
rownames(gamtab) <- T
colnames(gamtab) <- -3:3
# Make a list of tables to hold the activity coefficients, one for each IS
gamma <- rep(list(gamtab), length(IS))
names(gamma) <- IS
# Loop over the values of ionic strength
for(i in seq_along(IS)) {
# Calculate properties of species, including logarithm of activity coefficient
sres <- subcrt(species, T = T, IS = IS[i])
# Exponentiate to convert log10(gamma) to gamma, and put the values into the tables
for(j in seq_along(species)) gamma[[i]][, j] <- 10^sres$out[[j]]$loggam
}
# Print the output and make a plot
print(gamma)
matplot(T, gamma$`0.001`, type = "l")
title(main = "activity coefficients of -3, -2, -1, 0, +1, +2, +3 charged species")
## Alberty, 2003 p. 16 Table 1.3: adjusted pKa of acetic acid
## We use the 'IS' argument in subcrt() to calculate adjusted thermodynamic properties
# Set ideal.H to FALSE to calculate activity coefficients for the proton
# (needed for replication of the values in Alberty's book)
nonideal("Alberty")
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 of Alberty's pK values
Alberty_logK <- c(4.75, 4.54, 4.47)
# The maximum (absolute) pairwise difference between x and y
max(abs(Alberty_logK - sres$out$logK)) # 0.0072
# 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, add = TRUE, names = FALSE, 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 = FALSE, col = "red")
par(opar)
## Activity coefficients for monovalent ions at 700 degC, 10 kbar
# After Manning, 2013, Fig. 7
# Here we use the b_gamma equation
nonideal("bgamma")
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")
sres <- subcrt("Na+", T = 700, P = 10000, IS = IS)
water(oldwat)
# Compare the calculated activity coefficient to values from Manning's figure
gamma <- 10^sres$out[[1]]$loggam
Manning_gamma <- c(0.93, 0.82, 0.65, 0.76, 1.28, 2)
gamma - Manning_gamma
## Plot the data and splines used for calculating b_gamma
## (extended term parameter)
bgamma(showsplines = "T")
bgamma(showsplines = "P")
## A longer example, using nonideal() directly
# Alberty, 2003 p. 273-276: activity coefficient (gamma)
# as a function of ionic strength and temperature
nonideal("Alberty")
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)))
## Reset method to default
nonideal("Bdot") # or reset()