Hydraulic-radius correction to the equilibrium-width transport law (+ competence-limited self-arrest)

[awickert]

Summary

Make the hydraulic radius R_h the explicit variable of the equilibrium-width
theory, and add a linear cofactor that repairs the wide-channel
approximation in the one place it does harm — the sediment-discharge relation.
This costs nothing to the closed-form mathematics and gives the solver a
physically correct competence-limited self-arrest at the narrow-channel
limit, instead of an unphysical extrapolation or a crash.

The insight: R_h is the variable of interest

In the Wickert & Schildgen (2019) equilibrium-width derivation the symbol h
is really the hydraulic radius R_h everywhere except one relation. The
wide-channel approximation R_h ≈ h is invoked notationally throughout but is
only consequential in one place:

Relation	quantity that belongs	paper writes	consequence
Bed shear → threshold closure (Eq. 9)	R_h	h	none — the closure defines this length as R_h
Manning–Strickler u ∝ R_h^(2/3)	R_h	h	none — velocity genuinely wants R_h, and the closure supplies it
Continuity Q = u·b·h (area = b·h)	true flow depth	R_h	propagates into b (see below)

The closure pins the hydraulic radius,

$$R_h = \frac{C_h\,D}{S_c}, \qquad C_h = R\,(1+\varepsilon)\,\tau^*_c, \qquad R=\frac{\rho_s-\rho}{\rho},$$

and shear and velocity use it correctly.

Governing equation (corrected Eq. 20, W&S 2019 corrigendum)

$$\frac{\partial z}{\partial t} = \frac{k_{Qs}\,I}{\mathbb{S}^{7/6}(1-\lambda_p)}\left|\frac{\partial z}{\partial x}\right|^{1/6}\left[\frac{7}{6}\frac{Q}{B}\frac{\partial^2 z}{\partial x^2} + \frac{1}{B}\frac{\partial Q}{\partial x}\frac{\partial z}{\partial x}\right] + U$$

The cofactor: Eq. 16's b is the wetted perimeter

Width-side telling (the paper's own equations). Follow the chain:

• Eq. 2: Q_s = q_s·b — sediment discharge = (per-unit-width bed-load flux) × (the width it crosses, the bed).
• Eq. 8: q_s = k_qs·D^(3/2) — that flux, from MPM at threshold. It is per unit bed width.
• Eqs. 9, 12: h = C_h·D/S and u = 5.9·g^½·h^(2/3)·S^½/D^(1/6). The length here is R_h; both are correct in R_h.
• Eq. 13: q = u·h. This is discharge per unit width, and via Eq. 15 (Q = q·b) it is q = Q/b = u·h_true — per unit bed width. But the paper plugs in the same h = R_h. This is the single wide-channel substitution.
• Eqs. 15–16: b = Q/q = Q/(u·R_h). Since Q/u = A (flow area) and R_h ≡ A/P,

$$b = \frac{Q}{u\,R_h} = \frac{A}{R_h} = \frac{A}{A/P} = P = b_\text{true} + 2\,h_\text{true}.$$

So Eq. 16, b = k_b·Q·S^(7/6)/D^(3/2), is the wetted perimeter — not the bed width. (Dividing the flow area by the hydraulic radius returns the perimeter, by the definition of R_h. In the wide limit h ≪ b, P ≈ b_true, so the substitution is harmless; at finite aspect ratio it overshoots the bed by 2·h_true.)

• Eq. 17: Q_s = q_s·b = k_Qs·I·Q·S^(7/6). Here q_s is per bed width (Eq. 8) but the b it multiplies is the perimeter (Eq. 16), so Q_s overcounts the transporting width by P/b_true = b_wide/b_true = 1/f.

The fix — use the bed width in Eq. 2/17:

$$Q_s = q_s\,b_\text{true} = f\,k_{Qs}\,I\,Q\,S^{7/6}, \qquad f = \frac{b_\text{true}}{b_\text{wide}} = \frac{R_h}{h_\text{true}},$$

with the self-consistent geometry (wide-shallow branch):

$$h_\text{true} = \tfrac{1}{4}\left(b_\text{wide} - \sqrt{\,b_\text{wide}^2 - 8\,R_h\,b_\text{wide}\,}\right), \qquad b_\text{true} = b_\text{wide} - 2\,h_\text{true}.$$

(f = 1 − 2·R_h/b is exact only with the true bed width b_true; using b_wide there is a wide-limit shortcut that fails near the narrow limit.)

Depth-side telling (equivalent). From the other end: continuity used R_h where the true flow depth belongs; restoring the true depth gives the same f. The two are one correction seen from either end — "b is the perimeter, Q_s needs the bed" and "continuity used R_h for the depth."

Why only Q_s needs correcting — water routing is untouched. In continuity Q = u·b·h the perimeter-b and the R_h-h appear as a product, and

$$b\cdot h = P\cdot R_h = P\cdot\frac{A}{P} = A = b_\text{true}\cdot h_\text{true}.$$

The two errors carry the same product (the flow area A), so they cancel — the paper's Q = u·(perimeter)·(R_h) equals u·A, the correct discharge. The discrepancy surfaces only where b appears un-paired with h, i.e. Q_s = q_s·b. Hence exactly one correction, in the transport equation, and water routing was never wrong.

Linear, not (R_h/h)^(13/6). The h^(-13/6) that appears when the rate equation is rewritten in depth coordinates is a coordinate rewrite, not the site of the approximation — which lives in continuity, where it is linear.

Exact, because it iterates. f depends on the current geometry (S from z, and Q), so it is re-evaluated inside the Picard loop from the current iterate and converges with the profile and the |∂z/∂x|^(1/6) nonlinearity — it is not a one-shot post-hoc multiply (which would be only first-order). At convergence the wide-channel approximation is removed exactly — the exact rectangular threshold Q_s at any aspect ratio down to the b/h = 2 floor. The per-node geometry is closed-form (the quadratic above), so this costs essentially nothing; the iteration only handles the coupling to the evolving profile. Applied at the face flux (exact per face).

"Exact" here means the wide-channel approximation is gone. The model's other closures remain — notably the bed/bank shear partition (total boundary shear ρg·R_h·S is pinned to the bed, though banks carry a share in a narrow section) is a separate same-order O(h/b) term this cofactor does not remove.

Limit behaviour: competence-limited self-arrest

As a reach is driven deep and narrow, f → 1/2 at bed aspect ratio b/h = 2,
reached at a discharge floor

$$Q_\text{min} = 8\,u\,R_h^2 \qquad (\Leftrightarrow\ b_\text{wide} = 8\,R_h).$$

Below Q_min no rectangular threshold channel exists — which is not an error
but the onset of competence limitation: the flow can no longer hold the bed
at (1+ε)·τ*_c, so it drops below threshold and Meyer-Peter–Müller carries
transport smoothly to zero,

$$q_s = \phi\,R^{1/2} g^{1/2}\,(\tau^*_b - \tau^*_c)^{3/2}\,D^{3/2} \;\to\; 0 \quad\text{as}\quad \tau^*_b \to \tau^*_c, \qquad \text{then}\ \frac{\partial z}{\partial t}=0.$$

The two regimes join continuously at Q_min. This is exactly the switch
the fixed-width formulation already carries ("if τ*_b < τ*_c,
∂z/∂t = 0"); the equilibrium-width case never reached it because it assumed
threshold was always maintained. The sub-threshold regime needs an inherited
bed width b(x,t) — the dynamic-B work in #19; the two threads converge there.

Validity fences (set by different variables)

1. Aspect / discharge: f ≥ 1/2, i.e. b/h ≥ 2, i.e. Q ≥ Q_min.
2. Grain / slope: R_h ≳ a few D. Because R_h/D = C_h/S, this is a
   slope limit (S ≲ 0.02–0.03) — flow ceases to submerge the roughness and
   the channel becomes a boulder cascade. This is the paper's existing
   Lamb-slope / process-domain boundary and, on steep reaches, it bites first.

Known limitation: exact geometry, not exact closure

The cofactor removes the wide-channel geometry error but carries the
equilibrium-width shear closure τ_bed = (1+ε)·τ*_c (bed shear = 1.2× the bank
threshold) unchanged. That closure pins the bed shear to the reach-average
ρg·R_h·S, valid only while the bed dominates the perimeter. As the channel
narrows:

• Quantitatively, the banks take a growing share of the boundary shear, so the
  reach-average drops below the bed shear — but this drift is ε-suppressed
  (~8% at b/h=2, ~1.5% at b/h=20; about 1/6 the cofactor), sub-dominant in the
  range the cofactor matters.
• In premise, τ_bed = 1.2·τ_bank is Parker's (1978) self-formed
  mobile-bank near-threshold channel — inherently wide-ish. A b/h ≈ 2 deep slot
  is not a Parker channel; the premise fails.

The cofactor cannot rescue a closure whose premise is gone; it just carries it.
The self-arrest floor (b/h = 2) sits about where the premise gives out, so the
model stops there rather than extrapolating a broken closure (and f ≥ 1/2 keeps
it from pushing far in). So "exact" means exact geometry, not exact closure.
True narrow-channel fidelity would need Parker's lateral shear partition (a
bed-specific R_bed, modifying Eq. 9) — a deeper reformulation.

Magnitude

Negligible for substantial rivers: at the model's own equilibrium widths
b/R_h is in the hundreds and f ≈ 1 to sub-percent. The correction is
material (f ~ 0.5–0.9) only for small / steep / coarse reaches near the
theory's edge. This is a rigor + robustness improvement (correct limit
behaviour, no crashes on gentle / low-Q / waning reaches), not a
behaviour-changer for typical applications.

Naming / code change

compute_flow_depth() currently sets self.h from Eq. 9 — but that quantity is
the hydraulic radius, not the flow depth. Plan:

• Rename self.h → self.R_h and compute_flow_depth() → compute_hydraulic_radius().
• Introduce the true flow depth self.h only where it belongs — the
  continuity / cofactor step.
• lp.h is a public attribute (e.g. consumed as flow_depth in the
  characterization tests), so this is a minor breaking rename to sequence
  carefully (deprecation shim or coordinated bump).

Implementation plan

• Rename h → R_h (hydraulic radius); add true flow depth h.
• Compute the self-consistent geometry (h_true, b_true) and cofactor f
  per node from the current iterate.
• Apply f to the face sediment flux; fold into the Picard loop.
• Add the guards: geometric backstop b_wide ≥ 8·R_h, and the grain/slope
  fence R_h/D (warn near the boulder-cascade limit).
• Guard S ≈ 0: never evaluate C_h·D/S; drive the competence check from the
  actual shear ρg·R_h'·S (S in the numerator → 0) so a zero-slope reach
  freezes (Q_s = 0, ∂z/∂t = 0) instead of dividing by zero — this also
  makes flat / S = 0 initialization safe (the current code nans there).
• Sub-threshold regime below Q_min: MPM → 0, then ∂z/∂t = 0
  (competence-limited self-arrest); coordinate the inherited-width need with Valley realism: transient valley widening/narrowing + deposit (overbank) tracking #19.
• Feature switch (default off) so corrected and original wide-channel
  behaviour are both available; convergence test; check f → 1 as
  b/R_h → ∞.

Related

• Valley realism: transient valley widening/narrowing + deposit (overbank) tracking #19 — dynamic valley width B(x,t); the sub-threshold inherited-width need
  converges with that work.
• Aspect-ratio-dependent excess shear: ε(b/h) tapering the multiplier 1.2 → 1 (bank-stability closure) #21 — aspect-ratio-dependent excess shear ε(b/h): a continuous stress-multiplier
  taper (1.2 → 1) that would replace this issue's hard b/h=2 floor + self-arrest.
• Realistic hydraulic geometry: Parker's tapered stable-bank cross-section (beyond rectangular) #22 — realistic tapered (Parker) cross-section: the consistent geometry that
  dissolves the rectangular b/h=2 artifact and the bed/bank "Known limitation".
• Sand-bed conductance singularity at S→0 — may be capped by the hydraulic-radius correction + competence self-arrest (cf. GRLP #20) FluvTree#1 — sand-bed application.

🤖 Generated with Claude Code