feat(adv): add the ULTIMATE scheme of MT3DMS to the GWT and GWE ADV packages

[christianlangevin]

This PR adds a new ULTIMATE option to the SCHEME keyword of the GWT and GWE
ADV packages: an implementation of the third-order TVD (ULTIMATE) advection
scheme of MT3DMS. The face concentration is a QUICKEST reconstruction built
from the two cells straddling a face plus the second cell upstream, clamped by
the universal limiter to enforce the TVD condition. It resolves sharp fronts
with substantially less numerical dispersion than the existing TVD and UTVD
schemes and is exact for pure advection at a Courant number of one.

Formulation

ULTIMATE is explicit. The advective flux is evaluated entirely from the
concentration at the beginning of the time step, so ADV contributes nothing to
the coefficient matrix and everything to the right-hand side, giving an IMEX
scheme: advection explicit; storage, dispersion, reactions and sources
implicit as before. The face values depend only on the old dependent variable
and the step length, so they are computed once per time step through a new
prepare() hook on the interpolation scheme interface and stored per
connection, which guarantees both sides of a face use the same value and that
the mass budget closes. adv_cq reuses the same values.

The reconstruction is multidimensional: the face value is the average of the
field over the volume swept through the face during the time step, which
carries transverse-advection, transverse-curvature and twist (cross-derivative)
terms in addition to the QUICKEST terms normal to the face. The transverse
terms are what make the one-dimensional universal limiter usable in more than
one dimension: on uniform advection of a step front at 45 degrees to a 60x60
grid, the face-normal terms alone reach a maximum concentration of 2.03 while
satisfying the Courant check, whereas the full reconstruction stays within
[0, 1.04] across the reachable Courant range — within 0.01 of MT3DMS's own
overshoot — and is more accurate than MT3DMS at every step size measured.

Where sorption (GWT) or heat stored in the solid (GWE) retards the front, the
retardation factor enters the time step constraint and the reconstruction
alike, because the reconstruction traces the front backward over the step and
the front moves at the retarded velocity. ATS_PERCEL is therefore interpreted
as the fraction of a cell traversed by the retarded front. MST and EST
expose the retardation array and refresh it each time step. (MT3DMS
reconstructs from the unretarded velocity even though its update divides by
RETA; on a GWE problem with a thermal retardation of 5.3 this fix moves the
L1 error from 1.42, worse than TVD, to 1.05, better.)

Stability guard rails

The scheme is conditionally stable, so it is accepted only when

• the discretization is DIS (the stencil is two cells deep and axis-aligned),
• the simulation runs in serial (the scheme is programmed for serial
  simulations; parallel MPI simulations are rejected at input processing —
  an extension for parallel simulations could be implemented, by deepening
  the interface exchange to the two-cell stencil, but has not been
  programmed),
• ATS is active and ATS_PERCEL is specified with a value no greater than one.

Because ATS does not apply a submitted stability limit on the first step of a
stress period, prepare() also recomputes the maximum cell Courant number
from the step actually taken and terminates with the offending cell and a
suggested step length if it exceeds one.

Relation to MT3DMS

The transverse and twist terms are derived rather than transcribed. Expanded
as stencils, the derived form is algebraically identical to CFACE wherever
MT3DMS is self-consistent (positive transverse displacement, complete
stencils). It deliberately does not reproduce three defects in adv1.f, which
are catalogued in the source: locals with initializers carry an implicit SAVE
and reuse stale values next to inactive cells; the z-face branch assigns
TWISTY where it means TWISTZ, so the z-face twist term is whatever an earlier
face call left behind; and the upwinding correction to the transverse terms
does not follow the sign of the transverse flow, which is antidiffusive where
that flow is negative. MF6 zeroes any derivative whose stencil is incomplete,
computes the twist from the four diagonal neighbours in the plane of the face,
and derives the sign from the switch.

Accuracy

• 1-D advection of a step front: L1 error 1.31 against 3.14 for TVD/UTVD and
  6.30 for upstream; exact (L1 = 0) at C = 1.
• Against the Wexler analytical solutions (moc3d-p01a/b/c/d): lowest rms of
  the four schemes in all four cases, including the retardation case.
• Against MT3DMS ULTIMATE across the mt3dms example suite at identical time
  steps: closer than TVD in 9 of 11 problems, dramatically so where advection
  dominates (p01a 0.01% vs 7.12% rms of peak; p06 0.33% vs 9.17%). The two
  apparent exceptions are not advection results: p08 required the molecular
  diffusion the published MF6 model omitted (fixed in
  fix(ex-gwt-mt3dms-p08): use proper molecular diffusion modflow6-examples#346), after which ULTIMATE agrees to 0.7%,
  and p09's gap is a dispersion-formulation difference at a 1000x conductivity
  contrast — with dispersion off in both codes ULTIMATE reproduces MT3DMS to
  0.73%.

Documented limitations

• The scheme is not strictly monotone in more than one dimension; overshoot is
  bounded and comparable to MT3DMS (1.13 vs 1.11 worst-case on the oblique
  step front) and no undershoot is observed there. Making it rigorous means
  dimensional splitting, which is a different scheme.
• MF6 constrains the Courant number summed over a cell's outgoing faces, which
  it needs for the general case, so an n-dimensional problem runs at roughly
  1/n per direction and takes about n times as many steps as MT3DMS. Accuracy
  is flat in Courant number over the reachable range.
• A GWT/GWE model reading flows through FMI must be driven by a budget file
  with one time step per stress period, since ATS subdivides freely only then.
• Cells adjacent to a GWT-GWT/GWE-GWE exchange fall back to first-order
  upwind (the second upstream cell is outside the model).
• Vertical pass-through cells (IDOMAIN = -1) are supported via the
  six-direction neighbour map; a column containing one is geometrically
  identical to the same column without it and gives an identical answer.

Testing

• autotest/test_gwt_adv_ultimate.py: 1-D advection along each of the three
  grid axes at C = 0.5 and C = 1 (exact-shift), a two-cells-per-step variant,
  DTADJ = 0 (relies on the fix from fix(ats): honor a submitted stable time step when DTADJ is zero or one #2952), and a vertical pass-through case
  run at C = 0.25 so the upwind fallback would be visible if it were wrong.
• autotest/test_gwt_adv_ultimate_2d.py: uniform advection of a step front
  aligned and at 45 degrees to the grid, asserting error norms and bounds.
• The dfn changes ran through run-mf6ivar; release notes updated.

Checklist of items for pull request

• Replaced section above with description of pull request
• Referenced issue or pull request fix(ats): honor a submitted stable time step when DTADJ is zero or one #2952
• Added new test or modified an existing test
• Ran ruff on new and modified python scripts in .doc, autotests, doc, distribution, pymake, and utils subdirectories.
• Formatted new and modified Fortran source files with fprettify
• Added doxygen comments to new and modified procedures
• Updated meson files, makefiles, and Visual Studio project files for new source files
• Updated definition files
• Updated develop.toml with a plain-language description of the bug fix, change, feature; required for changes that may affect users

[Copilot]

Pull request overview

This PR adds a new ULTIMATE option to the SCHEME keyword of the GWT and GWE ADV packages — an explicit third-order TVD (QUICKEST + universal limiter) advection scheme ported from MT3DMS. It is a substantial addition (~700 lines of new Fortran implementing the scheme, plus a prepare() hook on the interpolation-scheme interface, retardation-factor support in MST/EST, and stability guard rails requiring DIS + serial + ATS + ATS_PERCEL ≤ 1). It builds on the ATS fix in #2952.

Changes:

• New UltimateSchemeType implementing an explicit, multidimensional QUICKEST reconstruction with a prepare() per-time-step hook, Courant verification, and universal limiter.
• Retardation-factor arrays and per-step refresh added to GWT MST and GWE EST; ADV now uses the retarded velocity for the ULTIMATE stability limit and reconstruction.
• Interface/enum/build/dfn/release-note plumbing plus two new autotests (1-D axis/pass-through/sorption and 2-D aligned/oblique).

Reviewed changes

Copilot reviewed 15 out of 15 changed files in this pull request and generated 1 comment.

Show a summary per file

File	Description
src/Model/TransportModel/InterpolationScheme/UltimateScheme.f90	New scheme: stencil map, derivatives, face-value reconstruction, limiter, Courant check.
src/Model/TransportModel/InterpolationScheme/InterpolationSchemeInterface.f90	Adds non-deferred prepare() default (no-op) hook.
src/Model/TransportModel/InterpolationScheme/AdvSchemeEnum.f90	Adds ADV_SCHEME_ULTIMATE = 4.
src/Model/TransportModel/tsp-adv.f90	New adv_ad/check_ultimate, adv_ar/adv_dt signature changes, retardation in time-step limit.
src/Model/GroundWaterTransport/gwt.f90 / gwe.f90	Pass porosity/retardation, call retardation refresh and adv_ad.
src/Model/GroundWaterTransport/gwt-mst.f90 / gwe-est.f90	New retardation array + *_calc_retardation and alloc/dealloc.
src/meson.build, msvs/mf6core.vfproj	Register new source file.
doc/ReleaseNotes/develop.toml, doc/mf6io/mf6ivar/dfn/gwt-adv.dfn, gwe-adv.dfn	Docs/dfn updates for the new scheme.
autotest/test_gwt_adv_ultimate.py, autotest/test_gwt_adv_ultimate_2d.py	New 1-D and 2-D tests.

I verified that all adv_ar/adv_dt callers were updated, interface models never select ULTIMATE (exchange scheme is independent), the retardation formulas are consistent with the MST/EST storage terms, check_ultimate terminates on failure, and this%xold matches the new pointer dummy argument. The one concrete issue found is the missing overshoot/undershoot assertion in the 2-D test (the property the test's own docstring says it verifies).

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.