 |
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)
|
______________________________________________________________________________
Chapters 3, 4 and 8 all solve one-dimensional flow along a conduit. When a sewer surcharges, the water that leaves the network spreads over the street, ponds behind a kerb, follows the terrain rather than the pipe, and re-enters the network somewhere else. Chapter 2's node-link model has no representation for any of that: SWMM 5 either discards the overflow or holds it in a fictitious ponded area above the node, to be returned to the same node later.
This chapter documents the optional two-dimensional overland-flow domain OpenSWMM provides for that water — a cell-centred finite-volume solver for the local-inertial shallow-water equations on an unstructured triangular mesh, coupled bidirectionally to the node-link network.
The module is activated by the presence of a mesh in the project ([2D_VERTICES] and [2D_TRIANGLES], inline or via [2D_MESH_FILE]) rather than by a routing option. IGNORE_2D YES disables it while leaving the mesh parsed and editable. The 1D network continues to route by whatever FLOW_ROUTING selects; the 2D domain supplements the network routing rather than replacing it.
A surface that is a domain in its own right. Ponded volume in the node-link model belongs to the node it came from. On the 2D mesh it belongs to the terrain: it flows downhill, splits at a crown, pools in a sag, and drains into whichever inlet it reaches. Which node receives the water is an outcome of the calculation rather than an input to it.
Figure 9-1 places the two domains side by side: the node–link network of Chapters 3, 4 and 8 with its routing and pressurisation alternatives, the overland surface of this chapter with its three momentum closures, and the exchange pathways that join them — vertex and cell coupling, street inlet capture, and the node–bed exchange of a mesh aquifer.
Figure 9-1 The one-dimensional network and the two-dimensional surface: exchange pathways
Flow paths the network does not contain. Overland routes — a road acting as a channel, flow over an embankment, a flow path between two otherwise unconnected catchments — exist in the terrain and nowhere in the pipe topology.
A flood extent. The depth field over a real surface is the quantity flood mapping needs. A ponded volume at a node is not one.
Several limitations should be noted. The method does not replace hydrology: unless RAINFALL_MODE says otherwise the mesh receives the same rainfall the subcatchments do, and both would deliver it (§9.8). It does not infiltrate — the mesh has no soil column, and water on it leaves only by flowing away, evaporating, or entering the network (§9.12). And its momentum equation is an approximation of the full shallow-water system: §9.2 states what is dropped and §9.10 measures what that costs.
Depth-averaged mass conservation over a surface of bed elevation \(z(x,y)\) and free surface \(\eta = z + h\), with unit-width discharge \(\mathbf{q} = h\mathbf{u}\) (m²/s):
\[\frac{\partial h}{\partial t} + \nabla \cdot \mathbf{q} = i - e + s\] | Continuity | (9-1) |
where \(i\) is rainfall intensity, \(e\) the evaporation rate and \(s\) the exchange with the 1D network, all as velocities normal to the surface.
Three momentum closures supply \(\mathbf{q}\) for the same continuity equation, selected by [2D_OPTIONS] MOMENTUM_EQUATION. They share one time-stepping engine — the tiered local time stepping, the flux-active set, the positivity share, the network coupling and the transport of §9.5–§9.7 — and differ only in the face law that supplies the discharge across each interior face and each boundary edge. The default, LOCAL_INERTIAL, reproduces the marcher as it stood before September 2026 bit for bit; the other two are dispatched by one branch per face firing and otherwise reuse every array, ledger and rule of §9.5.
The default momentum equation is the local-inertial (or inertial-wave) approximation of Bates et al. (2010):
\[\frac{\partial \mathbf{q}}{\partial t} + g\,h\,\nabla\eta + \frac{g\,n^{2}\,\lvert\mathbf{q}\rvert\,\mathbf{q}}{h^{7/3}} = 0\] | Momentum | (9-2) |
with \(g = 9.80665\) m/s² and \(n\) Manning's roughness. Compared with the full shallow-water momentum equation, the **convective acceleration term \(\nabla\cdot(\mathbf{q}\mathbf{q}/h)\) is dropped**. Everything else — local acceleration, the pressure gradient written as a free-surface slope, and bed friction — is retained.
This is the single most consequential modelling decision in the chapter, and it is a deliberate one. Retaining \(\partial\mathbf{q}/\partial t\) is what separates (9-2) from the diffusive wave: the diffusive wave has no inertia at all, so a flood front on a flat surface propagates at a speed set by the numerical step rather than by the physics, and the scheme becomes stiff exactly where floods are shallow. Dropping the convective term is what separates (9-2) from the full system, and it costs the Bernoulli terms: the scheme cannot represent a drawdown over a crest, a stable hydraulic jump position, or a fully supercritical profile. §9.10 shows both effects measured against closed-form solutions.
The practical justification is that urban flood flows are dominated by gravity and friction. Bates et al. (2010) and de Almeida and Bates (2013) delimit the regime: the local-inertial approximation is accurate for subcritical flows over gentle slopes at Froude numbers below about 0.5, and degrades progressively as \(Fr \to 1\). Flow that is persistently supercritical, or where the momentum flux through a contraction sets the answer, is outside it.
For steep faces the model would accelerate without bound, since nothing in (9-2) limits the velocity a slope can generate. A Froude clamp supplies that limit numerically (§9.5.1). The prognostic state is the cell volume \(V\) and one face-normal discharge \(q_f\) per interior face; the face law is (9-7)–(9-10).
FULL_SWE)MOMENTUM_EQUATION FULL_SWE restores the convective term and solves the shallow-water momentum equation in conservation form,
\[\frac{\partial \mathbf{q}}{\partial t} + \nabla\cdot\left( \frac{\mathbf{q}\,\mathbf{q}}{h} + \tfrac{1}{2}\,g\,h^{2}\,\mathbf{I} \right) = -\,g\,h\,\nabla z - \frac{g\,n^{2}\,\lvert\mathbf{q}\rvert\,\mathbf{q}}{h^{7/3}}\] | Momentum, conservative form | (9-31) |
which is (9-2) with the convective flux \(\mathbf{q}\mathbf{q}/h\) kept and the pressure gradient \(g h \nabla\eta\) split into the hydrostatic pressure flux \(\nabla(\tfrac{1}{2} g h^{2})\) and the bed-slope source \(-g h \nabla z\). The prognostic state changes with the equation: each cell carries its volume and a cell momentum \((hu, hv)\), and faces no longer hold a prognostic discharge — they evaluate a Godunov flux from the two adjacent cell states, after a hydrostatic reconstruction of the two depths against the higher of the two beds (Audusse et al., 2004), with a rotated HLLC Riemann solver (Toro, 2001) and a per-face bed-slope correction that cancels the pressure flux exactly at rest. Manning friction is applied per cell after the flux update, in the same semi-implicit division the local-inertial face law uses, so it can only shrink \(\lvert\mathbf{q}\rvert\) and never reverse it (the stopping condition of Liang and Marche, 2009, is satisfied by construction). §9.5.10 gives the face flux step by step.
What the term buys is every phenomenon §9.2.1 said the default cannot represent: transcritical control over a crest, a hydraulic jump held at its Rankine–Hugoniot position, dam-break fronts at the right speed, the Bernoulli drawdown of §9.10. The closure is valid at every Froude number and needs no Froude clamp; FROUDE_MAX is not consulted.
What it costs is a smaller step and more work per face. First-order Godunov with hydrostatic reconstruction is positivity-preserving at a Courant number of one half on triangles and quadrilaterals, so CFL_NUMBER is reduced to 0.5 with a notice when set higher, and the Courant length is the isoperimetric \(2A/P\) rather than the operator-derived \(L_{char}\) of (9-4); the \(\beta\) share of (9-15) remains as a backstop rather than the guarantee. Measured on the SWASHES strips of §9.10 the closure runs at roughly two to three times the wall-clock of the default. Two further switches belong to it: RECONSTRUCTION_ORDER 2 adds coupled MUSCL reconstruction of \((\eta,h,u,v)\): Green–Gauss in connected wet interiors and a conditioned wet-neighbour depth/surface fit at shorelines, with limited face depths, cell-centred shoreline velocity and SSP-RK2 time stepping in global-step mode (§9.5.6, §9.5.10), and is refused with a notice when overland transport is active; FRONT_REBUILD (§9.5.7) is switched on automatically, because a Godunov front outruns the fixed rebuild cadence. The Kokkos plugin backends serve only the all-triangle local-inertial scheme; a FULL_SWE model is announced and run on the CPU marcher (§9.11.1).
FULL_SWE keeps every property of §9.5: the lake at rest is exact — the Audusse correction cancels the \(\tfrac{1}{2} g h^{*2}\) pressure flux face by face, including at walls, which contribute a mirror-state Riemann flux — a dry higher neighbour is a wall, \(\sum V\) closes to round-off across tiers, and species ride the mass flux upwinded on the contact wave. Measured on the SWASHES strips of §9.10 (relative \(L^{1}\) depth error, local-inertial in parentheses): Stoker wet dam break 0.7 % (6.5 %), 0.3 % at second order; Ritter dry dam break 1.7 % (13 %), 1.1 % at second order; subcritical bump 0.5 % (6.8 %); transcritical bump 4.4 % (27 %); bump with shock 7.3 % (12 %). The transcritical and shock residuals are the 25-cell crest resolution of those decks — the 1D finite-volume solver of Chapter 8 sits at 7.5 % and 2.5 % on the same cases.
MOMENTUM_EQUATION DIFFUSIVE_WAVE drops both accelerations, so the momentum equation reduces to the quasi-steady Manning balance of the free-surface slope (Hunter et al., 2005; the explicit diffusive wave of LISFLOOD-FP):
\[\mathbf{q} = -\,\frac{h^{5/3}}{n}\,\frac{\nabla\eta}{\sqrt{\lvert\nabla\eta\rvert}}\] | Momentum, diffusive wave | (9-32) |
The prognostic state is the cell volume alone; a face discharge is recomputed from the two adjacent surfaces at every firing and carries no memory. Substituting (9-32) into (9-1) gives a nonlinear diffusion equation, \(\partial h/\partial t = \nabla\cdot(K\nabla\eta)\) with \(K = h^{5/3}/(n\sqrt{\lvert\nabla\eta\rvert})\), and that is both the closure's merit and its cost. Steady uniform flow is exact, slow floodplain inundation is well represented, and there is nothing to damp: no \(\theta\) blend and no Froude clamp are applied. But an explicit integration of a diffusion equation is bounded by \(\Delta t \le L^{2}/(4K)\) — the \(\Delta x^{2}\) restriction that the local-inertial term of §9.2.1 exists to remove — and the conductance \(K\) diverges on a flat surface, the classic explicit diffusive-wave failure. §9.5.11 gives the face law with its slope floor and the per-cell step bound. The closure is affordable here only because the tiered local time stepping of §9.5.6 lets each cell pay its own step.
Table 9-2 summarizes the three closures and the retired spelling that preceded FULL_SWE.
| Closure | Prognostic state | Face law | Valid regime | Cost | Status |
|---|---|---|---|---|---|
LOCAL_INERTIAL (default) | \(V\) per cell, \(q\) per face | de Almeida and Bates, (9-7)–(9-10) | \(Fr \lesssim 0.5\): ponds, streets, floodplains | 1× | Implemented |
FULL_SWE | \(V\), \(hu\), \(hv\) per cell | hydrostatic reconstruction (9-36), rotated HLLC (9-38), bed-slope correction (9-39), semi-implicit friction per cell (9-40) | all Froude numbers: transcritical control, jumps, dam breaks, drawdown over crests | ≈ 2–3× (Courant ≤ ½ on \(2A/P\), (9-41)); CPU marcher only | Implemented |
DIFFUSIVE_WAVE | \(V\) per cell | Manning quasi-steady balance (9-42) with a slope floor | slow floodplain inundation; steady uniform flow (exact) | \(\Delta x^{2}\)-bound steps (9-43), carried by the LTS tiers; CPU marcher only | Implemented |
ADVECTION YES | as LOCAL_INERTIAL | (9-7) plus a staggered upwind convective term (Stelling and Duinmeijer, 2003) | — | — | Retired |
Table 9-2 Momentum closures of the 2D marcher
ADVECTION YES is the deprecated (2026-09-06) spelling of the convective-term experiment that FULL_SWE superseded. A deck that still carries it is accepted with a warning naming the replacement and keeps the physics it asked for — the local-inertial law with the upwind convective difference added to (9-7) on wet–wet faces — rather than being remapped, so an old calibrated model does not change its answer silently. New models should not use it.
Unlike the internal discretization of Chapter 8, the mesh is an explicit part of the model. Its cells appear in results, carry per-cell parameters, and are addressable by the API.
Cells are triangles, listed in [2D_TRIANGLES] by three vertex indices from [2D_VERTICES], with Manning's \(n\) and an optional initial depth per cell:
Bed elevation is carried at the vertices. A cell's centroid elevation \(z_c\) is the mean of its three vertex elevations, so the cell bed is the plane through them — a fact the closure of §9.4 uses directly. Nothing in the solver reads a cell-constant bed.
Edge–neighbour adjacency is built by hashing each triangle's three sorted vertex pairs: an edge claimed twice is interior and the two triangles become neighbours; an edge claimed once is a domain boundary and carries a boundary condition (§9.6). The convention throughout is that **local edge \(e\) is opposite vertex \(e\)**, so its endpoints are the triangle's other two vertices. Both incident cells of an interior edge therefore compute the same endpoint elevations, which is what makes the face depth of §9.5.2 antisymmetric and the flux conservative.
Interior edges are then enumerated once each as unique faces, with the left cell \(L\), the right cell \(R\), the edge length \(\xi\), the outward normal of \(L\)'s slot (which points \(L \to R\)), the edge midpoint, and the face-normal centroid separation
\[d_{n} = \max\left( \lvert (\mathbf{x}_{R} - \mathbf{x}_{L})\cdot\hat{\mathbf{n}} \rvert,\ 0.3\,\lvert \mathbf{x}_{R} - \mathbf{x}_{L}\rvert \right)\] | (9-3) |
The floor keeps a near-degenerate sliver — where the centroid chord is nearly parallel to the shared edge — from producing an unbounded surface slope. Face roughness is the mean of the two cells' \(n\).
Each cell also carries a characteristic length derived from the discrete operator rather than from geometry, for the time-step bound of §9.5.5:
\[L_{char} = \sqrt{\frac{2A}{\sum_{f} \xi_{f}/d_{n,f}}}\] | (9-4) |
A cell none of whose edges is an interior face has an empty sum in (9-4) and keeps the altitude proxy \(2A/\xi_{max}\) instead; such a cell carries no interior flux until a neighbouring face opens, so the proxy never governs a coupled update.
Every quantity above is planimetric. Areas, lengths and normals are computed in plan rather than on the sloped surface. On terrain steep enough for that distinction to matter the shallow-water equations are themselves the wrong model, so this is a consistent rather than an incidental choice.
Per-edge conveyance factors in \([0,1]\) may be attached with [2D_EDGE_CONVEYANCE], multiplying the flux across that edge — the transmissivity \(\psi\) of the integral-porosity shallow-water literature (Sanders et al., 2008; Bruwier et al., 2017), used to represent sub-grid obstructions such as walls and fences without meshing them. Values are mirrored onto both slots of an interior edge, so the restriction stays symmetric and the flux stays conservative.
Initial conditions come from INIT_DEPTH in [2D_TRIANGLES] (a depth in mesh units, converted to a cell volume through the closure of §9.4) and optionally [2D_INITIAL_VELOCITY], which seeds face momentum from the cell velocities. Without the latter a model can only start from rest, which excludes solutions such as Thacker's oscillating basins whose initial state has \(\mathbf{u} \neq 0\).
For a finite local mesh bed, [2D_INITIAL_BED_QUALITY] seeds the bed porewater separately from the overland water. It uses the same *, TAG and one-based CELL scopes and * < TAG < CELL precedence as [2D_INITIAL_QUALITY]; each concentration is multiplied by that cell's bed pore volume (area × bed thickness × porosity). Pollutants use their declared concentration units, and __WATER_AGE__ is authored in hours. This section requires coupled heat with BED_HEAT_MODEL LOCAL_STORAGE. Temperature remains owned by the heat state. For groundwater-owned GROUND_COLUMN porewater, seed SAT or UNSAT through [GW_INITIAL_QUALITY] instead.
The mesh may mix triangles with convex quadrilaterals. Quads are listed in [2D_QUADS] after every triangle, by four vertex indices in cyclic order (either orientation — the builder orients every edge normal outward from the centroid), with the same roughness, initial-depth and tag columns as [2D_TRIANGLES]:
Cells are numbered triangles first, then quads, each in file order: the \(j\)-th quad row is cell \(n_{tri} + j\). Every cell-addressed section — [2D_BOUNDARY_CONDITIONS], [2D_BOUNDARY_QUALITY], [2D_INITIAL_VELOCITY], [2D_TRIANGLE_NODE_MAP], [2D_INFILTRATION] — uses this unified index, and their column headers keep the word TRI. A [2D_TRIANGLES] row appearing after a [2D_QUADS] row is a parse error, which is what keeps the numbering unambiguous. A non-convex or self-intersecting quad fails mesh validation with the cell index and vertex list; the GUI mesher must never emit one.
Every per-cell and per-edge array is padded to four slots per cell (kMaxCellVerts), so slot \(4t + k\) addresses local edge \(k\) of cell \(t\) whatever its vertex count \(n_v\); a padding slot has zero geometry and a neighbour index of \(-2\), so a stray read is inert. **Local edge \(k\) has endpoints \(v_{(k+1) \bmod n_v}\) and \(v_{(k+2) \bmod n_v}\).** For a triangle this is exactly the opposite-vertex rule of §9.3, so every existing TRI EDGE row keeps its meaning; for a quad it is a fixed rotation of the natural numbering — edge 0 is \((v_2, v_3)\), edge 1 \((v_3, v_4)\), edge 2 \((v_4, v_1)\) and edge 3 \((v_1, v_2)\) — and the EDGE column of a boundary row ranges 0..3. The adjacency hash of §9.3 is unchanged: an edge claimed twice is interior whichever shapes claim it, and a face between a triangle and a quad is an ordinary face. All of the face kernels — the three closures of §9.2, the positivity share, the transport — see only \((h_f, S, n_f, \xi, 1/d_n)\) and carry no shape branch.
Three quantities need care on a quad. Its area is the shoelace area, computed as the sum of the two sub-triangles \((v_0, v_1, v_2)\) and \((v_0, v_2, v_3)\). Its centroid is the area-weighted centroid of those two sub-triangles — a property of the polygon, independent of which diagonal is used — and not the vertex mean, which is not the centroid of a skewed quad and would bias the Perot arms of (9-16) and the ghost-cell distance of §9.6. Its mean bed \(z_c\), used by the FLAT closure, is the mean of the four vertex elevations. The centroid-to-edge distance that a triangle's stage boundary wrote as \(2A/(3L)\) is generalized to the normal distance from the centroid to the edge midpoint,
\[d_{e} = (\mathbf{m}_{e} - \mathbf{x}_{c})\cdot\hat{\mathbf{n}}_{e}\] | (9-33) |
which equals \(2A/(3L)\) for a triangle (that expression is the centroid-to-edge distance) and is exact for a quad; the triangle expression is retained verbatim in the solver so an all-triangle mesh stays bit-identical. Two other constants scale with the vertex count: the positivity share of (9-15) becomes \(\beta/n_v\) per exporting face, since a cell's at most \(n_v\) outgoing faces may take at most \(\beta V\) in total, and a quad none of whose edges is an interior face keeps \(L_{char} = 2\min_k d_{e,k}\) instead of the altitude proxy (which is \(\Delta x\) for a square and conservative for a skewed cell). The isoperimetric length \(2A/P\) that bounds the FULL_SWE step (§9.5.10) is computed over all \(n_v\) edges.
Under CELL_CLOSURE VFR a quad uses the two-plane storage model of §9.4.2. Results of a mixed mesh are written as UGRID mixed topology — Mesh2_face_nodes [nFace, 4] with _FillValue = -1, plus a Mesh2_face_nv vertex-count array — and the edge datasets and the bulk C-API edge arrays use stride 4; an all-triangle mesh keeps the [nFace, 3], stride-3 layout byte for byte. The Kokkos plugin backends serve all-triangle meshes only; a mixed mesh is announced and run on the CPU marcher (§9.11.1).
Figure 9-4 A quadrilateral and a triangle sharing a face: local edge numbering, the padded slot layout and the centroid-to-edge distance
The conserved variable is cell volume \(V\) rather than depth. Depth follows as the cell-mean \(\bar{h} = V/A\). The closure is the relation that recovers a free-surface elevation \(\eta\) from that volume, and it is where a triangular mesh with sloping cell beds differs sharply from a raster of flat cells.
CELL_CLOSURE FLAT (the default) uses
\[\eta = z_{c} + \bar{h}\] | (9-5) |
which is exact for a fully wetted cell and wrong for a partially wet one. On a cell whose bed spans a slope or a step, spreading the water uniformly over the whole cell raises the computed surface above the true waterline — by up to two-thirds of the cell's relief. That error is a head, and a head drives flux: the classic symptom is thin films creeping uphill at a shoreline, and a lake at rest that is not a steady state where it meets the bank.
CELL_CLOSURE VFR uses the exact stage–storage relation of the plane bed through the cell's three vertex elevations — the volume/free-surface relationships of Begnudelli and Sanders (2006, 2007). With sorted elevations \(z_1 \le z_2 \le z_3\) and \(\bar{z} = (z_1+z_2+z_3)/3\):
\[\bar{h}(\eta) = \frac{(\eta - z_{1})^{3}}{3(z_{2}-z_{1})(z_{3}-z_{1})}\] | \(z_{1} < \eta \le z_{2}\) | (9-6a) | |
\[\bar{h}(\eta) = (\eta - \bar{z}) + \frac{(z_{3}-\eta)^{3}}{3(z_{3}-z_{1})(z_{3}-z_{2})}\] | \(z_{2} < \eta \le z_{3}\) | (9-6b) | |
\[\bar{h}(\eta) = \eta - \bar{z}\] | \(\eta \ge z_{3}\) | (9-6c) |
The solver needs the inverse \(\eta(\bar{h})\). On the lower branch it is a closed-form cube root; on the upper branch a safeguarded Newton iteration bracketed by \([z_2, z_3]\), which converges unconditionally because \(d\bar{h}/d\eta = A_{wet}/A > 0\) there. Note that (9-6c) is identical to the flat closure (9-5), since \(z_c = \bar{z}\) — the two closures differ only on partially wet cells, which is precisely the claim.
As a cell dries, \(d\eta/dV = 1/A_{wet}\) diverges. VFR_MIN_WET_FRAC ( \(\varepsilon\), default 0.01) bounds it by continuing the exact relation below wetted fraction \(\varepsilon\) along its tangent line, keeping the closure \(C^1\) and monotone with \(d\eta/dV \le 1/(\varepsilon A)\). The regularized forward and inverse relations are exact inverses of each other for the same \(\varepsilon\), so seeding a state from heads and reading it back as volumes round-trips.
VFR is correct and is not the default. It restores the C-property at shorelines and removes the uphill-creep artifact, but in doing so it resolves a wetting and drying front that the flat closure freezes out — measured at three to eight times more solver substeps. For deep urban flooding, where partially wet cells are a thin fringe around a large wetted area, the flat closure is both faster and adequate. For shallow sheet flow on gentle slopes, where most cells are partially wet, VFR is the right choice and the cost is the price of the answer. Pair it with FACE_RECONSTRUCTION VFR_FACE (§9.5.2); the two halves address the same artifact from the cell side and the face side.
The relations (9-6) are evaluated by a single header-only routine shared by the serial marcher, the boundary path and the GPU kernels, so every backend reconstructs the identical surface. Three auxiliary quantities appear. The wetted-area fraction of the planar bed — which is also \(d\bar{h}/d\eta\) — with sorted elevations \(z_1 \le z_2 \le z_3\) and relief \(R = z_3 - z_1\):
\[w(\eta) = \frac{(\eta - z_{1})^{2}}{(z_{2}-z_{1})\,R}\] | \(z_{1} < \eta \le z_{2}\) | (9-20a) | |
\[w(\eta) = 1 - \frac{(z_{3} - \eta)^{2}}{R\,(z_{3}-z_{2})}\] | \(z_{2} < \eta < z_{3}\) | (9-20b) |
with \(w = 0\) below \(z_1\) and \(w = 1\) above \(z_3\). The stage at which the wetted fraction equals the regularization floor \(\varepsilon\):
\[\eta_{s} = z_{1} + \sqrt{\varepsilon\,(z_{2}-z_{1})\,R}\] | \(\varepsilon \le (z_{2}-z_{1})/R\) | (9-21a) | |
\[\eta_{s} = z_{3} - \sqrt{(1-\varepsilon)\,R\,(z_{3}-z_{2})}\] | otherwise | (9-21b) |
Below the switch depth \(\bar{h}_s = \bar{h}(\eta_s)\) from (9-6) the inverse continues along the tangent line,
\[\eta(\bar{h}) = \eta_{s} - \frac{\bar{h}_{s} - \bar{h}}{\varepsilon}\] | \(\bar{h} \le \bar{h}_{s}\) | (9-22) |
and on the lower branch it is closed-form,
\[\eta = z_{1} + \left\lbrack 3\,\bar{h}\,(z_{2}-z_{1})\,R \right\rbrack^{1/3}\] | \(\bar{h} \le (z_{2}-z_{1})^{2}/(3R)\) | (9-23) |
The inverse \(\eta(V)\) is evaluated by the following case ladder, in this order:
The forward relation \(\bar{h}(\eta)\) applies the same regularization — exact above \(\eta_s\), the tangent below, floored at zero — so forward and inverse are exact inverses for the same \(\varepsilon\) and head-seeded and volume-seeded states round-trip. A dry cell's head under VFR is seeded at \(\eta(0) = \eta_s - \bar{h}_s/\varepsilon \in (z_1, \eta_s)\), the value the closure itself returns at \(V = 0\), so that seeding the head back through the forward relation reproduces exactly zero volume. Figure 9-3 (§9.5.9) sketches the wetting cases.
Implementation: the closure lives in src/engine/2d/mesh/VfrClosure.hpp — openswmm::twoD::vfrSort3, openswmm::twoD::vfrWetFraction, openswmm::twoD::vfrStageAtWetFraction, openswmm::twoD::vfrMeanDepthFromEtaExact, openswmm::twoD::vfrMeanDepthFromEta, openswmm::twoD::vfrEtaFromMeanDepth and openswmm::twoD::vfrDryEta. The solver enters through openswmm::twoD::inertial::etaDepthScalar and openswmm::twoD::inertial::volumeFromEtaScalar in src/engine/2d/solver/InertialKernels.hpp, whose FLAT branch is (9-5); the flat-relief guard is kVfrFlatRelief ( \(10^{-9}\) m).
Four vertices do not in general lie on a plane, so the planar-bed relations (9-6) do not apply to a quad directly. Begnudelli and Sanders (2007) budget the storage of a quadrilateral cell by splitting it along one diagonal into two planar sub-triangles, and their piecewise cubic, quadratic and linear stage–storage relations (their Eqs. 8–13, with the coefficient tables of their appendix) are exactly the area-weighted sum of the two planar-triangle closures. The diagonal is chosen from the elevation ordering of the four vertices, \(n_1 \le n_2 \le n_3 \le n_4\), so that the split follows the fold of the surface:
| Case | Topology | Diagonal | Sub-triangles |
|---|---|---|---|
| 1 | \(n_1\) and \(n_4\) diagonally opposite | \(n_1\)– \(n_4\) | \((n_1, n_2, n_4)\) and \((n_1, n_3, n_4)\) |
| 2 | \(n_1\) and \(n_4\) adjacent, \(n_2\) adjacent to \(n_1\) | \(n_2\)– \(n_4\) | \((n_1, n_2, n_4)\) and \((n_2, n_3, n_4)\) |
| 3 | \(n_1\) and \(n_4\) adjacent, \(n_2\) adjacent to \(n_4\) | \(n_3\)– \(n_4\) | \((n_1, n_3, n_4)\) and \((n_2, n_3, n_4)\) |
Table 9-3 Diagonal selection of the two-plane quadrilateral closure (Begnudelli and Sanders, 2007, Cases 1–3)
Ties in elevation keep cyclic order (a stable sort), and any consistent choice is a valid split. With sub-triangle areas \(A_1, A_2\) and each sub-triangle's sorted elevations, the cell-mean depth at stage \(\eta\) is
\[\bar{h}(\eta) = \frac{A_{1}\,\bar{h}_{1}(\eta) + A_{2}\,\bar{h}_{2}(\eta)}{A_{1} + A_{2}}\] | (9-34) |
where \(\bar{h}_{1}, \bar{h}_{2}\) are the regularized planar-triangle relations of §9.4.1 — (9-6) above each sub-triangle's switch stage, the tangent (9-22) below it, floored at zero — evaluated with the same \(\varepsilon\). Quads and triangles therefore share one regularization and one round-trip contract, and the closed-form appendix coefficients are an optimization the engine does not need. The slope of the sum, which the Newton inverse uses, is
\[\frac{d\bar{h}}{d\eta} = \frac{A_{1}\,w_{1}^{*}(\eta) + A_{2}\,w_{2}^{*}(\eta)}{A_{1} + A_{2}}\] | (9-35) |
with \(w_{k}^{*}\) the wetted fraction (9-20) of sub-triangle \(k\) above its switch stage, \(\varepsilon\) inside its linear tail, and zero once that tail has reached zero depth.
The inverse \(\eta(\bar{h})\) follows the same case ladder as §9.4.1, with two differences in the fully wet and dry limits. Let \(\bar{z}_w\) be the area-weighted mean of the two sub-triangle mean beds — the true mean bed of the two-plane surface — and \(z_{top}\), \(z_{low}\) the highest and lowest of the six sorted elevations:
The derivative the coupling needs, \(d\eta/d\bar{h} = 1/\max(w^{*}, \varepsilon)\), is bounded by \(1/\varepsilon\) as on a triangle.
One consequence deserves stating. The FLAT closure (9-5) writes a quad's bed as the vertex mean \(z_c = \tfrac{1}{4}\sum z\), whereas the fully wet branch of the two-plane relation uses \(\bar{z}_w\), the mean bed of the folded surface. The two coincide for a planar quad and differ otherwise, so on a non-planar quad the two closures disagree even when the cell is fully wet, by the difference of the two bed means — the flat closure is exact for a fully wetted cell only when \(z_c\) is the cell's true mean bed, which for a quad the vertex mean is not. Wetting and drying classification is unchanged: the H_MOVE hysteresis of §9.5.7 acts on the cell-mean depth, and the paper's alternative vertex-wise criterion is not adopted — its threshold needed per-roughness tuning in the authors' own tests, and the flux-active set already plays that role.
Implementation: src/engine/2d/mesh/QuadVfr.hpp — quadVfrPrecompute selects the diagonal once per quad and stores the six sorted elevations and the two areas; quadMeanDepthFromEta, quadWetFraction, quadMeanDepthSlope, quadEtaFromMeanDepth and quadDEtaDMeanDepth are (9-34), (9-35) and the ladder above. The solver enters through openswmm::twoD::inertial::etaDepthQuadScalar and openswmm::twoD::inertial::volumeFromEtaQuadScalar in src/engine/2d/solver/InertialKernels.hpp; the quad branch of cellEtaDepth is kept out of line on purpose, because inlining its Newton loop stopped the compiler inlining the triangle path that every cell of every mesh runs. The MeshBuilder computes the per-quad data (quad_vfr_z, quad_vfr_a) after the cell geometry.
The discretization is a cell-centred finite-volume method with a staggered face variable: cells hold volume, faces hold the unit-width discharge \(q\) normal to the face, positive from \(L\) to \(R\). Time integration is explicit, with per-cell local time stepping.
Each face integrates (9-2) along its own normal over its own step \(\Delta t_f\):
\[q^{n+1} = \frac{\hat{q} - g\,h_{f}\,\Delta t_{f}\,S}{1 + g\,\Delta t_{f}\,n_{f}^{2}\,\lvert\mathbf{q}_{f}\rvert / h_{f}^{7/3}}\] | (9-7) |
with the free-surface slope
\[S = \frac{\eta_{R} - \eta_{L}}{d_{n}}\] | (9-8) |
Three details in (9-7) carry weight.
Friction is semi-implicit. Writing the friction term with \(q^{n+1}\) in the numerator and \(\lvert\mathbf{q}\rvert\) from the previous state puts it in the denominator, which is unconditionally stable: however large the friction coefficient, the update can only shrink \(q\) towards zero, never overshoot through it. An explicit friction term would impose a step limit that scales as \(h^{7/3}/n^{2}\) — unusable on thin films, which is where most of the cells are.
The friction magnitude is taken from the flow vector rather than the face-normal component. Manning friction acting on the normal component is \(n^{2} q_{n} \lvert\mathbf{q}\rvert / h^{7/3}\). Using \(\lvert q_n \rvert\) instead makes the damping a face applies depend on the face's orientation relative to the flow — a face at 45° to a uniform sheet under-damps by \(\sqrt{2}\) — so no smooth surface can satisfy every face of a triangulated slope simultaneously, and the steady state corrugates cell to cell. The vector at the face is the mean of the two incident cells' Perot-reconstructed discharge vectors (§9.5.4, Eq. 9-16), and its magnitude is floored at the face's own discharge:
\[\lvert\mathbf{q}_{f}\rvert = \max\!\left( \lvert q_{f} \rvert,\ \left\lvert \tfrac{1}{2}(\mathbf{q}_{L} + \mathbf{q}_{R}) \right\rvert \right)\] | (9-24) |
so a face whose cell reconstruction lags its own discharge — the first firing after a front arrives, or immediately after activation — never under-damps. With \(\theta = 1\) the cell discharge vectors are not allocated at all: the update uses \(\hat{q} = q_f\) and \(\lvert\mathbf{q}_f\rvert = \lvert q_f \rvert\), recovering the original Bates et al. (2010) scheme exactly. The friction exponent is evaluated as \(h^{7/3} = h^{2}\sqrt[3]{h}\) rather than through pow().
** \(\hat{q}\) is a lateral average, not \(q\) itself.** With THETA \(= \theta\),
\[\hat{q} = \theta\,q_{f} + (1-\theta)\,\tfrac{1}{2}\left( \mathbf{q}_{L} + \mathbf{q}_{R} \right)\cdot\hat{\mathbf{n}}\] | (9-9) |
\(\theta = 1\) recovers the original Bates et al. (2010) scheme, which carries no numerical diffusion and is prone to a checkerboard oscillation in thin films on steep faces. \(\theta < 1\) blends in the neighbouring cells' reconstructed discharge and damps it — the weighted formulation of de Almeida et al. (2012). The default \(\theta = 0.8\) applies enough diffusion to suppress the oscillation without visibly smearing fronts.
Smoothing over a fixed time interval. By default, THETA_REFERENCE_STEP = 0: the authored weight is applied on every face firing. Shortening a routing window can create more face firings per second and increase this numerical damping. A positive THETA_REFERENCE_STEP \(T_\theta\), in seconds, instead defines the retention over that interval:
\[\theta_f = \exp\!\left(\frac{\Delta t_f}{T_\theta}\ln\theta\right),\qquad \widehat q = q_f + (1-\theta_f)(q_{reconstructed}-q_f). \]
Here \(\Delta t_f\) is the actual face step, including local-timestep tiers and shortened window tails. For a fixed reconstructed target, subdividing an interval gives the same total relaxation. The implementation evaluates the blend fraction with expm1 to retain precision at small steps; CPU and Kokkos use the same scalar function. THETA = 1 still disables blending and retains the existing no-vector-friction branch described above.
For example, THETA 0.8 with THETA_REFERENCE_STEP 1 retains 80% of the momentum difference over one second, independently of the number of face steps. The equivalent relaxation timescale is \(-T_\theta/\ln\theta\). This is a numerical stabilization control, not a discharge coefficient or a calibration to measured flow. Select the interval explicitly and check timestep refinement, stability, heads and conservation. Omitting it preserves the historical calculation and restart identity; positive values are included in checkpoint compatibility checks. Surface plugins must be rebuilt for GPU ABI v9, which includes the new option layout.
Finally the result is clamped:
\[\lvert q^{n+1} \rvert \le Fr_{max}\,h_{f}\sqrt{g\,h_{f}}\] | (9-10) |
FROUDE_MAX defaults to 1.5. This is the steep-face guard: with no convective term there is nothing in (9-2) to arrest acceleration down a steep face, so the supercritical limit must be imposed rather than resolved. Raising it lets genuinely transcritical cases run (§9.10) at the cost of the guard.
As coded, one face firing over its step \(\Delta t_f\) proceeds in a fixed order: (i) evaluate the face depth (§9.5.2) and wall the face if \(h_f \le\) DRY_DEPTH, zeroing its momentum; (ii) form \(\hat{q}\) by (9-9) and the friction magnitude by (9-24); (iii) zero the free-surface difference if it lies below the \(10^{-12}\) m deadband (§9.5.3) and form the slope (9-8); (iv) apply (9-7) and then the clamp (9-10); (v) rescale by the positivity share (9-15) where it binds; (vi) book \(\pm\Delta M\) by (9-13). The stored face discharge is the post-clamp, post-rescale value, so the momentum a face carries always matches the mass it moved.
The depth in (9-7) is a property of the face, and how it is defined decides when water is allowed to cross. Under FACE_RECONSTRUCTION MEAN (the default):
\[h_{f} = \max(\eta_{L}, \eta_{R}) - \max(z_{c,L}, z_{c,R})\] | (9-11) |
\(h_f \le\) DRY_DEPTH makes the face a wall for that substep and its momentum is zeroed. This is the standard "flow depth above the higher
bed" rule, and it is what makes the scheme handle wetting and drying without regime-switching logic — but its bed is the higher centroid elevation. A thin crest resolved as a line of high vertices — a levee, a kerb, a road crown — has its height diluted by roughly a third when averaged into the flanking centroids, so water crosses it before it reaches it.
FACE_RECONSTRUCTION VFR_FACE uses instead the exact mean depth of the driving surface over the wetted portion of the shared edge, evaluated against the edge's true endpoint elevations \(z_{lo} \le z_{hi}\) (Begnudelli and Sanders, 2007, Eq. 14):
\[h_{f} = 0\] | \(\eta \le z_{lo}\) | (9-12a) | |
\[h_{f} = \frac{(\eta - z_{lo})^{2}}{2(z_{hi}-z_{lo})}\] | \(z_{lo} < \eta \le z_{hi}\) | (9-12b) | |
\[h_{f} = \eta - \tfrac{1}{2}(z_{lo}+z_{hi})\] | \(\eta > z_{hi}\) | (9-12c) |
with \(\eta = \max(\eta_L, \eta_R)\). The quadratic branch matches value and slope at both joins, so overtopping onset is \(C^1\) and the flux does not jump when the waterline crosses the edge. The gate (9-12a) is the substantive part: a cell holding water pooled below the whole shared edge conveys nothing across it. Embankments hold to their real crest, and drainage no longer strands water on slopes.
Both branches are single-sourced and used identically by the interior faces, the boundary edges and the GPU kernels, so all backends agree.
A body of water at rest over arbitrary bathymetry must stay at rest. This is the C-property, and it is not automatic — a scheme that discretizes the bed slope and the pressure gradient separately will generally produce a spurious flux from their imbalance.
Writing the pressure gradient as the free-surface slope (9-8) makes the property structural: at rest \(\eta_L = \eta_R\), so \(S = 0\) exactly, and (9-7) with \(\hat{q} = q = 0\) returns zero for any bed whatsoever. Likewise a dry neighbour standing higher gives \(h_f \le 0\) and the face is a wall — there is no uphill creep to suppress.
One numerical guard is needed. The closure round-trip \(V \to \eta\) introduces rounding noise of order 1 ulp, and the square-root character of the friction balance amplifies it: a persistent \(\Delta\eta \sim 10^{-16}\) m sustains \(q \sim 10^{-6}\) m²/s. A slope below \(10^{-12}\) m is therefore set to exactly zero, far below any physical head, after which the friction denominator decays \(q\) geometrically and rest states are exact rather than merely small. §9.10 measures the result at \(10^{-16}\) relative error on the SWASHES lake-at-rest cases.
A face firing books the identical volume transfer into a per-side accumulator:
\[\Delta M = q^{n+1}\,\xi\,\Delta t_{f}, \qquad \text{acc}_{L} \mathrel{-}= \Delta M, \quad \text{acc}_{R} \mathrel{+}= \Delta M\] | (9-13) |
and a cell firing gathers and clears its own side of each incident accumulator:
\[V^{n+1} = V^{n} + \sum_{f} \text{acc}_{f,i} + \Delta t_{c}\,A\,(i - e + s)\] | (9-14) |
after which \(\eta\) and \(\bar{h}\) are recomputed through the closure of §9.4. Because the two sides of a face are written from the same floating-point product, conservation is exact by construction rather than to within a tolerance — including across a local-time-stepping tier interface, where the two sides apply their halves at different times (§9.5.6). The sum of cell volumes plus pending accumulators is an invariant of the face phase, and the engine can assert it directly (OPENSWMM_2D_MARCHER_CHECK).
Positivity is enforced at face cadence rather than by a post-hoc clamp. A cell has at most three outgoing faces, so capping each exporting face at a share \(\beta/3\) of its exporting cell's volume bounds the total export at \(\beta V\) per cell step without any cross-face coordination:
\[\lvert q^{n+1}\rvert\,\xi\,\Delta t_{f} \le \frac{\beta}{3}\,\frac{V_{exp}}{2^{\,k_{exp} - k_{f}}}\] | (9-15) |
\(\beta\) is exchange_beta (0.8). The tier ratio in the denominator matters: an exporting cell republishes its volume only at its own firings, and a finer face fires \(2^{k_{exp}-k_f}\) times in between, so without dividing the share the repeated takes drain the cell. When a face is rescaled, the same rescaled flux updates both sides, so the cap costs nothing in conservation. A zero floor at the cell update remains as a backstop; with the face caps in place it does not engage.
The cell's discharge vector — needed for the friction magnitude and the \(\theta\) blend — is reconstructed at the cell's own cadence from its face fluxes by the Perot (2000) formula:
\[\mathbf{q}_{i} = \frac{1}{A_{i}}\sum_{f} s_{f}\,q_{f}\,\xi_{f}\,\left( \mathbf{x}_{f} - \mathbf{x}_{i} \right)\] | (9-16) |
with \(s_f = \pm 1\) the outward orientation of face \(f\) for cell \(i\) and \(\mathbf{x}_f\) the edge midpoint.
Each cell's stable step follows from the gravity-wave celerity and its own characteristic length:
\[\Delta t_{i} = \alpha \frac{L_{char,i}}{\sqrt{g h_{i}} + \lvert \mathbf{u}_{i}\rvert}\] | (9-17) |
with \(\alpha =\) CFL_NUMBER, default 0.7, and the base step of a macro cycle \(\Delta t_0 = \min_i \Delta t_i\), further capped by MAX_TIMESTEP. Only active cells wetter than DRY_DEPTH enter the census — a film the solver will not move imposes no constraint — and a fully quiescent active set falls back to \(\Delta t_0 =\) MAX_TIMESTEP. The advective augmentation \(\lvert\mathbf{u}_i\rvert = \lvert\mathbf{q}_i\rvert/h_i\) is evaluated from the Perot vector when the depth exceeds \(10^{-6}\) m and \(\theta < 1\), and taken as zero otherwise.
\(L_{char}\) is the operator-derived length of (9-4) rather than a geometric proxy, and the difference matters. The face update couples cells through \(g h \xi_f/(A\,d_{n,f})\); the worst (odd–even) mode of that operator has eigenvalue \(\lambda = 2(gh/A)\sum_f \xi_f/d_{n,f}\), and the explicit update is linearly stable for \(\Delta t \le 2/\sqrt{\lambda}\), which is exactly (9-4) divided by the celerity. Defining \(L_{char}\) this way makes \(\alpha\) a true Courant fraction: \(\alpha = 1\) is the linear stability limit on any mesh, and the default 0.7 is a uniform 30 % margin. A raster of squares recovers the classical \(c\,\Delta t/\Delta x \le 1/\sqrt{2}\); a union-jack pair of right triangles gets \(0.408\,\Delta x\). The obvious geometric proxy \(2A/\xi_{max}\) returns \(0.707\,\Delta x\) on that same union-jack mesh — an overstatement by \(\sqrt{3}\), and the reason frictionless basins seiched at nominal Courant numbers that looked conservative.
Between full rebuilds the tier lists are frozen while depths keep evolving, so \(\Delta t_0\) is re-minimized every macro cycle. It may be tightened at any time — every tier still satisfies \(\Delta t_i \ge 2^{k}\Delta t_0\) — but growing it requires reassigning tiers and therefore waits for a rebuild.
A flood mesh is heterogeneous by nature: a 0.5 m cell at a coupled manhole and a 20 m cell on a floodplain differ by two orders of magnitude in stable step. Marching the whole mesh at the smallest one wastes almost all of the work.
The solver instead assigns each cell a power-of-two tier \(k\) from the ratio \(\Delta t_i/\Delta t_0\), capped at LTS_TIERS (default 4, allowing an 8× spread; up to 8 tiers, 128×). A macro cycle is \(2^{K-1}\) base substeps; tier \(k\) fires every \(2^{k}\) substeps with \(\Delta t = 2^{k}\Delta t_0\). Within a substep all due faces fire first, then all due cells — faces read the surfaces their incident cells published at those cells' last firings.
A face belongs to the finer of its two incident cells' tiers. It therefore always integrates at the rate the sharper side needs, reading the coarser side's surface frozen since that cell last fired. This is what makes tier interfaces safe without interpolation, and (9-13) is what makes them conservative: the same \(\pm\Delta M\) is booked once and applied by each side at its own firing.
Cells whose forcing changes at the fastest cadence are pinned to tier 0 regardless of their Courant number — boundary cells and cells carrying a 1D coupling point. For these cells the forcing, rather than the celerity, sets the resolution requirement.
The per-cell step that feeds the tier assignment is closure-specific: (9-17) under the local-inertial law, the Courant bound (9-41) on the isoperimetric length under FULL_SWE, and the diffusion bound (9-43) under DIFFUSIVE_WAVE. Each is capped by MAX_TIMESTEP, and the tiers are built from whichever applies.
Global-step mode. RECONSTRUCTION_ORDER 2 under FULL_SWE (§9.5.10) integrates with a two-stage SSP-RK2 scheme over the global active lists, and a two-stage scheme cannot be tiered: LTS_TIERS is reduced to 1 with a notice and every active cell and face fires each substep at \(\Delta t_0\). LTS_TIERS 1 selects the same mode by hand for any closure; it is bit-identical to the tiered path on a mesh where tiering finds nothing to separate.
Front-cadence tiers. Under FRONT_REBUILD (§9.5.7) the dry cells of the widened halo would otherwise be assigned the coarsest tier — a dry cell's (9-17) is unbounded — so a front could enter them only every \(2^{K-1}\) substeps regardless of how wide the halo is. After the tier assignment, every dry active cell is therefore pulled down to the finest tier among its active neighbours, and the pull is propagated across the halo rings by a Gauss–Seidel sweep over the faces (at most five passes, stopping early when nothing changes). The sweep is kept serial on purpose: a Jacobi form reaches a different fixed point inside the ring budget, and on the Bellinge 30-minute slice that difference changed the DIFFUSIVE_WAVE inflow and spill volumes materially for a 4 % saving on the FULL_SWE deck alone. The receiving cells then update at the front's own cadence.
Most of a rain-on-grid mesh is not flowing. A cell is flux-active only above H_MOVE (default 3 mm), with hysteresis: entering cells need \(h_{move} + \delta\), active cells stay until \(h_{move} - \delta\), where \(\delta = \min(1\ \text{mm},\ h_{move}/2)\). Scaling the band with H_MOVE matters on shallow benchmarks — a fixed ±1 mm band made H_MOVE \(= 10^{-4}\) require 1.1 mm to activate, ten times the requested threshold, which freezes wetting fronts in place.
Two rules complete the set. A face flows only when both incident cells are active — a one-sided face would export volume into a cell whose update never runs, and measured as an 18 % basin loss when it was allowed. A one-ring halo around the active set therefore guarantees an advancing front always has an active receiving cell.
Inactive cells are not skipped, they are integrated lazily: rainfall and held coupling accumulate as pure storage over the whole interval since the last synchronization, in one pass, because a cell below H_MOVE has no face flux by construction. This is what makes rain-on-grid over a large dry mesh nearly free. A rainfall rate as such never activates a cell — activation follows only from the accumulated depth crossing the threshold at a rebuild, which is the point of the lazy tier — whereas a nonzero coupling flux activates its cell immediately, as do the pinned boundary and coupling-point cells of §9.5.6.
The active set and the tier assignment are rebuilt every four macro cycles rather than every substep; the cost of the rebuild is \(O(n_{cells})\) and dominated everything else when it ran per routing step on a large mesh.
Front-triggered rebuilds (FRONT_REBUILD). Between rebuilds the active set is frozen, so a wetting front can cross at most one cell ring per rebuild period. Under the local-inertial law that is not a practical restriction (§9.5.9): a fast front implies a small \(\Delta t_0\), so the four-cycle period shrinks with the front's own time scale. A Godunov front under FULL_SWE is different — it travels at about \(2\sqrt{gh}\), one cell per base substep at Courant ½, and crosses the one-ring halo inside a single macro cycle, after which it stalls until the next rebuild. That stall is what pinned the Ritter and Thacker results before the option existed. FRONT_REBUILD replaces the fixed cadence with a breach trigger. The halo is widened to five rings (kFrontHaloRings, enough for such a front to cross a macro cycle of \(2^{K-1}\) substeps), its outer ring — an active cell with an inactive neighbour — is marked as the frontier, and a cell firing that finds a frontier cell at or above the activation depth \(h_{on}\) of the last rebuild raises a breach flag; the next macro cycle then rebuilds immediately instead of waiting out the cadence. The dry cells of the widened halo are pulled to the front's tier cadence (§9.5.6), since a wide halo of cells that fire only at the coarsest tier would stall the front just as surely.
AUTO (the default) switches the trigger on under FULL_SWE only. It stays off under LOCAL_INERTIAL, which keeps that closure bit-identical to its earlier results, and off under DIFFUSIVE_WAVE, whose \(\Delta x^{2}\)-bounded step moves a front a fraction of a cell per macro cycle so the cheap cadence suffices. YES and NO force it either way for any closure.
Figure 9-5 The flux-active set at a wetting front: the default one-ring halo, the front-rebuild halo and the breach that forces an early rebuild
The solver's working set is three structure-of-arrays blocks, sized once at initialization.
Per unique interior face (built by the edge enumeration of §9.3): the incident cells \(c_L < c_R\); the edge length \(\xi\); the unit normal \(\hat{\mathbf{n}}\) oriented \(L \to R\) and the edge midpoint; the reciprocal face-normal centroid separation \(1/d_n\) of (9-3); the squared face roughness \(n_f^2 = \left(\tfrac{1}{2}(n_L + n_R)\right)^2\); the sorted true endpoint bed elevations \(z_{lo} \le z_{hi}\) of the shared edge (for VFR_FACE); and the two flat edge-slot indices \(3t + e\) used to publish fluxes back into the per-cell edge arrays. The prognostic state per face is the discharge \(q\) and the two pending-transfer accumulators \(\text{acc}_L\), \(\text{acc}_R\) of (9-13), plus a face tier.
Per cell: the conserved volume \(V\) and the reconstructed \(\eta\) and \(\bar{h}\); the source rates (rainfall, evaporation demand, coupling flux, all m/s); the Perot discharge vector \((q_{cx}, q_{cy})\), allocated only when \(\theta < 1\); the characteristic length \(L_{char}\) of (9-4); the tier, the active flag and the tier-0 pin flag; and a CSR incidence — for cell \(i\) the incident unique faces with orientation signs \(\pm 1\) — so the continuity gather is a race-free per-cell loop.
CPU transport and reconstruction storage. Private species face accumulators use [face * species_count + species]; cell gathering handles up to eight species in one incidence walk without changing each species' edge-summation order. Public cell-mass storage is unchanged. At second order, the existing gex_/gey_ vectors contain elevation gradients in their first cell-count block and depth gradients in their second block, alongside the four velocity gradient vectors and RK stage buffers. This adds two cell-count arrays of doubles to gradient storage relative to elevation/velocity-only reconstruction. The complete CPU strategy is documented in CHAPTER 7 — Two-Dimensional CPU Implementation and Qualification.
Per vertex: coordinates and bed elevation, plus the reconstruction stencils of §9.9. Per flat edge slot: the published volumetric flux and the boundary-condition tables of §9.6.
One solver advance over a window \([t,\ t + \Delta t]\) executes:
The face and cell passes are OpenMP-parallel with static scheduling and the project's THREADS setting. Each face is written by exactly one iteration and touches only its own accumulator slots, and each cell gathers only its own accumulator sides, so both passes are race-free and bit-identical to serial execution for any thread count. The boundary and coupling loops are serial; they are perimeter- and point-count-sized. The quantity \(\sum_i V_i + \sum_f (\text{acc}_{L,f} + \text{acc}_{R,f})\) is invariant under the face phase and is asserted directly when OPENSWMM_2D_MARCHER_CHECK is set.
Figure 9-2 assembles the co-advance batch of §9.7.3 and the marcher's substep loop into one workflow.
flowchart TD
A[1D routing step completes - node heads current] --> B{Pending span reaches the sync batch}
B -- no --> A
B -- yes --> C[Save surface diagnostic state]
C --> D[Outfall transfers and ledgers already applied each routing step]
D --> E[Refresh rainfall, forcing overrides and boundary values]
E --> F[Marcher advance over the batch span]
F --> G{Rebuild due}
G -- yes --> H[Settle accumulators, lazy sources, active set, tiers, dt0]
G -- no --> I[Tighten dt0 from live depths]
H --> J[Macro cycle of base substeps]
I --> J
J --> K[Fire due faces - depth, blend, update, clamp, positivity, book dM]
K --> L[Fire due cells - gather, sources, closure, Perot]
L --> M[Tier-0 firings - boundary edges and live junction exchange]
M --> N{Batch span filled}
N -- no --> G
N -- yes --> O[Book junction and boundary ledgers; outfalls already booked]
O --> P[Queue exchange volumes for uniform delivery to 1D lateral inflow]
P --> Q[Clear one-shot forcings, reset window accumulators]
Q --> A
Figure 9-2 One 1D–2D co-advance batch and the explicit marcher's substep loop within it (rendered diagram)
Implementation: openswmm::twoD::ExplicitInertialSolver::advance drives the loop; openswmm::twoD::ExplicitInertialSolver::fireFaces, openswmm::twoD::ExplicitInertialSolver::fireCells, openswmm::twoD::ExplicitInertialSolver::syncAndRebuild, openswmm::twoD::ExplicitInertialSolver::settleAccumulators and openswmm::twoD::ExplicitInertialSolver::refreshDt0 implement the numbered steps. The face layout and CSR incidence are built by openswmm::twoD::InertialEdges::build (src/engine/2d/solver/InertialEdges.cpp); the scalar kernels live in src/engine/2d/solver/InertialKernels.hpp; the cell state arrays are openswmm::twoD::SurfaceStateData (src/engine/2d/data/SurfaceStateData.hpp).
The scheme has no regime-switching logic; wetting and drying emerge from a small set of thresholds applied uniformly. Table 9-1 collects them.
| Constant | Value | Origin | Role |
|---|---|---|---|
DRY_DEPTH | 0.001 m | [2D_OPTIONS] | Face-wall depth, friction depth floor, evaporation taper scale, coupling ramp scale, CFL census cutoff |
H_MOVE | 0.003 m | [2D_OPTIONS] | Flux-activation depth |
| \(\delta\) | \(\min(0.001\ \text{m},\ h_{move}/2)\) | derived | Activation hysteresis half-band |
| \(\varepsilon\) | 0.01 | VFR_MIN_WET_FRAC | Wetted-fraction floor of the VFR closure |
| slope deadband | \(10^{-12}\) m | fixed | Free-surface differences treated as exactly zero |
| flat-relief guard | \(10^{-9}\) m | fixed | Cell or edge relief below which the geometry is flat |
| \(\beta\) | 0.8 | fixed | Positivity and exchange availability fraction |
| rebuild cadence | 4 macro cycles | fixed | Active-set and tier refresh period |
Table 9-1 Wetting and drying thresholds and guards of the 2D solver (SI units)
The face-wet criterion. A face conveys only when its flow depth exceeds DRY_DEPTH: under MEAN the depth (9-11), under VFR_FACE the wetted-edge depth (9-12) of the driving surface over the edge's true endpoint beds. A face that fails the test is a wall for that substep and its momentum is set to exactly zero — walls carry no stale discharge into their next wet substep. A face also requires both incident cells active (§9.5.7); faces bordering an inactive cell have their momentum zeroed at the rebuild.
How a dry cell wets. The activation thresholds carry hysteresis,
\[h_{on} = h_{move} + \delta, \qquad h_{off} = \max(0,\ h_{move} - \delta), \qquad \delta = \min(0.001,\ h_{move}/2)\] | (9-25) |
entering cells needing \(h_{on}\) and active cells persisting to \(h_{off}\). A dry cell gains water in one of three ways. Distributed sources (rainfall) accumulate lazily as pure storage; the cell joins the active set at the first rebuild whose census finds its depth at or above \(h_{on}\). A concentrated source — a nonzero coupling flux — activates the cell at the next rebuild regardless of depth, as does membership in the pinned set. A neighbouring active cell activates it through the one-ring halo, after which the shared face joins the face lists; the first term to act on the newly wet cell is then the mass transfer (9-13), driven by the slope term of (9-7) integrated from \(q = 0\) — the friction denominator is near unity at \(q = 0\), so the initial specific discharge after one face step is \(-g\,h_f\,\Delta t_f\,S\).
How a front advances. Between rebuilds the active set is frozen, so a wetting front can cross at most one cell ring per rebuild period (four macro cycles). This is not a practical restriction: a fast front implies a small \(\Delta t_0\) through (9-17), so the rebuild period shrinks with the front's own time scale, and the halo guarantees the front always finds an active receiving cell — a one-sided face is never allowed to fire (§9.5.7).
Shorelines and the lake at rest. A dry neighbour standing higher gives \(h_f \le 0\) under (9-11) and the face is a wall — there is no uphill creep to suppress. Under the VFR closure the partially wet shoreline cell's \(\eta\) is exact rather than biased high (§9.4), and under VFR_FACE the gate (9-12a) blocks conveyance across any edge whose low point stands above the driving surface, so an emerged bump holds a lake at rest exactly. The \(10^{-12}\) m slope deadband removes the closure round-trip noise that would otherwise sustain a \(\sim 10^{-6}\) m²/s residual discharge (§9.5.3); the emerged-bump verification case of §9.10 measures the combined result at \(1.4 \times 10^{-16}\) relative depth error.
Draining and positivity. Every exporting face is capped at the tier-scaled \(\beta/3\) share of its exporting cell's volume, Eq. (9-15), so a cell's at most three outgoing faces can remove at most \(\beta V\) per cell step with no cross-face coordination. Boundary edges and the coupling exchange clamp in volume space — the applied change is recomputed after flooring the provisional volume at zero, so booking matches application exactly. The cell update keeps a plain zero floor as a backstop; with the caps in place it does not engage (a debug build with OPENSWMM_2D_MARCHER_CHECK reports any clamp below \(-10^{-12}\) m³).
The friction depth floor. No separate floor exists: the face-wall test guarantees \(h_f >\) DRY_DEPTH \(= 1\) mm before (9-7) is evaluated, so \(h_f^{7/3} \ge 10^{-7}\) m \(^{7/3}\) and the semi-implicit denominator is always finite. Because friction enters only through the denominator, it can shrink \(q\) toward zero and never reverse it.
Rain and evaporation on dry cells. The evaporation sink is the demand rate tapered by a cubic Hermite ramp below DRY_DEPTH,
\[e_{eff} = e \cdot \sigma\!\left( h/h_{dry} \right), \qquad \sigma(t) = \min(1, t)^{2}\left( 3 - 2\min(1, t) \right)\] | (9-26) |
with \(e_{eff} = 0\) for \(h \le 0\) or \(e \le 0\) — a drying cell cannot evaporate more water than it holds, and negative demand is treated as zero rather than as a condensation source. Rainfall on an inactive cell integrates in a single lazy pass over the whole interval since the last synchronization; the result is floored at zero volume, so a forced evaporation override on a dry cell is harmless.
Figure 9-3 sketches the geometry the rules act on: the three wetting cases of the planar-bed cell (§9.4.1) and the wetted-edge face gate (§9.5.2).
Figure 9-3 Wetting cases of a planar-bed triangular cell and the wetted-edge face gate (placeholder — to be replaced by a final drawing)
Implementation: the face-wet tests are openswmm::twoD::inertial::faceFlowDepth and openswmm::twoD::inertial::faceFlowDepthVfr over openswmm::twoD::inertial::faceDepthFromEta; the positivity share is applied in openswmm::twoD::ExplicitInertialSolver::fireFaces with the scale factor of openswmm::twoD::inertial::positivityScale; the hysteretic activation, halo and lazy source passes are openswmm::twoD::ExplicitInertialSolver::syncAndRebuild and openswmm::twoD::ExplicitInertialSolver::lazySourcesOnly; the evaporation taper is openswmm::twoD::evapSink (src/engine/2d/solver/SurfaceFluxCalculator.hpp). The slope deadband is kEtaDeadband in src/engine/2d/solver/InertialKernels.hpp.
Under MOMENTUM_EQUATION FULL_SWE a face firing evaluates one Godunov flux from the two adjacent cell states — free surface \(\eta\), mean depth \(\bar{h}\) and cell momentum \(\mathbf{q} = (hu, hv)\) — and books mass and momentum on both sides from the same products. The steps, in the order coded:
Hydrostatic reconstruction. Each side's bed is taken consistent with its storage, \(z_L = \eta_L - \bar{h}_L\) and \(z_R = \eta_R - \bar{h}_R\) — the cell mean bed under FLAT, the flat-equivalent bed under VFR, which is what keeps a uniform \(\eta\) an exact rest state under either closure — and both depths are reconstructed against the higher of the two (Audusse et al., 2004):
\[z_{f} = \max(z_{L}, z_{R}), \qquad h_{L}^{*} = \max(0,\ \eta_{L} - z_{f}), \qquad h_{R}^{*} = \max(0,\ \eta_{R} - z_{f})\] | (9-36) |
A side whose cell depth is at or below DRY_DEPTH, or whose reconstructed depth vanishes, carries zero velocity, which removes the \(q/h\) blow-up at fronts; the reconstruction alone already removes the "dry cell standing higher" flux. A face dry on both sides carries nothing and books nothing. The velocities are rotated into the face frame, \(u_n = \mathbf{u}\cdot\hat{\mathbf{n}}\) and \(u_t = \mathbf{u}\cdot\hat{\mathbf{t}}\) with \(\hat{\mathbf{t}} = (-n_y, n_x)\).
Wave speeds. With \(c = \sqrt{g h^{*}}\) on each side, the two-rarefaction estimate (Toro, 2001, §10.4) for a wet–wet face and the exact dry-bed front speeds otherwise:
\[u^{\star} = \tfrac{1}{2}(u_{L} + u_{R}) + c_{L} - c_{R}, \qquad c^{\star} = \tfrac{1}{2}(c_{L} + c_{R}) + \tfrac{1}{4}(u_{L} - u_{R})\] | wet–wet | (9-37a) | |
\[S_{L} = \min(u_{L} - c_{L},\ u^{\star} - c^{\star}), \qquad S_{R} = \max(u_{R} + c_{R},\ u^{\star} + c^{\star})\] | wet–wet | (9-37b) | |
\[S_{L} = u_{L} - c_{L}, \qquad S_{R} = u_{L} + 2c_{L}\] | dry right | (9-37c) | |
\[S_{L} = u_{R} - 2c_{R}, \qquad S_{R} = u_{R} + c_{R}\] | dry left | (9-37d) |
The HLLC flux. On the state \(\mathbf{U} = [h,\ h u_n,\ h u_t]\) with physical flux \(\mathbf{F} = [h u_n,\ h u_n^{2} + \tfrac{1}{2} g h^{2},\ h u_n u_t]\),
\[\mathbf{F}_{f} = \mathbf{F}_{L}\ \ (S_{L} \ge 0); \qquad \mathbf{F}_{f} = \mathbf{F}_{R}\ \ (S_{R} \le 0); \qquad \text{else}\ \ F_{h,n} = \frac{S_{R}F_{L} - S_{L}F_{R} + S_{L}S_{R}(U_{R} - U_{L})}{S_{R} - S_{L}}\] | mass and normal momentum | (9-38a) | |
\[F_{t} = F_{h}\,u_{t}^{up}, \qquad S^{\star} = \frac{S_{L}\,d_{R} - S_{R}\,d_{L}}{d_{R} - d_{L}}, \quad d = h^{*}(u_{n} - S)\] | tangential momentum | (9-38b) |
The HLL middle state serves mass and normal momentum; the tangential momentum is a passive scalar upwinded on the contact speed \(S^{\star}\) (Toro, 2001, Eq. 10.70), which falls back to the mean normal velocity when both sides are at rest. Species ride the mass flux from the side \(S^{\star}\) points away from. The momentum flux is rotated back to \((x, y)\).
The bed-slope correction. Each side receives, per unit face length,
\[\mathbf{c}_{L} = \tfrac{1}{2}\,g\left( h_{L}^{*2} - \bar{h}_{L}^{2} \right)\hat{\mathbf{n}}, \qquad \mathbf{c}_{R} = -\,\tfrac{1}{2}\,g\left( h_{R}^{*2} - \bar{h}_{R}^{2} \right)\hat{\mathbf{n}}\] | (9-39) |
At rest the HLLC flux is \([0,\ \tfrac{1}{2} g h^{*2},\ 0]\) and (9-39) cancels it face by face, so a lake at rest over any bed is an exact steady state — the C-property of §9.5.3 in conservative form. Walls close the same balance: every WALL boundary slot of a wet cell contributes a mirror-state Riemann flux (same depth, reversed normal velocity) that carries no mass but \(\tfrac{1}{2} g h^{2}\) of normal momentum, booked at the cell's own cadence from a per-cell list of wall slots; without it a cell at rest beside a wall would accelerate away from it.
Booking. The per-edge conveyance factor of §9.3 scales the whole flux vector. The positivity share of (9-15), with \(\beta/n_v\) per face, scales mass and momentum together so the two stay consistent; under this closure it is a backstop rather than the guarantee (below). The stored face discharge is the limited mass flux, \(\pm\Delta M\) is booked into the volume accumulators exactly as in (9-13), and the momentum accumulators receive \((-\mathbf{F} + \mathbf{c}_{L})\,\xi\,\Delta t_f\) on the exporter side and \((+\mathbf{F} + \mathbf{c}_{R})\,\xi\,\Delta t_f\) on the receiver side.
The cell update. A cell firing gathers its volume by (9-14) and its momentum by
\[\mathbf{q}^{\,n+1} = \frac{A\,\mathbf{q}^{\,n} + \sum_{f}\text{macc}_{f,i}}{A}\ \Big/\ \left( 1 + \frac{g\,\Delta t_{c}\,n^{2}\,\lvert\mathbf{q}\rvert}{h_{new}^{7/3}} \right)\] | (9-40) |
— the gathered momentum, then semi-implicit Manning friction against the new depth, which can only shrink \(\lvert\mathbf{q}\rvert\). Sources (rain, coupling, infiltration, evaporation) move volume only, so the velocity adjusts with the depth; a cell that fell to DRY_DEPTH or below carries no momentum. The cell momentum stands in for the Perot vector (9-16) of the default closure wherever the marcher needs a cell velocity.
The step bound. First-order Godunov with hydrostatic reconstruction is positivity-preserving at a Courant number of one half on triangles and quadrilaterals, and the bound is enforced rather than assumed:
\[\Delta t_{i} = \alpha\,\frac{2A_{i}/P_{i}}{\sqrt{g h_{i}} + \lvert\mathbf{u}_{i}\rvert}, \qquad \alpha \le \tfrac{1}{2}\] | (9-41) |
with \(P_i\) the cell perimeter. CFL_NUMBER above ½ is reduced to ½ with a notice. The isoperimetric length replaces \(L_{char}\) of (9-4) because the positivity argument is about the volume a cell's faces can remove in one step, not about the odd–even mode of the staggered operator.
Second order (RECONSTRUCTION_ORDER 2). The CPU route reconstructs free-surface elevation and depth independently, as well as velocity where the wet stencil permits it. The bed used by this flux is the cell-equivalent bed \(z_i=\eta_i-h_i\). With VFR storage this is an algebraic equivalent of the stored mean depth, not a replacement for the geometric bed used by the volume–stage closure. At face midpoint \(\mathbf{x}_f\), with \(\mathbf{r}_{if}=\mathbf{x}_f-\mathbf{x}_i\),
\[\eta_{if}=\eta_i+\nabla\eta_i^{\ell}\cdot\mathbf{r}_{if},\qquad h_{if}=h_i+\nabla h_i^{\ell}\cdot\mathbf{r}_{if},\qquad z_{if}=\eta_{if}-h_{if}. \]
The superscript \(\ell\) denotes the limited gradient. Reconstructing only elevation over a constant cell bed gives the wrong face-depth variation on a slope. Reconstructing both quantities permits a smooth sloping bed and a smooth water surface while keeping the first-order limit when all gradients vanish.
Connected wet stencil. A neighbour contributes only when it is active, its depth exceeds DRY_DEPTH, and both cell surfaces stand more than DRY_DEPTH above the equivalent-bed sill \(\max(z_i,z_j)\). A wet cell behind a dry sill is not a sample of the same connected water surface. Boundary faces have no neighbour in the interior-face incidence list. Cells with \(h_i\le10\,\texttt{DRY_DEPTH}\) retain zero gradients.
For an interior without rejected neighbours, each of the four variables uses Green–Gauss integration, with the average of the adjacent cell values at an interior face and the cell's own value at a physical boundary. A rejected dry, inactive or sill-blocked neighbour selects the shoreline path instead. That path fits only \(\eta\) and \(h\) to connected wet neighbours. For centroid offsets \(\mathbf{d}_{ij}\),
\[M=\sum_j\frac{\mathbf{d}_{ij}\mathbf{d}_{ij}^{T}} {\lVert\mathbf{d}_{ij}\rVert^2},\qquad \mathbf{b}_w=\sum_j\frac{\mathbf{d}_{ij}(w_j-w_i)} {\lVert\mathbf{d}_{ij}\rVert^2},\qquad \nabla w=M^{-1}\mathbf{b}_w. \]
At least two connected neighbours are required, with \(\det M>0.01(\operatorname{tr}M)^2\). This scale-independent test rejects nearly collinear stencils (a spectral condition number of about 98 is the upper limit). An insufficient or rejected stencil stays first order. The shoreline velocity gradients remain zero even when the surface/depth fit succeeds. Extrapolating velocity from the remaining wet cells can export water faster than momentum and accelerate the residual film; keeping the exporting cell's velocity avoids that failure observed in the distorted-mesh review.
Limiting and face-depth positivity. Barth–Jespersen factors constrain each variable using connected wet neighbours and their face midpoints. At shoreline or physical-boundary cells, a further check includes every face midpoint, including faces omitted from the wet fit. For the already limited depth increment \(\delta h_{if}\),
\[\theta_i=\min\left(1,\min_{\delta h_{if}<-h_i} \frac{-h_i}{\delta h_{if}}\right),\qquad (\nabla h_i^{\ell},\nabla\eta_i^{\ell})\leftarrow \theta_i(\nabla h_i^{\ell},\nabla\eta_i^{\ell}). \]
An empty inner minimum imposes no additional reduction. Scaling both gradients together preserves their equality for a flat reconstructed bed. This is a slope restriction before the Riemann solve; simply clipping an extrapolated negative depth would define a different polynomial. The hydrostatic cut and available-volume flux bound remain separate safeguards.
Figure 9-7 Coupled surface/depth reconstruction and the connected-wet-stencil decision at a shoreline. The drawing is schematic, not model output.
Pressure and bed-force consistency. Hydrostatic cutting uses \(z_f^*=\max(z_{Lf},z_{Rf})\) and \(h_i^*=\max(0,\eta_{if}-z_f^*)\). The scalar correction for side \(i\), before multiplication by its outward normal, is
\[C_i=\frac{g}{2}\left[(h_i^*)^2-h_{if}^2 -(h_i+h_{if})(z_{if}-z_i)\right]. \]
The last term integrates the bed force with the mean cell/face depth along the centroid-to-face segment. It vanishes on a flat bed. At constant elevation, \(z_{if}-z_i=h_i-h_{if}\), so the expression reduces to \(g[(h_i^*)^2-h_i^2]/2\), the first-order hydrostatic correction. Thus a reconstructed cell can meet a first-order cell without losing the lake-at-rest cancellation. The common Riemann mass flux is still booked once with opposite signs; these side-specific pressure terms supply the bed force in momentum, not an additional water transfer.
Two stages and their ledgers. With volume and momentum state \(\mathbf U\), SSP-RK2 (Heun) advances \(\mathbf U^1=\mathbf U^0+\Delta t L(\mathbf U^0)\), \(\mathbf U^2=\mathbf U^1+\Delta t L(\mathbf U^1)\), then \(\mathbf U^{n+1}=(\mathbf U^0+\mathbf U^2)/2\) over the global active lists. Both stages recompute gradients and fluxes. Volume–stage closure is refreshed after the averaged state is formed. Source and exchange ledgers must receive the same trapezoidal weighting as the state: boundary/coupling volumes, infiltration and node spill budgets must not count two complete time intervals. Friction is applied through the cell momentum update.
Requested second order reduces LTS_TIERS to 1 with a notice. Active surface transport falls back to first order with a notice; the tier reduction occurs before that fallback, so request order 1 explicitly when transport should retain multiple tiers. Standalone CPU passive pollutant/age CLOSED_FORM/ENSLAVED/SIGMA groundwater supports the shared RK2 transaction, authored wells and HEAD/FLOW edges. Groundwater reactions/heat and network/catchment/runtime owners remain guarded; see Hydrology §9.11 and the groundwater transport design for the exact profile.
Implementation: src/engine/2d/solver/SweKernels.hpp — faceFlux (first order), faceFluxReconBed (second order; faceFluxRecon retains the constant-bed compatibility entry point), hllcFlux, waveSpeeds, bjLimiter, frictionUpdate; the face firing is openswmm::twoD::ExplicitInertialSolver::fireFacesSwe, the momentum gather and wall booking are the FULL_SWE branch of fireCells, the gradients are computeLimitedGradientsSwe and the two-stage step is runRk2Step in src/engine/2d/solver/ExplicitInertialSolver.cpp; the perimeter length is cell_lpos in InertialEdges.
Under MOMENTUM_EQUATION DIFFUSIVE_WAVE a face firing evaluates the face depth exactly as the default closure does — (9-11) under MEAN, (9-12) under VFR_FACE, a wall at or below DRY_DEPTH with the discharge zeroed — applies the \(10^{-12}\) m deadband of §9.5.3 to the free-surface difference, and then replaces (9-7) with the quasi-steady Manning balance,
\[q_{f} = -\,\frac{h_{f}^{5/3}\,S}{n_{f}\,\sqrt{\max(\lvert S \rvert,\ S_{\varepsilon})}}, \qquad S = \frac{\eta_{R} - \eta_{L}}{d_{n}}, \qquad S_{\varepsilon} = \frac{\text{FLUX\_DH\_EPS}}{d_{n}}\] | (9-42) |
with \(n_f = \tfrac{1}{2}(n_L + n_R)\) the face roughness. The floor \(S_{\varepsilon}\) is what makes an explicit diffusive wave usable: without it the conductance \(h^{5/3}/(n\sqrt{\lvert S\rvert})\) diverges as the surface flattens, and the flux across a nearly level pond grows without bound. Below the floor the discharge is linear in \(S\), which keeps the C-property — the flux still vanishes with the slope — while bounding the transmissivity. FLUX_DH_EPS is a head difference (default 4 mm), so the floor scales with the face spacing. The conveyance factor, the positivity share \(\beta/n_v\) of (9-15) and the booking (9-13) are those of the default closure, and the face discharge is recomputed from scratch at every firing — nothing is remembered from the previous one.
Substituting (9-42) into (9-1) makes the scheme a nonlinear diffusion \(\partial h/\partial t = \nabla\cdot(K\nabla\eta)\), \(K = h^{5/3}/(n\sqrt{\lvert S\rvert})\), and the explicit step of such an equation is bounded by \(\Delta t \le L^{2}/(4K)\). Per cell,
\[\Delta t_{i} = \alpha\,\frac{L_{char,i}^{2}\,n_{i}\,\sqrt{S_{max,i}}}{4\,h_{i}^{5/3}}, \qquad S_{max,i} = \max\!\left( \max_{f}\lvert S_{f}\rvert,\ \frac{\text{FLUX\_DH\_EPS}}{L_{char,i}} \right)\] | (9-43) |
with \(\alpha =\) CFL_NUMBER reused as the safety fraction, \(L_{char}\) the operator length of (9-4), and \(S_{max}\) the largest free-surface slope over the cell's faces at the rebuild (walled faces count as zero); a dry cell is unbounded, and MAX_TIMESTEP caps the result. The slope census is refreshed at every rebuild and every \(\Delta t_0\) re-minimization. Because the bound scales as \(L^{2}\) while (9-17) scales as \(L\), a fine cell at a coupling point is punished far more severely under this closure, and the tier spread of §9.5.6 is what keeps the rest of the mesh from paying for it. No Froude clamp and no \(\theta\) blend are applied — there is no inertia to damp.
Boundary edges under this closure use the collapsed-Manning conductance of §9.6 for every type, and FRONT_REBUILD stays off by default because a diffusive front cannot outrun the fixed rebuild cadence (§9.5.7).
Implementation: src/engine/2d/solver/DiffusiveKernels.hpp — faceDischarge is (9-42) and cellDiffusiveDt is (9-43); the face firing is openswmm::twoD::ExplicitInertialSolver::fireFacesDiffusive and the slope census refreshDiffusiveSlopes in src/engine/2d/solver/ExplicitInertialSolver.cpp.
Boundary edges — those claimed by only one cell — default to no-flux walls. [2D_BOUNDARY_CONDITIONS] assigns any of five types per edge:
| Type | Parameter | Meaning |
|---|---|---|
WALL | — | Zero flux (default) |
NORMAL_FLOW | signed bed slope \(S\) | Manning flow \(q = \mathrm{sign}(S)\,h^{5/3}\sqrt{\lvert S\rvert}/n\) per metre of edge; a positive slope falls away from the domain and drains it, a negative one feeds it with the interior depth (it cannot wet a dry cell) |
SPECIFIED_STAGE / TS_STAGE | head, or time series | Prescribed free-surface elevation |
SPECIFIED_FLOW / TS_FLOW | discharge per metre, or time series | Prescribed unit discharge, outward positive |
RATING_CURVE | curve name | Stage → unit discharge lookup, resolved each step from the boundary cell's stage |
TRI is the cell index of §9.3.1 and EDGE the local edge — 0..2 for a triangle, 0..3 for a quad. PARAM_2 is reserved and is written as *; it must be present whenever GROUP is. GROUP is an optional label, the sixth token of the row (* for none). It is stored with the row and written back by the .inp writer unchanged, and the engine never acts on it: it exists so that the edges of one shoreline, tidal reach or outfall line can be addressed together by a user interface or a script — assigned, re-typed or deleted as a set — without the engine having an opinion about which edges belong together. A row that writes a group also writes its * placeholder, so the two are never confused.
Time series and curve names are resolved to registry indices once, on the first advance, and evaluated every routing step thereafter. Prescribed stages share the mesh's vertical datum and prescribed flows the project's flow units, so both are converted to SI on the same terms as the mesh itself (§9.7.4). Under FACE_RECONSTRUCTION VFR_FACE every type conveys with the wetted-edge depth (9-12) of the driving surface over the boundary edge's own endpoint beds, so a cell whose water pools away from an outlet edge does not leak through it.
How a boundary edge is integrated depends on the momentum closure.
Local inertial. A stage boundary is integrated with the interior momentum law, not with a conductance. The ghost state holds \(\eta = \eta_{bc}\) with a zero-gradient discharge, sitting across the edge at the centroid-to-edge distance \(d_e\) of (9-33) — \(2A/(3L)\) for a triangle — and (9-7) is applied to it exactly as to an interior face, followed by the Froude clamp (9-10). The earlier treatment — a collapsed Manning flux toward the prescribed stage — was a diffusive-wave law grafted onto an inertial interior, and it showed: its conductance saturated the equilibrium clamp into a Dirichlet cell, and every boundary-driven steady case floated one head jump, of order \(v^{2}/2g\), above the stage it had been given. The bump cases of §9.10 pass because of this change. A per-substep equilibrium clamp remains as a backstop: one substep may move a cell at most to the prescribed stage, never past it. At the inertial law's gravity-scale fluxes it rarely binds; on a tiny cell it prevents overshoot. The other three types apply their per-metre discharge directly — the Manning law above for NORMAL_FLOW, the prescribed value for SPECIFIED_FLOW and RATING_CURVE — and the edge's contribution is added to the cell's Perot vector (9-16), without which a cell fed through its boundary carried a systematic \((1-\theta)\) drag on every face.
Full shallow water. Every boundary edge is a ghost-cell Riemann problem: the ghost state per type is chosen in the outward-normal frame and (9-36)–(9-39) are evaluated against it, giving both a mass flux and a momentum change for the cell.
WALL mirrors the normal velocity; the mass flux is zero and the pressure balance is closed (§9.5.10).NORMAL_FLOW prescribes the Manning mass flux above, with the ghost velocity \(q/h_i\); a zero slope, depth or roughness conveys nothing and the ghost is a mirror — a bare zero-gradient ghost is an absorbing boundary and drained a lake at rest through an inert outlet.SPECIFIED_STAGE holds \(\eta_{bc}\). For subcritical flow one characteristic leaves the domain, so its Riemann invariant fixes the ghost velocity for the prescribed depth,\[u_{n,g} = u_{n} + 2\left( \sqrt{g h_{i}} - \sqrt{g h_{g}} \right), \qquad h_{g} = \max(0,\ \eta_{bc} - z_{i})\] | (9-44) |
rather than a bare copy of the interior velocity; under supercritical outflow ( \(u_n \ge \sqrt{g h_i}\)) no characteristic enters and the ghost is transmissive — the prescribed stage cannot act.
SPECIFIED_FLOW and RATING_CURVE are flux boundaries: the mass flux is the prescribed value verbatim, because the Riemann solve's own mass flux is a wave-speed-weighted blend of the interior and prescribed states that under-delivered the SWASHES bump inflows by 7–21 % and reversed a supercritical inlet outright. The ghost still sets the momentum flux and keeps the Audusse pressure balance, so its depth is the physical boundary depth: the critical depth \((q^{2}/g)^{1/3}\) on a dry cell (dry-bed inflow), the interior depth for outflow or supercritical inflow, and for subcritical inflow the depth that the outgoing invariant \(r = u_n + 2\sqrt{g h_i}\) and the discharge together fix,\[2\sqrt{g}\,s^{3} - r\,s^{2} + q = 0, \qquad s = \sqrt{h_{g}}\] | (9-45) |
solved by a safeguarded Newton iteration from the interior depth and falling back to it when no physical root converges — a fallback costs accuracy in the momentum flux, never conservation.
No equilibrium clamp is applied under this closure — the Riemann flux is bounded by the two states — and when the availability clamp below shrinks the requested mass, the momentum already booked is rescaled in the same ratio.
Diffusive wave. Every type is served by the collapsed-Manning conductance the local-inertial stage boundary retired: for SPECIFIED_STAGE the per-metre flux toward the prescribed stage is \(-\,h_{up}^{5/3}\,\mathrm{sign}(\Delta\eta)\,\rho(\lvert\Delta\eta\rvert)/(n\sqrt{d_e})\) with \(h_{up}\) the upwind depth and \(\rho\) the regularized square root of (9-27) with FLUX_DH_EPS in place of \(\varepsilon_o\) — the same \(C^{1}\) quadratic below the floor, the bare root above it — so the transmissivity stays bounded as the two surfaces level; the other three types apply their per-metre discharge as above. Under the default closure FLUX_DH_EPS is not consulted at all: the inertial stage law replaced the only place it acted.
Whatever the closure, exchange is clamped in volume space and the booked flux re-derived from the applied change, so what is reported is exactly what was applied: the provisional volume is floored at zero, and under the local-inertial and diffusive laws a stage boundary is also held at the prescribed stage by the equilibrium clamp. Cumulative boundary volume is tracked per edge, outward positive, and enters the 2D mass balance of §9.9.
Coupling is bidirectional and mass-conservative, and it is where a 2D module earns or loses its credibility. Two mechanisms exist, one for junctions and one for outfalls.
[2D_VERTEX_NODE_MAP] and [2D_TRIANGLE_NODE_MAP] associate a mesh vertex or cell with a SWMM node, with a discharge coefficient \(C_d\) (default 0.65) and an exchange area. With the legacy COUPLING_LAW ORIFICE, exchange is an orifice law on the head difference:
\[Q = C_{d}\,A_{eff}\,\mathrm{sign}(\Delta h)\,\sqrt{2g}\ \varphi\!\left(\lvert\Delta h\rvert\right), \qquad \Delta h = h_{2D} - h_{1D}\] | (9-18) |
positive draining the surface into the network. Three regularizations turn (9-18) from a stiffness source into something an explicit solver can integrate:
A bounded square root. \(dQ/d\Delta h \to \infty\) as \(\Delta h \to 0\) is exactly the regime a fill-and-spill manhole hovers in. Below 2 cm, \(\varphi\) is a \(C^1\) quadratic matching \(\sqrt{x}\) in value and slope at the join and having finite slope at zero.
A capped-pipe gate. A manhole with its lid on exchanges through the network only when the higher of the two heads reaches the crown elevation \(z_{inv} + D_{full}\). A Hermite smoothstep over a 5 cm band above the crown opens the exchange, and the effective area transitions smoothly from the inlet area to twice it over the same scale as the node surcharges.
Source-side wet/dry ramps. \(Q\) is multiplied by a smoothstep on the source side's depth relative to DRY_DEPTH, so a drain self-limits to zero as the cell empties and a spill self-limits as the node empties. This replaces a held-flux availability cap and is what makes the exchange stable inside the solver's inner loop rather than only across a window.
The exchange is evaluated live, at tier-0 cadence, against the current 2D heads and the routing step's 1D heads, and \(\int Q\,dt\) is accumulated exactly per point. Two hard caps make the ledger unfalsifiable: a drain may take at most \(\beta\) of the source cell's volume per substep, and a spill draws against a per-node budget of the node's stored volume for the whole advance — so the same water cannot spill twice within a routing step.
The 2D head at a coupling point is not simply a cell head. Vertex-coupled points use a wet-masked, depth-weighted mean of the incident cells' free surfaces under the VFR closure: a manhole vertex is commonly carved below the surrounding terrain, and a geometric average over incident cells would read dry-cell bed elevations as a water surface and pin the exchange at a phantom head from the first step.
Under dynamic wave routing the exchange would otherwise be a zero-sensitivity explicit source in the node continuity equation, which churns the Picard iteration. The head sensitivity
\[G = -\frac{\partial Q}{\partial h_{1D}} = C_{d}A_{eff}\sqrt{2g}\ \varphi'\!\left(\lvert\Delta h\rvert\right)\cdot(\text{gate})\cdot(\text{ramp}) \ \ge\ 0\] | (9-19) |
is scattered into the node's \(\sum dQ/dH\) denominator each iteration. The gate and ramp derivatives are deliberately dropped so the term can only be positive — a pure damping contribution, never a destabilizing one.
Exchange volumes reach the 1D side through the lateral-inflow delivery queue, drained at a uniform rate over the batch span rather than as a single-step pulse.
COUPLING_LAW WEIR_ORIFICE combines free and submerged weir control with orifice control. Heads and geometry are evaluated in metres; Q > 0 denotes surface-to-network flow. With crest z, surface head h2, network head h1, surface depth hs = max(h2-z,0), crest perimeter L and opening area A:
(Qw^-16 + Qo^-16)^(-1/16). It is no greater than either control and is about 4.2% below either when the two are equal.The free/submerged transition is smoothed over a submergence ratio of 0.1. The regularised root phi equals sqrt(x) above its head threshold and is quadratic below, with a finite derivative at zero. The routed implementation uses a 0.02 m threshold. That is a numerical regularisation scale, and is large relative to the heads of the Rubinato laboratory rig. A comparison using a smaller threshold only in the isolated law test does not validate the routed model at that scale.
The dynamic-wave head sensitivity is a nonnegative damping approximation: blend-weight and availability-ramp derivatives are dropped. It should not be described as the exact derivative of the blended expression.
| Preset | Cw | Csw | Coi | Coo |
|---|---|---|---|---|
| DEFAULT | 0.540 | 0.360 | 0.167 | 0.160 |
| RUBINATO2017 | 0.540 | 0.056 | 0.167 | 0.160 |
| RUBINATO2017_LOWER | 0.493 | 0.055 | 0.160 | 0.159 |
| RUBINATO2017_UPPER | 0.587 | 0.057 | 0.174 | 0.161 |
These are implementation presets, not four independently calibrated physical regimes. Rubinato et al. calibrated free drainage (S1), submerged weir drainage (S2), wet outward orifice flow (S3, expected 0.167) and dry outward flow (S4, expected 0.160). The implementation uses Coo = 0.160 for both outward cases. Coi is borrowed by analogy; inward orifice flow was not calibrated in that paper. The default Csw = (2/3) Cw is a continuity choice, not the paper's submerged-weir calibration. Per-node coefficients override the selected preset; authored mapping-row Cd/Area can override the orifice terms and opening geometry.
COUPLING_HEAD CELL, VERTEX and STENCIL_WET_MEAN select different head representations. The wet stencil excludes dry bed elevations from a water surface estimate; it is not the same as averaging all numerical cells within a circular opening, as in the paper's HMho construction.
[2D_NODE_COUPLING] describes node eligibility and structure geometry before mesh mapping. NONE excludes a node; MANHOLE, GRATE, CURB, OPEN and OUTFALL distinguish the applicable geometry or boundary treatment. Geometry is authored in project length/area units and converted at routing.
For non-storage coupled junctions, COUPLING_POND_MODE FOOTPRINT, CHAMBER and TRACKING select the effective above-rim storage. The smoothed area transition lies above the rim. Storage nodes retain their storage curves. COUPLING_EQUALIZE bounds the 1D-side exchange over the interval for which its head is held; the surface state evolves within that interval. It is not a proof that a frozen two-reservoir equilibrium is never crossed. Availability budgets still prevent spending donor water twice.
A coupling configuration has three independent inputs: [2D_NODE_COUPLING] selects the eligible structure and optional node overrides; a vertex/cell map selects its position on the surface; [2D_OPTIONS] supplies global law and numerical defaults. A node setting alone creates no hydraulic exchange. A NONE node is excluded even if an old map still names it; initialization warns about ignored rows. An explicitly eligible node without a map also produces a warning. Node overrides persist independently of mesh regeneration.
| Node class | Default without an authored row | Allowed explicit structures |
|---|---|---|
| Ordinary junction, divider | MANHOLE | NONE, MANHOLE, GRATE, CURB |
| Inlet or virtual junction | NONE | NONE, MANHOLE, GRATE, CURB (explicit enabling requires a deliberate hydraulic representation) |
| Storage | OPEN | NONE, MANHOLE, GRATE, CURB, OPEN |
| Outfall | OUTFALL | NONE, OUTFALL |
An inlet junction already receives street drainage through its inlet representation; its default avoids adding a second surface route. A virtual junction is a conduit split point without a physical rim. OPEN uses storage geometry; OUTFALL follows the link/tailwater exchange described in §9.7.2, not the junction weir/orifice law.
Each row is NodeID TYPE KEY=value .... Omitted parameters inherit the defaults below. Clearing an override differs from entering zero. Node lengths are in project length units (m for SI flow units; ft for US flow units) and node areas in their square. An SI mesh tag (;; UNITS: SI (m)) does not convert these node inputs to metres on a US project. By contrast, map-row areas follow the mesh units. Coefficients and lid fractions are dimensionless in the exchange equations used in this chapter; tailwater relaxation is in seconds.
| Parameter | Interpretation, range and default |
|---|---|
D | Positive opening diameter; default 1.2 m or 4 ft. A manhole starts with \(A=\pi D^2/4\), \(L=\pi D\) |
AREA | Positive hydraulic opening area. For MANHOLE, replaces the diameter area; for GRATE/CURB, supplies opening area (fallback: authored map area, then diameter area). For OPEN, overrides storage surface area at full depth |
PERIMETER | Positive effective weir length. Overrides derived length. Without it, GRATE/CURB use \(4\sqrt A\), OPEN uses \(\sqrt{4\pi A}\), and MANHOLE uses \(\pi D\) |
CHAMBER_AREA | Positive chamber plan area. Under WEIR_ORIFICE, defaults to resolved opening area before lid scaling; under ORIFICE, defaults to the map area at structure resolution. Used by chamber/tracking storage treatment; storage-node curves remain authoritative |
LID_OPEN_FRACTION | In (0,1], default 1. Scales resolved opening area and weir length, not chamber area. A zero value is invalid; use NONE for no exchange |
CREST | NODE_RIM default: invert + full depth. MESH_BED: mapped vertex bed or cell-centroid bed. MAX or MIN: maximum/minimum of those elevations |
LAW | ORIFICE or WEIR_ORIFICE; absent inherits global COUPLING_LAW. Orifice mode retains map/global Cd and area; the weir geometry/crest/preset inputs do not redefine it |
COEFFS | DEFAULT, RUBINATO2017, RUBINATO2017_LOWER, RUBINATO2017_UPPER; absent inherits the global preset, subject to the OPEN rule below |
CW, CSW | Positive overrides for free-weir and submerged-weir drainage coefficients |
COI, COO | Positive overrides for surface-to-sewer and sewer-to-surface orifice coefficients |
SPREAD | Outfalls only: nonnegative discharge radius, default 0. Expands receivers to connected cells within the radius that are wet or below/equal to the highest original receiving bed; discharge is area weighted. Withdrawal uses original donor reservations |
TAILWATER_RELAX | Outfalls only: nonnegative relaxation time, default 0. Cached head updates by \(h_{tw}\leftarrow h_{tw}+\min(1,\Delta t/\tau)(h_{2D}-h_{tw})\) for positive τ; zero applies the current surface head immediately. First initialization uses the surface head |
For WEIR_ORIFICE, coefficient resolution uses the following values before individual overrides:
| Preset | CW | CSW | COI | COO |
|---|---|---|---|---|
DEFAULT | 0.540 | 2CW/3 = 0.360 | 0.167 | 0.160 |
RUBINATO2017 | 0.540 | 0.056 | 0.167 | 0.160 |
RUBINATO2017_LOWER | 0.493 | 0.055 | 0.160 | 0.159 |
RUBINATO2017_UPPER | 0.587 | 0.057 | 0.174 | 0.161 |
OPEN with no node COEFFS override | 1 | 2/3 | 1 | 1 |
OPEN applies its special coefficients even when the global preset differs; a node COEFFS override selects a preset explicitly. Individual CW/CSW/COI/COO then override the result. With DEFAULT or the special OPEN rule, CSW is computed as 2CW/3 after any CW override, unless CSW is explicitly supplied. The Rubinato ranges are laboratory coefficient alternatives, not confidence bounds on predictions at another site. The low experimental CSW is not the continuous-join DEFAULT value; see the comparison's limitations.
Finally, an authored map-row Cd overrides both COI and COO, and an authored map-row area overrides the resolved opening area. The perimeter is derived before this final area override; changing only map area does not resize an independently derived/authored perimeter. Chamber area defaults from the final opening area before lid scaling unless explicitly supplied. Thus a map inherited from a legacy model can mask changes to a node's coefficients or opening area. COUPLING_AREA AUTO applies only to ORIFICE: it derives 1.25 times the largest connected-conduit area, bounded to 0.05–2.0 m², where a positive conduit area is available. WEIR_ORIFICE retains its structure geometry.
In SWMMVis, set 2D Coupleable, the structure and its parameters in node properties or the node attribute table. NO authors NONE; YES restores the class default (or explicit MANHOLE when that default is NONE). Switching structure/eligibility clears prior structure parameters; Undo restores the full state. Blank editors restore inheritance; the table displays NULL for an unauthored numeric override. Node-table booleans display YES/NO and support single-value/selected-row paste; spreadsheet blocks can supply multiple geometry values. Set type before pasting parameters.
During mesh generation, 1D geometry influence controls Steiner points and constraints, while Map model nodes to the mesh after generation independently runs spatial registration. Map only coupleable nodes normally stays on; Set coupleable… edits the candidate set. Coincident nodes map to vertices; other in-domain nodes map to containing cells. Nodes not pinned because of minimum separation remain available for cell mapping. Outside nodes are reported.
Mesh 2D → Remap 1D↔2D previews Current/Proposed mappings and |rim − bed|. Preserve existing retains eligible mappings and fills missing ones; Clear & rebuild remaps by coordinate, retaining Cd/Area for a node returning to the same location. Apply records undoable edits. Use this after remeshing or changing coordinates; a node attribute edit alone does not relocate its map. The mapper leaves new Cd/Area fields unauthored when the global law is WEIR_ORIFICE; per-node law changes alone do not remove existing map overrides.
Neither a selected structure nor a fine mesh guarantees experimental accuracy. Verify opening dimensions and rim/terrain datum, inspect map overrides, compare routing timesteps and check continuity and exchange hydrographs. Head-sampling stencils and delivery cells are numerical supports; they do not automatically represent a spatial average over the physical opening footprint.
The accompanying Rubinato comparison compares all 36 steady measurements with the bare paper equations, a junction reconstruction and an actual-shaft-storage sensitivity. It preserves the fit, continuity, head residuals and successive-window drift separately. The paper's printed regression slopes have nonzero intercepts; matching the slopes does not establish absolute flow agreement.
The nine transient demonstration shapes use Table 1 peak ratios and durations. Their boundary forcing is synthetic and must not be labelled validation against the measured events or Figure 11. The 0.0125 m local mesh check demonstrates that the equalisation cap does not choke exchange; it does not establish mesh convergence or footprint-resolved delivery. The current comparison is a qualified numerical example. It does not support an experimental validation claim.
Reference: Rubinato et al. (2017), Journal of Hydrology 552, 421-432, doi:10.1016/j.jhydrol.2017.06.024.
Rubinato 2017: Sewer-Surface Exchange Comparison — measured steady comparisons, synthetic events, report and example downloads.
An outfall coupled to the mesh works in both directions.
Before each 1D routing step, surface cells reserve water for their coupled outfalls. If \(m_i\) outfall footprints share cell \(i\), each receives \(\max(V_i,0)/m_i\). The reservations of an outfall are summed to \(B_o\) and divided among its \(d_o\) incident hydraulic links. Each link's surface-to-pipe discharge is bounded by \(B_o/(d_o\Delta t)\) (converted from m³/s to ft³/s). The dynamic-wave solver applies this bound before node continuity; the finite-volume solver applies it to the boundary mass flux before either side integrates it. Positive pipe-to-surface discharge is not restricted by this donor bound.
After routing, the accepted net volume is transferred once. A withdrawal uses the same per-cell reservations, proportional to their available water; a discharge is distributed by cell area. The cell volumes and the outfall ledger receive identical increments. The transfer precedes surface routing, infiltration, evaporation and live junction exchange, so those processes cannot spend water already withdrawn by an outfall. Transported surface mass is removed in the same fraction as the donated water. Every routing step is booked in its actual direction, including when COUPLING_SYNC batches several surface steps. Reversals within a batch therefore retain both the incoming and outgoing totals; the current stored-volume total includes each transfer.
This is a conservative, first-order split exchange. Equal reservation shares can restrict one active outlet while another leaves its share unused; unused water remains on the surface and becomes available in the next routing step. Time-step refinement is still needed to assess exchange accuracy. Availability limiting is not a substitute for that accuracy check. Previously the code capped a stencil-wide request and then scattered a held sink with different weights; individual cells could reject that sink while the ledger booked the full request. That is not a conservative transfer.
Inward, the 2D surface acts as dynamic tailwater. The 2D stage at the coupling point is cached, and the outfall's boundary condition becomes \(\max(h_{standard}, h_{2D})\), applied inside the dynamic-wave iteration so it survives every Picard pass. Flap gates are honoured: the gate decision is made where the current \(h_{standard}\) is visible.
The wet/dry gate here is keyed on depth in excess of DRY_DEPTH, not on depth. The reason is specific: a draining cell comes to rest at a film at or just below DRY_DEPTH, which the solver treats as immovable. A ramp keyed on depth alone would read ≈ 1 at that resting film and pin the outfall at a tailwater it can never drain below — a deadlock in which the pipe cannot discharge and the cell cannot dry.
By default the two domains co-advance every routing step: the 2D solver advances over exactly \([t, t+\Delta t]\), and exchange volumes reach the 1D side with at most one routing step of lag. This keeps fill-and-spill coupling free of the batch-delay ringing that a longer exchange interval produces at weir and culvert ponds.
COUPLING_SYNC batches the 2D advance over a longer span (clamped to between one routing step and 60 s). It is a wall-clock lever for large meshes, where per-routing-step advances degenerate into the tail handling of §9.5.6, and it should be understood as trading accuracy for speed: the held-exchange error grows with the span.
The 2D solver runs internally in SI — metres, m³, m³/s, \(g = 9.80665\). The 1D engine always computes internally in feet, for every project, US or metric: its reader converts metric input to feet on load and converts back only at the display boundary. The 1D↔2D coupling factors are therefore always the feet–metres conversion, independent of FLOW_UNITS.
The mesh is different: it is authored in the project's display length units, so it is scaled to SI on load for US projects and left alone for metric ones. A mesh file may declare ;; UNITS: SI (m) to assert it is already metric regardless of FLOW_UNITS. Boundary stages share the mesh's datum and scale with it; boundary flows scale with FLOW_UNITS.
Getting this wrong is silent and severe — an earlier version tied the coupling factors to FLOW_UNITS, which collapsed them to 1.0 on metric projects and left every coupled head off by 3.28× and every exchanged volume off by 35×.
Parameters. Each coupling point carries the SWMM node index, the discharge coefficient \(C_d\) (default COUPLING_CD, 0.65), the exchange area \(A\) and, for outfalls, the flap-gate flag. A vertex row that authors no area defaults to 1.0 in mesh area units (scaled to m² with the mesh); with COUPLING_AREA AUTO, unauthored areas are derived at resolve time as \(\mathrm{clamp}(1.25 \times A_{conduit,max},\ 0.05,\ 2.0)\) m², where \(A_{conduit,max}\) is the full-flow area of the largest conduit connected to the node. A vertex-coupled point records the first triangle incident on its vertex as its host cell; that cell is pinned to tier 0 (§9.5.6).
The regularized exchange law. The pieces of (9-18), with all heads in the 2D metre frame ( \(h_{1D} = 0.3048 \times\) the node head, crown \(z_{cr} = 0.3048 \times (z_{inv} + D_{full})\), \(h_{max} = \max(h_{1D}, h_{2D})\)) and \(\varepsilon_o = 0.02\) m:
\[\varphi(x) = \sqrt{x}\] | \(x \ge \varepsilon_{o}\) | (9-27a) | |
\[\varphi(x) = \frac{3\,x}{2\sqrt{\varepsilon_{o}}} - \frac{x^{2}}{2\,\varepsilon_{o}^{3/2}}\] | \(0 \le x < \varepsilon_{o}\) | (9-27b) |
which matches \(\sqrt{x}\) in value and slope at \(\varepsilon_o\) and has the finite slope \(3/(2\sqrt{\varepsilon_o})\) at zero. The effective area grows linearly from the inlet area at the crown to twice it 5 cm above,
\[A_{eff} = A\left\lbrack 1 + \min\!\left( 1,\ \frac{h_{max} - z_{cr}}{0.05} \right) \right\rbrack\] | \(h_{max} \ge z_{cr}\); else \(A_{eff} = A\) | (9-28) |
and the capped-pipe gate and the wet/dry ramps multiply the orifice flow by cubic Hermite smoothsteps over the same scales,
\[Q \leftarrow Q\,\sigma(c), \qquad c = \mathrm{clamp}\!\left( \frac{h_{max} - z_{cr}}{0.05},\ 0,\ 1 \right)\] | gate | (9-29a) | |
\[Q \leftarrow Q\,\sigma\!\left( \mathrm{clamp}\!\left( d_{src}/h_{dry},\ 0,\ 1 \right) \right)\] | wet/dry ramp | (9-29b) |
with \(\sigma(t) = t^2(3 - 2t)\) and \(h_{dry} =\) DRY_DEPTH. The gate reads exactly zero at the crown, so a capped node exchanges nothing until one side surcharges past it. For a drain ( \(Q > 0\)) the source depth \(d_{src}\) is the maximum depth over the vertex stencil (or the single cell for centroid coupling); for a spill it is the node depth, converted. The driving head \(h_{2D}\) is the wet-masked, depth-weighted stencil mean under the VFR closure, the pseudo-Laplacian vertex head under FLAT, and the cell head for centroid coupling.
Caps. Inside the marcher, at each tier-0 substep of length \(\Delta t_0\):
\[Q_{drain} \le \frac{\beta \max(V_{c}, 0)}{\Delta t_{0}}\] | drain | (9-30a) | |
\[\lvert Q_{spill} \rvert\,\Delta t_{0} \le V_{node}\,f_{v} - D_{node}\] | spill | (9-30b) |
where \(\beta = 0.8\), \(V_c\) is the live coupling-cell volume, \(f_v\) the ft³-to-m³ factor and \(D_{node}\) the volume already drawn from that node during the current advance — the same stored water cannot spill twice within a routing step. The capped \(Q\,\Delta t_0\) is applied directly to the coupling cell's volume (floored at zero) and accumulated into the point's ledger \(\int Q\,dt\).
One routing step, in order:
DRY_DEPTH (§9.7.2). The outfall boundary logic applies \(\max(h_{standard}, h_{2D})\) inside every dynamic-wave iteration, blending by the cached factor and honouring flap gates. Reserve donor volumes and publish the per-link withdrawal limits described in §9.7.2.COUPLING_SYNC clamps to between one routing step and 60 s) the co-advance batch fires:coupling_flux sources or rebook their already applied ledger.Forcing override surface. The C API can replace or augment the computed exchange per cell: swmm_2d_force_coupling_flux(engine, idx,
value, mode, persist) prescribes a coupling rate (m/s, positive into the 2D domain) with mode selecting override or add and persist selecting one-shot or persistent application; swmm_2d_force_rainfall / swmm_2d_force_evap (and their _uniform variants) do the same for the meteorological sources, and swmm_2d_force_clear_all clears every prescription. One-shot prescriptions expire after the step they apply to; any change marks the forcing state dirty so it takes effect on the very next batch regardless of the 30 s refresh cadence.
Implementation: the exchange law and its sensitivity are openswmm::twoD::computeNodeCouplingQ and openswmm::twoD::computeNodeCouplingDQdh1d (src/engine/2d/coupling/NodeCoupling.cpp, which also holds the file-local orificePhi, effectiveArea, wetVertexEta, scatterCouplingFlux, budgetAvail/budgetDraw/budgetCredit helpers); outfall accumulation and tailwater caching are OutfallExchange and openswmm::twoD::updateOutfallBoundaries; the in-marcher exchange loop is the tier-0 tail of openswmm::twoD::ExplicitInertialSolver::fireCells; the batch orchestration and booking are openswmm::twoD::SurfaceRouter2D::advancePostRouting, openswmm::twoD::SurfaceRouter2D::coAdvanceStep and openswmm::twoD::SurfaceRouter2D::accumulateMassBalance; the forcing entry points are declared in include/openswmm/engine/openswmm_2d.h.
COUPLING_IN_FLOODING)The 1D routing continuity table of the status report carries a 2D Coupling Outflow row for the water that leaves a coupled junction onto the mesh, separate from Flooding Loss. That split is the accurate reporting: a coupling transfer is not flooding — the water is still in the model, on the surface, and returns through the drains and the coupled outfalls. COUPLING_IN_FLOODING YES books the spill into the Flooding Loss row instead, and exists for one reason: comparing a run against a build from before the split, where the spill sat inside Flooding Loss.
The key changes the report only. The spill is an outflow of the 1D system under either grouping, so the routing continuity error is identical; every routed volume is identical; and the per-step flooding rate that drives the flood reporting and the node flooding statistics keeps the spill in both cases, since those statistics are about water leaving the node upward, which a spill is whatever row the continuity table puts it in. The 2D → 1D drain is external inflow under either setting.
Figure 9-6 Where a node spill is booked under COUPLING_IN_FLOODING: the same totals and the same continuity error, in two different rows
Rainfall reaches the mesh from the project's rain gages, mapped by RAINFALL_MODE: NATURAL_NEIGHBOUR (the default) interpolates the located gages onto every cell centroid with natural-neighbour (Laplace) weights inside the gages' convex hull and inverse-distance weights outside it; NEAREST_NEIGHBOUR gives every cell the reading of its closest located gage; SYSTEM applies the arithmetic mean of all gages uniformly and is the automatic fallback when no gage has a map location; NONE applies no rain to the mesh. The interpolation — its weights, its Delaunay construction, its degenerate cases and the way the weights are built once and applied each step — is hydrology, and Hydrology Chapter 8 gives it. NONE is not an optimization, it is a modelling decision. If the project's subcatchments already convert the storm to runoff and deliver it to nodes, rain on the mesh double-counts the same storm. Use NONE whenever the surface is meant to receive only what the network gives it.
Evaporation on the mesh is governed by EVAPORATION:
YES (the default) makes evaporation a forcing-only sink: the demand rate is zero unless the C API prescribes one per cell or uniformly (swmm_2d_force_evap, §9.7.5). The project's [EVAPORATION] section does not reach the mesh.CLIMATE takes the project's climate evaporation rate — the same rate the subcatchments see, converted to m/s — as the demand on every cell not carrying a forcing override; an override replaces or adds to it per cell as the forcing API directs.NO holds the sink at zero and ignores forcing.Whatever supplies the demand, the sink actually applied is the demand tapered by the cubic Hermite ramp (9-26) below DRY_DEPTH, with \(e_{eff} = 0\) for a dry cell or a negative demand: a drying cell cannot evaporate more water than it holds, and negative demand is not a condensation source. The taper is applied in the same form at every point the sink enters — the lazy pass over inactive cells and the cell firing of §9.5.8 — and evaporation is taken first, with infiltration applied against what is left, so a surface-only project is bitwise unchanged by the infiltration machinery; both sinks feed a volume update floored at zero, so the pair cannot leave a cell with negative water.
Infiltration on the mesh — the per-cell methods, the destination of the infiltrated water (lost, the subcatchment aquifer, or the mesh aquifer) and its own step — is hydrology and is documented in Hydrology Chapter 8; the mesh aquifer it can feed is Hydrology Chapter 9. With INFILTRATION NO, or no infiltration rows, water on the surface leaves only by flowing away, evaporating, or entering the network.
The 2D domain carries its own mass balance, printed as a separate block of the status report because its sign conventions are opposite to the 1D routing balance:
Volumes are SI. Every term is booked from the volume actually applied — after every cap, clamp and rescale — so the balance reflects the volumes the solver actually applied. A mesh that infiltrates, or that carries an aquifer, gains the rows those processes book; they are described with the processes in Hydrology Chapter 8 and Hydrology Chapter 9. Where the spill appears in the 1D table is the subject of §9.7.6.
A 2D Solver Statistics block reports cumulative substeps, face-kernel evaluations, mean and last internal step, the minimum, mean and maximum active-cell fraction over the rebuild samples, and the occupancy share of each local-time-stepping tier. Those last two are the diagnostic for the two performance mechanisms of §9.5.6 and §9.5.7: a mesh with 100 % active cells is telling you the active set is not helping, and a tier histogram concentrated in tier 0 is telling you the same about local time stepping.
Per-cell results — depth, free-surface elevation, velocity, gradients, maximum-depth and maximum-velocity envelopes, cumulative volume and a per-cell continuity residual — are refreshed on a report-scale cadence rather than every routing step, because the six full-mesh passes cost several times the solver's own advance on a large mesh, and nothing consumes them faster than the reporting interval. Coupling and outfall heads read the solver's live state directly, never these derived fields. REPORT_2D NO skips the derived-gradient refresh; it does not suppress the results file. The reported gradients are Green–Gauss gradients of the free surface limited by the Jawahar–Kamath (2000) limiter, whose regularization is LIMITER_EPSILON; the limiter affects the output fields only and enters no solver update.
FULL_SWE cell velocities come from the accepted conserved cell momentum, including the initial snapshot and runs with all optional momentum processes disabled. Local-inertial and diffusive-wave cell velocities are output-only area-weighted least-squares fits to wet edge-normal velocities. For edge volume flux Q and hydraulic cross-sectional area A, the fit minimizes sum A (n dot u - Q/A)^2. A uses the selected MEAN/VFR face depth, edge length and interior transmissivity; stage boundaries use their own closure depths. Accumulating sum A n n^T and sum n Q avoids division by tiny edge sections and by the whole-cell mean depth. There is no diagnostic discharge clip or speed cap. The fit is a representative vector and need not reproduce every edge flux exactly; original edge fluxes remain the conserved hydraulic data.
Dry-cell vectors are zero. Wet-cell vectors are paired NaNs when wet edge normals cannot determine two components (relative determinant tolerance 1e-10 * trace^2), input is nonfinite, or an accepted averaged flow has no currently wet section. NaN means unavailable, not zero velocity. Velocity maxima use finite samples only, retain NaN if none exist, and cannot bound speeds during unavailable intervals. Depth and speed envelopes are sampled independently; their peaks need not occur at the same time. Good water continuity alone does not establish velocity accuracy.
Saved rendering uses available stored velocity arrays as authoritative, including their NaNs. Edge reconstruction is retained only for sources or frames without a stored pair. The bulk C API swmm_2d_get_velocities_bulk exposes the same physical Cartesian SI vectors and NaN convention. The GUI streams all three closures' vectors before the matching depth packet, including them in bounded history and sparse regional reads. Arrows and particles omit unavailable vectors and never borrow a neighboring vector into an unavailable cell. These output rules change no solver flux, momentum, wetting threshold, depth or water balance.
Two distinct vertex reconstructions exist and must not be confused. The solver field is the pseudo-Laplacian stencil of Kumar et al. (2009), built once from cell-centre geometry with Lagrange multipliers enforcing linear exactness; dry cells contribute their bed elevations, which the solver relies on. The rendering field is a wet-masked, depth-weighted mean of the incident wet cells' free surfaces, with a wetted-contact gate so that a cell votes at a corner only where its water actually reaches it. Interpolating the solver field for display would drag water surfaces up dry banks and down into thin films; interpolating the rendering field in the solver would break the active-set logic.
The results file. With OUTPUT_FILE set in [2D_OPTIONS], results are written to an HDF5 file following the CF-1.11 and UGRID-1.0 conventions for unstructured meshes — mesh topology (mixed topology for a mesh with quads, §9.3.1), node and face coordinates, bed elevations and roughness written once, then the selected time-varying fields appended. The file opens directly in ParaView, QGIS, or any CF/UGRID-aware reader. Four keys shape it:
OUTPUT_PRECISION — FLOAT32 (the default, halving the file at about seven significant digits) or FLOAT64 for bit-for-bit parity tooling. Mesh geometry always stays float64: projected coordinates are of order \(10^{6}\) m and single precision would lose centimetres there, whereas depth and velocity fields carry no such offset.OUTPUT_COMPRESSION — a deflate level 0–9 applied, with the byte-shuffle filter, to the chunked time-varying datasets; 0 writes them uncompressed.REPORT_2D_STEP — the cadence of the time-varying datasets, in HH:MM:SS or numeric seconds. Fractional values survive API queries and INP save/reopen; numeric seconds also apply to INFIL_STEP. It must be a positive multiple of REPORT_STEP; unset, the 2D file follows the report step.REPORT_2D_VARIABLES — which dataset groups are written. The value is a preset or a list of the tokens of Table 9-4; DEPTH is always included, since every reader keys the time axis on it, and a group not selected is not created in the file, so readers treat every time-varying dataset as optional. DEFAULT is everything a user renders or plots with the solver diagnostics off; MINIMAL is the depth map, the render reconstruction and the envelopes; ALL is every group; NONE is DEPTH alone.| Token | Datasets | In DEFAULT |
|---|---|---|
DEPTH | Mesh2_face_depth, Mesh2_face_head | always |
VELOCITY | Mesh2_face_vx, Mesh2_face_vy | yes |
EDGE_FLUX | Mesh2_edge_flux (GUI velocity reconstruction, profile flux) | yes |
NODE_HEAD | Mesh2_node_head, Mesh2_node_depth (the rendering reconstruction) | yes |
SPECIES | Mesh2_face_species_conc, when transport rows exist | yes |
RAINFALL | Mesh2_face_rainfall, Mesh2_face_rain_cum | yes |
INFILTRATION | Mesh2_face_infil_rate, Mesh2_face_infil_cum | yes |
COUPLING | Mesh2_face_coupling_flux, Mesh2_face_net_source | no |
GRADIENTS | Mesh2_face_grad_hx, Mesh2_face_grad_hy and their limited pair (solver diagnostics) | no |
CONTINUITY | Mesh2_face_continuity_err (solver diagnostic) | no |
ENVELOPES | Mesh2_face_max_depth, Mesh2_face_max_velocity, Mesh2_face_max_continuity_err | yes |
BUILDUP | Mesh2_face_buildup, when [2D_COVERAGES] rows resolved | yes |
GROUNDWATER | the per-cell Mesh2_face_gw_* fields, the groundwater ledger and the node-exchange cumulative, when a mesh aquifer resolved | yes |
GW_DETAILED | Mesh2_face_gw_theta_sigma per layer and face — \(m\) layers × nFace per step | no |
HEAT | Water, mesh bed, soil, and mesh-connected LID temperature/energy/capacity fields; fixed-grid soil layers; signed cumulative thermal receipts and residuals | no |
Table 9-4 The fifteen REPORT_2D_VARIABLES tokens
HEAT is opt-in and writes fields only when the corresponding thermal owner is initialized. Face fields use the [time, face] layout; soil-layer fields use [time, layer, face], with the top layer first and a static center-depth field in metres. Mesh-connected LID physical cells are reported on a separate [time, lid_layer] axis. The static Mesh2_lid_layer_* mapping datasets identify each LID layer's network node, mesh face, physical cell, control, authored layer, elevation bounds, and area. A zero-capacity physical layer has zero energy and capacity and a NaN temperature; this is distinct from an unowned mesh face, which is NaN in face-based fields. Mesh2_lid_layer_ground_contact_energy is cumulative signed energy in joules, positive into that LID layer from its mesh-ground contact. It is the equal-and-opposite per-layer counterpart to ground_heat_budget's lid_contact_j term. mesh_heat_budget and ground_heat_budget use [time, face, term], units of joules, and carry named term/sign-convention attributes. Their receipts are cumulative from initialization; internal water/bed and water/channel contacts are signed and cancel in the combined owner balance. The groundwater budget residual is stored energy minus initial energy and the listed cumulative receipts.
When BED_HEAT_MODEL LOCAL_STORAGE is active, HEAT writes the local bed's initial and current sensible energy, heat capacity, temperature, and fixed porewater volume. With transported species it also writes porewater concentrations, initial species amounts, cumulative surface-to-bed species exchange, and generated bed age. Species dimensions follow each dataset's species_names and species_units attributes; the temperature row is excluded, age concentration is in hours, amount fields multiply the species units by m³, and exchange is positive from surface water into the bed. Initial energy and initial species amounts are constant since-initialization fields; cells without local bed porewater are NaN. These fields use the model's local finite-bed storage; they do not include groundwater-owned channel contacts. The thermal result schema version is 6. HEAT is excluded from DEFAULT to preserve existing results-file contents; add it explicitly or select ALL to write it.
REPORT_2D_SPECIES narrows the SPECIES group to a list of species names (ALL, the default, writes every transported species). The species themselves — pollutants, MSX species, water age and temperature on the mesh — are the subject of Quality Chapter 10; the groundwater datasets and the aquifer block of the status report are described in Hydrology Chapter 9.
The solver is verified against the SWASHES compilation of analytic shallow-water solutions (Delestre et al., 2013), with the reference formulas implemented independently of the engine. The figures below are the relative \(L^1\) depth error and the mass-balance error over the run; steady cases are graded on the time mean of the final half of the simulation.
| Case | SWASHES § | rel. \(L^1\) depth error | mass error | graded against |
|---|---|---|---|---|
| Lake at rest, immersed bump | 3.1.1 | \(6.3\times10^{-11}\) | \(-7\times10^{-14}\) % | analytic |
| Lake at rest, emerged bump | 3.1.2 | \(1.4\times10^{-16}\) | 0 | analytic |
| Subcritical flow over a bump | 3.1.3 | 0.67 % | \(-6\times10^{-13}\) % | analytic |
| MacDonald 1000 m, subcritical | 3.2.1 | 2.0 % | \(2\times10^{-11}\) % | analytic |
| Transcritical, no shock | 3.1.4 | 27 % | \(-7\times10^{-13}\) % | baseline |
| Transcritical with shock | 3.1.5 | 12 % | \(-8\times10^{-13}\) % | baseline |
| Stoker wet-bed dam break | — | 6.2 % | \(-1\times10^{-12}\) % | baseline |
| Ritter dry-bed dam break | — | 10 % | \(-3\times10^{-12}\) % | baseline |
| Thacker planar, 1D | — | 78 % | \(7\times10^{-13}\) % | baseline |
| Thacker radial, 2D | — | 29 % | \(-7\times10^{-14}\) % | baseline |
| Thacker planar, 2D | — | 43 % | \(7\times10^{-14}\) % | baseline |
| MacDonald 1000 m, supercritical | 3.2.1 | 19 % | \(-2\times10^{-13}\) % | expected failure |
Four observations follow.
Well-balancedness is exact. Both lake-at-rest cases sit at rounding, including the emerged bump, which is a wetting and drying problem. The C-property of §9.5.3 holds exactly.
Mass conservation is exact. Every case closes to \(10^{-11}\) % or better, which is round-off for the volumes involved. This is the structural guarantee of (9-13) and (9-14), and it holds through wetting, drying, positivity rescaling and boundary clamping.
Accuracy is good where the approximation holds and degrades where it does not. The subcritical bump and the subcritical MacDonald channel meet analytic tolerances. The transcritical, dam-break and oscillating cases are graded against recorded baselines instead, because the missing convective term is a physics limit rather than a discretization error that a finer mesh would remove.
The failure mode is specific and worth recognizing. On the subcritical bump, the computed free surface is dead flat over the crest, where the analytic solution dips. That is not a defect: a flat \(\eta\) is the exact frictionless steady state of (9-2), because the dip is \(\Delta(v^{2}/2g)\) and there is no \(q^{2}/h\) term to produce it. The residual 0.7 % error is that dip and nothing else. The supercritical MacDonald channel is recorded as an expected failure for the same reason, one order more severely: with no convective inertia the fully supercritical profile never steadies at all, and develops a roll-wave-like unsteadiness. If your problem looks like that case, this is not the solver for it.
The October 2026 review separated numerical corrections from arithmetic-preserving CPU optimization. The corrected implementation uses connected wet stencils, a conditioned shoreline fit, a positive face-depth polynomial and a consistent within-cell bed force. Resting lakes alone cannot expose all defects: a method may preserve a constant surface but still give the wrong mass flux or bed force for a sloping surface. The focused tests therefore also check linear wet-face mass flux, cell-to-face pressure integration, bounded thin-depth forces, phase/depth/velocity in planar Thacker oscillations, and distorted triangular, quadrilateral and mixed shorelines.
The final qualification rebuilt the retained CPU solver and passed 169 of 170 tests across 15 suites; one unavailable Kokkos OpenMP plugin test was skipped. Twelve fresh paired full-engine runs matched API states, water accounting, normalized reports and all HDF5 contents exactly for local-inertial rain/infiltration/drain/pipe workloads with zero, eight and 32 transported species. The preceding optimization qualification also has 70 exact analytical pairs and 12 exact coupled pairs, including FULL_SWE, FLAT/VFR storage and a wet Bellinge interval. Those results were carried forward after verifying source and evidence hashes, rather than counted as newly executed tests.
For accuracy per CPU cost, the reference is the phase-7 reconstruction at 128 divisions and the final implementation at 96 divisions, over three Thacker periods, one CPU thread and second-order reconstruction. The final mesh has 43.75% fewer cells. Both endpoint and time-mean depth errors improve in all four cases:
| Bowl / mesh | Reference endpoint depth L1 | Final endpoint depth L1 | Median paired CPU reduction |
|---|---|---|---|
| Radial / triangles | 3.88% | 2.97% | 61.1% |
| Radial / quadrilaterals | 7.29% | 4.90% | 60.1% |
| Planar / triangles | 7.84% | 6.30% | 59.9% |
| Planar / quadrilaterals | 10.70% | 9.51% | 60.9% |
The depth error is the area-weighted absolute depth error divided by the area-weighted exact depth. A separate history run computes the arithmetic mean of the sampled relative errors. All eight comparison runs keep nonnegative depth, with maximum absolute relative volume error \(1.34\times10^{-14}\). Three randomized timing pairs after one discarded warm-up pair per case provide the CPU ratios; history diagnostics are excluded from the timing executable. This is a depth-accuracy comparison, not equal velocity error, same-grid speedup or a general convergence order.
Figure 9-8 Final CPU cost and depth error relative to the earlier reconstruction, at different mesh resolutions. Points are measured pairs; bars are paired medians.
The wet Bellinge control uses 25,600 loaded cells, a cold-start ten-minute storm interval, VFR storage and second-order FULL_SWE. It advances 2,728 internal steps and 104,574,000 face evaluations, with 99.5% mean active cells. Its exact before/after agreement validates the CPU optimization, not catchment calibration or agreement with observed flooding. The separately timed early Bellinge interval is dry and must not be used as wet-storm performance evidence.
Long-period damping remains, especially on coarse meshes. Wet-front velocity safeguards and the reported depth improvements do not constitute an arbitrary-mesh stability proof. GPU execution was unavailable and is not validated by CPU agreement. The engine manual's CHAPTER 7 — Two-Dimensional CPU Implementation and Qualification describes the implementation, benchmark controls and reproduction prerequisites. Detailed source manifests, tests and measurements are in reviews/cpu_2d_2026_10_03/accuracy_phase9, cpu_phase10, cpu_phase11 and cpu_phase12 in the engine repository. Full-engine review binaries use a consistent frozen supporting-engine build and do not qualify unrelated concurrent 1D changes.
All keys live in [2D_OPTIONS]. This section is an index, not a reference: the theory each key selects or tunes is in the section named beside it, and the accepted values and defaults are documented once, in the engine manual (Section: [2D_OPTIONS]), so that a default changed there is not contradicted here. The 2D grammar as a whole is CHAPTER 2 - Input File Reference, with the boundary and mesh sections at Section: [2D_BOUNDARY_CONDITIONS], Section: [2D_QUADS] and Section: [2D_MESH_FILE].
| Key | Theory section |
|---|---|
MAX_TIMESTEP | §9.5.5; also caps the co-advance span of §9.7.3 |
DRY_DEPTH | §9.5.9, Table 9-1 |
COUPLING_CD | §9.7.1, §9.7.5 |
COUPLING_SYNC | §9.7.3 |
LIMITER_EPSILON | §9.9 |
FLUX_DH_EPS | §9.5.11, §9.6 |
RAINFALL_MODE | §9.8; Chapter 8: Distributed Surface Processes on the 2D Mesh |
REPORT_2D | §9.9 |
CELL_CLOSURE | §9.4, §9.4.2 |
FACE_RECONSTRUCTION | §9.5.2 |
VFR_MIN_WET_FRAC | §9.4.1 |
OUTPUT_FILE | §9.9 |
INTEGRATOR | §9.5; EXPLICIT is the only value, the others are retired (below) |
MOMENTUM_EQUATION | §9.2, Table 9-2 |
RECONSTRUCTION_ORDER | §9.5.10, §9.5.6 |
FRONT_REBUILD | §9.5.7 |
THETA | §9.5.1 |
THETA_REFERENCE_STEP | §9.5.1; zero keeps the per-firing weight, positive seconds define physical-time retention |
CFL_NUMBER | §9.5.5; clamped to ½ under FULL_SWE, §9.5.10; the diffusion safety fraction of §9.5.11 |
H_MOVE | §9.5.7 |
LTS_TIERS | §9.5.6 |
FROUDE_MAX | §9.5.1 |
ADVECTION | deprecated spelling of the convective term — §9.2.3, Table 9-2, and the retired keys below |
COUPLING_AREA | §9.7.5 |
COUPLING_IN_FLOODING | §9.7.6 |
BACKEND | §9.11.1 |
OUTPUT_PRECISION | §9.9 |
OUTPUT_COMPRESSION | §9.9 |
REPORT_2D_VARIABLES | §9.9, Table 9-4 |
REPORT_2D_SPECIES | §9.9; Chapter 10: Surface Quality and Transport on the 2D Mesh |
REPORT_2D_STEP | §9.9 |
INFILTRATION | Chapter 8: Distributed Surface Processes on the 2D Mesh |
INFIL_STEP | Chapter 8: Distributed Surface Processes on the 2D Mesh |
INFIL_DEFAULT_METHOD | Chapter 8: Distributed Surface Processes on the 2D Mesh |
INFIL_DESTINATION | Chapter 8: Distributed Surface Processes on the 2D Mesh; the AQUIFER_2D destination is Chapter 9: Spatially Explicit Groundwater on the 2D Mesh |
EVAPORATION | §9.8 |
TRANSPORT_POLLUTANTS | Chapter 10: Surface Quality and Transport on the 2D Mesh |
TRANSPORT_MSX | Chapter 10: Surface Quality and Transport on the 2D Mesh |
TRANSPORT_AGE | Chapter 10: Surface Quality and Transport on the 2D Mesh |
TRANSPORT_TEMPERATURE | Chapter 10: Surface Quality and Transport on the 2D Mesh |
GROUNDWATER | Chapter 9: Spatially Explicit Groundwater on the 2D Mesh |
GW_ET | Chapter 9: Spatially Explicit Groundwater on the 2D Mesh |
DISPERSION | Chapter 10: Surface Quality and Transport on the 2D Mesh |
Table 9-5 The forty-two live [2D_OPTIONS] keys and the theory each selects
The computational backend is chosen by BACKEND. The built-in CPU marcher is OpenMP-threaded on the host and uses the project's THREADS setting; OMP, CUDA, HIP and SYCL name a Kokkos plugin and load it outright (a plugin that is not installed, or reports no usable device, falls back to CPU with a notice — never a hard failure). AUTO, the default, prefers a device plugin when one is installed and the mesh is above the device floor (OPENSWMM_2D_MIN_PARALLEL_CELLS_DEVICE, default 10 000 cells), then the OpenMP plugin above its own, much higher floor (OPENSWMM_2D_MIN_PARALLEL_CELLS, default 50 000), else CPU. The OPENSWMM_2D_BACKEND environment variable overrides the option when set (the same precedence OPENSWMM_FV_BACKEND has over FV_BACKEND). The OpenMP floor is deliberately high — on a 25 000-cell coupled model the OpenMP plugin measured an order of magnitude slower than the built-in marcher, because the plugin pays a launch cost on every substep while the built-in path is already threaded. Whether a discrete GPU pays off is likewise problem-dependent (§9.7.3: every routing step is a host↔device round trip, and the junction-exchange and boundary passes run serially), so a model that measures slower on AUTO should pin BACKEND CPU.
The plugins implement the all-triangle local-inertial scheme without transport. A model that mixes in quads (§9.3.1), selects FULL_SWE or DIFFUSIVE_WAVE (§9.2), or carries species, age or temperature rows is refused by the plugin loudly — under AUTO as well, since a model that silently lost its plugin would look like a performance regression — and served by the CPU marcher. There is no silent physics substitution.
Keys retired with the earlier implicit integrator are accepted with a warning and ignored on file load, so legacy models still open; on the programmatic option-set path they are errors. Each row names where the behaviour it once configured is now described.
| Key | Theory section |
|---|---|
MIN_TIMESTEP | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9; the marcher's step is §9.5.5 |
REL_TOLERANCE | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
ABS_TOLERANCE | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
MAX_CVODE_STEPS | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
MAX_KRYLOV_DIM | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
LINEAR_SOLVER | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
PRECONDITIONER | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
JACOBIAN | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
ATOL_AREA_REF | retired with the CVODE/ARKODE integrators — Chapter 10: Planned and Retired Formulations §10.9 |
COUPLING_INTERVAL | superseded by the per-routing-step co-advance and COUPLING_SYNC — §9.7.3 |
COUPLING_WINDOW | superseded by the per-routing-step co-advance and COUPLING_SYNC — §9.7.3 |
ACTIVE_SET | the flux-active set is always on — §9.5.7 |
ACTIVE_SET_HALO | the halo is the one-ring (five rings under FRONT_REBUILD) — §9.5.7 |
MOMENTUM | superseded by MOMENTUM_EQUATION — §9.2; the IMEX blend it selected is Chapter 10: Planned and Retired Formulations §10.9 |
INTEGRATOR other than EXPLICIT | the CVODE/ARKODE selections — Chapter 10: Planned and Retired Formulations §10.9 |
ADVECTION | deprecated 2026-09-06 in favour of MOMENTUM_EQUATION FULL_SWE — §9.2.2, Table 9-2; still honoured with a warning so an old deck keeps its physics |
Table 9-6 The sixteen retired [2D_OPTIONS] keys and values
MOMENTUM_EQUATION FULL_SWE (§9.2.2) restores the term, at the cost of a half-Courant step and CPU-only execution.DIFFUSIVE_WAVE the penalty is quadratic in the cell size (9-43).RECONSTRUCTION_ORDER 2 runs in global-step mode and falls back to first order when overland transport is active. Groundwater RK2 has a standalone CPU water-only profile; reactions/heat and borrowed owners remain guarded (§9.5.10). Long-period damping remains, and the shoreline fallback does not establish uniform second-order convergence on arbitrary meshes.COUPLING_SYNC > 0 the 1D heads a 2D advance sees are held for the batch. The error grows with the span.SWMMEngine::saveCheckpoint and SWMMEngine::loadCheckpoint, with the same attached 2D common-source stage kept alive for the session. Restore is same-platform and requires a freshly initialized matching model and mesh/contact identity. The checkpoint profile accepts supported in-memory climate sources; active runtime forcing overrides and file-backed climate sources are refused. The input-file SAVE HOTSTART / USE HOTSTART commands do not create or restore this joint checkpoint. For the qualified attached 2D session, the engine checkpoint preserves shared 1D node/link temperatures and bed state together with the mesh soil columns and pending thermal receipts. Mapped open conduits using MESH_CONTACTS can therefore resume their channel-to-ground hyporheic heat exchange with the same mesh/contact identity. Restore is same-platform into a freshly initialized matching engine and attached mesh session. Direct hydraulic-component checkpoint files do not include thermal state.Audusse, E., Bouchut, F., Bristeau, M.-O., Klein, R., and Perthame, B. (2004). "A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows." SIAM Journal on Scientific Computing, 25(6), 2050–2065.
Barth, T. J., and Jespersen, D. C. (1989). "The design and application of upwind schemes on unstructured meshes." AIAA Paper 89-0366, 27th Aerospace Sciences Meeting, Reno.
Bates, P. D., Horritt, M. S., and Fewtrell, T. J. (2010). "A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling." Journal of Hydrology, 387(1–2), 33–45.
Begnudelli, L., and Sanders, B. F. (2006). "Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying." Journal of Hydraulic Engineering, 132(4), 371–384.
Begnudelli, L., and Sanders, B. F. (2007). "Conservative wetting and drying methodology for quadrilateral grid finite-volume models." Journal of Hydraulic Engineering, 133(3), 312–322.
Belikov, V. V., Ivanov, V. D., Kontorovich, V. K., Korytnik, S. A., and Semenov, A. Y. (1997). "The non-Sibsonian interpolation: A new method of interpolation of the values of a function on an arbitrary set of points." Computational Mathematics and Mathematical Physics, 37(1), 9–15.
Bowyer, A. (1981). "Computing Dirichlet tessellations." The Computer Journal, 24(2), 162–166.
Bruwier, M., Archambeau, P., Erpicum, S., Pirotton, M., and Dewals, B. (2017). "Shallow-water models with anisotropic porosity and merging for flood modelling on Cartesian grids." Journal of Hydrology, 554, 693–709.
de Almeida, G. A. M., Bates, P., Freer, J. E., and Souvignet, M. (2012). "Improving the stability of a simple formulation of the shallow water equations for 2-D flood modeling." Water Resources Research, 48, W05528.
de Almeida, G. A. M., and Bates, P. (2013). "Applicability of the local inertial approximation of the shallow water equations to flood modeling." Water Resources Research, 49(8), 4833–4844.
Delestre, O., Lucas, C., Ksinant, P.-A., Darboux, F., Laguerre, C., Vo, T.-N.-T., James, F., and Cordier, S. (2013). "SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies." International Journal for Numerical Methods in Fluids, 72(3), 269–300.
Gottlieb, S., Shu, C.-W., and Tadmor, E. (2001). "Strong stability-preserving high-order time discretization methods." SIAM Review, 43(1), 89–112.
Hunter, N. M., Horritt, M. S., Bates, P. D., Wilson, M. D., and Werner, M. G. F. (2005). "An adaptive time step solution for raster-based storage cell modelling of floodplain inundation." Advances in Water Resources, 28(9), 975–991.
Jawahar, P., and Kamath, H. (2000). "A high-resolution procedure for Euler and Navier–Stokes computations on unstructured grids." Journal of Computational Physics, 164(1), 165–203.
Kumar, M., Duffy, C. J., and Salvage, K. M. (2009). "A second-order accurate, finite volume-based, integrated hydrologic modeling (FIHM) framework for simulation of surface and subsurface flow." Vadose Zone Journal, 8(4), 873–890.
Liang, Q., and Marche, F. (2009). "Numerical resolution of well-balanced shallow water equations with complex source terms." Advances in Water Resources, 32(6), 873–884.
Perot, B. (2000). "Conservation properties of unstructured staggered mesh schemes." Journal of Computational Physics, 159(1), 58–89.
Sanders, B. F., Schubert, J. E., and Gallegos, H. A. (2008). "Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling." Journal of Hydrology, 362(1–2), 19–38.
Stelling, G. S., and Duinmeijer, S. P. A. (2003). "A staggered conservative scheme for every Froude number in rapidly varied shallow water flows." International Journal for Numerical Methods in Fluids, 43(12), 1329–1354.
Thacker, W. C. (1981). "Some exact solutions to the nonlinear shallow-water wave equations." Journal of Fluid Mechanics, 107, 499–508.
Toro, E. F. (2001). Shock-Capturing Methods for Free-Surface Shallow Flows. Wiley, Chichester.
Watson, D. F. (1981). "Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes." The Computer Journal, 24(2), 167–172.
The CPU full shallow-water solver accepts the following [2D_OPTIONS] keys. Defaults preserve the existing solver path.
| Key | Values and units | Default |
|---|---|---|
CORIOLIS | OFF, LATITUDE, CONSTANT_F | OFF |
CORIOLIS_LATITUDE | Finite degrees north, −90 to 90; required in latitude mode | Unset |
CORIOLIS_F | Finite signed s⁻¹; required in direct-parameter mode | Unset |
TURBULENCE_MODEL | NONE, SMAGORINSKY, MIXING_LENGTH | NONE |
SMAGORINSKY_CS | Finite nonnegative dimensionless value | 0.10 |
EDDY_VISCOSITY_BACKGROUND | Finite nonnegative m²/s; inactive with NONE | 0 |
MIXING_LENGTH_COEFFICIENT | Finite nonnegative dimensionless value; active only with MIXING_LENGTH | 0.5 |
UNSET clears either optional Coriolis input through the extended option API. Inactive values can remain authored; the active mode selects the parameter. The input writer omits default process switches and preserves authored parameters. Generic extended-option access also makes the controls available to programmatic clients without changing C ABI option structures.
Qualification is limited to CPU FULL_SWE, LTS_TIERS 1, first-order or second-order reconstruction. Unsupported momentum modes and coupled surface/source profiles produce a diagnostic. Inertial and diffusive-wave closures and accelerator kernels do not acquire these processes implicitly. A selected accelerator follows the established full-SWE CPU fallback policy.
In a right-handed horizontal frame, Coriolis adds (f hv, -f hu) to the two momentum equations, with f = 2 Ω sin(latitude) and Ω = 7.2921158553 × 10⁻⁵ s⁻¹. Exact local momentum rotations are symmetrically composed around the hydrodynamic step. Isolated rotation preserves kinetic energy up to floating-point roundoff. This source integration is not a claim of exact geostrophic well-balancing; rotating pressure equilibria must be checked for convergence.
For Smagorinsky, Δ = sqrt(A), νt = ν0 + (Cs Δ)² |S|, and
|S| = sqrt(2 ux² + 2 vy² + (uy + vx)²).
Here ν0 is EDDY_VISCOSITY_BACKGROUND, a constant horizontal eddy viscosity. With Cs = 0 it gives constant-viscosity mixing; under Froude similarity a physically equivalent value scales as length^1.5, so one value does not transfer between model scales.
TURBULENCE_MODEL MIXING_LENGTH uses the same stress operator with a shear-based mixing length: νt = ν0 + l² |ω| with l = min(c h, κ d), where h is the cell depth, c is MIXING_LENGTH_COEFFICIENT, κ = 0.41, d is the distance from the cell centroid to the nearest mesh boundary edge and |ω| = |vx − uy| uses the same least-squares gradients as |S|.
c = 0.5.|ω| instead of |S| gives no eddy viscosity in irrotational strain, such as converging flow into a contraction, and keeps it in shear layers and curved jumps. A strain-based depth length put up to several m²/s of viscosity into a deep reservoir and its outlet during UK Test 6A development.κ d limits the length near walls, as in common depth-averaged mixing-length models. Here d is measured to every exterior edge, including open boundaries; there the cap only acts within about c h/κ (≈ 1.2 depths) of the edge.νt depends on depth, distance and the flow, not on cell size, so refining a mesh converges toward a fixed eddy viscosity. The Smagorinsky value falls with cell area.The model has no transported state and no bed-generated term. It is intended for shallow flows, where horizontal eddies are limited by depth.
The added term is the conservative depth-integrated stress divergence div(2 h νt S). Unrestricted least-squares velocity gradients are separate from shock-limited MUSCL reconstruction. Shared-face gradients include an over-relaxed normal correction reproducing the centroid velocity difference on nonorthogonal cells. Face coefficients use interpolated cell viscosity, connected wet depth and conveyance. Each face transfers equal/opposite momentum; dry sills, holes and zero-conveyance edges transmit no turbulent momentum. Rank-deficient corner gradients borrow a connected complete stencil where available. Otherwise their cell gradient remains zero. Walls and open boundaries have zero turbulent traction; no wall-friction or wall-damping model is included.
The mixing operator is integrated with SSP-RK2 at fixed water storage inside symmetric source splitting. Its timestep uses absolute row sums of the reconstructed face operator, including cross stresses. Nonlinear stage coefficients are reevaluated; a tighter bound or increasing kinetic energy retries the unpublished momentum substep. Viscosity is never clipped to meet a desired timestep. A mesh/state for which no dissipative admissible substep can be found is refused. Both stage checks and energy checks are necessary: conservation alone does not prove a skew-mesh stress discretization is accurate or dissipative.
The optional global mode includes all cells in the schedule and evaluates exterior Riemann states at the same beginning-of-stage time as interior states. This prevents boundary splitting from generating a nonuniform flow in the uniform inertial-oscillation test. The disabled path retains its prior numerical behavior.
REPORT_2D_VARIABLES TURBULENCE adds Mesh2_face_eddy_viscosity (m²/s) and Mesh2_face_strain_rate (s⁻¹) when mixing is enabled. They are computed from the accepted reported state. Root attributes record the effective f, turbulence model, the Smagorinsky coefficient and filter or the mixing-length coefficient, and the background viscosity (eddy_viscosity_background_m2ps). Optional-process reported velocity comes from accepted cell momentum; edge fluxes continue to represent hydraulic transfers and boundary window means. No plugin snapshot or C ABI structure is enlarged by this output extension.
A depth-averaged closure does not resolve vertical turbulence, nonhydrostatic pressure or three-dimensional obstacle wakes. Bed-generated mixing, dynamic Smagorinsky, scalar turbulent diffusivity and wall damping are separate future closures. Suggested coefficient sensitivity is 0.05, 0.10 and 0.20; mesh and timestep convergence remain necessary. Adding mixing does not automatically resolve the UK Test 6 gauge discrepancy.
WALL remains impermeable and free-slip. Select a solid exterior edge with WALL_DRAG Cw in [2D_BOUNDARY_CONDITIONS] to add tangential traction:
The example coefficient is a diagnostic value, not a material recommendation. The convention is tau_t/rho = -Cw |u_t| u_t (no factor of one half; Darcy friction factor f corresponds to Cw=f/8). Cw must be finite and nonnegative; zero is free-slip. It is independent of bed Manning roughness. The finite-volume force uses edge length times wetted wall height; VFR integrates the edge endpoint elevations against the current water surface. Only wall-parallel momentum is changed. Exact scalar drag decay, composed symmetrically over wall faces, removes kinetic energy without reversing the wall-parallel component in a single wall update. Normal impermeability and resolved pressure/form drag are retained. Do not also assign a duplicate empirical contraction/form-loss term for the same resolved geometry.
Initial support: CPU MOMENTUM_EQUATION FULL_SWE, LTS_TIERS 1, orders 1/2, uncoupled surface hydraulics. Smagorinsky may be on or off. Its interior stress has zero exterior traction; WALL_DRAG supplies the selected wall traction once. Unsupported closures, coupling and wall checkpoint profiles are refused. Wall coefficients can be authored through swmm_2d_get_edge_wall_drag and swmm_2d_set_edge_wall_drag (OPENED/BUILDING only), and Python Surface2D.get_edge_wall_drag / set_edge_wall_drag. INP, mesh sidecars and GeoPackage preserve the prescription. Runtime wall edits are unsupported.
The run log reports selected face count and cumulative wall energy removal divided by density, in m5/s2. A constant prescribed Cw is a partial-slip drag closure; it is not a smooth/rough log-law wall function or a resolved no-slip boundary layer. Mesh, timestep and coefficient sensitivity remain necessary.
For authoring APIs, first call swmm_2d_prepare_for_edit (Python Surface2D.prepare_for_edit()) after opening the model. Positive wall drag and nonempty groups are checked against authored edge connectivity, before runtime neighbours have been built. Triangle INP/sidecar/GeoPackage wall round trips are qualified. For quad wall edge 3, use INP: the existing GeoPackage boundary table restricts local edge indices to 0–2.
For subsecond reports, use an explicit clock such as REPORT_STEP 0:00:00.1. Fractional clock seconds are preserved for REPORT_STEP, WET_STEP and DRY_STEP, including save/reopen. Bare numeric values for these options retain legacy decimal-hour semantics; for example, REPORT_STEP 0.1 means six minutes. ROUTING_STEP retains its numeric-second convention. The Test 6 sensitivity runner verifies the actual native HDF5 report grid before starting the long cases.
1.13.2