Chemistry that happens on a surface

Before this page

Standard states, Solid solutions and Oxidation state: a site balance is written like the charge row, and site mixing like a solid solution.

Most of this package treats a solid as a bulk phase: so many moles of calcite, with a composition and a chemical potential, and nothing about its shape. That is enough for dissolution and precipitation. It is not enough for the class of phenomena where what matters is the interface: a proton that attaches to an oxide's outer layer, an alkali that a C-S-H holds without ever forming a separate phase, a chloride that a cement binds without becoming Friedel's salt.

This page explains what is added to describe those, in the order the ideas depend on one another. No system is built and nothing is solved here; the syntax and a worked case are elsewhere.

1. The vocabulary, because the words are not interchangeable

A grain of mineral is the support. Some of the atomic groups exposed on its outside can react with dissolved species; those are binding sites. Sites that behave the same way and share one capacity form a site family. A dissolved species that attaches to one becomes a surface complex.

The canonical reaction is not an exotic one — it is an acid-base reaction:

The ≡X is not an element to look up in a periodic table: it means attached to the support. Read the reaction carefully, because two things happen and one does not:

• the charge of the surface changes, and so does the distribution of hydrogen;

• the number of sites does not. One site on the left, one on the right.

That second fact is the whole reason a surface needs a formalism of its own. Together with its mirror image, ≡XOH ⇌ ≡XO⁻ + H⁺, it is what makes an oxide amphoteric: positively charged in acid, negatively charged in base, and neutral somewhere between.

Two more words, because they are easy to conflate. Adsorption is accumulation at an interface; precipitation is the appearance of a new phase with its own stability. They compete for the same dissolved matter, which is exactly why a surface belongs inside the general equilibrium calculation rather than beside it. And monodentate means one site per adsorbed molecule, multidentate several — a distinction that follows from how a site is defined, and cannot be settled by a good curve fit.

2. A site is a conserved quantity that is not an element

Every calculation in this package rests on conservation: whatever redistributes, the atoms are still there. That is the matrix equation

with one row per component and one column per species. Sites obey a law of exactly the same shape — every site is either free or occupied, and the total is fixed while the support is — but a site is not a chemical element, and its row is not an element row.

The package already had one conserved quantity that is not an element: electric charge. It rides in a species' formula as a placeholder symbol, and the matrix assembly turns it into a row without knowing it is special. Sites are carried the same way, each family with its own placeholder. For a single family with a free site and two occupied states, the row reads

with N the family's site budget in moles. Nothing was added to the matrix machinery to obtain it.

Three states of the same site family, drawn schematically. Protonation changes the charge and hydrogen inventory of a surface group, while each state still counts once in the site balance. The fixed shown here assumes a fixed support capacity; §10 treats a capacity that follows the host mineral.

Where N comes from

Three ways, because published data comes in three shapes and converting between them needs a quantity that is usually not measured: a site density per unit area times an area, a capacity per kilogram of dry support times its mass, or a prescribed total. The second needs no area at all, which matters — turning a clay's exchange capacity into the first would mean inventing a BET area to divide by.

3. The activity of something bound to a surface

A chemical potential is always written the same way,

and the modeling is entirely in what a_j means. For a dissolved species it is built from a concentration; for a pure solid it is one; for an end-member of a solid solution it is a mole fraction.

For a species occupying a site it is a site fraction:

its share of the family's budget. The free site is a species like any other and has an activity of its own — and that is the part worth pausing on, because it is where saturation comes from. As sites fill, n of the free site falls, its activity falls with it, and the reaction that consumes it becomes harder to drive. Nothing in the code says "stop at saturation"; saturation is what running out of a reactant looks like.

Reference state and capacity are two different things

The literature distinguishes a reference density Γ°, which fixes where the activity scale has its zero, from a capacity density Γ_C, which is how many sites physically exist (Kulik, 2002). Changing the first means changing the standard potentials to match, and changes no physics. Changing the second changes the physics. Conflating them is the most common way a set of published constants stops meaning what it meant.

4. Langmuir is a consequence, not an ingredient

Here is the part that surprises people who have met adsorption through isotherms. No isotherm is written anywhere in this package.

Take one family, a free site ≡X, and adsorbates j binding one site each. Each binding reaction has its law of mass action,

— the two N cancel, which is why the result does not depend on how the budget was specified. Write β_j = K_j a_j, the drive of adsorbate j at the dissolved activity it actually has. Then n_{≡Xj} = β_j n_{≡X}, and the site balance n_{≡X} + Σ_j n_{≡Xj} = N closes it:

That is the competitive Langmuir isotherm. It was not assumed; it fell out of a conservation law and a mass-action law, which is all the package ever applies.

Read the denominator, because it carries the physics:

• a weak drive, β ≪ 1, leaves the sites mostly free and gives n_j/N ≈ β_j — the linear regime, where sorption looks like a partition coefficient;

• a strong drive, β ≫ 1, gives n_j/N → 1: saturation;

• the sum is shared. Every adsorbate appears in every other's denominator, because they are drawing on one finite capacity. That shared denominator is competition, and it is the thing no single-species isotherm can express.

5. The surface activity coefficient, and why it is absent here

The literature often writes the same physics with the free site eliminated, carrying only the occupied species. That formulation then needs a correction factor to reproduce what the free site was doing, and Kulik (2006) derives it: for monodentate species on one family,

with θ_Σ the occupied fraction. It grows without bound as the sites fill, which is exactly the resistance to further binding that an explicit free site provides by simply becoming scarce.

This package keeps the free site, so γ_L is not applied. Applying it on top would count saturation twice. The identity between the two formulations is worth knowing, and it is verified rather than asserted: at equilibrium, 1/(1 − θ_Σ) equals N divided by the free-site amount, which is the same statement written twice.

What about several sites per molecule?

The site balance holds for any denticity — a bidentate species contributes twice to it. Ideal mixing of occupied and free sites does not: counting the ways a molecule can straddle two neighboring sites is a different combinatorial problem, treated by the quasi-chemical approximation (Kulik, 2006), whose own derivation excludes mixed adsorbates. This release therefore refuses a multidentate species rather than returning a number that looks like an isotherm.

6. What the certificate still means

This package does not only solve; it proves that what it returns is a constrained minimum. A new kind of species could break that proof, so it is worth saying why this one does not.

Ideal site mixing contributes, per family,

and the entropy term Σ n ln(n/Σn) is convex — the same function, and the same argument, as an ideal solid solution. The site balance is linear. So the surface adds a convex term under a linear constraint, and the certificate covers it unchanged.

That is true of this model and is not a general statement about surfaces. Two additions reopen the question rather than inherit the answer. An electrostatic term changes the energy: §8 shows that a constant capacitance keeps it convex, and §9 that a diffuse layer, whose potential depends on the ionic strength of the aqueous phase, leaves no Gibbs energy behind the activities at all. The certificate then reports a self-consistent speciation rather than a minimum, and says so in its scope. An attractive lateral-interaction term, which can make the mixing energy non-convex near half coverage, is not implemented (§12).

7. Cation exchange: the same machinery, counting charge

A permanent-charge clay is not an oxide. Its charge does not come and go with pH — it is built into the mineral, from substitutions in the lattice — and there is no such thing as an unoccupied site: every unit of charge is compensated by some cation, and what varies is which.

That sounds like a different model. It is the same one, with the sites counted differently.

The budget is in charge, and it says so in the formula

Write the exchanger species as Na-X and Ca-X₂. The second carries two of the family's pseudo-elements because a calcium neutralizes two units of charge, and the conservation row that comes out of the matrix assembly is therefore

— the cation exchange capacity, in moles of charge. Nothing was declared to make that happen: the coefficient is in the formula, and the row counts what the formula says.

There is also no free site to serve as the reference of the mixing, so one of the forms is chosen as the reference — the abundant monovalent one, usually.

Two conventions, and no factor between them

Here the two literatures part company, and the package refuses to guess which one a number came from.

convention	activity of an exchanger species	counts
Vanselow		particles
Gaines-Thomas		charges

One calcium and one sodium are one particle each and two charges against one. For a homovalent exchange, Na⁺/K⁺, every is 1 and the two fractions are the same number — which is why an oxide surface never has to choose. For a heterovalent one they diverge, and so do the selectivity coefficients fitted under each.

The literature quotes conversion factors of 2, 3 or 4 between them. Those are trace-composition limits, not constant offsets: the exact relation depends on the exchanger's composition, which is what the calculation is solving for (Marinich et al., 2025). So this package converts nothing implicitly. The convention is part of the declaration, and a constant fitted under one and used under the other is a different model — not a rescaled one.

Both conventions are checked against Reaktoro (Leal et al., 2017), which carries both, on a heterovalent Na/K/Ca exchange, and agree to 4 parts in 10⁹. On that composition they differ from each other by 48 % on the sodium, which is the point.

8. A charged surface, and why it needs no new unknown

Everything above ignored one thing: a surface that binds protons becomes charged, and the work of putting one more charge on an already charged object is not zero. The chemical potential of a member acquires an electrical term,

with the electrostatic potential of the surface plane (V), the charge number of the species, the Faraday constant, the gas constant and the temperature. What relates to the charge the surface carries, per unit area (C/m²), is a closure, and the simplest is a capacitor:   , with the capacitance per unit area (F/m²).

The elimination, and what it buys

The literature presents an electrostatic surface model as one extra unknown per surface with one extra equation. For the diffuse layer that is unavoidable: there depends on the ionic strength through a relation with no closed inverse. For a constant capacitance it is not. The closure inverts, so

where is the area of the surface (m²), the amount of the surface species (mol) and its charge number, so that   is the charge of the whole surface. This is an explicit function of the composition — and a composition-dependent term in a chemical potential is exactly what an activity coefficient is. So the model belongs with the mixing, and the solver needs no new machinery at all.

The certificate survives it, and that is a proof

The electrical work of charging the surface is   per unit area, running over the charge density from zero to . Over the area , and with  , it integrates to

a quadratic form whose Hessian   is positive semi-definite for any positive capacitance. Convex term, convex mixing, linear constraint: the certificate of §6 covers the sum unchanged. Its gradient is , which is how one knows the energy is the right one.

That is this model. It does not transfer: the diffuse layer's Ψ depends on the ionic strength, hence on the aqueous composition, so convexity at fixed ionic strength is not convexity in the composition, and that question has to be reopened rather than inherited.

What it does, and where it stops

A surface already charged resists charging further, so a titration curve flattens: the transitions spread over more pH units than the constants alone would give. On hydrous ferric oxide at pH 5, protonation falls from 0.995 without the term to 0.81 at C = 3 F/m².

Convexity makes the minimum unique, so any failure to find it is numerical. There is one, and its scale is

the potential the surface would reach fully charged. Below about 5 the solve is found to machine precision; above it the Newton loses it. On that same oxide that is C ≳ 3 F/m², which puts the usual oxide range of 1–3 at the edge.

The remedy is the formulation, not the tolerance. Carrying Ψ as an unknown with   as its equation is the same problem — the elimination proved that — but the Newton then controls the potential directly instead of meeting it through a stiff exponential. The extra unknown of the textbooks is a preconditioner.

9. The diffuse layer, and the price of eliminating a potential

§8 removed an unknown by inverting a closure. The diffuse layer is where that trick is usually said to stop, and it is worth being precise about what stops and what does not.

The closure inverts too

Gouy and Chapman's solution of the Poisson-Boltzmann equation beside a flat surface, for a symmetric 1:1 electrolyte, relates the charge the surface carries to the potential it raises and to how well the solution screens it:

with the ionic strength of the bulk solution, the relative permittivity of water at the temperature of the solve unless a value is given (DiffuseLayer), and   kg/m³ the factor that makes it a volumetric concentration. This is transcendental in — and monotone in it, which is a different thing. A monotone relation inverts, and this one inverts in closed form:

asinh is smooth, its derivative is bounded by one, and it costs the solver a single call. So the diffuse layer, like the constant capacitance, is an activity coefficient: ln a_k = ln x_k + z_k ψ̃, with the same convention and a different closure. No unknown is added.

What is actually lost is not the unknown

Two different questions are usually merged into one word, "consistency", and this is the model that separates them.

	is it a gradient?	is that gradient extensive?
ideal site mixing	yes	yes
constant capacitance	yes	no
diffuse layer	no	no

The first column is the one the certificate needs: a minimization has to be minimizing something. Symmetry of the activity Jacobian,  , is exactly the statement that such a something exists.

A constant capacitance passes it — its charging work   is a genuine potential — and fails the second column, because that expression is homogeneous of degree two in the amounts while a Gibbs energy is homogeneous of degree one. That failure is not a defect. The area is a parameter of the problem, fixed from outside like a volume; scaling the amounts without scaling the surface is not the extensive scaling Gibbs-Duhem is about.

A diffuse layer fails the first column, and that one is a defect. depends on , which depends on the aqueous composition, while the activity of an aqueous ion does not depend in return on what is bound to the surface. The coupling is one-way, so the Jacobian is asymmetric, and an asymmetric Jacobian is the Hessian of nothing.

This is a property of the model, not of any implementation of it. Treating the bulk as a reservoir whose ionic strength is a parameter is precisely the approximation that lets a diffuse layer be written without carrying its own inventory of counter-ions — the approximation Dzombak and Morel (1990) make, and the one the default SURFACE block of PHREEQC (Parkhurst and Appelo, 2013) makes. What comes back from such a solve is a self-consistent speciation, mass action and conservation satisfied together. It is not a certified minimum, and the certificate reports it as a self-consistent speciation (scope = :self_consistent).

And a second price, which is the one that bites

Writing a potential as an activity coefficient means the solver reaches it by successive substitution: a composition implies activities, which imply a composition. That iteration converges while the activity moves less than the composition does — while the dimensionless sensitivity

is small. Above roughly five it is not, and what the solve returns is not a second root but a point that violates mass action outright.

Three things follow, and all three are measured rather than argued:

• The physics is never the problem. Charging a surface always opposes further charging, so the equilibrium is unique and stable at every composition. Only the elimination fails.

• A diffuse layer is hard in the middle and easy at the edges. Far from the point of zero charge is large, asinh flattens, and the same system that will not solve at pH 7 solves at pH 4. It is also harder in a dilute background, as , which is the opposite of the usual intuition about difficult chemistry.

• A constant capacitance is judged once. Its sensitivity, , does not depend on the composition at all.

The forecast is electrostatic_stiffness and the threshold ELECTROSTATIC_STIFFNESS_LIMIT; the verdict is the stationarity residual of the certificate, which separates the two regimes by fourteen orders of magnitude.

Which is why is carried as an unknown

The resolution is to stop eliminating it. Writing the potential as an activity coefficient is what creates the fixed point; declaring it an unknown of the solve, with Gouy-Chapman's relation as its equation,

moves the coupling into the outer Newton, which has a Jacobian for it. The activity model is then told its potential rather than computing one, so during the inner solve is a constant and that loop converges in a single pass — exactly as it does for ideal mixing.

This is what solve does by default for a family whose model needs it, and nothing else changes: the residual is already in units of , so it sits in the Newton system unscaled, and the starting guess   is the uncharged surface a solve without electrostatics would return.

The elimination remains reachable, by surface_potential = :eliminated, and it is kept because it is what makes the statement above a measurement: over the same eighteen points, the eliminated route meets the optimality conditions at three and the unknown at all eighteen — in the :self_consistent sense, as the next paragraph says.

None of this repairs the first price. The model is still not the gradient of a Gibbs energy, whichever way the potential is obtained. An unknown makes the system solvable; it does not make an asymmetric Jacobian symmetric, and what comes back is still a self-consistent speciation rather than a certified minimum.

10. A support that appears and disappears

Everything above holds the sorbent still. In a hydrating cement it does not: the C-S-H that would carry the sites precipitates, and a budget posed once is then a statement about a quantity that no longer exists.

A capacity measured per unit mass, or per unit specific area, is linear in the host's amount,

so the whole dependence is one coefficient, , and the site balance becomes

That is equation (30) of Kulik (2002) read as a coefficient. Whether a capacity really is homogeneous of degree one is measured rather than assumed: sites_per_host evaluates it at two scaled host amounts and requires the budget to scale with them — with held fixed, since scaling it too makes a ShrinkingCoreArea ratio equal one everywhere and hides the nonlinearity the probe exists to find.

The free sites belong to the host

How the coupling enters the constraint matrix follows from what the host's database formula contains. A solid such as is described by one formula and one energy, measured on real grains whose surface exposes the hydroxyls that the sites are made of: the free sites are part of the host, and counting them on top of it would count their atoms twice. The constraint matrix therefore loses, in the host's column, times the column of the free site,

which is the bookkeeping of Kulik (2002) for surface groups that belong to their sorbent. On the site row it reads   ; on the element rows it removes from the host the atoms carried by its free sites, so that every atom is counted once; and since the free site is neutral, the charge of the system is conserved while the host dissolves or grows. The subtraction holds in any basis of primaries, the free site or a bare site component carrying the site row.

The free site then carries zero standard energy. An intact grain, every site free, has the composition of the database formula, and it has the database energy only if the free sites add nothing to it. With a fixed budget the free site's energy cancels from every surface reaction and is a gauge; with a budget that follows its host it enters the host's energy, and a value away from zero moves the host's solubility by  , which host_coupling_bias reports and the construction of the system refuses above 0.05 log units.

What the coupling costs, and what it does not

The host's saturation index includes the energy of its surface, which is the thermodynamic statement that a sorbing surface is not the same phase as a bare one; saturation_indices forms the index with the same matrix the solve was constrained with. The certificate is not affected: the coupled matrix is linear, so the feasible set is a polyhedron and the objective is unchanged.

Measured

On hydrous ferric oxide at the weak-site density of Dzombak and Morel (1990),  , titrated by hydrochloric acid until half of it has dissolved, every point certifies, the charge of the system is conserved to mol, and the site total is times the host amount to ; the free site and a bare site component as primaries give the same host amount to . Against PHREEQC, whose SURFACE can be coupled to an EQUILIBRIUM_PHASES mineral, the site total is the declared coefficient times the phase amount to   over five partially dissolved states.

A diffuse layer on a support that grows

A diffuse layer of fixed area on a support that grows would spread a growing charge over an unchanged area, so its charge density, and its potential, would follow the amount of the phase rather than its state. The charge per unit area is a property of the phase: DiffuseLayer can be given a specific area of the host, per mole or per kilogram, instead of an area. When the system is built it becomes an area per mole of sites of the support,   with summed over the families on it, and the area is read from the site amounts at every solve,

the sums over the members of the support, each occupying one site. At the solution   , so   , and is homogeneous of degree zero in the amounts: a function of the occupancy alone. Reading the area from the sites rather than from the host keeps it valid when the host is not a species of the solve at all — a kinetic host, outside the equilibrium partition, whose budget the run carries (partial-equilibrium kinetics).

The eliminated potential then reacts to the composition through the area as well as through the charge,    , which ChemistryLab.electrostatic_stiffness reports. On the weak sites of hydrous ferric oxide in acid NaCl, multiplying every amount by three leaves and unchanged to , where a layer held at its first area raises by 90 % and from 6.33 to 7.61.

11. Three planes of charge

A capacitor and a diffuse layer describe a surface whose charge sits in one plane. On an oxide it does not. The protons and the inner-sphere complexes sit on the surface itself; the ions of the background electrolyte that pair with the charged groups keep part of their water of hydration and sit farther out; an oxyanion bound to the surface by one oxygen keeps its other oxygens, and most of its charge, on the side of the solution. Davis et al. (1978) placed the ion pairs in a plane between the surface and the diffuse layer, the triple-layer model; Westall and Hohl (1980) compared it with the basic Stern model, in which the ion pairs sit where the diffuse layer begins; Hiemstra and Van Riemsdijk (1996) shared the charge of a complex between two planes, the charge distribution model.

ChargePlanes holds the three. A surface has three planes: plane 0, the surface; plane 1, the head of the Stern layer; plane 2, where the diffuse layer begins. Each member of a family puts a charge on each plane , and the charge densities, summed over every family on the surface,

raise three potentials. Gauss's law across each capacitor gives the field between two planes from the charge inside them, and the diffuse layer balances the whole:

with  , and the capacitances of the inner and outer Stern layers, and as in §9. The electrical work of member is the sum over the planes of its charge times their potential,

which for a reaction between members is the of Hiemstra and Van Riemsdijk (1996). The models of the literature differ by where the charges are put and by the capacitances:

model			ion pairs	inner-sphere complex
basic Stern (Westall and Hohl, 1980)			plane 2 (= 1)	plane 0
triple layer (Davis et al., 1978)			plane 1	plane 0
charge distribution (Hiemstra and Van Riemsdijk, 1996)			plane 2	shared between 0 and 1

  merges planes 1 and 2. The surface with one plane and no diffuse layer is ConstantCapacitance (§8), and the one with a diffuse layer alone is DiffuseLayer (§9), which ChargePlanes reproduces when both capacitances are infinite and every charge is on plane 0.

The charge of a group, and the charge it is counted with

The surface groups of a mineral carry fractional charges. On goethite a singly coordinated group and a triply coordinated one each carry half a unit, the share of the bonds of their oxygen that the iron atoms do not satisfy (Hiemstra and Van Riemsdijk, 1996). That half unit raises the potential and belongs on plane 0. A formula, which counts the conservation of charge, carries whole numbers. A member therefore has two charges: , its formula's, and , its planes'. Their difference is the intrinsic charge of the group, the same for every member of a family, so that a reaction between members moves the same charge in both counts; the construction of a family refuses members that do not agree on it. The two conservation statements, and , then differ by the intrinsic charge times the number of sites, which the site balance holds constant.

The potentials are unknowns

As for the diffuse layer (§9), the three potentials are carried as unknowns of the solve, one per plane and per surface, closed by the three equations above; the families of one surface share them, and must agree on its area, its capacitances and its permittivity. The two capacitors alone would be the gradient of a quadratic form, like the constant capacitance; the diffuse layer is not, for the reason §9 gives, so the certificate of a solve on such a surface reports a self-consistent speciation.

Against PHREEQC's CD-MUSIC surface, with the activity coefficients set to one on both sides so that the surface models alone are compared, the goethite of Hiemstra and Van Riemsdijk (1996) in the three placements of the table, 56 points from pH 4 to 11 and 0.001 to 0.1 mol/kg of NaCl, agrees to   in the site fractions and 0.05 mV in the potentials, except at 0.001 mol/kg in the basic Stern model, where the diffuse layer carries 0.18 V at pH 4 and the two codes part by   and 0.3 mV (test/charge_planes.jl).

12. What this page does not cover

Saying what is absent is part of describing what is present.

• No spatial profile of the diffuse layer. Its ions are either left in the solution (§9) or counted in a layer of water of fixed thickness at the average enrichment that balances the surface charge, the Donnan approximation of DonnanLayer; the Poisson–Boltzmann profile across the layer is not computed, and neither is the overlap of two layers in a pore narrower than they are.

• No enthalpy of a surface reaction. The constant of each reaction of a family built by site_family is held at every temperature: the standard energy of a complex follows the aqueous species of its reaction, so that does not move, which is a zero reaction enthalpy. A published constant carries none, and PHREEQC holds those of phreeqc.dat the same way: on the weak sites of hydrous ferric oxide, its site fractions at 50 °C are those at 25 °C at the same proton activity, and the calculation here agrees with them to without a diffuse layer and within   with one (test/surface_temperature.jl). A measured reaction enthalpy has nowhere to go.

• No more than three charge planes, and no plane-resolved surface without a diffuse layer other than the constant capacitance; stacking two electrostatic models is refused rather than summed (§11).

• No multidentate species, per the note above.

• No lateral interactions. Neighbors on a surface affect each other's binding energy, and the models that describe it — Frumkin's interaction term, the quasi-chemical approximation — can make the mixing energy non-convex. Neither is here.

See also

• Surface areas — the syntax, and how a site family is declared.

• Adsorption on a single site family — the smallest complete case, computed, against the closed form above.

• Proving that an answer is the answer — what the certificate checks, and what it does not.

• Solid solutions — the other mixing phase, whose arithmetic is the same and whose closure is not.