Water table and generic duricrust#

goSPL can evolve a water table and a generic duricrust — near-surface hydrology coupled to a chemical armoring of the erodibility. Each step a shallow groundwater head is solved on the mesh from the infiltrated rainfall; where the water-table depth lies in a capillary fringe an indurated crust precipitates and resists erosion, producing relief inversion and — with stratigraphy on — stacked, buried and re-exhumed duricrusts. The feature is enabled by a groundwater section in the input file and is a strict no-op when absent. The full design (algorithm, invariants, caveats) is in docs/DESIGN_WATERTABLE_DURICRUST.md.

The whole update runs once per step, after flow accumulation and before erosion, so each process reads the previous process’s state (soil → groundwater → duricrust → erosion). Every step of that loop is slow (\(10^4\)–\(10^6\) yr) and the coupling is explicit/sequential within a step, so there is no stiff intra-step feedback — the same stability argument as goSPL’s other explicit couplings.

Water-table head — Dupuit–Boussinesq#

The saturated groundwater flow under a free surface is described by the unconfined (Dupuit–Boussinesq) approximation. The water-table head \(h\) (elevation of the saturated surface) obeys

\[S\,\frac{\partial h}{\partial t} = \nabla\!\cdot\!\big(T(h)\,\nabla h\big) + R, \qquad T(h) = K_h\,\max(h - z_{bed},\, b_{min}),\]

with \(S\) the specific yield (drainable porosity), \(T\) the transmissivity, \(K_h\) the hydraulic conductivity, \(z_{bed}\) the impermeable base, \(b_{min}\) a minimum saturated thickness (so \(T>0\)), and \(R\) the net recharge. Recharge is the fraction of effective rainfall that infiltrates,

\[R = f_{infil}\,\max(0,\, \mathrm{rain} - \mathrm{evap}),\]

held at zero under standing water (sea, ponded lakes) and under ice, where the head is pinned to the surface or precipitation does not infiltrate.

The infiltration fraction \(f_{infil}\) is a scalar or a per-vertex map, and three optional (default-off) modulations refine it: subglacial recharge lets a fraction of the glacial meltwater (the ice model’s iceMeltRiverL) infiltrate under ice — the one recharge path allowed there; lithology scales it by the exposed coarse fraction so coarse infiltrates more than fine (dual lithology); and slope reduces it on steeper terrain as \(f/(1+\mathrm{slope}/\mathrm{slope}_{ref})\).

Implicit discretisation (and why)#

The groundwater response time \(\tau_{gw}\) is short relative to goSPL’s century-to-millennium steps \(\Delta t\), so an explicit update would be severely time-step limited. goSPL instead takes one implicit (backward-Euler) solve per step on the mesh’s finite-volume Laplacian \(L(T)\) (the same operator the marine diffusion uses):

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

As \(\Delta t/S\) grows this tends to the steady elliptic limit \(L(T)\,h = R\), so the solve is quasi-steady at long steps — exactly the regime goSPL runs in. The operator is a stiff 2-D elliptic system: it is solved with a Krylov method preconditioned by algebraic multigrid (fgmres + hypre BoomerAMG). A simple block-Jacobi/ILU preconditioner does not converge on it and silently drifts the water table to the surface, so multigrid is required — the same lesson as the flexure biharmonic.

Note

The implicit solve was validated against the exact steady 1-D Dupuit parabola \(s(x)^2 = s_d^2 + (R/K)(2Wx - x^2)\) on a drained strip (the benchmarks/test_dupuit.py fixture): the steady head matches the analytic profile to a root-mean-square error of 0.00 %.

Two non-linearities are handled with cheap inner iterations:

  • Unconfined transmissivity \(T(h)\): a few Picard sweeps (picard_its) lag \(T\) on the previous iterate.

  • Seepage free boundary \(h \le z\): the head is pinned to the surface (Dirichlet) at the base seepage set — sea, ponded lakes and open outlets — then, after each Picard block, any cell whose head over-topped the surface is added to that set and the solve repeated, up to seepage_passes passes. A dry-aquifer floor \(h \ge z_{bed}\) is imposed by clipping.

The aquifer base \(z_{bed}\) is prescribed (aquifer_base scalar or per-vertex map), or — with aquifer_base: from_soil — tied to the bedrock elevation \(z_{bed} = l_{Hbed} - d_{bedrock}\) (the permeable regolith over impermeable bedrock model of laterite terrains), which requires soil tracking. In a depositional basin the porous sediment fill is itself the aquifer, so the base is taken as the deeper of \(l_{Hbed} - d_{bedrock}\) and the bottom of the stratigraphic pile \(z - \sum H\) — uplands keep the thin-regolith base, basins deepen to the fill bottom.

Baseflow closure#

With conserve_baseflow on, the recharge that does not go into aquifer storage returns to the surface-water network as baseflow at the seepage nodes,

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

(written as the baseflow output). It is also re-injected into the surface-flow source: the runoff becomes \(b = \mathrm{rain}\,A - R\,A + Q_{seep}\), so the infiltrated recharge leaves surface runoff and re-emerges at the seepage nodes — rivers are physically baseflow-fed and total discharge stays \(\approx \mathrm{rain} - \mathrm{evap}\) (net-neutral globally; in the steady limit \(\sum \mathrm{baseflow} \approx \sum \mathrm{recharge}\)). It is applied once per step (the single runoff-source reset) so the two per-step flow-accumulation solves do not double-count it.

Lake ↔ aquifer volume coupling#

Lakes, rivers and the sea are fixed-head boundaries (\(h = z\)) for the head solve. With the opt-in lake_exchange the lake volume budget is additionally debited/credited by the across-bed groundwater flux. After the head solve the signed per-node exchange \(\nabla\!\cdot\!(T\nabla h)\,A = -(L h)\,A\) is formed (positive = the aquifer discharges into the lake, negative = the lake leaks into the aquifer), aggregated per lake, and fed into the lake’s inflow — gains added, leakage debited and clamped at the available water (mirroring the lake-surface evaporation). Off by default (the conventional fixed-head treatment).

Duricrust formation and armoring#

A generic capillary-fringe crust grows from the water-table depth \(wt = z - h\). Its thickness \(duriH\) evolves as

\[\frac{\mathrm{d}\,duriH}{\mathrm{d}t} = k_{form}\,\Phi\,\Psi\,\Big(1 - \frac{duriH}{duriH_{max}}\Big), \qquad \Phi = \exp\!\left[-\Big(\frac{wt - d_0}{w}\Big)^2\right],\]

self-limiting to \(duriH_{max}\). The fringe favourability \(\Phi\) is a Gaussian band on the water-table depth centred on the mean fringe depth \(d_0\) (half-width \(w\)): the crust forms only where the table sits at the fringe. The weathering supply \(\Psi\) has three modes: a default climate/temperature proxy \(\max(0, \mathrm{rain}-\mathrm{evap})^{p}\), an explicit rate (a Maher–Chamberlain kinetic–thermodynamic law driven by the recharge), or prodsoil (reuse of the soil production rate). When soil is tracked, formation is additionally regolith-limited — chemical crust growth cannot outpace physical regolith production. A breakdown term strips the crust under surface incision and relaxes it away from the fringe.

By default the crust indurates wherever the water table is shallow — the relative-accumulation (in-situ) style that forms the plateau / bowal cuirasses blanketing flat, low-relief laterite uplands. The optional discharge_gate switches to the absolute-accumulation (lateral) style of valley / footslope ferricrete: the favourability \(\Phi\) is multiplied by a groundwater-discharge weight

\[G = \frac{(-\nabla\!\cdot\mathbf{q})^{+}}{(-\nabla\!\cdot\mathbf{q})^{+} + R} \in [0,1],\]

the fraction of a cell’s upward discharge that is imported by lateral convergence of the groundwater flux \(\mathbf{q} = -T\nabla h\) rather than supplied by the local recharge \(R\) (both area-normalised rates, so \(G\) is dimensionless and needs no tuning constant). \(G \to 1\) at a convergent valley floor or seepage face and \(G \to 0\) at a divergent recharge rise, so the crust tracks the drainage network. The gate is off by default (\(G \equiv 1\), unchanged).

The induration degree \(duriF = duriH/duriH_{max} \in [0,1]\) is the single control on erodibility. It enters as a multiplicative armoring factor at the shared surface-erodibility hook,

\[K_{eff} = K \cdot \mathrm{surfLithoK} \cdot \big(1 - armor_{max}\,duriF\big),\]

so a fully indurated cell (\(duriF = 1\)) is \(armor_{max}\) less erodible (e.g. 0.9 ⇒ 10× more resistant). This modulates all three eroders (SPL, non-linear SPL and soil-aware SPL) with no branching in them. Armoring is rate-only: it changes no elevation the step it forms, so flow routing and mass balance are untouched — the crusted cells simply erode more slowly. This produces the classic dissected cuirasse: a laterite sheet is stripped from the incising valleys and survives on the armored divides as flat-topped capped mesas (differential-erosion relief growth). Note the armoring adds no geometry (no valley aggradation), so a valley-fill cap gets no elevation head-start; the standing landforms are the armored highs, not inverted valley floors.

Stratigraphic induration record#

When stratigraphy is recorded (a positive strat interval), the induration is archived per layer as stratDuri — an intensive property in \([0,1]\), distinct from the depositional erodibility stratK (they multiply at the exposed surface). Each step the live crust is written down across its depth range — every layer whose top lies within the crust thickness \(duriH\) below the surface is raised toward the live induration, so a thick crust indurates the several thin layers it physically spans (formation); a fresh deposit buries it with an uncemented (0) layer (preservation); and when erosion exhumes a previously buried indurated layer its preserved value re-arms the surface, which therefore resists incision over the crust’s full thickness (exhumation). This is the multi-cycle, stacked-duricrust / relief-inversion behaviour of cratonic (e.g. Australian laterite) landscapes. stratDuri advects with the pile like the porosity, is compaction-neutral, and is written to / restored from the stratal output; it is exposed per layer as the induration field of gospl-strata-volume.

Conservative geochemistry (Level B, opt-in)#

By default the duricrust supply is a local, non-conservative rate (§ above). Adding a groundwater: geochem: block turns on a conservative solute cycle: one or more lumped tracers are dissolved, transported with the groundwater flux, precipitated at the fringe (feeding the crust) and exported to the rivers, with a closed mass budget. See docs/DESIGN_WATERTABLE_GEOCHEM.md.

Each tracer c (an n_species array — a single tracer by default, several for calcrete / silcrete / ferricrete typing) obeys a steady advection–reaction balance, solved once per step (the transport equilibrates within a century step, just like the head):

\[\nabla\!\cdot\!(q\,c) = D - P, \qquad q = -T\nabla h,\]

  • Dissolution \(D\) — a chemical-weathering source on subaerial land (the climate/temperature supply scaled per tracer by weatherability), drawing from a conserved per-node source pool seeded as source_pool × cell area (a weatherable-rock reservoir, default 1e6). Dissolution debits the pool, so a pool too small for the run’s weathering rate exhausts and the tracer goes inert (a one-time warning is printed) — set it small to model an exhaustible weathering front or a wet-phase pulse, or large to keep weathering rate-limited over the whole run. The weatherability may vary in space (a per-vertex map, or a per-(class, species) table gathered by a standalone lithology map or the provenance class — the static source_class bedrock or the dynamic surface_class, the top stratigraphic layer that tracks exhumation/burial), so lithology sets which species each region yields — mafic rock → Fe/silica, a carbonate platform → carbonate.

  • Transport — first-order upwind advection by the groundwater flux q, built from the head operator’s face conductances (geometry-correct on flat and global meshes). A diagonal discharge-to-surface sink max(0, R − ∇·q) (the recharge that cannot be transmitted laterally leaves at the surface) anchors the operator as a well-posed M-matrix; the recharge term is essential where the groundwater carries no lateral flow — a fully saturated seepage node or a dry aquifer node — which would otherwise have a zero diagonal and spike the concentration. The steady transport is solved with a direct LU factorisation (MUMPS in parallel), which is exact and robust to the weak diagonal dominance at those poorly-drained nodes (an iterative solve can stall there and silently break the mass balance).

  • Precipitation \(P = k_p\,\Phi\,c\) — a linear sink at the capillary fringe Φ that feeds the crust duriH (replacing the local proxy supply when geochem is on). The saturation threshold c_sat is a future refinement.

  • Export — the seepage sink removes the solute discharging to the surface; because the upwind operator is exactly conservative, each step ``dissolved = precipitated + exported`` (machine precision). The dissolved export is the soluteflux output and a per-tracer flux to the ocean — the weathering-derived alkalinity/solute delivery relevant to long-term carbonate and climate budgets.

With several tracers, the dominant crust former per node is written as the crust_type field (different tracers win in different settings). With in-model provenance on, the solute is additionally transported per source-rock class (by linearity of the transport operator), so the downstream crust is attributed to the upgradient region where its solute dissolved — the dominant source is the crust_source field (solute-source provenance). The crust’s dominant species and dominant source region are additionally archived per stratigraphic layer (stratCrustType / stratCrustSource, the categorical companions to the stratDuri induration degree) — frozen on burial and exposed by gospl-strata-volume — so a cross-section shows what each buried crust is and where its chemistry came from. Outputs: solute (concentration), soluteflux (export), crust_type, crust_source. Coupled multi-species aqueous equilibrium (speciation, pH) is out of scope by design — the lumped multi-tracer model is the scale-appropriate choice at km / My.

With geochem: river_load on, the exported solute does not simply vanish to the ocean: it is routed down the surface drainage network, per species (the river dissolved load — the dominant natural pathway for weathering products to the sea). For each tracer an implicit accumulation solve reuses the flow-accumulation matrix with the per-node seepage export \(s\) as the source, so the dissolved load \(L\) (riverSolute) grows downstream and is delivered at the shoreline; solute routed into a closed continental basin is trapped there (evaporite behaviour). No deposition or cascade is needed — unlike sediment, dissolved load neither settles nor fills pits. An optional per-species river_decay adds a first-order in-transit loss (a diagonal term \((I - W^\mathsf{T} + \kappa I)\,L = s\) — in-channel precipitation/uptake), and marine_coupling accumulates the delivered coastal flux into a per-species marine reservoir. The coastal-delivered total cross-checks the lumped per-tracer ocean flux.

Compatibility#

The water table and duricrust compose with the other stratigraphic features: stratDuri is one more independent, intensive per-layer field alongside the dual-lithology fine pile (stratHf/phiF) and the provenance partition (stratP), and the armoring is a multiplicative factor at the same erodibility hook the others already use. A run with groundwater, duricrust, dual lithology and provenance all enabled conserves every invariant.