Overview

The electrolyte CPA (Cubic Plus Association) model in NeqSim extends the standard CPA equation of state to handle aqueous electrolyte solutions. The model is based on the work of Solbraa (2002) and combines:

  1. CPA equation of state for non-electrolyte interactions (van der Waals + association)
  2. Fürst electrostatic contribution for ion-ion and ion-solvent interactions
  3. Short-range Wij parameters for specific ion-solvent and ion-ion correlations

Implementation Classes

The Statoil (now Equinor) implementation of the electrolyte CPA model:

import neqsim.thermo.system.SystemElectrolyteCPAstatoil;

SystemElectrolyteCPAstatoil system = new SystemElectrolyteCPAstatoil(298.15, 1.01325);
system.addComponent("water", 55.5);
system.addComponent("Na+", 1.0);
system.addComponent("Cl-", 1.0);
system.chemicalReactionInit();
system.createDatabase(true);
system.setMixingRule(10);  // Required: CPA mixing rule with temperature/composition dependency

Reaction species and mixing-rule initialization order

chemicalReactionInit() can add species such as H3O+, OH-, HCO3-, and CO3--. Select the mixing rule after reaction initialization and createDatabase(true), as in the example above. The interaction matrices must cover the final component set.

If the rule was selected first and the component count subsequently changes, Fürst/electrolyte-CPA init(...) throws an IllegalStateException reporting the configured and current counts and the recovery steps. Direct electrolyte calcWij(...) calls also reject this stale state. Previously the reported four-to-eight-component case failed with an ArrayIndexOutOfBoundsException. The check covers species addition and removal in both Fürst phase implementations and their electrolyte CPA subclasses. It is not a global guard on ordinary cubic-EOS initialization: non-electrolyte component projection, characterization, and stream-union setup retain their existing initialization sequence.

To recover, complete reaction initialization, call createDatabase(true), and reapply the chosen mixing rule (10 for electrolyte CPA). Apply any custom interaction parameters after rebuilding the rule, then call init(0) before property calculations or a flash. The check does not silently reload fitted parameters or fill missing interactions with zeros. Changing only component amounts, or repeating chemistry initialization without adding species, does not require a mixing-rule refresh.

This count check does not detect same-count component replacement or renaming; rebuild the database and mixing rule after those changes as well. The separate chemical-reaction topology guard still requires chemicalReactionInit() again when component identities change after chemistry was initialized.

Regression coverage: MixingRuleComponentCountTest checks all five Fürst/CPA electrolyte variants, species removal, explicit recovery, unchanged initialized properties and inventories relative to a correctly ordered setup, and preservation of custom parameters when stale state is rejected.

Key Features of the Statoil Implementation

Feature Description
Model Name Electrolyte-CPA-EOS-statoil
Base Class Extends SystemFurstElectrolyteEos
Phase Class PhaseElectrolyteCPAstatoil
Component Class ComponentElectrolyteCPAstatoil
Attractive Term Term 15 (Mathias-Copeman alpha function)
Volume Correction Enabled by default
Fürst Parameters Uses electrolyteCPA parameter set

Ionic covolume selection and compatibility

CPA ionic components read covolume coefficients [0] and [1] directly from FurstElectrolyteConstants.furstParamsCPA, including construction and initFurstParam() reinitialization. This also applies to the Statoil and Advanced subclasses. ScRK components retain their separate furstParams defaults. Constructing an electrolyte CPA system no longer reassigns those shared ScRK defaults. Interleaving systems, adding ions later, cloning and serialization must therefore preserve each model’s covolume identity.

With the default tables, Na+ covolume is approximately 2.58461732544 for ScRK and 3.87214760778 for CPA in NeqSim internal units. Earlier versions could give the CPA value to ScRK solely because a CPA system had been constructed first. Removing that construction-order dependence can change affected ScRK results; it is not a new parameter fit or an estimate of fluid-property error.

For deliberate customization, use setFurstParamCPA(index, value) for the CPA table and setFurstParam(index, value) for the ScRK table, then reinitialize the affected components/phases. These remain process-wide mutable customization APIs, not per-system or thread-safe parameter stores. The legacy setFurstParams("electrolyteCPA") explicitly aliases the ScRK table to CPA and should not be used to select a system model. Code that relied on a CPA constructor making setFurstParam target CPA must migrate to setFurstParamCPA.

The covolume correction is independent of the reference-interaction correction described below. Neither changes fitted database pairs.

Class Hierarchy

SystemThermo
  └── SystemSrkEos
        └── SystemFurstElectrolyteEos
              └── SystemElectrolyteCPAstatoil
                    ├── PhaseElectrolyteCPAstatoil (phase calculations)
                    └── ComponentElectrolyteCPAstatoil (component properties)

Mixing Rule 10 - CPA with Temperature/Composition Dependency

Overview

Mixing rule 10 is the recommended mixing rule for all CPA and electrolyte CPA systems. It automatically selects the appropriate sub-type based on the binary interaction parameters:

system.setMixingRule(10);  // Automatically selects optimal sub-type

Automatic Sub-Type Selection

The mixing rule analyzes the binary interaction parameter matrices and selects:

Condition Sub-Type Class Description
Symmetric kij, no T-dependency classic-CPA ClassicSRK Simple symmetric mixing
Symmetric kij, with T-dependency classic-CPA_T ClassicSRKT2 Temperature-dependent symmetric
Asymmetric kij (kij ≠ kji) classic-CPA_Tx ClassicSRKT2x Full asymmetric + T-dependent

Mathematical Formulation

The a parameter mixing rule:

\[a = \sum_i \sum_j x_i x_j \sqrt{a_i a_j} (1 - k_{ij})\]

For asymmetric mixing (ClassicSRKT2x):

\[k_{ij} \neq k_{ji}\]

Temperature dependency:

\[k_{ij}(T) = k_{ij,0} + k_{ij,T} \cdot T\]

Why Mixing Rule 10?

  1. Automatic optimization: Selects the simplest sufficient mixing rule
  2. Temperature dependency: Captures T-dependent phase behavior
  3. Asymmetric parameters: Handles non-symmetric ion-solvent interactions
  4. Database integration: Uses binary parameters from NeqSim database

Theoretical Background

Helmholtz Energy Decomposition

The total residual Helmholtz energy is decomposed as:

\[A^{res} = A^{CPA} + A^{elec}\]

Where:

Fürst Electrostatic Model

The electrostatic contribution follows the Fürst model, which combines:

  1. Mean Spherical Approximation (MSA) for ion-ion interactions
  2. Born solvation term for ion-solvent interactions
  3. Short-range Wij terms for specific ion interactions
\[A^{elec} = A^{MSA} + A^{Born} + A^{SR}\]

Mean Spherical Approximation (MSA)

The MSA term accounts for the electrostatic screening between ions:

\[\frac{A^{MSA}}{RT} = -\frac{V}{3\pi} \left[ \Gamma^3 + \frac{3\Gamma\sigma_+ \sigma_-}{1 + \Gamma\sigma_{+-}} \right]\]

Where:

Born Solvation Term

The Born term accounts for the solvation energy of ions in the dielectric medium:

\[\frac{A^{Born}}{RT} = -\frac{e^2 N_A}{4\pi\varepsilon_0 RT} \sum_i n_i \frac{z_i^2}{\sigma_i} \left(1 - \frac{1}{\varepsilon_r}\right)\]

Where:

Here $n_i$ is in mol and $\sigma_i$ is the diameter in metres used by the Furst-based implementation. The equivalent radius convention uses $8\pi\varepsilon_0 RT$ and $r_i=\sigma_i/2$. Do not combine a diameter with the radius prefactor, or use $N_A/k_B$ in place of $N_A/R$ for this molar expression.

The mole derivative includes both the ionic contribution and the derivative of solvent permittivity. The latter must scale inversely with phase amount so that the chemical potential stays intensive. See Born derivatives and phase-size invariance for the chain rule, the legacy Mod2004 correction, the related EOS audit, and the literature-based development recommendations.

The base Furst/electrolyte CPA derivative implementation includes the implicit response of the stationary shielding parameter and the complete dielectric/Born second-derivative chain rules. See the same guide for finite-difference qualification and the remaining nitrogen, n-butane and n-hexane dielectric-data limitations. Those components share malformed default coefficients; replacing them requires species-valid property data and separate salting-out validation. A converged flash or consistent Jacobian does not by itself qualify the model’s salting-out or mixed-inhibitor accuracy.

Call system.init(3) before inspecting composition derivatives or running their consistency checks. A flash followed only by initProperties() may leave those arrays uninitialized. Derivative initialization refreshes electrolyte and CPA association caches at the final returned volume; repeated initialization and detached copies are covered by the final-volume regression tests. Component clones preserve fugacity composition-derivative caches in independently owned arrays. Electrolyte CPA component clones also own their association-composition derivative rows, so initializing one liquid phase cannot overwrite another phase’s component derivatives. This does not qualify concurrent access to all cloned phase workspaces.

Ionic shielding requires finite positive phase permittivity. If an extrapolated dielectric correlation violates that domain, the calculation raises an explicit error; ion-free states retain exactly zero shielding. In particular, a 1000 K endpoint for aqueous-phase disappearance in gas/oil/brine is invalid with the current solvent correlation. These model corrections do not qualify the legacy TP-flash boundary solver: brine boiling and gas/oil/brine pressure boundaries require the general equilibrium solver work for strict conservation, ion-confinement and fugacity-closure checks.

Short-Range Interaction Parameters (Wij)

Overview

The short-range Wij parameters capture specific ion-solvent and ion-ion interactions not described by the electrostatic terms. These are fitted to experimental activity coefficient and osmotic coefficient data.

Parameter Correlations

The reference coefficient wij[0] for calculated cation-water and ordinary cation-anion pairs follows the phase’s EOS family. SystemFurstElectrolyteEos and SystemFurstElectrolyteEosMod2004 use the existing ScRK furstParams table; the three electrolyte CPA variants use furstParamsCPA. Selection follows the phase type, including CPA subclasses, and is independent of construction order. Database pairs marked fitted (CalcWij != 0) retain their reference and temperature coefficients. The named-solvent tables, the Piperazine approximation, the MDEA+ ion-ion override, and gas-ion/organic-inhibitor overrides are unchanged.

For ScRK, the existing six-coefficient table gives the same correlation for monovalent and divalent cations, with diameters in angstroms:

Wij(cation-water, 298.15 K) = 6.99219e-5 * stokesDiameter + 4.3984e-6
Wij(cation-anion, 298.15 K) = -6.06e-8 * (d_cat + d_an)^4 - 2.1795e-5

These are the current ScRK table values, not a new fit. The separate ScRK/CPA model identity follows the correlations discussed by Solbraa (2002), equations 8.14-8.15 and table 8-12; the mutable repository tables remain the exact numerical source. The CPA-specific valence fits are described below and are unchanged.

Numerical compatibility: correcting the unconditional CPA dispatch introduced in #1787 changes calculated ScRK pairs. With the default Na+ Stokes diameter of 5.68 angstroms, wij[0](Na+, water) changes from 1.6297556353584585e-4 to 4.01554792e-4. ScRK phase equilibria and properties that depend on these pairs can therefore change. This is a model-selection regression correction (#3850), not evidence that the original or corrected calibration meets a given experimental accuracy target. CPA reference interactions retain their values.

The current shared furstParamsCPA_TDep temperature correction is deliberately retained for both families; the correction does not restore a complete historical ScRK calibration. Regressions check the reference coefficients, symmetric/reversed pairs, both valences, repeated initialization, CPA interleaving, fitted precedence, and the retained temperature formula at 298.15 and 323.15 K. Qualification of the combined ScRK correlation against experimental data, including temperature dependence, remains part of #3144.

The CPA Wij values are calculated using the following ionic-diameter correlations:

Monovalent (1+) Cations

Wij(cation-water) = furstParamsCPA[2] × stokesDiameter + furstParamsCPA[3]
Wij(cation-anion) = furstParamsCPA[4] × (d_cat + d_an)^4 + furstParamsCPA[5]

Current fitted values (2024):

Divalent (2+) Cations

Wij(2+ cation-water) = furstParamsCPA[6] × stokesDiameter + furstParamsCPA[7]
Wij(2+ cation-anion) = furstParamsCPA[8] × (d_cat + d_an)^4 + furstParamsCPA[9]

Current fitted values (refitted December 2024):

Why Separate Parameters for Divalent Cations?

The divalent cation parameters differ significantly from monovalent:

Using unified parameters would give dramatically wrong Wij values for divalent cations (sometimes even wrong sign), making separate parameters essential for accuracy.

Parameter Monovalent Divalent Ratio
Slope 4.98e-05 5.40e-05 1.08
Intercept -1.22e-04 -1.72e-04 1.42

Validation Against Experimental Data

Robinson & Stokes Data (1965)

The model has been validated against Robinson & Stokes experimental data for mean activity coefficients (γ±) and osmotic coefficients (φ) at 25°C.

Monovalent Salts (1:1 electrolytes)

Salt Type γ± Error φ Error
NaCl 1-1 2.4% 1.6%
KCl 1-1 4.3% 1.0%
LiCl 1-1 3.4% 2.5%
NaBr 1-1 2.8% 2.0%
KBr 1-1 1.4% 2.0%

Divalent Cation Salts (2:1 electrolytes)

Salt Type γ± Error φ Error
CaCl₂ 2-1 7.0% 4.2%
MgCl₂ 2-1 9.6% 4.6%
BaCl₂ 2-1 2.3% 1.5%

Divalent Anion Salts (1:2 electrolytes)

Salt Type γ± Error φ Error
Na₂SO₄ 1-2 20.0% 19.7%
K₂SO₄ 1-2 2.9% 1.6%

Overall Performance

Dielectric Constant Mixing Rules

The model supports three dielectric constant mixing rules:

\[\varepsilon_{mix} = \sum_i x_i \varepsilon_i\]

2. VOLUME_AVERAGE

\[\varepsilon_{mix} = \sum_i \phi_i \varepsilon_i\]

Where $\phi_i$ is the volume fraction.

3. LOOYENGA

\[\varepsilon_{mix}^{1/3} = \sum_i \phi_i \varepsilon_i^{1/3}\]

Thermodynamic Consistency

The model has been verified for thermodynamic consistency using built-in checks:

Fugacity Coefficient Identities

✅ $\sum_i x_i \ln\phi_i = \frac{G^{res}}{RT}$ - PASSED

Derivative Consistency

✅ $\left(\frac{\partial \ln\phi_i}{\partial P}\right)T = \frac{\bar{V}_i - V{ig}}{RT}$ - PASSED

✅ $\left(\frac{\partial \ln\phi_i}{\partial T}\right)P = \frac{H{ig} - \bar{H}_i}{RT^2}$ - PASSED

✅ $\left(\frac{\partial \ln\phi_i}{\partial n_j}\right){T,P,n{k\neq j}} = \left(\frac{\partial \ln\phi_j}{\partial n_i}\right){T,P,n{k\neq i}}$ (symmetry) - PASSED

Usage Example

// Create electrolyte CPA system
SystemInterface system = new SystemElectrolyteCPAstatoil(298.15, 1.01325);

// Add water (solvent)
system.addComponent("water", 55.5); // mol

// Add electrolyte
system.addComponent("Na+", 1.0);
system.addComponent("Cl-", 1.0);

// Initialize
system.chemicalReactionInit();
system.createDatabase(true);
system.setMixingRule(10); // Electrolyte CPA mixing rule

// Run flash calculation
ThermodynamicOperations ops = new ThermodynamicOperations(system);
ops.TPflash();

// Get activity coefficients
int aq = system.getPhaseNumberOfPhase("aqueous");
int naIdx = system.getPhase(aq).getComponent("Na+").getComponentNumber();
int clIdx = system.getPhase(aq).getComponent("Cl-").getComponentNumber();
int waterIdx = system.getPhase(aq).getComponent("water").getComponentNumber();

double gammaNa = system.getPhase(aq).getActivityCoefficient(naIdx, waterIdx);
double gammaCl = system.getPhase(aq).getActivityCoefficient(clIdx, waterIdx);
double meanGamma = Math.sqrt(gammaNa * gammaCl); // Mean activity coefficient

double phi = system.getPhase(aq).getOsmoticCoefficientOfWater();

Mixed Solvent Systems

The model supports mixed solvent systems including:

Separate Wij parameters are available for each solvent system.

For calculated cation-glycol pairs, TEG and the mixing-rule alias triethylene glycol use furstParamsCPA_TEG; MEG and ethylene glycol use furstParamsCPA_MEG. Each set uses indices [2]/[3] for the monovalent slope/intercept and [6]/[7] for divalent cations. The TEG parameters are initial estimates based on water parameters and have not been fitted to experimental TEG-water-electrolyte data.

The correction in issue #3846 removes an earlier MEG fallback that shadowed the TEG-specific branch. With the default parameters and the Na+ Stokes diameter of 5.68, the calculated Na+-TEG reference parameter is now 0.000160864, instead of the MEG value 0.0003394. TEG-containing calculations using the shared Furst short-range mixing rule can therefore change. Explicitly fitted pair parameters retain precedence, and the existing temperature-dependent coefficients are unchanged. This corrects parameter selection; it does not validate the TEG estimates against experimental data.

Gas-Ion Interaction Parameters

The electrolyte CPA model includes legacy gas-ion short-range interaction parameters. They are model coefficients in m3/mol, not measured Setschenow coefficients, and the original generic parameter provenance has not been independently reconstructed. Gas-ion behavior therefore requires a source-specific check rather than relying on the generic labels alone.

Setschenow Equation

\[\log\left(\frac{S_0}{S}\right) = k_s \cdot c_{salt}\]

Where:

SystemElectrolyteCPAstatoil has one source-reconstructed specialization: its default CO2/K+ coefficient reproduces the Onda et al. (1970) KCl ion-contribution correlation at 298.15 K. The calibration, coefficient units, override priority, independent holdouts, and unqualified regions are documented in CO2–KCl interaction provenance.

Usage Example - Methane Solubility with Salt:

SystemInterface fluid = new SystemElectrolyteCPAstatoil(298.15, 50.0);
fluid.addComponent("methane", 0.1);
fluid.addComponent("water", 10.0);
fluid.addComponent("Na+", 0.5);
fluid.addComponent("Cl-", 0.5);
fluid.setMixingRule(10);

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.TPflash();

// Methane solubility decreases with salt (salting-out effect)
double ch4InAqueous = fluid.getPhase(PhaseType.AQUEOUS).getComponent("methane").getx();

Organic Inhibitor-Ion Interaction Parameters

Hu-Lee-Sum Correlation for Hydrate Inhibition

For hydrate equilibrium calculations with combined salt + organic inhibitor systems, the Hu-Lee-Sum universal correlation (AIChE Journal 2017, 2018) states that water activity effects should be additive:

\[\ln(a_w^{combined}) = \ln(a_w^{salt}) + \ln(a_w^{OI})\]

This is critical for accurate hydrate inhibitor dosing calculations when both thermodynamic inhibitors (MEG, methanol) and formation water salts are present.

OI-Ion Interaction Parameters

Without explicit organic inhibitor-ion (OI-ion) parameters, the combined effect may not be additive. The following parameters ensure correct additive behavior:

Inhibitor W_MeOH-cation W_MeOH-anion Purpose
Methanol 1.5e-4 1.5e-4 Ensures MEG+salt gives more inhibition
MEG 0.0 0.0 Default calculation works correctly
Ethanol 1.3e-4 1.3e-4 Interpolated value

Validation Results

System Hydrate T (°C) Validation
Pure gas + water +20.0 Baseline
With NaCl only +10.8 ✅ Salt inhibition
With MEG only -2.3 ✅ MEG inhibition
With methanol only -4.6 ✅ Methanol inhibition
MEG + NaCl -18.9 ✅ ~16°C additional depression
Methanol + NaCl -5.4 ✅ ~0.8°C additional depression

Example - Combined Inhibitor Hydrate Calculation:

SystemInterface fluid = new SystemElectrolyteCPAstatoil(273.15 + 10.0, 100.0);
fluid.addComponent("methane", 0.85);
fluid.addComponent("water", 0.12);
fluid.addComponent("methanol", 0.03);
fluid.addComponent("Na+", 0.01);
fluid.addComponent("Cl-", 0.01);
fluid.setMixingRule(10);
fluid.setHydrateCheck(true);

ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.hydrateFormationTemperature();

// Combined effect is additive per Hu-Lee-Sum correlation
System.out.println("Hydrate T: " + fluid.getTemperature("C") + " °C");

Component conservation in hydrate-temperature calculations

Non-reactive brines retain the supplied molecular and ionic inventories. The multiphase solver temporarily uses a normalized, ion-free molecular feed for phase discovery. Stripped ions are excluded from its phase-fraction equations; their very small fugacity coefficients must not enter the molecular Newton matrix. Restoring the ions transforms both phase fractions and aqueous composition back to the full feed basis. Charge balance alone does not establish component conservation.

The final ionic gas/aqueous refinement compares Gibbs energies only when the reference state conserves the feed. A converged, normalized, conservative candidate must not be rejected because a state with a different inventory has a lower extensive Gibbs energy.

At every fluid evaluation, the hydrate-temperature operation checks phase and composition normalization, each component’s recovered inventory against the input, and aqueous ion confinement. An invalid non-reactive electrolyte state raises IllegalStateException with the failed diagnostic instead of returning a hydrate temperature. The caller’s multiphase-check setting is restored on success and failure. Reactive calculations retain their existing species and element-balance handling.

Regression coverage includes CO2 with NaCl–CaCl2 and NaCl–KCl mixtures on a 1 kg water / 10 mol CO2 basis, nearby pressures and salt concentrations, and equivalent salt-addition orders and repeated calculations. These are numerical conservation checks, not experimental qualification of mixed-salt hydrate temperatures.

For the CO2/water/explicit-ion subset, hydrate-temperature calculations also compare independently initialized vapour and liquid CO2 trials with a conserved aqueous feed. The constrained solver keeps ions in water, checks molecular fugacity equality for both gas/aqueous and liquid-CO2/aqueous splits, and accepts a single aqueous phase only after CO2 stability testing. This addresses the phase-state failures in issue #3584. It is a hydrate-temperature fluid-evaluation path; the generic TPflash() API and mixed-inhibitor/chemical paths are unchanged.

Use the operation’s getDiagnostics() snapshot to distinguish a numerically converged finite-inventory result from a saturated CO2 boundary. The hydrate phase-state guide documents the public API and the separate saturated/undersaturated experimental assessment. The selected public data give maximum absolute temperature errors of 1.12 K for saturated systems and 4.10 K for undersaturated systems, so this numerical correction does not establish full-range experimental accuracy or high-pressure drilling-fluid applicability.

Known Limitations

  1. Divalent anions (SO₄²⁻): Higher errors for 1:2 electrolytes like Na₂SO₄ (~20%)
  2. High concentrations: Accuracy decreases above ~2 mol/kg for some salts
  3. VOLUME_AVERAGE/LOOYENGA mixing rules: Incomplete composition derivatives
  4. Temperature range: Parameters fitted at 25°C, extrapolation accuracy may vary

Comparison: CPA vs Electrolyte CPA

Feature SystemSrkCPAstatoil SystemElectrolyteCPAstatoil
Use Case Non-ionic associating systems Aqueous electrolyte solutions
Electrostatics None MSA + Born solvation
Ions Not supported Na+, K+, Ca++, Mg++, Cl-, etc.
Mixing Rule 10 (recommended) 10 (required)
Chemical Reactions Optional Recommended (pH, speciation)
Phase Class PhaseSrkCPAs PhaseElectrolyteCPAstatoil

Quick Start Examples

Simple NaCl Solution

SystemInterface system = new SystemElectrolyteCPAstatoil(298.15, 1.01325);
system.addComponent("water", 55.5);
system.addComponent("Na+", 0.5);
system.addComponent("Cl-", 0.5);
system.createDatabase(true);
system.setMixingRule(10);
system.init(0);
system.init(1);

// Get mean activity coefficient
double gammaMean = system.getPhase(0).getMeanIonicActivityCoefficient("Na+", "Cl-");

With Chemical Reactions (pH Calculation)

SystemInterface system = new SystemElectrolyteCPAstatoil(298.15, 1.01325);
system.addComponent("water", 55.5);
system.addComponent("CO2", 0.01);
system.addComponent("Na+", 0.1);
system.addComponent("Cl-", 0.1);
system.chemicalReactionInit();  // Enable pH and speciation
system.createDatabase(true);
system.setMixingRule(10);

ThermodynamicOperations ops = new ThermodynamicOperations(system);
ops.TPflash();

// Access aqueous phase
int aq = system.getPhaseNumberOfPhase("aqueous");
double pH = -Math.log10(system.getPhase(aq).getComponent("H3O+").getx() * 55.5);

In-situ pH — getpH() and the acid-gas fallback

The recommended way to read the aqueous pH is getpH() on the aqueous phase (or SystemInterface.getpH(), which delegates to it). It supports several scales:

Method Definition Notes
getpH() / getpH("activity") -log10(γ_x · x_H3O+) default; mole-fraction activity, consistent with NeqSim’s reaction K’s
getpH("molality") -log10(γ_m · m_H3O+) IUPAC standard (mol H₃O⁺ / kg water)
getpH("molarity") -log10([H3O+]) mol/L, ignores activity coefficient
getpH("acidgas") carbonic/hydrosulfuric acid dissociation estimate screening in-situ pH from dissolved CO₂/H₂S

Acid-gas fallback (no chemicalReactionInit() required). The rigorous activity/molality/molarity scales need explicit H3O+ species, which only exist after chemicalReactionInit() has solved the water/carbonic dissociation equilibria. If those reactions were not initialised (or the electrolyte solver is numerically unstable and yields no/NaN H3O+, e.g. at low pressure), getpH() now falls back to an acid-gas dissociation estimate whenever the aqueous phase contains water together with dissolved CO₂ and/or H₂S — instead of silently returning a flat, unphysical 7.0.

The fallback uses the first-dissociation equilibria of carbonic and hydrosulfuric acid and an electroneutrality balance dominated by the conjugate bases (HCO₃⁻, HS⁻):

\[[\mathrm{H}^+] = \sqrt{K_1^{\mathrm{CO_2}}\,C_{\mathrm{CO_2}} + K_1^{\mathrm{H_2S}}\,C_{\mathrm{H_2S}} + K_w}, \qquad \mathrm{pH} = -\log_{10}[\mathrm{H}^+]\]

where $C_i$ are aqueous molar concentrations (mol/L) obtained from the aqueous phase density and molar mass. The constants are temperature-dependent: $K_1^{\mathrm{CO_2}}$ from Plummer & Busenberg (1982), $K_1^{\mathrm{H_2S}}$ from a van’t Hoff fit anchored at pK₁ = 7.05 (25 °C), and $K_w$ from a van’t Hoff fit anchored at 1.0×10⁻¹⁴ (25 °C). For CO₂-saturated water at ambient conditions this reproduces the textbook value pH ≈ 3.9; acid-gas-free water returns ≈ 7.0.

// Plain electrolyte fluid — NO chemicalReactionInit(), excess CO2 over water
SystemInterface system = new SystemElectrolyteCPAstatoil(298.15, 1.01325);
system.addComponent("CO2", 5.0);
system.addComponent("water", 10.0);
system.createDatabase(true);
system.setMixingRule(10);
system.setMultiPhaseCheck(true);

ThermodynamicOperations ops = new ThermodynamicOperations(system);
ops.TPflash();
system.initProperties();

double pH = system.getpH();   // ≈ 3.9 for CO2-saturated water (acidic)

Screening scope. The acid-gas fallback ignores alkalinity from other ions (bicarbonate buffering, dissolved salts) and the second dissociation steps. It is intended for acid-gas corrosion screening (e.g. NORSOK M-506), not for buffered-brine speciation. For a fully rigorous speciated pH in a buffered brine, run chemicalReactionInit() so explicit H3O+/HCO3-/OH- species are solved. For corrosion work that must always return a finite, source-tagged value, see neqsim.process.corrosion.RobustAqueousPH.

Gas-Liquid Equilibrium with Electrolytes

SystemInterface system = new SystemElectrolyteCPAstatoil(323.15, 50.0);
system.addComponent("methane", 10.0);
system.addComponent("water", 100.0);
system.addComponent("MEG", 20.0);
system.addComponent("Na+", 1.0);
system.addComponent("Cl-", 1.0);
system.createDatabase(true);
system.setMixingRule(10);

ThermodynamicOperations ops = new ThermodynamicOperations(system);
ops.TPflash();

// Check phase compositions
for (int i = 0; i < system.getNumberOfPhases(); i++) {
    System.out.println("Phase " + i + ": " + system.getPhase(i).getType());
}

References

  1. Solbraa, E. (2002). “Measurement and Modelling of Absorption of Carbon Dioxide into Methyldiethanolamine Solutions at High Pressures.” PhD Thesis, Norwegian University of Science and Technology.

  2. Fürst, W., & Renon, H. (1993). “Representation of excess properties of electrolyte solutions using a new equation of state.” AIChE Journal, 39(2), 335-343.

  3. Robinson, R.A., & Stokes, R.H. (1965). “Electrolyte Solutions.” 2nd Edition, Butterworths, London.

  4. Kontogeorgis, G.M., & Folas, G.K. (2010). “Thermodynamic Models for Industrial Applications.” Wiley.

  5. Michelsen, M.L., & Mollerup, J.M. (2007). “Thermodynamic Models: Fundamentals & Computational Aspects.” Tie-Line Publications.

  6. Hu, Y., Lee, B.R., Sum, A.K. (2017). “Universal correlation for gas hydrates suppression temperature of inhibited systems: I. Single salts.” AIChE Journal, 63(11), 5111-5124. DOI: 10.1002/aic.15868

  7. Hu, Y., Lee, B.R., Sum, A.K. (2018). “Universal correlation for gas hydrates suppression temperature of inhibited systems: II. Mixed salts and structure type.” AIChE Journal, 64(6), 2240-2250. DOI: 10.1002/aic.16generalized

Parameter History

Date Change Impact
2002 Initial parameters from Solbraa thesis Baseline model
2024 Refitted monovalent parameters to Robinson & Stokes γ± error: 2.8%
Dec 2024 Refitted divalent cation parameters [6-9] CaCl₂: 16%→7%, MgCl₂: 22%→10%
Dec 2024 Updated chemical equilibrium solver Improved pH accuracy
Dec 2024 Added gas-ion parameters for C2-C5+, N₂, H₂S, H₂ Correct salting-out for all gases
Feb 2026 Added OI-ion parameters for Hu-Lee-Sum compliance Additive hydrate inhibition with combined inhibitors
Jul 2026 getpH() acid-gas dissociation fallback Acidic in-situ pH for CO₂/H₂S water without chemicalReactionInit(); no more silent 7.0

Source Code References

System Classes

Phase Classes

Component Classes

Mixing Rules

Parameters

Tests

Last updated: December 27, 2024