Class GWMesh#

class flow.gwplex.GWMesh(*args, **kwargs)[source]#

Water table (groundwater) + generic duricrust — opt-in near-surface hydrology and chemical armoring. See docs/DESIGN_WATERTABLE_DURICRUST.md and the technical guide docs/tech_guide/groundwater.rst.

Enabled by the YAML groundwater: block (self.gwOn / self.duriOn); when absent, every path here is a no-op and goSPL is byte-identical to a run without it. The whole update runs once per step in updateGroundwater (after flow accumulation, before erosion), so each process reads the previous process’s state — a slow, stable, explicit/sequential coupling.

Water table. An implicit (backward-Euler) Dupuit-Boussinesq head solve on the DMPlex (_solveHead): (I + (Δt/S)·L(T))·h = h_old + (Δt/S)·R with transmissivity T = Ksat·max(h − z_bed, b_min), a seepage free boundary (h ≤ z at sea / lakes / open outlets) and a dry-aquifer floor (h ≥ z_bed). z_bed is prescribed (aquifer_base scalar/map) or tied to the bedrock (aquifer_base: from_soil ⇒ lHbed − bedrock_depth). Recharge R = f_infil·max(0, rain − evap) is zeroed under water and ice. Opt-in baseflow closure returns the seepage discharge to the rivers (conserve_baseflow).

Duricrust. A generic capillary-fringe crust (_updateDuricrust) grows where the water-table depth sits in the fringe band (Φ Gaussian on wt = z − h), fed by a weathering supply Ψ (climate proxy, an explicit Maher-Chamberlain rate, or the soil-production congruency; regolith-limited when soil is tracked). Its induration duriF ∈ [0,1] armors erodibility through the single _surfaceArmoringK hook (1 − armor_max·duriF), which slows erosion of crusted cells and drives relief inversion. When stratigraphy is on, the induration is archived per layer (stratDuri) so buried crusts are preserved and re-arm on exhumation (stacked-duricrust / cratonic laterite behaviour).

State & lifecycle. Persistent Vecs (headL/headG, duriHL/duriHG, rechargeL, baseflowL) and the cached elliptic solver (_ksp_gw, fgmres + hypre BoomerAMG) are registered in destroy_DMPlex; head and duriH are model memory written to / restored from the output HDF5 on restart. Compatible with dual lithology and provenance (independent per-layer fields; armoring composes at the shared K hook).

Methods

updateGroundwater()

Per-step groundwater / duricrust update.

Initialise

__init__(*args, **kwargs)

Initialise the groundwater / duricrust state.

Public Methods

updateGroundwater()

Per-step groundwater / duricrust update.

Private Methods

`_makeGWKSP`()	Cached KSP for the elliptic head solve: fgmres + hypre BoomerAMG (algebraic multigrid).
`_solveHead`()	Implicit (backward-Euler) Dupuit-Boussinesq water-table head solve (DESIGN_WATERTABLE_DURICRUST.md §2/§3).
`_gwZbed`(z)	Aquifer-base elevation `z_bed` (m).
`_surfaceSlope`()	Per-node steepest-descent slope (m/m, ≥ 0) from the flow-direction receivers built by `flowAccumulation` — `(z − z_rcv)/dist` to the primary receiver `rcvID[:,0]`.
`_baseflowClosure`(hold, hnew, seep)	Baseflow (seepage-return) accounting (DESIGN_WATERTABLE_DURICRUST.md §3 step 7), opt-in via `conserve_baseflow`.
`_lakeExchangeFlux`(hloc, zbed, bmin, zeroKp)	Signed groundwater flux exchanged with the surface at every node (m³/yr), stored in `self.gwLakeFlux` for the lake volume coupling (§15).
`_arrhenius`()	Optional Arrhenius temperature factor for the weathering supply, `exp(Ea/Rg · (1/T_ref − 1/T))` (dimensionless, 1 at `T = T_ref`).
`_weatheringSupply`()	Solute-supply rate `Ψ` feeding in-situ fringe precipitation (DESIGN_WATERTABLE_DURICRUST.md §3a).
`_updateDuricrust`()	Evolve the generic capillary-fringe duricrust over `self.dt` (DESIGN_WATERTABLE_DURICRUST.md §3, step 5).
`_regolithSupplyRate`()	Rate (m/yr) at which weathering-produced regolith can feed the duricrust, the cap on chemical crust formation when soil is tracked (§8): the soil production rate `prodSoil · rain` (chemical∝physical weathering congruency).
`_recordInduration`()	Sync the live per-node induration `duriF` with the per-layer stratigraphic archive `stratDuri` (DESIGN_WATERTABLE_DURICRUST.md §9).
`_surfaceSourceClass`()	Per-node current surface lithology class — the dominant provenance class of the top non-empty stratigraphic layer (from `stratP`), falling back to the static bedrock `source_class` where a column has no sediment.
`_resolveGeoWeather`()	Resolve the per-species weatherability to a list of scalar or ``(lpoints,)`` array entries (Level-B extension 1 — lithology → chemistry).
`_makeSoluteKSP`()	Cached KSP for the solute-transport solve: a direct LU factorisation (`preonly` + `lu`; MUMPS in parallel), `gw_solute_` prefix.
`_soluteAdvecCoeffs`(h, T)	Upwind finite-volume coefficients for the steady solute advection `∇·(q c)` by the groundwater flux `q = −T∇h` (area-normalised, per cell).
`_dischargeWeight`(divq, R)	Absolute-accumulation gate (opt-in `discharge_gate`): a per-node weight `G ∈ [0,1]` that restricts crust formation to groundwater discharge zones.
`_solveSoluteTransport`(source, dmask, dval)	Solve one steady tracer transport `M c = source` with Dirichlet nodes `dmask` pinned to `dval` (`M` = upwind advection of §``_soluteAdvecCoeffs`` at the current head).
`_soluteSolveRHS`(M, rhs, guess)	Solve `M c = rhs` for the (pre-assembled) transport operator `M` with a warm-start `guess`, returning the local concentration (clipped ≥ 0).
`_updateSolute`()	Level-B per-step solute update (G2), per tracer: dissolve → transport → precipitate → export, conservatively accounted.
`_routeRiverSolute`()	Route the groundwater-exported (seepage/baseflow) solute down the surface drainage network to the shoreline — the river dissolved load (Level-B extension 2, DESIGN_GEOCHEM_EXTENSIONS.md §2), per species.

Public functions#

GWMesh.updateGroundwater()[source]#

Per-step groundwater / duricrust update.

Recharge R = f_infil * max(0, rain - evap) (m/yr), held at 0 under standing water (marine seaID or a ponded continental lake pitIDs>-1 & lFill>hl) and under ice, where rain does not infiltrate. Stored in self.rechargeL for output.
Head solve (_solveHead): the implicit Dupuit-Boussinesq water table from that recharge; sets self.wtDepth = z - h.
Duricrust (_updateDuricrust, only when duriOn): evolve the capillary-fringe crust duriH and the induration duriF / armor multiplier duriKarmor from the new water-table depth.
Geochemistry (_updateSolute, only when gwGeochemOn — Level B): per tracer, dissolve → transport (down q = -T∇h) → precipitate at the fringe (feeding duriH, replacing the Level-A supply) → export; a closed mass budget with per-species crust typing, optional solute-source provenance and a dissolved ocean flux.
River routing (_routeRiverSolute, only when gwRiverLoad): route the exported solute down the surface drainage network to the coast (per species; optional in-transit loss + marine coupling).
Stratigraphic archive (_recordInduration, when duriOn and stratigraphy is on): record the induration degree and, with geochem, the crust’s dominant species / source per layer (exhumation re-arm + formation write-down).

No-op when gwOn is off. The head, solute-transport and river-routing solves are collective (KSP); recharge, duricrust and the archive are purely rank-local (per-node).

Private functions#

GWMesh._makeGWKSP()[source]#: Cached KSP for the elliptic head solve: fgmres + hypre BoomerAMG (algebraic multigrid). The head operator I + (Δt/S)·L(T) is a stiff, high-condition 2-D elliptic operator; block-Jacobi/ILU does NOT converge on it (it stalls at DIVERGED_ITS and returns a garbage iterate that drifts to the surface) — an elliptic solve needs a multigrid PC. Same rationale as the flexure biharmonic wanting a strong solver. Options prefix gw_ scopes any override (env -gw_pc_type ...).

GWMesh._solveHead()[source]#

Implicit (backward-Euler) Dupuit-Boussinesq water-table head solve (DESIGN_WATERTABLE_DURICRUST.md §2/§3). One elliptic solve per step:

\[(I + \tfrac{\Delta t}{S}\,L(T))\,h = h_{old} + \tfrac{\Delta t}{S}\,R\]

with transmissivity \(T = K_h\max(h - z_{bed}, b_{min})\) (unconfined), \(L(T)=-\nabla\cdot(T\nabla)\) the FV neg-Laplacian (from jacobiancoeff, area-normalised — the same operator the marine Picard solver uses). Reduces to the steady elliptic limit as \(\Delta t/S\) grows, so it is quasi-steady at century steps.

Unconfined non-linearity T(h): gwPicardIts Picard sweeps (lag T on the previous iterate; operator rebuilt per sweep).
Seepage free boundary h <= z: Dirichlet h = z at the base seepage set (marine seaID + ponded lakes + open outletIDs) via zeroRowsLocal; after each Picard block, clip h = min(h, z) and add any newly-over-topped cell to the Dirichlet set, up to gwSeepagePasses passes (break on an Allreduce’d new-seepage count — collective-safe). aquifer_base is prescribed (scalar/map) OR, with aquifer_base: from_soil, tied to the bedrock elevation z_bed = lHbed − bedrock_depth (Phase 5 soil coupling; _gwZbed).
Baseflow closure (opt-in conserve_baseflow): the seepage-return discharge is accounted into self.baseflowL (_baseflowClosure).

Scratch: uses self.tmp (global rhs); leaves it defined. Writes self.headL/headG and self.wtDepth = z - h.

GWMesh._gwZbed(z)[source]#

Aquifer-base elevation z_bed (m). Prescribed by aquifer_base (scalar or per-vertex map, z_bed = z − aquifer_base), OR — with aquifer_base: from_soil (Phase-5 soil coupling) — tied to the bedrock elevation z_bed = lHbed − bedrock_depth, the “permeable regolith over impermeable bedrock” model. from_soil needs soilSPL (lHbed); it falls back to the surface (a self-consistent no-aquifer floor) with a one-time warning if soil is off.

Depositional basins (`from_soil` + stratigraphy). Under Option-2.5 lHbed = z − Lsoil is only the base of the thin weathering regolith; deposited sediment lives in the stratigraphy, not Lsoil. So in a filled basin the porous sediment column is the aquifer and its base is the bottom of the (non-sentinel) stratigraphic pile, deeper than lHbed. The base is therefore taken as the deeper (lower) of lHbed − bedrock_depth and z − Σ sediment thickness — so uplands (thin pile) keep the regolith base while basins deepen to the fill base. Rank-local.

GWMesh._surfaceSlope()[source]#: Per-node steepest-descent slope (m/m, ≥ 0) from the flow-direction receivers built by flowAccumulation — (z − z_rcv)/dist to the primary receiver rcvID[:,0]. Flats / sinks / outlets (dist=0 or a self-receiver) return 0. A cheap proxy for the optional slope-dependent infiltration (goSPL stores no explicit slope). Rank-local.

GWMesh._baseflowClosure(hold, hnew, seep)[source]#

Baseflow (seepage-return) accounting (DESIGN_WATERTABLE_DURICRUST.md §3 step 7), opt-in via conserve_baseflow. Over the quasi-steady step the net recharge that does not go into aquifer storage discharges back to the surface-water network at the seepage nodes:

\[Q_{seep} = \sum_{owned}(R\,A) - \sum_{owned} S\,(h-h_{old})\,A/\Delta t\]

(m³/yr). The global budget is reduced across ranks (Allreduce) and distributed over the owned seepage nodes weighted by cell area, so Σ_owned baseflowL ≈ Σ_owned recharge in the steady limit (ΔS→0) — river discharge stays ≈ rain − evap. Stored in self.baseflowL (diagnostic/output); re-injection into the surface-flow source (making rivers physically baseflow-fed, which redistributes erosion) is the next increment. Collective (two Allreduce); no rank-local collective gate.

GWMesh._lakeExchangeFlux(hloc, zbed, bmin, zeroKp)[source]#: Signed groundwater flux exchanged with the surface at every node (m³/yr), stored in self.gwLakeFlux for the lake volume coupling (§15). It is the FV divergence of the lateral groundwater flow, ∇·(T∇h)·A = −(L·h)·A (L = the area-normalised neg-Laplacian from jacobiancoeff, so L·h = −∇·(T∇h)): positive = the aquifer discharges INTO the surface (a gaining lake), negative = the surface leaks INTO the aquifer (a losing lake). One extra operator assemble + mat-vec, only when lake_exchange is on. Uses the un-zeroed operator (the seepage-Dirichlet rows would otherwise null the flux exactly at the lake nodes we need).

GWMesh._arrhenius()[source]#

Optional Arrhenius temperature factor for the weathering supply, exp(Ea/Rg · (1/T_ref − 1/T)) (dimensionless, 1 at T = T_ref).

Returns 1.0 (no temperature dependence) when weather_Ea <= 0 OR no temperature map is available — so duricrust runs soil-free by default. Reuses the soil tempMap (self.tempFile/tempData/tempRef) when present; loaded once and cached. Rank-local.

GWMesh._weatheringSupply()[source]#

Solute-supply rate Ψ feeding in-situ fringe precipitation (DESIGN_WATERTABLE_DURICRUST.md §3a). Three modes (duriWeatherMode); all plug in at the same place, so the formation ODE is mode-agnostic. Supply-only (no mass debit) — see §3a Level-A caveat. Rank-local.

proxy (default): Ψ = max(0, rain − evap)^p · arrhenius(T) — a climate/temperature stand-in, no new inputs.
rate (Level A): Maher-Chamberlain kinetic×thermodynamic law W = R·C_eq·(1 − exp(−Dw/(R·L)))·arrhenius·weatherability driven by the recharge R the head solve already computes; L = Lsoil when soilSPL is on (capped by the regolith production supply), else the prescribed path_length.
prodsoil: reuse the temperature-scaled soilSPL.prodSoil directly, scaled by water availability (the chemical∝physical congruency). Falls back to the proxy when soil is not tracked.

GWMesh._updateDuricrust()[source]#

Evolve the generic capillary-fringe duricrust over self.dt (DESIGN_WATERTABLE_DURICRUST.md §3, step 5). Per-node, rank-local — no collective; localToGlobal on duriH at the end for the halo.

Fringe favourability Φ = exp(−((wt − d0)/w)²) — a Gaussian band on the water-table depth wt = z − h centred on the mean fringe depth d0 (half-width w); →1 at the fringe, →0 far above/below. With the opt-in discharge_gate (_dischargeWeight) Φ is additionally restricted to groundwater discharge zones (absolute-accumulation crust).
Formation dduriH/dt = k_form·Φ·Ψ·(1 − duriH/duriH_max) (self-limiting to duriH_max); Ψ from _weatheringSupply. When soil is tracked (cptSoil) the formation supply is additionally regolith-limited — capped by _regolithSupplyRate (chemical crust growth cannot outpace physical regolith production; §8 soil coupling).
Breakdown: the per-step surface incision (z_last − z > 0) strips the crust top (−k_break·incision), plus a slow disequilibrium decay away from the fringe (−k_decay·(1−Φ)·duriH). Clipped to [0, duriH_max] (a crust eroded through resets to 0).

Writes duriH (duriHL/duriHG), the induration duriF = duriH/duriH_max and the armor multiplier duriKarmor = 1 − armor_max·duriF (consumed by the Phase-4 erodibility hook).

GWMesh._regolithSupplyRate()[source]#: Rate (m/yr) at which weathering-produced regolith can feed the duricrust, the cap on chemical crust formation when soil is tracked (§8): the soil production rate prodSoil · rain (chemical∝physical weathering congruency). Returns +inf (no cap) when prodSoil is unavailable, so a soil-off run is unaffected. Rank-local.

GWMesh._recordInduration()[source]#

Sync the live per-node induration duriF with the per-layer stratigraphic archive stratDuri (DESIGN_WATERTABLE_DURICRUST.md §9). Two directions:

Exhumation (read-up) — a previously buried, indurated layer now at the surface (its overburden eroded through last step) re-arms the live crust: duriF = max(duriF, stratDuri[top]) at the top non-empty layer (found as in _surfaceComposition), refreshing duriKarmor. The stacked-duricrust / relief-inversion behaviour of cratonic laterites.
Formation (write-down) — the live crust has a real thickness duriH that spans a depth range, not just the surface layer, so the induration is written into every layer whose top lies within ``duriH`` below the surface: stratDuri[k] = max(stratDuri[k], duriF) for all k with depth_above(k) < duriH (depth_above = the summed thickness of the layers above k). A thick crust therefore indurates several thin layers; each is preserved when later buried (deposeStrat fresh layers start at 0) and re-arms the surface when re-exhumed — so an exhumed crust resists incision over its full thickness, not one layer’s worth.

With Level-B geochemistry on, the same write-down also stamps the crust’s dominant solute species into stratCrustType (and, with in-model provenance, its dominant source region into stratCrustSource) across the crust’s depth range, so a stratigraphic section preserves what each crust layer is and where its chemistry came from — the categorical companions to the stratDuri degree. These are archive-only (no exhumation read-up).

No-op (surface-only duriF, no archive) when stratDuri is unallocated (stratNb == 0). Composition-only — no geometry change. Rank-local (per-node); no collective.

GWMesh._surfaceSourceClass()[source]#

Per-node current surface lithology class — the dominant provenance class of the top non-empty stratigraphic layer (from stratP), falling back to the static bedrock source_class where a column has no sediment. As erosion exhumes deeper layers (or deposition buries the surface under transported sediment of a different provenance), the top layer — and hence this label — changes, so a surface_class-driven weatherability tracks the rock actually exposed each step.

Uses the same top-non-empty-layer scan as _recordInduration / _surfaceComposition. Needs provenance (stratP); returns the plain source_class (or None) when the per-layer class record is absent. Rank-local; no collective.

GWMesh._resolveGeoWeather()[source]#

Resolve the per-species weatherability to a list of scalar or ``(lpoints,)`` array entries (Level-B extension 1 — lithology → chemistry). Lets lithology control which species each region yields:

(a) per-vertex map — a species’ weatherability: [file, key] is loaded and subset to the local partition (self.locIDs), exactly like duriWeatherability / gwInfiltration.
(b) table by a standalone lithology map — a species’ weatherability_by_class: [...] is gathered by a per-vertex integer lithology: [file, key] map (independent of provenance).
(c) table by the provenance label — the same weatherability_by_class gathered by a per-vertex class label: weatherability_from: source_class uses the static bedrock class, while weatherability_from: surface_class uses the dynamic current surface class (the top stratigraphic layer, recomputed each step, so weatherability tracks exhumation/burial). Both need no new input beyond provenance.
otherwise the scalar (the default; byte-identical to the old path).

The lithology label for (b)/(c) is the lithology: map when given, else the provenance class (static source_class or dynamic surface_class). Resolved lazily (first _updateSolute) so it runs after provenance seeds source_class; the dynamic surface_class form is re-resolved every step. Rank-local, partition-exact, no collective.

GWMesh._makeSoluteKSP()[source]#

Cached KSP for the solute-transport solve: a direct LU factorisation (preonly + lu; MUMPS in parallel), gw_solute_ prefix.

The upwind advection operator is a non-symmetric M-matrix, but its diagonal is only weakly dominant: at a flat, poorly-drained saturated node the lateral outflow, the seepage sink and the recharge floor are all small, so the operator is badly conditioned. A Krylov + block-Jacobi solve then STALLS there for the strong (high-weatherability) tracers — it returns a non-converged iterate with a concentration ~10× too large, which silently breaks the per-step dissolved = precipitated + exported mass balance (the ocean/marine budget then drifts by an order of magnitude). A direct solve is exact, so it is both robust to the conditioning and conservative to round-off; the operator is the same size / sparsity as the flexure biharmonic that already uses this path, and the per-species factorisation is reused across the provenance right-hand sides. Env-overridable via the gw_solute_ prefix (e.g. request an iterative solver for meshes too large to factorise).

GWMesh._soluteAdvecCoeffs(h, T)[source]#: Upwind finite-volume coefficients for the steady solute advection ∇·(q c) by the groundwater flux q = −T∇h (area-normalised, per cell). Reuses jacobiancoeff — its off-diagonals are the area-normalised face conductances C_ik/A_i (correct for flat and global meshes alike), so the signed face flux i→k is f_ik = (C_ik/A_i)·(h_i − h_k). First-order upwinding puts the outflow (f>0) on the diagonal (carries c_i) and the inflow (f<0) on the neighbour column (carries c_k). Returns (adv, divq): adv is the (lpoints, 1+maxnb) array for _assembleDiffMatCSR (col 0 = diagonal), and divq = Σ_k f_ik is the lateral divergence (used for the seepage sink / export).

GWMesh._dischargeWeight(divq, R)[source]#

Absolute-accumulation gate (opt-in discharge_gate): a per-node weight G ∈ [0,1] that restricts crust formation to groundwater discharge zones. The lateral groundwater flux q = −T∇h converges into a discharge cell (div q < 0), delivering dissolved load that emerges at the surface and precipitates — the absolute-accumulation (lateral) style of valley/footslope ferricrete, as opposed to the in-situ / relative accumulation that indurates wherever the water table is shallow.

G = (−div q)⁺ / ((−div q)⁺ + R) — the fraction of the cell’s upward discharge that was imported laterally (vs supplied by local recharge R). Both terms are area-normalised rates (m/yr), so G is dimensionless and self-scaling (no new tuning constant): G→1 at a strongly convergent valley floor, G→0 at a recharge (divergent) rise. Returns all-ones when the gate is off (backwards-compatible). Rank-local.

GWMesh._solveSoluteTransport(source, dmask, dval)[source]#: Solve one steady tracer transport M c = source with Dirichlet nodes dmask pinned to dval (M = upwind advection of §``_soluteAdvecCoeffs`` at the current head). The seepage set is pinned (c leaves the aquifer there — the export sink), which anchors the M-matrix. G1 machinery: the physical dissolution source / fringe precipitation / baseflow export are added in G2/G3. Collective (KSP); returns the local concentration array.

GWMesh._soluteSolveRHS(M, rhs, guess)[source]#: Solve M c = rhs for the (pre-assembled) transport operator M with a warm-start guess, returning the local concentration (clipped ≥ 0). A thin wrapper so _updateSolute can reuse one assembled operator for several right-hand sides (the per-source-class provenance solves, G5). Collective (KSP).

GWMesh._updateSolute()[source]#

Level-B per-step solute update (G2), per tracer: dissolve → transport → precipitate → export, conservatively accounted.

Dissolution — a chemical-weathering source (the Level-A _weatheringSupply scaled per tracer by weatherability), on subaerial land only, debiting the conserved source pool.
Transport — the steady M c = D solve of _soluteAdvecCoeffs (upwind advection by q = −T∇h + the vertical seepage sink), plus a precipitation sink at the capillary fringe added to the diagonal (k_p·Φ where the tracer is super-saturated — a one-step-lagged c > c_sat gate).
Precipitation — the sink mass k_p·Φ·c feeds the duricrust duriH (Level-B thus replaces the Level-A proxy supply as the crust source; the induration/armoring then follow unchanged).
Export — by domain mass balance, the solute that is neither precipitated nor left in solution has discharged to the surface network (→ ocean; the flux is formalised in G3).

Budget per tracer (owned nodes): dissolved = precipitated + exported + Δ(in solution) — accumulated in gwDissolved/gwPrecip/gwOceanFlux for the conservation guard. Rank-local accounting (KSP solve collective).

GWMesh._routeRiverSolute()[source]#

Route the groundwater-exported (seepage/baseflow) solute down the surface drainage network to the shoreline — the river dissolved load (Level-B extension 2, DESIGN_GEOCHEM_EXTENSIONS.md §2), per species.

For each tracer an implicit accumulation solve reuses the cached flow-accumulation matrix fMati (the _getSedFlux pattern) with that species’ per-node seepage export gwSoluteFluxSp[:, k] (already an extensive m³/yr flux — no area weighting) as the source:

conservative (river_decay = 0): (I − Wᵀ) L = s — a passive tracer; solute reaching a coast/outlet leaves the continent, solute into a closed basin is trapped (evaporite), both from the filled-topo matrix with no special handling.
in-transit loss (river_decay = κ > 0): a first-order removal along the network adds a diagonal sink, (I − Wᵀ + κ I) L = s (built as fMati.shift(κ)); the lost mass κ·Σ L is a per-species diagnostic (in-channel precipitation / uptake).

Per-species fields riverSoluteSp / riverSoluteToOceanSp / riverSoluteLost are summed into the totals riverSolute / riverSoluteToOcean. With marine_coupling on, the delivered coastal flux accumulates into a per-species marine reservoir marineSolute (m³, integrated) and the per-node coastal input marineSoluteInput.

Called at the end of updateGroundwater — flowAccumulation (which caches fMati) runs earlier in the step, so the matrix is fresh and no reordering is needed. No-op unless river_load is on; the flow matrix must exist (guarded for the first call). Collective (KSP per tracer).