# subcrt: Properties of Species and Reactions In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry

## 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

 ```1 2 3 4 5``` ``` 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, exceed.Ttr = FALSE, exceed.rhomin = FALSE, logact = NULL, autobalance = TRUE, 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 `T` numeric, temperature(s) of the calculation `P` numeric, pressure(s) of the calculation, or character, Psat `grid` character, type of `P`x`T` grid to produce (NULL, the default, means no gridding) `exceed.Ttr` logical, calculate Gibbs energies of mineral phases and other species beyond their transition temperatures? `exceed.rhomin` logical, return properties of species in the HKF model below 0.35 g cm\S-3? `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`)? `autobalance` logical, attempt to automatically balance reaction with basis species? `IS` numeric, ionic strength(s) at which to calculate adjusted molal properties, mol kg^-1

## 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 \H2O 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 \H2O (`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`x`P` or `P`x`T` grid. The `property`s 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 `coeff`icients, 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 `autobalance` is `TRUE`, the `basis` species of a system were previously defined, and all elements in the reaction are represented by the basis species, an unbalanced reaction given in the arguments to `subcrt` will be balanced automatically. The auto balancing doesn't change the reaction coefficients of any species that are do not correspond to basis species.

Minerals with polymorphic transitions (denoted by having states cr (lowest T phase), cr2, cr3 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 cr, 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` to calculate the properties of mineral polymorphs (i.e., using cr) the function will identify the stable polymorph using the calculated Gibbs energies of the phase species instead of the tabulated transition temperatures. 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 adjusted 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 adjusted 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 \H2O instead of pressure, or calculations of temperatures of univariant curves (i.e. for which `logK` is zero).

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.

## 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

Although SUPCRT92 prohibits calculations above 350 \degC at \Psat (“beyond range of applicability of aqueous species equations”), CHNOSZ does not impose this limitation, and allows calculations up to the critical temperature (373.917 \degC) at \Psat. Interpret calculations between 350 \degC and the critical temperature at \Psat at your own risk. The discontinuity in the value of \logK at \Psat that is apparent in `demos("NaCl")` demonstrates one unexpected result.

`NA`s are produced for calculations at Psat when the temperature exceeds the critical temperature of \H2O. In addition, properties of species using the revised HKF equations are set to `NA` wherever the density of \H2O < 0.35 g/cm\S3 (threshold just above the critical isochore; Johnson et al., 1992). Both of these situations produce warnings, which are stored in the warnings element of the return value.

`NA`s are also output if the T, P conditions are otherwise beyond the capabilities of the water equations of state derived from SUPCRT92 (H2O92D.f), but the messages about this are produced by `water.SUPCRT92` rather than `subcrt`.

## 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\degC. Comp. Geosci. 18, 899–947. https://doi.org/10.1016/0098-3004(92)90029-Q

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. Geochim. Cosmochim. Acta 62, 985–1081. https://doi.org/10.1016/S0016-7037(97)00219-6

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 \degC and 5 kbar. J. Chem. Soc. Faraday Trans. 88, 803–826. https://doi.org/10.1039/FT9928800803

`info` can be used to find species in the thermodynamic database. `makeup` is used by `subcrt` for parsing formulas to check mass balance of reactions. `demo("ORP")` and `nonideal` for examples using the `IS` argument.
 ``` 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 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192``` ```## 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)) # because O2,aq is in the basis, this is a non-reaction # (O2,aq to O2,aq) subcrt("O2", 1, "aq") # but this one auto-balances into a reaction # (O2,aq to O2,gas) subcrt("O2", 1, "gas") # properties of a species and a formation # reaction for that species subcrt("C2H5OH") # species basis("CHNOS") subcrt("C2H5OH", 1) # reaction ## mineral polymorphs # properties of the stable polymorph subcrt("pyrrhotite") # properties of just the high-T phase subcrt(c("pyrrhotite"), state="cr2") # polymorphic 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 # above the T, P limits for the H2O equations of state subcrt("alanine", T=c(2250, 2251), P=c(30000, 30001), grid="T") # Psat is not defined above the critical point ## Not run: ## (also gives a warning) subcrt("alanine", T=seq(0, 5000, by=1000)) ## End(Not run) ## 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[]\$T, s\$out[]\$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, as CHNOSZ has updated # parameters for species e.g. Cu+ from Shock et al., 1997. # Check the calculated phase transitions for chalcopyrite stopifnot(all.equal(s\$polymorphs\$chalcopyrite[1:11], c(1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 3))) ## 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 #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 opar <- 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))) par(opar) ## Calculations using a variable-pressure standard state thermo("opt\$varP" = TRUE) # Calculate the boiling point of octane at 2 and 20 bar # We need exceed.Ttr=TRUE because the liquid is metastable # at high temperatures (also, the gas is metastable at low # temperatures, but that doesn't produce NA in the output) sout2 <- subcrt(rep("octane", 2), c("liq", "gas"), c(-1, 1), T=seq(-50, 300, 0.1), P=2, exceed.Ttr=TRUE)\$out sout20 <- subcrt(rep("octane", 2), c("liq", "gas"), c(-1, 1), T=seq(-50, 300, 0.1), P=20, exceed.Ttr=TRUE)\$out # find T with the Gibbs energy of reaction that is closest to zero Tvap2 <- sout2\$T[which.min(abs(sout2\$G))] Tvap20 <- sout20\$T[which.min(abs(sout20\$G))] # the boiling point increases with pressure stopifnot(Tvap20 > Tvap2) # more precisely, the calculated boiling points should be near the # empirical values (digitized from Fig. 1 of Helgeson et al., 1998) Tvap_2bar <- 156 Tvap_20bar <- 276 stopifnot(abs(Tvap2 - Tvap_2bar) < 6) stopifnot(abs(Tvap20 - Tvap_20bar) < 25) # those comparisons would fail if varP were FALSE (the default) thermo("opt\$varP" = FALSE) ```