The shallow water equations (SWE) describe the depth-integrated flow of a thin layer of fluid under gravity when the horizontal scale greatly exceeds the depth. They govern tidal and flood waves, storm surges, tsunamis racing across ocean basins, bore formation in rivers, and geostrophic eddies in the atmosphere. The same conservation structure appears in gas dynamics — with depth h playing the role of density ρ and the shallow-water Froude number Fr = u/csw analogous to the Mach number.

1. Derivation and Conservation Form

Starting from the 3-D incompressible Navier-Stokes equations, depth-averaging under the hydrostatic and kinematic surface conditions yields the 1-D SWE in conservation form:

∂U/∂t + ∂F(U)/∂x = S(U, x) U = [h, hu]ᵀ (conserved variables) F(U) = [hu, hu² + ½gh²]ᵀ (flux vector) S = [0, −gh ∂b/∂x]ᵀ (topography source) g = 9.81 m/s², b(x) = bed elevation, η = h + b = surface elevation

The eigenvalues of the Jacobian ∂F/∂U are the wave speeds:

λ₁ = u − c, λ₂ = u + c, c = √(gh) (shallow-water celerity) Froude number: Fr = |u|/c Fr < 1: subcritical (tranquil) flow Fr = 1: critical (hydraulic jump transition) Fr > 1: supercritical (shooting) flow

2. Roe's Approximate Riemann Solver

The cell-centred finite volume method updates each cell average Ui from interface fluxes. At each interface i+½, an exact Riemann problem is expensive. Philip Roe (1981) proposed a linearised solver: replace the nonlinear problem with a linear one having a specially constructed Roe average state Û that satisfies exact conservation across shocks:

Roe-averaged velocity and celerity: û = (√h_L · u_L + √h_R · u_R) / (√h_L + √h_R) ĉ = √( g(h_L + h_R)/2 ) Eigenvalues: λ̂₁ = û − ĉ, λ̂₂ = û + ĉ Wave strengths (characteristic decomposition): α₁ = Δ(hu)/2ĉ − Δh/2 · (û/ĉ − 1) α₂ = Δh − α₁ − α₃ … (full 2-wave split) Interface flux: F_{i+½} = ½(F_L + F_R) − ½ Σ |λ̂ₖ| αₖ rₖ

The term ½ Σ |λ̂ₖ| αₖ rₖ is the numerical diffusion — it damps oscillations but smears shocks over a few cells. High-order MUSCL reconstruction before the Roe solve recovers second-order accuracy while keeping the scheme conservative.

3. MUSCL Reconstruction and Limiters

First-order Godunov is stable but diffusive. MUSCL (Monotone Upstream-Centred Schemes for Conservation Laws, van Leer 1979) extends to second order by reconstructing a linear profile within each cell:

U_L at i+½ : U_i + ½ φ(r) (U_i − U_{i-1}) U_R at i+½ : U_{i+1} − ½ φ(1/r) (U_{i+2} − U_{i+1}) r = (U_i − U_{i-1}) / (U_{i+1} − U_i) (slope ratio) Slope limiters φ(r): minmod: max(0, min(1, r)) van Leer: (r + |r|)/(1 + |r|) superbee: max(0, min(2r, 1), min(r, 2))

The minmod limiter is the most diffusive but guaranteed monotone. Van Leer is a good default balancing accuracy and stability. Superbee is the least diffusive but can steepen smooth extrema.

4. CFL Stability Condition

Explicit time integration requires the Courant-Friedrichs-Lewy condition. For the SWE, the fastest characteristic wave speed is |u| + c:

Δt ≤ CFL · Δx / max_i(|uᵢ| + cᵢ) Safe CFL ≤ 0.9 for first-order Godunov Safe CFL ≤ 0.5 for MUSCL with Runge-Kutta-2 CFL > 1: information travels more than one cell per step → instability

5. Wetting and Drying

Flood simulations frequently involve dry beds (h = 0). Naive Roe solvers produce negative depths and blow up. Standard fixes:

Technique Description
Depth threshold if h < h_dry (e.g. 10⁻³ m), set u = 0 and treat as dry wall boundary
Positivity-preserving limiter After each step, clip h ≥ 0; correct hu accordingly
Hydrostatic reconstruction Audusse et al. 2004 — reconstruct h at interface accounting for bed slope to guarantee h_L, h_R ≥ 0
Interface celerity Replace ĉ with max(ĉ, √(g · max(h_L, h_R) / 2)) in Roe to avoid division by zero

6. Source Terms: Topography and Friction

The well-balanced property ensures that a lake at rest (u = 0, η = const) remains exactly stationary despite non-zero bed gradients. A naive operator-splitting fails this: the flux difference and the source term must cancel exactly. The hydrostatic reconstruction of Audusse et al. achieves this automatically.

Bottom friction is modelled by a Manning-type term added to the momentum source:

S_friction = −g n² |u| u / h^(4/3) n = Manning roughness coefficient Open water: n ≈ 0.013–0.025 m^(−1/3) s Floodplain grass: n ≈ 0.030–0.050 Dense vegetation: n ≈ 0.060–0.120 Implicit treatment of friction avoids timestep restriction for shallow h: uⁿ⁺¹ = (uⁿ + Δt·advection) / (1 + Δt·g n² |uⁿ| h^(−4/3))

7. 2-D Extension and Coriolis

Extending to 2-D adds a y-momentum equation and, for geophysical flows, the Coriolis force:

∂h/∂t + ∂(hu)/∂x + ∂(hv)/∂y = 0 ∂(hu)/∂t + ∂(hu²+½gh²)/∂x + ∂(huv)/∂y = −gh∂b/∂x + fhv ∂(hv)/∂t + ∂(huv)/∂x + ∂(hv²+½gh²)/∂y = −gh∂b/∂y − fhu f = 2Ω sin(φ) (Coriolis parameter, φ = latitude, Ω = 7.27×10⁻⁵ rad/s) f·c/g ≪ 1: gravity waves dominate (coastal/fluvial) f·L/c ≳ 1: Coriolis-dominated (oceanic mesoscale eddies, Rossby waves)
Phenomenon Dominant balance Typical scale
Dam break Inertia + pressure L ~ 10–1000 m, Fr ~ 1–3
Tidal bore Inertia + pressure + friction L ~ 10 km, Fr ~ 1
Ocean tsunami Pressure + dispersion correction L ~ 1000 km, c ~ 200 m/s
Geostrophic eddy Coriolis + pressure L ~ 100 km, Ro ≪ 1
Hurricane storm surge Wind stress + Coriolis + friction L ~ 500 km, h ~ 1–5 m surge

8. Interactive: 1-D Dam Break

This simulation solves the 1-D SWE on a 400-cell grid using first-order Roe + CFL-adaptive timestepping. Adjust the upstream depth ratio and watch the rarefaction wave propagate left and the bore propagate right. The exact analytical solution (Stoker 1957) is overlaid in amber.

t = 0.00 s