Chemical Kinetics

This tutorial covers the kinetics module of ChemistryLab.jl, which implements mineral dissolution and precipitation kinetics coupled to fast aqueous speciation, following the methodology of (Leal et al., 2017).

Background

The kinetics algorithm solves:

\[\frac{d n_k}{dt} = r_k(t), \quad k \in \text{kinetic minerals}\]

where $r_k$ [mol/s] is the net rate of reaction $k$ (positive = dissolution, negative stoichiometric coefficient for the mineral so dn/dt < 0). At each evaluation of this ODE right-hand side, the aqueous speciation is optionally re-equilibrated, providing accurate activity coefficients for the saturation ratio $\Omega = \text{IAP}/K$.

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

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](@cite))
k_acid = arrhenius_rate_constant(5.012e-1, 14400.0)   # k₀ [mol/(m²s)], Ea [J/mol]

k_acid(; T = 298.15)    # → 0.5012 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: (Parrot and Killoh, 1984) model

parrot_killoh is the factory for the (Parrot 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
#   PK_PARAMS_C3S, PK_PARAMS_C2S, PK_PARAMS_C3A, PK_PARAMS_C4AF
# Each has keys: K₁, N₁, K₂, N₂, K₃, N₃, B, Ea, T_ref (with units)

pk_C3S = parrot_killoh(PK_PARAMS_C3S, "C3S")    # → KineticFunc

# Custom α_max ([Powers1948](@cite) limit for w/c = 0.40)
WC    = 0.40
α_max = min(1.0, WC / 0.42)
pk_C3S_wc = parrot_killoh(PK_PARAMS_C3S, "C3S"; α_max = α_max)

# Custom parameters — accepts Quantity or plain SI Real
pk_ggbs = parrot_killoh(
    (K₁=0.15u"1/d", N₁=2.0, K₂=0.003u"1/d", N₂=2.0,
     K₃=0.0015u"1/d", N₃=3.5, B=0.2, Ea=46_000.0u"J/mol", T_ref=293.15u"K"),
    "GGBS";
    α_max = 0.9,
)

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:

k_neutral = arrhenius_rate_constant(1.549e-6, 23500.0)
k_acid    = arrhenius_rate_constant(5.012e-1, 14400.0)

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 parrot_killoh
kr = KineticReaction(cs, cs.dict_reactions["calcite dissolution"], tst)

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.

Defining kinetic reactions

KineticReaction associates a Reaction with a KineticFunc:

pk = parrot_killoh(PK_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] = parrot_killoh(PK_PARAMS_C3S, "C3S")
kr = KineticReaction(cs, rxn)

Reaction-centric workflow

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

Properties table

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 = parrot_killoh(PK_PARAMS_C3S, "C3S"; α_max = min(1.0, WC / 0.42))

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.

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]

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

The semi-adiabatic calorimeter solves:

\[\frac{dT}{dt} = \frac{\dot{q}(t) - \varphi(T(t) - T_{\rm env})}{C_p + \sum_i n_i C^\circ_{p,i}(T)}\]

The denominator uses the variable total heat capacity Cp_total = Cp + Σᵢ nᵢ Cp°ᵢ(T), where Cp°ᵢ(T) are the molar heat capacities from the thermodynamic database ((Lavergne et al., 2018)).

SemiAdiabaticCalorimeter bundles hardware parameters and initial temperature:

using ChemistryLab, DynamicQuantities

# Quadratic heat loss ([Lavergne2018](@cite): φ = a·ΔT + b·ΔT²)
cal = SemiAdiabaticCalorimeter(;
    Cp        = (1.0 * 800.0 + WC * 4186.0 + 1.0 * 900.0) * u"J/K",   # ≈ 3449 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]

Full workflow: OPC clinker hydration with KineticsProblem

Complete chain — ChemicalSystem → ChemicalState → KineticReaction → KineticsProblem → integrate:

using ChemistryLab, OrdinaryDiffEq, DynamicQuantities, Printf

# ── 1. ChemicalSystem from CEMDATA18 ────────────────────────────────────────
DATA_FILE = joinpath(pkgdir(ChemistryLab), "data", "cemdata18-thermofun.json")
substances = build_species(DATA_FILE)

input_species = split(
    "C3S C2S C3A C4AF " *
    "Portlandite Jennite ettringite monosulphate12 C3AH6 C3FH6 H2O@",
)
species = speciation(substances, input_species; aggregate_state = [AS_AQUEOUS])
cs = ChemicalSystem(species, CEMDATA_PRIMARIES)

# ── 2. Initial state: 1 kg OPC (CEM I 52.5 R), w/c = 0.40 ─────────────────
WC          = 0.40
COMPOSITION = (C3S=0.619, C2S=0.165, C3A=0.080, C4AF=0.087)

state0 = ChemicalState(cs)
for (name, frac) in pairs(COMPOSITION)
    set_quantity!(state0, string(name), frac * u"kg")
end
set_quantity!(state0, "H2O@", WC * u"kg")

# ── 3. [ParrotKilloh1984](@cite) rate functions with Powers α_max ───────────────────────
α_max   = min(1.0, WC / 0.42)
pk_C3S  = parrot_killoh(PK_PARAMS_C3S,  "C3S";  α_max)
pk_C2S  = parrot_killoh(PK_PARAMS_C2S,  "C2S";  α_max)
pk_C3A  = parrot_killoh(PK_PARAMS_C3A,  "C3A";  α_max)
pk_C4AF = parrot_killoh(PK_PARAMS_C4AF, "C4AF"; α_max)

# ── 4. Kinetic reactions (reaction-centric) ─────────────────────────────────
# Reactions follow [LothenbachWinnefeld2006](@cite) — Jennite = Ca₉Si₆O₁₈(OH)₆·8H₂O
# Balanced hydration reactions — ΔᵣH⁰ computed from species ΔₐH⁰.
sp(name) = cs[name]

rxn_C3S = Reaction(
    OrderedDict(sp("C3S") => 1.0, sp("H2O@") => 103/30),
    OrderedDict(sp("Jennite") => 1.0, sp("Portlandite") => 4/3);
    symbol = "C₃S hydration",
)
rxn_C3S[:rate] = pk_C3S

rxn_C2S = Reaction(
    OrderedDict(sp("C2S") => 1.0, sp("H2O@") => 73/30),
    OrderedDict(sp("Jennite") => 1.0, sp("Portlandite") => 1/3);
    symbol = "C₂S hydration",
)
rxn_C2S[:rate] = pk_C2S

rxn_C3A = Reaction(
    OrderedDict(sp("C3A") => 1.0, sp("H2O@") => 6.0),
    OrderedDict(sp("C3AH6") => 1.0);
    symbol = "C₃A hydration",
)
rxn_C3A[:rate] = pk_C3A

rxn_C4AF = Reaction(
    OrderedDict(sp("C4AF") => 1.0, sp("Portlandite") => 2.0, sp("H2O@") => 10.0),
    OrderedDict(sp("C3AH6") => 1.0, sp("C3FH6") => 1.0);
    symbol = "C₄AF hydration",
)
rxn_C4AF[:rate] = pk_C4AF

# ── 5. Problem + semi-adiabatic calorimeter ─────────────────────────────────
cal = SemiAdiabaticCalorimeter(;
    Cp        = (1.0 * 800.0 + WC * 4186.0 + 1.0 * 900.0) * u"J/K",
    T_env     = 293.15u"K",
    heat_loss = ΔT -> 0.30 * ΔT + 0.003 * ΔT^2,
    T0        = 293.15u"K",
)

kp = KineticsProblem(
    cs, [rxn_C3S, rxn_C2S, rxn_C3A, rxn_C4AF], state0, (0.0, 7.0 * 86400.0);
    calorimeter = cal,
    equilibrium_solver = nothing,
)

# ── 6. Integrate and post-process ───────────────────────────────────────────
ks  = KineticsSolver(; ode_solver = Rodas5P(), reltol = 1e-6, abstol = 1e-9)
sol = integrate(kp, ks)

_, T_vec = temperature_profile(sol, cal)
_, Q_vec = cumulative_heat(sol, cal)

n0_kin = [sol.prob.p.n_initial_full[i] for i in kp.idx_kinetic]
n_kin  = [[u[i] for u in sol.u] for i in eachindex(n0_kin)]

function phase_alpha(cs, kp, n0_kin, n_kin, name)
    sp_idx = findfirst(s -> ChemistryLab.symbol(s) == name, cs.species)
    pos    = findfirst(==(sp_idx), kp.idx_kinetic)
    isnothing(pos) && return fill(NaN, length(sol.t))
    return 1.0 .- n_kin[pos] ./ n0_kin[pos]
end

α_C3S = phase_alpha(cs, kp, n0_kin, n_kin, "C3S")

@printf "ΔT_max = %.2f °C   Q_7d = %.1f kJ/kg   α(C3S) = %.4f\n" \
    maximum(T_vec .- 273.15) - 20.0  Q_vec[end]/1000  α_C3S[end]

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 = parrot_killoh(PK_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)

Parameter sensitivity and optimization

The entire chain is ForwardDiff-compatible:

using ForwardDiff

function n_C3S_final(K₁)
    params = merge(PK_PARAMS_C3S, (K₁ = K₁,))   # override K₁ (plain Real, SI = 1/s)
    pk  = parrot_killoh(params, "C3S")
    kr_ = KineticReaction(cs, "C3S", pk)
    kp_ = KineticsProblem(cs, [kr_], state0, (0.0, 86400.0))
    sol = integrate(kp_, ks)
    return sol.u[end][1]
end

dn_dK₁ = ForwardDiff.derivative(n_C3S_final, safe_ustrip(us"1/s", PK_PARAMS_C3S.K₁))