Activity models

Before this page

Thermochemistry and Standard states, which define the activity this page 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 five built-in models — DiluteSolutionModel, DaviesActivityModel, HKFActivityModel, SITActivityModel and PitzerActivityModel — differ in both; §5 measures by how much for the first three, and §6a and §6 treat the other two.

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 relative permittivity 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 at 25 °C are derived as well as tabulated: they agree to 0.01 % with the values of the LLNL aqueous model that Parkhurst and Appelo (2013) tabulate, from 0 to 100 °C — compared at three temperatures in What the choice of activity model costs. Helgeson et al. (1981), Table 1, computed from the water properties of the time, gives 0.5091 and 0.3283 at 25 °C.

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  – mol/kg, 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 (1969) 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.

Away from 25 °C the model keeps at its 25 °C value. With temperature_dependent = true, and are recomputed at the temperature and pressure of the solution from the water properties of §1, which the equation of state of water carries to 1000 °C; , fitted to NaCl at 25 °C, has no temperature function here, and a calculation far from 25 °C in a solution concentrated enough for the term to matter should supply its own value of it, a keyword away.

The Davies equation (Davies, 1962) makes the same move with one more simplification. It drops the ion size, replacing   by  , and fixes the linear term once for every ion:

No per-species datum is needed, and every ion of a given charge has the same activity coefficient. The term plays the part of , with a weight proportional to , and its status is the same: a correction fitted to data, not derived from screening.

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. 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. Once the activities of the solutes are chosen, it leaves the solvent none to choose. Integrating it over the solutes' terms gives the osmotic coefficient ,

and that is how the solvent row of the ideal model and of the two extended Debye-Hückel models is built. The ideal model has  , the exact partner of . Davies, whose ions all carry one function of times , has a closed form for (see DaviesActivityModel). HKFActivityModel integrates its , and . Raoult's mole fraction agrees with the ideal to first order, the two differing by , but breaks the relation by  ; a model that corrected its solutes and left its solvent at would set the two halves of one model against each other. TruesdellJonesActivityModel and SITActivityModel (§6a) keep it.

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	SITActivityModel	PitzerActivityModel
solute scale	molarity	molality	molality	molality	molality
		Davies	extended D-H +	D-H +	virial expansion
		osmotic coefficient, closed form	osmotic coefficient	Raoult (§6a)	osmotic coefficient
per-species data	none	none	ion radii (tabulated, overridable)	one parameter per ion pair — caller input	a parameter per ion pair and per triplet — caller input
scalar inputs	none	, ,	, , , ,	and the convention	the shape constants , ,
, dependence	none	on request	, on request	on request; the set is fitted at one temperature	only; the set is fitted at one temperature
returns	for every species	same	same	same	same
useful to	(1 % on a monovalent ion)			to	the range its set was fitted over, a few mol/kg
Gibbs-Duhem consistent	exactly	exactly on ions; not for a neutral with   (§5b)	to (§5)	see §6a	exactly, by construction (§6)

All of them 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. activity_coefficients evaluates the per-ion formula of the first three directly. SIT and Pitzer couple the solutes and have none; their coefficient is read back from the closure, which writes    with built on its own floored amounts, so that    is exact for a trace solute as for a major one.

One limit is common to all of them and is not one of accuracy. A cation cannot be added to a solution without an anion, so a measurement gives the mean coefficient of a neutral combination,  , and the ratios that appear in balanced reactions, never the coefficient of one ion (Richet, 2001; Secs. 11.6b–c, pp. 273–275). How the mean is shared among the ions is the model's, and so is any quantity that rests on one ion alone, the pH, , among them (Richet, 2001; Sec. 12.1b, p. 287).

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 above 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 below a molal, and in a single salt the Gibbs-Duhem relation holds for every model to rounding: each builds its solvent row from its own solutes' terms, and two ions of one charge and one size leave nothing to break it. What breaks it shows in a mixed solution, through the symmetry measured in §5b: ions of different sizes or charges under a linear term, and neutral species carrying a salting-out coefficient the ions' coefficients do not return. The cost of such a break is not a large error in : an equilibrium is set by the values of the potentials, and two water activities within one percent of each other shift by less than  . It is that the solvent and the solutes no longer draw their potentials from one Gibbs energy, so that there is no function whose minimum the equilibrium is; §5b tests this exactly, and the certificate depends on it.

The B-dot model's own approximation — a single charge-weighted mean ion size in its osmotic coefficient — costs nothing in NaCl, whose two ions carry the ion size of the salt, and appears only in a solution whose ions differ in size.

5b. Does a model come from one excess Gibbs energy at all?

§5 measures a Gibbs-Duhem residual along one path. A small residual there is necessary and not sufficient: a model can cancel along a dissolution and fail badly in another direction. The question underneath is whether the activities are the gradient of any function of , and that has an exact test.

If   for some , then second derivatives commute, so

and, being homogeneous of degree zero,

The first is the sharper: it fails exactly when no has these activities for its gradient, whatever might be. Both are read off the Jacobian of the log-activities, the closure activity_model returns, by automatic differentiation — no is constructed, no solve runs, nothing is differenced against a second solve.

What is measured

On H2O@ / Na⁺ / Cl⁻ / Ca²⁺ / SO₄²⁻ at 0.1, 0.12, 0.02 and 0.01 mol, worst relative asymmetry over pairs, ion/ion and solvent/ion separated:

model	ion/ion	solvent/ion
Pitzer (on NaCl)
ideal dilute
Davies (no neutral species)
B-dot, per ion,   (default)
B-dot, common ,
B-dot, per ion,
B-dot, common and
Debye-Hückel limiting law ( ,  )
,   (one GEM-Selektor run (Kulik et al., 2013), not the prescription of Cemdata18 (Lothenbach et al., 2019))

Two readings, and the second is the one that was not expected.

The limiting law is exact. It is the classical result, and it is the anchor: that model does come from one excess energy, to machine precision.

It takes both corrections to break it, and removing either alone does not mend it. Dropping while keeping per-ion radii makes the ion/ion asymmetry worse.

Why, in two lines

Every of §1 depends on composition only through , and, at a fixed amount of water,   . So

and Maxwell demands that not depend on .

• With a common and no , is proportional to exactly, and the ratio is the same for every ion. Hence the zeros in the table.

• A per-ion puts inside the denominator  , where it does not factor out.

• The term is added with the same coefficient to every ion, so it contributes and not . An ion of charge 2 and an ion of charge 1 then disagree by construction.

Neither is a defect of this implementation: both are the published extended form, which every geochemical code uses. §2 already says what is — a deviation function fitted to one salt — and this is the same fact seen from the thermodynamics rather than from the fit. The model is an approximation with a stated domain; outside a fairly dilute solution it does not have a potential behind it, and saying so is not a criticism of the model or of the codes built on it.

What it means for a certificate

equilibrate_certified proves a KKT point of the problem it is handed. When the activities are not the gradient of a single , that problem is not the minimization of a Gibbs energy, and optimal = true says less than its name suggests: a stationary point of the residuals, not a minimum of a potential.

In practice this is the regime every geochemical calculation on extended Debye-Hückel has always worked in, and the table says how far from a potential one is standing. Where the water activity carries the answer — a sealed paste, a concentrated brine — PitzerActivityModel is the model that has one, and §6 is about what that costs.

test/activities.jl pins the zeros and the orders of magnitude above, so a change to any shows up here rather than in a result.

The two lines of the derivation are checked as such rather than only through their consequence. log10_gamma_expression returns the kernel as a symbolic expression, so can be differentiated and compared between ions: with a common and no the two agree to , and with either correction in place they do not. The same expression is asserted to agree with the compiled kernel, which is what makes the formula written in §1 the formula that runs.

6a. SIT: one parameter per pair, and no radius

Between the corrected Debye-Hückel laws above and the virial expansion below sits the Specific ion Interaction Theory (SIT). It keeps the screening term and replaces the empirical deviation term with a sum over pairs:

the sum running over ions of charge opposite to .

Read that against §2. Davies and the form both write the deviation as one number times the ionic strength, which says the correction depends on how much salt there is and not on which salt. SIT says it does depend on which: and are and , of opposite sign. That is the whole content of the model, and it is why a sodium sulfate solution and a sodium chloride solution of the same ionic strength do not have the same activity coefficients.

The is a convention, not a fit. Every published was obtained with it, so changing it makes a compilation inconsistent with the model that produced it; SITActivityModel exposes it as a field only to make that visible.

What it buys, and where it stops

under SIT does not fall forever: the pairwise term is positive for many pairs and eventually overcomes the screening term, so the curve turns around. In sodium chloride at 25 °C, reaches near   and has come back to at  . A one-parameter deviation term cannot do that and match the dilute end too, which is why the range is the test rather than a single point.

The model is a truncation: the pairwise term is the first order of the same virial expansion §6 carries to second order, so at high ionic strength the terms it drops stop being small. The range usually quoted,  – mol/kg (the table of §4), is a rule of thumb from practice, not a bound derived from the model. It also expects genuine ion pairs to appear in the speciation and not inside .

Why it matters here, beyond being one more option

The thermodynamic databases of the NEA reviews and of ANDRA's ThermoChimie are calibrated in SIT. A taken from those and used under Davies is not the constant that was fitted, and the difference is not a rounding: it is the same category of error as using a surface constant fitted with a diffuse layer in a model without one. Reproducing published work from those compilations needs the model they were written in.

What it costs

The water activity here is Raoult's mole fraction, as in TruesdellJonesActivityModel, not an osmotic coefficient: a departure from a full SIT treatment, in which the solvent follows from the terms of the solutes through    (Richet, 2001; Sec. 11.6d, Eq. (11.116), p. 276), and one that does not enter a comparison made at prescribed proton activity.

And an unlisted is taken as zero, which is the literature's convention and never a measurement. missing_epsilon_pairs reports which pairs of a given system are resting on it, because a convention silently applied is indistinguishable from a coefficient somebody determined.

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.

where   is the mass of water in kilograms, the basis of the molalities.

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. The ThermoFun databases read by this package carry none; PHREEQC distributes a set in its pitzer.dat, and the package ships the cement set of Reardon (1990) in data/pitzer-reardon1990.toml, read by build_pitzer_parameters. What the package refuses is to invent a missing one. 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 further terms belong to the model, and both are carried:

• the higher-order electrostatic terms and of Pitzer (1975), for two ions of like sign and unlike charge, Na⁺ with Ca²⁺ in a cement pore solution. They vanish for a symmetrical pair, enter and the osmotic coefficient from the same excess energy, and use Pitzer's function by the Chebyshev approximation PHREEQC uses (Plummer et al., 1988). A set whose were fitted without them is used with etheta = false;

• temperature terms of the coefficients, in the form PHREEQC's databases give them, used when the model is built with temperature_dependent = true. A published set fitted at one temperature, as the shipped one is, carries none, and is used at its own values whatever the temperature.

A set in PHREEQC's format is read by build_pitzer_parameters with format = :phreeqc, from the file the caller has.

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 et al. (1984) fit the 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 set of Reardon (1990) 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.

Where to go next

The chapter continues with Solid solutions, the other place where a mixture departs from ideality. The expressions of this page are evaluated side by side in What the choice of activity model costs, the ion-interaction model is loaded and checked in The Pitzer model, and what it can be used with, and the syntax for choosing a model is in Activity models.