Equilibrium

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.

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.

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 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.

DaviesActivityModel(; A=0.5114, b=0.3, bₙ=0.1, temperature_dependent=false)

Construct a DaviesActivityModel.

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())

HKFActivityModel(; A=0.5114, B=0.3288, Ḃ=0.041, Kₙ=0.1, å_default=3.72,
                   temperature_dependent=false) -> HKFActivityModel

Construct an HKFActivityModel with the given parameters.

Default values are from Helgeson et al. (1981), Table 1, at 25 °C / 1 bar.

_excess_ln_gamma(model, k, x, T) -> Real

Return the excess log-activity coefficient ln γₖ for end-member k (1-based index) of a solid solution with mole-fraction vector x at temperature T (K).

AD-compatible: all branches preserve ForwardDiff.Dual through computations on x.

Methods:

• IdealSolidSolutionModel: returns zero(eltype(x)).
• RedlichKisterModel: binary Redlich-Kister formula (requires length(x) == 2).

_hkf_sigma(x) -> Real

Compute the σ function used in the osmotic coefficient formula (Helgeson et al. 1981, Eq. 132–137):

σ(x) = (3/x³)(x − 2 ln(1+x) − 1/(1+x) + 1)

For |x| < 1e-3, a Taylor series 1 − (3/2)x + (9/5)x² is used to avoid catastrophic cancellation. Both branches agree to O(x³) at the threshold, so the gradient is continuous.

AD-compatible: branching is on ForwardDiff.value(x), not on the Dual itself.

_solid_solution_lna!(out, _n, ss_groups, ss_models, T, ϵ)

Fill out[i] with ln aᵢ = ln xᵢ + ln γᵢ for all solid-solution end-members.

ss_groups[k] and ss_models[k] describe the k-th solid-solution phase. T is the temperature in K (only relevant for non-ideal models). ϵ is a regularisation floor to avoid log(0).

ForwardDiff-compatible.

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.

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.

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

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).

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 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

RedlichKisterModel(; a0=0.0, a1=0.0, a2=0.0)

Keyword constructor for RedlichKisterModel. Parameters are promoted to a common type.

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

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

Construct and validate a SolidSolutionPhase.

End-members whose class is not already SC_SSENDMEMBER are automatically requalified via with_class. Passing database species with SC_COMPONENT therefore works directly, without a prior call to with_class.

end_members(ss::SolidSolutionPhase) -> Vector{<:AbstractSpecies}

Return the end-member species of the solid solution phase.

model(ss::SolidSolutionPhase) -> AbstractSolidSolutionModel

Return the activity model of the solid solution phase.

name(ss::SolidSolutionPhase) -> String

Return the name of the solid solution phase.

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: nothing).
• 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.

EquilibriumProblem(A, μ, u0; b=A*u0, p=nothing, lb=fill(Tu(1e-16), length(u0)), ub=maximum(abs.(A))/minimum(abs.(A[.!iszero.(A)]))*sum(u0)*one.(u0))

Construct an EquilibriumProblem with the given stoichiometric matrix A, chemical potential function μ, and initial guess u0.

Arguments

• A: stoichiometric matrix (conservation matrix).
• μ: chemical potential function μ(n, p).
• u0: initial guess for species amounts.
• b: conservation vector (elemental abundances). Defaults to A * u0.
• p: coefficients for the potential function. Defaults to nothing.
• lb: lower bounds for species amounts. Defaults to fill(Tu(1e-16), length(u0)).
• ub: upper bounds for species amounts. Defaults to maximum(abs.(A))/minimum(abs.(A[.!iszero.(A)]))*sum(u0)*one.(u0).

Returns

An EquilibriumProblem instance.

_DEFAULT_SOLVER_FACTORY

Internal Ref{Union{Nothing, Function}} — populated by extension __init__ functions to register a default solver factory.

• OptimizationIpoptExt.__init__: registers only if nothing is set (low priority).
• OptimaSolverExt.__init__: always overrides (high priority).

Result: OptimaSolver wins whenever loaded, regardless of load order.

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

EquilibriumSolver(cs, model, solver; variable_space=Val(:linear), kwargs...)

Construct an EquilibriumSolver from a ChemicalSystem, an activity model, and a SciML solver.

The potential function is built once at construction time from cs and model. Repeated calls to solve with different ChemicalState inputs reuse it.

Arguments

• cs: the ChemicalSystem defining species and conservation matrix.
• model: an AbstractActivityModel (e.g. DiluteSolutionModel()).
• solver: any Optimization.jl solver.
• variable_space: Val(:linear) (default) or Val(:log).
• kwargs...: forwarded to the underlying solve call (tolerances, verbosity...).

_build_n0(state::ChemicalState) -> Vector

Extract the dimensionless mole vector from a ChemicalState. Type is inferred from the state — compatible with ForwardDiff dual numbers.

_build_params(state::ChemicalState; ϵ=1e-16) -> NamedTuple

Extract dimensionless parameters from a ChemicalState. ΔₐG⁰overT is evaluated at the current T and P of the state. Units are stripped — compatible with ForwardDiff dual numbers.

Returned fields

• ΔₐG⁰overT: vector of standard Gibbs energies divided by RT (dimensionless).
• T: temperature in K (plain number, Dual-safe).
• P: pressure in Pa (plain number, Dual-safe).
• ϵ: regularisation floor (default 1e-16).

T and P are included so that temperature-dependent activity models (e.g. HKFActivityModel with temperature_dependent=true) can recompute their parameters inside the potential closure.

equilibrate(state::ChemicalState, solver; model=..., variable_space=..., ϵ=...) -> ChemicalState
equilibrate(state::ChemicalState; kwargs...) -> ChemicalState

Compute the chemical equilibrium state by minimising the Gibbs free energy.

Two-argument form (solver explicit, always available once an extension is loaded):

using Optimization, OptimizationIpopt
state_eq = equilibrate(state, IpoptOptimizer())

using OptimaSolver
state_eq = equilibrate(state, OptimaOptimizer())

One-argument form (uses the default solver of the loaded extension):

state_eq = equilibrate(state)

When both extensions are loaded, OptimaSolverExt always takes priority.

Arguments

• state: initial ChemicalState — defines the system, T, P, and composition.
• solver: any SciML-compatible solver (e.g. IpoptOptimizer(), OptimaOptimizer()).
• model: activity model (default: DiluteSolutionModel()).
• variable_space: Val(:linear) (default) or Val(:log).
• ϵ: regularisation floor for mole amounts (default: 1e-16).
• kwargs...: forwarded to the underlying solver.