Equilibrium

abstract type AbstractActivityModel end

Base type for all activity models. Each concrete subtype must implement activity_model(cs::ChemicalSystem, model::AbstractActivityModel) which returns a closure a(n, p) -> Vector{Float64} of log-activities.

struct DiluteSolutionModel <: AbstractActivityModel

Ideal dilute solution model:

• Solvent: Raoult's law — ln a = ln(x_solvent)
• Solutes: Henry's law — ln a = ln(c_i / c°) where c° = 1 mol/L
• Crystals: pure solid — ln a = 0
• Gas: ideal mixture — ln a = ln(x_i)
• SS end-members: ideal mixing — ln a = ln(xᵢ) within the solid-solution phase

struct HKFActivityModel{T<:Real} <: AbstractActivityModel

Extended Debye-Hückel activity model with B-dot term (Helgeson 1969 / Helgeson, Kirkham & Flowers 1981). This is the model used by PHREEQC and EQ3/6.

Activity coefficient formulas

For ionic species (charge z ≠ 0):

log₁₀ γᵢ = −A zᵢ² √I / (1 + B åᵢ √I)  +  Ḃ I

For neutral aqueous species (charge z = 0):

log₁₀ γᵢ = Kₙ I

Log-activity (molality convention, standard state m° = 1 mol/kg):

ln aᵢ = ln(10) × log₁₀(γᵢ) + ln(mᵢ)

Water activity is computed from the osmotic coefficient φ (Gibbs-Duhem):

ln a_w = −Mw × Σmⱼ × φ

Ionic strength: I = ½ Σ mⱼ zⱼ²

Fields

• A: Debye-Hückel A parameter [(kg/mol)^(1/2)]. Default 0.5114 at 25 °C/1 bar.
• B: Debye-Hückel B parameter [Å⁻¹(kg/mol)^(1/2)]. Default 0.3288.
• Ḃ: B-dot extended term [kg/mol]. Default 0.041.
• Kₙ: salting-out coefficient for neutral species [kg/mol]. Default 0.1.
• å_default: global fallback effective ionic radius [Å]. Default 3.72.
• temperature_dependent: if true, recompute A and B from p.T, p.P at each call to the activity closure (requires T and P in p). Default false.

Ionic radius lookup

The effective radius åᵢ is resolved in order:

1. sp[:å] — explicit value in the species properties dict.
2. REJ_HKF — Helgeson et al. (1981) Table 3, keyed by PHREEQC formula.
3. REJ_CHARGE_DEFAULT — fallback by formal charge.
4. model.å_default — global fallback.

Valid range

I ≲ 1 mol/kg. Beyond ~2 mol/kg, the Pitzer model is recommended.

References

• Helgeson, H.C. (1969). Am. J. Sci. 267, 729–804.
• Helgeson, H.C., Kirkham, D.H. & Flowers, G.C. (1981). Am. J. Sci. 281, 1249–1516.
• Parkhurst, D.L. & Appelo, C.A.J. (2013). USGS Techniques Methods, Book 6, ch. A43.

Examples

# Default model at 25 °C / 1 bar (fixed A, B)
model = HKFActivityModel()

# Temperature-dependent A and B (recomputed at each solve)
model_tdep = HKFActivityModel(temperature_dependent=true)

state_eq = equilibrate(state; model=HKFActivityModel())

struct DaviesActivityModel{T<:Real} <: AbstractActivityModel

Davies (1962) activity model — a simplified Debye-Hückel without ionic radii.

Activity coefficient formula

For ionic species (charge z ≠ 0):

log₁₀ γᵢ = −A zᵢ² (√I / (1 + √I)  −  b I)

For neutral aqueous species (charge z = 0):

log₁₀ γᵢ = bₙ I

Log-activity (molality convention): ln aᵢ = ln(10) × log₁₀(γᵢ) + ln(mᵢ)

Water activity: Raoult approximation ln a_w = ln(x_w) (V1; accurate for I < 0.1).

Fields

• A: Debye-Hückel A parameter [(kg/mol)^(1/2)]. Default 0.5114 at 25 °C/1 bar.
• b: Davies empirical constant (default 0.3).
• bₙ: salting-out coefficient for neutral species (default 0.1).
• temperature_dependent: if true, recompute A from p.T, p.P at each call.

References

• Davies, C.W. (1962). Ion Association. Butterworths, London.

activity_model(cs::ChemicalSystem, ::DiluteSolutionModel) -> Function

Return a closure lna(n, p) -> Vector{Float64} computing the vector of log-activities for the dilute ideal solution model.

The returned function has signature lna(n, p) where:

• n: dimensionless mole vector (same indexing as cs.species)
• p: NamedTuple containing at least ϵ (floor value to avoid log(0))

Solid-solution end-members (class SC_SSENDMEMBER) receive ln aᵢ = ln xᵢ where xᵢ = nᵢ / Σnⱼ within the same solid-solution phase.

All quantities are dimensionless — units are stripped at construction time.

activity_model(cs::ChemicalSystem, model::HKFActivityModel) -> Function

Return a closure lna(n, p) -> Vector computing log-activities for the extended Debye-Hückel (B-dot) model of Helgeson (1969).

The closure captures all species indices and ionic radii at construction time. Inside lna:

• Solutes: molality convention, B-dot formula for ions, salting-out for neutrals.
• Solvent: osmotic coefficient from Gibbs-Duhem (σ-function).
• Crystals: ln a = 0 (pure solid).
• Gas: ideal mixture ln a = ln(xᵢ).

If model.temperature_dependent=true, p must contain T (K) and P (Pa) — both are provided automatically by _build_params.

AD-compatible: all closure computations accept ForwardDiff.Dual inputs.

activity_model(cs::ChemicalSystem, model::DaviesActivityModel) -> Function

Return a closure lna(n, p) -> Vector computing log-activities for the Davies (1962) model. No species-specific ionic radii are required.

AD-compatible: all closure computations accept ForwardDiff.Dual inputs.

build_potentials(cs::ChemicalSystem, model::AbstractActivityModel) -> Function

Return a closure μ(n, p) -> Vector{Float64} computing dimensionless chemical potentials μ_i / RT for all species.

$\mu_i / RT = \Delta_a G_i^0 / RT + \ln a_i$

The returned function is compatible with SciML solvers:

• n: dimensionless mole vector
• p: NamedTuple containing:
  • ΔₐG⁰overT: vector of standard Gibbs energies of formation divided by RT
  • ϵ: regularization floor (e.g. 1e-30)

All quantities are dimensionless — caller is responsible for stripping units from ΔₐG⁰overT before passing them in p.

Examples

julia> cs = ChemicalSystem([
           Species("H2O"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("Na+"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE),
       ]);

julia> μ = build_potentials(cs, DiluteSolutionModel());

julia> n = [55.5, 0.1];

julia> p = (ΔₐG⁰overT = [-95.6, -105.6], ϵ = 1e-30);

julia> length(μ(n, p)) == 2
true

REJ_HKF::Dict{String,Float64}

Effective electrostatic radii åᵢ [Å] for aqueous ions from Helgeson, Kirkham & Flowers (1981), Am. J. Sci. 281, Table 3.

Keys are PHREEQC-format formula strings (e.g. "Na+", "Ca+2", "SO4-2"). Used by HKFActivityModel with priority 2 in the radius lookup chain: sp[:å] > REJ_HKF > REJ_CHARGE_DEFAULT > model.å_default.

See also: REJ_CHARGE_DEFAULT, HKFActivityModel.

REJ_CHARGE_DEFAULT::Dict{Int,Float64}

Fallback effective electrostatic radii åᵢ [Å] indexed by formal charge, from ToughReact V2 (Xu et al. 2011, Table A2; after Helgeson et al. 1981).

Used by HKFActivityModel with priority 3 in the radius lookup chain: sp[:å] > REJ_HKF > REJ_CHARGE_DEFAULT > model.å_default.

See also: REJ_HKF, HKFActivityModel.

hkf_debye_huckel_params(T_K, P_Pa) -> NamedTuple{(:A, :B)}

Compute the Debye-Hückel A and B parameters from the water density ρ [g/cm³] and dielectric constant εᵣ at temperature T_K (K) and pressure P_Pa (Pa).

Formulas (Helgeson et al. 1981):

A(T,P) = 1.824829238×10⁶ × √ρ / (εᵣ T)^(3/2)    [(kg/mol)^(1/2)]
B(T,P) = 50.29158649      × √ρ / √(εᵣ T)          [Å⁻¹ (kg/mol)^(1/2)]

where ρ is in g/cm³.

Uses water_thermo_props (HGK equation of state) and water_electro_props_jn (Johnson-Norton dielectric constant).

AD-compatible (ForwardDiff-safe). Returns a NamedTuple (A=..., B=...).

Examples

julia> p = hkf_debye_huckel_params(298.15, 1e5);

julia> isapprox(p.A, 0.5114; rtol=1e-3)
true

julia> isapprox(p.B, 0.3288; rtol=1e-3)
true

abstract type AbstractSolidSolutionModel end

Base type for activity models within a solid solution phase. Each concrete subtype implements _excess_ln_gamma(model, k, x, T).

struct IdealSolidSolutionModel <: AbstractSolidSolutionModel

Ideal (Temkin) solid solution mixing: ln γᵢ = 0, hence ln aᵢ = ln xᵢ.

Valid for any number of end-members.

Example

julia> m = IdealSolidSolutionModel()
IdealSolidSolutionModel()

struct RedlichKisterModel{T<:Real} <: AbstractSolidSolutionModel

Binary Redlich-Kister (asymmetric Margules) model for non-ideal solid solutions. Requires exactly 2 end-members per solid solution phase.

Parameters a0, a1, a2 are in J/mol and divided by RT inside the activity closure.

Activity coefficients (Guggenheim / ThermoCalc convention):

ln γ₁ = (x₂²/RT)[a₀ + a₁(3x₁ − x₂) + a₂(x₁ − x₂)(5x₁ − x₂)]
ln γ₂ = (x₁²/RT)[a₀ − a₁(3x₂ − x₁) + a₂(x₂ − x₁)(5x₂ − x₁)]

AD-compatible: all computations propagate ForwardDiff.Dual numbers.

Examples

julia> m = RedlichKisterModel(a0 = 4000.0, a1 = 500.0)
RedlichKisterModel{Float64}(4000.0, 500.0, 0.0)

julia> m.a0
4000.0

abstract type AbstractSolidSolutionPhase end

Base type for solid solution phases. Concrete subtypes group a set of end-member AbstractSpecies and associate them with an AbstractSolidSolutionModel.

struct SolidSolutionPhase{T<:AbstractSpecies, M<:AbstractSolidSolutionModel}
        <: AbstractSolidSolutionPhase

A solid-solution phase consisting of end_members (species with AS_CRYSTAL aggregate state) mixing according to model.

End-members are automatically requalified to SC_SSENDMEMBER at construction time, so database species with SC_COMPONENT can be passed directly.

Construction

Use the keyword constructor:

SolidSolutionPhase(name, end_members; model = IdealSolidSolutionModel())

Validation at construction time:

• All end-members must have aggregate_state == AS_CRYSTAL.
• RedlichKisterModel requires exactly 2 end-members.

Example

julia> em1 = Species("Ca2SiO4"; aggregate_state=AS_CRYSTAL, class=SC_COMPONENT);

julia> em2 = Species("Ca3Si2O7"; aggregate_state=AS_CRYSTAL, class=SC_COMPONENT);

julia> ss = SolidSolutionPhase("CSH", [em1, em2])
SolidSolutionPhase{Species{Int64}, IdealSolidSolutionModel}
  name: CSH
  end-members (2): Ca2SiO4, Ca3Si2O7
  model: IdealSolidSolutionModel

julia> class(end_members(ss)[1])
SC_SSENDMEMBER::Class = 5

EquilibriumProblem

Definition of a chemical equilibrium problem.

Fields

• b: conservation vector (elemental abundances).
• A: stoichiometric matrix (conservation matrix).
• μ: chemical potential function μ(n, p).
• u0: initial guess for species amounts.
• p: coefficients for the potential function (default: NullParameters).
• lb: lower bounds for species amounts.
• ub: upper bounds for species amounts.

The problem solves for the species distribution that minimizes the Gibbs energy subject to mass conservation constraints A * n = b.

struct EquilibriumSolver{F<:Function, S, V<:Val}

Encapsulates all fixed ingredients of a chemical equilibrium calculation: the potential function, the SciML solver, and the variable space.

Construct once, call repeatedly with different ChemicalState inputs.

Fields

• μ: chemical potential closure μ(n, p) -> Vector{Float64}.
• solver: any Optimization.jl-compatible solver (e.g. IpoptOptimizer()).
• variable_space: variable space — Val(:linear) or Val(:log).
• kwargs: solver keyword arguments forwarded to solve.

Examples

julia> cs = ChemicalSystem([
           Species("H2O"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLVENT),
           Species("H+";  aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE),
           Species("OH-"; aggregate_state=AS_AQUEOUS, class=SC_AQSOLUTE),
       ]);

julia> solver = EquilibriumSolver(cs, DiluteSolutionModel(), IpoptOptimizer());

julia> solver isa EquilibriumSolver
true

equilibrate(state::ChemicalState;
            model    = DiluteSolutionModel(),
            solver   = IpoptOptimizer(...),
            variable_space  = Val(:log),
            ϵ        = 1e-16,
            kwargs...) -> ChemicalState

Compute the chemical equilibrium state from an initial ChemicalState.

This is a high-level convenience function that wraps EquilibriumSolver and solve with sensible defaults. It is intended for users who do not need to fine-tune the solver.

Arguments

• state: initial ChemicalState — defines the system, T, P, and composition.
• model: activity model (default: DiluteSolutionModel()).
• solver: SciML-compatible solver (default: IpoptOptimizer with tight tolerances).
• variable_space: variable space — Val(:linear) or Val(:log) (default). Val(:log) is more robust for systems spanning many orders of magnitude.
• ϵ: regularization floor for mole amounts (default: 1e-16).
• kwargs...: additional keyword arguments forwarded to the solver.

Returns

A new ChemicalState at thermodynamic equilibrium, with all derived quantities (pH, pOH, volumes, porosity, saturation) already computed.

Examples

# Minimal usage — all defaults
state_eq = equilibrate(state)

# Custom temperature-dependent run
set_temperature!(state, 350.0u"K")
state_eq = equilibrate(state)

# Custom solver tolerance
state_eq = equilibrate(state; abstol=1e-12, reltol=1e-12)

# Custom activity model (future)
state_eq = equilibrate(state; model=DebyeHuckelModel(A=0.51, B=3.28))