Theory and governing equations

GRLP evolves the long profile of a transport-limited gravel-bed river. This page gives the governing equations as amended by the 2020 corrigendum to Wickert & Schildgen (2019) — the forms the code actually solves. The papers are open access (CC-BY); see References. Fuller, source-linked notes on all of the literature live in docs/literature/README.md in the repository.

Symbols

Symbol	Meaning	Code	Typical
\(z\)	valley-floor (bed) elevation	z	—
\(x\)	down-valley distance	x	—
\(Q_s\)	bed-load sediment discharge	—	—
\(Q\)	channel-forming water discharge	Q	—
\(B\)	valley width	B	—
\(U\)	uplift / subsidence (source–sink)	set_uplift_rate	—
\(\lambda_p\)	sediment porosity	lambda_p	0.35
\(I\)	intermittency	intermittency	0.01–1
\(\mathbb{S}\)	sinuosity	sinuosity	\(\ge 1\)
\(k_{Q_s}\)	lumped transport coefficient	k_Qs	\(\approx 0.041\)

Sediment transport

For a near-threshold, self-formed (Parker, 1978) equilibrium-width gravel channel, the bed-load sediment discharge is

\[ Q_s = -\frac{k_{Q_s}\, I}{\mathbb{S}^{7/6}}\; Q\; \frac{\partial z}{\partial x}\, \left|\frac{\partial z}{\partial x}\right|^{1/6}. \]

A defining feature of the equilibrium-width assumption is that grain size cancels: the specific transport rate scales as \(q_s \propto D^{3/2}\) while the channel width scales as \(b \propto D^{-3/2}\), so no grain size enters \(Q_s\). The lumped coefficient \(k_{Q_s} = k_{q_s}\,k_b \approx 0.041\) is built from first principles (Wong & Parker, 2006; Parker, 1978) in Segment.bedload_lumped_constants().

Mass conservation (Exner)

\[ \frac{\partial z}{\partial t} = -\frac{1}{B(1-\lambda_p)}\frac{\partial Q_s}{\partial x} + U. \]

Valley width \(B\) enters only as a local prefactor \(1/B\), not through its downstream derivative. (Keeping \(B\) inside the derivative was the error fixed by the corrigendum — see below.)

Long-profile evolution equation

Combining the two gives a nonlinear diffusion equation:

\[ \frac{\partial z}{\partial t} = \frac{k_{Q_s}\, 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 first bracketed term is diffusive; the second is an advective term arising from downstream discharge growth. There is no valley-width-derivative term.

Boundary conditions

Upstream (Neumann, sediment flux)

A prescribed input sediment discharge sets the boundary slope, \(S_0 = -\,\mathrm{sgn}(Q_s)\,\mathbb{S}\, \bigl(Q_s /(k_{Q_s} I Q)\bigr)^{6/7}\). A higher sediment-to-water supply ratio steepens the boundary (Lane’s balance). Set with set_Qs_input_upstream(Q_s_0), or slope-first with set_S0(S0).

Downstream (Dirichlet, base level)

Elevation is fixed at the outlet with set_z_bl(z_bl). The outlet may also be moved horizontally with set_x_bl — for example to follow a shoreline migrating along a continental-shelf gradient (Ruby et al., 2026).

Numerical method

The equation is solved by a semi-implicit finite-difference scheme. The diffusive part is a tridiagonal system solved directly; the weak nonlinearity from \(|\partial z/\partial x|^{1/6}\) (and any nonlinear \(Q(x)\), \(B(x)\)) is resolved by Picard iteration (set_niter). On a network, one global sparse system is assembled over all nodes by walking the channel topology (grlp.solver).

Steady-state and slope–area

The analytical steady-state profile between endpoints \((x_0, z_0)\) and \((x_1, z_1)\), useful as a verification benchmark, is

\[ z = (z_1 - z_0)\, \frac{x^{a} - x_0^{a}}{x_1^{a} - x_0^{a}} + z_0, \qquad a = 1 - \tfrac{6}{7}P_{x,Q}, \]

with the exponent depending on the distance→discharge power \(P_{x,Q}\) alone. The corresponding slope–area relation is \(S = S_0 (A_0/A)^{(6/7)P_{A,Q}}\), so channel concavity \(\theta = \tfrac{6}{7}P_{A,Q}\) is set by the drainage-area→discharge exponent.

Linearized response: gain and lag

Under small sinusoidal forcing, the equation linearizes to a constant-coefficient diffusion equation for the elevation perturbation, with a sediment-transport diffusivity

\[ \kappa = \frac{7}{6}\, \frac{k_{Q_s}\, I\, \overline{Q}} {\mathbb{S}^{7/6} B(1-\lambda_p)}\, \left|\frac{\partial\bar z}{\partial x}\right|^{1/6}, \qquad T_{eq} = \frac{L^2}{\kappa}. \]

The equilibration time \(T_{eq}\) controls the response: fast forcing (\(P \ll T_{eq}\)) is strongly damped, slow forcing (\(P \gg T_{eq}\)) is nearly in phase with gain \(\to 6/7\). GRLP computes these directly (compute_diffusivity, compute_equilibration_time, compute_z_gain, compute_z_lag, compute_Qs_gain, compute_Qs_lag); the framework is from McNab et al. (2023).

The 2020 corrigendum

The original Wickert & Schildgen (2019) equations mistakenly kept valley width \(B\) inside the spatial derivative of the Exner equation, which incorrectly implied that downstream valley widening or narrowing could by itself drive aggradation or incision. The corrigendum removes every valley-width-derivative term (so \(B\) is a \(1/B\) prefactor), restores an omitted \(\mathbb{S}^{7/6}\) sinuosity factor in the numerical scheme, corrects the upstream boundary condition’s sign, and redefines concavity as \(\theta = \tfrac{6}{7}P_{A,Q}\). The GRLP release that brought the code into line with the corrected equations was v1.4.1. Use the corrected forms on this page.

Networks

A network is a set of long-profile segments joined at confluences, each obeying the evolution equation above, coupled by three junction rules (McNab et al., 2025): (1) a segment’s upstream sediment supply is the sum of its upstream segments’ sediment discharges; (2) elevation is continuous across every junction; (3) water discharge steps up at each tributary junction. Random network topologies are generated as Shreve (1966, 1974) random binary trees (generate_random_network). The integrated network response resembles a single profile whose length is the network’s mean inlet-to-outlet distance, not its longest path.