Providing your own data and parameters¶

Every constant used in a calculation is a plain attribute you can set. There are five ways to supply data, from quickest to most structured:

Python DSL: declare species and reactions directly.
YAML or JSON specification: the same content as a file, validated strictly.
Database plus changes: start from a PHREEQC database, then add or override.
Your own PHREEQC-format file: reuse existing .dat data.
Fitting: estimate constants from measurements (see fitting.md).

1. Species and reactions in Python¶

from discopt.speciation import System

s = System("my system", activity="truesdell-jones")
Ca  = s.aqueous("Ca+2", a0=5.0, b=0.165)       # charge from the name; ion size a0 (Å) and b
SO4 = s.aqueous("SO4-2", a0=5.0, b=-0.04)
CaSO4 = s.aqueous("CaSO4")                     # neutral: log γ = b I, b defaults to 0.1
H2O = s.solvent("H2O")                         # water as a reacting species (optional)
gyp = s.solid("Gypsum")

s.reaction(Ca + SO4 >> CaSO4, logK=2.25, dH=5544.0, name="CaSO4 pair")          # J/mol
s.reaction("Gypsum = Ca+2 + SO4-2 + 2 H2O", name="Gypsum",
           analytic=[68.2401, 0.0, -3221.51, -25.0627])                            # PHREEQC A1..A6

Reaction forms. A reaction can be written with operators, as an equation string, or as a dict of net coefficients ({"Gypsum": -1, "Ca+2": 1, ...}).

Constraints. - Every solid and gas needs exactly one reaction linking it to non-phase species. - Reactions must be linearly independent. - A dependent reaction is rejected, because its logK would be an independent number that the thermodynamics does not allow.

Temperature and pressure dependence¶

\(\log K(T, P)\) uses the first form that is present:

data	formula
`analytic=[A1..A6]`	\(A_1 + A_2T + A_3/T + A_4\log_{10}T + A_5/T^2 + A_6T^2\)
`logK`, `dH`, `dCp`	\(\log K_r - \frac{\Delta H}{R\ln 10}(\frac1T-\frac1{T_r}) + \frac{\Delta C_p}{R\ln10}(\frac{T_r}{T}-1+\ln\frac{T}{T_r})\)
`dV`	adds \(-\Delta V (P - P_r)/(RT\ln 10)\)

Molar volumes. Above 1 atm, species molar volumes add the same term with \(\Delta V_r = \sum_i \nu_{ri} V_i(T, P, I)\), as in PHREEQC (see Concepts). Give them as PHREEQC's -Vm:

s.aqueous("Ca+2", Vm=[-0.3456, -7.252, 6.149, -2.479, 1.239, 5, 1.60, -57.1, -6.12e-3, 1])  # a1..a4 W b i1..i4
s.solid("Calcite", Vm=36.9)                                                                 # cm3/mol

Water's molar volume (\(18.016/\rho_w\)) counts once any species has Vm. A reaction's own dV is added to the volume term.

Tabulated data: fit analytic coefficients to your \(\log K(T)\) table with numpy.polyfit-style least squares, or fit logK and dH directly with fitting.md.

Other phases¶

s.gas("CO2(g)")
s.reaction("CO2(g) = CO2", logK=-1.468, dH=-19983.0)
s.gas_phase(["CO2(g)"], pressure=1.0)              # or volume=1.0 (L); or fix log10 P in equilibrate

s.exchange_site("X")                                # returns the virtual master X-
s.exchange("NaX"); s.exchange("CaX2")               # sites per formula from the name
s.reaction("Na+ + X- = NaX", logK=0.0)
s.reaction("Ca+2 + 2 X- = CaX2", logK=0.8)

s.surface_type("Hfo", ["Hfo_w", "Hfo_s"], area=600, mass=0.09)   # diffuse double layer
s.surface("Hfo_wOH"); s.surface("Hfo_wOH2+"); ...

s.solid_solution("(Ca,Sr)CO3", ["Calcite", "Strontianite"])       # ideal
s.solid_solution("Arag-Stront", ["Aragonite", "Strontianite"], guggenheim=[3.43, -1.82])  # non-ideal: stability.md
s.electron()                                         # e- for redox; then equilibrate(..., pe=4)

Gases with critical constants¶

A gas with a critical temperature Tc (K), critical pressure Pc (atm) and acentric factor omega follows the Peng–Robinson equation of state; one without them is ideal (see Concepts).

s.gas("Ar(g)", Tc=150.8, Pc=48.0, omega=0.0)          # when declaring it
s.critical_constants("CO2(g)", Tc=304.2, Pc=72.86, omega=0.225)   # or later, e.g. for a database gas
s.gas_phase(["CO2(g)", "Ar(g)"], pressure=100.0,
            kij={("CO2(g)", "Ar(g)"): 0.1})            # optional binary parameters (default 0)
s.gas_phase(["CO2(g)"], volume=1.0, eos="ideal")      # per phase: "ideal" or "peng-robinson"
System("name", gas_eos="ideal")                       # the whole system, fixed pressures included

Mixing. All members of a Peng–Robinson gas phase need constants. A phase that mixes gases with and without them raises an error, as in PHREEQC.
Defaults for kij. Pairs not given use the database's GAS_BINARY_PARAMETERS, then PHREEQC's built-in values for H₂O(g) with CO₂, H₂S, CH₄ and N₂ (Søreide and Whitson 1992).
Roots. A fixed-pressure gas takes the largest (vapor) root of the cubic. A fixed-volume gas is followed from the ideal gas, so it can become dense and liquid-like (CO₂ above ~65 atm at 25 °C), as in PHREEQC; a root with negative pressure is excluded.
Limits. Condensation into a separate liquid phase is not modelled. Inside the liquid–vapor loop PHREEQC caps the pressure at the loop maximum and the package does not (validation). The gibbs method and optimize do not support Peng–Robinson gases.

2. YAML / JSON specification¶

The same content as a strictly validated document (unknown keys are errors):

name: gypsum
activity: truesdell-jones
temperature: 298.15
species:
  - {name: Ca+2, a0: 5.0, b: 0.165}
  - {name: SO4-2, a0: 5.0, b: -0.04}
  - {name: CaSO4}
  - {name: Gypsum, phase: solid}
reactions:
  - {equation: "Ca+2 + SO4-2 = CaSO4", logK: 2.25, dH: 5544}
  - {equation: "Gypsum = Ca+2 + SO4-2", logK: -4.58, name: Gypsum}
problem:
  amounts: {Gypsum: 1.0}

from discopt.speciation import from_yaml, solve_spec
s = from_yaml("gypsum.yaml")      # a System
out = solve_spec("gypsum.yaml")   # JSON-able result dict (agents, scripts)
s.to_spec()                       # any System back to a spec dict

Sections. The spec also takes these sections:

exchange_sites
surfaces
solid_solutions
gas_phases
masters: input aliases such as {Ca: Ca+2}
species phase: solvent, gas, exchange, surface or electron
species params: extra activity parameters; for a gas, {Tc: 304.2, Pc: 72.86, omega: 0.225}
gas_eos (auto, ideal, peng-robinson), and per gas phase eos and kij: [[CO2(g), CH4(g), 0.1]]
reaction fields analytic, dCp and dV

Activity models. An activity model is given as a name or a dict, for example:

activity: {model: davies, A: 0.509}
activity: {model: sit, eps: [[Na+, Cl-, 0.03]]}

3. Start from a database, then add or override¶

db = PhreeqcDatabase.load("phreeqc.dat")
s = db.system(["Ca", "S(6)", "Na", "Cl"], phases=["Gypsum"])

gyp = next(r for r in s.reactions if r.name == "Gypsum")
gyp.analytic = None; gyp.logK = -4.60; gyp.dH = -460.0     # your own value

s.aqueous("CaCl+", a0=5.0, b=0.0)        # a species phreeqc.dat lacks (check "CaCl+" in s.species first)
s.reaction("Ca+2 + Cl- = CaCl+", logK=0.40, name="CaCl+")

s["Ca+2"].params["a0"] = 5.2             # change an ion size

Override rules. - Setting analytic = None is required when you override a reaction that has analytic coefficients, since those take precedence over logK. - Parameter changes take effect at the next equilibrate.

In YAML:

database: {file: phreeqc.dat, components: [Ca, "S(6)"], phases: [Gypsum]}
overrides:
  Gypsum: {logK: -4.60, dH: -460}      # clears analytic unless you give one
species:
  - {name: MyComplex}
reactions:
  - {equation: "Ca+2 + SO4-2 = MyComplex", logK: 0.3, name: MyComplex}

4. Your own PHREEQC-format database¶

Any file in PHREEQC format works. You can write a small one with only the species you need:

SOLUTION_MASTER_SPECIES
H       H+      -1.0    H       1.008
H(1)    H+      -1.0    0
E       e-      0       0.0     0
O       H2O     0       O       16.0
O(-2)   H2O     0       0
Ca      Ca+2    0       Ca      40.08
S(6)    SO4-2   0       SO4     32.064
SOLUTION_SPECIES
H+ = H+
    log_k 0
e- = e-
    log_k 0
H2O = H2O
    log_k 0
Ca+2 = Ca+2
    log_k 0
    -gamma 5.0 0.165
SO4-2 = SO4-2
    log_k 0
    -gamma 5.0 -0.04
H2O = OH- + H+
    log_k -14.0
    delta_h 13.362 kcal
Ca+2 + SO4-2 = CaSO4
    log_k 2.25
    delta_h 1.325 kcal
PHASES
Gypsum
    CaSO4:2H2O = Ca+2 + SO4-2 + 2 H2O
    log_k -4.58
    delta_h -0.109 kcal
END

db = PhreeqcDatabase.from_file("mine.dat")    # or PhreeqcDatabase(text)
s = db.system(["Ca", "S(6)"], phases=["Gypsum"])

Units. delta_h is read in kJ/mol unless followed by kcal (or J).

Gases. Give critical constants as in PHREEQC, and binary parameters in a GAS_BINARY_PARAMETERS block:

PHASES
Ar(g)
    Ar = Ar
    log_k -2.85
    -T_c 150.8; -P_c 48.0; -Omega 0.0
GAS_BINARY_PARAMETERS
CO2(g)  Ar(g)  0.1

Pitzer data. Add a PITZER block to your file, or build the parameters in Python:

from discopt.speciation import Pitzer, PitzerParameters

p = PitzerParameters()
p.add("B0", ["Na+", "Cl-"], [0.0765, -777.03, -4.4706, 0.008946, -3.3158e-6])  # a0..a5 (T dependence)
p.add("B1", ["Na+", "Cl-"], 0.2664)
p.add("C0", ["Na+", "Cl-"], 0.00127)
p.add("THETA", ["Na+", "K+"], -0.012)
p.add("PSI", ["Na+", "K+", "Cl-"], -0.0018)
s = System("brine", activity=Pitzer(p))

Temperature form. The Pitzer coefficient temperature form is PHREEQC's: \(P(T) = a_0 + a_1(1/T - 1/T_r) + a_2\ln(T/T_r) + a_3(T-T_r) + a_4(T^2-T_r^2) + a_5(1/T^2-1/T_r^2)\).

Other models. SIT takes SIT(eps={("Na+", "Cl-"): 0.03}). B-dot takes per-species a0 and a common bdot. See activity-models.md.

Checking your data¶

Structure: s.validate() checks for undeclared species, reactions without a constant, and solids or gases that do not have exactly one reaction. Linearly dependent reactions are rejected when equilibrating.
Result state: r.residuals reports the mass-action, balance and complementarity residuals of the result; r.converged is true when all are below 1e-8.
Sensitivity: perturb a constant and re-solve. fit.FitResult.jacobian holds these sensitivities for a set of observations.