10.4 Shallow water equations

Derivation of shallow water equations

Shallow water equations model fluid flow in rivers, coastal areas, and atmospheric systems. They simplify the full 3D fluid dynamics by assuming vertical motion is negligible compared to horizontal flow. The result is a depth-averaged system that captures waves, shocks, and other key features while remaining computationally tractable. These equations are the backbone of flood prediction, tsunami modeling, and large-scale weather forecasting.

Conservation laws in shallow water flows

The shallow water equations rest on two fundamental conservation principles:

• Mass conservation requires that the total fluid mass is preserved, producing the continuity equation.
• Momentum conservation balances the forces acting on the fluid (pressure gradients, gravity, friction), producing the momentum equations.

Together, these two laws fully determine how the water depth and velocity evolve over time.

Assumptions and simplifications

Four key assumptions let you reduce the 3D Navier-Stokes equations to the 2D shallow water system:

1. Small aspect ratio: The vertical length scale is much smaller than the horizontal length scale ($H \ll L$). This is what makes the flow "shallow."

2. Incompressibility: Fluid density $\rho$ stays constant throughout the flow.

3. Hydrostatic pressure: Vertical acceleration is negligible compared to gravitational acceleration, so pressure at any depth is simply the weight of the water above it: $p = \rho g (h - z)$.

4. Slowly varying bathymetry: The bottom topography changes gradually enough that depth-averaging remains valid.

When any of these assumptions breaks down, the standard shallow water equations lose accuracy and you need extensions (covered later in this guide).

Depth-averaged continuity equation

You obtain the continuity equation by integrating the 3D incompressibility condition over the fluid depth. The result relates changes in water depth to the divergence of the depth-averaged velocity:

$\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} + \frac{\partial (hv)}{\partial y} = 0$

• $h$ is the water depth
• $u$ and $v$ are the depth-averaged velocity components in the $x$ and $y$ directions

In words: if more water flows into a region than flows out, the depth $h$ at that location increases.

Depth-averaged momentum equations

Similarly, integrating the 3D momentum equations over the depth gives two equations (one per horizontal direction) that balance pressure gradients, bottom slope effects, and friction:

$\frac{\partial (hu)}{\partial t} + \frac{\partial (hu^2 + \frac{1}{2}gh^2)}{\partial x} + \frac{\partial (huv)}{\partial y} = -gh\frac{\partial z_b}{\partial x} - \frac{\tau_{bx}}{\rho}$

$\frac{\partial (hv)}{\partial t} + \frac{\partial (huv)}{\partial x} + \frac{\partial (hv^2 + \frac{1}{2}gh^2)}{\partial y} = -gh\frac{\partial z_b}{\partial y} - \frac{\tau_{by}}{\rho}$

• $g$ is gravitational acceleration
• $z_b$ is the bottom elevation (bathymetry)
• $\tau_{bx}$, $\tau_{by}$ are bottom friction stresses
• $\rho$ is fluid density

The $\frac{1}{2}gh^2$ terms come from integrating the hydrostatic pressure over the depth. The right-hand side captures forcing from bottom slopes and frictional drag.

Properties of shallow water equations

Hyperbolic system of equations

The shallow water equations form a hyperbolic system of PDEs. This has direct physical consequences:

• Information propagates at finite speeds (the wave speeds).
• The system admits discontinuous solutions like hydraulic jumps and bores.
• The mathematical structure closely parallels the compressible Euler equations in gas dynamics, with water depth $h$ playing a role analogous to density.

Eigenvalues and eigenvectors

The eigenvalues of the system give the characteristic speeds at which information travels. For the 1D shallow water equations, these are:

$\lambda_1 = u - \sqrt{gh}, \quad \lambda_2 = u + \sqrt{gh}$

The quantity $\sqrt{gh}$ is the shallow water wave speed (analogous to the speed of sound in gas dynamics). So $\lambda_1$ and $\lambda_2$ represent waves traveling upstream and downstream relative to the flow.

The corresponding eigenvectors define the directions in state space along which these waves carry information. This eigenstructure is the foundation for building numerical Riemann solvers and classifying flow regimes:

• Subcritical flow ($|u| < \sqrt{gh}$): Both characteristics point in different directions, so information can travel both upstream and downstream.
• Supercritical flow ($|u| > \sqrt{gh}$): Both characteristics point downstream, and upstream information cannot propagate.

Riemann invariants

Riemann invariants are quantities that stay constant along characteristic curves. For the 1D shallow water equations:

• $R_1 = u - 2\sqrt{gh}$ is constant along $\frac{dx}{dt} = u - \sqrt{gh}$
• $R_2 = u + 2\sqrt{gh}$ is constant along $\frac{dx}{dt} = u + \sqrt{gh}$

These are powerful because they convert the PDEs into ordinary differential equations along characteristics. You can use them to construct exact solutions (like the dam-break problem) and to set boundary conditions correctly based on which characteristics enter or leave the domain.

Characteristic curves

Characteristic curves are the paths in the $x$-$t$ plane along which information propagates. In 1D:

$\frac{dx}{dt} = u - \sqrt{gh} \quad \text{and} \quad \frac{dx}{dt} = u + \sqrt{gh}$

These curves carry the Riemann invariants described above. Where characteristics converge, shocks form. Where they diverge, rarefaction waves develop. Tracing characteristics is one of the most direct ways to understand how a shallow water flow evolves.

Numerical methods for shallow water equations

Finite difference schemes

Finite difference methods approximate derivatives on a structured grid using Taylor expansions. Common time-stepping approaches include:

• Explicit (e.g., Forward Euler): Simple to implement, but restricted by the CFL stability condition $\Delta t \leq \frac{\Delta x}{\max|\lambda|}$.
• Implicit (e.g., Backward Euler): Unconditionally stable for linear problems, but requires solving a system of equations at each time step.
• Semi-implicit (e.g., Crank-Nicolson): Balances accuracy and stability by averaging explicit and implicit contributions.

Finite difference schemes are straightforward to code on regular grids but struggle with complex geometries and can smear discontinuities.

Finite volume methods

Finite volume methods work with the integral (conservative) form of the equations. They divide the domain into control volumes and track fluxes across cell interfaces. This approach:

• Guarantees conservation of mass and momentum at the discrete level
• Handles complex geometries through unstructured meshes
• Naturally captures shocks and discontinuities

The key step is computing the numerical flux at each cell interface, which is where Riemann solvers come in.

Godunov's scheme

Godunov's scheme is the prototypical finite volume method for hyperbolic systems:

1. Assume the solution is piecewise constant within each cell.
2. At each cell interface, solve the Riemann problem between the two adjacent states.
3. Use the Riemann solution to compute the flux across the interface.
4. Update cell averages using these fluxes.

The scheme is conservative and robust at capturing shocks, but it's only first-order accurate. This means it introduces numerical diffusion that smears out sharp gradients and smooth features alike.

High-resolution schemes

To get better than first-order accuracy without introducing spurious oscillations near discontinuities, high-resolution schemes use limiters to control the reconstruction:

• TVD (Total Variation Diminishing) schemes prevent the total variation from increasing, suppressing oscillations.
• ENO (Essentially Non-Oscillatory) schemes adaptively choose the smoothest stencil for reconstruction.
• WENO (Weighted ENO) schemes use a weighted combination of stencils for even better accuracy in smooth regions.

These methods achieve second-order (or higher) accuracy in smooth regions while remaining non-oscillatory near shocks. The trade-off is increased computational cost and implementation complexity.

Applications of shallow water equations

River and channel flows

Shallow water equations are the standard tool for modeling open-channel hydraulics. Typical applications include:

• Flood inundation mapping: Predicting which areas get flooded and to what depth during extreme rainfall events.
• Flow-structure interaction: Simulating how bridges, weirs, and levees affect water levels and velocities.
• Sediment transport: Coupling the flow equations with sediment continuity to study erosion, deposition, and river morphology changes.

The equations capture the essential quantities (depth, velocity, discharge) that engineers need for design and risk assessment.

Tidal flows and storm surges

In coastal regions, the shallow water equations model two important phenomena:

• Tidal flows are driven by gravitational forces from the moon and sun, producing periodic variations in water levels and currents. These are relatively predictable.
• Storm surges are driven by extreme weather (hurricanes, cyclones) that push water onshore through wind stress and low atmospheric pressure. A storm surge can raise water levels by several meters.

Shallow water models predict the extent and timing of coastal inundation, which is critical for evacuation planning and infrastructure design.

Tsunami modeling

Tsunamis are long-wavelength waves triggered by underwater earthquakes, landslides, or volcanic eruptions. Because their wavelength (often hundreds of kilometers) vastly exceeds the ocean depth (a few kilometers), the shallow water approximation applies well during open-ocean propagation.

The equations capture refraction (bending due to depth changes), shoaling (wave height increase as depth decreases), and runup onto land. Operational tsunami warning systems like NOAA's MOST model are built on shallow water solvers.

Atmospheric flows and weather prediction

The shallow water equations serve as a simplified model for large-scale atmospheric dynamics, capturing:

• Rossby waves: Large-scale planetary waves driven by the variation of the Coriolis parameter with latitude.
• Gravity waves: Waves driven by buoyancy, analogous to surface water waves.

Atmospheric shallow water models include the Coriolis force and treat orography (mountain ranges) as bottom topography. They're used as testbeds for numerical weather prediction algorithms and for understanding fundamental atmospheric dynamics before moving to full 3D models.

Limitations and extensions

Validity of shallow water assumptions

Each assumption behind the equations defines a regime where they can fail:

• Hydrostatic assumption breaks down when vertical accelerations become significant (steep waves, flow over abrupt steps).
• Depth-averaging loses information about vertical structure, which matters in stratified flows (e.g., salt wedges in estuaries).
• Slowly varying bathymetry assumption fails near steep underwater cliffs, sharp channel contractions, or abrupt depth changes.

Always check whether your problem satisfies $H/L \ll 1$ before trusting a shallow water model.

Inclusion of bottom topography

Bottom topography enters through the source terms $-gh \frac{\partial z_b}{\partial x}$ and $-gh \frac{\partial z_b}{\partial y}$ on the right-hand side of the momentum equations. Handling these terms correctly in numerical schemes is surprisingly tricky:

• The scheme must preserve the lake-at-rest steady state ($u = v = 0$, $h + z_b = \text{const}$) exactly. Schemes that do this are called well-balanced.
• Bottom steps and obstacles require special treatment to avoid spurious waves.

Coriolis force and rotating frames

For large-scale flows (oceans, atmosphere), Earth's rotation matters. Adding the Coriolis force introduces terms $fhv$ and $-fhu$ to the $x$ and $y$ momentum equations, where $f = 2\Omega \sin\phi$ is the Coriolis parameter ($\Omega$ is Earth's angular velocity, $\phi$ is latitude).

The Coriolis force produces:

• Geostrophic balance: A steady state where pressure gradients balance the Coriolis force, driving flow along (not across) pressure contours.
• Rossby waves: Slow, large-scale waves that arise from the variation of $f$ with latitude.

Non-hydrostatic pressure effects

When vertical accelerations are no longer negligible, you need non-hydrostatic corrections. This happens with:

• Steep or short-wavelength waves
• Rapid changes in bottom topography
• Undular bores (oscillatory wave trains behind a bore front)

Non-hydrostatic shallow water models add a correction pressure term that accounts for vertical acceleration. They're more expensive to solve but capture dispersive wave behavior that the standard equations miss entirely.

Boundary conditions and initial conditions

Inflow and outflow boundaries

Setting boundary conditions correctly depends on the flow regime, which you determine from the characteristic speeds:

• Subcritical inflow (one characteristic entering): Specify one quantity (e.g., discharge) and extrapolate the other from the interior.
• Supercritical inflow (both characteristics entering): Specify both depth and velocity.
• Subcritical outflow (one characteristic leaving): Specify one quantity (e.g., water level) and extrapolate the other.
• Supercritical outflow (both characteristics leaving): Extrapolate everything from the interior.

Mismatching the number of specified conditions with the number of incoming characteristics is a common source of instability.

Solid wall boundaries

Solid walls (coastlines, river banks, structures) enforce the no-penetration condition: the velocity component normal to the wall is zero.

For the tangential component, you choose between:

• Free-slip: Zero shear stress at the wall (tangential velocity is unconstrained). Appropriate when wall friction is negligible.
• No-slip: Tangential velocity is also zero. More realistic but requires resolving the boundary layer.

Corners and sharp edges need careful treatment to avoid numerical artifacts.

Free surface and bottom boundaries

• Free surface: Pressure equals atmospheric pressure, and the kinematic condition requires fluid particles on the surface to stay on the surface. In the depth-averaged framework, this condition is already built into the continuity equation.
• Bottom boundary: The no-penetration condition applies. Bottom friction is typically modeled using empirical formulas:
  • Manning's equation: $\tau_b = \rho g n^2 \frac{|\mathbf{u}|\mathbf{u}}{h^{1/3}}$, where $n$ is Manning's roughness coefficient.
  • Chézy's equation: $\tau_b = \rho g \frac{|\mathbf{u}|\mathbf{u}}{C^2}$, where $C$ is the Chézy coefficient.

Initial conditions for specific problems

Initial conditions specify $h$, $u$, and $v$ everywhere at $t = 0$. Common choices include:

• Quiescent state: Constant depth $h_0$, zero velocity. Used as a baseline or for perturbation studies.
• Steady-state solution: Obtained analytically or from a prior simulation. Useful for studying the response to perturbations.
• Dam-break configuration: A discontinuity in $h$ at some location, with zero velocity everywhere. The classic test problem.

Poor initialization can trigger transient oscillations or instabilities, so it's worth ensuring your initial state is physically consistent (e.g., satisfying the well-balanced condition over variable bathymetry).

Conservation properties and energy

Mass and momentum conservation

The shallow water equations are written in conservation form, meaning each equation has the structure:

$\frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{f}}{\partial x} + \frac{\partial \mathbf{g}}{\partial y} = \mathbf{s}$

where $\mathbf{q}$ is the vector of conserved quantities, $\mathbf{f}$ and $\mathbf{g}$ are flux vectors, and $\mathbf{s}$ is the source term. This structure guarantees that mass and momentum are conserved globally (what enters a region through fluxes must be accounted for).

Numerical schemes should preserve this conservation at the discrete level. Finite volume methods do this by construction; finite difference methods on non-conservative forms can violate conservation near shocks, producing incorrect jump conditions.

Energy conservation and dissipation

The total energy of the shallow water system is the sum of kinetic and potential energy:

$E = \frac{1}{2}\rho h(u^2 + v^2) + \frac{1}{2}\rho g h^2$

Without friction or other dissipation, total energy is conserved. Bottom friction acts as an energy sink, converting kinetic energy into heat. Numerical schemes that artificially create or destroy energy can produce unphysical results, especially in long-duration simulations.

Entropy and entropy production

In the context of hyperbolic conservation laws, entropy serves as a selection criterion for physically correct solutions. Multiple weak solutions can satisfy the conservation equations across a discontinuity, but only one satisfies the entropy condition (the second law of thermodynamics, loosely speaking).

For shallow water equations, entropy production occurs at hydraulic jumps and bores: mechanical energy is irreversibly converted to turbulent dissipation. A numerical scheme that respects the entropy condition will:

• Select the correct (entropy-satisfying) weak solution at shocks
• Prevent non-physical expansion shocks

Symmetries and conservation laws

The shallow water equations exhibit several symmetries:

• Translation invariance in space and time (leading to conservation of momentum and energy)
• Rotation invariance in the absence of Coriolis force (leading to conservation of angular momentum)

By Noether's theorem, each continuous symmetry corresponds to a conservation law. Numerical schemes that preserve these symmetries (sometimes called structure-preserving or geometric integrators) tend to have better long-term stability and accuracy, particularly for problems involving wave propagation over many periods.

Analytical solutions and test cases

Dam-break problem

The dam-break problem is the most widely used test case for shallow water solvers. The setup is simple: a wall separates water at depth $h_L$ (left) from depth $h_R$ (right), and the wall is instantaneously removed at $t = 0$.

The exact solution (on a flat, frictionless bed) consists of:

1. A rarefaction wave propagating upstream into the deeper water
2. A shock wave (or bore) propagating downstream into the shallower water
3. A contact region of constant depth between them

This solution tests a scheme's ability to handle both smooth (rarefaction) and discontinuous (shock) features simultaneously. The dry-bed case ($h_R = 0$) is an especially challenging variant.

Soliton solutions

Solitons are waves that maintain their shape and speed due to a balance between nonlinearity (which steepens the wave) and dispersion (which spreads it out). Strictly speaking, the standard shallow water equations are non-dispersive, so soliton solutions belong to dispersive extensions like the Korteweg-de Vries (KdV) equation or the Kadomtsev-Petviashvili (KP) equation.

These solutions are valuable benchmarks for testing numerical schemes that include dispersive corrections, verifying that the scheme preserves the soliton's shape and speed over long propagation distances.

Riemann problem and Riemann solvers

The Riemann problem generalizes the dam-break problem: given two constant states separated by a discontinuity, find the resulting wave pattern. The solution consists of combinations of shocks, rarefactions, and contact discontinuities.

Common Riemann solvers used in shallow water codes:

• Exact solver: Iteratively solves for the exact intermediate state. Accurate but computationally expensive.
• Roe solver: Linearizes the problem around an averaged state. Fast and accurate for most flows, but can fail at sonic points or near dry beds.
• HLL solver: Estimates the wave speeds bounding the Riemann fan and assumes a single intermediate state. Robust but more diffusive than Roe.

These solvers provide the numerical fluxes in Godunov-type finite volume schemes.

Benchmark test cases for validation

Standard benchmarks for validating shallow water codes include:

• Steady-state problems: Flow over a bump, hydraulic jump in a channel. Tests whether the scheme maintains steady states accurately.
• Transient problems: Dam-break (with and without dry bed), solitary wave propagation. Tests shock-capturing and wave-propagation accuracy.
• Experimental comparisons: Laboratory flume data and field measurements provide ground truth that accounts for effects (turbulence, 3D flow) not captured by the equations themselves.

Comparing your numerical results against these benchmarks is the standard way to verify that a code is implemented correctly and to quantify its accuracy and robustness.