Writing a kinetic model

Before this page

The tutorial Chemical Kinetics, which runs a complete model end to end; the snippets below are the pieces of that run, shown rather than executed.

A kinetic model is assembled from a small number of objects: rate functions that read the state, rate constants, kinetic reactions that attach a rate to a stoichiometry, a problem that gathers them with an initial state, and optionally a calorimeter. This page describes each of them, after the choice between the two routes that advance a model in time, and ends with the features needed by multi-step hydration schemes.

Which route: two ways to advance in time

There are two, they answer different questions, and choosing wrongly is the most common mistake. Both take the same KineticReaction objects and the same rate laws.

	implicit step — kinetic_step	ODE integration — integrate
what advances	one Gibbs minimization per step, extents included	dn/dt = ν' r(n), handed to SciML
kinetics–equilibrium coupling	strong: the aqueous phase is at equilibrium at the end of the step	with an equilibrium solver, the equilibrium partition is re-speciated once per accepted step and held within it; without one, none
the hydrate assemblage	thermodynamic — whatever minimizes G	thermodynamic for the species of the equilibrium partition (with an equilibrium solver); otherwise imposed by the stoichiometry you wrote
time stepping	Δt is yours	adaptive, with a choice of stiff solvers
overshooting saturation	impossible, at any Δt	controlled by the integrator's tolerances
several reactions per mineral	yes	yes

Use the implicit step when the products should be decided by thermodynamics — a mineral dissolving into a solution, carbonation, an assemblage you do not want to prescribe. It is the step of Leal et al. (2017), the one Reaktoro implements, and it is unconditionally stable (The implicit step); Coupling kinetics and equilibrium runs it against the ODE route on calcite.

Use the ODE route when the reactions themselves are the model — a stoichiometric hydration scheme in the form of Lavergne et al. (2018), where C₃S + H₂O → C-S-H + CH is written out and the products are prescribed, not proposed. Every product then follows the coefficients you gave, several pathways may share a mineral, and the integrator chooses the step.

Prescribed products AND strong coupling: coupling = :species

The third combination, which the ODE route does not offer: a solid-to-solid scheme whose products are imposed by the stoichiometry, solved inside the same Gibbs minimization as the aqueous equilibrium.

coupling = :reactions (the default) puts one constraint per reaction, so the minimization still rearranges the solids and the assemblage stays thermodynamic. coupling = :species puts one per kinetic species, nᵢ = nᵢ(0) + Σⱼ νᵢⱼ Δξⱼ, so every product follows the coefficients you wrote while the aqueous phase still minimizes G under element conservation.

kss = KineticStepSolver(cs, DiluteSolutionModel(), reactions; coupling = :species)

Measured on two C₃A pathways — ettringite and monosulfoaluminate, one thousand seconds, constant rates — against what the stoichiometry demands:

	:reactions	:species
extents Δξ	Δt·M·r, exactly	Δt·r, exactly
C₃A, gypsum, AFt, AFm	rearranged by G; AFm at zero	the stoichiometry, to 8–9 digits
element balance	2×10⁻¹¹ mol	2×10⁻¹⁴ mol
certificate	proved optimal	proved optimal

The pinned species are eliminated rather than constrained, which The implicit step explains.

Use :reactions when the assemblage should come from thermodynamics, the ODE route when no aqueous coupling is needed, and :species when you want a stoichiometric scheme with the pore solution still at equilibrium.

The implicit step, in full

using ChemistryLab, DynamicQuantities, OptimaSolver, OrderedCollections

rxn  = Reaction(OrderedDict(cs["Cal"] => 1.0),
                OrderedDict(cs["Ca+2"] => 1.0, cs["CO3-2"] => 1.0);
                symbol = "calcite dissolution")
rate = KineticFunc((T, P, t, n, lna, n0) -> 1.0e-6, NamedTuple(), u"mol/s")
kr   = KineticReaction(cs, rxn, rate)

kss = KineticStepSolver(cs, DiluteSolutionModel(), [kr])

Δξ  = Ref(Float64[])
st1 = kinetic_step(kss, st0, 100.0u"s"; parameters = Δξ)
Δξ[]        # the reaction extents the step found, in mol

The step is one minimization with the extents among its unknowns and the rate evaluated at the end of the step, written in The implicit step. On calcite its reactivity matrix is K = [−1], so that Δξ = Δt·r: the calcite left after 1000 s at 1 µmol/s is n₀ − kΔt.

Choosing the step length

kinetic_step_adaptive takes one step of Δt and two of Δt/2, compares the extents, and accepts the step when their difference is within reltol of the amount each reaction acts on and the step is certified; it halves the step otherwise, and returns the step taken and a suggestion for the next. Three implicit solves per accepted step is the price.

cert = Ref{Any}(nothing)
st1 = kinetic_step(kss, st0, 1.0e5u"s"; certificate = cert)
cert[].optimal                     # read it before trusting a step much longer than 1/k
st2, Δt_used, Δt_next = kinetic_step_adaptive(kss, st0, 1.0e5u"s")   # one accepted step

A single step far beyond the relaxation time of a law that vanishes at equilibrium can land on a second root, the mineral wholly dissolved, which the certificate refuses and the adaptive march never accepts (The implicit step).

integrate decides by itself whether the partition must be solved inside the right-hand side (speciation = :auto); speciation = :frozen or :rhs forces either, and a trajectory that reaches amounts the system cannot hold is returned with retcode = Unstable. Why freezing the partition under a law that reads it overshoots the equilibrium is in When the partition may be frozen within a step; Coupling kinetics and equilibrium shows the three routes on one mineral.

Solid solutions in the step, and the two settings the certificate decides

A solid solution takes part like any other phase — the products of a kinetic dissolution may be a C-S-H solution, decided by thermodynamics at the end of each step. On C₃S dissolving at 20 µmol/s into a four-member CSHQ solution, the step comes out certified at a stationarity of 2×10⁻¹³, with Δξ = Δt·r exactly, all four end-members present, and ten times the step giving ten times the C-S-H.

Two settings are decided by the certificate rather than by the caller, for the reasons The implicit step gives:

• warm_start (default true) equilibrates the starting guess of a system with solid solutions before the step, leaving the component totals untouched;

• pin_minerals (default :auto) tries the step with and without the kinetic minerals held in the active set and keeps the answer the certificate proves.

Pass certificate = Ref{Any}(nothing) to see the proof. It is taken on the augmented problem, so a kinetically held mineral is not reported as supersaturated.

Rate functions: KineticFunc and StateView

All rate functions are encapsulated as KineticFunc objects. They share the positional six-argument calling convention:

kf(T, P, t, n::StateView, lna::StateView, n_initial::StateView) -> Real [mol/s]

where

Argument	Unit	Description
T	K	Temperature (plain Real or ForwardDiff.Dual)
P	Pa	Pressure
t	s	Current integration time
n	mol	Moles of all species — named access via StateView
lna	—	Log-activities — named access via StateView
n_initial	mol	Initial moles — named access via StateView

What a rate law may depend on

Everything in that table, by species name. A rate law is an ordinary Julia function, so there is no expression language to learn and no restriction on what it may read:

# 1. a constant rate
KineticFunc((T, P, t, n, lna, n0) -> 1.0e-6, NamedTuple(), u"mol/s")

# 2. depending on how much is left — first order in the mineral
KineticFunc((T, P, t, n, lna, n0) -> k * n["C3S"], NamedTuple(), u"mol/s")

# 3. depending on the SOLUTION — pH, or any activity
KineticFunc(
    (T, P, t, n, lna, n0) -> begin
        a_H = exp(lna["H+"])            # `lna` is the natural log of the activity
        k * a_H^0.5
    end, NamedTuple(), u"mol/s")

# 4. depending on time explicitly
KineticFunc((T, P, t, n, lna, n0) -> k * exp(-t / 86400), NamedTuple(), u"mol/s")

# 5. depending on the degree of reaction, since `n0` is there too
KineticFunc(
    (T, P, t, n, lna, n0) -> begin
        α = 1 - n["C3S"] / n0["C3S"]
        k * (1 - α)^(2 / 3)
    end, NamedTuple(), u"mol/s")

A reaction that waits for another species to run out

This is the pattern Lavergne et al. (2018) use for the aluminate sequence, and it needs nothing special: read the species you are waiting on and multiply by a gate.

In that model C₃A takes three successive routes, and which one runs depends on what is left:

stage	while	C₃A reacts with	forming
1	gypsum remains	gypsum	ettringite (AFt)
2	gypsum gone, ettringite remains	the ettringite already formed	monosulfoaluminate (AFm)
3	both gone	water alone	C₃AH₆

Three KineticReaction objects on the same mineral, each gated on the reactant its stage consumes:

# Stage 1 — runs while there is gypsum
r_aft = KineticFunc(
    (T, P, t, n, lna, n0) -> begin
        Gp = n["Gp"]
        Gp <= 0 && return zero(Gp)
        k1 * n["C3A"] * _gate(Gp, n0["Gp"])
    end, NamedTuple(), u"mol/s")

# Stage 2 — takes over once the gypsum is exhausted, and consumes the AFt
r_afm = KineticFunc(
    (T, P, t, n, lna, n0) -> begin
        ett = n["ettringite"]
        ett <= 0 && return zero(ett)
        k2 * n["C3A"] * (1 - _gate(n["Gp"], n0["Gp"])) * _gate(ett, n0["Gp"])
    end, NamedTuple(), u"mol/s")

# Stage 3 — only once BOTH are gone
r_hydrogarnet = KineticFunc(
    (T, P, t, n, lna, n0) -> begin
        k3 * n["C3A"] *
            (1 - _gate(n["Gp"], n0["Gp"])) * (1 - _gate(n["ettringite"], n0["Gp"]))
    end, NamedTuple(), u"mol/s")

# A smooth gate: 1 while the reactant is plentiful, 0 once it is spent.
_gate(x, xref) = tanh(max(x, zero(x)) / (1e-3 * xref))

Note that the gates are written from the amounts, and that stage 2 consumes a species stage 1 produced. Both work because the reaction extents are carried in the state, so a non-kinetic amount moves during the step — see the two rules below.

Two rules make gates like this behave:

• Return a zero of the right type, zero(CH) rather than 0.0. The rate is differentiated for the stiff solver's Jacobian and for parameter sensitivity, and a hard Float64 zero breaks the dual-number chain.

• Prefer a smooth gate where you can, as _gate above does. A min or a hard cut-off is a kink, and an adaptive integrator will shorten its step to crawl over it; tanh costs nothing and integrates cleanly. Keep the hard <= 0 && return zero(x) guard as well — it is what keeps a rate from being computed on a negative amount the integrator proposed and rejected.

A gate on the saturation of the reaction itself — a phase that dissolves and never reprecipitates, or one that only grows — need not be written by hand: dissolution_only and growth_only wrap any law with a factor that is zero on the forbidden side of   and continuous with its slope there (What a rate law may be):

law = first_order_rate(k, cs, rxn, BlaineSurfaceArea(blaine))
kr  = KineticReaction(cs, rxn, dissolution_only(law, cs, rxn))

A gate on a species the reactions consume only works because the reaction extents are carried in the ODE state, which is what lets non-kinetic amounts move during the step. Before that they were pinned at their initial values inside the residual, so a gate never closed and the reaction ran past depletion — while the kinetic mass balance stayed exact and the solver reported success. That failure mode is pinned by a test; see Rate laws that depend on a consumed reactant.

Several reactions on one mineral

Write one KineticReaction per pathway, each with its own rate law, and pass them all. The ODE state holds one entry per unique mineral and the contributions accumulate, so C₃A reacting through both an ettringite and a monosulfoaluminate pathway is two reactions sharing a species — see Multiple kinetic reactions per mineral. Gating one pathway on the exhaustion of the sulfate is how the switch between them is written.

StateView provides O(1) named access to a species vector via a pre-built dictionary (sv["C3S"] — no per-step dict allocation):

index = Dict("C3S" => 1, "C2S" => 2)
data  = [2.5, 1.8]          # moles
sv    = StateView(data, index)

sv["C3S"]          # → 2.5
haskey(sv, "C3S")  # → true

The index dict is built once at KineticsProblem construction; the same dict is shared by all StateView wrappers inside the ODE hot path.

Setting up rate constants

Rate constants are built as NumericFunc objects using the Arrhenius factory:

using ChemistryLab, ForwardDiff

# Arrhenius rate constant for calcite acid dissolution [PalandriKharaka2004; Table 33](@cite): log k at 25 °C and the activation energy, read from the data file
calcite = literature_row("PalandriKharaka2004", "carbonate_rates", "calcite")
k_acid = arrhenius_rate_constant(10.0^calcite.acid_log_k, calcite.acid_E)

k_acid(; T = 298.15)    # → 10^log k, in mol/(m²s)
k_acid(; T = 310.0)     # → higher value at elevated T

ForwardDiff.derivative(T -> k_acid(; T = T), 298.15)  # AD-compatible

Cement clinker hydration: Parrott and Killoh (1984) model

parrott_killoh_avrami is the factory for the (Parrott and Killoh, 1984) kinetic model. It returns a KineticFunc that uses the StateView-based calling convention and is fully AD-compatible.

using ChemistryLab, DynamicQuantities

# Predefined NamedTuple parameters for the four main clinker phases
#   PK84_PARAMS_C3S, PK84_PARAMS_C2S, PK84_PARAMS_C3A, PK84_PARAMS_C4AF
# Each has keys: k₁, n₁, k₂, k₃, n₃, Ea, T_ref (with units)

pk_C3S = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S")    # → KineticFunc

# Water availability limit of [Powers1948](@citet) for w/c = 0.40
WC    = 0.40
α_max = powers_alpha_max(WC)          # 0.952 at w/c = 0.40
pk_C3S_wc = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S"; α_max = α_max)

# Fineness correction — the rate scales as B / 385 m²/kg
pk_C3S_fine = parrott_killoh_avrami(
    PK84_PARAMS_C3S, "C3S"; α_max = α_max, blaine = 380u"m^2/kg",
)

parrott_killoh is deprecated — use parrott_killoh_avrami

ChemistryLab also exports parrott_killoh, a smoothed variant with its own parameter sets PK_PARAMS_*. Its formulas are not those of (Parrott and Killoh, 1984) and its parameters match no published set, so the attribution has been withdrawn. With PK_PARAMS_* the diffusion branch takes over within the first few percent of hydration, and C₃S, C₂S and C₃A then all land on α(7 d) = 0.2386 whatever their K₁ — against the ≈ 0.61 the literature reports for a CEM I at w/c = 0.40. See Two Parrott–Killoh variants. Supplementary cementitious materials do not follow Parrott–Killoh at all: use waller.

Rate evaluation:

index = Dict("C3S" => 1)
n_sv  = StateView([0.5], index)    # 0.5 mol C3S remaining
n0_sv = StateView([1.0], index)    # 1.0 mol initial
lna   = StateView([0.0], index)

r = pk_C3S(293.15, 1.0e5, 0.0, n_sv, lna, n0_sv)   # mol/s

AD smoke-test:

using ForwardDiff

drdT = ForwardDiff.derivative(T -> pk_C3S(T, 1.0e5, 0.0, n_sv, lna, n0_sv), 293.15)
@assert isfinite(drdT) && drdT > 0

Transition-state theory rate model

For models based on solution chemistry (calcite, quartz, …), transition_state builds a multi-mechanism TST KineticFunc:

calcite   = literature_row("PalandriKharaka2004", "carbonate_rates", "calcite")
k_neutral = arrhenius_rate_constant(10.0^calcite.neutral_log_k, calcite.neutral_E)
k_acid    = arrhenius_rate_constant(10.0^calcite.acid_log_k, calcite.acid_E)

surface = BETSurfaceArea(90.0)   # 90 m²/kg

tst = transition_state(
    [
        RateMechanism(k_neutral, 1.0, 1.0),
        RateMechanism(k_acid,    1.0, 1.0, [RateModelCatalyst("H+", 1.0)]),
    ],
    cs,
    cs.dict_reactions["calcite dissolution"],
    surface,
)

# tst is a KineticFunc — same calling convention as parrott_killoh
kr = KineticReaction(cs, cs.dict_reactions["calcite dissolution"], tst)

cs.dict_reactions holds only the kinetic reactions

The lookup above works only if calcite was declared kinetic on the ChemicalSystem. For any other reaction, build the Reaction yourself and pass it directly — dict_reactions is not a catalog of the database.

The transition_state factory captures the ΔₐG⁰ callables from the reaction and recomputes ln K(T) at every ODE step, so the saturation ratio Ω is correct even when the temperature changes (semi-adiabatic calorimetry).

For a single mechanism without catalysts, use first_order_rate as a convenience wrapper.

The mechanisms of a mineral, from Palandri and Kharaka

palandri_kharaka builds the list above from the tables of Palandri and Kharaka (2004) transcribed in data/literature/PalandriKharaka2004.json: each mechanism the report gives the mineral a rate constant for, with its activation energy and its order in H⁺, negative for the base mechanism. palandri_kharaka_minerals lists the minerals, among which quartz and amorphous silica, the hydroxides (brucite, gibbsite), the sulfates (gypsum, anhydrite), hematite and wollastonite.

mech = palandri_kharaka("gibbsite")              # acid, neutral and base
tst = transition_state(mech, cs, rxn_gibbsite, BETSurfaceArea(10.0))

# Quartz is tabulated against its geometric and its BET area: the row says which
qtz = palandri_kharaka("quartz, BET surface area")

A mechanism the report does not give is simply absent, and a rate constant without an activation energy (the neutral mechanism of gypsum) is refused unless assume_Ea supplies one. The carbonate mechanism reads the partial pressure of CO₂ and is built only when asked for, with the gas species that carries it: palandri_kharaka("calcite"; mechanisms = (:acid, :neutral, :carbonate), pco2 = "CO2").

Defining kinetic reactions

KineticReaction associates a Reaction with a KineticFunc:

pk = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S")

# Convenience: look up "C3S" in cs by symbol/formula, build minimal dissolution Reaction
kr = KineticReaction(cs, "C3S", pk)

# Explicit Reaction object (e.g. from cs.dict_reactions)
rxn = cs.dict_reactions["calcite dissolution"]
kr  = KineticReaction(cs, rxn, tst_calcite)

# Reaction-centric (kinetics stored in rxn.properties)
rxn[:rate] = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S")
kr = KineticReaction(cs, rxn)

Reaction-centric workflow

Kinetics can be attached directly to Reaction objects via their properties dict:

Properties table

Key	Type	Required	Description
:rate	KineticFunc or (T,P,t,n,lna,n₀)->Real	✓	Net dissolution rate [mol/s]
:ΔᵣH⁰	AbstractFunc or Number		Reaction enthalpy [J/mol] (thermodynamic convention: negative = exothermic). Computed automatically from species ΔₐH⁰ for balanced reactions; set explicitly for custom species without ΔₐH⁰.

When :rate is a plain callable (not a KineticFunc), it is automatically wrapped in a KineticFunc with empty refs.

Example: C₃S hydration (reaction-centric)

using ChemistryLab, OrdinaryDiffEq, DynamicQuantities

sp(name) = cs[name]

pk_C3S = parrott_killoh_avrami(PK84_PARAMS_C3S, "C3S"; α_max = powers_alpha_max(WC))

rxn_C3S = Reaction(
    OrderedDict(sp("C3S") => 1.0, sp("H2O@") => 3.33),
    OrderedDict(sp("Jennite") => 0.167, sp("Portlandite") => 1.5);
    symbol = "C₃S hydration",
)
rxn_C3S[:rate] = pk_C3S   # heat computed from species ΔₐH⁰ (CEMDATA18)

kp = KineticsProblem(cs, [rxn_C3S, rxn_C2S, rxn_C3A, rxn_C4AF], state0, tspan)

Unit-aware API

All constructors support DynamicQuantities Quantity values. Plain Real → SI assumed.

Constructor	Parameter	Default SI unit	Example
arrhenius_rate_constant	k₀	mol/(m²·s)	0.5u"mmol/(m^2*s)"
arrhenius_rate_constant	Ea	J/mol	62.0u"kJ/mol"
parrott_killoh_avrami params	k₁, k₂, k₃	s⁻¹	1.5u"1/d"
parrott_killoh_avrami params	Ea	J/mol	42.0u"kJ/mol"
parrott_killoh_avrami params	T_ref	K	293.15u"K"
FixedSurfaceArea	A	m²	500.0u"cm^2"
BETSurfaceArea	A_specific	m²/kg	0.09u"m^2/g"
Reaction	:ΔᵣH⁰ (custom species)	J/mol	NumericFunc((T,) -> -36100.0, (:T,), u"J/mol")
KineticsProblem	tspan	s	(0.0u"s", 7.0u"d")
IsothermalCalorimeter	T	K	298.15u"K"
SemiAdiabaticCalorimeter	Cp	J/K	4.0u"kJ/K"
SemiAdiabaticCalorimeter	T_env	K	293.15u"K"
SemiAdiabaticCalorimeter	L	W/K	500.0u"mW/K"
SemiAdiabaticCalorimeter	T0	K	293.15u"K"

Running a simulation

using OrdinaryDiffEq   # activates KineticsOrdinaryDiffEqExt

kp  = KineticsProblem(cs, [kr], state0, (0.0, 7200.0))
sol = integrate(kp)   # default Rodas5P

ks  = KineticsSolver(; ode_solver = Rodas5P(), reltol = 1e-8, abstol = 1e-10)
sol = integrate(kp, ks)

Isothermal calorimetry

cal = IsothermalCalorimeter(298.15u"K")   # T = 25 °C
kp = KineticsProblem(cs, reactions, state0, tspan; calorimeter = cal)
sol = integrate(kp, ks)

t, Q    = cumulative_heat(sol, cal)   # Q(t) [J]
t, qdot = heat_flow(sol, cal)         # q̇(t) [W]

Calorimetry under partial equilibrium

With an equilibrium_solver attached, both calorimeters balance the enthalpy of the whole composition, the partition included, and a semi-adiabatic cell solves its temperature with the partition at every evaluation (The calorimeters). Every species then needs an enthalpy of formation, and a system where one lacks it is refused.

kp  = KineticsProblem(cs, reactions, state0, tspan;
                      calorimeter = IsothermalCalorimeter(293.15u"K"),
                      equilibrium_solver = EquilibriumSolver(cs, model, OptimaOptimizer()))
sol = integrate(kp, ks)

t, Q = cumulative_heat(sol, kp.calorimeter)           # J, H₀ − H at each saved state
t, Q, q̇ = heat_release(sol, kp; times = my_times)     # J and W, from certified states

heat_release replays each instant through the certifying solver and gives the same heat as the run: both solve the same partition.

enthalpy and heat_capacity give the same sums for a single state, and missing_enthalpy lists the species that carry no ΔₐH⁰ and are therefore absent from the balance — worth checking, because a heat curve missing one hydrate is not visibly wrong.

Semi-adiabatic calorimetry (Lavergne et al., 2018)

Without an equilibrium partition the cell integrates its temperature from the heat of the kinetic reactions and the total heat capacity of the paste and the vessel; under partial equilibrium it solves it from the enthalpy of the cell (The calorimeters).

SemiAdiabaticCalorimeter bundles hardware parameters and initial temperature:

using ChemistryLab, DynamicQuantities

# Quadratic heat loss ([Lavergne2018](@cite): φ = a·ΔT + b·ΔT²)
# `Cp` is what the ODE does not compute: the vessel, and anything inert in it.
# The sample's Σᵢ nᵢ Cp°ᵢ(T) is added at every step, so it is not counted here.
cal = SemiAdiabaticCalorimeter(;
    Cp        = 900.0u"J/K",
    T_env     = 293.15u"K",
    heat_loss = ΔT -> 0.3 * ΔT + 0.003 * ΔT^2,
    T0        = 293.15u"K",
)

# Shorthand for linear Newton cooling (L keyword)
cal_lin = SemiAdiabaticCalorimeter(;
    Cp = 4000.0u"J/K", T_env = 293.15u"K", L = 0.5u"W/K", T0 = 293.15u"K",
)

kp = KineticsProblem(cs, reactions, state0, tspan; calorimeter = cal)
sol = integrate(kp, ks)

t, T_profile = temperature_profile(sol, cal)   # T(t) [K]
t, qdot      = heat_flow(sol, cal)             # q̇(t) [W]
t, Q         = cumulative_heat(sol, cal)       # Q(t) [J]

Multiple kinetic reactions per mineral

Following Leal et al. (2017), reactions — not species — carry kinetics. A single mineral can appear in multiple KineticReaction objects (e.g. C₃A via an early ettringite pathway and a late monosulphate pathway). The ODE state contains one entry per unique mineral; contributions accumulate in du[j].

pk_c3a = parrott_killoh_avrami(PK84_PARAMS_C3A, "C3A"; α_max)

kr_C3A_ett  = KineticReaction(cs, rxn_C3A_ettringite,   pk_c3a)
kr_C3A_mono = KineticReaction(cs, rxn_C3A_monosulphate, pk_c3a)

kp = KineticsProblem(
    cs, [kr_C3S, kr_C2S, kr_C3A_ett, kr_C3A_mono, kr_C4AF], state0, tspan,
)
# length(build_u0(kp)) == 4  (one entry per unique mineral: C3S, C2S, C3A, C4AF)

Species classification Leal2017

KineticsProblem partitions the species as Leal et al. (2017) do:

• Kinetic species (tracked in the ODE state u): the species each reaction controls, which is the species the system declares kinetic, otherwise the reaction's first crystal reactant, otherwise its first reactant (see KineticReaction); and, when that species occupies a surface site, every member of its site family (see sorption_rate).

• Equilibrium species: all the others. With an equilibrium_solver they are re-speciated from their element totals; without one, they follow the stoichiometry of the reactions, n = n(0) + νᵀξ.

Parameter sensitivity and optimization

The rate laws themselves are ForwardDiff-clean, and tested to be: no Float64 casts in the evaluation path, branch guards through _primal, and arrhenius_rate_constant differentiable in k₀, Ea and T_ref.

using ForwardDiff

pk(Ea) = parrott_killoh_avrami(merge(PK84_PARAMS_C3S, (Ea = Ea,)), "C3S")
idx = Dict("C3S" => 1)
f(Ea) = pk(Ea)(300.0, 1.0e5, 3600.0, StateView([0.9], idx), StateView([0.0], idx),
    StateView([1.0], idx))
ForwardDiff.derivative(f, 42_000.0)

Differentiating a run

A run whose rate laws, initial amounts, temperature, calorimeter or explicit heat_per_mol carry ForwardDiff dual numbers is integrated on them. The state and every buffer of the right-hand side take the number type of what the problem is given, the parameters a rate closure captures included, and under partial equilibrium the partition is the certified equilibrium of each step, lifted to the duals by the implicit-function theorem at its answer. Nested differentiations are taken one level at a time, so second derivatives are exact as well.

using ChemistryLab, DynamicQuantities, ForwardDiff, OrdinaryDiffEq

subs = build_species(datapath("cemdata18-thermofun.json"))
cs = ChemicalSystem(speciation(subs, ["Cal", "Portlandite"]; aggregate_state = [AS_AQUEOUS]), CEMDATA_PRIMARIES)

# First-order dissolution of calcite, `n(t) = n₀ e^{−kt}`: the derivative with
# respect to `k` is `−t n₀ e^{−kt}`.
function calcite_left(k)
    st = ChemicalState(cs, [zero(k) * u"mol" for _ in cs.species])
    set_quantity!(st, "H2O@", 1.0u"kg")
    set_quantity!(st, "Cal", 0.5u"mol")
    rxn = Reaction([cs["Cal"]], [cs[p] for p in ("Ca+2", "CO3-2")]; symbol = "dissolution")
    rxn[:rate] = KineticFunc(
        (T, P, t, n, lna, n0) -> k * n["Cal"], (T = 298.15u"K", P = 1.0e5u"Pa"), u"mol/s",
    )
    kp = KineticsProblem(cs, [rxn], st, (0.0, 1.0e4); equilibrium_solver = nothing)
    return integrate(kp, KineticsSolver(; reltol = 1.0e-12, abstol = 1.0e-14)).u[end][1]
end
ForwardDiff.derivative(calcite_left, 1.0e-4)     # −1e4 × 0.5 × e⁻¹ = −1839.4

The derivative is that of the numerical solution, the steps of the integrator included: it departs from the sensitivity of the exact trajectory by the integrator's tolerance, which is the tolerance to give a run that is differentiated. The observables read off it are exact derivatives of that solution: heat_flow from the derivative of the interpolant, heat_release from −dH/dt at each certified state.

Two Parrott–Killoh variants

ChemistryLab ships two implementations of the clinker hydration model of Parrott and Killoh (1984). They are not interchangeable, and their parameter sets are not transferable between them.

	parrott_killoh (deprecated)	parrott_killoh_avrami
Nucleation–growth	(K₁/N₁)(1-ξ)^N₁ / (1 + B·ξ^N₃)	(k₁/n₁)(1-ξ)(-ln(1-ξ))^(1-n₁) (Avrami)
Second mechanism	K₂(1-ξ)^N₂	k₂(1-ξ)^(2/3) / (1-(1-ξ)^(1/3)) (Jander)
Third mechanism	3K₃(1-ξ)^(2/3) / (N₃(1-(1-ξ)^(1/3)))	k₃(1-ξ)^n₃ (power law)
Combination	min(max(r_NG, r_I), r_D)	min of all three
Parameters	PK_PARAMS_C3S …, provenance unestablished	PK84_PARAMS_C3S …
Status	deprecated, warns once	use this one

parrott_killoh_avrami is the canonical 1984 formulation, with the parameters of Parrott and Killoh (1984) as reported by Lothenbach and Winnefeld (2006) and used by Lavergne et al. (2018). Use it when you want the α(t) curves of the cement literature. A useful signature of a correct transcription: with those parameters C₂S has no nucleation–growth stage and C₃S no diffusion-controlled stage.

parrott_killoh is kept only so that existing scripts keep running. Its nucleation–growth term carries no Avrami logarithm, K₃ sits where the canonical form has k₂, N₁ = 3.3 is the canonical n₃, and k₃ = 1.1 has no counterpart at all — so the two are different models, not two parameterizations of one. The primary source, Parrott and Killoh (1984), has no DOI and could not be consulted, so the smoothed variant is not attributed to it. With PK_PARAMS_* the diffusion branch takes over at α ≈ 0.003 (C₂S), 0.013 (C₃S) and 0.057 (C₃A); those three then follow the closed form α(t) = α_max·[1 − (1 − √(2·K₃·t / N₃))³], giving α(7 d) = 0.2386 for all three, while C₄AF is limited by its own nucleation branch and reaches only 0.193. Converting the four clinker demos to the canonical variant moved a CEM I paste at w/c = 0.40 from a mean degree of 0.234 to 0.628, an adiabatic temperature rise from 2.0 to 14.2 °C, and the heat released at seven days from 115 to 308 kJ/kg of cement.

pk = parrott_killoh_avrami(
    PK84_PARAMS_C3S, "C3S";
    α_max    = powers_alpha_max(0.5),      # Powers water availability limit
    blaine   = 380u"m^2/kg",               # rate ∝ B/385 m²/kg
    humidity = 0.95,                       # β_h reduction, 0 below 80 % RH
)

The three corrections are exported separately — powers_alpha_max, blaine_factor and humidity_factor — and multiply the rate. humidity also accepts a callable t -> h(t) for a drying history. A fourth, w_c, is the water/cement factor of Parrott and Killoh (1984), pk_wc_factor: it slows a phase once its degree passes 1.333 w/c, or H w/c with the critical degree H that Lothenbach et al. (2008) fit for each clinker phase, passed as H. In a blend, w/c is the water over the cement alone, so the factor is what lets the clinker of a slag cement hydrate further than that of the plain one.

Supplementary cementitious materials do not follow the law of Parrott and Killoh (1984) at all. Their pozzolanic or latent-hydraulic reaction follows waller, a sigmoid in log-time, with WALLER_PARAMS_FLY_ASH or WALLER_PARAMS_SILICA_FUME. No set is shipped for a slag: the 100-day time this package attributed to Waller (1999) is not in that thesis, so a slag's characteristic time is the caller's, with its source, merge(WALLER_PARAMS_FLY_ASH, (τ = τ_slag,)).

Rate laws that depend on a consumed reactant

The ODE state carries the extents of reaction ξ alongside the kinetic moles, integrating dξ/dt = r. The residual therefore reconstructs every species from n = n(0) + νᵀ ξ before evaluating the rates, so a rate closure may read any amount — including species that are not themselves kinetic.

That is what makes reaction sequencing expressible. Sulfate available for ettringite, portlandite available for a pozzolanic reaction, a product inhibiting its own formation: all are written directly.

# C3A reacts with gypsum while sulfate lasts, then by the sulfate-free route.
# "Gp" is consumed but is not a kinetic species — reading it here is correct.
ε = 1.0e-3
sulfated  = (T, P, t, n, lna, n0) -> pk(T, P, t, n, lna, n0) * n["Gp"] / (n["Gp"] + ε)
unsulfated = (T, P, t, n, lna, n0) -> pk(T, P, t, n, lna, n0) * ε / (n["Gp"] + ε)

Keep such a gate smooth. A hard n["Gp"] > 0 ? r : 0 is a discontinuity in the residual, which a stiff solver will either step over or grind against; the x/(x+ε) form above costs nothing and stays differentiable.

Left without a gate, a reaction that consumes a species it does not read runs on after that species is gone, and the reconstruction takes it below zero. The run then comes back with the retcode Unstable and a warning that names the species, the time it runs out and the shortfall, once the shortfall exceeds the relative tolerance of the integration times the most of the species the element totals allow. Only the species whose amount is an amount of substance are judged, a phase, a gas, the solvent: without an equilibrium partition a solute stands for the total of its component, and H+ below zero is hydroxide.

With an equilibrium solver, the non-kinetic partition is piecewise constant

When re-speciation is active the equilibrium partition is owned by the equilibrium solve, which runs once per accepted step by operator splitting. Those amounts are therefore refreshed between steps but held constant within a step, and are not reconstructed from ξ. A gate on such a species still works; it simply lags by at most one step.

Where to go next

The rate laws and the provenance of their parameters are discussed in Rate laws, and every parameter in them, and the coupling with equilibrium in Coupling kinetics and equilibrium. The application pages Hydration kinetics of a CEM I 52.5 R clinker and A complete CEM I 52.5 N, through its pore solution use these objects on complete cements.