Thermochemistry: the quantities and how they connect

This page is the shortest path from "what is a chemical potential" to the expressions this package actually evaluates, with the notation of the code rather than of a textbook. Everything here is either a definition or an identity; nothing is fitted.

1. One potential, and the two pieces it is made of

The Gibbs energy of a system holding nᵢ moles of each species is  , and the chemical potential of a species is what one more mole of it costs:

Two pieces, and they come from two different places in the package:

• , the standard potential, is thermodynamic data — it belongs to the species alone and comes from the database (§3);

• , the activity, is where the mixture enters, and it comes from the activity model.

Divide by and the code appears verbatim. build_potentials returns the closure

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

so throughout this documentation

Working in is not cosmetic: it is dimensionless, it is what the minimizer differentiates, and it makes and the standard potential directly comparable in size.

The first equality,  , is not a definition but a consequence of the extensivity of . At fixed temperature and pressure, multiplying every amount by multiplies by , and Euler's theorem on functions homogeneous of degree one then expresses as the sum of the amounts weighted by the partial derivatives, which are the partial molar Gibbs energies (Anderson and Crerar, 1993) (§9.2). Being intensive, the cannot vary independently of one another: differentiating   and subtracting   , valid at fixed and , yields the Gibbs-Duhem relation

which holds within each phase. The activities of the species of one phase are thus bound together, and a model cannot prescribe the activity coefficients of the solutes and the activity of the solvent separately; it is by this relation that Activity models §3 obtains the activity of water.

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

An activity is always a ratio to a standard state, and the choice of standard state is a convention that must be stated or the numbers are meaningless. The conventions in use here 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	(exact at )	the pure ideal gas at   bar	mole fraction
solid-solution end-member		the pure end-member	mole fraction (§5)

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 therefore not a property of a species but of a species in a solution, and   in the standard state's own limit — infinite dilution for a solute, purity for a solid.

3. Where comes from

The database stores, per species, the coefficients of a heat-capacity polynomial, and the package integrates it. That is the whole of the temperature dependence, and it is exact rather than a series truncation: THERMO_MODELS (src/thermodynamics/thermo_models.jl) holds and the analytic integrals of the same coefficients for , and , so the three identities

hold by construction rather than numerically. Aqueous solutes use a second 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 electrostatic part is the Born term evaluated with the dielectric constant of water — the same water model that supplies and to the activity models.

The HKF standard Gibbs energy of a solute adds to the reference value two contributions of different origins (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 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.

Two consequences worth carrying:

• a reaction's temperature dependence is not an input; it follows from the of its reactants and products, which is the Gibbs-Helmholtz relation    (van 't Hoff, when is taken constant);

• a species whose database entry has no model cannot be moved off at all, and the package says so rather than extrapolating.

Apparent and formation Gibbs energies

The standard potential is not an absolute energy, since only differences of energy are accessible to measurement. A database stores instead the energy of formation of each species from the chemical elements, the latter being taken in their stable form at the reference conditions   K and   bar, and it is the meaning of this difference away from that deserves attention: two definitions are in use, and they coincide at only (Anderson and Crerar, 1993) (§7.4).

The standard Gibbs energy of formation compares the species with its elements at the same temperature . Its evaluation then requires the heat capacity of every element over the whole range, jumps at their own phase transitions included, although all these contributions disappear from any balanced reaction. The apparent Gibbs energy of formation, in the convention associated with the names of Benson and Helgeson, avoids the detour by leaving the elements at the reference conditions

where denotes the coefficient of element , in its reference form, in the reaction that forms species . The sum being a constant of the species, it cancels from every balanced reaction, so that

holds with apparent quantities exactly as it would with absolute ones, and nothing is required from the elements besides their properties at .

The temperature dependence follows from    and from the identities recalled above. Integrating from , where reduces to the tabulated , yields

which is the expression assembled by build_thermo_functions (src/thermodynamics/thermo_models.jl), the two integrals being the analytic antiderivatives stored in THERMO_MODELS. The third-law entropy of the species multiplies  , and not its entropy of formation, precisely because the elements, frozen at , bring no temperature dependence of their own. The Helgeson-Kirkham-Flowers model of aqueous solutes is written in the same convention, its reference term being and its other terms depending on the species alone, so that a solute and a mineral enter one reaction without any conversion.

For a heat capacity linear in temperature,   , both integrals are elementary, and the closed form can be set against the function the package builds from the same data. With the parameters of gaseous CO₂ used in Thermodynamic Functions, at 500 K:

using DynamicQuantities
Tr, T = 298.15, 500.0
ΔfG, S, a₀, a₁ = -394373.0, 213.785, 33.98, 23.88e-3

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

dtf = build_thermo_functions(
    :cp_ft_equation,
    Dict(
        :S⁰ => S * u"J/K/mol", :ΔₐH⁰ => -393510.0u"J/mol", :ΔₐG⁰ => ΔfG * u"J/mol",
        :a₀ => a₀ * u"J/K/mol", :a₁ => a₁ * u"J/(mol*K^2)", :T => Tr * u"K",
    ),
)
(code = dtf[:ΔₐG⁰](T = T), closed_form = closed_form)

(code = -439937.13910932385, closed_form = -439937.1391093238)

This is the quantity every species carries as ΔₐG⁰, and the one the minimizer reads, divided by , as ΔₐG⁰overRT; the subscript is thus 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, due to Berman and Brown, also removes the entropies of the elements at ; 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 Berman-Brown values, 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. "The matter available" is a linear constraint, because atoms and charge are conserved whatever the reactions do:

where is the conservation matrix cs.SM.A — one row per component, one column per species — and the element (and charge) budget. This is the formulation of (Leal et al., 2017), and it is why no reaction list is needed: reactions are whatever moves inside the null space of .

The Lagrangian multipliers of the equality constraints are the interesting output. Writing for the multiplier of row , the first-order conditions are

The are component potentials (element potentials, when the components are elements): one number per component, from which the chemical potential of any species is a dot product. The two lines above are the complementarity conditions of the KKT system, and they are exactly what optimality_certificate checks — see Proving that an answer is the answer for the certificate itself and why convexity makes it sufficient rather than merely necessary.

5. Reaction quotients, , and the saturation index

Take a species and write its formation from the system's primary species,    — the column of is that reaction. For any reaction,

so  : the sign of says which way the reaction runs, and   is equilibrium. For a solid dissolving into its ions, is the ion activity product IAP and is the solubility product, and the logarithm of their ratio is the saturation index.

Because the rows of are labeled by primary species, a component potential is that primary's , and the index becomes a difference of potentials with no equilibrium constant to look up:

Negative means undersaturated, zero means in equilibrium with the solution, positive means the phase should have precipitated. This is what saturation_indices returns.

Two consequences of the identity are worth carrying, because they are what make the index trustworthy rather than merely conventional:

• no equilibrium constant is looked up. never appears in the computation; it is implied by the standard potentials, so an index cannot disagree with the the rest of the calculation used;

• every phase present at an equilibrium must come out at exactly zero. That is not a property of the answer, it is the first-order condition of §4 restated, so it is a check on the solve: measured on a CEM I paste, the twelve present solids land within  . If they do not, the state is not an equilibrium and no other index in the result means anything.

The identity itself — that the difference of component potentials above equals what saturation_indices computes — is asserted in the test suite (test/aqueous_properties.jl), where it is recomputed by hand from this formula and compared to the function.

6. Volume, porosity and chemical shrinkage

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

with no excess volume of mixing. That is an assumption, and it is the one behind every porosity this package reports. 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 — hydrates occupying less than the water and clinker they were made from, which is the mechanism behind 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 not assumed, and what is

The two lists, side by side:

assumed	not assumed
one well-mixed phase per aggregate state	a reaction list, a reaction path, a sequence
ideal molar volumes, no excess volume	ideal activities — that is the model's business
the activity model's own domain of validity	that the answer is a local minimum: the certificate proves global (see below)
that   is the whole of the conservation	which phases appear: they are a result

Convexity is what makes the last cell true. is convex in for the built-in activity models, so a KKT point is the minimum and the certificate is sufficient rather than necessary — the one documented exception being CapillaryWater, where a composition-dependent shift of breaks the convexity argument and the certificate falls back to proving a KKT point.

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 half of

• Solid solutions — the mole-fraction half

• Chemical Equilibrium — driving the solver