Coupling kinetics and equilibrium

Some reactions are fast enough to be treated as instantaneous, others are not. A hydrating cement paste has both: aqueous speciation — protonation, complexation, water autoprotolysis — reaches equilibrium in microseconds, while alite dissolves over days. Integrating everything with rate laws would demand kinetic parameters nobody measures for the fast reactions and would force the integrator down to their timescale. Equilibrating everything would dissolve the clinker instantly.

The way out is (Leal et al., 2017)'s partition, which is also what Reaktoro implements. This page derives it; The hydrating paste, end to end puts it to work.

The partition

Split the species into two sets:

• the kinetic partition — phases whose transformation is rate-limited (clinker, most minerals);

• the equilibrium partition — everything assumed instantaneously equilibrated (the aqueous species, and any phase you accept as equilibrated).

Each kinetic reaction has a stoichiometric row that touches both, so the full stoichiometric matrix splits column-wise as    .

Why the state is and not

The obvious choice — integrate every species — does not work, and the reason is worth stating because it explains the whole construction.

Take  : this advances the equilibrium species along the kinetic reactions without re-equilibrating them. After one step the aqueous phase is no longer at equilibrium, so the assumption that justified the partition has been abandoned.

Re-equilibrating after each step does not rescue it either, because the natural constraint for the equilibrium problem is not a composition — it is a set of element totals. Along the way an individual species genuinely wants to go negative: the dissolution reactions generated from the nullspace are written in , and a cement paste contains no acid. What has physical meaning is not "how much was consumed" but "how much hydrogen is in the equilibrium partition", and it is the minimizer, not the caller, that redistributes those elements over a feasible composition.

So the state carries the element amounts of the equilibrium partition,

where is the conservation matrix restricted to that partition. These are conserved by every fast reaction by construction, and only the kinetic reactions move them.

Aₑ must be the matrix the solve is posed on

The formula matrix over the canonical elements and the matrix over the system's primary species are different matrices. Building Aₑ from one while the minimization is posed on the other makes every equilibrium solve infeasible — in one intermediate version of this package, 89 failures out of 89 steps. Aₑ is taken from the equilibrium sub-system, which inherits the parent system's primaries.

The system

With that state, the coupled problem is ((Leal et al., 2017), Eqs. 54–65):

closed by the equilibrium map

Two ODEs and one constrained minimization. The rates depend on the full composition — a dissolution rate needs the saturation index, hence the aqueous activities — so feeds back into the right-hand side.

Solving it: operator splitting

is not evaluated inside the ODE right-hand side. It is applied once per accepted step, as a DiscreteCallback.

This is a deliberate choice, and both reasons matter:

• a stiff integrator evaluates the right-hand side many times per step and differentiates it to build a Jacobian. An optimization solve in there makes the cost unpredictable, and the Jacobian would be taken through an active-set-dependent map;

• with the solve outside, the residual and its Jacobian see the same speciation, so the Jacobian is consistent with the model actually being integrated.

The initial state is equilibrated before the first step, so the trajectory starts on the constraint manifold rather than drifting onto it.

    u = (bₑ, nₖ)  ──▶  ODE step  ──▶  accepted?  ──▶  respeciate!  ──▶  u
                          ▲                              │
                          └───── frozen nₑ ◀─────────────┘

What to check when it runs

Two diagnostics tell you whether the coupling is doing its job:

check	meaning	healthy value
‖Aₑ nₑ − bₑ‖∞	the speciation carries the integrated element totals	~10⁻⁸ mol
re-speciation failure count	steps that kept a frozen composition	0

The second is reported automatically at the end of integrate. It exists because an earlier version of this package swallowed solve failures in a bare catch: a run in which re-speciation never once succeeded looked exactly like a healthy one. Silence is not success — if the count is nonzero, the trajectory is not the one the equations describe.

Validation status

Both halves are checked against Reaktoro: the equilibrium species by species, and the coupling along a constant-rate dissolution trajectory, where the two agree to 4.3 % or better — see Validation against Reaktoro.

See also

• The hydrating paste, end to end — a worked example

• Validation against Reaktoro — what agrees and what does not

• Chemical Kinetics — rate laws and the KineticFunc interface