Kinetics under partial equilibrium: the right-hand side and the calorimeters
Before this page
Rate laws, and every parameter in them, which writes the rates
A hydrating paste holds reactions on two time scales. The aqueous speciation, protonation, complexation and the autoprotolysis of water, reaches equilibrium in microseconds; alite dissolves over days. Integrating every reaction with a rate law would need kinetic constants that are not measured for the fast ones and would bring the integrator down to their time scale; equilibrating every phase would dissolve the clinker at once. Partial equilibrium integrates the slow reactions and solves the fast ones by a minimization of the Gibbs energy at each instant (Leal et al., 2015; Leal et al., 2017).
The equilibria computed this way are conditional on the kinetic amounts of the instant. They are not jumps over activation barriers, and they do not say that the paste has reached its stable equilibrium: they are minima of
Section 1 sets out the matrix every later equation is written with, the stoichiometric matrix, over the elements and over the primary species. The package then advances such a system in two ways, which answer two questions. Sections 2 to 4 write the differential system of Leal et al. (2015), whose right-hand side holds the minimization and whose Jacobian is exact; Section 5 the implicit step of Leal et al. (2017), where the step itself is one minimization; Section 6 extends the first to the energy balance of a calorimeter.
1. The stoichiometric matrix
Everything below rests on one matrix: what each species is made of. The partition, the state of the integration and the reactions it advances are all written with it, and the two ways of writing it, over the elements or over a few chosen species, are the source of the one error a kinetic scheme makes most easily, a budget kept in one basis and spent in the other.
Over the elements
A species is a column of numbers, the atoms of each element in its formula and its charge. Placed side by side, the columns of the Zz in the package, and the amounts of the elements in a composition
the balance every equilibrium and every reaction conserves (Leal et al., 2017). The carbonate system, water with dissolved carbon dioxide and its ions, is the smallest one that shows everything that follows:
using ChemistryLab, LinearAlgebra
aq(f, c = SC_AQSOLUTE) = Species(f; aggregate_state = AS_AQUEOUS, class = c)
H2O = aq("H2O", SC_AQSOLVENT)
Hp, OHm, CO2, HCO3m, CO3mm = aq("H+"), aq("OH-"), aq("CO2"), aq("HCO3-"), aq("CO3-2")
carbonate_species = [H2O, Hp, OHm, CO2, HCO3m, CO3mm]
W = CanonicalStoichMatrix(carbonate_species)
pprint(W; label = :symbol)┌────┬─────┬────┬─────┬─────┬───────┬───────┐
│ │ H2O │ H+ │ OH- │ CO2 │ HCO3- │ CO3-2 │
├────┼─────┼────┼─────┼─────┼───────┼───────┤
│ C │ │ │ │ 1 │ 1 │ 1 │
│ H │ 2 │ 1 │ 1 │ │ 1 │ │
│ O │ 1 │ │ 1 │ 2 │ 3 │ 3 │
│ Zz │ │ 1 │ -1 │ │ -1 │ -2 │
└────┴─────┴────┴─────┴─────┴───────┴───────┘Four rows, six columns, and yet only three of the rows are independent:
println("rank W = ", rank(Float64.(W.A)))rank W = 3The fourth row is a combination of the other three. Weighting the hydrogen by
Primary species: a basis of the columns
Linear algebra says more. The
and each decomposition is a balanced reaction, the formation of the secondary species
The same decomposition rewrites the whole formula matrix. With
SM.A in the package, of
SM = StoichMatrix(carbonate_species, [H2O, Hp, CO2])
pprint(SM; label = :symbol)┌─────┬─────┬────┬─────┬─────┬───────┬───────┐
│ │ H2O │ H+ │ OH- │ CO2 │ HCO3- │ CO3-2 │
├─────┼─────┼────┼─────┼─────┼───────┼───────┤
│ H2O │ 1 │ │ 1 │ │ 1 │ 1 │
│ H+ │ │ 1 │ -1 │ │ -1 │ -2 │
│ CO2 │ │ │ │ 1 │ 1 │ 1 │
└─────┴─────┴────┴─────┴─────┴───────┴───────┘The column of the carbonate ion reads
pprint(StoichMatrix(carbonate_species, [H2O, Hp, HCO3m]); label = :symbol)┌───────┬─────┬────┬─────┬─────┬───────┬───────┐
│ │ H2O │ H+ │ OH- │ CO2 │ HCO3- │ CO3-2 │
├───────┼─────┼────┼─────┼─────┼───────┼───────┤
│ H2O │ 1 │ │ 1 │ -1 │ │ │
│ H+ │ │ 1 │ -1 │ 1 │ │ -1 │
│ HCO3- │ │ │ │ 1 │ 1 │ 1 │
└───────┴─────┴────┴─────┴─────┴───────┴───────┘while the budgets it describes are the same compositions:
# mol, in the order of `carbonate`; the proton balances the anions
n = [55.5, 1e-4 + 1e-10 + 2e-9, 1e-10, 1e-2, 1e-4, 1e-9]
A₁ = Float64.(SM.A)
A₂ = Float64.(StoichMatrix(carbonate_species, [H2O, Hp, HCO3m]).A)
println("over H₂O, H⁺, CO₂: b = ", round.(A₁ * n; digits = 6))
println("over H₂O, H⁺, HCO₃⁻: b = ", round.(A₂ * n; digits = 6))
println("over the elements: b = ", round.(Float64.(W.A) * n; digits = 6))over H₂O, H⁺, CO₂: b = [55.5001, 0.0, 0.0101]
over H₂O, H⁺, HCO₃⁻: b = [55.49, 0.0101, 0.0101]
over the elements: b = [0.0101, 111.0002, 55.5203, 0.0]Three budgets, one composition. In the first basis the row of the proton is the row of the charge, so its coordinate is the charge of the solution, zero here, and it would be negative for a solution holding a base: a coordinate, not an amount. Which basis is chosen changes the numbers, not the constraints; what it may not do is change in the middle of a calculation. The package uses the primaries of the system throughout, a choice made once (Formation from primary species for the energies that go with it).
Reactions: the null space
A reaction is a vector of coefficients
and any other balanced reaction is a combination of these. A minimization of the Gibbs energy needs no list of reactions at all: it moves
for r in reactions(SM)
println(r.equation)
endH₂O = OH⁻ + H⁺
H₂O + CO₂ = HCO₃⁻ + H⁺
H₂O + CO₂ = CO₃²⁻ + 2H⁺ν = [r[s] for r in reactions(SM), s in carbonate_species] # one row per reaction, r[s] < 0 for a reactant
println("A νᵀ = 0: ", iszero(SM.A * transpose(ν)))A νᵀ = 0: true2. The partition, and why the state carries element amounts
The species are split into a kinetic partition, of amounts
The obvious state, every amount
with
On calcite dissolving in the carbonate solution, calcite alone is kinetic. Its canonical reaction on the primaries is its formation; turned so that calcite carries
Cal = Species("CaCO3"; aggregate_state = AS_CRYSTAL, class = SC_COMPONENT)
Cap = aq("Ca+2")
species = [carbonate_species; Cap; Cal]
SMc = StoichMatrix(species, [H2O, Hp, CO2, Cap])
formation = only(r for r in reactions(SMc) if r[Cal] != 0)
println(formation.equation)
ν = [formation[s] for s in species]
ν = ν ./ -ν[end] # the dissolution: calcite carries −1
e, k = 1:7, 8:8 # the two partitions
show(b) = join(("$(symbol(p)) $(x)" for (p, x) in zip(SMc.primaries, b)), ", ")
println("A_e ν_eᵀ = ", show(SMc.A[:, e] * ν[e])) # per mole dissolved
println("−A_k ν_kᵀ = ", show(-SMc.A[:, k] * ν[k]))H₂O + CO₂ + Ca²⁺ = CaCO₃ + 2H⁺
A_e ν_eᵀ = H2O 1.0, H+ -2.0, CO2 1.0, Ca+2 1.0
−A_k ν_kᵀ = H2O 1.0, H+ -2.0, CO2 1.0, Ca+2 1.0The budget moves by the calcite's column, the two expressions are one, and no amount of any equilibrium species is touched: the minimization decides where the released calcium and carbon go. Advancing the amounts themselves instead would remove two protons per mole from a solution that holds almost none. In a cement paste at pH 13 a liter holds
n_H = 1e-13 # mol of H⁺ in the solution
ξ = 1e-3 # mol of calcite dissolved
println("H⁺ advanced along the reaction: ", round(n_H + ν[2] * ξ; sigdigits = 3), " mol")H⁺ advanced along the reaction: -0.002 molThe coordinate of the budget on the proton may well be negative, it is a coordinate; the amounts the minimization returns for it are not, the protons "consumed" being taken from the hydroxide and the water, as the equilibrium says.
Three properties of the construction follow, each of which an implementation can get wrong while conserving matter:
is the matrix the minimization is posed on. The formula matrix over the elements and the matrix over the primary species of the system are different matrices, the same conservation in two bases (§1); a budget built in one basis and a minimization posed in the other is infeasible at every step. The equilibrium sub-system inherits the primary species of the parent system, so that , its derivative and the minimization share one basis. The minimization runs over the equilibrium partition only. Posed on the whole system it would equilibrate the kinetic minerals instantaneously, which is what a kinetic description exists to prevent. The sub-system keeps the phases of the parent whose members are all in the partition, its solid solutions and its site families, on an all-or-nothing rule: a family split between the two partitions is refused, its members sharing one budget.
A site family is fast or slow as a whole. Its budget is a conservation row of the minimization, and moving one member across would take the row with it. A family no reaction controls belongs to the equilibrium partition, its states redistributing as fast as the aqueous speciation. A family one of whose members a kinetic reaction controls goes, whole, to the kinetic partition: A slow surface.
A slow surface
When equilibrium is the right description. Partial equilibrium asks the reactions it treats as equilibrated to be fast compared with what changes the composition, a criterion that is pragmatic, set by the error one accepts and by the time scales of the other processes (Rubin, 1983; p. 1232). For a surface, the time scales reported are short. Greskowiak et al. (2011) describe the surface complexation of metals as reaching equilibrium within milliseconds, citing pressure-jump relaxation measurements, and attribute the nonequilibrium observed in sediments to diffusion-limited exchange (p. 1). Furrer and Stumm (1986) take the protonation of an oxide surface as equilibrated with the solution while the detachment of a metal center limits the dissolution (p. 1850). On calcium silicate hydrates, Tits et al. (2006) find the uptake of strontium fast, with equilibrium reached within one day (Sec. 3.2, p. 81). A cement paste changes over hours to years, so equilibrium is the default here.
What a slow uptake usually is. Where uptake goes on for days or months, the authors who analyzed it attribute it to transport into the solid rather than to a slow reaction at the site.
Axe and Trivedi (2002) describe the slow uptake of metals by amorphous Al, Fe and Mn oxides as diffusion along the walls of micropores, which hold 40 to 90 % of the sites, every site being at local equilibrium (pp. 259–263).
Spiesz et al. (2012) describe chloride binding during a migration test as a mass transfer toward the binding isotherm,
between 1.0 and 10.5 × 10⁻⁶ s⁻¹, smaller for longer tests (Table 3). Liu et al. (2008) relax each class of sites toward the equilibrium of a surface complexation model at a rate drawn from a lognormal distribution (Eqs. 5 and 6), rates that may stand for diffusion as well as for chemistry (paragraph 23).
Greskowiak et al. (2011) compare that chemical description with a physical one, diffusion into immobile water with equilibrium inside it: the two are equivalent for a linear isotherm (p. 2) and differ when the partition coefficient depends on the composition of the water (Eq. 11, p. 7).
A rate constant fitted on a slow uptake is therefore an apparent parameter. It lumps a transport that a zero-dimensional calculation does not resolve, and it carries the conditions of the experiment it was fitted on.
The model. A family one of whose members a kinetic reaction controls is moved, whole, to the kinetic partition, its free site included. Its site row leaves the minimization and is conserved by the stoichiometry of the reactions, each of which moves one site from one state to another. sorption_rate gives such a reaction the law
a forward and a backward mass action whose ratio is the equilibrium constant of the database: the kinetics of a Langmuir site (Rubin, 1983; Eqs. (66) and (67), p. 1244). Section 7 of Rate laws gives the properties that follow. In a solution whose activities stay constant, a site family starting bare fills as
which the test suite checks against an integrated run, to
Two sites, and several. A fast family at equilibrium beside a slow one is the two-site description: the fast family follows the solution at every instant, the slow one only at the end, and both reach the equilibrium of the whole system, which is checked. A distribution of rates, as Liu et al. (2008) use, can be approximated by several slow families with different constants; the package does not draw the distribution itself.
A slow surface on a host that dissolves or grows. When the budget follows a host (SITES_FOLLOW_HOST) that a rate law controls,
with
What is refused, and what is assumed.
A slow family whose budget follows a host at equilibrium is refused: its sites would appear or vanish with an equilibrium phase while their occupancy is integrated.
An electrostatic family sharing its support with a family left at equilibrium is refused: both would feel one potential, from a charge the equilibrium solve does not hold. On a support of its own, the potential of a slow family is the eliminated one, a function of its own composition.
A state of a slow family that no reaction forms is refused, since its amount would never change.
The activity coefficients of the family, an exchange convention or a surface potential, are in
, hence in the backward term. Thermodynamics fixes the ratio of the two terms, not how a non-ideality is shared between them: this is the choice made.
3. The state and the right-hand side
The state is
the amounts of the equilibrium partition being those of the minimum of the Gibbs energy at the budget and the temperature of the state,
The right-hand side
When the partition may be frozen within a step
When no rate law reads the equilibrium partition, as for a law of Parrott and Killoh (1984) or of Waller (1999), speciation = :frozen). speciated_states certifies the partition at any instant afterwards.
When a rate law reads it, through an activity, a saturation ratio or the amount of an equilibrium species, the partition has to be solved where the integrator evaluates speciation = :rhs). Frozen within a step instead, the rate is constant over the step, and the stiff method integrates the extent explicitly however implicit it is. A rate integrate tells the two cases apart from the rate laws themselves (speciation = :auto), and a trajectory that reaches amounts the system cannot hold is returned as a failure rather than a success.
A semi-adiabatic calorimeter takes the second route whatever its laws read: its temperature is solved with the partition (Section 6).
4. The Jacobian
A stiff method needs
in which only
and the same with
whose value part vanishes and whose dual part is exactly the derivative above, in every direction the dual numbers carry. The Jacobian of
5. The implicit step
The differential system of Sections 3 and 4 lets the integrator choose the step. The implicit step of Leal et al. (2017) makes the step itself one minimization, in which the extents are unknowns beside the amounts:
Backward Euler makes the step unconditionally stable. Under
Two constraints on the products. One reactivity row per reaction (coupling = :reactions) leaves the solids free to rearrange, the assemblage being whatever minimizes coupling = :species),
Choosing the step, and why the certificate decides. For a rate that vanishes at equilibrium the step has two solutions: the equilibrium root, and the composition with the mineral wholly dissolved, which satisfies the element balance and the reactivity row kinetic_step_adaptive estimates the error of a step of
What the certificate of a step proves. It is taken on the augmented problem, whose reactivity rows carry the multipliers that make the stationarity of a kinetically held mineral satisfiable. The same composition tested against the unconstrained equilibrium fails, because a mineral held back by a rate law is supersaturated by construction, which is what being held back means. Two choices of the step are left to the certificate rather than to the user because neither works on every problem: whether the starting guess of a system with solid solutions is equilibrated first (a mixing phase is admitted by a tangent-plane test, which a composition without the phase cannot pass), and whether the kinetic minerals are held in the active set (a mineral exhausted at equilibrium is otherwise removed by the active-set rule, after which nothing enforces its reactivity row). The problem being convex, the certificate decides between the candidates exactly.
6. The calorimeters
The enthalpy of the paste is the sum of the standard enthalpies of its species,
the excess enthalpies of the activity models being left out, as What the package counts as heat states. Under partial equilibrium the kinetic reactions only release ions, the hydrates being precipitated by the minimization, so that the heat of those reactions is not the heat of the paste; the heat is the change of
The isothermal cell
The bath holds the paste at
and requires no integration. Its rate along the trajectory is the derivative of the same function,
evaluated by lifting the partition along heat_flow), and the certified replay (heat_release) computes both from the same partition.
The semi-adiabatic cell
The vessel and what it holds besides the paste, of heat capacity SemiAdiabaticCalorimeter). The enthalpy of the cell,
and the temperature is the root of the energy balance
The temperature is solved for, not integrated. Its derivative with respect to
is the heat capacity of the cell at equilibrium, that of the vessel and of the paste at fixed composition, plus the heat the partition takes up as it shifts with the temperature. It is positive wherever the balance has a unique root, and the temperature is found by Newton's method,
The implicit-function theorem applied to
and the package obtains them, as in Section 4, by a single correction lifted into the dual numbers of the state,
That is the reason for carrying
The heat and the temperature the accessors report are functions of the state: temperature_profile solves the root at each instant, cumulative_heat returns heat_flow its rate,
Without an equilibrium partition
When the kinetic reactions produce the hydrates themselves and no minimization is solved, every amount follows from the extents,
whose right-hand side holds no derivative of a minimization; this formulation integrates the temperature, and an isothermal cell the heat
7. What is checked
The formulation is held to identities that a route sharing nothing with it computes:
The partition carries the integrated budget. At every accepted step
integratereports both. A run in which the minimization never succeeded would otherwise be indistinguishable from a healthy one.The coupling reaches the equilibrium it claims. A solid releasing calcium at a constant rate onto a sorbent (
test/kinetics/test_surface_coupling.jl) holds the site budget at every instant, accounts for every mole released, and satisfies the mass action of the binding reaction on the integrated state; the equilibrium and the coupling along a constant-rate dissolution agree with Reaktoro (Validation against Reaktoro).The Jacobian holds the derivative of the partition. Its columns in
equal computed through the certified solver ( test/kinetics/test_rhs_speciation.jl).The calorimeters conserve energy (
test/kinetics/test_calorimetry.jl): the enthalpy of an adiabatic cell,to within of the heat released; computed with the certified equilibrium differentiated with respect to its temperature, and the same with respect to the budget; the heat of an isothermal cell equals the enthalpy difference of the certified replay, and its rate the rate of the replay; a Rosenbrock, a BDF and an explicit Runge–Kutta method integrate the same temperature, to K, with losses through the walls equal to .
8. Limits
The excess enthalpies of the activity models are not in
, whereas the minimization carries the temperature dependence of the activity coefficients; then holds a term that the excess enthalpies would have made a quadratic form, and its positivity is not guaranteed by the stability of the equilibrium alone. A balance with a non-positive is refused at the evaluation rather than corrected. is smooth only while the assemblage does not change. Where a phase appears or vanishes it has a kink, which the integrator meets as a loss of accuracy and answers by shortening its step; no event is located there. The partition is the answer of the certified solver to its tolerance, and the temperature is solved to
K. The implicit step is first order in time; its accuracy is controlled by the adaptive march, not by the step alone.
See also
Coupling kinetics and equilibrium, the three routes run on one mineral, with their trajectories.
The silicates of a CEM I clinker, hydrating end to end, the partition on a clinker.
A CEM I 52.5 N mortar in a semi-adiabatic calorimeter, inside the kinetics, the semi-adiabatic cell on a cement.
Writing a kinetic model, the syntax of the routes and of the calorimeters.