 |
OpenSWMM Engine
6.0.0-alpha.4
Data-oriented, plugin-extensible SWMM Engine (6.0.0-alpha.4)
|
|
OpenSWMM Engine
6.0.0-alpha.4
Data-oriented, plugin-extensible SWMM Engine (6.0.0-alpha.4)
|
The two-zone aquifer of Chapter 5 is integrated per subcatchment: one unsaturated moisture content and one water-table depth each, lateral outflow by the empirical power function of §5.3.5, no communication between neighbouring aquifers, and one feedback to the rest of the engine — the end-of-step infiltration cap. Four physically real processes are therefore not representable in it:
This chapter documents the opt-in kernel that supplies them: a vertically integrated two-zone subsurface on the cells of the 2D mesh of Hydraulics Chapter 9, following the semi-discrete finite-volume formulation of Qu and Duffy (2007). Each cell carries a saturated thickness and an unsaturated store joined at a moving water table; neighbouring saturated zones exchange by Darcy's law across the mesh's own edges; nodes on the mesh exchange with the aquifer beneath them by a conductance; and water a saturated column cannot hold returns to the surface cell above it, where the shallow-water solver routes it.
It is opt-in and changes nothing else. The kernel runs only when a [2D_AQUIFER] section is authored (or [2D_OPTIONS] GROUNDWATER YES says so); every deck without one is unchanged and the Chapter 5 aquifers keep their own path. The two cannot own the same water: a cell whose infiltration is routed to a subcatchment aquifer (Chapter 8 §8.4) cannot coexist with a [2D_AQUIFER], and a storage node with a bed on the mesh has its [STORAGE] exfiltration replaced by the node exchange of §9.6.
MODE PER_SUBCATCH. The design intends a second vehicle for the same kernel — one degenerate cell per subcatchment, no lateral flux, exchange with the outlet node. The code implements less: the token is parsed, written, and removes the lateral edge topology, but the kernel still runs on the mesh cells (the per-subcatchment cell list and outlet mapping are declared and never populated). Today it means "the mesh kernel without lateral flow" (§9.11).
Everything in the kernel is SI; the authored [2D_AQUIFER] rows stay in the project's units and are converted once, into the per-cell state (§9.12).
Figure 9-1 One aquifer column and its fluxes, with the unsaturated zone held as a bulk store (closure A) and as sigma layers (closure B)
Per cell \(i\) of planimetric area \(A_i\) and bed elevation \(z_{c,i}\):
| Symbol | Meaning | Source |
|---|---|---|
| \(z_s\) | Soil column thickness, bottom to surface (m) | ZS |
| \(z_{bed} = z_c - z_s\) | Aquifer bottom elevation (m) | derived |
| \(h_g \in [0, z_s]\) | Saturated thickness above the bottom (m); the saturated state | HG0 seeds it; default 0 |
| \(L = z_s - h_g\) | Unsaturated column thickness (m) | derived |
| \(H = z_{bed} + h_g\) | Water-table elevation, the head for lateral flow and node exchange (m) | derived |
| \(h_u\) | Closure A: bulk unsaturated storage as an equivalent depth of water (m) | seeded hydrostatic |
| \(\theta_j,\ j = 0..m-1\) | Closure B: water content of σ layer \(j\), layer 0 at the surface | seeded hydrostatic |
| \(K_s, \theta_s, \theta_r, \alpha, \psi_b, \lambda, n, L_M, c_{loss}\) | Soil parameters (§9.3, §9.5) | row columns and keywords |
Water in the cell — the continuity ledger's storage term (§9.10) — is
\[W_i = \bigl(\theta_s\, h_g + S_{col}\bigr)\, A_i, \qquad S_{col} = \begin{cases} h_u & \text{closures A and ENSLAVED} \\ \sum_j \theta_j\, L\, \Delta\sigma & \text{closure B} \end{cases} \]
(9-1)
Rows resolve onto cells by scope, * then TAG then CELL, and by authored order within a scope; a residual content at or above the porosity is clamped to half the porosity with a warning. The closure-B layer count \(m\) is global (M_LAYERS, default 8, range 2..128; a per-row value is parsed but unused).
Four laws are selectable, globally through [2D_AQUIFER_OPTIONS] SOIL_CHAR or per row. Each supplies an effective saturation \(S_e(\psi)\), a water content \(\theta(\psi) = \theta_r + S_e(\theta_s - \theta_r)\) and a relative conductivity \(K_r(\psi)\), with suction \(\psi \ge 0\) measured upward from the water table ( \(\psi = L\) at the surface) and \(K(\psi) = K_s K_r(\psi)\):
| Law | \(S_e(\psi)\) | \(K_r(\psi)\) | Beyond \(K_s, \theta_s, \theta_r\) |
|---|---|---|---|
GARDNER (1958) | \(e^{-\alpha\psi}\) | \(e^{-\alpha\psi}\) | \(\alpha\) |
RUSSO (1988), default | \(\bigl[(1 + \tfrac12\alpha\psi)\, e^{-\alpha\psi/2}\bigr]^{2/(2+L_M)}\) | \(\bigl[(1 + \tfrac12\alpha\psi)\, e^{-\alpha\psi/2}\bigr]^{2}\) | \(\alpha\), Mualem \(L_M\) (0.5) |
BROOKS_COREY (1964) | 1 for \(\psi \le \psi_b\), else \((\psi_b/\psi)^{\lambda}\) | 1 for \(\psi \le \psi_b\), else \((\psi_b/\psi)^{2 + 3\lambda}\) | \(\psi_b, \lambda\) |
VAN_GENUCHTEN (1980) | \(\bigl[1 + (\alpha\psi)^{n}\bigr]^{-m},\ m = 1 - 1/n\) | \(S_e^{L_M}\bigl[1 - (1 - S_e^{1/m})^{m}\bigr]^{2}\) | \(\alpha, n, L_M\) |
Russo's retention is the Mualem-consistent inverse of his conductivity, so it costs the same as Gardner and is the default; Brooks–Corey's air-entry branch is the one non-smooth law; van Genuchten opens the Rosetta and UNSODA databases. \(S_e\) is kept in \([10^{-8}, 1]\).
Figure 9-2 Retention and relative conductivity of the four soil-characteristic laws at their default parameters
Hydrostatic-equilibrium storage. The water a column of thickness \(L\) holds when hydrostatic above the table is
\[ h_u^{\ast}(L) = \int_0^{L} \theta(\psi)\, d\psi = \theta_r L + (\theta_s - \theta_r)\, \sigma(L), \qquad \sigma(L) = \int_0^{L} S_e(\psi)\, d\psi \]
(9-2)
with \(\sigma(L)\) closed-form for Gardner, \((1 - e^{-\alpha L})/\alpha\), and for Brooks–Corey, \(\psi_b + \psi_b\bigl[(L/\psi_b)^{1-\lambda} - 1\bigr]/(1 - \lambda)\) for \(L > \psi_b\) (else \(L\); \(\lambda = 1\) takes the logarithm), and by four-point Gauss–Legendre quadrature for Russo and van Genuchten; \(L\) is floored at 1 mm. This integral is the ENSLAVED closure's whole unsaturated state and the initial condition of every column.
The closed-form recharge. For Gardner, the quasi-steady recharge across the table is Qu and Duffy's equation 22 verbatim, on the saturation scale \(\tilde h_u = (h_u - \theta_r L)/(\theta_s - \theta_r)\):
\[ q_0(h_u, h_g) = K_s\, \frac{1 - e^{-\alpha L} - \alpha\, \tilde h_u}{\alpha L\,\bigl(1 - e^{-\alpha L}\bigr)}, \qquad L = z_s - h_g \]
(9-3)
positive downward, negative for capillary rise, smooth in both states. It is identical to a relaxation of the stored water toward its hydrostatic equilibrium at a conductivity-set rate,
\[ q_0 = C(L)\,\bigl[\tilde\sigma^{\ast}(L) - \tilde h_u\bigr], \qquad \tilde\sigma^{\ast}(L) = \frac{1 - e^{-\alpha L}}{\alpha}, \quad C(L) = \frac{K_s}{L\,\bigl(1 - e^{-\alpha L}\bigr)} \]
(9-4)
The other three laws generalise (9-4): the equilibrium is the law's own \(\sigma(L)\) from (9-2) and the rate is the conductivity at the column's mean suction over that integral,
\[ q_0 = \frac{\bar K}{\sigma(L)}\,\bigl[\sigma(L) - \tilde h_u\bigr], \qquad \bar K = K\!\left(\psi = \tfrac12 L\right) \]
(9-5)
which for Gardner gives the same order and limits as (9-4).
The Experimental caveat. The code states, and this manual repeats: (9-5) is the engine's generalisation, not the literature's. Gardner reproduces the published closed form exactly and is verified against it. Russo, Brooks–Corey and van Genuchten reproduce the correct equilibrium and sign, but their relaxation rate is a modelling choice rather than a published result. The closure-ladder benchmark and HYDRUS-1D comparison the plan schedules would license them; until those run the three carry Experimental in §9.11, Gardner or Russo are preferred for anything quantitative, and closure B — which integrates the real Richards flux and needs none of this — is preferred where the answer matters. The caveat concerns the recharge closure only.
The dimensionless group read by the closure selection is \(\alpha L\), the column thickness in capillary lengths; for Brooks–Corey the law's own inverse length replaces \(\alpha\) ( \(\alpha L \equiv L/\psi_b\)).
The wilting stress reducing evaporative demand is a smooth Feddes-style multiplier of the suction,
\[ S(\psi) = \begin{cases} 1 & \psi \le 0 \\ s^{2}(3 - 2s),\ s = 1 - \psi/\psi_w & 0 < \psi < \psi_w \\ 0 & \psi \ge \psi_w \end{cases} \]
(9-6)
deliberately not the legacy if (infil <= 0) gate and not piecewise-linear: a kink here shows up as a kink in the ET ledger that reads like a solver artefact.
The unsaturated store is represented per cell by CLOSURE (in [2D_AQUIFER_OPTIONS] or on a row), resolved from AUTO at initialisation.
AUTO reads \(\alpha L\) of the cell's initial column:
| \(\alpha L\) | Closure | Reasoning |
|---|---|---|
| \(< 1\) | ENSLAVED | The column equilibrates far faster than the table moves; its storage is a function of the table |
\(1 \le \alpha L < 5\), or FORCE_CLOSED_FORM YES | CLOSED_FORM | The quasi-steady recharge (9-3)/(9-5) holds |
| \(\ge 5\) | SIGMA | The quasi-steady assumption fails: travel time, profile memory and a localised capillary fringe are the physics, and only a resolved column carries them |
A warning reports how many cells resolved to each. The design record also describes a soft warning near \(\alpha L \approx 3\) and a refusal of ENSLAVED above \(\alpha L = 2\); neither is in the code, which accepts CLOSURE ENSLAVED on any column.
CLOSED_FORM (closure A)One bulk state \(h_u\) driven by the quasi-steady recharge. At a firing of step \(\Delta t\), with the table moving from \(L_0\) to \(L_1\) (§9.5),
\[ h_u^{\,1} = h_u^{\,0} + \bigl(q^{+} - q_{ET} - q_0\bigr)\,\Delta t + \theta_{bot}\,(L_1 - L_0) \]
(9-7)
where \(q^{+}\) is the infiltration delivered from the surface (Chapter 8), \(q_{ET}\) the demand of §9.8, \(q_0\) from (9-3) or (9-5), and the last term the handover — the water the moving table carries across the zone boundary. The slab content is the bulk store's own mean, \(\theta_{bot} = \mathrm{clamp}(h_u^{\,0}/L_0,\ \theta_r,\ \theta_s)\), and it is the same number the saturated update pairs with its storage coefficient \(S_y = \theta_s - \theta_{bot}\); if the two zones disagreed on how much water a moving table carries, the cell would leak at a rate proportional to \(\lvert \dot h_g \rvert\). Drained, it reduces to \(\theta_s - \theta_r\).
Two caps keep the books. Above \(\theta_s L_1\) the excess is rejected to the surface (§9.7). Below \(\theta_r L_1\) — a table rising faster than the fluxes refill the column — the shortfall is drawn back up out of the saturated zone in the same firing, lowering the table by \(\delta = (\theta_r L_1 - h_u^{\,1})/(\theta_s - \theta_r)\) (bounded by its thickness); a shortfall no water anywhere can meet is booked as ET that was never available rather than invented. The recharge is also capped by availability: downward, \(q_0 \le (h_u - \theta_r L_0)/\Delta t\); upward, \(-q_0 \le h_g(\theta_s - \theta_r)/\Delta t\). Without the first cap a cell that cannot reach equilibrium drains past residual and keeps "delivering" at \(\bar K\) forever.
ENSLAVEDThe shallow-water-table reduction of Qu and Duffy's equation 39: \(h_u\) is algebraic in \(h_g\), \(h_u = h_u^{\ast}(L)\) from (9-2), and the cell has one state. Since \(d h_u^{\ast}/dL = \theta(L)\), the column contributes its surface water content as a storage debit, \(S_y = \theta_s - \theta(L)\), and infiltration drives the table with no lag: \(q_0 = q^{+} - q_{ET}\).
The table is not updated by (9-10) alone. The balance \(W(h) = \theta_s h + h_u^{\ast}(z_s - h)\) is monotone, so the new table solves \(W(h_g^{\,1}) = W(h_g^{\,0}) + \Delta V_{sat}/A\) on \([0, z_s]\) by a safeguarded Newton with bisection fallback on the volume residual ( \(W\) is flat across the capillary fringe, where only bisection works). What the interval cannot hold is booked exactly: above \(z_s\) as Dunne excess, below 0 as a refund to the sinks (§9.5). The unguarded Newton this replaced created 23.5 m³ on a two-cell exfiltration deck.
SIGMA (closure B)An explicit column of \(m\) layers on \(\sigma = z/L \in [0, 1]\), \(z\) downward from the ground, \(\Delta\sigma = 1/m\). Layers stretch and compress as the table moves with no regridding and a fixed SoA stride, so a layer sweep vectorises across columns. The conserved quantity is the water depth per layer \(w_j = \theta_j L \Delta\sigma\); \(\theta\) is re-derived from \(w\) and the new \(L\) at the end of the sweep, which makes the compression exact. With \(\dot L = (L_1 - L_0)/\Delta t\) and face \(j + \tfrac12\) at \(\sigma_f = (j+1)/m\), the downward-positive flux is
\[ F_{j+1/2} = K\!\bigl(\psi(\theta_j)\bigr) \;-\; \theta_{donor}\,\sigma_f\,\dot L \;-\; \frac{\bar D_{j+1/2}}{L_0\,\Delta\sigma}\,\bigl(\theta_{j+1} - \theta_j\bigr) \]
(9-8)
Gravity is downward, so the advective donor is the layer above (first-order upwind; MUSCL is a recorded follow-up). The grid term is upwinded on its own sign: with \(\dot L > 0\) (table falling) the donor is the layer below. The capillary term, \(D = K\, d\psi/d\theta\) harmonic-mean across the face, is included only under CAPILLARY_DIFF YES (default off). The boundaries:
\[ F_{top} = \min\!\bigl(q^{+},\ (\theta_s - \theta_0)\, L_0 \Delta\sigma / \Delta t\bigr) - q_{ET}^{taken}, \qquad F_{bot} = q_{0,phys} - \theta_{m-1}\,\dot L \]
(9-9)
At the top, the first layer's headroom caps the infiltration accepted and the remainder is rejected to the surface (§9.7); the ET taken is the demand reduced by (9-6) at the top layer's own suction (wilting suction at 1 % effective saturation) and limited to what the layer holds above residual. At the bottom, \(q_{0,phys}\) is the bottom layer's gravity drainage \(K(\psi(\theta_{m-1}))\) and the second term is the handover swept across the moving boundary. Both together are what the saturated zone gains, which is why (9-10) books \(q_{0,phys}\) with \(S_y = \theta_s - \theta_{m-1}\) rather than \(F_{bot}\) with \(\theta_s\): the same equation, but only the first conserves once the column retains water above a falling table. The design record writes \(\theta_s\); \(\theta_s - \theta_{bot}\) is the textbook specific yield and what the code uses.
Positivity is a per-layer availability share: a layer may not export more than it holds above residual in one step. After the update \(\theta\) is re-derived from the new geometry. A rising table compresses every layer, and a layer whose water then exceeds \(\theta_s L_1 \Delta\sigma\) is saturated: the surplus cascades downward and whatever reaches the bottom is reported as an overflow to the saturated zone applied to \(h_g\) in the same firing; a deficit below residual is drawn upward the same way. Clamping the surplus away inside the sweep is the classic σ-grid failure and lost 0.39 m of water on the closed-column oscillation gate. Columns are seeded hydrostatic above their initial table, layer centres at \(\psi = L(1 - \sigma)\).
At a cell's firing every saturated source shares one storage coefficient:
\[ \Delta V_{sat} = \bigl(q_0 - q_{deep}\bigr)\, A\, \Delta t \;+\; V_{lat} \;-\; V_{node} \;+\; V_{link}, \qquad h_g^{\,1} = h_g^{\,0} + \frac{\Delta V_{sat}}{S_y\, A} \]
(9-10)
where \(V_{lat}\) is the lateral volume gathered from the cell's edge accumulators since its last firing (positive in), \(V_{node}\) the node exchange gathered from its bed accumulators (positive out, §9.6), \(V_{link}\) the conduit seepage delivered (§9.6), and \(S_y\) the closure's specific yield of §9.4, floored at \(10^{-3}\). The split is: recharge from the current column state, then (9-10), then the clamp, then the column sweep at the clamped \(\dot L\), then the column's overflow or deficit applied to \(h_g\), then Dunne. Deriving \(\dot L\) from the clamped table is what makes the clamp harmless: no water is created at \(h_g = z_s\) or destroyed at \(h_g = 0\).
The clamp keeps its books. A table above \(z_s\) becomes Dunne volume \((h_g^{\,1} - z_s)\, S_y A\) (§9.7). A negative table means the sinks asked for more than the aquifer holds; the shortfall \(-h_g^{\,1} S_y A\) is refunded to the sinks that overdrew — deep loss first, then the node exchange — and the refunded amounts are what the ledger records.
Per unique interior edge \(e\) between cells \(l\) and \(r\), of length \(\xi_e\) and centroid distance \(d_e\),
\[ Q_e = \bar K_e\, T_e\, \xi_e\, \frac{H_l - H_r}{d_e}, \qquad \bar K_e = \frac{2 K_{s,l} K_{s,r}}{K_{s,l} + K_{s,r}}, \qquad T_e = h_{g,\mathrm{donor}} \]
(9-11)
positive from \(l\) to \(r\). The conductivity is the harmonic mean (zero if either side is) and the transmissivity is the donor cell's saturated thickness, upwind on the head difference — the mean here is the classic way to make a drying aquifer pump water it does not have. Positivity is a per-face share of the donor's drainable water, split across the donor's faces so a cell with many wet neighbours cannot be emptied several times over:
\[ \lvert Q_e \rvert \le \frac{\tfrac12\, T_e\, S_{y,\mathrm{donor}}\, A_{\mathrm{donor}}}{n_{f,\mathrm{donor}}\, \Delta t} \]
(9-12)
The volume \(Q_e \Delta t\) is booked \(\mp\) into the edge's two side accumulators at the face's cadence and gathered by each cell at its own firing — the surface marcher's strategy, which makes conservation across tiers inherited rather than re-proved (§9.9). Under MODE PER_SUBCATCH the edge topology is empty and (9-11) is never evaluated.
\[ q_{deep} = c_{loss}\, \frac{h_g}{z_s} \]
(9-13)
a linear sink at full-saturation rate \(c_{loss}\) (C_LOSS, default 0) — exactly Chapter 5's deep percolation (5-21), with \(c_{loss}\) as \(DP\) and \(h_g/z_s\) as \(d_L/(E_G - E_B)\). The design record also allows a [GWF] DEEP expression; the code implements the linear form only.
A 1D node on the mesh exchanges with the aquifer of the cell containing it through a MODFLOW-River conductance, sign-determined by heads. With the pipe head \(H_{pipe} = z_{inv} + d\) (converted to metres) and the aquifer head \(H_{gw} = z_{bed} + h_g\),
\[Q_{node} = C\,\bigl(H_{gw} - H_{pipe}\bigr), \qquad C = \begin{cases} \dfrac{K_c\, A_b}{d_C} & \text{semi-confining bed given (config b)} \\[1.2em] \dfrac{K_s\, A_b}{\tfrac12 z_s} & \text{no bed: direct Darcy over half the column (config a)} \end{cases} \]
(9-14)
positive out of the aquifer into the pipe. \(K_c\) and \(d_C\) are the bed's conductivity and thickness (KC, DC); \(A_b\) is AREA, or the cell's area when omitted. Table above pipe head is groundwater inflow, pipe head above table is exfiltration: one equation, no fitted coefficients.
Figure 9-3 Node-aquifer exchange through a semi-confining bed, in both directions, with the caps that bound it
Caps. A drain ( \(Q > 0\)) may take at most half the cell's drainable water per step, \(\tfrac12 h_g S_y A / \Delta t\). A recharge ( \(Q < 0\)) is guarded three ways: a column whose table is within \(10^{-4} z_s\) of the ground takes nothing (it would come straight back as Dunne excess and ping-pong against the frozen node head); the saturated zone's headroom \((z_s - h_g)\, S_y A\) is taken in the same half-per-step share; and the node cannot give more than it holds for the batch plus what flows through it — a junction stores nothing below its rim but can lose water at the rate it is fed.
Cadence and booking. The exchange is sampled at the marcher's tier-0 cadence against the batch-frozen 1D heads, like the surface's junction exchange, and booked into the bed's accumulator; the aquifer cell gathers its share when it fires, and neither side is pinned to the other's tier. Per node, the batch's volume is added to the node's coupling volume for the next routing step's lateral-inflow assembly — without entering the surface ledger; its home is the aquifer's led_node.
Which nodes exchange. NODE_ENROLMENT AUTO (default) enrols every node whose [COORDINATES] fall in a mesh cell with a direct-Darcy bed over the cell's area unless a [2D_AQUIFER_NODE] row names it; a row overrides with its own bed (cell given explicitly, 1-based, or as CELL AUTO) or opts out with EXCHANGE NO; ROWS limits exchange to the authored rows. One bed per node is enforced; a bed off the mesh is a warning and the node exchanges with nothing. A storage node with a bed has its [STORAGE] exfiltration switched off, since both would lose the same water twice.
Conduit seepage. With LINK_SEEPAGE AUTO (default), every conduit with a [LOSSES] seepage rate whose polyline crosses the mesh delivers its seepage volume into the cells it crosses, weighted by the length inside each (exact Cyrus–Beck clipping). The 1D side keeps booking the loss; the volume lands in led_link and rides (9-10) as \(V_{link}\). NONE keeps it a system loss.
Two physically distinct routes carry water back to the surface cell above, and both are booked to the same per-cell accumulator:
ENSLAVED, the volume above \(W(z_s)\)) and the table is held at \(z_s\).Rejected infiltration — water the top of the column would not take: under closure A the excess of (9-7) above \(\theta_s L_1\), under closure B \(q^{+} - F_{top}\) from (9-9). ENSLAVED has no column to reject from.
\[ V_{\to surf} = \bigl[\,\mathrm{DUNNE}\,\bigr]\, V_{Dunne} + q_{rej}\, A\, \Delta t \]
(9-15)
Together they are the mass balance, which is why DUNNE defaults on; the flag exists for compatibility comparisons and disables only the first route. The volume waits in the cell's surface-bound accumulator until the surface cell fires and drains it as a source in its own budget (Chapter 8, equation 8-6). A dry surface cell above a saturating column would otherwise hold that water forever, so the marcher pins every cell with pending surface-bound water active at each rebuild — the one seed of the active set that does not come from the surface's own state, and what lets return flow re-emerge in a valley bottom with no rain falling. Rainfall on such a cell simply lands on the surface, where the column accepts none of it and the surface routes it.
GW_ET — in [2D_AQUIFER_OPTIONS], or as a [2D_OPTIONS] alias folded into the same option at open — selects how the subsurface meets evaporative demand:
| Mode | Capillary rise (negative \(q_0\)) | Demand at the column top |
|---|---|---|
NONE (default) | Suppressed: a negative \(q_0\) is zeroed | None |
CAPILLARY_RISE | Allowed | None |
BOUNDARY_ET | Suppressed | Yes |
BOTH | Allowed | Yes |
The demand is the surface cell's evaporation rate of Chapter 8 §8.5 — zero under the default [2D_OPTIONS] EVAPORATION YES with no forcing, the project's [EVAPORATION] under CLIMATE — reduced by (9-6) at the column's mean suction with the wilting suction at 150 m (about \(-15\) bar):
\[ q_{ET} = e_i\; S\bigl(\psi(\bar S_e)\bigr), \qquad \bar S_e = \mathrm{clamp}\!\left(\frac{h_u/L_0 - \theta_r}{\theta_s - \theta_r},\ 10^{-6},\ 1\right), \quad \psi_w = 150\ \mathrm{m} \]
(9-16)
with \(\psi(\bar S_e)\) the law's retention inverse. Under closure A and ENSLAVED this \(q_{ET}\) is removed from the column (or, enslaved, from the table) and booked. Under closure B it is offered to the column sweep, which applies its own stress at the top layer's suction with the wilting suction at 1 % effective saturation and limits the take to what the top layer holds; the amount actually removed is booked. The two stresses therefore compound under SIGMA — an implementation detail of the column rather than a documented design choice.
This replaces Chapter 5's three-parameter split — upper-zone fraction, lower-zone depth and the if (infil <= 0) gate — with one boundary condition at the top of the unsaturated zone and one physical path from below: there is no separate lower-zone ET term, water reaches the top by capillary rise (the negative branch of (9-3)/(9-5), or the σ column's upward propagation), and the closure that governs recharge downward governs capillary rise upward. The wilting point survives; the lower-zone depth does not. Capillary rise is capped by what the table holds (§9.4.2) and counted as the negative share of recharge.
The subsurface has no loop of its own. It joins the surface marcher's halving-order tier ladder through four hooks — tier assignment at every rebuild, a face phase and a cell phase per due tier, and a settle before any re-tier — under two guarantees enforced by construction. G-A, independence: a cell's tier comes only from its own bound \(\min(\Delta t_g, \Delta t_u)\); the surface's pin-to-tier-0 rule for coupling cells is not applied, so a cell with node exchange, infiltration or Dunne transfer keeps its tier and the cross-domain volume accumulates until it fires. G-B, conservation: every flux channel books \(\pm\Delta V\) into a side accumulator at the producer's cadence and is gathered by the owner at its own firing — lateral Darcy on the edges, node exchange on the beds, infiltration and seepage per cell, Dunne on the surface-bound accumulator.
The saturated bound is that of explicit diffusion on the head:
\[ \Delta t_g \le C_{gw}\, \frac{A\, S_y}{\displaystyle \sum_e \bar K_e\, T_e\, \xi_e / d_e \;+\; c_{loss} A / z_s \;+\; \sum_{beds} C_b}, \qquad T_e = \max\bigl(h_{g,i},\ h_{g,j},\ 0.01\, z_s\bigr) \]
(9-17)
with \(S_y = \theta_s - \theta_r\) (floored), the wetter of the pair with a floor for the transmissivity so an empty aquifer does not report an infinite step and then take one, and the deep-loss and node conductances folded in as linear sinks. C_GW defaults to 0.5, range (0, 1]. For realistic conductivities the bound is minutes to hours — orders above the surface tiers, which is what local time stepping exploits.
The unsaturated bound depends on the closure. A σ column's own bound is
\[ \Delta t_u \le C_{col}\, \min_j \frac{L\, \Delta\sigma}{c(\theta_j) + \lvert \sigma_j \dot L \rvert + 2\bar D_j / (L\Delta\sigma)}, \qquad c(\theta) = \frac{dK}{d\theta} \]
(9-18)
with the celerity by a centred difference through the retention inverse and \(\dot L\) estimated from the previous firing, \(-q_{0,last}/S_y\). Closure A is a relaxation whose stiffness \(C = -\partial q_0/\partial h_u\) is taken numerically, so \(\Delta t_u \le C_{col}/C\) covers every soil law; ENSLAVED has no unsaturated bound. C_COL defaults to 0.9.
Tiers and ordering. A cell's rung is \(\lfloor \log_2(\Delta t_{cell}/dt_0) \rfloor\) by the same ilogb arithmetic the surface uses, clamped to the ladder. The base step is the minimum of the surface's and the aquifer's finest bound, so a groundwater cell finer than every wet surface cell (or a whole domain with a dry surface) does not fire beyond its bound on rung 0. The ladder height is the surface's LTS_TIERS (default 4, at most 8); a slower cell lands on the coarsest rung and fires more often than it needs — the design record's runtime tier count (D-N1) is deliberately not taken, and requiredTiers() is unused. A face fires at the finer of its two cells' rungs. Within a base substep the order is surface faces, subsurface faces, surface cells, subsurface cells (so a cell that infiltrated hands the water down within the substep), then node sampling at tier 0. A routing window that does not fit a whole macro cycle ends with global tail steps in which every aquifer rung fires once.
Settling does not apply pending volumes: turning a volume into a table rise also moves the unsaturated store, so all volume goes through the cell firing. Nothing strands, since accumulators are gathered through incidences a re-tier does not change; continuity audits between firings must add them (§9.10).
The kernel keeps cumulative domain totals in m³:
| Term | Meaning | Sign |
|---|---|---|
infil_in | \(q^{+}\) delivered from the surface | in |
link | Conduit seepage delivered | in |
lateral | Net lateral Darcy across the domain edge | net in |
deep | Deep percolation | out |
node | Node exchange | out of the aquifer |
et | Subsurface ET | out |
dunne | Saturation excess and rejection returned to the surface | out |
recharge, caprise | Unsaturated to saturated transfer and its negative share | internal |
init_storage, storage | Water at the start; water now | — |
Three storage views exist. storage() is (9-1) summed over cells. The live storage adds every side accumulator (lateral edge pairs, water in flight to or from the surface, seepage in flight) and is what the results file and the C API report. The ledgered storage adds only the lateral accumulators, since the surface accumulators are always one firing out of phase with infil_in and dunne. The continuity residual is
\[ \mathcal{R} = S_{ledgered} - S_0 - \bigl[(F_{in} + Q_{lat} + F_{link}) - (F_{deep} + F_{node} + F_{ET} + F_{Dunne})\bigr] \]
(9-19)
zero to machine precision for a conserving kernel; recharge and capillary rise are internal to a cell and do not appear.
When [REPORT] CONTINUITY is on and the kernel ran, a 2D Aquifer Continuity block follows the 2D Surface Routing Continuity block, in cubic metres: Initial Stored Volume, Infiltration Inflow, Conduit Seepage Inflow, Lateral Net Inflow, Deep Percolation, Node Exchange Outflow, Subsurface ET, Saturation Excess Return, Final Stored Volume (the live storage), then Recharge and Capillary Rise listed as internal. Its continuity error is (9-19) as a percentage of \(S_0 + F_{in} + F_{link} + \max(0, Q_{lat})\). By design the surface block's Infiltration Loss (Chapter 8 §8.6) is this block's Infiltration Inflow, and this block's Saturation Excess Return is the surface's aquifer inflow term.
With the GROUNDWATER group of REPORT_2D_VARIABLES (in DEFAULT), the results file carries per-cell fields held from each cell's last firing:
| Dataset | Content | Units |
|---|---|---|
Mesh2_face_gw_table_elev | Water-table elevation \(z_{bed} + h_g\) | m |
Mesh2_face_gw_hg | Saturated thickness | m |
Mesh2_face_gw_hu | Closure A unsaturated storage (0 under closure B) | m |
Mesh2_face_gw_recharge | \(q_0\), positive down, negative capillary rise | m/s |
Mesh2_face_gw_lateral | Net lateral Darcy into the cell | m³/s |
Mesh2_face_gw_node_exchange | Exchange with the node bed, positive into the pipe | m³/s |
Mesh2_face_gw_deep | Deep percolation | m/s |
Mesh2_face_gw_et | Subsurface ET | m/s |
Mesh2_face_gw_dunne | Saturation excess returned to the surface | m³/s |
Mesh2_face_gw_infil_in | Infiltration delivered from the surface | m/s |
Mesh2_face_gw_link_seepage | Conduit seepage delivered | m³/s |
plus two static descriptors, Mesh2_face_gw_bed_elev and Mesh2_face_gw_closure (0 closed form, 1 enslaved, 2 sigma); the domain series groundwater_ledger of shape [time, 12] in the order recharge, lateral, deep, node, dunne, caprise, et, infil_in, init_storage, storage, link, continuity_residual (read its terms attribute — the seepage term was appended and moved the residual); and groundwater_node_exchange_cum, the cumulative exchange per bed, positive out of the aquifer, with node names as an attribute. The GW_DETAILED group adds Mesh2_face_gw_theta_sigma of shape [time, layer, face], layer 0 at the ground, zero for a non-sigma cell; it is \(m \times\) the cell count per step and off by default.
| Family | Alternative | Behaviour | Status |
|---|---|---|---|
| Soil characteristic | GARDNER | Exponential retention and conductivity; recharge by equation 22 verbatim, verified against it | Implemented |
| Soil characteristic | RUSSO (default) | Mualem-consistent retention at Gardner's cost; recharge by (9-5), whose rate is the engine's | Experimental |
| Soil characteristic | BROOKS_COREY | Air-entry branch, closed-form equilibrium storage; recharge by (9-5) | Experimental |
| Soil characteristic | VAN_GENUCHTEN | Rosetta/UNSODA parameters, quadrature equilibrium; recharge by (9-5) | Experimental |
| Closure | AUTO | Per-cell selection from \(\alpha L\) at 1 and 5 | Implemented |
| Closure | CLOSED_FORM | Bulk \(h_u\) with quasi-steady recharge and the handover term | Implemented |
| Closure | ENSLAVED | \(h_u\) algebraic in \(h_g\); bracketed volume solve | Implemented |
| Closure | SIGMA | Explicit ALE σ column, first-order upwind | Implemented |
| Closure | SIGMA with MUSCL reconstruction | Recorded follow-up in the column's contract | Planned |
| Column option | CAPILLARY_DIFF YES | Harmonic-mean diffusive term in (9-8) and its bound | Implemented |
| Mode | MODE MESH (default) | Kernel on the mesh cells with lateral Darcy | Implemented |
| Mode | MODE PER_SUBCATCH | Today: the mesh kernel with lateral topology removed; the one-cell-per-subcatchment vehicle is not in the code | Experimental |
| Deep loss | Linear \(c_{loss} h_g / z_s\) | Chapter 5's (5-21) form | Implemented |
| Deep loss | [GWF] DEEP expression | Named in the design record only | Planned |
| Node exchange | NODE_ENROLMENT AUTO | Every node inside the mesh, direct Darcy over the cell area | Implemented |
| Node exchange | NODE_ENROLMENT ROWS | Authored [2D_AQUIFER_NODE] rows only | Implemented |
| Node exchange | Semi-confining bed (KC, DC) | MODFLOW-River config b conductance | Implemented |
| Conduit seepage | LINK_SEEPAGE AUTO | Length-weighted delivery of [LOSSES] seepage into crossed cells | Implemented |
| Conduit seepage | LINK_SEEPAGE NONE | Seepage stays a system loss | Implemented |
| Return flow | DUNNE YES (default) | Table-at-ground excess to the surface | Implemented |
| Return flow | DUNNE NO | Compatibility only; rejection still returns | Implemented |
| ET | GW_ET NONE, CAPILLARY_RISE, BOUNDARY_ET, BOTH | §9.8 | Implemented |
| Time stepping | Standalone CLOSED_FORM/ENSLAVED/SIGMA transport/age/water-temperature FULL_SWE SSP-RK2 | Global tier, conserved equilibrium-water blends for ENSLAVED and layer-water blends with starting-state ALE velocity for SIGMA, knot-aligned wells and HEAD/FLOW quality edges, canonical bulk/layer species and age ledgers, symmetric groundwater/surface reactions, rollback and exact component checkpoint; a porous profile supports stored soil energy, midpoint/lateral conduction, chemical/age inventories and phase-aware reactions with sensible/latent soil ET and constant, project-series or held AIR thermal boundaries; a fully covered COUPLED/GROUND_COLUMN finite-contact profile also supports RK2 heat, chemical/age contact and reactions with nonlinear radiation, sensible atmosphere and accepted surface latent cooling; copied surface/groundwater sources, runtime aquifer HEAD/FLOW edges and rain/quality/ET/paired infiltration and authored surface STAGE/FLOW, in-memory series/rating, normal-flow boundaries and copied overrides use bounded stages, shared expiry and component continuation, with accepted rain reports; managed standalone RK2 sessions preserve native weather/calendar state, pending batches, copied prescriptions and numerical, binary or combined mesh-HDF5 history; disjoint water-only or qualified native LEGACY pollutant networks also preserve their actual hydraulic/quality state; independent owned degree-day/UEB snow supports water delivery with native or cached weather, full point/ledger restart and all four output modes; zero transport rows, heat disabled, no soil contact or copied weather are required; further heat, wider network quality and shared RK2 network owners remain guarded | Experimental |
| Time stepping | Shared tier ladder, height LTS_TIERS | §9.9 | Implemented |
| Time stepping | Runtime tier count reaching the slowest cell (D-N1) | requiredTiers() exists; the marcher does not grow the ladder | Planned |
| Verification law | Broadbridge–White | Exact Richards benchmark in the test harness; not selectable | Retired |
| Backend | Host OpenMP raw CLOSED_FORM recharge | Optional ABI v1 and separate C/Python operator status; all stores, transport, reactions and heat remain CPU-owned. Other operators and CUDA/HIP/SYCL remain unqualified | Experimental |
Native LEGACY network pollutants now support numerical/binary checkpoint continuation, including signed native loads, current/old concentrations, load statistics and pollutant balances. First-order matched conservative mesh/aquifer rows can share actual surface/node exchange; an independent RK2 mesh can retain the native network quality owner alongside its own integrator. Pollutant order, units, decay, defaults and per-element initial quality participate in model identity. This profile excludes catchments, LIDs, treatment expressions, native heat/age and multispecies reactions. Shared network RK2, ARD/LARD continuation and broader quality owners remain guarded.
Managed global-RK2 qualification also includes skewed all-quad and mixed meshes with heterogeneous aquifers, bulk or independent chemical/thermal grids, groundwater HEAD/FLOW and surface stage/flow/normal/rating boundaries, reactions, age and optional latent heat/cached weather. Independent water/species/energy budgets and numerical/combined-output continuation pass together. These tests qualify the stated integration and restart profile, without extending its borrowed-owner or device restrictions or claiming general spatial order.
A [2D_AQUIFER] row is *, TAG name or CELL n, then five positional columns KS ZS THETA_S THETA_R ALPHA, then keyword pairs (so a fifth law never renumbers a file). Values are authored in project units, converted once:
| Column or keyword | Symbol | US units | SI units | Default |
|---|---|---|---|---|
KS | \(K_s\) | in/hr | mm/hr | 1e-5 m/s internally |
ZS | \(z_s\) | ft | m | 5 m |
THETA_S, THETA_R | \(\theta_s, \theta_r\) | — | — | 0.45, 0.10 |
ALPHA | \(\alpha\) | 1/ft | 1/m | 2 /m |
PSI_B, LAMBDA | \(\psi_b, \lambda\) | ft, — | m, — | 0.2 m, 0.4 |
N, L | \(n, L_M\) | — | — | 1.6, 0.5 |
C_LOSS | \(c_{loss}\) | in/hr | mm/hr | 0 |
HG0 | initial \(h_g\) | ft | m | 0 (empty aquifer) |
KC, DC, AREA (node rows) | \(K_c, d_C, A_b\) | in/hr, ft, ft² | mm/hr, m, m² | 0, 0, cell area |
Validation requires \(K_s, z_s, \alpha, \lambda > 0\), \(0 < \theta_s \le 1\), \(0 \le \theta_r < \theta_s\), \(n > 1\), \(\psi_b, c_{loss} \ge 0\), and a positive DC wherever KC is given.
Shared with Chapter 5. The tables of §5.5 apply directly where the parameter is the same physical quantity:
| This chapter | Chapter 5 | Where |
|---|---|---|
| \(\theta_s\) | porosity \(\varphi\) | §5.5.1 and its soil-texture tables and pedotransfer relations |
| \(\theta_r\) | nearest analogue is the wilting point \(\theta_{WP}\) (an upper bound on residual content) | §5.5.1 |
| \(K_s\) | saturated hydraulic conductivity \(K_S\) | §5.5.2 |
| \(z_s\) | \(E_G - E_B\), ground minus aquifer-bottom elevation | §5.2 |
HG0 | initial water-table depth \(d_L\) | §5.2 |
| \(c_{loss}\) | the deep-percolation coefficient \(DP\) of (5-21) | §5.3.4 |
Not carried over. The percolation equation's conductivity and tension slopes (§5.5.2) are replaced by the soil law's own \(K(\psi)\); the ET coefficients \(UEF\) and \(DEL\) (§5.5.3) by the boundary condition of §9.8; the discharge constants \(A1, B1, A2, B2, A3\) (§5.5.4) by the head-driven exchange of §9.6. Calibrated values of these do not apply.
New parameters. \(\alpha\) is the Gardner/Russo sorptive number, the inverse capillary length (order 1 to 10 /m for sands, below 1 /m for clays; a 5 m column with \(\alpha = 2\) /m starts at \(\alpha L = 10\) and resolves to SIGMA under AUTO). Brooks–Corey's \(\psi_b, \lambda\) come from engineering soil tests, van Genuchten's \(\alpha, n\) from Rosetta or UNSODA, and \(K_c, d_C\) from a MODFLOW River package of the same area.
| Concern | Where |
|---|---|
| State, options, rows, ledger | src/engine/2d/subsurface/SubsurfaceData.hpp, .cpp (SubsurfaceState, GwOptions, GwAquiferRow, GwNodeBed, GwLinkShare; the three storage views and the residual) |
| Soil laws, equilibrium storage, recharge, stress | src/engine/2d/subsurface/SoilCharacteristic.hpp, .cpp |
| The σ column | src/engine/2d/subsurface/SigmaColumn.hpp, .cpp (advanceColumn, columnDtLimit, seedHydrostatic, columnStorage) |
| The kernel | src/engine/2d/subsurface/SubsurfaceSolver.hpp, .cpp — initialize(); the hooks refreshDtCell(), assignTiers(), fireGwFaces(), fireGwCells(), settle(); sampleNodeExchange(); bookInfiltrationFromSurface(), bookLinkSeepage(), takeToSurface(); the per-cell fireCell() |
| Grammar, units, resolution, writer | src/engine/2d/subsurface/SubsurfaceSections.cpp — [2D_AQUIFER_OPTIONS], [2D_AQUIFER], [2D_AQUIFER_NODE]; gwUnitFactors(), resolveSubsurface(), resolveLinkSeepage(); documented in [2D_AQUIFER_OPTIONS], [2D_AQUIFER] and [2D_AQUIFER_NODE], with the GROUNDWATER and GW_ET keys in [2D_OPTIONS] of the input-file chapter |
| Wiring into the router and marcher | src/engine/2d/SurfaceRouter2D.cpp initialize() (enable, resolve, setSubsurface, one-owner checks), coAdvanceStep() (node volumes to coupling_volume), advancePostRouting() (seepage booking); src/engine/2d/solver/ExplicitInertialSolver.cpp syncAndRebuild() (dt0 fold, tiers, pending-cell pinning), runMacroCycle() (hook order), fireCellsImpl() (infiltration down, return up) |
| Report and results | src/engine/plugins/DefaultReportPlugin.cpp (2D Aquifer Continuity); src/engine/2d/output/Default2DOutputPlugin.cpp (Mesh2_face_gw_*, groundwater_ledger, groundwater_node_exchange_cum) |
| C API | include/openswmm/engine/openswmm_gw2d.h (swmm_gw2d_*, SWMM_GW2D_LED_*); src/engine/2d/api/ApiGw2D.cpp |
The design records are plans/TWO_ZONE_GROUNDWATER_FV_INTEGRATION_PLAN.md (manuscript summary, the four feedbacks) and plans/TWO_ZONE_GROUNDWATER_EXPLICIT_LTS_PLAN_2026-08-15.md (the explicit formulation, the σ column, the LTS guarantees). Where this chapter and those documents differ, the chapter follows the code, and the differences are named where they arise: the specific yield \(\theta_s - \theta_{bot}\) in place of \(\theta_s\) (§9.4.4), the closure thresholds without the \(\alpha L \approx 3\) warning or the ENSLAVED refusal (§9.4.1), the linear deep loss without a [GWF] DEEP expression (§9.5.3), the fixed ladder height in place of the runtime tier count (§9.9), and MODE PER_SUBCATCH as implemented (§9.1).
1.13.2