title: "Chapter 5 — Diffusive-Wave Routing" chapter: 5

Chapter 5 — Diffusive-Wave Routing

Chapter 4 builds the kinematic wave from the Saint-Venant equations and shows it translating a flood pulse downstream without ever changing its shape. That assumption — friction balances gravity, full stop — is exactly right on a steep mountain reach and exactly wrong on a flat valley floor. This chapter is entirely about the gap between those two cases: where the kinematic wave actually gives the wrong answer, what extra physics fixes it, how that physics is written into OPM's production router, and what it costs to run all of it under one single, global time step.

5.1 — Where Kinematic Wave Fails

Drop a stone in a long, still canal and watch the ripple travel away. It does not stay a sharp spike — it spreads out and gets shorter the farther it goes. A flood wave does the same: the sharp crest that leaves the mountains arrives at the city lower and broader, smeared over a longer time. This flattening is called peak attenuation, and the kinematic wave cannot reproduce it — it translates the crest without ever lowering it, so it tends to over-predict the peak discharge at a gauge, exactly the number engineers care most about.

Now pour water steadily into one end of a bathtub with the drain open. Near the drain the water is shallow and quick; but if the drain is small, water piles up behind it and the level rises upstream. Rivers do this constantly — a pond, a flat reach, a narrow culvert, or a downstream lake raises the water level and slows the flow coming in from upstream. This is the backwater effect, and the kinematic wave is blind to it too: it uses only the bed slope, which always points downhill and knows nothing about how deep the water already is in the next cell. A kinematic cell keeps draining at full Manning's rate even when the cell below it is brim-full and physically cannot accept more.

Both pictures correspond to the kinematic-wave number from Chapter 2, $K = S_0\,g\,L / V_0^2$, falling below the safe threshold of about 20:

1. Flat slopes ($S_0 < 10^{-3}$) — the bed slope no longer dominates; the water's own depth differences matter as much as gravity.
2. Peak attenuation — the kinematic wave keeps the crest sharp; reality flattens it.
3. Backwater and adverse gradients — a downstream pool can be deeper than the upstream cell, opposing the flow; the kinematic wave never sees it.

The widget below runs the same two-cell setup three ways — Normal, Backwater, and Adverse — holding the upstream cell's depth fixed and only changing how full the downstream cell is. Watch $Q_\text{kin}$: it never moves. It cannot, because the kinematic formula never looks at $h_\text{ds}$.

In the Adverse lane, the downstream water surface is higher than the upstream one — flow should stop. The diffusive scheme correctly clamps to $Q_\text{diff} = 0$. The kinematic scheme keeps draining at $0.02024\ \text{m}^3/\text{s}$ regardless, into an already-fuller cell — not a rounding error, a physically wrong answer.

5.2 — Follow the Water Surface, Not the Bed

The fix is one idea: drive flow with the slope of the water surface, not the slope of the bed. Define the water-surface elevation as bed plus depth, $\text{WSE} = z + h$, and replace the bed slope $S_0$ with

Three cases fall out immediately:

• $h_i > h_\text{ds}$ — the upstream cell is fuller than the bed slope alone would suggest. $S_w > S_0$: a small boost.
• $h_i = h_\text{ds}$ — both cells equally full. $S_w = S_0$ exactly: identical to kinematic.
• $h_i < h_\text{ds}$ — the downstream water surface is higher. $S_w < S_0$, and once $h_\text{ds} - h_i$ exceeds $S_0\,\Delta x$, $S_w$ goes negative. OPM clamps it: $S_\text{eff} = \max(S_w, 0)$ — flow halts rather than running backward.

Drag $h_i$ and $h_\text{ds}$ below and watch the water-surface line tilt, flatten, and reverse. With $S_0 = 0.001$ and $\Delta x = 100\,\text{m}$, the exact halt threshold is $h_\text{ds} - h_i = 0.10\,\text{m}$ — verify it on the slider.

5.3 — Deriving the Diffusive-Wave Equation From Scratch

Chapter 4 showed the full Saint-Venant momentum equation has five terms, and that the kinematic wave keeps only friction and gravity:

Kinematic drops ①②③, leaving $S_f = S_0$. The diffusive wave makes a smaller cut: it drops only the inertial terms ①② — local and convective acceleration, the two terms responsible for genuine flood-wave dynamics like surging and reflection — but it keeps the pressure-gradient term ③. What survives is

which is exactly the continuous-space version of the discrete $S_w$ from §5.2. Linearizing around a reference depth $h_0$ and velocity $V_0$ (small perturbation $h' = h - h_0$) turns the diffusive-wave momentum balance into an advection–diffusion equation:

The celerity $c_k$ is identical to Chapter 4's kinematic celerity — the wave still travels at $\tfrac{5}{3}$ the flow velocity. The new term is $D$, a genuine diffusivity: it spreads and lowers the pulse exactly the way heat spreads through a rod. Because $S_0$ sits in the denominator, flat slopes ($S_0\to 0$) blow $D$ up — flat reaches diffuse strongly, steep reaches barely diffuse at all:

The widget below walks through the full derivation step by step, then lets you turn a $D$ slider and watch a pulse spread on a 1-D grid — the same equation, made visible.

5.4 — The Real OPM Code: GSSHA-Style Conveyance Depth

$S_w$ alone is not quite enough. On a normal downhill bed the upstream cell's own depth $h_i$ is the right depth to carry the flow. But at a D8 saddle or a local DEM step, the downstream cell's ground can sit higher than the upstream cell's ($z_\text{ds} > z_i$) even while its water surface is still lower — water is conveyed across that higher lip, and the relevant depth is the depth of water over that higher bed, not over the upstream one.

OPM's diffusive_wave_discharge() in routing_utils.py borrows the answer from the CASC2D family of models (Julien et al., 1995; later GSSHA, and the LISFLOOD-FP diffusion wave): use the depth above whichever of the two beds is higher,

A single blend weight $\theta$ ties the slope and the conveyance depth together as a coherent pair: $\theta=0$ reduces to kinematic exactly (own depth, bed slope); $\theta=1$ — OPM's production default — is the full diffusion wave (water-over-higher-bed, water-surface slope). The widget below renders the actual, unabridged function body from the production code, with six numbered callouts, plus the ds_safe trick that lets every cell — including the outlet, which has no downstream neighbor — run the identical branch-free formula on the GPU.

def diffusive_wave_discharge(depth, dem, dist, slope_bnd, n, ds_safe, valid_ds,
                             theta, cell_size, xp, min_depth):
    depth_ds = depth[ds_safe]
    dem_ds   = dem[ds_safe]

    # ① water-surface slope = bed slope + theta * depth gradient
    S_w = slope_bnd + theta * (depth - depth_ds) / dist
    # ② outlet cells fall back to bed slope; clamp adverse gradient to 0
    S_eff = xp.where(valid_ds, S_w, slope_bnd)
    S_eff = xp.maximum(S_eff, 0.0)

    # ③ water-surface elevations (this cell and its downstream neighbour)
    wse      = dem + depth
    wse_ds   = wse[ds_safe]
    # ④ CASC2D / GSSHA conveyance depth = water over the higher bed
    h_higher = xp.maximum(wse, wse_ds) - xp.maximum(dem, dem_ds)
    # ⑤ blend own depth (kinematic) with h_higher (diffusive) by theta
    h_flow   = (1.0 - theta) * depth + theta * h_higher
    h_flow   = xp.where(valid_ds, h_flow, depth)
    h_flow   = xp.maximum(h_flow, min_depth)

    # ⑥ Manning-style discharge on (h_flow, S_eff)
    return (1.0 / n) * (h_flow ** (5.0 / 3.0)) * (S_eff ** 0.5) * cell_size

valid_ds = ds_idx >= 0                       # True for interior cells
ds_safe  = xp.where(valid_ds, ds_idx, 0)     # outlet's -1 -> 0 (any valid index)

5.5 — One Global Δt, Many Slopes: The Hidden Cost

Every explicit scheme needs a stable $\Delta t$, and OPM picks exactly one for the whole watershed. Look at the real self-check in kinematic_wave_router.py:

max_slope         = float(slope_1d.max().item())
n_min             = float(n.min().item()) if hasattr(n, 'min') else float(n)
V_at_1m_max_slope = (1.0 / n_min) * (1.0 ** (2.0 / 3.0)) * (max_slope ** 0.5)
C_indicator       = V_at_1m_max_slope * dt / dx
safe_dt           = dx / V_at_1m_max_slope

This runs once, on the static bed-slope grid, using the single steepest cell and the smallest roughness anywhere in the basin — and it runs identically whether ROUTING_SCHEME is 'kinematic' or 'diffusive'. Diffusive routing inherits the kinematic CFL's global $\Delta t$; it does not get its own.

That single choice is sized for the steepest cell, which means every gentler cell is running far below its own stability limit — and an under-sized Courant number is exactly what manufactures numerical diffusion ($D_\text{num} \approx \tfrac{c\,\Delta x}{2}(1-C)$, from Chapter 4 §4.4; the global-vs-local trade-off is Chapter 4 §4.5). The widget below puts OPM's real $\Delta x \approx 100\,\text{m}$ and default $\Delta t = 0.7\,\text{s}$ into three illustrative zones spanning a real watershed's slope range:

max_slope         = float(slope_1d.max().item())
n_min             = float(n.min().item()) if hasattr(n, 'min') else float(n)
V_at_1m_max_slope = (1.0 / n_min) * (1.0 ** (2.0 / 3.0)) * (max_slope ** 0.5)
C_indicator       = V_at_1m_max_slope * dt / dx
safe_dt           = dx / V_at_1m_max_slope

At the default settings, the flat valley-floor cell's Courant number is 164× smaller than the steep headwater cell that set $\Delta t$ in the first place — and it is 131× smaller than Chapter 4's already-failing "Too Small Δt" demo lane. The consequence: at that flat cell, the scheme's own numerical smoothing is already ≈28% the size of the physical backwater diffusivity it was built to capture. $D_\text{num}$ and $D_\text{phys}$ enter the exact same equation, so OPM has no way to report which fraction of a flat-valley hydrograph's attenuation is real physics and which is a grid artifact. (At the steep cell the ratio looks huge too, but both numbers are negligible there in absolute terms — consistent with §5.3's table that steep slopes barely diffuse regardless.)

5.6 — Should OPM Have Used Implicit?

An implicit scheme removes the CFL limit entirely — but it is not free:

Chapter 4 §4.6 already showed a real implicit scheme — the Preissmann θ-method solved by the Thomas algorithm — and noted it works because a single 1-D channel is a chain: each interior cell has exactly one upstream and one downstream neighbor, which collapses the linear system to a lower-bidiagonal matrix solvable in one forward sweep.

OPM's domain is not a chain. It is a watershed-wide D8 flow-direction graph, assembled and routed with scatter_add, where a confluence — two or more upstream cells feeding one downstream cell — is the normal case, not the exception.

A confluence breaks the tridiagonal/bidiagonal structure a Thomas algorithm needs: the downstream row in the matrix now has two or more nonzero upstream entries, not one. Solving that requires a general sparse solver, not a cheap single-sweep elimination — architecturally different, and far more expensive, at the $10^4$–$10^6$-cell scale OPM targets.

Verdict: explicit was the right call for OPM's actual architecture and scale — the per-step cost of a sparse implicit solve over a branching DAG would dwarf the savings, and the flux limiter already neutralizes the one thing implicit buys (robustness to large $C$). But that decision leaves the §5.5 cost sitting on the table: a watershed with both steep and flat cells pays needless numerical diffusion on every flat cell, every step. The fix that would address it without a full implicit rewrite is local/block time stepping, a TVD flux, or Muskingum–Cunge matched diffusion (Chapter 4 §4.5 lays out all four routes and the GPU trade-off) — none implemented in OPM today, and left here as clearly scoped future work rather than a flaw to patch immediately.

5.7 — Side-by-Side & Mass Balance

This is why OPM defaults to ROUTING_SCHEME = 'diffusive': a basin with both steep hillslopes and flat valley floors is handled correctly by one scheme, because the diffusive formula automatically collapses to the kinematic answer wherever $h_\text{higher} \approx h_i$ and $S_w \approx S_0$.

Adding backwater and flow-halting does not threaten mass conservation. The volume update from Chapter 2 is completely unchanged —

— only the formula for $Q_\text{out}$ differs between schemes, and both are capped by the identical flux limiter and tracked by the identical mass-balance accumulators. When $S_\text{eff} = 0$, $Q_\text{diff} = 0$: no water leaves the cell, the limiter is trivially satisfied, and the volume is simply unchanged — stopping flow does not lose water.

Drag $\theta_\text{diff}$ below to separate the kinematic, blended, and full-diffusive curves on the same hydrograph:

5.8 — Chapter Summary

Key takeaways:

1. The kinematic wave fails on flat slopes, cannot attenuate peaks, and cannot feel backwater — all because it uses the static bed slope.
2. The diffusive wave's one idea: drive flow with the water-surface slope $S_w$ instead of $S_0$.
3. Three cases — $S_w > S_0$ (normal boost), $S_w = S_0$ (identical to kinematic), $S_w < 0$ (clamped to 0, flow halts).
4. The $\theta$-blend makes $\theta=0$ exactly kinematic and $\theta=1$ fully diffusive; OPM uses $\theta=1$.
5. The CASC2D/GSSHA conveyance depth $h_\text{higher}$ correctly handles adverse slopes and reduces to $h_i$ for normal downhill flow.
6. The ds_safe trick replaces the outlet's $-1$ index with 0 and masks the result — a branch-free, GPU-friendly downstream lookup.
7. OPM's single global $\Delta t$, sized for the steepest cell, leaves flat cells running at a tiny fraction of their own stability limit — manufacturing numerical diffusion that is indistinguishable, inside the model, from the real physical diffusivity it was built to capture.
8. Explicit was the right architectural call for OPM's branching D8 domain — a true implicit solve would need a general sparse solver, not the cheap single-sweep Thomas algorithm a 1-D channel allows — but a per-zone or flow-adaptive $\Delta t$ remains the clean fix for the cost identified in §5.5, and is not yet implemented.