Default IMBNUM to SATNUM when not in deck

[dlmorenob]

Summary

When a deck activates relperm hysteresis (SATOPTS HYSTER + EHYSTR) but omits the IMBNUM keyword, OPM initialises every cell's imbibition region to 1. In decks with multiple SATNUM regions this silently selects the wrong imbibition tables for every cell whose SATNUM is not 1, producing a non-zero deltaSwImbKrn_ in the hysteresis scanning curve and a systematic mass-balance drift.

The sign and magnitude of the drift depend on how much of the grid has SATNUM != 1, how much of that volume is swept by gas, the drainage/imbibition endpoint spread, and the simulation horizon. Single-SATNUM decks are unaffected. Multi-SATNUM decks fail silently: no warning, no error, just a running mismatch in FGMIP vs WGMIT that grows over time.

Fix: default IMBNUM = SATNUM when the deck never mentions IMBNUM — the Eclipse-documented idiom for "switch hysteresis off in specific blocks", applied globally. A single 15-line change in FieldProps::FieldProps(...) fixes every downstream consumer (init_satfunc, EclEpsGridProperties, EclMaterialLawInitParams, ...) without touching any of them.

How to use hysteresis after this fix

Three patterns are all valid. IMBNUM can now be omitted in every case.

1. Full I-prefixed saturation tables. Classic setup: one SGFN/SWFN per SATNUM region plus matching ISGFN/ISWFN. Unchanged by this PR.

2. ENDSCALE + scaled imbibition endpoints (no I-tables). One drainage set per region, with the imbibition curve built analytically from the per-cell residual endpoint. Recommended when only the residual-Sg endpoint is known — typical of CO2-storage studies.

RUNSPEC
  SATOPTS
    HYSTER /
  ENDSCALE
  /

PROPS
  SGFN
    -- region 1 drainage (e.g. seal)
    ... /
    -- region 2 drainage (e.g. reservoir)
    ... /

  SWFN
    ... /
    ... /

  EHYSTR
    -- item 2: 0 or 1 = Carlson, 2 or 3 = Killough
    0.1  0  0.1  1*  KR  1*  1*  1* /

  -- Per-region imbibition critical gas saturation.
  -- Must exceed the drainage Sgcr of the same region.
  EQUALS
    'ISGCR' 0.30  1 80  1 80  1  5 /
    'ISGCR' 0.15  1 80  1 80  5 14 /
  /

REGIONS
  SATNUM
    ... /
  -- IMBNUM intentionally omitted; post-fix, it defaults to SATNUM.

3. No imbibition data at all. IMBNUM = SATNUM collapses the scanning curve onto the drainage curve — the simulator behaves as if SATOPTS HYSTER were absent. Mass is conserved.

For Pattern 2, additional scaled endpoints (ISGU, ISWCR, ISWU, ISWL) may be supplied per region via the same EQUALS mechanism if finer control is needed, but ISGCR alone is sufficient for residual gas trapping.

Worked example (not a universal symptom)

Reproducer: 80x80x14 CO2STORE deck with two SATNUM regions (shale seal over sand reservoir, ~72% of cells in region 2), CO2 injection into region 2 for 30 years, then 370 years of post-injection.

Config on this deck	FGMIP vs WGMIT
Pre-fix	+109.34% @ y=170
Post-fix, no I-data (hysteresis-off control)	-0.04% @ y=399
Post-fix, Carlson + ISGCR=0.25 uniform	+2.23% @ y=399
Post-fix, Killough + ISGCR=0.25 uniform	+0.60% @ y=399
Post-fix, Carlson + per-region ISGCR (R1=0.30, R2=0.15)	+0.37% @ y=399

Other deck geometries and injection schedules will land at different magnitudes; the pre-fix blowup on this deck is illustrative, not a universal figure.

Known remaining drift (out of scope for this PR)

A residual sub-percent drift remains on the reproducer deck when hysteresis is active alongside DISGASW + DRSDTCON. Between end-of-injection (y=30) and end-of-simulation (y=399), with WGMIT flat in that window:

Delta(FGMGP) + Delta(FGMDS) > 0

e.g. Carlson + ISGCR=0.25: Delta(FGMGP) = -2.50e+10, Delta(FGMDS) = +2.65e+10, net +1.46e+09 lb of gas mass appears in a closed system. Magnitude on this deck: +0.4% (per-region ISGCR) to +2.2% (uniform ISGCR) over 369 years.

Bisection controls on the same reproducer deck:

Configuration	Drift @ y=399
Hysteresis effectively off + DISGASW on	-0.04% (no drift)
Hysteresis on + DISGASW + VAPWAT off	+0.017% (no drift; numerical noise)
Hysteresis on + DISGASW + VAPWAT on	+0.4% to +2.2%

The drift is therefore localised to the hysteresis-dissolution coupling, not to the scanning-curve math (hysteresis on its own is mass-conservative to within numerical noise). Suspected culprit: dissolution of trapped (residual) gas into the water phase adds dissolved mass without a matching debit on the gas-phase side of the accounting. Flagged here for follow-up; not addressed by this PR.

Tests

Two new cases in tests/parser/FieldPropsTests.cpp:

• IMBNUM_DEFAULTS_TO_SATNUM — multi-SATNUM deck without IMBNUM must produce IMBNUM == SATNUM cell-by-cell.
• IMBNUM_EXPLICIT_PRESERVED — deck-supplied IMBNUM must not be overwritten, even when SATNUM differs.

All 62 FieldPropsTests + 32 test_hysteresis cases pass.

Test plan

• ctest -R FieldPropsTests
• ctest -R test_hysteresis
• Reproducer deck verified (see table above)
• Hysteresis-alone conservation verified (Case E: +0.017% over 400 years)
• Broader opm-tests regression suite (none of the existing EHYSTR decks combine multi-SATNUM with absent IMBNUM, so results are unchanged by construction)

[totto82]

Thanks for the PR. Defaulting IMBNUM to SATNUM seems logical. I will ask @bska to review the code since I am not 100% sure if amending the props as you suggest has unintended consequences for restart and output.

[bska]

Thank you for your work.

I appreciate the sentiment, but I don't think this is the right solution. If a model wants hysteresis, then it's (supposed to be) a modelling error if IMBNUM is not provided in the model description. If it were me, I'd rather remove the default initial value of 1 from the IMBNUM region array descriptor (i.e., remove .init(1)) instead of defaulting IMBNUM to SATNUM (i.e., no hysteresis).

Then you'd get a runtime diagnostic saying that IMBNUM does not exist when the simulator tries to access the array. That diagnostic will likely be somewhat confusing, I don't remember the exact wording right now, but it will at least state that the IMBNUM array does not exist in the input. Then the modeller gets a choice for how to handle that situation, which could be something simple like

COPY
  SATNUM IMBNUM /
/

or defining IMBNUM in some other way.

[dlmorenob]

Reworked as suggested — force-pushed to the same branch:

• FieldProps.hpp: removed .init(1) from IMBNUM
• FieldPropsTests.cpp: kept the "explicit IMBNUM preserved" test, dropped the "defaults to SATNUM" case
• Reverted the FieldProps.cpp copy-logic

One follow-up thought on documentation: in a CO2 storage workflow you often move between hysteresis patterns over a field's lifecycle —early on, using just endpoints (ISGCR etc.) to estimate residual trapping for uncertainty quantification; later, once the field is better characterized, switching to a fully specified IMBNUM with full I-tables is customary. Working out how the mechanics worked was not obvious to me. Would it be useful to document the common hysteresis patterns in the manual?

[bska]

Thanks a lot for the updates. This is starting to look good to me now, save for one minor detail in the new unit test.

[bska]

Just for your information, the Boost.Test library has a helper macro, BOOST_CHECK_EQUAL_COLLECTIONS, that runs this loop for you and also checks that the collections have the same size. Using that macro you could write

BOOST_CHECK_EQUAL_COLLECTIONS(imbnum.begin(), imbnum.end(), expected.begin(), expected.end());

Another benefit of using that macro is that it will also report the corresponding index if one or more elements differ between the two collections.

[bska]

jenkins build this please

[dlmorenob]

Thanks for the tip...been learning a lot. :)

[bska]

Thanks for the tip...been learning a lot. :)

Happy to help!

That said, I now see that OPM Flow user manual documents the default value of IMBNUM = 1 and that there are three existing tests that depend on this behaviour and therefore fail as a result of this PR.

I personally would prefer to not have a default value for IMBNUM–i.e., to require that models explicitly define IMBNUM when using hysteresis–but as this would require changing documented behaviour I can't make the decision without broader discussion.

[dlmorenob]

Just to add practical context: for CO2 storage workflows using analytical hysteresis via EHYSTR + ENDSCALE with per-region endpoints (no full imbibition tables), the current IMBNUM=1 default assigns the wrong imbibition region to every cell where SATNUM != 1. On my test case this produced a silent mass-balance drift where trapped-gas mass more than doubled.

Making IMBNUM required would mean every multi-SATNUM deck using analytical hysteresis needs COPY SATNUM IMBNUM / just to get correct behavior — This invalids the purpose of having analytical models.

I understand that changing documented behavior needs broader discussion, but I'd argue the current default is actively harmful for multi-region hysteresis decks: it's a silent physics-level error rather than a deliberate design choice. Defaulting IMBNUM to SATNUM when the keyword is omitted fixes the drift, keeps analytical hysteresis usable out of the box, and still allows explicit IMBNUM for users who need different imbibition regions.