20b. 2D methods: equations, numerics and validation

This chapter is for the engineer who has to defend a 2D result. It gives the equations the solver integrates, the exact update formulas it applies, how the time step, wet/dry cells, boundaries, rain and infiltration are handled, the 1D-2D interface formulas with their coefficients and sources, the two file formats the 2D engine owns, the validation cases with the numbers the test suite asserts, and the limitations. The user-facing side — how to set up a grid, draw interfaces and read the results in the app — is in chapter 20.

The solver is written from the published papers cited at the end; no code was taken from any other flood model.

1. Governing equations

Depth-averaged flow over a fixed bed obeys the shallow-water equations. With h the depth, z the bed elevation, η = z + h the water surface, q = (qx, qy) the unit discharge (m²/s or ft²/s) and u = q / h the velocity:

∂h/∂t + ∂qx/∂x + ∂qy/∂y = R − I                                (continuity)
∂q/∂t + ∂(q u)/∂x + ∂(q v)/∂y + g h ∇η + g h S_f = 0            (momentum)

R is the rain rate on the cell, I the infiltration rate, g gravity and S_f the friction slope from Manning's equation,

S_f = n² |q| q / (k² h^{10/3})

with n Manning's roughness and k the unit constant: k = 1 in metric units (g = 9.81 m/s²) and k = 1.486 in US customary units (g = 32.174 ft/s²). The system is chosen by the model's FLOW_UNITS: CMS, LPS and MLD are metric, everything else is US. All lengths on the grid are the DEM's units, which must be the model's length units.

The solver offers two momentum treatments:

2. Discretisation

2.1 Grid

The DEM's cells are the control volumes: square, side Δx, row-major with the top row first (the raster convention). Depth h, bed z and roughness n live at cell centres. For the local-inertial scheme the unit discharges live on the faces, staggered: qx on the (ncols+1) × nrows vertical faces (positive to +x, east) and qy on the (nrows+1) × ncols horizontal faces (positive towards increasing row index, i.e. south; the results file converts to world +y). A cell whose DEM value is no-data is a wall: every face it touches carries zero flux and it is never wet.

2.2 Local-inertial flux update

For the face between cells L and R (west and east, or north and south), Bates et al. (2010) eq. 11–12 with de Almeida et al. (2012) eq. 12:

h_f  = max(η_L, η_R) − max(z_L, z_R)                    (effective face depth)
q_c  = θ q^n + ½ (1 − θ) (q^n_prev + q^n_next)          (q-centred average)
q^{n+1} = [ q_c − g h_f Δt (η_R − η_L) / Δx ]
          / [ 1 + g Δt n_f² |q^n| / (k² h_f^{7/3}) ]

where q^n_prev and q^n_next are the fluxes on the two neighbouring faces along the same axis, n_f = ½(n_L + n_R), and θ = 0.7 (θ = 1 recovers the original Bates scheme). The friction term is evaluated implicitly in q^{n+1}, which is what keeps thin sheets of water on rough ground stable. If h_f ≤ dry_depth the flux is zero.

The scheme is well balanced by construction: the driving term is the free-surface gradient, which vanishes on a lake at rest whatever the bed does (§6.1).

2.3 Continuity update and the outflow limiter

h^{n+1} = h^n + (Δt/Δx) (q_W − q_E + q_N − q_S) + R Δt − I Δt

with the four face fluxes taken as inflow-positive. Before the update, every cell's total outflow over the step is compared with what it holds; if Δt/Δx · Σ q_out > h, all of that cell's outgoing fluxes are scaled by h / (Δt/Δx · Σ q_out). A face is scaled by its upwind cell only, so the volume leaving one cell is exactly the volume entering its neighbour and the limiter conserves mass. The scaled fluxes are the ones carried into the next step's momentum equation. The limiter is a safety net for the rare cell that empties within a step; it is never active in the validation runs except at a dam-break tip.

2.4 Full-dynamic (HLL) update

The conserved variables (h, hu, hv) live at cell centres. The update comes in two orders ([RUN] ORDER, §5.1; Second order in 2D Setup): first order, described in this section, where each face sees the two cell averages and one forward-Euler step is taken; and second order (the default), which feeds the same flux and source terms with reconstructed face states and takes a two-stage step (§2.4a). For each face the two bed elevations are reconciled by hydrostatic reconstruction (Audusse et al. 2004):

z_f = max(z_L, z_R),  h_L* = max(0, η_L − z_f),  h_R* = max(0, η_R − z_f)

and the HLL flux F = (F_h, F_hu, F_hv) is formed from (h_L*, u_L, v_L) and (h_R*, u_R, v_R) with wave speeds S_L = min(u_L − c_L, u_R − c_R), S_R = max(u_L + c_L, u_R + c_R), c = √(g h*), and Toro's dry-bed estimates (S_L = u_R − 2c_R, S_R = u_R + c_R when the left is dry, mirrored when the right is dry):

F = F_L               if S_L ≥ 0
F = F_R               if S_R ≤ 0
F = [S_R F_L − S_L F_R + S_L S_R (U_R − U_L)] / (S_R − S_L)   otherwise

The transverse momentum is carried as F_hv = F_h · v_upwind. The cell update is

h^{n+1}   = h  − (Δt/Δx) Σ F_h (outward)
hu^{n+1}  = hu − (Δt/Δx) (F_hu,E − F_hu,W + F_hv,S − F_hv,N)
            + (Δt/Δx) ½ g [(h_E,L*)² − (h_W,R*)²]

(likewise for hv with the north/south faces), where h_E,L* is this cell's reconstructed depth on its east face and h_W,R* on its west face; the pressure difference of the reconstructed depths is what balances the bed slope exactly at rest. Manning friction is applied semi-implicitly after the update, hu ← hu / (1 + Δt g n² |u| / (k² h^{4/3})). Walls (closed edges and no-data neighbours) use the reflecting flux (0, ½ g h², 0). The same outflow limiter as §2.3 is applied to the mass flux and, with the same factor, to the momentum fluxes.

2.4a Second-order reconstruction and time stepping

A first-order Godunov scheme puts the cell average on both sides of every face, which adds a numerical diffusion proportional to the cell size: a front or a hydraulic jump spreads over several cells, and a smooth wave loses amplitude as it travels. The second-order option removes the leading error term in both space and time, with the established ingredients for well-balanced, positivity-preserving shallow-water schemes:

Reconstruction (MUSCL, van Leer 1979). Along each axis every cell carries a limited linear profile of four quantities: the free surface η = z + h, the depth h, and the velocities u and v. The face states are the profile's values at the two faces. Reconstructing η rather than h alone is the surface-gradient method of Zhou et al. (2001): a flat surface has zero slope over any bed, so still water stays still. The face bed is then z_face = η_face − h_face (Audusse et al. 2004, §3), so the face depth is the reconstructed depth, which the limiter keeps non-negative. Velocities, not discharges, are reconstructed: a limited velocity stays between its neighbours' even where the depth is small, where a reconstructed discharge divided by a small face depth can give a spurious large speed.

Limiter. The monotonised-central limiter (MC; van Leer 1977), written for unequal cells. For a cell of extent d with centre gaps g⁻, g⁺ to its neighbours and differences Δ⁻ = w − w⁻, Δ⁺ = w⁺ − w, the change from the centre to the positive face is

δ = sign(Δ⁻) · min( |a Δ⁻ + b Δ⁺|, |Δ⁻|, |Δ⁺| )    if Δ⁻ Δ⁺ > 0, else 0
a = d g⁺ / (2 g⁻ (g⁻ + g⁺)),   b = d g⁻ / (2 g⁺ (g⁻ + g⁺))
w_face⁺ = w + δ,   w_face⁻ = w − δ

a Δ⁻ + b Δ⁺ is d/2 times the central slope (Δ⁻/g⁻ · g⁺ + Δ⁺/g⁺ · g⁻)/(g⁻ + g⁺), which is exact for a quadratic on any spacing; on a uniform grid a = b = ¼ and this is the textbook MC limiter. The bounds |Δ⁻| and |Δ⁺| are MC's "twice the one-sided slope", stated as the property they guarantee on any spacing: no face value passes the neighbour's average on its side. That is what keeps face depths positive and the scheme free of new extrema on a graded mesh, where the uniform-grid factor of 2 would let a face overshoot a neighbour that is closer than the cell is wide. The gaps and widths are the mesh's true ones (§3), computed once per column and row.

MC was chosen over minmod on the measured errors (§6.3, §6.3a–c): on the same cases minmod's errors were 1.6 times MC's on Ritter's dam break (L1), 1.5 times on the wet-bed bore, 3 times on the smooth wave at the finest grid (where minmod's clipping of smooth extrema also lowered the observed order to 1.73–1.91), and 1.1 times on Thacker's bowl. Both are positive and conservative by the same argument; minmod would only buy extra dissipation, which the wet/dry fallback below already supplies where it is needed.

Wet/dry fallback. Along an axis, a cell is reconstructed only when it and both of its neighbours along that axis are deeper than 10 × the dry depth (3 cm at the default 3 mm). Anywhere else — at a shoreline, in thin sheet flow, beside a wall, a blocked cell or the grid edge — its face states are its averages, i.e. first order there. The threshold is in swmm::twod::solver::SECOND_ORDER_DEPTH_FACTOR.

Fluxes and the bed slope. The face states go through the hydrostatic reconstruction and the HLL flux of §2.4 unchanged. With face states that differ from the cell averages, the bed-slope source needs Audusse et al.'s (2004) second-order form: each face's flux is corrected on the cell's side by ½ g (h_face² − h*²), and the cell gains a centred term from its own two face states,

S_x = ½ g (h_W⁺ + h_E⁻) (z_W⁺ − z_E⁻)

where h_W⁺, z_W⁺ are the cell's depth and bed at its west face and h_E⁻, z_E⁻ at its east face. For still water (h + z equal on both faces) the corrected fluxes and S_x cancel exactly; at first order S_x is zero and the correction reduces to §2.4's pressure difference.

Time stepping. Heun's method in its strong-stability-preserving form (SSP-RK2; Shu & Osher 1988, Gottlieb, Shu & Tadmor 2001):

U¹    = U  + Δt L(U)
U²    = U¹ + Δt L(U¹)
Uⁿ⁺¹  = ½ (U + U²)

Each stage is a complete first-order-style update — the face fluxes from that stage's state, the outflow limiter (§2.3) for that stage, the conserved update and the semi-implicit friction. Because each stage is a forward-Euler step under the same CFL condition and the limiter keeps it non-negative, the average is non-negative too, and since both stages are in flux form the step conserves mass to round-off. Rain, infiltration, evaporation, sources and sinks are applied once, after the average, with the step's full Δt; the boundary volumes are the mean of the two stages'.

Cost. Two stages and the reconstruction make a step about 2–2.5 times the first-order cost at the same Δt (§6.9). The time step itself is the same (§2.5).

First order stays available, bit for bit as it was (the same flux, source and friction arithmetic in the same order), for comparison and for runs that must reproduce earlier results. A sidecar written before the option existed reads as first order (§5.1).

2.5 Time step

Adaptive, from the CFL condition on the fastest wave:

Δt = α Δx / √(2 g h_max)                     (local inertial)
Δt = ½ α Δx / max(|u| + √(g h))              (HLL)

with α the Courant number from the sidecar (COURANT, default 0.7). The HLL step carries the ½ because the unsplit two-dimensional update is stable for ν_x + ν_y ≤ 1; without it round-off on a lake at rest grows exponentially (that is how the factor was found). The local-inertial step carries a √2 for the same reason: its staggered forward-backward update of gravity waves is stable for c Δt √(1/Δx² + 1/Δy²) ≤ 1, which on square cells is c Δt / Δx ≤ 1/√2 ≈ 0.707. Bates et al. (2010) state the step as α Δx / √(g h) with α up to 0.7, which in two dimensions sits on that limit. The EA benchmark Test 2 (§6.11) showed what that costs: at the old step a still pond 0.3 m deep on 20 m cells locked into a checkerboard of ±0.3 m that was still there after 48 h, and the run settled correctly with an effective α of 0.65 or less. The √2 makes local-inertial runs take about 1.4 times as many steps as before. Both HLL orders use the same step: each SSP-RK2 stage is a forward-Euler step and is stable and positive under the same condition as the first-order scheme (Gottlieb, Shu & Tadmor 2001), and the outflow limiter of §2.3 guards positivity in every stage regardless. h_max is the largest depth on the grid after the previous step. MAX_DT caps the step; while the whole grid is dry the step is 1 s (capped the same way). Steps are shortened to land exactly on every output frame and every 1D synchronisation time, so the frames are at the times they claim.

2.6 Wet and dry cells

A cell is wet when h > dry_depth (DRY_DEPTH, default 0.003 in the DEM's units). Faces whose effective depth h_f is at or below it carry no flux; velocities are reported as zero on dry cells; the arrival-time grid records the first step at which a cell exceeds the threshold. Water below the threshold is still stored and still counted in the balance; it just does not move. The second-order HLL scheme falls back to first order wherever a cell or one of its neighbours along an axis is within ten dry depths of dry (§2.4a), so its reconstruction never reaches into a shoreline.

2.7 Boundaries

Every face on the grid's outer edge has a condition. [BOUNDARY] sets one for all of them; boundary segments (SEGMENT rows, §5.1) give chosen parts of the edge their own. A face no segment covers keeps the [BOUNDARY] setting, and where segments overlap the face goes to the one listed first (the setup warns, and warns of a segment that covers no face). A file without segments is the all-edges case exactly as it was: its runs are bit for bit those of the solver before segments existed, under both schemes and both HLL orders (checked on a bumpy basin with a source under every all-edges setting).

The all-edges settings.

Where a segment lies. EDGE side [from to] takes the faces of one side whose midpoints lie between two coordinates measured along it (y on the west and east sides, x on the north and south), or the whole side. LINE x1 y1 … takes each face whose midpoint a polyline segment spans, measured along the side, where that polyline segment runs along the side (no more than 45° off it) and passes within one cell of the edge there — the extent of the cell behind the face, normal to the edge. A line drawn down the edge, or through the edge cells' centres, so takes the faces beside it and not the first face of the adjoining side at a corner. Faces of no-data cells are left out.

Values. A constant, a [TIMESERIES] (dated points placed on the model clock, undated ones hours from the start, as for [SOURCES]), or a file of time value rows. A stage is an elevation in the DEM's units, linear between points and held at its first and last values outside them — a level is a state, and dropping it to zero past the last point would drain the grid. An inflow is in the model's flow units, converted as for sources, linear between points and zero outside them; negative values count as zero.

Stage. Local inertial: the FIXED_HEAD ghost — the edge cell's bed, depth max(0, H − z) — with the level H(t) sampled at the middle of the step, and the face discharge from the same momentum equation as any interior face (Bates et al. 2010 eq. 11 with de Almeida et al. 2012's θ-weighting). This is how LISFLOOD-FP's fixed- and varying-head boundaries work. HLL: the ghost has depth h_b = max(0, H − z) on the edge cell's bed, and its normal velocity comes from the Riemann invariant carried out of the grid by the characteristic that reaches the edge from inside, w + 2√(g h) with w the cell's outward velocity:

w_b = w + 2 (√(g h) − √(g h_b)),   held to |w_b| ≤ √(g h_b)

A subcritical boundary may impose one quantity and must take the other from the outgoing invariant; this is the standard treatment of a level boundary for a Godunov scheme (Toro 2001). The bound is there because a level cannot control supercritical flow: water leaving faster than critical passes the edge at critical depth (the free-overfall control), and water enters at most at critical. Where the edge cell already leaves supercritically, nothing comes back from outside and the ghost copies the cell. The tangential velocity is the cell's (zero gradient). The face flux is then the HLL flux between the ghost and the cell with the hydrostatic reconstruction, so still water at the boundary level is exactly at rest. The all-edges FIXED_HEAD keeps its ghost at rest, which treats the outside as a still reservoir: inflow then drops the level at the face below H, and a basin fills more slowly than its boundary level says. At second order each of Heun's stages sees the boundary at its own time — the predictor the level at the step's start, the corrector at its end — with the ghost made from the state that stage starts from.

Inflow. Over each step the hydrograph's volume V = ∫ Q dt is integrated exactly (trapezoids on the piecewise-linear series, as for a source) and shared among the segment's faces by conveyance,

w_i = K_i / Σ K_j,   K_i = L_i h_i^{5/3} / n_i,   q_i = V w_i / (Δt L_i)

with L_i the face length and h_i, n_i the depth and roughness of the cell behind it. Conveyance puts the water where a channel crossing the edge would carry it — the deep, smooth part — which is also how HEC-RAS shares a 2D flow-hydrograph boundary along its line (USACE HEC 2016). An equal share per length (the report's "vertical flow" spread, EA Test 3) puts as much on dry banks as in the channel and starts sheet flow there that the river would not; a share by depth concentrates less than the flow does. When no cell behind the segment holds water, the volume goes to its lowest cells (bed within the dry depth of the lowest), by face length, so the first water enters where a channel would. Local inertial: the face carries q_i into the grid, and its momentum reaches the next face through the q-centred θ-average. HLL: the mass flux is exactly q_i, and the normal momentum flux q_i²/h_b + g h_b²/2 (no tangential momentum: the flow enters normal to the edge) is that of the ghost depth satisfying the same outgoing invariant,

−q/h_b + 2 √(g h_b) = w + 2 √(g h)

or critical depth (q²/g)^{1/3} where that root would be supercritical (a dry edge cell, or one draining hard towards the edge). The left side is increasing and concave in h_b, so Newton's method from critical depth (left of the root) climbs onto it without overshooting. The segment books exactly V; at second order both of Heun's stages carry the same q_i, each over half the step.

Normal depth S. Local inertial: a ghost one gap beyond the edge whose bed and water surface are both lower by S Δ, so the face sees the edge cell's depth and a surface slope S; at steady state the momentum equation then gives Manning's q = (k/n) h^{5/3} S^{1/2}, with the scheme's inertia and implicit friction on the way there. HLL: a ghost on a bed lower by S Δ, with the cell's depth, moving outwards at the Manning velocity (k/n) h^{2/3} S^{1/2} — in uniform flow exactly the neighbour an interior cell has. Neither ever lets water in. Unlike OPEN, the slope is the one given (the channel's beyond the edge), not the local one, and under HLL it replaces a zero-gradient outflow that backs water up at a sloping edge (§6.7c).

Closed and open segments are CLOSED and OPEN above, on just their faces.

Time step. The water a stage segment holds over its edge cells counts in the CFL depth even while the grid is still dry, and under HLL the previous step's fastest ghost wave (at most 2√(g h_b)) counts in the wave speed. Neither applies without segments.

Books. Boundary fluxes are booked in the mass balance as inflow or outflow. Each edge face's volume — q Δt L from the limited fluxes the update applied, per Heun stage over half the step at second order — is also booked to its segment, or to (other edges). These are what the run summary lists per segment; they add up to the balance's boundary share to round-off.

2.8 Rain on the grid, infiltration, sources and sinks

2.9 Mass balance

The solver keeps five volumes in the DEM's cubic units: inflow (rain, sources, boundary inflow, exchange onto the surface), outflow (boundary outflow, sinks, exchange off the surface), infiltrated, evaporated, and the storage at the start. The reported error is

ε = (inflow − outflow − infiltrated − evaporated − (stored − stored₀)) / max(inflow, stored₀) × 100 %

Because every transfer is an addition to one cell and a subtraction from another (or from a booked total), the error is round-off: below 1e-9 % in every validation run. Boundary inflow and outflow are also kept per boundary segment (§2.7); their sums agree with the balance's boundary share to round-off.

With groundwater on (§2.11) the surface gains a sixth volume, the water exfiltrated back onto it (seepage), counted with the inflow:

ε_surface = (inflow + exfiltrated − outflow − infiltrated − evaporated − Δstored) / max(inflow + exfiltrated, stored₀)

and infiltration under a store is no longer a loss but a transfer, which the groundwater balance (§2.11.5) takes up.

2.10 Threads

The four passes of a step (fluxes, limiter, update, limiter applied) each run over bands of rows with std::thread::scope, every band writing a disjoint slice and reading the previous pass's arrays. Results do not depend on the thread count (THREADS, 0 = all cores). Every band costs a thread spawn per pass, a few hundred microseconds on Windows, so the grid is cut into at most one band per 10 000 cells, and grids under 40 000 cells run inline on the calling thread with no spawn at all. Before this was so, a 710-cell grid (EA Test 1) spent nearly all its time creating threads (40 s for a run that now takes 2 s), and an 80 000-cell grid stepped twice as fast on eight bands as on one band per core (24).

2.11 Groundwater

What it is. SWMM's own lumped, two-zone groundwater model (Rossman & Huber 2016, ch. 5) put under the 2D surface. Each store is one SWMM aquifer: an unsaturated upper zone of uniform moisture content θ above a saturated lower zone of depth d_L, between the aquifer bottom E_b and the ground surface E_s (total depth D = E_s − E_b). The 2D cells above a store recharge it with their infiltration, a full store refuses infiltration so the rain runs off, a high water table seeps back onto the low cells, and the store drains to a network node (SWMM's GWF equation) and into sewer pipes below the water table (GWI).

What it is not. It is not a 2D or 3D groundwater flow model. A store has one water table; there is no lateral groundwater flow between stores or between cells, no Darcy flow field, no anisotropy, no wells. Water moves sideways only as the GWF outflow to a node and as GWI into pipes. That is exactly SWMM's model, chosen because it is what a SWMM model's [AQUIFERS]/[GROUNDWATER] already parameterise; coupling to a groundwater flow code (MODFLOW-type) is a possible later option.

2.11.1 Stores and their parameters

2.11.2 Equations (SWMM's, unchanged)

Rates are per unit area of the store (length/s):

dθ/dt   = (f_I − e_U − f_U) / (D − d_L)
dd_L/dt = (f_U − e_L − f_L − f_G − f_GWI − f_S) / (φ − θ)

The pair is integrated over each groundwater step with SWMM's own solver: the fifth-order Cash–Karp Runge–Kutta scheme with Numerical Recipes' adaptive step control, relative tolerance 1e-4, first step = Δt — the same code path as EPA's odesolve.c — followed by gwater.c's bounds (θ ≥ WP; θ < φ and d_L < D by 0.001).

2.11.3 Coupling to the surface

The groundwater step (sidecar STEP, 60 s by default) is the synchronisation interval: a run's sync steps are cut to it, so in coupled runs the network sees the groundwater at the exchange times.

2.11.4 Coupling to the network

2.11.5 Conservation

SWMM's groundwater continuity is approximate (its report carries a groundwater continuity error): the states come from the ODE and the bounds, while the outflow is booked by the trapezoid of the end rates. Here every flux is booked once, as the volume actually moved — infiltration as the surface tallied it, seepage as the cells received it, lateral and GWI as the network is given them, ET and deep percolation as the Runge–Kutta integrals of their rates (fifth-order weights over the accepted steps) — and the state is then projected onto the stored depth the books imply. Stored depth is V = (D − d_L) θ + d_L φ, linear in d_L at fixed θ, so d_L = (V − Dθ)/(φ − θ) exactly; if that leaves [0, D) the upper zone takes the rest; water above the store's limits waits to seep back (above), and an overdraft (impossible with the limits above, but recorded as deficit if it ever happened) is never hidden. The correction is the size of the ODE tolerance — the water table stays within 5e-5 ft of SWMM's own (§6.12).

The groundwater balance is

ε_gw = initial + infiltrated + deficit − exfiltrated − e_U − e_L − deep − lateral − GWI − stored

reported as a percent of the larger of the infiltration and the initial storage; the run summary carries each term. The surface and groundwater balances together form the site's: rain, sources and boundary inflow in; boundary outflow, sinks, surface and subsurface ET, deep percolation, lateral flow and GWI out; the infiltration and exfiltration terms cancel. The boundary terms are the sums of the per-segment books (§2.7), so each segment's inflow and outflow is part of the combined balance as it is of the surface's.

2.11.6 Who owns a subcatchment's groundwater

A subcatchment that is a SWMM runoff subcatchment with a [GROUNDWATER] row and is covered by 2D cells could have its groundwater computed twice. The run refuses to start until the sidecar says which (OWNER; 2D Setup Owner):

Without such a subcatchment the owner need not be chosen.

2.11.7 Continuous runs: the state file

Every run with groundwater writes <model>.2d.gwstate (§5.3) at its end; pointing the next run's STATE (2D Setup Start state) at a copy of it starts each store where the last run left it, so a forecast chain carries the water table forward. Stores are matched by name (the subcatchment, or *DEFAULT); a store not in the file starts from its aquifer's own initial state. On the command line aquaswmm-swmm twod run ... --gw-state-in FILE --gw-state-out FILE does the same, like SWMM's USE HOTSTART and SAVE HOTSTART, and the Python package's aquaswmm.forecast.gwstate reads and writes the file and moves moisture and water table to and from a SWMM hotstart's groundwater records (chapter 25 §25.2), so the 2D state rides beside a hotstart chain.

3. Grid setup

Setup::build resolves the sidecar against the model:

  1. The DEM is read (.asc; .tif through the GIS chapter's GeoTIFF reader), cropped to WINDOW (cells whose centre lies inside) and resampled to CELL by bilinear interpolation between cell centres (nearest for the roughness-class raster, which must not be interpolated).
  2. The roughness grid is UNIFORM n, a GRID sampled onto the DEM cells (missing or non-positive values become 0.05 with a warning), or a CLASSES raster looked up in the CLASS id n table (unknown classes become 0.05 with a warning).
  3. Every node in [JUNCTIONS], [STORAGE], [OUTFALLS] and [DIVIDERS] with [COORDINATES] is placed in its cell. Its interface is the [NODES] row if there is one, SEALED if it is listed under [SEALED], else a default MANHOLE. The rim is the node's invert plus its full depth exactly as the engine computes it (EPA SWMM 5.2 link.c, link_validate): MaxDepth, raised to the crown of every link that meets the node — upstream ends of every link but pumps and bottom orifices, downstream ends of conduits only — except at storage units, which keep their MaxDepth. It is always a maximum, so a MaxDepth of 0 gives the highest crown, and a written depth below a crown is raised too. Nodes outside the DEM or on a no-data cell are warnings, not errors. Before the nodes are placed, each [BUILDINGS] BLOCK layer (.shp or .geojson polygons) removes from the flow field every cell whose centre lies inside a footprint, holes excepted (chapter 20, section 20.9b). A blocked cell is a no-data cell to the solver, so the wall rule of §2.7 applies at a building's edge and the mass balance needs no new term; a node under a footprint gets no interface and is named in the warnings.
  4. Each [BANKS] polyline (or the conduit's own centreline when no points are given) is rasterised with Bresenham's line algorithm between the cells of consecutive vertices, segments clipped to the grid; the crest at each cell is the given crest or the DEM.

4. 1D-2D interfaces

4.1 Where the formulas come from

The exchange formulas are those of Chen, Djordjević, Leandro & Savić (2007) and Leandro, Chen, Djordjević & Savić (2009): a free weir while the receiving side is below the crest, an orifice once the opening is submerged, and Villemonte's (1947) reduction for a submerged weir. Written with the constants explicit so one dimensionless coefficient c serves both unit systems:

free weir      Q = c · (2/3) · √(2g) · P · Δh^{3/2}
orifice        Q = c · A · √(2g Δh)
submerged weir Q = Q_free · [1 − (h₂/h₁)^{3/2}]^{0.385}     (Villemonte)

The weir forms take the weir coefficient c_w, the orifice form the orifice coefficient c_o. Each is set per interface (weir_coeff, orifice_coeff; a bank has only a weir) or, where an interface leaves it unset, comes from the model-wide [EXCHANGE] defaults, which are 0.6 for both. The setup resolves the defaults into every located interface, so the solver sees one number per form per interface. In metric units c_w (2/3) √(2g) = 1.77 at c_w = 0.6, the familiar broad-crested constant; in US units it is 3.21. The defaults are the literature values of Chen et al. (2007) and Leandro et al. (2009), not calibrated ones; a grate or a lid with a known rating is the place to change them. P is the opening's perimeter and A its clear area: the inlet's own for an INLET; for a MANHOLE, a round lid of 0.6 m (metric) or 2 ft (US), P = π D, A = π D² / 4.

4.2 Node interfaces

At each synchronisation time both models are at the same t. With H the node's head, z the DEM ground at its cell, d the surface depth there and η = z + d:

A sealed node is one whose lid is bolted or that is not on the surface at all (an underground storage). Outfall interfaces are sinks (§2.8).

4.3 Bank interfaces

For every cell under a bank line, the channel water surface is the head at the conduit's from-node and to-node interpolated by the cell's fraction of the line, η₁ = H_from + f (H_to − H_from). With η₂ the cell's surface and z_c the crest, h₁ = max(η₁, η₂) − z_c and h₂ = min(η₁, η₂) − z_c:

Q = c (2/3) √(2g) Δx h₁^{3/2} · [1 − (h₂/h₁)^{3/2}]^{0.385}   if h₂ > 0, else no reduction

positive from channel to surface when η₁ > η₂. Each cell's volume is split between the two end nodes by (1 − f, f) and set as (negative) lateral inflow in tight mode. Flow back into the channel is bounded by what the cell holds.

Spill onto the surface is capped the way a node's formula surcharge is. In tight mode every node has one allowance per exchange interval — the one formula surcharge uses (§4.2) — shared by everything that withdraws from it: its own lid first, then the bank lines on its links in the order they are listed. A bank's spill is charged to its end nodes by (1 − f, f), so each end's share is capped by what is left of that node's allowance; one scale factor, set by the tighter end, applies to the whole line, so the spill keeps its shape along the bank and the split between the ends stays the chainage split. A reach whose upstream node holds and receives nothing therefore spills nothing, however full the downstream node is. In iterative mode bank spill is not carried between passes the way a lid overflow is (§4.2), so there banks only take water back from the surface.

A bank line is placed on the solver's own cells: on a refined mesh the line is walked in steps of a quarter of the smallest cell and every mesh cell it enters is used, with the crest taken from that cell's bed where the bank has no crest of its own.

4.4 Modes

Both modes write the surface results file and return a CoupledSummary with the 1D files, the iterations run, and the total volumes surcharged and captured.

5. File formats

5.1 The sidecar <model>.2d

Plain text beside the .inp, which stays EPA's format untouched. Sections and keywords are case-insensitive, ; starts a comment, paths or names with spaces are double-quoted, a missing file is the default configuration.

[FORMAT]
VERSION      2                        ; absent in files from before version 2

[GRID]
DEM          "C:/data/site dem.asc"   ; .asc or .tif, in the model's map units
CELL         2                        ; resample to this cell size (optional)
WINDOW       xmin ymin xmax ymax      ; restrict to this window (optional)

[ROUGHNESS]
UNIFORM      0.05                     ; or:
GRID         n.asc                    ; a raster of Manning's n
CLASSES      landcover.asc            ; a class raster with a table:
CLASS        11 0.03
CLASS        22 0.1

[BUILDINGS]                           ; optional; footprints block cells (§3)
BLOCK        "buildings.shp"          ; .shp or .geojson, one layer per line

[REFINE]                              ; optional; finer cells near nodes/boxes
CELL         0.5                      ; fine cell size
RADIUS       10                       ; how far the fine cells extend
RATIO        2                        ; largest size ratio between neighbours
NODES        YES                      ; refine around interfaced nodes
BOX          x0 y0 x1 y1              ; and/or inside boxes

[RAIN]
NONE | GAGE RainGage | CONSTANT 25    ; in/hr or mm/hr
SUBCATCHMENTS [FallbackGage]          ; each cell its subcatchment's gage
GRIDDED RATE 300                      ; or DEPTH; seconds each frame lasts
FRAME 0:00 "radar/r0000.asc"          ; start (H:MM[:SS] or s) and raster
FRAME 0:05 "radar/r0005.asc"

[INFILTRATION]
NONE | CONSTANT 5 | HORTON f0 fc k    ; rain units; k in 1/hr
GREEN_AMPT 3.5 0.5 0.26 SUBCATCHMENTS ; suction (in|mm) Ksat (rain units) IMD;
                                      ; SUBCATCHMENTS is optional

[BOUNDARY]
CLOSED | OPEN | FIXED_HEAD 98.25      ; every edge face no segment covers
SEGMENT Tide STAGE TideLevel EDGE WEST          ; optional, any number
SEGMENT Creek INFLOW FILE "q.csv" EDGE NORTH 120 180
SEGMENT Out NORMAL_DEPTH 0.002 LINE 500 0 500 40
SEGMENT Wall CLOSED EDGE SOUTH 0 50

[EXCHANGE]                            ; written only when not the defaults
WEIR         0.6                      ; weir coefficient for interfaces on auto
ORIFICE      0.6                      ; orifice coefficient for interfaces on auto

[RUN]
DURATION     7200                     ; seconds (default: the model's)
OUTPUT_STEP  300                      ; seconds between frames
DRY_DEPTH    0.003                    ; DEM units
COURANT      0.7                      ; 0 < α ≤ 1
MAX_DT       2                        ; seconds (optional)
THREADS      0                        ; 0 = all cores
SCHEME       INERTIAL                 ; or HLL: the full-dynamic scheme (§2.4)
ORDER        2                        ; HLL only: 1 or 2 (§2.4a)

[SOURCES]
;;Name  X  Y  Series
hose    12.5 40 Q1

[NODES]
;;Node  MANHOLE [weir|*] [ORIFICE c|*] [OPEN] | INLET perimeter area [weir|*] [ORIFICE c|*] | SEALED
J1      MANHOLE 0.5 ORIFICE * OPEN
J2      INLET 3.0 0.2 * ORIFICE 0.67
J3      SEALED

[BANKS]
;;Link  LEFT|RIGHT  crest|DEM  coeff|*  x1 y1 x2 y2 ...
C1      RIGHT DEM 0.7 0 0 10 5

[SEALED]
O1

[GROUNDWATER]                          ; optional; written only when set (§2.11)
ON                                     ; or OFF (the default)
OWNER        2D                        ; or SWMM; needed only when a subcatchment is both (§2.11.6)
STEP         60                        ; groundwater step, seconds
GWI          0.002                     ; GWI leakance C, 1/day; 0 = none (§2.11.4)
SEEPAGE      YES                       ; seepage onto cells below the water table
DEFAULT      Loam J1 * 0.001 1.5 0 0 0 0   ; aquifer node ground|* A1 B1 A2 B2 A3 fixedDepth
STATE        "runs/last.2d.gwstate"    ; start from this state (§5.3)

Config::to_text writes exactly this layout and Config::parse reads it; the round trip is tested. In a [NODES] row, * is "the model default". A row with a weir coefficient and no ORIFICE — every row written before the two coefficients were separate — applies that one coefficient to both forms, as it always did; the writer therefore spells out ORIFICE whenever a row sets either coefficient. Sections added later ([EXCHANGE] and the others marked as optional) are written only when they differ from the defaults, so a sidecar that uses none of them is the same file it was.

Versions. [FORMAT] VERSION says which format a file is; a file without it is version 1. Version 2 added SCHEME and ORDER, and writes both on every save. A version-1 file reads as it ran: the local-inertial scheme, and first order should the command line's --scheme hll select the full-dynamic one, so opening an old model changes nothing about how it runs. In a version-2 file a missing ORDER is the current default (second). A file of a newer version than the program reads is refused with a message rather than misread. Older AquaSWMM builds refuse version-2 files (they do not know [FORMAT]), which is the point of the marker.

[GROUNDWATER] is new in this release and optional: a file without it reads exactly as before, with groundwater off, and a configuration with groundwater at its defaults writes no section, so the format version is still 2. An older AquaSWMM refuses a file that has the section (unknown section) rather than running it without groundwater.

5.2 The results file <model>.2d.out

Little-endian binary, written frame by frame so a stopped run still opens.

header   "SS2D"  u32 version  u32 ncols  u32 nrows  f64 x0  f64 y0  f64 cell
         u8 metric  f64 frame_step_s  u32 n_frames  f64 dry_depth
         u32 n_nodes  { u16 len, utf-8 name } × n_nodes
         [version 2] f64 x_edges[ncols+1]  f64 y_edges[nrows+1]
frame    f64 time_s  f32 depth[n]  f32 vx[n]  f32 vy[n]  f32 node[n_nodes]
trailer  f32 max_depth[n]  f32 max_speed[n]  f32 max_hazard[n]  f32 arrival[n]
         u32 n_frames  "2DSS"

Version 1 is a uniform grid, fully described by cell. Version 2 adds the column and row edges for a refined mesh, whose cells are not all one size; x_edges ascends from the left, y_edges descends from the top, and cell is then the smallest cell, for readers that need one number. A uniform run is still written as version 1, so a reader from before refinement keeps opening every file it could already open.

n = ncols × nrows, row-major top row first, x0, y0 the lower-left corner as in an ESRI ASCII grid. Depth is NaN on no-data cells; vx, vy are world axes (x east, y north); node[i] is the average exchange rate at interface i over the frame interval, positive onto the surface, in cubic DEM units per second. The trailer holds the maxima over every step (not just the frames): depth, speed, depth × speed (the hazard index), and the arrival time in seconds (NaN where never wet). The header's n_frames is 0 until the file is sealed; a reader that finds no 2DSS tail counts frames from the file size and computes the maxima from the frames it has. The CLI writes the maximum-depth grid or one frame as .asc (twod max, twod frame). An .asc holds one cell size, so a refined run's grid exports are resampled onto the uniform grid its finest cell implies — nearest cell, no interpolation, since smearing depths across a refinement boundary would invent water the solver never routed.

5.3 The groundwater state <model>.2d.gwstate

Plain text, written at the end of every run with groundwater on, read when a sidecar's STATE names it. It is SWMM's own groundwater hot-start state (gwater_getState: moisture, water-table elevation, lateral flow) per store, plus the two things only the 2D stores have:

; AquaSWMM groundwater state (see manual chapter 20b, groundwater)
[GWSTATE]
VERSION 1
TIME_S 21600                 ; seconds from the run's start when written
UNITS SI                     ; US (ft, cfs) or SI (m, cms)
;;Store  Theta  WaterTable  LateralFlow  Pending  [GWI per pipe ...]
S1 0.287000000 9.412000000 3.1e-4 0e0 1.2e-5 1.1e-5
*DEFAULT 0.250000000 9.000000000 0e0 0e0

Store is the subcatchment (or *DEFAULT); Theta the upper-zone moisture; WaterTable the water-table elevation (so a state survives a change of the aquifer bottom); LateralFlow the store's flow to its node, in cubic length units per second (ft³/s or m³/s); Pending water above the store's capacity still to seep back, cubic length units; then the GWI rate of each of the store's pipes, in [CONDUITS] order. A store absent from the file starts from its aquifer's initial state; rows for stores the run does not have are ignored. The format is deliberately simple so a forecast chain (for instance one driving SWMM hot-start files) can read, edit or write it.

5.4 The groundwater results <model>.2d.gw

Little-endian binary beside <model>.2d.out, one frame per saved 2D frame:

header  "AQGW"  u32 version (1)  u32 n_cells  u32 n_stores
        { u16 len, utf-8 name, f64 ground, f64 bottom, f64 area } × n_stores
        i32 store[n_cells]   (−1: no store)   f32 ground[n_cells]
frame   f64 time_s  { f64 θ, f64 water_table, f64 lateral, f64 gwi } × n_stores

The map's Water table depth layer is ground − water_table per cell (negative where the water table stands above the ground).

6. Validation

The cases below are the integration tests in swmm/tests/twod_validation.rs and swmm/tests/twod_coupling.rs; the thresholds are what the tests assert on every build.

6.1 Lake at rest on a random bed

50 × 50 cells of 2 m, bed uniformly random in 100–103 m (a fixed LCG seed), one no-data cell as an island, filled still to 104 m, 1000 steps, both schemes and both HLL orders; HLL again filled only to 102 m, so a third of the bed stands out of the water as dry islands with shorelines between. Asserted: every cell velocity ≤ 1e-12, depth drift ≤ 1e-12, balance error < 1e-9 %. The local-inertial scheme is exactly balanced because its driving term is the surface slope; HLL is balanced by the hydrostatic reconstruction (§2.4), and at second order by reconstructing the surface rather than the depth together with Audusse's centred source term (§2.4a). Measured (largest speed, largest depth drift): local inertial 0 and 0; HLL first order 8e-15 m/s and 9e-15 m (1e-13 and 2e-14 with the dry islands); second order 2e-13 and 1e-14 (3e-13 and 1e-14) — round-off.

6.2 Steady uniform flow (Manning)

A 200 m plane of 1 m cells at S = 0.01, n = 0.03, q = 0.1 m²/s fed at the top edge, open downstream boundary, run to 2000 s. Manning's normal depth is h_n = (q n / (k √S))^{3/5} = 0.03^{0.6} = 0.12203 m. Asserted: depth at 80, 100, 120 and 150 m within 1 % of h_n for both schemes; unit discharge within 1 % of q (local inertial), 2 % (HLL first order, whose numerical diffusion costs 1.49 % here) and 0.1 % (HLL second order: measured 0.003 %); no cross-slope difference (< 1e-6). The HLL depth errors are 0.99 % (first order) and 0.64 % (second); what remains at second order is the semi-implicit friction, whose factor is evaluated from the momentum before friction and so carries a first-order-in-Δt bias (§7).

6.3 Dam break on a dry frictionless bed (Ritter)

Ritter (1892), in the form given by Stoker (1957, Water Waves, §10.8): with c₀ = √(g h₀) the depth at time t is

h = h₀                          x < −c₀ t
h = (2c₀ − x/t)² / (9g)         −c₀ t ≤ x ≤ 2c₀ t
h = 0                           x > 2c₀ t

1200 × 3 cells of 1 m, h₀ = 1 m behind a dam at 600 m, n = 0, checked at t = 40 s for −0.8 c₀t ≤ x ≤ 1.5 c₀t, i.e. away from the negative wave and the dry tip where any first-order scheme smears the solution, and as an L1 error over the whole channel, fronts included (normalised by h₀ and the wave's extent 3 c₀ t).

worst away from the fronts L1, whole wave
HLL first order 0.894 % of h₀ 0.00553
HLL second order 0.101 % of h₀ 0.00086

Asserted: first order within 5 % of h₀ in the range; second order within 0.3 %, under a fifth of the first-order worst error, and an L1 error under a quarter of first order's (measured: 8.9 and 6.4 times smaller); the still water and the dry bed far from the dam untouched (1e-9); balance < 1e-6 %; no negative depth. This case runs the HLL scheme. It is not a test of the local-inertial scheme: Ritter's solution has Froude numbers up to 2 at the tip and depends on the convective term that scheme omits, which is exactly the regime Bates et al. (2010) exclude; on this case the local-inertial profile is wrong by tens of percent, as expected, so a dam break must be run with Scheme::Hll.

6.3a Order of accuracy on a smooth wave

The convergence test of the second-order scheme (§2.4a) needs an exact solution that stays smooth. On a flat frictionless bed a state with u − 2c = −2c₀ everywhere (c = √(g h)) is a right-going simple wave: c is constant along the characteristics dx/dt = u + c = 3c − 2c₀, so

c(x, t) = c(ξ, 0)  where  x = ξ + (3 c(ξ, 0) − 2c₀) t

(Stoker 1957, §10.2; Toro 2001, §3.1), solved here for ξ by bisection. The initial celerity is a Gaussian hump, c = c₀ + 0.2 exp(−((x − 40)/10)²) on 1 m of still water, u = 2(c − c₀), in a 120 m channel. The wave steepens as it travels and would break at about 19.4 s; it is compared at 8 s, on 120, 240, 480 and 960 cells (the time step follows the cell, so the observed order is that of the whole space–time scheme).

cells L1 depth error, first order second order
120 3.70e-1 4.06e-2
240 2.03e-1 1.09e-2
480 1.07e-1 2.86e-3
960 5.54e-2 7.37e-4
observed order 0.86, 0.92, 0.95 1.89, 1.93, 1.96

Asserted: first order between 0.7 and 1.3; second order above 1.5 on every pair and above 1.8 on the finest, and more accurate on every grid. The second-order error on 120 cells is below the first-order error on 960.

6.3b A moving hydraulic jump: dam break onto still water

Stoker's (1957, §10.7) wet-bed dam break: h_L = 1 m, h_R = 0.1 m, flat and frictionless. The exact solution is a rarefaction, a uniform middle state and a bore (a moving hydraulic jump) into the still water. The middle depth comes from u_m = 2(c_L − c_m) across the rarefaction and the Rankine–Hugoniot conditions u_m = S (1 − h_R/h_m), S² = g h_m (h_m + h_R)/(2 h_R) across the jump, solved by bisection. 800 cells of 0.25 m, t = 10 s.

cells across the jump (10–90 % of its height) L1 error
HLL first order 3 0.00731
HLL second order 2 0.00138

Asserted: at most 3 cells across the jump at second order and fewer than at first order, an L1 error under 0.3 of first order's (measured 5.3 times smaller), and the balance to 1e-9 %. A shock steepens itself, so even first order keeps it to a few cells; the gain is in the rarefaction, the plateau's level and the jump's position.

6.3c Wet/dry: Thacker's oscillation in a paraboloid

Thacker's (1981) exact solution for water sloshing in a frictionless paraboloidal bowl (the curved-surface, axisymmetric case; also in the SWASHES collection, Delestre et al. 2013): bed z = h₀ (r²/a² − 1), and with ω = √(8 g h₀)/a, A = (a² − r₀²)/(a² + r₀²),

η = h₀ [ √(1−A²)/(1 − A cos ωt) − 1 − (r²/a²) ((1−A²)/(1 − A cos ωt)² − 1) ]
h = max(η − z, 0)

a lens of water whose shoreline moves in and out every period. a = 1 m, r₀ = 0.8 m, h₀ = 0.1 m on a 4 m square of 4 cm cells (100 × 100), dry depth 1e-5 m, one period (2.24 s). Every step is checked for a negative depth.

L1 depth error / volume after one period balance
HLL first order 0.0998 4e-13 %
HLL second order 0.0291 3e-13 %

Asserted: no negative depth in any step, the balance to 1e-9 % at both orders, and the second-order error under half the first-order one (measured 3.4 times smaller). Most of what remains is at the moving shoreline, where both orders are first order by construction (§2.4a).

6.4 Mass conservation in a closed basin

A 60 × 60 parabolic bowl, closed edges, a triangular source hydrograph peaking at 2 m³/s over 600 s (600 m³ in all), run to 900 s. Asserted: inflow booked = 600 m³ to 1e-9, balance error < 0.1 % (achieved: round-off), stored volume within 0.1 % of 600 m³, zero boundary outflow.

6.5 Symmetry

A constant 0.5 m³/s source at the centre of a flat 41 × 41 plate, 120 s. Asserted: the depth field is mirror-symmetric in x, in y and under the diagonal swap to 1e-9, and the water has spread past 5 cells.

6.6 Time step

A 10 × 10 plate filled to 0.25 m then 4 m. Asserted: Δt = 0.7 Δx / √(2 g h) to 1e-12, the ratio of the two steps is exactly 4 (√16), and MAX_DT caps it.

A second case guards the two-dimensional limit of §2.5: a 25 × 25 grid of 20 m cells shaped like one depression of the EA Test 2 egg box (0.5 m deep, n = 0.03), started with a 0.1 m tilt on the water surface and run for 10 h at the default α = 0.7. Asserted: every wet cell's level within 1 mm of the mean at the end. At the old step (equivalent to α = 0.99 here) the residual is 0.09 m, a standing checkerboard.

6.7 Rain, infiltration and a fixed-head edge

A 30 × 30 plate under 36 mm/h with 7.2 mm/h infiltration for 600 s, then one edge opened to a head below the bed. Asserted: rain and infiltration volumes booked exactly (1e-9), water leaves through the edge, balance < 1e-9 % throughout.

6.7c Boundary segments

swmm/tests/twod_boundaries.rs, every case under the local-inertial scheme and HLL at first and second order (LI, HLL1, HLL2), on a uniform mesh and on a graded one (a band of quarter-size cells across the middle third, graded by halves):

swmm/tests/twod_setup.rs resolves segments from a sidecar against the Site Drainage model: a [TIMESERIES] stage, an inflow file (H:MM and decimal-hour times), a constant inflow converted from CFS, faces by side range and by line, the overlap and empty-segment warnings, and errors that name the segment (a missing series, a zero slope); a 900 s run books each inflow segment's volume exactly. It also reads a version-1 sidecar as the writer at 42664ed wrote it and a version-2 one from 09cd840, writing the latter back byte for byte. Compared bit for bit with the solver before segments (a 23 × 17 basin with a source, CLOSED, OPEN and two FIXED_HEAD levels, LI and HLL at both orders), every depth is identical.

6.7a Rain that varies over the grid

6.7b Green–Ampt infiltration

A flat closed 6 × 6 plate of 5 m cells under 50 mm/h on a soil with Ks = 10 mm/h, ψ = 110 mm, θd = 0.3, against the hand calculation of Rossman & Huber (2016, ch. 4) with the Mein & Larson (1973) ponding time: all rain infiltrates until F_s = Ks C / (i − Ks) = 8.25 mm at t_s = F_s / i = 594 s, after which F − F_s − C ln[(F + C)/(F_s + C)] = Ks (t − t_s). Asserted: before ponding F = i t to 1e-12; at 1200, 2400 and 3600 s every cell's F is within 0.2 % of the hand value and the ponded depth is i t − F; the rate over the next minute is Ks (1 + C/F); infiltration booked equals F times the plate area and the balance is below 1e-9 % throughout. Running the same storm with steps capped at 2 s and at 60 s gives the same F at 2400 s within 0.2 %. The capacity solver satisfies the integrated equation to 1e-12 relative over C from 0 to 0.5 m, F₁ from 0 to 5 m and Δt from 0.01 s to an hour, is never below Ks Δt, and from a dry start tends to √(2 Ks C Δt). On the graded mesh of §6.10 with Ks varying cell to cell, the infiltrated volume booked equals Σ F · area and the balance closes to 1e-9 %. A setup test checks that cells take their own subcatchment's [INFILTRATION] row, and that a subcatchment on Horton and the ground outside every outline take the sidecar's values.

6.8 Interfaces

6.9 Performance

300 × 300 cells with 1 m of water over a third of the grid, 60 s of flow (several hundred CFL steps): asserted under 5 s in a release build (ignored in debug builds, where it is about ten times slower). This runs the default scheme, local inertial.

The same case with HLL compares the two orders (an #[ignore]d test, so it does not run beside the 5 s one and slow it). Release build, this machine, with other builds running on it:

per step 1 thread 8 threads all 24 cores
local inertial (367 steps) 5.5 ms 3.7 ms 4.1 ms
HLL first order (646 steps) 9.2–11.8 ms 3.4–4.8 ms 4.7 ms
HLL second order (654 steps) 20.7–31.2 ms 8.2–13.1 ms 10.5 ms

(Ranges are repeat runs at different loads.) Second order costs 2.2–2.7 times first order per step, at the same number of steps; the 300 × 300 HLL run takes 3.0–3.3 s at first order and 6.9–7.7 s at second on all cores. A coupled run (the EPA Site Drainage sample with inlets at every junction, 3 in/hr on a DEM through the nodes, iterative, 2 passes, 30 min) on 205 × 117 cells of 5 ft: the 2D-only surface took 19.9 s (HLL first order) and 37.2 s (second), 1.9 times; the coupled run's surface 19.1 s and 52.3 s with the engine's runs in between on a loaded machine. On 54 × 32 cells of 20 ft: 1.05 s and 1.96 s, and 2.3 s and 4.1 s wall for the coupled run. The second-order surface captured 20,903 ft³ through the inlets against 19,437 ft³ at first order (20 ft cells), with the balance at round-off in all of them.

6.10 A graded mesh

The cases that would break first if the variable-cell arithmetic were wrong, re-run on graded meshes.

The first two use a 10 × 10 grid of 4 m refined to 1 m around a point, graded at a ratio of 2 — the shape a refinement actually leaves behind.

The resolver is tested separately (swmm/tests/twod_setup.rs): a 50 ft grid refined to 12.5 ft around the interfaced nodes keeps the DEM's extent and area exactly, holds every cell between the target and the base size, keeps neighbouring columns within the grading ratio, puts every node in a fine cell with real terrain under it, and gives the base grid back cell for cell when refinement is switched off.

6.11 EA 2D benchmarks

The UK Environment Agency's benchmark of 2D packages (Néelz & Pender 2013, report SC120002) gives eight tests with prescribed geometry, roughness, boundary data and output points. Nineteen packages ran them, among them TUFLOW (classic, FV and GPU), MIKE FLOOD, InfoWorks ICM, ISIS 2D, SOBEK, JFLOW+, XPSTORM, ANUGA, LISFLOOD-FP and Ceasg. The cases are in swmm/tests/twod_ea_benchmarks.rs. Their inputs come from the authors' dataset deposit (Zenodo record 15256842); scripts/fetch-ea-benchmarks.ps1 downloads it, checks each zip's MD5 and unpacks it into target/ea2d.

How the comparison is made. The report publishes no tables of results, only plots of every package's time series at each output point. Every reference range below was read off those figures by sampling the plotted curves' pixels at a given time. The tests quote each figure and page, and state each axis's reading precision (about four pixels: 0.003 m on most level axes, 0.02–0.05 m on the taller ones). A range is the spread of the packages the report treats as a group. A package the report itself sets apart is left out, and the test comment says which and why: for example, the 0- and 2-term models in Test 1, or Flowroute run at n = 0.03 in Test 3. "In range" means inside that spread, widened by the reading precision. Where a result falls outside, the test's margin is set to the observed excess, so a regression still fails, and the excess is reported below rather than hidden.

Inputs as specified, and where they could not be.

Results. LI is the local-inertial scheme (§2.2); HLL 1st and HLL 2nd are the full-dynamic one (§2.4) at first and second order (§2.4a, the default). "Out" or a signed excess marks a value outside the spread. The first-order column was produced before the second-order option existed; AQUASWMM_EA2D_HLL_ORDER=1 (or 2) runs the HLL cases at that order.

Test Quantity (report figure, page) Reference spread LI HLL 1st HLL 2nd
1 Flat pond, level BC (Fig 4.3, p. 20) P1 level at 2.5 h 10.165–10.255 m 10.187 10.219 10.171
time to 10.34 m 3.03–3.91 h 3.55 3.32 3.67
P1/P2 level at 13 h / at 19.9 h 10.262–10.277 / 10.249–10.263 10.274 / 10.256 10.267 / 10.251 10.268 / 10.251
2 Depressions (Fig 4.5, pp. 24–29) final levels in 12 depressions; 5 dry per depression all in all in (P11 8.548 vs 8.515–8.545, Ceasg 8.572) all in but P12 8.900 vs 8.858–8.895 (+5 mm); P11 8.557
final volume (Table 4.3) 97,200 m³ 97,200 97,200 97,200
3 Obstruction (Fig 4.8, pp. 37–38) P1 level at 300 s 9.976–9.990 m 9.928 (out) 9.986 9.979 at 5 m
P1 peak speed / time 1.6–2.2 m/s / 65–95 s 1.68 / 94 s 1.79 / 76 s 2.16 / 72 s at 5 m (2.35 at 2.5 m, +0.15)
P2 level at 300 s 9.803–9.827 m 9.928 (out) 9.791 at 5 m (−12 mm) / 9.809 at 2.5 m 9.812 at 5 m / 9.821 at 2.5 m
P2 reaches 9.765 m 148–197 s 158 s 202 s at 5 m / 160 s at 2.5 m 156 s at 5 m / 142 s at 2.5 m (−6 s)
4 Flat plain (Figs 4.10–4.11, pp. 42–45) P1 level at 121 / 199 min 0.313–0.331 / 0.315–0.333 m 0.306 / 0.309 (−7 mm) 0.319 / 0.322 0.320 / 0.323
P5 level at 199 / 251 min 0.165–0.178 / 0.172–0.183 m 0.164 / 0.170 (−2 mm) 0.171 / 0.177 0.175 / 0.181
P1 / P5 / P6 speed at 199 min 0.332–0.509 / 0.079–0.089 / 0.084–0.099 m/s 0.375 / 0.084 / 0.097 0.428 / 0.087 / 0.084 0.376 / 0.089 / 0.086
arrival P5 / P6 (2 cm deep) 50.0–63.2 / 56.1–66.7 min 64.0 / 61.0 min (P6 first) 59.0 / 62.7 min 61.7 / 65.7 min
5 Valley (Figs 4.16–4.17, pp. 53–56) P1 / P3 / P7 peak level 174.53–175.06 / 155.64–156.00 / 155.50–155.83 m 174.63 / 155.81 / 155.64 174.76 / 155.99 / 155.75 174.85 / 155.91 / 155.71
P5 level at 29.9 h 148.51–148.85 m 148.51 148.58 148.55
P5 reaches 145.5 m 2.77–3.48 h (ANUGA, UIM, LISFLOOD-FP ~4 h) 3.50 3.57 (+5 min) 3.52 (+2 min)
P4 / P7 peak speed (§4.6.4 text) 1.5–2.2 / 1.2–1.6 m/s 1.87 / 1.67 (+0.07) 1.69 / 1.45 1.86 / 1.40
P1 peak speed 1.9–2.7 m/s 2.76 (+0.06) 1.98 2.27
6A Dam break, lab (Fig 4.25, p. 71) G6 at 9.9 / 19.9 / 29.9 s 0.268–0.290 / 0.203–0.225 / 0.156–0.178 m — 0.288 / 0.223 / 0.176 0.283 / 0.219 / 0.172
6B Dam break, ×20 (Figs 4.27–4.28, pp. 77–78) G6 at 51 / 101 / 199 / 301 / 590 s 5.25–5.61 / 3.91–4.32 / 2.31–2.67 / 1.75–1.91 / 1.27–1.35 m — 5.60 / 4.30 / 2.62 / 1.87 / 1.30 5.56 / 4.26 / 2.61 / 1.82 / 1.30
G2 / G4 / G5 peak level 2.80–3.25 / 2.05–2.53 / 1.90–2.43 m — 2.96 / 2.14 / 1.94 2.95 / 2.24 / 1.90
G2 / G4 peak speed 4.7–5.4 / 3.9–4.9 m/s — 5.18 / 4.63 5.02 / 4.59
8A Urban rain + inflow (Figs 4.42–4.43, pp. 101–102) P1 peaks (rain / inflow) 27.600–27.640 / 27.655–27.695 m 27.612 / 27.695 27.628 / 27.675 27.629 / 27.672
P2 peaks 28.700–28.745 / 28.760–28.820 m 28.723 / 28.778 28.740 / 28.829 (+9 mm) 28.740 / 28.834 (+14 mm)
P2 peak speeds 0.48–0.82 / 0.61–0.92 m/s 0.64 / 0.89 0.53 / 0.80 0.52 / 0.80
P3 (pond) level at 150 min 24.340–24.400 m 24.337 (−3 mm) 24.368 24.369

The tests compare more points than the table shows: 16 depressions in Test 2, four points and their arrival times in Test 4, G2/G4/G5 at three times in Test 6B, and others. The mass balance closes to round-off (< 1e-10 %) in every run.

Reading the results.

Test 8B: a surcharging culvert, coupled. This is the one test of the 1D-2D coupling (§4), and the one with a table of results: Table 4.11 gives each package's total manhole flow volume. The surface is Test 8A's (2 m, n = 0.02 on roads and pavements and 0.05 elsewhere, closed edges, dry start, 5 h), now with the building footprints blocking their cells (8466 cells) and no rain. Under it runs a 1400 mm culvert at n = 0.017, with inverts from the dataset's geometry table and a free outfall. It carries a 1.6 m³/s baseflow, which is the initial flow in every conduit, and the hydrograph Test8B-pipe-inflow.csv, which peaks at 6.28 m³/s at 2 h. The water surcharges through a 1 m² vertical manhole 467 m down the culvert, with its rim at 31.46 m.

The test builds the culvert as an EPA SWMM .inp and the surface as its .2d sidecar, both from the published files and both at test time. The manhole MH is the only node with coordinates, so it is the only node on the surface. The other six junctions are sealed joints with a 100 m surcharge depth. MIN_SURFAREA 1 m² gives the manhole its plan area. The routing step is 1 s; 0.5 s and 5 s change the flood volume by under 2 m³. Two discrepancies in the inputs:

The culvert runs with the slot surcharge method, not EPA's default EXTRAN. Under EXTRAN a pressurised junction has no storage and its head is set from the flows, so a withdrawal from it does not come out of the network's volume. With the Opening lid, EXTRAN reported −16.5 % flow routing continuity at a 1 s exchange (−3.0 % at 10 s), and the manhole flow changed sign every few minutes. The slot method closes to +0.05 %. With the Spec lid both methods are clean, but they flood different volumes: EXTRAN gives 5951 m³ at a 2.288 m³/s peak, and the slot method gives 5763 m³ at 2.267 m³/s.

Results: tight coupling (1 s exchange) unless marked. The reference ranges cover the full shallow-water packages plus UIM, which the report finds "generally in agreement" with them (§4.10.3). ISIS Fast and ISIS Fast Dynamic are left out. Each figure was read like the others: four pixels are 0.016 m³/s, 0.004 / 0.002 / 0.006 m at points 1 / 2 / 3, and 1.1 min. Manhole flows are 60 s means. Signed values are excesses beyond the range.

Quantity (report table/figure, page) Reference Spec, LI Spec, LI iterative Spec, HLL 1st / 2nd Opening, LI Opening, HLL 1st / 2nd
net manhole volume (Table 4.11, p. 110) 4954–5840 m³ 5764 5762 5764 / 5764 4618 (−336) 4557 / 4574 (−380)
peak manhole flow (Fig 4.46, p. 110) 2.00–2.23 m³/s 2.267 (+0.037) 2.267 2.267 1.909 (−0.09) 1.880 / 1.885 (−0.11)
its time 118.5–121 min 119.5 119.5 119.5 119.5 118.5
outflow starts / ends 81–88 / 155–157 min 91.5 (+3.5) / 154.5 91.5 / 154.5 91.5 / 154.5 91.5 / 152.5 (−2.5) 91.5 / 153.5
flow back into the pipe 0 to −0.35 m³/s none none none −0.32 m³/s, 130 m³ −0.31 m³/s, 92 / 93 m³
P1 peak level (Fig 4.47, p. 111) 27.612–27.646 m 27.659 (+13 mm) 27.658 27.641 / 27.651 (+5 mm) 27.650 (+4 mm) 27.625 / 27.636
P1 level at 299 min 27.517–27.555 m 27.558 27.557 27.557 / 27.557 27.557 27.556 / 27.556
P1 arrival (5 cm) 107–112 min 110 110 112 / 112 112 114 / 114
P2 peak level (Fig 4.47, p. 111) 28.724–28.753 m 28.739 28.743 28.767 (+14 mm) / 28.763 (+10 mm) 28.730 28.751 / 28.748
P3 level at 149 min (Fig 4.47, p. 112) 24.055–24.145 m 24.104 24.110 24.098 / 24.102 24.038 (−17 mm) 24.035 / 24.042 (−13 mm)
P3 level at 299 min 24.220–24.279 m 24.261 24.260 24.265 / 24.267 24.191 (−29 mm) 24.191 / 24.193 (−27 mm)
P3 arrival (5 cm) 120.5–131 min 128 127 128 / 127 131 130 / 130
P3 peak speed (Fig 4.48, p. 112) 0.34–0.75 m/s 0.40 0.38 0.58 / 0.58 0.40 0.59 / 0.74
engine flow-routing continuity +0.048 % +0.048 % +0.048 % +0.049 % +0.049 %
surface mass error 1e-13 % 2e-14 % 3e-13 % 3e-13 % 2e-13 %
coupling: surface − engine +1.0 m³ of 5763 flooded −1.2 m³ +1.0 m³ −1.9 m³ withdrawn, −6.5 m³ returned −0.9 / −1.0 m³, −5.9 / −6.4 m³

With the Opening lid, iterative coupling (local-inertial, four passes) was run once to see how it behaves; there is no test for it.

Use tight coupling for a two-way manhole (§7).

Reading Test 8B.

Two bugs Test 8B found, both in the coupler:

A bug the benchmarks found. At the old local-inertial time step, Δt = α Δx / √(g h), Test 2 never settled. Depression 4 held a standing ±0.3 m checkerboard for all 48 h, and ended at 9.625 m against 9.534–9.545 m. §2.5 explains the cause and the fix: the two-dimensional stability limit is a factor √2 tighter. A regression case is in §6.6. After the fix, every depression is inside the spread.

Two performance fixes came from the same runs (§2.10):

Run times. Release build, this machine, with other runs sharing it:

Test LI HLL
Test 2 25 s 60 s
Test 3, 5 m 0.1 s 0.2 s
Test 4 3 min 4 min
Test 5 4 min 3.5 min
Test 6A — 50 s
Test 6B — 90 s
Test 8A 3–9 min 6 min
Test 8B, coupled tight 5–10 min 7 min (1st order), 12–18 min (2nd)
Test 8B, coupled iterative 6–8 min (1 pass) —

The report's packages took from seconds (GPU) to hours on the same tests. Run times are not compared further, because the hardware differs.

Not run, and why.

The long cases are #[ignore]d, and they skip with a message when the data has not been fetched. The Test 8B cases also skip without an EPA SWMM engine and, for tight coupling, without the 32-bit bridge (cargo build --release -p aquaswmm-swmm-bridge --target i686-pc-windows-msvc; a release test binary finds it in target/i686-pc-windows-msvc/release or through AQUASWMM_SWMM_BRIDGE). To run them:

pwsh scripts/fetch-ea-benchmarks.ps1
cargo test -p aquaswmm-swmm --release --test twod_ea_benchmarks -- --include-ignored --nocapture

Test 3 at 5 m (both HLL orders), and the check of the embedded profiles, run on every build. Test 1 and Test 3 at 2.5 m run on every release build. At second order the HLL cases take about twice as long as the run times above.

6.12 Groundwater

swmm/tests/twod_groundwater.rs.

Against EPA SWMM. A 10-acre, fully pervious Horton subcatchment over an aquifer (φ 0.46, WP 0.13, FC 0.23, K_s 2 in/hr, HCO 10, HCU 15, ET_u 0.35, lower ET depth 4 ft, deep percolation 0.002 in/hr, bottom 0, water table 5 ft, θ₀ 0.24, ground 10 ft, node invert 2 ft), three days of 1-minute steps with two storms and evaporation, is run in EPA SWMM 5.2.4. The store is then driven with SWMM's own reported infiltration and surface evaporation, step by step, and compared with SWMM's reported water table, upper-zone moisture and groundwater flow at every step (4 320 steps):

Case max |Δ water table| max |Δθ| max |Δ GW flow| (peak)
A1 0.0015, B1 1.6; deep percolation; ET 0.15 in/day 2.8e-5 ft 6.8e-8 3.2e-6 cfs (0.117 cfs)
A1 0.004, B1 1.2, A2 0.001, B2 1, A3 2e-5, fixed depth 1.5 ft, plus [GWF] LATERAL 0.0005*STEP(Hgw-Hcb)*(Hgw-Hcb); ET 0.1 in/day 4.8e-5 ft 1.1e-7 5.6e-6 cfs (0.193 cfs)

The tests assert 1e-4 ft, 1e-5 and 1e-4 of the peak flow. What is left is the projection (§2.11.5) — SWMM's own states are not quite conservative — and SWMM's single-precision output.

Saturation excess. 100 m × 100 m flat at 10 m, the water table 5 cm below the ground, 30 mm/hr of rain for 4 h on a 30 mm/hr constant infiltration capacity. Without groundwater all 1 200 m³ infiltrate; with it, 82.1 m³ — the 75 m³ of room above the water table plus the 7.2 m³ that drained to the node meanwhile — and the rest stays on the surface. Groundwater balance 1.5e-12 %, surface 5.7e-13 %, and the infiltration the store received equals what the surface lost to it.

Seepage. A grid sloping 9.5 → 10.5 m over a store with its water table at 9.9 m: 280.6 m³ seeps onto the cells below 9.9 m over 4 h; the high cells stay dry; the Water table depth layer is negative at the low corner and positive at the high one; both balances below 5e-13 %.

GWI and lateral flow. The same site with the water table 1 m below the ground, two 0.6 m pipes 2.5 m down and GWI 0.05: driving the simulator interval by interval as a coupler does, the network is given 9.895 657 m³ and the stores booked 4.454 772 m³ lateral + 5.440 885 m³ GWI — equal to 1e-12. Coupled through the engine: tight, the lateral inflows set equal the books to 1e-9 and the engine's own reported lateral inflow at the three nodes integrates to 9.875 m³ (its report samples a step-shaped input once a minute); iterative (two passes), the written series carry exactly the books (9.895 380 m³) and the engine's record integrates to 9.885 m³.

With boundary segments. A site sloping 10.2 → 9.82 m over a store whose water table is at 9.9 m, a 0.2 m³/s inflow segment on the west edge and a normal-depth segment on the east, 20 mm/hr infiltration and GWI: 2 880 m³ in and 2 518.8 m³ out through the segments, 471.3 m³ recharged, 312.9 m³ seeping back on the low ground, 7.1 m³ lateral, 2.7 m³ GWI. The segment books equal the surface's boundary terms, and the combined surface-plus-groundwater residual is 4e-11 m³ (surface 3e-13 %, groundwater 5e-13 %).

Ownership and state. A subcatchment with SWMM groundwater under 2D cells stops the run until the owner is chosen; with OWNER 2D it has its own store and its row is gone from the engine's copy, with OWNER SWMM its cells (half the grid) have none. A 2-hour run's state, read by the next run, starts its water table where the first ended.

7. Limitations

8. References