Note

This page was generated from a jupyter notebook.

Unit Tests for the Landlab GravelBedrockEroder Component

G.E. Tucker, CIRES and Department of Geological Sciences, University of Colorado Boulder

This notebook describes a series of basic unit tests for the GravelBedrockEroder component. These tests are implemented in the file test_gravel_bedrock_eroder.py, and are in addition to the unit tests implemented as doctests in the source code.

Types of test

The theory starts with a representation of a layer of sediment overlying rock. Each grid cell is assumed to contain a primary channel that drains to an adjacent cell. The cell may also receive inflow of water and sediment from neighboring cells.

Processes and effects include:

1. Dynamic adjustment of channel width, based on near-threshold theory (implicit in derivation; not calculated explicitly unless requested via function call)

2. Transport of coarse (assumed gravel-size) sediment as bed load

3. Abrasion of sediment, which turns coarse sediment into wash load (not tracked)

4. Abrasion of underlying bedrock

5. Plucking erosion of underlying bedrock, which supplies coarse sediment

Ideally, each of these elements should be tested, both separately and in combination. Two types of test are used: instantaneous tests, which are “single iteration” comparisons between predicted and computed values, and equilibrium tests, in which a small terrain is run with baselevel forcing (“uplift”) until it reaches an equilibrium that is then compared with an independently calculated solution.

Parameters, their default values, and mathematical symbols are:

• intermittency_factor=0.01 (\(I\)) dimensionless

• transport_coefficient=0.041 (\(k_Q\)) dimensionless

• abrasion_coefficient=0.0 (\(\beta\)) [1/L]

• sediment_porosity=0.35 (\(\phi\)) dimensionless

• depth_decay_scale=1.0 (\(H_*\)) [L]

• plucking_coefficient=1.0e-6 (\(k_p\)) [1/L]

• coarse_fraction_from_plucking=1.0 (\(\gamma\)) dimensionless

Mathematical symbols for the GravelBedrockEroder state variables and their corresponding at-node fields are:

• \(\eta\) = channel elevation = topographic__elevation

• \(H\) = sediment thickness = soil__depth

• \(\eta_b\) = bedrock surface elevation = bedrock__elevation

• \(S\) = flow-wise slope gradient = topographic__steepest_slope

• \(Q\) = bankfull discharge [L\(^3/T\)] = surface_water__discharge (externally calculated)

• \(Q_{in}\) = sediment flux entering the cell at a node [L\(^3\)/T] = bedload_sediment__volume_influx

• \(Q_{out} = Q_s\) = sediment flux leaving a cell along its flow link [L\(^3\)/T] = bedload_sediment__volume_outflux

• \(\alpha\) = bedrock exposure fraction [-] = bedrock__exposure_fraction

Mathematical symbols for other variables used in the calculations below:

• \(U\) = uplift rate relative to baselevel [L/T]

• \(\Lambda\) = grid cell (projected) surface area [L\(^2\)]

• \(\lambda\) = grid link length [L]

• \(t\) = time (in the underlying governing equations)

In addition, runoff rate \(r\) is provided indirectly, through the field surface_water__discharge. Its default is considered to be unity (1 m/y). Bankfull discharge \(Q\) in these tests is equal to \(r\) times drainage area \(A\), which is equal to cell area \(\Lambda\) times the number of cells drained. The local flow-directed slope gradient, \(S\), is calculated at each grid node. Rock uplift (or sediment uplift, in the unlimited sediment case) is applied as a boundary condition.

Channel width is adequately tested by the doctests, so additional tests of width are not included here.

Test setup

The doctests in the code use a RasterModelGrid. For these external tests, we will use a HexModelGrid with three core nodes and a single open boundary node. This configuration is small enough for tests to be quick and efficient, but large enough to include flow convergence (two cells feel flow into a third, which then drains to the open boundary).

Imports

import numpy as np
from numpy.testing import assert_almost_equal

from landlab import HexModelGrid
from landlab.components import FlowAccumulator, GravelBedrockEroder

Instantaneous tests

Transport rate

Test condition: 3-core-node hex grid with gradient of 0.01 along flow paths. Sediment cover: ample (100 m), limited (0.5 m, with depth decay scale also set to 0.5 m), and none.

Predicted sediment transport rate under ample cover (not limited by bedrock exposure):

\[Q_s = k_Q I Q S^{7/6}\]

Here the default value \(k_Q=0.041\) is used, but the intermittency factor is set to 0.02. The discharge by default will be one meter per year times drainage area. The drainage area of one cell is the cell’s area, here about 866,025 m\(^2\), and the drainage area of the cell that receives flow is three times this. Therefore the predicted sediment transport rate for the two “upstream” cells and for the single “downstream” cell is:

Q = (3.0**0.5 / 2.0) * 1e6 * np.array([3, 1, 1])
Qs = 0.041 * 0.02 * Q * (0.01) ** (7.0 / 6.0)
print("Predicted transport rate:", Qs)

For the case with limiting sediment cover, when the cover thickness is equal to the depth decay scale (set to 0.5 m), the transport rate should be reduced by a factor of \(1 - 1/e\). This works out to:

print("Predicted transport rate:", Qs * (1.0 - np.exp(-1.0)))

Finally, with no sediment at all, the transport rate should be zero.

Note that in order to give the grid nodes a gradient of 0.01, the elevations need to rise in the y-direction at a rate equal to 0.01 / cos 30\(^\circ\).

def test_transport_rate():
    grid = HexModelGrid((4, 2), spacing=1000.0)
    grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
    grid.status_at_node[0] = grid.BC_NODE_IS_FIXED_VALUE

    elev = grid.add_zeros("topographic__elevation", at="node")
    elev[:] = (0.01 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed = grid.add_zeros("soil__depth", at="node")
    sed[:] = 100.0

    fa = FlowAccumulator(grid)
    fa.run_one_step()
    gbe = GravelBedrockEroder(grid, intermittency_factor=0.02, depth_decay_scale=0.5)
    rock = grid.at_node["bedrock__elevation"]
    qs_out = grid.at_node["bedload_sediment__volume_outflux"]

    gbe.run_one_step(1.0e-6)  # using dt=0 prevents change to sed, rock, or elev
    assert_almost_equal(qs_out[grid.core_nodes], [9.88854526, 3.29618175, 3.29618175])

    elev[:] = (0.01 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed[:] = 0.5
    rock[:] = elev - sed

    gbe.run_one_step(1.0e-6)
    assert_almost_equal(qs_out[grid.core_nodes], [6.25075275, 2.08358425, 2.08358425])

    elev[:] = (0.01 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed[:] = 0.0
    rock[:] = elev

    gbe.run_one_step(1.0e-6)
    assert_almost_equal(qs_out[grid.core_nodes], [0.0, 0.0, 0.0])

test_transport_rate()

Sediment abrasion rate

Consider the first of the cases above, in which the transport rate is about 3.3 m\(^3\)/y for the upstream cells and 9.9 m\(^3\)/y for the downstream ones. If the abrasion coefficient is 10\(^{-4}\) m\(^{-1}\), then we can calculate the resulting lowering rate of the thickness of sediment as:

\[\frac{dH}{dt} = -\beta (Q_{in} + Q_{out}) \lambda / 2 \Lambda\]

where \(\beta\) is the abrasion coefficient, \(Q_{in}\) is incoming sediment flux (m\(^3\)/y), \(Q_{out}\) is the outgoing sediment flux, \(\lambda\) is the length of the flow path (distance from the node to its downstream neighbor), and \(\Lambda\) is the surface area of the cell. The factor of two is there to average between \(Q_{in}\) and \(Q_{out}\). The flow length \(\lambda\) appears because the abrasion rate \(\beta Q_s\) is the rate per length, so we need to multiply by length to get the total volume rate. Finally, cell area \(\Lambda\) appears in the denominator in order to convert a volume rate to a lowering rate.

In this case, the numbers are as follows (here sediment flux is half the above case because we are using the default intermittency factor):

beta = 1.0e-4  # abrasion coefficient, 1/m
Qout = 0.5 * 3.29618175  # transport rate, m3/y
path_length = 1000.0  # node spacing, m
cell_area = 1000.0 * 1000.0 * 0.5 * 3.0**0.5
print(
    "Rate of thickness loss from sediment abrasion (upstream):",
    beta * 0.5 * (0.0 + Qout) * path_length / cell_area,
)

Qin = 2 * Qout
Qout = 0.5 * 9.88854526
print(
    "Rate of thickness loss from sediment abrasion (downstream):",
    beta * 0.5 * (Qin + Qout) * path_length / cell_area,
)

def test_sediment_abrasion_rate():
    grid = HexModelGrid((4, 2), spacing=1000.0)
    grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
    grid.status_at_node[0] = grid.BC_NODE_IS_FIXED_VALUE

    elev = grid.add_zeros("topographic__elevation", at="node")
    elev[:] = (0.01 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed = grid.add_zeros("soil__depth", at="node")
    sed[:] = 100.0

    fa = FlowAccumulator(grid)
    fa.run_one_step()
    gbe = GravelBedrockEroder(grid, abrasion_coefficient=1.0e-4)
    gbe.run_one_step(1.0)

    assert_almost_equal(
        gbe._abrasion[grid.core_nodes],
        [4.7576285545378313e-07, 9.515257103302159e-08, 9.515257103302159e-08],
    )

test_sediment_abrasion_rate()

Bedrock abrasion

Here we test the calculation of bedrock abrasion rate. We need a test that has some sediment, but not so much that the bed is totally shielded. We’ll use 1 m thick sediment. That reduces the transport capacity, which should be equal to the above transport rates times the fractional alluvial cover, which is 1 minus the bedrock exposure fraction:

beta = 1.0e-4
path_length = 1000.0
frac_bed_exposed = np.exp(-1.0)
cell_area = 1.0e6 * 0.5 * 3.0**0.5
Qs_out = (
    0.041
    * 0.01
    * 0.01 ** (7.0 / 6.0)
    * cell_area
    * np.array([3, 1, 1])
    * (1.0 - frac_bed_exposed)
)
Qs_in = np.array([Qs_out[1] + Qs_out[2], 0.0, 0.0])

print("Sed outflux:", Qs_out)
print("Sed influx:", Qs_in)

sed_abr_rate = beta * 0.5 * (Qs_in + Qs_out) * path_length / cell_area
print("Sediment abrasion rate:", sed_abr_rate)

rock_abr_rate = sed_abr_rate * frac_bed_exposed
print("Bed abrasion rate:", rock_abr_rate)

def test_rock_abrasion_rate():
    grid = HexModelGrid((4, 2), spacing=1000.0)
    grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
    grid.status_at_node[0] = grid.BC_NODE_IS_FIXED_VALUE

    elev = grid.add_zeros("topographic__elevation", at="node")
    elev[:] = (0.01 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed = grid.add_zeros("soil__depth", at="node")
    sed[:] = 1.0

    fa = FlowAccumulator(grid)
    fa.run_one_step()
    gbe = GravelBedrockEroder(grid, abrasion_coefficient=1.0e-4)
    gbe.run_one_step(1.0)

    assert_almost_equal(
        gbe._sediment_outflux[grid.core_nodes], [3.12537638, 1.04179213, 1.04179213]
    )
    assert_almost_equal(
        gbe._rock_abrasion_rate[grid.core_nodes],
        [1.10635873e-07, 2.21271745e-08, 2.21271745e-08],
    )

test_rock_abrasion_rate()

Plucking erosion

Here we test the calculation of bedrock plucking rate.

plucking_coef = 1.0e-6
intermittency_factor = 0.01
frac_bed_exposed = np.exp(-1.0)
flow_link_length = 1000
cell_area = (3.0**0.5 / 2.0) * 1e6
Q = cell_area * np.array([3, 1, 1])
slope = 0.01

pluck_rate = (
    plucking_coef
    * intermittency_factor
    * Q
    * slope ** (7.0 / 6.0)
    * frac_bed_exposed
    * (flow_link_length / cell_area)
)
print("Plucking rate:", pluck_rate)

grid = HexModelGrid((4, 2), spacing=1000.0)
grid.length_of_link[0]
grid.area_of_cell[0]

def test_rock_plucking_rate():
    grid = HexModelGrid((4, 2), spacing=1000.0)
    grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
    grid.status_at_node[0] = grid.BC_NODE_IS_FIXED_VALUE

    elev = grid.add_zeros("topographic__elevation", at="node")
    elev[:] = (0.01 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed = grid.add_zeros("soil__depth", at="node")
    sed[:] = 1.0

    fa = FlowAccumulator(grid)
    fa.run_one_step()
    gbe = GravelBedrockEroder(grid, plucking_coefficient=1.0e-4)
    gbe.run_one_step(1.0)

    assert_almost_equal(
        gbe._pluck_rate[grid.core_nodes],
        [5.12263532e-06, 1.70754511e-06, 1.70754511e-06],
    )

test_rock_plucking_rate()

Equilibrium tests

Case of unlimited sediment

We start with an all-sediment case, using the same 3-cell configuration as above. We will impose a rate of uplift relative to the fixed baselevel node of \(U=10^{-4}\) m/y. Once the system reaches steady state, each of the upper cells should show a balance between the rate of sediment input via uplift \(U\Lambda\), the rate of loss to abrasion \(\beta Q_{out}^0 / 2\), and the rate of output downstream, \(Q_{out}^0\) (the superscript zero notation means we’re talking about the upper two cells):

\[Q_{out}^0 + \beta \lambda Q_{out}^0/2 = (1-\phi)U\Lambda\]

Solving for \(Q_{out}^0\),

\[Q_{out}^0 = \frac{(1-\phi)U\Lambda}{1+\beta\lambda /2}\]

The balance for the downstream cell is similar, except that there are two inputs as well. Denoting outflux from the downstream cell as \(Q_{out}^1\),

\[Q_{out}^1 + \beta\lambda (Q_{out}^1 + 2 Q_{out}^0)/2 = (1-\phi)U\Lambda + 2Q_{out}^0\]

or

\[Q_{out}^1 = \frac{(1-\phi)U\Lambda + (2 - \beta\lambda) Q_{out}^0}{1+\beta\lambda /2}\]

Using \(U=10^{-4}\) m/y, \(\Lambda = 866,025\) m\(^2\), \(\lambda = 1000\) m, \(\phi = 0.35\), and \(\beta = 0.0005\), here are the predicted values:

beta = 0.0005
U = 0.0001
Lambda = 866025.0
length = 1000.0  # "little lambda"
phi = 0.35  # porosity

Qout0 = (1 - phi) * U * Lambda / (1 + beta * length / 2)
print("Qout0:", Qout0, "m3/y")

Qout1 = ((1 - phi) * U * Lambda + (2 - beta * length) * Qout0) / (1 + beta * length / 2)
print("Qout1:", Qout1, "m3/y")

We can find the corresponding slope from the flux laws. For the upper nodes,

\[Q_{out}^0 = k_Q I r \Lambda S^{7/6}\]

\[S_0 = \left(\frac{Q_{out}^0}{k_Q I r \Lambda}\right)^{6/7}\]

For the lower node,

\[S_1 = \left( \frac{Q_{out}^1}{3 k_Q I r \Lambda} \right)^{6/7}\]

kQ = 0.041
intermittency = 0.01
r = 1.0

S_pred0 = (Qout0 / (kQ * intermittency * r * Lambda)) ** (6.0 / 7.0)
print("Predicted slope at upper nodes:", S_pred0)

S_pred1 = (Qout1 / (3 * kQ * intermittency * r * Lambda)) ** (6.0 / 7.0)
print("Predicted slope at lower node:", S_pred1)

The test code follows.

def test_steady_unlimited_sediment():
    grid = HexModelGrid((4, 2), spacing=1000.0)
    grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
    grid.status_at_node[0] = grid.BC_NODE_IS_FIXED_VALUE

    elev = grid.add_zeros("topographic__elevation", at="node")
    elev[:] = (0.13 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed = grid.add_zeros("soil__depth", at="node")
    sed[:] = 10000.0
    rock = grid.add_zeros("bedrock__elevation", at="node")
    rock[:] = elev - sed

    fa = FlowAccumulator(grid)
    fa.run_one_step()
    gbe = GravelBedrockEroder(grid, abrasion_coefficient=0.0005)

    dt = 4.0e4
    uplift_rate = 0.0001
    nsteps = 500
    for i in range(nsteps):
        elev[grid.core_nodes] += uplift_rate * dt
        sed[grid.core_nodes] += uplift_rate * dt
        gbe.run_one_step(dt)

    assert_almost_equal(
        grid.at_node["bedload_sediment__volume_outflux"][grid.core_nodes],
        [99.073, 45.033, 45.033],
        decimal=2,
    )
    assert_almost_equal(
        grid.at_node["topographic__steepest_slope"][grid.core_nodes],
        [0.130579, 0.170346, 0.170346],
        decimal=5,
    )

test_steady_unlimited_sediment()

Case of general equilibrium

For this test, we will use a grid with a single core node.

Volume rate of gravel sediment produced by plucking is the areal average rate times cell area,

\[Q_p = k_p I r S^{7/6} \gamma \Lambda \lambda \alpha\]

Equilibrium \(Q_s\) equals sediment produced from plucking minus sediment lost to abrasion:

\[Q_s = k_p I r S^{7/6} \gamma \Lambda \lambda \alpha - \beta\lambda Q_s / 2\]

Solve for \(Q_s\):

\[\boxed{Q_s = \frac{k_p I r S^{7/6} \gamma \Lambda \lambda \alpha}{1 + \beta\lambda / 2}}\]

Equate this with transport rate:

\[k_Q I r \Lambda S^{7/6} (1 - \alpha) = \frac{k_p I r S^{7/6} \gamma \Lambda \lambda \alpha}{1 + \beta\lambda / 2}\]

Cancel common factors:

\[k_Q (1 - \alpha) = \frac{k_p \gamma \lambda \alpha}{1 + \beta\lambda / 2}\]

Rearrange to isolate \(\alpha\):

\[\frac{1 - \alpha}{\alpha} = \frac{k_p \gamma \lambda}{k_Q(1 + \beta\lambda / 2)}\]

\[\frac{1}{\alpha} - 1 = \frac{k_p \gamma \lambda}{k_Q(1 + \beta\lambda / 2)}\]

\[\alpha = \frac{1}{1 + \frac{k_p \gamma \lambda}{k_Q(1 + \beta\lambda / 2)}}\]

\[\boxed{\alpha = \frac{k_Q(1 + \beta\lambda / 2)}{k_Q(1 + \beta\lambda / 2) + k_p \gamma \lambda}}\]

To understand this equation for \(\alpha\), note that if plucking produces no coarse sediment at all—that is, if \(\gamma=0\), then the bedrock exposure fraction is unity. One the other hand, if the plucking process is highly efficient, such that the second term in the denominator is much bigger than the first, then \(\alpha\rightarrow 0\), meaning that the bed is mostly covered by sediment and only the most miniscule bedrock exposure fraction suffices to allow downcutting.

Equilibrium bedrock lowering rate:

\[U = k_p I r S^{7/6} \lambda \alpha + \beta\lambda Q_s \alpha / 2 \Lambda\]

Rearrange to isolate \(Q_s\):

\[\beta\lambda Q_s \alpha / 2 \Lambda = U - k_p I r S^{7/6} \lambda \alpha\]

\[Q_s = 2\Lambda\frac{U/\alpha - k_p I r S^{7/6} \lambda}{\beta\lambda}\]

Plug in transport rate,

\[k_Q I r \Lambda S^{7/6} (1 - \alpha) = 2\Lambda\frac{U/\alpha - k_p I r S^{7/6} \lambda}{\beta\lambda}\]

Rearrange to solve for slope,

\[k_Q I r \Lambda S^{7/6} (1 - \alpha) = \frac{2\Lambda U/\alpha}{\beta\lambda} - \frac{2\Lambda k_p I r S^{7/6} \lambda}{\beta\lambda}\]

\[k_Q I r \Lambda S^{7/6} (1 - \alpha) + \frac{2\Lambda k_p I r S^{7/6} \lambda}{\beta\lambda}= \frac{2\Lambda U/\alpha}{\beta\lambda}\]

\[\left(k_Q I r \Lambda (1 - \alpha) + \frac{2\Lambda k_p I r \lambda}{\beta\lambda}\right) S^{7/6} = \frac{2\Lambda U/\alpha}{\beta\lambda}\]

\[S^{7/6} = \frac{2\Lambda U/\alpha}{\beta\lambda\left(k_Q I r \Lambda (1 - \alpha) + \frac{2\Lambda k_p I r \lambda}{\beta\lambda}\right)}\]

\[S^{7/6} = \frac{2\Lambda U/\alpha}{\beta\lambda k_Q I r \Lambda (1 - \alpha) + 2\Lambda k_p I r \lambda}\]

\[S^{7/6} = \frac{2 U}{\beta\lambda k_Q I r \alpha (1 - \alpha) + 2k_p I r \lambda\alpha}\]

\[\boxed{S = \left( \frac{2U}{\lambda I r \alpha (\beta k_Q (1 - \alpha ) + 2k_p)} \right)^{6/7}}\]

This has some familiar pieces to it in terms of a slope-area (or in this case, slope-\(\lambda\)) relationship. Slope gradient is directly proportional to uplift rate, and inversely proportional to rainfall rate (\(rI\)). Slope gets smaller with more efficient transport (higher \(k_Q\)), more efficient plucking (\(k_p\)), or more efficient abrasion (\(\beta\)).

kQ = 0.041
beta = 0.0005
length = 1000.0
kp = 1.0e-4
gamma = 0.5
U = 1.0e-4
Lambda = 866025.404
intermittency = 0.01
r = 1.0

# alpha (bedrock exposure fraction)
term1 = kQ * (1.0 + beta * length / 2.0)
term2 = kp * gamma * length
alpha = term1 / (term1 + term2)
print("Predicted alpha:", alpha)

# Qs (gravel sediment transport rate out of cell)

# S (slope gradient)
S = (
    2 * U / (length * intermittency * r * alpha * (beta * kQ * (1.0 - alpha) + 2 * kp))
) ** (6.0 / 7.0)
print("Predicted S:", S)

# Qs (gravel sediment transport rate out of cell)
Qs = (kp * intermittency * r * S ** (7 / 6) * gamma * Lambda * length * alpha) / (
    1 + beta * length / 2
)
print("Predicted Qs:", Qs)

def test_steady_general():
    grid = HexModelGrid((3, 2), spacing=1000.0)
    grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
    grid.status_at_node[0] = grid.BC_NODE_IS_FIXED_VALUE

    elev = grid.add_zeros("topographic__elevation", at="node")
    elev[:] = (0.2 * grid.y_of_node) / np.cos(np.radians(30.0))
    sed = grid.add_zeros("soil__depth", at="node")
    sed[:] = 1.0
    rock = grid.add_zeros("bedrock__elevation", at="node")
    rock[:] = elev - sed

    fa = FlowAccumulator(grid)
    fa.run_one_step()
    gbe = GravelBedrockEroder(
        grid, abrasion_coefficient=0.0005, coarse_fraction_from_plucking=0.5
    )

    dt = 7500.0
    uplift_rate = 0.0001
    nsteps = 3000
    for i in range(nsteps):
        elev[grid.core_nodes] += uplift_rate * dt
        rock[grid.core_nodes] += uplift_rate * dt
        gbe.run_one_step(dt)

    assert_almost_equal(
        grid.at_node["bedrock__exposure_fraction"][grid.core_nodes], 0.5062, decimal=4
    )
    assert_almost_equal(
        grid.at_node["topographic__steepest_slope"][grid.core_nodes], 0.2387, decimal=4
    )
    assert_almost_equal(
        grid.at_node["bedload_sediment__volume_outflux"][grid.core_nodes],
        32.972,
        decimal=3,
    )

import time

start = time.time()
test_steady_general()
print(time.time() - start)

References and further reading

Attal, M., & Lavé, J. (2006). Changes of bedload characteristics along the Marsyandi River (central Nepal): Implications for understanding hillslope sediment supply, sediment load evolution along fluvial networks, and denudation in active orogenic belts. Geol. Soc. Am. Spec. Pap, 398, 143-171.

Attal, M., & Lavé, J. (2009). Pebble abrasion during fluvial transport: Experimental results and implications for the evolution of the sediment load along rivers. Journal of Geophysical Research: Earth Surface, 114(F4).

Meyer-Peter, E., & Müller, R. (1948). Formulas for bed-load transport. In IAHSR 2nd meeting, Stockholm, appendix 2. IAHR.

Parker, G. (1978). Self-formed straight rivers with equilibrium banks and mobile bed. Part 2. The gravel river. Journal of Fluid mechanics, 89(1), 127-146.

Phillips, C. B., & Jerolmack, D. J. (2016). Self-organization of river channels as a critical filter on climate signals. Science, 352(6286), 694-697.

Wickert, A. D., & Schildgen, T. F. (2019). Long-profile evolution of transport-limited gravel-bed rivers. Earth Surface Dynamics, 7(1), 17-43.

Willgoose, G., Bras, R. L., & Rodriguez‐Iturbe, I. (1991). A physical explanation of an observed link area‐slope relationship. Water Resources Research, 27(7), 1697-1702.

Willgoose, G. (1994). A physical explanation for an observed area‐slope‐elevation relationship for catchments with declining relief. Water Resources Research, 30(2), 151-159.

Wong, M., & Parker, G. (2006). Reanalysis and correction of bed-load relation of Meyer-Peter and Müller using their own database. Journal of Hydraulic Engineering, 132(11), 1159-1168.

---

Generated by nbsphinx from a Jupyter notebook.