Thermochemistry: the quantities and how they connect

Formation quantities and the database ends with the chemical potential of a species written through the apparent Gibbs energy a database stores. This page follows that quantity into the code: the potential the solver minimizes, the models by which its standard part is evaluated at any temperature and pressure, the minimization and the multipliers it returns, and the quantities read from them.

1. The potential the solver minimizes

The terms of the elements and of the charge that the chemical potential contains besides the apparent Gibbs energy change neither an equilibrium composition nor a heat computed between two states (The chemical potential and the apparent Gibbs energy), and the package leaves them out. It works with the remaining potential divided by ,

which the rest of the chapter calls the chemical potential divided by , the difference shifting the component potentials of §4 by constants and leaving compositions and heats unchanged. The two parts come from different places: the standard part belongs to the species and is read from the database (§3), whereas the logarithm of the activity depends on the composition of the phase and is computed by the activity model. build_potentials returns their sum as a closure,

μ(n, p) = p.ΔₐG⁰overRT .+ lna(n, p)      #  μᵢ/RT  =  ΔₐG⁰ᵢ/RT  +  ln aᵢ

The division by makes the potential dimensionless, which is the form the minimizer differentiates, and brings its standard part and the logarithm of the activity to the same scale.

2. Activities, standard states, and the conventions in use

An activity is a ratio to a standard state, and the choice of the standard state is a convention that must be stated for the numbers to have a meaning. The conventions of the package are summarized below and argued, class by class, in Standard states:

class	activity	standard state	convention
aqueous solute		mol/kg of solvent	molality
aqueous solute (ideal model)		mol/L	molarity
solvent (water)		pure water	mole fraction
pure solid, pure phase	⇒	the pure substance	—
gas in an ideal mixture		the pure ideal gas at   bar	mole fraction
solid-solution end-member		the pure end-member	mole fraction

concentration_scale is how a model declares which of the first two it uses, and the aqueous accessors read it rather than guessing:

using ChemistryLab
concentration_scale(DiluteSolutionModel()), concentration_scale(HKFActivityModel())

(:molarity, :molality)

An activity coefficient is thus a property of a species in a given solution rather than of the species alone, and it tends to one in the limit that defines the standard state, infinite dilution for a solute and purity for a solid.

3. How is evaluated

The apparent Gibbs energy of a species is its value of formation at plus integrals of its heat capacity and of its volume (Apparent and formation Gibbs energies). A database supplies, for each species, the reference values and a model for these integrals. Three models cover the species of the databases the package reads: a heat-capacity polynomial for the minerals and the gases, the HKF equations for the solutes, and the equation of state of water for the solvent.

Minerals and gases

Minerals and gases carry a heat-capacity polynomial, the model :cp_ft_equation, of up to eleven terms in powers of , and . THERMO_MODELS (src/thermodynamics/thermo_models.jl) stores together with the analytic antiderivatives of the same coefficients for , and , so that the relations

hold by construction rather than numerically. The model has no pressure term of its own; a record that declares a constant molar volume adds   to its Gibbs energy and enthalpy, as Standard states §1 explains. A heat capacity given on several temperature intervals separated by phase transitions is followed into each of them: the functions are anchored on the interval that contains , carried continuously to the next, and a transition the record places at a boundary adds its enthalpy and its entropy , the Gibbs energy staying continuous; the molar volume stays the record's at . An entry that gives a single heat capacity at is extrapolated with the same model reduced to its constant term; in CEMDATA18 (Lothenbach et al., 2019), the zeolites and the magnesium silicate hydrates are among them, and the solvent's record is of this kind too, though the solvent follows the equation of state of water instead (below). An entry that gives none, but gives the entropy, is extrapolated with a zero heat capacity, so that its apparent Gibbs energy still decreases as with temperature. That zero is a placeholder, heat capacities being positive (Richet, 2001; Sec. 1.2, p. 8). Leaving out a constant raises by      and shifts by  : the enthalpy, and the heats computed with it, are off at first order in   where the Gibbs energy is off at second order only. The additivity of the heat capacities of the constituents, the rule of Neumann and Kopp, is the usual estimate above room temperature (Richet, 2001; Sec. 4.7b, Eq. (4.39), p. 83).

For the Maier-Kelley heat capacity,    , both temperature integrals are elementary, and the closed form can be set against the function the package builds from the same data. With the parameters of calcite used in Thermodynamic Functions, at 500 K:

using DynamicQuantities
Tr, T = 298.15, 500.0
# Thermoddem's calcite, as transcribed in data/literature/Blanc2012.json (SI units)
cal = literature_row("Blanc2012", "individual_properties", "Calcite")
mk = literature_row("Blanc2012", "maier_kelley", "Calcite")
ΔfG, S = ustrip(cal.dfG), ustrip(cal.S)
a₀, a₁, a₂ = ustrip(mk.a), ustrip(mk.b), ustrip(mk.c)

∫Cp = a₀ * (T - Tr) + a₁ / 2 * (T^2 - Tr^2) - a₂ * (1 / T - 1 / Tr)        # ∫ Cp dT
∫Cp_T = a₀ * log(T / Tr) + a₁ * (T - Tr) - a₂ / 2 * (1 / T^2 - 1 / Tr^2)   # ∫ Cp/T dT
closed_form = ΔfG - S * (T - Tr) + ∫Cp - T * ∫Cp_T

dtf = build_thermo_functions(
    :cp_ft_equation,
    Dict(
        :S⁰ => cal.S, :ΔₐH⁰ => cal.dfH, :ΔₐG⁰ => cal.dfG,
        :a₀ => mk.a, :a₁ => mk.b, :a₂ => mk.c, :T => Tr * u"K",
    ),
)
(code = dtf[:ΔₐG⁰](T = T), closed_form = closed_form)

(code = -1.1528438070232493e6, closed_form = -1.1528438070232493e6)

Aqueous solutes

Aqueous solutes use the model :solute_hkf88_reaktoro, the Helgeson-Kirkham-Flowers (HKF) equation of state (Helgeson et al., 1981) in its revised form (Tanger and Helgeson, 1988), whose standard Gibbs energy adds two contributions of different origins to the reference value (Anderson and Crerar, 1993) (§17.9). The first, called nonsolvation, describes the species itself through a heat capacity (coefficients , ) and a volume (coefficients to ) that depend on and , the pressure entering through   and    with   bar. The second is the energy of solvation of a charge in a dielectric continuum, given by the Born equation

where is the relative permittivity of water at and and the Born coefficient of the species, itself a function of and for an ion. The code evaluates it through the Born function    as  , measured from its value at the reference conditions, a term that follows the permittivity of water, which falls from about 78 at 25 °C to about 55 at 100 °C; the same model of water supplies the parameters of the activity models. The standard state of a solute is thus taken at the pressure of the system, and the parameters of an ion are conventional, those of H⁺ being zero at every temperature and pressure.

Where the HKF equations hold

The HKF equations describe a solute in water dense enough to solvate it. The solvent function , which corrects the Born coefficients of the charged species for the compressibility of water, was retrieved by Shock et al. (1992) from dissociation constants of NaCl measured down to a density of 0.35 g/cm³, and SUPCRT92 restricts its calculations to that density (Johnson et al., 1992). A solute's standard properties are computed down to 350 kg/m³ and refused below, with a DomainError, rather than extrapolated: at 350 °C and 1 bar water is a vapor. SUPCRT92 also withholds the properties of charged species between 350 and 400 °C below 500 bar, close to the critical point, where their uncertainties grow large; the package computes them there, and the comparison with the ThermoFun library in the next section leaves out the states above 360 °C and below 600 bar. The solvent and the minerals are computed at every state.

The solvent, from 0 to 1000 °C and up to 5000 bar

The solvent of the ThermoFun databases follows the equation of state of water of Haar et al. (1984) at every temperature and pressure. Written in the steam-table convention, where the entropy and the internal energy of liquid water vanish at its triple point, the equation gives a specific enthalpy and entropy , which SUPCRT92 refers to the properties of liquid water at the triple point (Johnson et al., 1992):

Those constants, which SUPCRT92 takes from Helgeson and Kirkham (1974), are not those of the other records of a database. At 25 °C and 1 bar they put the Gibbs energy 1.35 J/mol and the enthalpy 48.6 J/mol above the values the solvent's record tabulates; and   , which the element entropies of the file fix at for every species so that the enthalpy of a reaction is the one its Gibbs energy implies, comes out 45 J/mol off, where the tabulated values come within 2.4. The solubility products of Cemdata18, from which the energies of its hydrates were derived, close to in against the tabulated water and only to 0.007 against the equation's (ettringite, whose dissolution releases thirty waters). The solvent is therefore built as every other record is, from its tabulated values at and increments from there, the increments being the equation's:

the three constants making the three functions equal to the record's at the reference. The heat capacity is and the volume , so that the five functions are those of one Gibbs energy. At 25 °C and 1 bar the solvent is the record; elsewhere it follows water, where 0.33.0 extrapolated the heat capacity of 25 °C and held the volume at 1 bar constant: the Gibbs energy moves by 0.2 J/mol at most from 0 to 100 °C at 1 bar, the volume by the thermal expansion of water, 0.12 % at 20 °C and 2.6 % at 80 °C, and the Gibbs energy by 0.53 kJ/mol at 300 °C and 1 kbar.

The ThermoFun library uses the equation's functions without the shift, with   K in the constant , a convention that places its Gibbs energy of water 0.633 J/mol below the one computed here. Its values and this package's differ by these constants and by the molar mass of water, 18.015268 g/mol in the library, 18.015 g/mol from the atomic masses here; with the three recomposed, the two agree to at 652 states from 0 to 1000 °C and from 1 to 5000 bar (test/water_eos_reference.jl).

The density is the root of   that the state calls for: the liquid where it exists, the largest root, which a Newton iteration reaches by descending from a density above any the equation attains, the vapor or the supercritical fluid otherwise. The iteration stops on the pressure, to of its scale, never on the step alone, which a stationary point of the isotherm would satisfy as well. The derivatives of the density with respect to and are those of the implicit-function theorem, carried by automatic differentiation.

The apparent convention

This is the quantity every species carries as ΔₐG⁰, and the one the minimizer reads, divided by , as ΔₐG⁰overRT; the subscript is not a variant spelling of . A value of read from a table built in the traditional convention cannot be combined with the apparent energies of a database, since the two differ by    and this difference no longer cancels between species taken from different sources. A second apparent convention also removes the elemental entropies at (Anderson and Crerar, 1993) (§7.4.2); its values differ from the former by the constant  , which is consistent within one database and inconsistent across two. Which convention the code implements can be read on the formula itself: the anchor   rules out the values of that second convention, and the absolute entropy in the linear term rules out the traditional ones.

4. Equilibrium is a constrained minimization, and its dual is the useful part

At fixed and , equilibrium is the composition that minimizes subject to the matter available, which is a linear constraint, since the elements and the charge are conserved whatever the reactions:

where is the conservation matrix cs.SM.A, with one row per component and one column per species, and the budget of the components. The components are the primary species of the system, or the elements with the charge, and the two choices express the same constraints (Formation from primary species). This is the formulation of Leal et al. (2017), the classical minimization of the Gibbs energy with one multiplier per conserved component (Richet, 2001; Sec. 8.6d, Eq. (8.68), p. 194), in which no list of reactions is needed: the reactions are the moves of within the null space of .

The Lagrange multipliers of the equality constraints are the useful output. Writing for the multiplier of row , and     for the potential they give species (the components of   ), the first-order conditions are

The are the component potentials, or element potentials when the components are elements: one number per component, from which the chemical potential of any species follows as a scalar product,   . The two lines above are the complementarity conditions of the Karush-Kuhn-Tucker (KKT) system, and they are what optimality_certificate checks; the certificate itself, and the conditions under which it proves a global minimum, are the subject of Proving that an answer is the answer.

5. The saturation index

The column of for a species is its reaction of formation from the primary species,   , whose Gibbs energy is    and vanishes at equilibrium (Reactions and equilibrium constants). For a solid, the reverse reaction is its dissolution into the primaries, of that dissolution is the ion activity product IAP, the solubility product, and the logarithm of their ratio the saturation index. Since the rows of are labeled by primary species, a component potential is the potential of the corresponding primary, and the index becomes a difference of potentials, with no equilibrium constant to look up,

being that of the reaction of formation. A negative index denotes an undersaturated phase, a zero index a phase in equilibrium with the solution, and a positive index a phase that should have precipitated; this is what saturation_indices returns. A positive index states a direction, not a rate: a precipitate starts from a nucleus whose interface costs energy, so that a solution can stay supersaturated, metastable, for a long time (Richet, 2001; Sec. 12.3d, p. 298). The index also depends on the solid it is computed for, a metastable polymorph or a poorly crystallized form being more soluble than the stable crystal (Richet, 2001; Sec. 12.3d, p. 299).

Two consequences of the identity make the index a check rather than a convention. Since never enters the computation, being implied by the standard potentials, an index cannot disagree with the used by the rest of the calculation. And since every phase present at an equilibrium satisfies the first line of the conditions of §4, its index must vanish to the tolerance of the solve, failing which the state is not an equilibrium and no other index of the result has a meaning. The test suite (test/aqueous_properties.jl) recomputes the index by hand from this formula and compares it with the function.

6. Volume, porosity and chemical shrinkage

Volumes are treated as ideal: the volume of a phase is the sum of the standard molar volumes of its species,

with no excess volume of mixing. This assumption underlies every porosity the package reports. For a solute, is a standard partial molar volume, which can be negative, that of OH⁻ for one (Richet, 2001; Sec. 12.3b, p. 296), and the sum gives the volume of a solution only in the dilute limit; an activity coefficient that depends on pressure carries an excess volume the sum leaves out. volume returns the split by aggregate state, porosity the void fraction relative to a reference state, and chemical_shrinkage the volume the reaction itself consumes, the hydrates occupying less than the water and the clinker they were made from, which is the mechanism of self-desiccation (Self-desiccation).

A species carrying no contributes zero and would corrupt a porosity in silence, which is why the volume machinery reports missing_molar_volumes(state) and why CapillaryWater refuses to be constructed when that list is not empty.

7. What is assumed, and what is not

assumed	not assumed
well-mixed declared phases, each with one composition per instance	a list of reactions, a reaction path, a sequence
ideal molar volumes, no excess volume	ideal activities, which are the business of the activity model
the domain of validity of the activity model	the phases present, which are a result
that   is the whole of the conservation	that a point satisfying the optimality conditions is the minimum

The last cell rests on the certificate. The optimality conditions prove a global minimum when the potentials derive from one convex Gibbs energy, which holds for some activity models and not for others, and fails as well for a solid solution declared inside a miscibility gap or for a constraint that shifts the activity of water; the certificate then reports a KKT point or a composition consistent with its own activities, as stated in What the certificate proves, and when.

Where to go next

Standard states is the next page of the chapter and states, class by class, what each activity of §2 is measured from. The other pages this one leads to are:

• Proving that an answer is the answer, the certificate built on the conditions of §4;

• Activity models, the part of ;

• Solid solutions, the part due to mixing in a solid;

• Chemical Equilibrium, which drives the solver.

Water and quartz from 0 to 1000 °C compares the ionization constant of water and the solubility of quartz computed with the solvent and the solutes of §3 with formulations fitted to the measurements over the whole range.