EOS Parameter Regression Workflow

Guide for fitting equation-of-state parameters to experimental data using NeqSim.

When to Use This Skill

• Fitting binary interaction parameters (kij) to VLE/LLE data
• Matching saturation pressure (bubble/dew point) to lab measurements
• Tuning C7+ characterization to PVT experiments (CME, CVD, DL)
• Validating EOS predictions against NIST or experimental data
• Selecting the best EOS for a specific fluid system

Regression Strategy

For reservoir material-balance handoff, use CalibratedPvtModel to freeze original and independently tuned fluids with sample/model/laboratory revision references. Compare complete held-out CCE/CVD/DLE and separator experiments before accepting the tuned model. Never split rows from one depletion experiment between tuning and validation. MaterialBalancePvtHistory requires an explicitly attributed depletion and separator path; constant-composition flash tables cannot be passed as CVD/DLE history. It checks alignment and provenance fields, not the external evidence's truth or experimental accuracy. See docs/pvtsimulation/reservoir_material_balance.md and hand off measured-history volume estimation to @field-development, keeping laboratory EOS calibration separate from reservoir/aquifer fitting.

For explicit finite uncertainty scenarios, hand the frozen laboratory model to ReservoirHistoryEnsemble rather than fitting EOS parameters against field volume. Its posterior parameterQuantiles(q) reports marginal lower-tail in-place volume and support assumptions in getParameterNames() order; pressureQuantiles(q) reports model-pressure uncertainty. Zero/one select the positive-mass support endpoints, including a tiny positive tail or explicit infinite capacity. Marginal quantiles are not one joint scenario, continuous posterior coverage, petroleum exceedance P10 or recoverable reserves. Retain grid/prior/EOS/history sensitivity and the independent evidence boundary in docs/pvtsimulation/reservoir_aquifer_uncertainty.md.

Step 1: Select Base EOS

Step 2: Gather Experimental Data

For reported gas-phase p/T/Z/molar-mass/density tuples, apply the optional PvtGasDensityConsistency gate before regression. All quantities must describe the same phase/state and molar-volume basis; pressure is absolute Pa, temperature K, molar mass kg/mol, density kg/m3. Freeze justified absolute bounds first and retain failing rows and source/page/table/row identities. Do not treat printed precision as measurement sigma or widen bounds to manufacture agreement. Passing same-sample evidence can be frozen in CalibratedPvtModel via its overload. See docs/pvtsimulation/laboratory_gas_density_consistency.md and its executable regression; the check is arithmetic compatibility, not independent accuracy.

For saturated-vapor-reference gas-condensate CVD, screen the complete observed experiment with CvdMaterialBalance before fitting. Use absolute Pa, constant K, liquid volume per initial saturated cell volume and cumulative moles per initial mole. Each withdrawal uses that endpoint vapor composition. Preserve source row IDs, smoothed/adjusted versus raw analyses, and separately attributed final liquid analysis bounds. Absent/unresolved/inconsistent liquid has no available composition; never normalize or clip failed inventories. Optional checkpoint attachment uses CalibratedPvtModel(..., gasDensityEvidence, cvdEvidence) and rejects a failed or wrong-sample report. Arithmetic consistency is not independent measured accuracy. See docs/pvtsimulation/laboratory_cvd_material_balance.md; the primary NS-1 calculated liquid is a numerical reference, while experimental CH4/C7+ remain incompatible with printed-rounding-only bounds. Wet, critical and other reference protocols are outside this gate, not silently extrapolated.

Required data types by priority:

1. Saturation pressure (bubble/dew point) at reservoir temperature — most critical
2. Liquid density at reservoir conditions — validates volume translation
3. GOR/Rs from separator tests — validates phase split
4. Viscosity — validates transport property correlations
5. Compositional data from CVD/DL — validates K-values

Step 3: Set Up Regression in NeqSim

// Create fluid
SystemInterface fluid = new SystemPrEos(273.15 + 100.0, 200.0);
fluid.addComponent("methane", 0.70);
fluid.addComponent("ethane", 0.10);
fluid.addComponent("propane", 0.05);
fluid.addComponent("n-heptane", 0.15);
fluid.setMixingRule("classic");

// Access binary interaction parameters
fluid.setMixingRule("classic");

// Get current kij
double kij = fluid.getPhase(0).getMixingRule()
    .getBinaryInteractionParameter(0, 3);  // methane-nC7

// Set modified kij
fluid.getPhase(0).getMixingRule()
    .setBinaryInteractionParameter(0, 3, 0.02);
fluid.getPhase(1).getMixingRule()
    .setBinaryInteractionParameter(0, 3, 0.02);

Step 4: Define Objective Function

// Calculate saturation pressure
ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.bubblePointPressureFlash(false);
double calcPsat = fluid.getPressure("bara");

// Compare with experimental
double expPsat = 250.0;  // bara (from lab)
double error = Math.abs(calcPsat - expPsat) / expPsat * 100.0;

Step 5: Iterate and Optimize

# Python regression loop using scipy
from scipy.optimize import minimize_scalar
from neqsim import jneqsim

def objective(kij_value):
    fluid = jneqsim.thermo.system.SystemPrEos(373.15, 200.0)
    fluid.addComponent("methane", 0.70)
    fluid.addComponent("n-heptane", 0.30)
    fluid.setMixingRule("classic")
    fluid.getPhase(0).getMixingRule().setBinaryInteractionParameter(0, 1, kij_value)
    fluid.getPhase(1).getMixingRule().setBinaryInteractionParameter(0, 1, kij_value)

    ops = jneqsim.thermodynamicoperations.ThermodynamicOperations(fluid)
    ops.bubblePointPressureFlash(False)
    calc_psat = fluid.getPressure("bara")

    exp_psat = 250.0
    return (calc_psat - exp_psat) ** 2

result = minimize_scalar(objective, bounds=(-0.05, 0.10), method='bounded')
print(f"Optimized kij = {result.x:.4f}")

PVT Experiment Matching

Constant Mass Expansion (CME)

// Match relative volume vs pressure
SystemInterface fluid = createCharacterizedFluid();
ThermodynamicOperations ops = new ThermodynamicOperations(fluid);

double[] pressures = {400, 350, 300, 250, 200, 150, 100};  // bara
double[] expRelVol = {0.95, 0.97, 0.99, 1.00, 1.05, 1.15, 1.35};

for (int i = 0; i < pressures.length; i++) {
    fluid.setPressure(pressures[i], "bara");
    ops.TPflash();
    fluid.initProperties();
    double calcVol = fluid.getVolume("m3");
    // Compare with expRelVol[i]
}

Constant Volume Depletion (CVD)

// Match gas Z-factor and liquid dropout
SystemInterface fluid = createCharacterizedFluid();
neqsim.pvtsimulation.simulation.ConstantVolumeDepletion cvd =
    new neqsim.pvtsimulation.simulation.ConstantVolumeDepletion(fluid);
cvd.setPressures(new double[]{350, 300, 250, 200, 150, 100});
cvd.setTemperature(373.15);
cvd.runCalc();
double[] calcZgas = cvd.getZgas();
double[] calcLiqVol = cvd.getRelativeVolume();

C7+ Characterization Tuning

TBP Fraction Setup

SystemInterface fluid = new SystemPrEos(273.15 + 100.0, 200.0);
fluid.addComponent("methane", 0.60);
fluid.addComponent("ethane", 0.08);
fluid.addComponent("propane", 0.05);
fluid.addComponent("n-butane", 0.03);
fluid.addComponent("n-pentane", 0.02);

// TBP fractions: (name, moleFrac, MW_kg/mol, density_g/cm3)
fluid.addTBPfraction("C7", 0.05, 92.0 / 1000, 0.727);
fluid.addTBPfraction("C8", 0.04, 104.0 / 1000, 0.749);
fluid.addTBPfraction("C9", 0.03, 119.0 / 1000, 0.768);
fluid.addPlusFraction("C10", 0.10, 200.0 / 1000, 0.830); // no "+" in the name -- breaks parsing

fluid.getCharacterization().configureLumping().model("standard").totalPseudoComponents(6).build();
fluid.getCharacterization().characterisePlusFraction();
fluid.setMixingRule("classic"); // MUST come AFTER characterisePlusFraction() -- see below

Ordering note (confirmed 2026-10-12 against a real NeqSim 3.24.0 build): characterisePlusFraction() can SPLIT a plus/TBP fraction into several lumped pseudo-components, changing the total component count. Setting the mixing rule BEFORE that split leaves the EOS's kij matrix sized for the pre-split count; the first flash that actually needs the matrix then throws ArrayIndexOutOfBoundsException: Index N out of bounds for length N inside EosMixingRuleHandler$ClassicSRK.getkij. Always characterize first, then set the mixing rule.

Lumping count note (confirmed 2026-10-12, PR review on #4303, reproduced with a failing JUnit assertion): do NOT call fluid.getCharacterization().getLumpingModel().setNumberOfLumpedComponents(n) directly on the default ("standard") lumping model to control the total pseudo-component count -- StandardLumpingModel.generateLumpedComposition() overwrites numberOfLumpedComponents with numberOfPseudocomponents (both default to 7) at generation time, so setNumberOfLumpedComponents(6) alone still produces 7 pseudo-components (12 total components here: 5 real + 4 TBP/plus heavy fractions lumped into 7), not 6 (11 total). Use configureLumping().model("standard").totalPseudoComponents(n).build(), which sets numberOfPseudocomponents directly -- the field the generator actually reads.

Tuning Parameters

Validation Checklist

• Saturation pressure within ±2% of experimental
• Liquid density within ±3% of measured
• GOR within ±5% of separator test
• Z-factor within ±2% of CVD data
• Oil formation volume factor within ±3%
• Viscosity within ±20% (transport properties have larger uncertainty)
• Gas condensate: dew point within 10 % and C7+ density within 3 %; CME/CVD not degraded by the tuning
• Report the result as a pvt_tuning_quality block (before = untuned template, after = tuned, parameters with bounds and at-bound flags, excluded data with reasons): see neqsim-professional-reporting Principle 9a
• Lab data checked for physical validity first: a gas "CME" that stays near 1.0 was run on separator liquid

Common Pitfalls

1. Overfitting: Don't tune more parameters than you have independent data points
2. Non-unique solutions: Multiple kij sets can match Psat — validate with additional data
3. Temperature extrapolation: Fitted parameters may not work at different temperatures
4. Compositional shift: Parameters fitted to one composition may fail for depleted fluids
5. Forgetting setMixingRule: Must call before accessing kij parameters
6. Wrong units: MW in kg/mol (not g/mol) for addTBPfraction

Mechanics of Applying Tuning (read before writing a regression loop)

Iterate getPhases(), never getPhase(i)

getPhase(i) resolves through the phase-index map. On a single-phase system getPhase(0), getPhase(1) and getPhase(2) can all return the same phase object, while other phase objects the flash will later use are never touched. Always walk the raw array:

for phase in fluid.getPhases():
    if phase is None:          # the array has 6 slots, several are null
        continue
    c = phase.getComponent(i)
    c.setTC(c.getTC() * tc_mult)
    c.setPC(c.getPC() * pc_mult)

Regress in the representation you export

Applying parameters fitted on a 40-component characterization to a 22-component lumped fluid moved the dew point 12%. Run the regression on the lumped fluid in its own right and the lumping error disappears by construction. Round-trip the exported E300 file and re-check Psat and density — that is the only proof the delivered file is the model you fitted.

Do not regress a parameter the data cannot see

On a gas condensate, every CME density point is single-phase gas where the heavy pseudo-components are ~1 mol%; their Péneloux volume shift moves mixture density by <0.02%. A regression drives it to whichever bound it is given while barely changing the objective — and that arbitrary value then dominates the condensate liquid density in the reservoir model. Check the objective across the bound range before including a parameter; if it is monotonic to the bound, drop it and keep the correlation value.

Plus-fraction back-out fails on lean gas

Backing MW(C36+) out of the C7+ mass balance works for a black oil but not when the plus fraction is 0.001–0.006 mol%: the residual is smaller than the reporting precision and the molecular weight comes out negative. Instead keep the Katz SCN shape and rescale it with one MW factor and one SG factor to reproduce the measured mw7p and den7p — both factors are simply lab/Katz.

Saturation-pressure flashes on near-critical fluids

dewPointPressureFlash() and dewPointPressureFlashHC() can return the initial guess plus a small offset, so the "saturation pressure" silently tracks whatever start value the regression supplied. Verify against a pressure scan before trusting them. A reliable fallback is bisection on the pressure where TPflash first produces a second phase — but clone the fluid for every evaluation, because a reused fluid object carries its previous converged K-values into the next flash and the answer then depends on the order the bracket is walked. Cross-check the result against calcPTphaseEnvelope.

Quantify the data noise floor before chasing residuals

Compare samples that are compositionally near-identical. If two such samples report saturation pressures 24% apart, no EOS will fit both; tuning harder just fits noise. Report bias and scatter separately, and stop when the scatter matches the data's own inconsistency.

Explicit initial-volume depletion handoff

DepletionPvtHistoryGenerator.replayNativeFixedCell preserves original/tuned checkpoints and fixed separator provenance for independently declared initial occupied-volume experiments. Its immutable reference is not a dew-point estimate or measured saturation normalization. The bounded gas-bearing phase family must remain unchanged; critical and other-phase states fail. Keep the ordinary replayNative laboratory protocol separate. See docs/pvtsimulation/depletion_pvt_generation.md and FixedCellDepletionHandoffTest; synthetic EOS consistency does not establish measured calibration, held-out experiment accuracy or recoverable reserves.

Explicit-reference fixed-cell replay refines accepted endpoints through the existing multiphase beta solver to 1e-12 equilibrium and relative volume targets. Preserve all phase-family, Gibbs and component gates; failed refinement rejects the transaction. The Gibbs trial must use the same stored component inventory: the caller reconstructs a feasible homogeneous-phase allocation on copies before comparing energies. Keep FixedCellFeasibleGibbsTest alongside the independent fixed-cell references; this fix does not relax energy/fugacity/balance bounds or qualify global stability. Its108 SRK/PR/volume-correction controls check analytic extensivity and fresh homogeneous construction only; do not treat them as new depletion-path admission or independent EOS/measurement accuracy. Preserve the frozen EOS_CONTROL_SCOPE.md. Precise restoration reuses its caller-qualified endpoint before the first volume check and reflashes after withdrawal; preserve FixedCellInitialStateTest's named zero-inventory and independent binary controls. Do not widen volume tolerances to conceal a repeated-flash inventory perturbation. Legacy CVD keeps its sequence. Do not weaken absolute inventory bounds or claim arbitrary-scale/measured accuracy. The original 216-state and intermediate 480-state matrices use identical frozen bounds.

For independently declared dry initial-volume experiments with oil appearance or loss, use DepletionPvtHistoryGenerator.replayNativeDryPhaseChangingFixedCell (core ConstantVolumeDepletion.runDryPhaseChangingFixedCellDepletion). Preserve FixedCellResult.allowsDryPhaseChanges() with the immutable experiment reference. The existing fixed-cell method retains its unchanged-family contract. Gas must remain; water/aqueous/other phases, critical-degenerate splits and unclosed paths fail transactionally. Every local refinement retains its phase-family/Gibbs/ component gates; between-endpoint dry oil appearance/disappearance is permitted. DryFixedCellPhaseTransitionTest and the fixed-cell-transitions reference qualify 27dry binary SRK paths at three scales/both checkpoint selections, including the gas-only pressure-grid limit and every separator inventory. This is independent numerical consistency, not measured calibration, reservoir prediction or Tran RCVD. Read docs/pvtsimulation/depletion_pvt_generation.md#dry-phase-changing-fixed-cell-experiment.

Explicit wet-gas/aqueous fixed-cell path

Use DepletionPvtHistoryGenerator.replayNativeWetGasFixedCell (core ConstantVolumeDepletion.runWetGasFixedCellDepletion) only for a separately accepted initial-occupied-volume gas/aqueous experiment. Aqueous liquid may disappear; aqueous appearance/reappearance is rejected; gas must remain and oil/other phases fail transactionally. Water carried in produced gas is conserved through the unchanged separator. Preserve checkpoint selection, whole-experiment identity and actual final-product basis. allowsWetGasPhaseChanges() is distinct from the dry protocol flag. No wet-dew inference, saturation-normalized CVD, black-oil extrapolation, reservoir-volume fitting or measured-accuracy claim follows from this route. Read WetGasFixedCellHandoffTest and the frozen independent reference at src/test/resources/neqsim/pvtsimulation/wet-gas-fixed-cell/. Read OFFGRID_SCOPE.md and OFFGRID_VALIDATION.md in that reference folder for 27 additional 20-pressure paths, unchanged bounds and separately labeled current-master compatibility. Numerical case counts are not independent experimental measurements or laboratory/reservoir holdout evidence. Guide: docs/pvtsimulation/depletion_pvt_generation.md#wet-gas-fixed-cell-aqueous-transitions.