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

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.

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 (Salting-Out Effect)

The electrolyte CPA model includes gas-ion interaction parameters to correctly predict the salting-out effect. The salting-out coefficient (k_s) determines how much dissolved gas solubility decreases with salt concentration.

Setschenow Equation

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

Where:

Fitted Gas-Ion Parameters

Gas k_s (L/mol) W_cation W_anion Notes
CO₂ 0.10 1.05e-4 1.05e-4 Polar, acidic gas
CH₄ 0.12 1.10e-4 1.10e-4 Reference hydrocarbon
C₂H₆ 0.13 1.13e-4 1.13e-4 Slightly larger
C₃H₈ 0.14 1.15e-4 1.15e-4 Larger hydrocarbon
C₄ (butanes) 0.15 1.20e-4 1.20e-4 Heavier hydrocarbon
C₅+ 0.16 1.25e-4 1.25e-4 Heaviest fractions
N₂ 0.10-0.12 1.05e-4 1.05e-4 Similar to CH₄
H₂S 0.06-0.08 1.10e-4 1.10e-4 Polar, acidic
H₂ 0.10 1.30e-4 1.30e-4 Very small molecule

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