Self-desiccation: where Powers' 0.42 comes from
Before this page
The tutorial Chemical Equilibrium, and The water budget of a hydrating paste for the theory this page applies.
A sealed cement paste stops hydrating before it runs out of cement. The rule of thumb of Powers (1948) caps the degree of hydration at
and powers_alpha_max supplies it to the rate laws. The coefficient is empirical, and this page takes it apart: it writes the water budget that produces it, fills each term from an independent source — the chemistry from a Gibbs energy minimization, the arrest saturation from a published desorption isotherm — and reports what the two together imply.
Nothing below is fitted to 0.42. The companion script runs the whole calculation outside the documentation.
Read this first: the Kelvin term does not arrest hydration
It is tempting to expect the arrest from thermodynamics — water in a fine pore is held at a reduced activity, so hydrates that consume water become less stable. The magnitude is wrong by two orders of magnitude. At
A real paste stops at 75–80 % relative humidity because transport and nucleation stop. So the humidity belongs in the rate law, through humidity_factor and PoreHumidity, and CapillaryWater is what makes the water activity of the equilibrium state mean the pore water rather than the mole fraction. This page uses the thermodynamics for the water budget and the humidity only as the arrest criterion.
1. The water budget
Take one gram of cement and
Some of it is bound into the hydrate formulae — the
What is left is pore solution, and its volume follows:
And the products occupy less space than the reactants they consumed — the Le Chatelier contraction. In a sealed specimen that deficit cannot be filled from outside, so it becomes empty porosity,
The pore space is therefore liquid plus void, and its degree of saturation is
Now suppose hydration stops when the saturation falls to some value
Two things follow, and the first is a warning about what this page can prove.
Powers' form is structural; only the coefficient is predicted
What is genuinely predicted is
2. The chemistry:
The cement is the one whose isotherm this page uses further down — taking the isotherm from one paper and the clinker from another would compare two materials. Baroghel-Bouny et al. (1999) give its mineral composition in their Table 2.
using ChemistryLab
using DynamicQuantities
using Logging
using OptimaSolver
using Printf
# Baroghel-Bouny et al. (1999): the cement of their Table 2, in mass percent, and
# the water-to-cement ratio of their mix CO, Table 1. The balance of Table 2 is
# free lime and alkalis, which this species list does not carry.
pct(c) = literature_row("BaroghelBouny1999", "cement_composition", c).content_percent / 100
compo = ["C3S" => pct("C3S"), "C2S" => pct("C2S"), "C3A" => pct("C3A"),
"C4AF" => pct("C4AF"), "Gp" => pct("Gypsum"), "Cal" => pct("CaCO3")]
cmass = sum(last.(compo))
wc = literature_row("BaroghelBouny1999", "mixes", "CO").W_C
phases = split("C3S C2S C3A C4AF Gp Anh Cal Portlandite Jennite H2O@ " *
"ettringite monosulphate12 C3AH6 C3FH6 C4FH13 monocarbonate")
substances = build_species(datapath("cemdata18-thermofun.json"); verbose = false)
cs = ChemicalSystem(
speciation(substances, phases; aggregate_state = [AS_AQUEOUS]), CEMDATA_PRIMARIES
)
iw = only(cs.idx_solvent)
ic = [findfirst(s -> symbol(s) == sym, cs.species) for (sym, _) in compo]
M_H2O = ustrip(us"kg/mol", cs.species[iw][:M])An equilibrium calculation has no notion of a reaction that has not happened, so the degree of hydration is imposed: a fraction
function paste(α)
mtot = cmass + wc * cmass
st = ChemicalState(cs)
for (sym, mfrac) in compo
set_quantity!(st, sym, α * mfrac / mtot * u"kg")
end
set_quantity!(st, "H2O@", wc * cmass / mtot * u"kg")
V = volume(st)
set_quantity!(st, "H+", 1e-7u"mol/L" * V.liquid)
set_quantity!(st, "OH-", 1e-7u"mol/L" * V.liquid)
return st
end
cement_mass(st) = sum(ustrip(us"kg", st.n[i] * cs.species[i][:M]) for i in ic)
function paste_budget(α)
fresh = paste(1.0)
mc = cement_mass(fresh)
w_tot = ustrip(us"mol", fresh.n[iw]) * M_H2O
eq, cert = equilibrate_certified(paste(α))
n = collect(eq.n) # put the unreacted cement back
for i in ic
n[i] += (1 - α) * fresh.n[i]
end
ϕ = porosity(ChemicalState(cs, n), fresh)
w_free = ustrip(us"mol", eq.n[iw]) * M_H2O
V_ref = ustrip(us"m^3", volume(fresh).total)
return (; b = (w_tot - w_free) / mc / α,
s = ϕ.void * V_ref / mc / α * 1e3, # m³/kg → cm³/g
w_free = w_free / mc, porosity = ϕ.total, certified = cert.optimal)
end
println(" α b (g/g) s (cm³/g) free water (g/g) porosity certified")
for α in (0.55, 0.60, 0.65, 0.70, 0.80)
r = paste_budget(α)
@printf("%5.2f %9.4f %10.4f %16.4f %9.4f %s\n",
α, r.b, r.s, r.w_free, r.porosity, r.certified)
end α b (g/g) s (cm³/g) free water (g/g) porosity certified
0.55 0.3095 0.0640 0.1698 0.3118 true
0.60 0.3095 0.0639 0.1543 0.2931 true
0.65 0.3095 0.0639 0.1389 0.2743 true
0.70 0.3095 0.0639 0.1234 0.2556 true
0.80 0.3095 0.0638 0.0924 0.2182 trueTwo things to read off that table.
Every solve is certified. For a convex problem the KKT conditions are sufficient, so certified = true is a proof that the composition is the Gibbs minimum and not the point an iteration stopped at — see Proving that an answer is the answer.
ref = paste_budget(0.65)
b_model, s_shrink = ref.b, ref.s
@printf("b = %.4f g/g (formula water) s = %.4f cm³/g (chemical shrinkage)\n",
b_model, s_shrink)b = 0.3095 g/g (formula water) s = 0.0639 cm³/g (chemical shrinkage)For scale: a chemical shrinkage of 0.064 cm³ per gram of cement is the value cement chemistry reports for an ordinary Portland cement, and it is here a consequence of the standard molar volumes in the database rather than an input.
3. The measurement:
The retention curve is not chemistry. It says how tightly a particular material holds the water still in it, it is measured, and it is the one thing on this page that cannot come out of a thermodynamic database.
Baroghel-Bouny et al. (1999) measured water-vapor desorption isotherms on two pastes and two concretes and fitted each with
which is the van Genuchten form written with VanGenuchten as
Baroghel-Bouny et al. (1999), Table 5, with the water-to-cement ratios of their Table 1. CO and CH are cement pastes, BO and BH concretes, and CH and BH contain 10 % silica fume:
fits = literature_table("BaroghelBouny1999", "retention_fit")
println("mix W/C a (MPa) b m = 1/b")
for (mix, a, b) in zip(fits.mix, fits.a, fits.b)
w_c = literature_row("BaroghelBouny1999", "mixes", mix).W_C
@printf("%-4s %5.2f %9.4f %6.4f %9.5f\n", mix, w_c, ustrip(a) / 1.0e6, b, 1 / b)
endmix W/C a (MPa) b m = 1/b
BO 0.48 18.6237 2.2748 0.43960
BH 0.26 46.9364 2.0601 0.48541
CO 0.34 37.5479 2.1684 0.46117
CH 0.19 96.2837 1.9540 0.51177The capillary pressure becomes a water activity through Kelvin, water_activity does:
fit_co = literature_row("BaroghelBouny1999", "retention_fit", "CO")
co = VanGenuchten(; a = fit_co.a, m = 1 / fit_co.b) # their mix CO
T_K = 298.15
# Liquid water at T_K: its molar volume from the HGK equation of state, its
# surface tension from the IAPWS equation.
V_m = ustrip(us"kg/mol", Species("H2O").M) / water_thermo_props(T_K, 1.0e5).D # m³/mol
γ_w = water_surface_tension(T_K) # N/m
a_w(S) = water_activity(co, S; V_m = V_m, T = T_K)
function saturation_at(rh) # the form is not invertible in closed form
lo, hi = 1e-4, 1 - 1e-12
for _ in 1:60
mid = (lo + hi) / 2
a_w(mid) > rh ? (hi = mid) : (lo = mid)
end
return (lo + hi) / 2
end
println(" RH S* p_c (MPa) Kelvin radius (nm)")
for rh in (0.75, 0.80, 0.85, 0.90, 0.95)
S = saturation_at(rh)
@printf(" %4.2f %6.4f %10.2f %18.2f\n", rh, S,
capillary_pressure(co, S) / 1e6,
kelvin_radius(rh; γ = γ_w, V_m = V_m, T = T_K) * 1e9)
end RH S* p_c (MPa) Kelvin radius (nm)
0.75 0.7107 39.47 3.65
0.80 0.7862 30.62 4.70
0.85 0.8619 22.30 6.46
0.90 0.9301 14.46 9.96
0.95 0.9800 7.04 20.45The last column sets the scale: 80 % relative humidity corresponds to a meniscus of radius 4.7 nm. That is the gel-pore scale, which is why the water Powers (1948) assigns to gel pores and the water a sealed paste cannot use are the same water.
using Plots
Ss = range(0.40, 0.999; length = 200)
p1 = plot(Ss, a_w.(Ss); xlabel = "degree of saturation S", ylabel = "water activity a_w",
label = "mix CO, measured fit", linewidth = 2, color = :steelblue,
title = "Desorption isotherm (Baroghel-Bouny et al. 1999)", legend = :bottomright)
hline!(p1, [0.80]; label = "RH = 0.80", linestyle = :dash, color = :firebrick)
vline!(p1, [saturation_at(0.80)]; label = "S* = $(round(saturation_at(0.80); digits = 4))",
linestyle = :dash, color = :seagreen)
plot(p1; size = (700, 420), left_margin = 8Plots.mm, bottom_margin = 8Plots.mm)
4. Closing the budget
Everything the boxed formula needs is now in hand, each from its own source.
powers_k(b, s, S★) = b + s * S★ / (1 - S★)
k_powers = literature_value("Powers1948", "w_c_sealed")
S80 = saturation_at(0.80)
@printf("S* = %.4f so S*/(1-S*) = %.4f\n\n", S80, S80 / (1 - S80))
@printf("with the model's formula water, b = %.4f:\n", b_model)
@printf(" k = %.4f α_max(w/c=%.2f) = %.4f (Powers: %.2f and %.4f)\n",
powers_k(b_model, s_shrink, S80), wc, wc / powers_k(b_model, s_shrink, S80),
k_powers, min(1.0, wc / k_powers))S* = 0.7862 so S*/(1-S*) = 3.6774
with the model's formula water, b = 0.3095:
k = 0.5444 α_max(w/c=0.34) = 0.6245 (Powers: 0.42 and 0.8095)30 % above Powers (1948). Before calling that a disagreement, look at what
The 0.42 of Powers (1948) is itself a sum: about 0.23 g/g of non-evaporable water plus about 0.19 g/g of gel water. And his 0.23 is an operational quantity — the water that survives D-drying, over
b_powers = literature_value("Powers1948", "non_evaporable_water") # his own number
@printf("model formula water %.4f g/g\n", b_model)
@printf("Powers' w_n %.4f g/g (defined by D-drying)\n", b_powers)
@printf("difference %.4f g/g — interlayer water, counted differently\n",
b_model - b_powers)model formula water 0.3095 g/g
Powers' w_n 0.2300 g/g (defined by D-drying)
difference 0.0795 g/g — interlayer water, counted differentlySubstituting the non-evaporable water of Powers (1948), and leaving the chemical shrinkage and the isotherm exactly as measured:
@printf("with Powers' w_n = %.2f:\n", b_powers)
@printf(" k = %.4f α_max(w/c=%.2f) = %.4f (Powers: %.2f and %.4f)\n",
powers_k(b_powers, s_shrink, S80), wc, wc / powers_k(b_powers, s_shrink, S80),
k_powers, min(1.0, wc / k_powers))with Powers' w_n = 0.23:
k = 0.4650 α_max(w/c=0.34) = 0.7312 (Powers: 0.42 and 0.8095)11 % above. The remaining gap is a fifth of a point of saturation, and §6 shows how little that is.
5. The result: inverting Powers
The comparison above runs one way — assume the arrest is at 80 % and see what
function invert_k(b, s, k_target)
lo, hi = 1e-6, 1 - 1e-12
for _ in 1:80
mid = (lo + hi) / 2
powers_k(b, s, mid) > k_target ? (hi = mid) : (lo = mid)
end
return (lo + hi) / 2
end
@printf("k = %.2f implies:\n", k_powers)
for (label, b) in (("the model's formula water", b_model), ("Powers' own w_n", b_powers))
S = invert_k(b, s_shrink, k_powers)
@printf(" b = %.4f (%-26s) S* = %.4f RH = %.4f\n", b, label, S, a_w(S))
endk = 0.42 implies:
b = 0.3095 (the model's formula water ) S* = 0.6337 RH = 0.6956
b = 0.2300 (Powers' own w_n ) S* = 0.7483 RH = 0.7752With the non-evaporable water of Powers (1948), his 0.42 corresponds to a hydration arrest at 77.5 % relative humidity.
That number was not put in anywhere. It came out of a chemical shrinkage computed from a thermodynamic database, a desorption isotherm measured for a drying study, and the water split of Powers (1948) — three sources, none of which knows what the others are for. And 75–80 % is the window in which sealed pastes are independently reported to stop hydrating; it is the same window humidity_factor implements as its cut-off.
A corroboration from the same paper's own measurements
(Baroghel-Bouny et al., 1999) also measured the internal relative humidity of their sealed specimens at 28 days, in their Table 3: 97 % for the paste CO (W/C 0.34) and 88.5 % for the paste CH (W/C 0.19); 97 % for the concrete BO (W/C 0.48) and 77.5 % for the concrete BH (W/C 0.26).
Monotone in W/C within each material class, as self-desiccation requires. And the mixes Powers (1948) says should have arrested by 28 days — the low-W/C ones — are the ones sitting in the 77.5–88.5 % band that the inversion above points at, while the two mixes with water to spare sit at 97 %. Those specimens are not arrested pastes at equilibrium, so this is corroboration of a range and not a fourth decimal; it is worth stating because the numbers come from the same table as the isotherm.
6. How much of this is the isotherm?
@printf("dk/dS* = %.3f at S* = %.4f\n\n", s_shrink / (1 - S80)^2, S80)
@printf(" RH S* k (model b) k (Powers w_n) α_max(%.2f)\n", wc)
for rh in (0.70, 0.75, 0.80, 0.85, 0.90)
S = saturation_at(rh)
@printf(" %4.2f %6.4f %13.4f %16.4f %12.4f\n", rh, S,
powers_k(b_model, s_shrink, S), powers_k(b_powers, s_shrink, S),
min(1.0, wc / powers_k(b_powers, s_shrink, S)))
enddk/dS* = 1.398 at S* = 0.7862
RH S* k (model b) k (Powers w_n) α_max(0.34)
0.70 0.6397 0.4229 0.3434 0.9900
0.75 0.7107 0.4664 0.3870 0.8786
0.80 0.7862 0.5444 0.4650 0.7312
0.85 0.8619 0.7083 0.6289 0.5406
0.90 0.9301 1.1600 1.0806 0.3146rhs = range(0.65, 0.92; length = 120)
ks_m = [powers_k(b_model, s_shrink, saturation_at(r)) for r in rhs]
ks_p = [powers_k(b_powers, s_shrink, saturation_at(r)) for r in rhs]
p2 = plot(rhs, ks_m; xlabel = "assumed arrest humidity", ylabel = "k = (w/c) / α_max",
label = "b = model formula water", linewidth = 2, color = :steelblue,
title = "What the coefficient depends on", legend = :topleft, ylims = (0.25, 1.0))
plot!(p2, rhs, ks_p; label = "b = Powers' w_n = $(b_powers)", linewidth = 2, color = :firebrick)
k_saturated = literature_value("Powers1948", "w_c_saturated")
hline!(p2, [k_powers]; label = "Powers $(k_powers)", linestyle = :dash, color = :black)
hline!(p2, [k_saturated]; label = "Powers $(k_saturated), saturated curing", linestyle = :dot, color = :gray)
plot(p2; size = (700, 420), left_margin = 8Plots.mm, bottom_margin = 8Plots.mm)
A tenth of a point of saturation moves
7. Two negative controls
Two claims carry the rest of the page: that the arrest criterion belongs in the rate law and not in the Gibbs energy, and that the proportional form of Powers (1948) is no evidence for his coefficient. Both are testable on the material at hand, so both are measured here rather than asserted.
The thermodynamic route, measured
CapillaryWater imposes a water activity on the equilibrium itself. If self-desiccation arrested hydration thermodynamically, imposing the arrest humidity on a paste with all its cement available would leave some of that cement unreacted. It does not:
fresh = paste(1.0)
ip = findfirst(s -> symbol(s) == "Portlandite", cs.species)
ij = findfirst(s -> symbol(s) == "Jennite", cs.species)
println(" a_w imposed n(Portlandite) n(Jennite) free water (mol) certified")
for aw in (1.00, 0.90, 0.80, 0.50)
# The last column is the certificate; the warning of a refusal would repeat it.
eq, cert = with_logger(NullLogger()) do
equilibrate_certified(
paste(1.0); constraint = CapillaryWater(_ -> aw; reference = fresh)
)
end
@printf(" %11.2f %14.6f %10.6f %16.5f %s\n", aw,
ustrip(us"mol", eq.n[ip]), ustrip(us"mol", eq.n[ij]),
ustrip(us"mol", eq.n[iw]), cert.optimal)
end a_w imposed n(Portlandite) n(Jennite) free water (mol) certified
1.00 2.659505 2.967346 1.26544 true
0.90 2.659505 2.967345 1.26544 true
0.80 2.659505 2.967345 1.26545 true
0.50 2.659210 2.967342 2.42921 trueFrom saturation down to the arrest humidity the assemblage is identical to six digits — the imposed shift moves the answer by less than the printing precision — and each of those answers is certified. The last row is there to show the other half of that: further down, no route certifies, the last column says so, and the answer equilibrate_certified returns is a KKT candidate rather than a proof. The arithmetic in the opening admonition says why nothing is expected to happen in that range anyway.
So a Gibbs minimization under CapillaryWater arrests where it ran out of water stoichiometrically, not where Powers (1948) says. What the constraint is good for is what §3–§5 use it for: making the water activity of a state mean the water in the pores rather than a mole fraction, and handing that humidity to a rate law through PoreHumidity.
The proportional form, measured
§1 argued that
fit_bh = literature_row("BaroghelBouny1999", "retention_fit", "BH")
bh = VanGenuchten(; a = fit_bh.a, m = 1 / fit_bh.b) # their mix BH
function saturation_of(law, rh)
lo, hi = 1e-4, 1 - 1e-12
for _ in 1:60
mid = (lo + hi) / 2
water_activity(law, mid; V_m = V_m, T = T_K) > rh ? (hi = mid) : (lo = mid)
end
return (lo + hi) / 2
end
for (name, law) in (("CO", co), ("BH", bh))
S = saturation_of(law, 0.80)
k = powers_k(b_model, s_shrink, S)
@printf("\ncurve %s: S* = %.4f k = %.4f\n", name, S, k)
for w in (0.25, 0.30, 0.35, 0.40)
@printf(" w/c %.2f → α_max = %.6f α_max/(w/c) = %.6f\n", w, w / k, 1 / k)
end
end
curve CO: S* = 0.7862 k = 0.5444
w/c 0.25 → α_max = 0.459188 α_max/(w/c) = 1.836754
w/c 0.30 → α_max = 0.551026 α_max/(w/c) = 1.836754
w/c 0.35 → α_max = 0.642864 α_max/(w/c) = 1.836754
w/c 0.40 → α_max = 0.734702 α_max/(w/c) = 1.836754
curve BH: S* = 0.8389 k = 0.6423
w/c 0.25 → α_max = 0.389229 α_max/(w/c) = 1.556918
w/c 0.30 → α_max = 0.467075 α_max/(w/c) = 1.556918
w/c 0.35 → α_max = 0.544921 α_max/(w/c) = 1.556918
w/c 0.40 → α_max = 0.622767 α_max/(w/c) = 1.556918The last column is the same number at every
Neither curve gives 0.42 at RH 0.80 with this system's formula water. Recovering the coefficient of Powers (1948) takes his own water split as well, and §5 does that explicitly and says so.
8. What went in, and what came out
| input | value | source | is that source about Powers? |
|---|---|---|---|
| clinker composition | Table 2 | Baroghel-Bouny et al. (1999) | no |
| thermodynamic data | Cemdata18 | Lothenbach et al. (2019) | no |
| 0.3095 g/g, 0.0639 cm³/g | computed here, certified | no | |
| retention curve | mix CO | Baroghel-Bouny et al. (1999) Table 5 | no |
| 0.0720 N/m, 1.807e-5 m³/mol, 298.15 K | water at 25 °C: IAPWS (2014), and the HGK equation of state (Haar et al., 1984) | no | |
| Powers' own split of his 0.42 | Powers (1948) | yes |
output: the arrest humidity, 77.5 %.
Only the last row knows about Powers (1948), and it contributes his decomposition, not his coefficient. Nothing on this page was adjusted to improve the agreement, and the one number that could have been — the retention curve — was fitted by Baroghel-Bouny et al. (1999) to a drying experiment two decades before this calculation existed.
9. Assumptions, and where each one bites
The degree of hydration is imposed, not predicted. §2's construction reacts a fraction of the cement and leaves the rest inert. What this page computes is the
at which the arrest criterion is met, not a trajectory in time. For that, hand PoreHumiditytoparrott_killoh_avramias itshumidityand integrate.One retention curve for an evolving pore structure. The measured isotherm belongs to a mature paste; the page applies it at every
. The pore structure of a young paste is coarser, so its true is higher and the arrest earlier. The published fit excludes the humidity range this page works in. Baroghel-Bouny et al. (1999) state that "the experimental data corresponding to the lowest capillary pressures (corresponding to the highest RH) have not been accounted [for] due to their weak reliability". Above about 90 % RH the curve used here is extrapolation into a region its authors declined to fit — which is exactly where an unarrested paste sits, and a reason the corroboration in §5 is stated as a range.
The saturations are normalized differently. The
of Baroghel-Bouny et al. (1999) is relative to their measured total porosity, 30.3 % for mix CO; this model gives 27.4 % at , and by §6 worth a few hundredths of . The C-S-H model matters as much as the isotherm.
one for one. Measured on this system: Jennitegives= 0.3095 g/g and the CSHQsolid solution gives 0.3684, so it moves away from the value of Powers (1948), not toward it. Neither is wrong; they draw the formula-water line in different places.No alkalis. This species list carries none, so the pore solution's osmotic depression of
is absent. At 0.1–0.5 mol/kg it is worth one to two points of relative humidity, in the same direction as the capillary term. Ideal molar volumes, and a closed species list — the standing assumptions of every equilibrium page here.
10. Reproducing this
julia --project=scripts scripts/self_desiccation_powers.jlEvery block above is executed when this page is built, so the numbers printed are whatever the code produced. The script computes the same quantities outside Documenter and prints the intermediate ones this page summarizes.
Where to go next
The degree of hydration is an input on this page; Self-desiccation in time makes it the result of a rate law that reads the internal humidity. What a Gibbs minimization does without any arrest criterion is shown in the w/c example. The API entries are CapillaryWater, PoreHumidity, WaterRetention and powers_alpha_max.