subcrt: Properties of Species and Reactions

Description

Calculate the standard molal thermodynamic properties of one or more species or a reaction between species as a function of temperature and pressure.

Usage

subcrt(species, coeff = 1, state = NULL,
    property = c("logK","G","H","S","V","Cp"),
    T = seq(273.15,623.15,25), P = "Psat", grid = NULL, 
    convert = TRUE, check.Ttr = TRUE, exceed.Ttr = FALSE,
    logact = NULL, action.unbalanced = "warn", IS = 0)

Arguments

species

character, name or formula of species, or numeric, rownumber of species in thermo$obigt

coeff

numeric, reaction coefficients on species

state

character, state(s) of species

property

character, property(s) to calculate

numeric, temperature(s) of the calculation

numeric, pressure(s) of the calculation, or character, Psat

grid

character, type of P$\times$T grid to produce (NULL, the default, means no gridding)

check.Ttr

logical, check if temperatures are beyond transition temperatures?

exceed.Ttr

logical, calculate Gibbs energies of mineral phases and other species beyond their transition temperatures?

logact

numeric, logarithms of activities of species in reaction

convert

logical, are input and output units of T and P those of the user (TRUE) (see T.units), or are they Kelvin and bar (FALSE)?

action.unbalanced

character warn or NULL, what action to take if unbalanced reaction is provided

numeric, ionic strength(s) at which to calculate apparent molal properties, mol kg$^{-1}$

Value

For subcrt, a list of length two or three. If the properties of a reaction were calculated, the first element of the list (named reaction) contains a dataframe with the reaction parameters; the second element, named out, is a dataframe containing the calculated properties. Otherwise, the properties of species (not reactions) are returned: the first element, named species, contains a dataframe with the species identification; the second element, named out, is itself a list, each element of which is a dataframe of properties for a given species. If minerals with phase transitions are present, a third element (a dataframe) in the list indicates for all such minerals the stable phase at each grid point.

Warning

Comparison of the output of subcrt with that of SUPCRT92 indicates the two give similar values of properties for neutral aqueous species up to 1000 $^{\circ}$C and 5000 bar. Changes were made in CHNOSZ 0.9-9 to improve the calculation of the $g$- and $f$-functions (Shock et al., 1992) for the partial derivatives of the omega parameter which are used by the hkf function, so thermodynamic properties of charged aqueous species over the temperature range 0 to 1000 $^{\circ}$C are now mostly consistent with SUPCRT92. Note, however, that while SUPCRT92 limits calculations at Psat to 350 $^{\circ}$C (“beyond range of applicability of aqueous species equations”), there is no corresponding limit in place in subcrt (or hkf), so that inapplicable calculations will be performed up to the critical temperature (373.917 $^{\circ}$C) without any warning. It is probably for this reason that there is a noticeable jump in the value of logK at Psat shown in the one of the examples (demo("NaCl")).

A known bug is misidentification of the stable polymorph of some minerals at high temperature; an example of this bug is shown in the stopifnot at the end of the skarn example below.

Details

subcrt calculates the standard molal thermodynamic properties of species and reactions as a function of temperature and pressure. For each of the species (a formula or name), optionally identified in a given state, the standard molal thermodynamic properties and equations-of-state parameters are retrieved via info (except for liquid). The standard molal properties of the species are computed using equations-of-state functions for aqueous species (hkf), crystalline, gas, and liquid species (cgl) and liquid or supercritical (water).

T and P denote the temperature and pressure conditions for the calculations and should generally be of the same length, unless P is Psat (the default; see water). Argument grid if present can be one of T or P to perform the computation of a T$\times$P or P$\times$T grid. The propertys to be calculated can be one or more of those shown below:

`rho`	Density of water	g cm$^{-3}$
`logK`	Logarithm of equilibrium constant	dimensionless
`G`	Gibbs energy	(cal \| J) mol$^{-1}$
`H`	Enthalpy	(cal \| J) mol$^{-1}$
`S`	Entropy	(cal \| J) K$^{-1}$ mol$^{-1}$
`V`	Volume	cm$^3$ mol$^{-1}$
`Cp`	Heat capacity	(cal \| J) K$^{-1}$ mol$^{-1}$
`E`	Exapansibility	cm$^3$ K$^{-1}$
`kT`	Isothermal compressibility	cm$^3$ bar$^{-1}$

Note that E and kT can only be calculated for aqueous species and only if the option (thermo$opt$water) for calculations of properties using water is set to IAPWS. On the other hand, if the water option is SUPCRT (the default), E and kT can be calculated for water but not for aqueous species. (This is not an inherent limitation in either formulation, but it is just a matter of implementation.)

Depending on the units currently defined (E.units) the values of G, H, S and Cp are returned using calories or Joules as the unit of energy, but only if convert is TRUE. Likewise, the input values of T and P are interpreted to have the units specified through T.units and P.units, but setting convert to FALSE forces subcrt to treat them as Kelvin and bar, respectively.

A chemical reaction is defined if coeff is given. In this mode the standard molal properties of species are summed according to the stoichiometric coefficients, where negative values denote reactants. Reactions that do not conserve elements are permitted; subcrt prints the missing composition needed to balance the reaction and produces a warning but computes anyway. Alternatively, if the basis species of a system were previously defined, and if the species being considered are within the compositional range of the basis species, an unbalanced reaction given in the arguments to subcrt will be balanced automatically, without altering the coefficients on the species identified in the arguments (unless perhaps they correspond to basis species), and without a warning. However, if a reaction is unbalanced and action.unbalanced is set to NULL, no warnings are generated and no attempt is made to balance the reaction.

Minerals with phase transitions (denoted by having states cr1, cr2 etc.) can be defined generically by cr in the state argument with a character value for species. subcrt in this case simultaneously calculates the requested properties of all the phases of each such mineral (phase species) and, using the values of the transition temperatures calculated from those at P = 1 bar given in the thermodynamic database together with functions of the entropies and volumes of transitions (see dPdTtr), determines the stable phase of the mineral at any grid point and substitutes the properties of this phase at that point for further calculations. If individual phase species of minerals are specified (by cr1, cr2 etc. in state), and exceed.Ttr is FALSE (the default), the Gibbs energies of these minerals are assigned values of NA at conditions beyond their transition temperature, where another of the phases is stable. If you set exceed.Ttr to TRUE while investigating the properties of phases of minerals identified generically (by cr) the function will identify the stable phase on the basis not of the transition temperatures but of the calculated Gibbs energies of the phase species. This is not generally advised, since the computed standard molal properties of a phase species lose their physical meaning beyond the transition temperatures of the phase.

If logact is provided, the chemical affinities of reactions are calculated. logact indicates the logarithms of activities (fugacities for gases) of species in the reaction; if there are fewer values of logact than number of species those values are repeated as necessary. If the reaction was unbalanced to start, the logarithms of activities of any basis species added to the reaction are taken from the current definition of the basis species. Columns appended to the output are logQ for the log10 of the activity product of the reaction, and A for the chemical affinity, in the units set by E.units. Note that affinity provides related functionality but is geared toward the properties of formation reactions of species from the basis species and can be performed in more dimensions. Calculations of chemical affinity in subcrt can be performed for any reaction of interest; however, they are currently limited to constant values of the logarithms of activities of species in the reactions, and hence of logQ, across the computational range.

If IS is set to a single value other than zero, nonideal is used to calculate the apparent properties (G, H, S and Cp) of charged aqueous species at the given ionic strength. To perform calculations at a single P and T and for multiple values of ionic strength, supply these values in IS. Calculations can also be performed on a P-IS, T-IS or P,T-IS grid. Values of logK of reactions calculated for IS not equal to zero are consistent with the apparent Gibbs energies of the charged aqueous species.

subcrt is modeled after the functionality of the SUPCRT92 package (Johnson et al., 1992). Certain features of SUPCRT92 are not available here, for example, calculations as a function of density of instead of pressure, or calculations of temperatures of univariant curves (i.e. for which logK is zero), or calculations of the molar volumes of quartz and coesite as a function of temperature (taken here to be constant). The informative messages produced by SUPCRT92 when temperature or pressure limits of the equations of state are exceeded generally are not reproduced here. However, NAs may be produced in the output of subcrt if the requisite thermodynamic or electrostatic properties of water can not be calculated at given conditions. Specifically, NAs are produced for calculations at Psat when the temperature exceeds the critical temperature of .

For calculations below 273.16 K, the pressure should be set to 1, as Psat is not defined in these conditions.

If thermo$opt$varP is TRUE, standard Gibbs energies of gases will be converted from a standard state at 1 bar (as used in SUPCRT) to a variable pressure standard state (see chapter 12 in Anderson and Crerar, 1993). This is useful for constructing e.g. boiling curves for organic compounds.

References

Anderson, G. M. and Crerar, D. A. (1993) Thermodynamics in Geochemistry: The Equilibrium Model, Oxford University Press. http://www.worldcat.org/oclc/803272549

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$^{\circ}$C. Comp. Geosci. 18, 899--947. https://doi.org/10.1016/0098-3004(92)90029-Q

LaRowe, D. E. and Helgeson, H. C. (2007) Quantifying the energetics of metabolic reactions in diverse biogeochemical systems: electron flow and ATP synthesis. Geobiology 5, 153--168. https://doi.org/10.1111/j.1472-4669.2007.00099.x

Schulte, M. D. and Shock, E. L. (1995) Thermodynamics of Strecker synthesis in hydrothermal systems. Orig. Life Evol. Biosph. 25, 161--173. https://doi.org/10.1007/BF01581580

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 $^{\circ}$C and 5 kbar. J. Chem. Soc. Faraday Trans. 88, 803--826. https://doi.org/10.1039/FT9928800803

Examples

Run this code

# NOT RUN {
## properties of species
subcrt("water")
# calculating at different temperatures
subcrt("water", T=seq(0, 100, 10))
# calculating at even increments  	
subcrt("water", T=seq(500, 1000, length.out=10),
  P=seq(5000, 10000, length.out=10))
# calculating on a temperature-pressure grid
subcrt("water", T=c(500, 1000), P=c(5000, 10000), grid="P")
# to calculate only selected properties
subcrt("water", property=c("G", "E"))	
# the properties of multiple species
subcrt(c("glucose", "ethanol", "CO2"))

## properties of reactions
subcrt(c("H2O", "H+", "k-feldspar", "kaolinite", "K+", "SiO2"),
  c(-1, -2, -2, 1, 2, 4))
subcrt(c("glucose", "ethanol", "CO2"), c(-1, 2, 2))
# to specify the states
subcrt(c("glucose", "ethanol", "CO2"), c(-1, 2, 2), c("aq", "aq", "gas")) 

## auto balancing reactions
# the basis species must first be defined
basis(c("CO2", "H2O", "NH3", "H2S", "O2"))
subcrt(c("glucose", "ethanol"), c(-1, 3))
# a bug in CHNOSZ <0.9 caused the following
# to initiate an infinite loop
basis(c("H2O", "H2S", "O2", "H+"))
subcrt(c("HS-", "SO4-2"), c(-1, 1))
# note the next one is a non-reaction
# (products same as reactants)
subcrt("H2O", 1)
# but this one auto-balances into a reaction
# (water to steam)
subcrt("H2O", 1, "gas")
# properties of a species and a formation 
# reaction for that species
subcrt("C2H5OH")    # species
basis("CHNOS")
subcrt("C2H5OH",1)  # reaction

## properties of mineral phases
# properties of phase species
subcrt(c("pyrrhotite", "pyrrhotite"), state=c("cr1", "cr2"))
# properties of the stable phases
subcrt("pyrrhotite")
# phase transitions in a reaction
subcrt(c("pyrite", "pyrrhotite", "H2O", "H2S", "O2"), c(-1, 1, -1, 1, 0.5))

## these produce messages about problems with the calculation
# Psat is not defined above the critical point
subcrt("alanine", T=seq(0, 5000, by=1000))
# above the T, P limits for the H2O equations of state
subcrt("alanine", T=c(2250, 2251), P=c(30000, 30001), grid="T")

## minerals with phase transitions
# compare calculated values of heat capacity of iron with
# values from Robie and Hemingway, 1995
T.units("K")
E.units("J")
# we set pressure here otherwise subcrt uses Psat (saturation 
# vapor pressure of H2O above 100 degrees C) which can not be 
# calculated above the critical point of H2O (~647 K)
s <- subcrt("Fe", T=seq(300, 1800, 20), P=1)
plot(s$out[[1]]$T, s$out[[1]]$Cp, type="l",
  xlab=axis.label("T"), ylab=axis.label("Cp"))
# add points from RH95
RH95 <- read.csv(system.file("extdata/cpetc/RH95.csv", package="CHNOSZ"))
points(RH95[,1], RH95[,2])
title(main=paste("Heat capacity of Fe(cr)\n",
  "(points - Robie and Hemingway, 1995)"))
# reset the units to default values
T.units("C")
E.units("cal")

## Skarn example after Johnson et al., 1992
P <- seq(500, 5000, 500)
# this is like the temperature specification used 
# in the example by Johnson et al., 1992
# T <- seq(0, 1000, 100)
# we use this one to avoid warnings at 0 deg C, 5000 bar
T <- c(2, seq(100, 1000, 100))
s <- subcrt(c("andradite", "carbon dioxide", "H2S", "Cu+", "quartz", 
  "calcite", "chalcopyrite", "H+", "H2O"), 
  coeff=c(-1, -3, -4, -2, 3, 3, 2, 2, 3),
  state=c("cr", "g", "aq", "aq", "cr", "cr", "cr", "aq", "liq"),
  P=P, T=T, grid="P")
# the results are not identical to SUPCRT92, at least because CHNOSZ
# doesn't have volume changes for quartz, and has updated
# parameters for species e.g. Cu+ from Shock et al., 1997
# the following is to help detect unintended changes to the code
# across revisions; the code *should* be fixed sometime so that 
# the last 1 becomes a 3
stopifnot(all.equal(s$state$chalcopyrite[1:11], 
  c(1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 1)))

## Standard Gibbs energy of reactions with HCN and 
## formaldehyde, after Schulte and Shock, 1995 Fig. 1
rxn1 <- subcrt(c("formaldehyde","HCN","H2O","glycolic acid","NH3"),
  c(-1,-1,-2,1,1),P=300)
rxn2 <- subcrt(c("formaldehyde","HCN","H2O","glycine"),
  c(-1,-1,-1,1),P=300)
plot(x=rxn1$out$T,rxn1$out$G/1000,type="l",ylim=c(-40,-10),
  xlab=axis.label("T"),ylab=axis.label("DG0r","k"))
lines(rxn1$out$T,rxn2$out$G/1000)
# write the reactions on the plot
text(150, -14, describe.reaction(rxn1$reaction, iname=c(1,2,4)))
text(200, -35, describe.reaction(rxn2$reaction, iname=c(1,2)))
title(main=paste("Standard Gibbs energy of reactions",
  "after Schulte and Shock, 1995",sep="\n"))

## Calculation of chemical affinities
# after LaRowe and Helgeson, 2007, Fig. 3 (a): reduction of nicotinamide 
# adenine dinucleotide (NAD) coupled to oxidation of glucose
# list the available NAD species
info("NAD ")
T <- seq(0, 120, 10)
# oxidation of glucose (C6H12O6)
basis(c("glucose", "H2O", "NH3", "CO2", "H+"), c(-3, 0, 999, -3, -7))
s <- subcrt(c("NAD(ox)-", "NAD(red)-2"), c(-12, 12), logact=c(0, 0), T=T)
# LH07's diagrams are shown per mole of electron (24 e- per 12 NAD)
A <- s$out$A/24/1000
plot(x=T, y=A, xlim=range(T), ylim=c(1.4, 5.4),
  xlab=axis.label("T"), ylab=axis.label("A", prefix="k"), type="l")
text("NAD(ox)-/NAD(red)-2 = 1", x=53, y=median(A), srt=21)
# different activity ratio
s <- subcrt(c("NAD(ox)-","NAD(red)-2"), c(-12, 12), logact=c(1, 0), T=T)
A <- s$out$A/24/1000
lines(x=T, y=A)
text("NAD(ox)-/NAD(red)-2 = 10", x=55, y=median(A), srt=24)
# different activity ratio
s <- subcrt(c("NAD(ox)-","NAD(red)-2"), c(-12, 12), logact=c(0, 1), T=T)
A <- s$out$A/24/1000
lines(x=T, y=A)
text("NAD(ox)-/NAD(red)-2 = 0.1", x=52, y=median(A), srt=18)
# print the reaction and chemical conditions on the plot
text(0, 5.3, describe.reaction(s$reaction, iname=c(1, 2)), adj=0)
text(0, 5.1, describe.basis(oneline=TRUE, ibasis=c(1, 2, 4, 5)), adj=0)
# label the plot
title(main=paste("Reduction of NAD coupled to oxidation of glucose",
 "after LaRowe and Helgeson, 2007", sep="\n"))

## Subzero (degrees C) calculations
# uncomment the following to try IAPWS95 instead of SUPCRT92
#thermo$opt$water <<- "IAPWS95"
# the limit for H2O92D.f (from SUPCRT92) is currently -20 deg C
# but we go to -30 knowing properties will become NA
sb <- subcrt(c("H2O", "Na+"), T=seq(-30, 10), P=1)$out
# start plot with extra room on right
par(mar=c(5, 4, 4, 4))
# plot G
plot(sb$water$T, sb$water$G, ylim=c(-63000, -56000), xlab=axis.label("T"),
  ylab=axis.label("DG0"))
points(sb$`Na+`$T, sb$`Na+`$G, pch=2)
# add Cp
# change y-axis
par("usr"=c(par("usr")[1:2], -100, 25))
axis(4)
mtext(axis.label("Cp0"), side=4, line=3)
points(sb$water$T, sb$water$Cp, pch=16)
points(sb$`Na+`$T, sb$`Na+`$Cp, pch=17)
legend("topleft", pch=c(16, 1, 17, 2), legend=c("H2O Cp", "H2O G", "Na+ Cp", "Na+ G"))
H2O <- expr.species("H2O")
Na <- expr.species("Na+")
degC <- expr.units("T")
title(main=substitute(H2O~and~Na~to~-20~degC, list(H2O=H2O, Na=Na, degC=degC)))
# }

Run the code above in your browser using DataLab

`rho`	Density of water	g cm\(^{-3}\)
`logK`	Logarithm of equilibrium constant	dimensionless
`G`	Gibbs energy	(cal \| J) mol\(^{-1}\)
`H`	Enthalpy	(cal \| J) mol\(^{-1}\)
`S`	Entropy	(cal \| J) K\(^{-1}\) mol\(^{-1}\)
`V`	Volume	cm\(^3\) mol\(^{-1}\)
`Cp`	Heat capacity	(cal \| J) K\(^{-1}\) mol\(^{-1}\)
`E`	Exapansibility	cm\(^3\) K\(^{-1}\)
`kT`	Isothermal compressibility	cm\(^3\) bar\(^{-1}\)