Activity models

Every equilibrium this package computes rests on

and on nothing else about the solution. The standard potential comes from the thermodynamic database; the activity comes from the activity model, which is therefore the single place where a real solution stops being ideal. This page is about what those models are, where they come from, what they need, and where they stop being true.

A model has to answer two questions, not one:

1. the solute activity  , which sets every solubility and every saturation index;

2. the solvent activity , which is what a hydrate reaction consumes, and which in a cement paste short of mixing water is the quantity that decides how far the reaction can go.

The three built-in models — DiluteSolutionModel, DaviesActivityModel, HKFActivityModel — differ in both, and §5 measures by how much.

1. Where the comes from

An ion in a solution of ions is not in the same state as one alone at the same molality, and the reason is electrostatic. Around a cation, anions are slightly more likely to be found than cations — not because of any structure, only because their energy there is lower — so the ion sits at the center of a diffuse ionic cloud of opposite net charge. That cloud screens the ion's own field and lowers its energy, which is why activity coefficients of ions are below 1.

The classical treatment takes three steps, and the assumptions it makes on the way are exactly the assumptions that later fail:

1. Poisson's equation relates the mean potential to the mean charge density,   , with the solvent entering only as a continuum of bulk permittivity .

2. Boltzmann statistics give the local ion densities,   , treating the ions as independent in the mean field of all the others.

3. Linearization,  , turns the pair into  , whose solution is a screened Coulomb potential with a single length scale

the Debye length. Charging the central ion reversibly inside its own cloud gives the excess chemical potential, and hence

exact as  . The is not fitted: it is , and  .

Two ions cannot approach closer than the sum of their radii, so the cloud is excluded from a shell of radius around the central ion. Carrying that through gives the extended law, which is what HKFActivityModel implements:

and are properties of water, not fitting constants

Collecting the constants of step 3 gives, with in g/cm³ and the dielectric constant of water,

which is what hkf_debye_huckel_params evaluates from this package's own equation of state for water. So the   and   that the models carry as defaults are derived, not adopted, and they agree with (Helgeson et al., 1981) Table 1 — evaluated at three temperatures in What the choice of activity model costs.

Both rise with temperature, because water's dielectric constant falls faster than rises: hot water screens worse, so the same ionic strength costs more.

The screening length is the size of a gel pore

The Debye length is   in ångström when is in Å⁻¹(kg/mol)^½. At 25 °C it runs from 9.6 nm at   mol/kg to 0.30 nm at 1 mol/kg, passing 0.55 nm at   mol/kg — the values are tabulated in What the choice of activity model costs.

A cement pore solution sits around  –, so its screening length is a few ångström — the thickness of two or three water molecules. That is the same scale as the water films in the gel pores of C-S-H, and both assumptions of step 1 and step 2 above are strained there: a continuum of bulk permittivity, and ions independent in a mean field. Nothing in the formulas announces it. It is the physical reason to treat an extended Debye-Hückel model as a correlation valid in bulk solution rather than as a theory of confined water.

2. The term is a deviation function, not a physical term

Measured activity coefficients turn back upwards at high ionic strength, which no screening argument produces. The B-dot model adds a linear term for it,

Its status is that of an empirical correlation. (Anderson and Crerar, 1993) (§17.7.1, pp. 445–446) record that Helgeson defined as a deviation function: the difference between the observed activity coefficient of an electrolyte — NaCl — and what the extended Debye-Hückel expression predicts for it. So it carries short-range ion-solvent and ion-ion interaction and whatever the first two terms failed to capture, together, in one number fitted to one salt. They add that (Helgeson et al., 1981) later split it into a hydration term from the Born equation and a residual short-range term.

That is why this model has a ceiling rather than an asymptote, and why the ceiling is quoted vaguely as "about a molal": the term is not wrong so much as it is standing in for physics it does not contain.

Neutral species get the Setschenow form,  : water engaged in the solvation shells of ions is water unavailable to solvate a neutral molecule, so its activity rises with ionic strength and its solubility falls. This is salting out, and CO₂(aq) is the case that matters for carbonation.

3. The water activity, and why it cannot be assumed separately

The solvent is not a solute and its activity is not obtained by the same formula. Two routes exist in this package.

Raoult —  , the mole fraction — is what DiluteSolutionModel and DaviesActivityModel use. It counts molecules and knows nothing about what they are.

The osmotic coefficient is what HKFActivityModel uses:

with obtained by integrating the Gibbs-Duhem relation over the same , and that produced the . That is the whole point of it. At constant and ,

which is not an optional refinement: it is the statement that the solvent and the solutes are parts of one thermodynamic system. A model that corrects its solutes and leaves its solvent ideal violates it by construction, and §5 measures by how much.

The one approximation in the B-dot route is that uses a single charge-weighted mean radius,  , where the use per-ion radii. §5 measures that too.

4. Inputs and outputs, model by model

What each model needs and what it returns — the full parameter tables, with the provenance of every default, are in the docstrings (HKFActivityModel, DaviesActivityModel); this is the summary that lets you choose.

	DiluteSolutionModel	DaviesActivityModel	HKFActivityModel	PitzerActivityModel
solute scale	molarity	molality	molality	molality
		Davies	extended D-H +	virial expansion
	Raoult	Raoult	osmotic coefficient	osmotic coefficient
per-species data	none	none	ion radii (tabulated, overridable)	a parameter per ion pair and per triplet — caller input
scalar inputs	none	, ,	, , , ,	the shape constants , ,
, dependence	none	on request	, on request	only; the set is fitted at one temperature
returns	for every species	same	same	same
useful to				the range its set was fitted over, a few mol/kg
Gibbs-Duhem consistent	approximately	no (§5)	to (§5)	exactly, by construction (§6)

All three return the same object — a vector of indexed like cs.species, covering solutes, solvent, pure crystals (), gases and solid-solution end-members — so they are interchangeable at every call site, and concentration_scale tells the accessors which convention was used.

5. What the difference is worth, measured

The models are compared, on an imposed NaCl composition and without solving anything, in What the choice of activity model costs. Three results from that page belong here, because they are about the theory rather than about the numbers:

The ideal model is already several percent off at a millimolal. Ideality is not a safe default that degrades gracefully; it is exact only in a limit.

Davies and the B-dot model part company around a tenth molal, and by 3 mol/kg Davies returns   while the B-dot model is still below 1. That is the term of §2 taking over from the screening term — the ceiling of a deviation function arriving, visible in the numbers.

The values of barely separate at all — Raoult and the osmotic route differ by a few parts in a thousand even at 3 mol/kg — while their derivatives differ by four orders of magnitude. Measured as the Gibbs-Duhem residual   along a dissolution at 1 mol/kg:   for Davies against   for the B-dot model, with the ideal model at   in between.

The third result follows from §3: Davies is less thermodynamically consistent than assuming ideality. Correcting the solutes while leaving the solvent at   sets the two halves of one model against each other, and a model can be more wrong for being partly corrected. Since equilibrium is set by derivatives and not by values, a disagreement invisible in is decisive in .

The B-dot model's own defect — the single charge-weighted mean radius in its osmotic coefficient — appears only along a composition change that holds and fixed, where it measures a few parts in a thousand. That is the direction test/activities.jl uses, which is why its tolerance is 5e-3 rather than solver tolerance.

6. The ion-interaction model: a different kind of object

Everything above is a corrected Debye-Hückel law: one screening term derived from electrostatics, one size correction, and one empirical term standing in for the rest. PitzerActivityModel is not that. It starts from the excess Gibbs energy and expands it as a virial series in the molalities, exactly as the pressure of a non-ideal gas is expanded in its density — (Anderson and Crerar, 1993) (§17.8) derive it as a cluster expansion with osmotic pressure in place of pressure, which is where the analogy is exact:

the first term in a virial equation always represents ideal behavior; in the second term represents the non-ideal contribution from pairwise interactions of molecules; gives the interactions of triples.

Hence the shape of the parameter set: one coefficient per ion pair, one per triplet, and no per-species radius at all.

The long-range term keeps the Debye-Hückel physics — it must, since the limiting law is exact — in the form

and the pair term carries the ionic-strength dependence that a single constant cannot:

is the same Debye-Hückel slope as §1, on the osmotic basis:  , so it comes from the water model rather than from a second implementation. , and are not fitted — they fix the functional form, and a published table is only valid with the values it was fitted against.

Why the whole construction matters here

and are both partial derivatives of one function. So the Gibbs-Duhem relation between the solutes and the solvent is an identity of the algebra, not an approximation — and that is the defect §3 measured in the other two models. It is verified rather than asserted: with the derivative taken analytically,   comes out exactly zero for this model along three composition directions at 0.1, 1 and 3 mol/kg, where the B-dot model's residual is not small.

What it costs

The parameters are caller input, and there is no way around it: no thermodynamic database ships them, and this package refuses to invent them. A PitzerParameters takes every table as a keyword without a default, and completeness is checked against the species list rather than against the set alone — every cation-anion pair present must have a , and the error names those that do not.

Two limitations belong here rather than in a footnote:

• the higher-order electrostatic terms , are not implemented. They vanish identically for a symmetrical pair, so a single 1-1 or 2-2 electrolyte is unaffected; in a mixture of Na⁺ with Ca²⁺ — a cement pore solution — they are a real omission;

• a published set is fitted at one temperature, and nothing here extrapolates the interaction parameters away from it.

The speciation a set assumes is part of the set

This one is easy to miss and changes results. A Pitzer set absorbs ion association into its coefficients: Harvie, Møller and Weare fit Ca–SO₄ interaction rather than postulating a CaSO₄⁰ complex. A species list that carries both the free ions and the ion pairs therefore counts the same association twice, and CEMDATA18 does carry them — Ca(SO4)@, CaOH+, Na(SO4)-, NaOH@ among others. It also names the silica species differently (HSiO3-, SiO2@) from the set's H3SiO4- and H2SiO4-2.

So the shipped Reardon set does not drop into a CEMDATA18 cement calculation, and the completeness check refuses such a system rather than returning a number. That refusal is the correct outcome and not a limitation of the implementation: combining a dissociated parameterization with an associated speciation is not a defensible calculation in any code.

7. Outside the domain

None of this survives arbitrary concentration, and the failure is silent: the formulas go on returning finite, plausible numbers. Two guards exist rather than one warning.

solvent_fraction reports the mole fraction of water within the aqueous phase, and SOLVENT_FRACTION_FLOOR is where equilibrate_certified refuses to call the answer a solution at all. The regime it catches is real: a cement paste below  , on the species lists used here, drives the free water to the solver's floor and the ionic strength to hundreds of mol/kg — a composition at which every molality is per kilogram of a solvent that is no longer there. See the w/c example and the water budget.

A Pitzer model does not rescue that regime. Its fitted range is a few mol/kg, not hundreds, and no activity model repairs a composition that has left the physical picture of a solution.