The Pitzer model, and what it can be used with
Before this page
Activity models §6.
Activity models §6 derives the ion-interaction model and says why its Gibbs-Duhem consistency is exact rather than approximate. This page loads the shipped parameter set, checks that claim, and establishes what the set may and may not be combined with — which turns out to be the practical question.
using ChemistryLab
using DynamicQuantities
using ForwardDiff
using LinearAlgebra
using Printf
pars = build_pitzer_parameters(datapath("pitzer-reardon1990.toml"))
model = PitzerActivityModel(; parameters = pars)
substances = build_species(datapath("slop98-inorganic-thermofun.json"); verbose = false)
dict = Dict(symbol(s) => s for s in substances)
nacl = ChemicalSystem([dict[s] for s in split("H2O@ Na+ Cl-")], ["H2O@", "Na+", "Cl-"])
# the molar mass the closures themselves use, so a molality formed here is the
# one they form
M_w = ustrip(us"kg/mol", nacl.species[only(nacl.idx_solvent)][:M])
n_w = 1 / M_w
prm(k) = (ΔₐG⁰overRT = zeros(k), T = 298.15, P = 1.0e5, ϵ = 1.0e-30)1. The set says which of its numbers were measured
The tables of Reardon (1990) are those of Harvie et al. (1984), except that the silicate, aluminate and ferrate parameters had no published data at all and were estimated by analogy — he names the analog in each case. Those are exactly the ions a cement assemblage needs, so the set carries the distinction and pitzer_origin reads it back:
for (c, a) in (("Na+", "Cl-"), ("Ca+2", "SO4-2"), ("Ca+2", "OH-"),
("Ca+2", "Al(OH)4-"), ("Na+", "H3SiO4-"), ("Na+", "H2SiO4-2"))
@printf(" %-5s / %-9s β⁰ = %+8.4f β¹ = %+8.4f β² = %+8.2f %s\n",
c, a, pars.beta0[(c, a)], pars.beta1[(c, a)], pars.beta2[(c, a)],
pitzer_origin(pars, c, a))
end Na+ / Cl- β⁰ = +0.0765 β¹ = +0.2664 β² = +0.00 fitted
Ca+2 / SO4-2 β⁰ = +0.2000 β¹ = +3.1973 β² = -54.24 fitted
Ca+2 / OH- β⁰ = -0.1747 β¹ = -0.2303 β² = -5.72 fitted
Ca+2 / Al(OH)4- β⁰ = +0.2145 β¹ = +2.5300 β² = +0.00 estimated:HSO4-
Na+ / H3SiO4- β⁰ = +0.0454 β¹ = +0.3980 β² = +0.00 estimated:HSO4-
Na+ / H2SiO4-2 β⁰ = +0.0196 β¹ = +1.1130 β² = +0.00 estimated:SO4-2A
2. Against measurement, over four decades of molality
Internal consistency cannot tell a correct parameter set from a self-consistent wrong one. Hamer and Wu (1972) can: their Table 16 is a critical compilation of the osmotic and mean activity coefficients of NaCl at 25 °C, and it gives both, so each half of the model is checked separately.
# Hamer & Wu (1972), Table 16, read from the data file: m [mol/kg], φ, γ±.
hw = literature_table("HamerWu1972", "nacl")
HW_ALL = collect(zip(ustrip.(us"mol/kg", hw.m), hw.phi, hw.gamma))
# Ten of its 29 rows, a decade at a time and then every molal, keep the table short.
HW = filter(r -> r[1] in (0.001, 0.01, 0.1, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0), HW_ALL)
bdot = HKFActivityModel()
# γ± of NaCl at m mol/kg, the geometric mean of the coefficients of its two ions
γ_of(mdl, m) = let γ = activity_coefficients(ChemicalState(nacl, [n_w, m, m] .* u"mol"), mdl)
sqrt(γ["Na+"] * γ["Cl-"])
end
# φ from the solvent row of the closure, ln a_w = −M_w Σm φ
lna_pz = activity_model(nacl, model)
φ_of(f, m) = let o = f([n_w, m, m], prm(3)); -o[1] / (M_w * 2m) end
println(" m γ± measured Pitzer dev B-dot dev φ meas φ Pitzer")
for (m, φm, γm) in HW
g1, g2, f1 = γ_of(model, m), γ_of(bdot, m), φ_of(lna_pz, m)
@printf("%7.3f %11.4f %8.4f %6.2f%% %8.4f %6.1f%% %8.4f %9.4f\n",
m, γm, g1, 100 * (g1 - γm) / γm, g2, 100 * (g2 - γm) / γm, φm, f1)
end m γ± measured Pitzer dev B-dot dev φ meas φ Pitzer
0.001 0.9650 0.9650 -0.00% 0.9649 -0.0% 0.9880 0.9884
0.010 0.9030 0.9020 -0.11% 0.9012 -0.2% 0.9680 0.9679
0.100 0.7790 0.7762 -0.35% 0.7718 -0.9% 0.9330 0.9318
0.500 0.6810 0.6786 -0.35% 0.6708 -1.5% 0.9210 0.9208
1.000 0.6570 0.6544 -0.40% 0.6471 -1.5% 0.9360 0.9354
2.000 0.6680 0.6660 -0.31% 0.6562 -1.8% 0.9840 0.9839
3.000 0.7140 0.7116 -0.34% 0.6902 -3.3% 1.0450 1.0454
4.000 0.7830 0.7806 -0.31% 0.7366 -5.9% 1.1160 1.1154
5.000 0.8740 0.8716 -0.27% 0.7922 -9.4% 1.1910 1.1919
6.000 0.9860 0.9866 0.06% 0.8561 -13.2% 1.2700 1.2737The Pitzer column (Pitzer, 1975) follows the measurement to better than half a percent from a millimolal to six molal — through the minimum near 1 mol/kg and the climb back to 0.99 at six molal, neither of which any Debye-Hückel form can produce, since both require a term that grows faster than
The B-dot column (Helgeson, 1969) follows the measurement to 1.5 % up to one molal, with the ion size Helgeson et al. (1981) give NaCl. Its range is nonetheless real — 3 % out at three molal, 13 % at six — and nothing in its output announces the exit.
using Plots
ms = exp10.(range(-3, log10(6.0); length = 120))
p1 = plot(ms, [γ_of(model, m) for m in ms];
xscale = :log10, xlabel = "molality m (mol/kg)", ylabel = "γ±",
label = "Pitzer (Reardon set)", linewidth = 2, color = :steelblue,
title = "NaCl mean activity coefficient at 25 °C", legend = :bottomleft)
plot!(p1, ms, [γ_of(bdot, m) for m in ms];
label = "B-dot", linewidth = 2, color = :firebrick)
A25 = hkf_debye_huckel_params(298.15, 1.0e5).A
plot!(p1, ms, [exp(-A25 * sqrt(m) * log(10)) for m in ms];
label = "Debye-Hückel limiting law", linestyle = :dot, color = :gray, linewidth = 2)
scatter!(p1, [h[1] for h in HW_ALL], [h[3] for h in HW_ALL];
label = "Hamer & Wu (1972), measured", color = :black, markersize = 5)
plot(p1; size = (720, 430), left_margin = 8Plots.mm, bottom_margin = 8Plots.mm)
Read on a logarithmic axis, the three models are one curve below a millimolal — they must be, the limiting law is exact there — and separate irreversibly above a hundredth molal. The measured points sit on the Pitzer curve throughout.
p2 = plot(ms, [100 * (γ_of(model, m) - 1) * 0 for m in ms];
label = "", color = :black, linewidth = 1, linestyle = :dash,
xscale = :log10, xlabel = "molality m (mol/kg)",
ylabel = "deviation from measurement (%)",
title = "Where each model leaves the data", legend = :bottomleft)
scatter!(p2, [h[1] for h in HW], [100 * (γ_of(model, h[1]) - h[3]) / h[3] for h in HW];
label = "Pitzer", color = :steelblue, markersize = 5)
scatter!(p2, [h[1] for h in HW], [100 * (γ_of(bdot, h[1]) - h[3]) / h[3] for h in HW];
label = "B-dot", color = :firebrick, markersize = 5, markershape = :diamond)
hline!(p2, [-1, 1]; label = "± 1 %", color = :seagreen, linestyle = :dot)
plot(p2; size = (720, 430), left_margin = 10Plots.mm, bottom_margin = 8Plots.mm)
This is also, incidentally, a check on the transcription of the parameter file: the Na/Cl coefficients were read off a scanned table, and nothing mistyped reproduces a measured curve over four decades.
3. The dilute limit, exactly
A model that expands around ideality has no freedom as the solution empties:
println(" m ln γ± (Pitzer) limiting law ratio")
for m in (1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3)
lnγ = log(γ_of(model, m))
law = -A25 * sqrt(m) * log(10)
@printf(" %8.1e %14.6f %14.6f %8.4f\n", m, lnγ, law, lnγ / law)
end m ln γ± (Pitzer) limiting law ratio
1.0e-06 -0.001176 -0.001177 0.9987
1.0e-05 -0.003707 -0.003724 0.9957
1.0e-04 -0.011614 -0.011775 0.9864
1.0e-03 -0.035665 -0.037236 0.95784. Gibbs-Duhem, exactly
The reason to prefer this model.
μ_pz = build_potentials(nacl, model)
μ_bd = build_potentials(nacl, HKFActivityModel())
function gd(mu, m, dn)
n0 = [n_w, m, m]
dμ = ForwardDiff.jacobian(nn -> mu(nn, prm(3)), n0) * dn
return abs(sum(n0 .* dμ)) / max(norm(n0 .* abs.(dμ)), 1.0)
end
println(" Pitzer B-dot")
for (name, dn) in ("dissolution " => [0.0, 1.0, 1.0],
"ion exchange " => [0.0, 1.0, -1.0],
"water removal" => [-1.0, 0.0, 0.0])
for m in (0.1, 1.0, 3.0)
@printf(" %s m = %4.1f %10.3e %10.3e\n",
name, m, gd(μ_pz, m, dn), gd(μ_bd, m, dn))
end
end Pitzer B-dot
dissolution m = 0.1 9.890e-17 2.980e-16
dissolution m = 1.0 0.000e+00 0.000e+00
dissolution m = 3.0 0.000e+00 0.000e+00
ion exchange m = 0.1 0.000e+00 0.000e+00
ion exchange m = 1.0 0.000e+00 0.000e+00
ion exchange m = 3.0 0.000e+00 0.000e+00
water removal m = 0.1 4.337e-19 1.301e-18
water removal m = 1.0 1.388e-17 0.000e+00
water removal m = 3.0 2.776e-17 2.776e-17Zero, at machine precision, in every direction and at every molality — for both models, on NaCl. The B-dot model holds the identity too as long as its ions share one ion size, which Na⁺ and Cl⁻ do: both carry that of the salt. It loses it in a solution whose ions differ in size, where its osmotic coefficient rests on one mean ion size and its activity coefficients on one each; the Pitzer construction holds it whatever the composition, because both halves come from one excess Gibbs energy.
5. What the set refuses, and why that is correct
A Pitzer model cannot fall back on ideal behavior for one pair, so completeness is checked when the model meets a species list, and the error names what is missing:
with_br = ChemicalSystem(
[dict[s] for s in split("H2O@ Na+ Cl- Br-")], ["H2O@", "Na+", "Cl-", "Br-"]
)
try
activity_model(with_br, model)
catch err
println(first(split(err.msg, ". ")), ".")
endPitzerActivityModel: the parameter set has no beta0 for 1 cation-anion pair(s) this system contains: Na+/Br-.The same mechanism refuses something more important. A Pitzer set assumes the speciation it was fitted with, and that of Reardon (1990) is fully dissociated: the association of Ca with SO₄, of Na with OH, is inside the
cem = build_species(datapath("cemdata18-thermofun.json"); verbose = false)
pairs_present = [
symbol(s) for s in cem
if symbol(s) in ("Ca(SO4)@", "CaOH+", "Na(SO4)-", "K(SO4)-", "NaOH@", "KOH@",
"Na(CO3)-", "MgSO4@", "HSiO3-", "SiO2@")
]
println("ion pairs and silica species CEMDATA18 carries: ", join(pairs_present, ", "))ion pairs and silica species CEMDATA18 carries: SiO2@, Na(CO3)-, MgSO4@, Na(SO4)-, CaOH+, KOH@, NaOH@, Ca(SO4)@, HSiO3-, K(SO4)-So the shipped set is not a drop-in for a CEMDATA18 cement calculation, and the refusal above is the right outcome rather than a shortcoming of the implementation: a dissociated parameterization combined with an associated speciation is not a defensible calculation in any code. Using this model on a cement pore solution means building the species list the set describes — the free ions and the solids — or obtaining a parameter set fitted for the speciation at hand.
6. Where it does apply
A dissociated alkali-hydroxide-sulfate solution of the kind a cement pore solution approximates, at an ionic strength where the B-dot model is already out of its range:
pore = ChemicalSystem(
[dict[s] for s in split("H2O@ Na+ K+ Ca+2 OH- SO4-2")],
["H2O@", "Na+", "K+", "Ca+2", "OH-", "SO4-2"],
)
# a mixture in the proportions an alkali-rich pore solution reaches
st = ChemicalState(pore, [n_w, 0.30, 0.15, 0.02, 0.40, 0.045] .* u"mol")
γ_pz, γ_bd = activity_coefficients(st, model), activity_coefficients(st, bdot)
println(" γ (Pitzer) γ (B-dot) ratio")
for nm in ("Na+", "K+", "Ca+2", "OH-", "SO4-2")
@printf(" %-6s %10.4f %10.4f %8.3f\n", nm, γ_pz[nm], γ_bd[nm], γ_pz[nm] / γ_bd[nm])
end
@printf("\n a_w: %10.5f %10.5f\n", water_activity(st, model), water_activity(st, bdot)) γ (Pitzer) γ (B-dot) ratio
Na+ 0.6689 0.6659 1.004
K+ 0.6967 0.6795 1.025
Ca+2 0.0643 0.1918 0.335
OH- 0.6503 0.6492 1.002
SO4-2 0.1218 0.2042 0.597
a_w: 0.98559 0.98530The two models disagree by tens of percent on the ions, and the disagreement is not a refinement: at this ionic strength the B-dot term is doing work it was never fitted for, while the Pitzer set was fitted over exactly this range. Which is right is a question for measurement, not for a docstring — but only one of the two is thermodynamically consistent with its own water activity, and §3 says which.
7. An equilibrium: the solubility of halite
Everything above evaluates the model at a given composition. An equilibrium asks the reverse: the composition at which halite, the solid NaCl, stops dissolving. The certified solver recovers the dissolved ions from their potentials by Newton's method on the model's own Jacobian (the inner level), and Hamer and Wu (1972) give the molality of the saturated solution at 25 °C.
sat = ChemicalSystem([dict[s] for s in split("H2O@ Na+ Cl- Hl")], ["H2O@", "Na+", "Cl-"])
st = ChemicalState(sat)
set_quantity!(st, "H2O@", 1.0u"kg")
set_quantity!(st, "Hl", 8.0u"mol") # more halite than water can take
eq, cert = equilibrate_certified(st; model)
n = Dict(symbol(s) => ustrip(us"mol", x) for (s, x) in zip(sat.species, eq.n))
m_sat = n["Na+"] / (n["H2O@"] * M_w)
measured = ustrip(us"mol/kg", literature_value("HamerWu1972", "nacl_saturated_molality"))
@printf("certified: %s saturated molality %.4f mol/kg, measured %.3f (%+.2f %%)\n",
cert.optimal, m_sat, measured, 100 * (m_sat - measured) / measured)certified: true saturated molality 6.1605 mol/kg, measured 6.144 (+0.27 %)The two agree to a quarter of a percent, from the standard Gibbs energy of halite in slop98 (Johnson et al., 1992) and the Na–Cl parameters of the set, neither of them fitted to this number, at a molality six times the 1 mol/kg up to which the Debye–Hückel forms are stated valid.
Where to go next
The derivation is in Activity models, and the API entries are PitzerParameters, build_pitzer_parameters and pitzer_origin. The next group of applications moves to surfaces, beginning with Adsorption on a single site family.