ChemistryLab.jl

Appearance

CO₂ Dissolution and the Carbonate System ​

Before this page

The tutorial Chemical Equilibrium.

The carbonate system — CO₂(aq) / HCO₃⁻ / CO₃²⁻ — is one of the most important chemical equilibria in natural water chemistry, geochemistry, and environmental science. It controls the pH of rainwater, rivers, oceans, and aquifers, and governs the dissolution and precipitation of carbonate minerals such as calcite and dolomite.

This example shows how to:

  1. Load species from a thermodynamic database and build a ChemicalSystem.

  2. Extract dissociation constants (pKₐ₁, pKₐ₂) directly from standard Gibbs energies.

  3. Simulate CO₂ dissolution in pure water over several orders of magnitude.

  4. Construct the analytical speciation diagram (fraction of each carbonate species vs pH).

System setup ​

All species are loaded from the SLOP98 inorganic database. The system contains only H, O, C and charge-bearing species; no other elements are needed.

SymbolSpeciesPhase
H2O@H₂O — wateraqueous solvent
H+H⁺ — protonaqueous solute
OH-OH⁻ — hydroxideaqueous solute
CO2@CO₂(aq) — dissolved CO₂aqueous solute
HCO3-HCO₃⁻ — bicarbonateaqueous solute
CO3-2CO₃²⁻ — carbonateaqueous solute

The three primary species H2O@, H+, CO3-2 form a basis for the whole system: every other aqueous species is a linear combination of these three.

julia
using Optimization, OptimizationIpopt
using ChemistryLab
using DynamicQuantities

substances = build_species(datapath("slop98-inorganic-thermofun.json"))

# Select only the six carbonate-system species
dict = Dict(symbol(s) => s for s in substances)
species = [dict[sym] for sym in split("H2O@ H+ OH- CO2@ HCO3- CO3-2")]

cs = ChemicalSystem(species, ["H2O@", "H+", "CO3-2"])
The chemical system in full
julia
cs
6-element ChemicalSystem{Species{Int64}, AbstractReaction, StoichMatrix{Int64, Symbol, Vector{Symbol}, Matrix{Int64}, Species{Int64}}, StoichMatrix{Int64, Species{Int64}, Vector{Species{Int64}}, Matrix{Int64}, Species{Int64}}, Nothing, Nothing}:
 H2O@ {Water HGK} [H2O@ ◆ H₂O@]
 H+ {H+} [H+ ◆ H⁺]
 OH- {OH- hydroXyl ion} [OH- ◆ OH⁻]
 CO2@ {CO2,aq (+ H2O = H2CO3,aq )} [CO2@ ◆ CO₂@]
 HCO3- {HCO3- bicarbonate ion} [HCO3- ◆ HCO₃⁻]
 CO3-2 {CO3-2 carbonate ion} [CO3-2 ◆ CO₃²⁻]

Larger carbonate systems

If calcium is relevant (e.g. hard water, calcite dissolution), add Ca+2 and its complexes to the species list and extend the primaries to ["H2O@", "H+", "CO3-2", "Ca+2"]. The SLOP98 inorganic database also contains Ca(HCO3)+ and Ca(CO3)@.

Dissociation constants from thermodynamic data ​

The two acid–base equilibria of the carbonate system are:

Both constants are computed from standard Gibbs energies of formation at 25 °C, exactly as in the acetic acid titration example:

julia
R  = ustrip(Constants.R)
T  = 298.15   # K
RT = R * T

sp = Dict(symbol(s) => s for s in cs.species)

# pKa1 : CO₂(aq) + H2O ⇌  HCO₃⁻ + H⁺
Ka1  = exp(-(sp["HCO3-"].ΔₐG⁰(T = T) + sp["H+"].ΔₐG⁰(T = T) - sp["H2O@"].ΔₐG⁰(T = T) - sp["CO2@"].ΔₐG⁰(T = T)) / RT)
pKa1 = -log10(Ka1)

# pKa2 : HCO₃⁻ ⇌ H⁺ + CO₃²⁻
Ka2  = exp(-(sp["CO3-2"].ΔₐG⁰(T = T) + sp["H+"].ΔₐG⁰(T = T) - sp["HCO3-"].ΔₐG⁰(T = T)) / RT)
pKa2 = -log10(Ka2)

# pKw : H₂O ⇌ H⁺ + OH⁻
Kw   = exp(-(sp["H+"].ΔₐG⁰(T = T) + sp["OH-"].ΔₐG⁰(T = T) - sp["H2O@"].ΔₐG⁰(T = T)) / RT)
pKw  = -log10(Kw)

println("pKa1 (CO₂(aq) / HCO₃⁻) = ", round(pKa1, digits = 2), "   (lit. 6.35)")
println("pKa2 (HCO₃⁻  / CO₃²⁻)  = ", round(pKa2, digits = 2), "   (lit. 10.33)")
println("pKw  (H₂O    / H⁺+OH⁻) = ", round(pKw,  digits = 2), "   (lit. 14.00)")
pKa1 (CO₂(aq) / HCO₃⁻) = 6.34   (lit. 6.35)
pKa2 (HCO₃⁻  / CO₃²⁻)  = 10.33   (lit. 10.33)
pKw  (H₂O    / H⁺+OH⁻) = 14.0   (lit. 14.00)

Building the equilibrium solver ​

A single EquilibriumSolver is compiled once and reused for every point of the scan:

julia
using Optimization, OptimizationIpopt

solver = EquilibriumSolver(
    cs,
    DiluteSolutionModel(),
    IpoptOptimizer(
        acceptable_tol        = 1e-10,
        dual_inf_tol          = 1e-10,
        acceptable_iter       = 1000,
        constr_viol_tol       = 1e-10,
        warm_start_init_point = "no",
        additional_options    = Dict{String, Any}("print_level" => 0, "sb" => "yes"),
    );
    variable_space = Val(:linear),
    reltol  = 1e-8,
)
EquilibriumSolver{ChemistryLab.var"#μ#build_potentials##0"{ChemistryLab.var"#lna#activity_model##0"{Type{Float64}, ChemistryLab._MixingTerms{Vector{Vector{Int64}}, Nothing, Vector{Vector{Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing}, Nothing, Bool, Float64, Vector{Int64}, Vector{Int64}, Int64, Bool}}, IpoptOptimizer, Val{:linear}, DiluteSolutionModel}(ChemistryLab.var"#μ#build_potentials##0"{ChemistryLab.var"#lna#activity_model##0"{Type{Float64}, ChemistryLab._MixingTerms{Vector{Vector{Int64}}, Nothing, Vector{Vector{Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing}, Nothing, Bool, Float64, Vector{Int64}, Vector{Int64}, Int64, Bool}}(ChemistryLab.var"#lna#activity_model##0"{Type{Float64}, ChemistryLab._MixingTerms{Vector{Vector{Int64}}, Nothing, Vector{Vector{Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing}, Nothing, Bool, Float64, Vector{Int64}, Vector{Int64}, Int64, Bool}(Float64, ChemistryLab._MixingTerms{Vector{Vector{Int64}}, Nothing, Vector{Vector{Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing}(false, Vector{Int64}[], nothing, false, Vector{Int64}[], nothing, nothing, nothing, false, 1, Int64[], Float64[], 1.0, nothing, nothing), nothing, false, 4.016550535127659, Int64[], [2, 3, 4, 5, 6], 1, true)), IpoptOptimizer(1.0e-10, 1000, 1.0e-10, 1.0e-10, 0.0001, "no", "no", "mumps", "none", "", "", "yes", "gradient-based", 100.0, "no", "no", "monotone", "quality-function", 0.1, "obj-constr-filter", "no", "exact", 6, "bfgs", "no", "filter", "no", Dict{String, Any}("sb" => "yes", "print_level" => 0)), Val{:linear}(), Base.Pairs(:reltol => 1.0e-8), DiluteSolutionModel())

CO₂ dissolution scan ​

The experiment: 1 kg of pure water receives an increasing amount of dissolved CO₂ (CO₂(aq)), from trace atmospheric levels (~10⁻⁵ mol) up to a highly carbonated solution (~0.1 mol). Three representative concentrations illustrate the pH trend:

julia
s = ChemicalState(cs)
for (n_CO2, label) in [
    (1e-5, "atmospheric, ~400 ppm"),
    (1e-3, "lightly carbonated"),
    (1e-1, "strongly carbonated"),
]
    set_quantity!(s, "H2O@", 1.0u"kg")
    set_quantity!(s, "CO2@", n_CO2 * u"mol")
    V_liq = volume(s).liquid
    set_quantity!(s, "H+",  1e-7u"mol/L" * V_liq)
    set_quantity!(s, "OH-", 1e-7u"mol/L" * V_liq)
    s_eq = solve(solver, s)
    println("pH at [CO₂]₀ = ", n_CO2, " mol (", label, ") : ",
            round(pH(s_eq), digits = 2))
end
pH at [CO₂]₀ = 1.0e-5 mol (atmospheric, ~400 ppm) : 5.72
pH at [CO₂]₀ = 0.001 mol (lightly carbonated) : 4.68
pH at [CO₂]₀ = 0.1 mol (strongly carbonated) : 3.68

For the full scan used to generate the plots below, loop over 60 log-spaced points:

julia
sp_idx = Dict(symbol(s) => i for (i, s) in enumerate(cs.species))

n_CO2_range = exp.(range(log(1e-5), log(0.1); length = 60))   # mol, in 1 kg H₂O

pH_vals     = Float64[]
f_CO2_vals  = Float64[]
f_HCO3_vals = Float64[]
f_CO3_vals  = Float64[]

s = ChemicalState(cs)
for n_CO2 in n_CO2_range
    set_quantity!(s, "H2O@",  1.0u"kg")
    set_quantity!(s, "CO2@",  n_CO2 * u"mol")

    V_liq = volume(s).liquid
    set_quantity!(s, "H+",  1e-7u"mol/L" * V_liq)   # pH-neutral seed
    set_quantity!(s, "OH-", 1e-7u"mol/L" * V_liq)

    s_eq = solve(solver, s)
    push!(pH_vals, pH(s_eq))

    n_a = ustrip(s_eq.n[sp_idx["CO2@"]])
    n_b = ustrip(s_eq.n[sp_idx["HCO3-"]])
    n_c = ustrip(s_eq.n[sp_idx["CO3-2"]])
    n_C = n_a + n_b + n_c

    push!(f_CO2_vals,  n_a / n_C)
    push!(f_HCO3_vals, n_b / n_C)
    push!(f_CO3_vals,  n_c / n_C)
end

Result: pH as a function of dissolved CO₂ ​

julia
using Plots

p1 = plot(
    log10.(n_CO2_range), pH_vals;
    xlabel    = "log₁₀([CO₂]₀) (mol/L)",
    ylabel    = "Equilibrium pH",
    label     = "Numerical (ChemistryLab)",
    linewidth = 2,
    marker    = :circle,
    markersize = 3,
    color     = :steelblue,
    title     = "CO₂ dissolution in pure water (25 °C)",
    ylims     = (3, 7),
    legend    = :right,
)
hline!(p1, [pKa1]; linestyle = :dot, color = :orange, label = "pKa1 = $(round(pKa1, digits=2))")
vline!(p1, [log10(1.4e-5)]; linestyle = :dash, color = :grey,
       label = "Atmospheric CO₂ (~400 ppm)")
p1

Speciation diagram ​

The distribution of carbonate species as a function of pH follows from the two pKₐ values. Using the computed constants, the analytical mole fractions are:

julia
Ka1_val = 10^(-pKa1)
Ka2_val = 10^(-pKa2)

function carbonate_fractions(ph)
    H     = 10^(-ph)
    denom = H^2 + Ka1_val * H + Ka1_val * Ka2_val
    return H^2 / denom, Ka1_val * H / denom, Ka1_val * Ka2_val / denom
end

pH_axis = range(2, 14; length = 1000)
fracs   = carbonate_fractions.(pH_axis)
f0 = [f[1] for f in fracs]
f1 = [f[2] for f in fracs]
f2 = [f[3] for f in fracs]

p2 = plot(
    collect(pH_axis), [f0 f1 f2];
    xlabel    = "pH",
    ylabel    = "Mole fraction",
    label     = ["CO₂(aq)" "HCO₃⁻" "CO₃²⁻"],
    linewidth = 2,
    color     = [:steelblue :orange :green],
    title     = "Carbonate speciation (25 °C, pKₐ from database)",
    legend    = :right,
)
vline!(p2, [pKa1]; linestyle = :dash, color = :grey,   label = "pKa1 = $(round(pKa1, digits=2))")
vline!(p2, [pKa2]; linestyle = :dash, color = :black,  label = "pKa2 = $(round(pKa2, digits=2))")

# Overlay the numerical points from the CO₂ scan
scatter!(p2, pH_vals, f_CO2_vals;  color = :steelblue, markersize = 3, label = "")
scatter!(p2, pH_vals, f_HCO3_vals; color = :orange,    markersize = 3, label = "")
p2

Analysis ​

ZonepH rangeDominant species
Acidic / high CO₂pH < 6.35 (pKₐ₁)CO₂(aq)
Buffer zone 1pH ≈ pKₐ₁ = 6.35CO₂(aq) / HCO₃⁻ (equal amounts)
Mid-range6.35 < pH < 10.33HCO₃⁻
Buffer zone 2pH ≈ pKₐ₂ = 10.33HCO₃⁻ / CO₃²⁻ (equal amounts)
AlkalinepH > 10.33CO₃²⁻

  • Atmospheric CO₂ (~10⁻⁵ mol/L) produces a pH of ≈ 5.6 in pure water — slightly acidic rain even without any pollutants.

  • Carbonated water (industrial, ~0.05 mol/L CO₂) reaches pH ≈ 4.0–4.3.

  • Because dissolved CO₂ always pushes the system into the CO₂(aq)-dominated zone (pH < pKₐ₁), the HCO₃⁻ and CO₃²⁻ domains are only accessible by adding a base (e.g. NaOH, CaCO₃ dissolution) — which is the essence of ocean buffering and concrete alkalinity.

Link to cement chemistry

In the pore solution of hydrated cement (pH ≈ 12.5–13), the dominant carbonate species is CO₃²⁻. When atmospheric CO₂ penetrates the concrete cover, the carbonate system shifts leftward until portlandite (Ca(OH)₂) is consumed and the pH drops below 9. This process — carbonation — is the main cause of corrosion initiation in reinforced concrete, and is explored in the cement carbonation example.

Where to go next ​

The same equilibria in a cement pore solution, where they drive the carbonation of the hydrates, are computed in Carbonation of a CEM I paste. The first calculation on a complete cement is A CEM I from its clinker phases.