GravelRiverTransporter: bed-load transport and downstream abrasion in a network of gravel rivers#

class GravelRiverTransporter(*args, **kwds)[source]#

Bases: Component

Model drainage network evolution for a network of transport-limited gravel-bed rivers with downstream abrasion.

GravelRiverTransporter is designed to operate together with a flow-routing component such as PriorityFloodFlowRouter, so that each grid node has a defined flow direction toward one of its neighbor nodes. Each core node is assumed to contain one outgoing fluvial channel, and (depending on the drainage structure) zero, one, or more incoming channels. These channels are treated as effectively sub-grid-scale features that are embedded in valleys that have a width of one grid cell. The rate of gravel transport out of a given node is calculated as the product of bankfull discharge, channel gradient (to the 7/6 power), a dimensionless transport coefficient, and an intermittency factor that represents the fraction of time that bankfull flow occurs. The derivation of the transport law is given by Wickert & Schildgen (2019), and it derives from the assumption that channels are gravel-bedded and that they “instantaneously” adjust their width such that bankfull bed shear stress is just slightly higher than the threshold for grain motion. The substrate is assumed to consist entirely of gravel-size material with a given bulk porosity. The component calculates the loss of gravel-sized material to abrasion (i.e., conversion to finer sediment, which is not explicitly tracked) as a function of the volumetric transport rate, an abrasion coefficient with units of inverse length, and the local transport distance (for example, if a grid node is carrying a gravel load Qs to a neighboring node dx meters downstream, the rate of gravel loss in volume per time per area at the node will be beta Qs dx, where beta is the abrasion coefficient). Sediment mass conservation is calculated across each entire grid cell. For example, if a cell has surface area A, a total volume influx Qin, and downstream transport rate Qs, the resulting rate of change of elevation will be (Qin - Qs / (A (1 - phi)), where phi is porosity.

Parameters:
  • grid (ModelGrid) – A Landlab model grid object

  • intermittency_factor (float (default 0.01)) – Fraction of time that bankfull flow occurs

  • transport_coefficient (float (default 0.041)) – Dimensionless transport efficiency factor (see Wickert & Schildgen 2019)

  • abrasion_coefficient (float (default 0.0 1/m)) – Abrasion coefficient with units of inverse length

  • sediment_porosity (float (default 0.35)) – Bulk porosity of bed sediment

  • solver (string (default "explicit")) – Solver type (currently only “explicit” is tested and operational)

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 3), xy_spacing=1000.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[5] = grid.BC_NODE_IS_FIXED_VALUE
>>> fa = FlowAccumulator(grid, runoff_rate=10.0)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid, abrasion_coefficient=0.0005)
>>> for _ in range(200):
...     fa.run_one_step()
...     elev[grid.core_nodes] += 1.0
...     transporter.run_one_step(10000.0)
...
>>> int(elev[4] * 100)
2366

Initialize GravelRiverTransporter.

__init__(grid, intermittency_factor=0.01, transport_coefficient=0.041, abrasion_coefficient=0.0, sediment_porosity=0.35, solver='explicit')[source]#

Initialize GravelRiverTransporter.

calc_abrasion_rate()[source]#

Update the rate of bedload loss to abrasion, per unit area.

Here we use the average of incoming and outgoing sediment flux to calculate the loss rate to abrasion.

The factor dx (node spacing) appears in the denominator to represent flow segment length (i.e., length of the link along which water is flowing in the cell) divided by cell area. This would need to be updated to handle non-raster and/or non-uniform grids.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 3), xy_spacing=1000.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[3:] = 10.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid, abrasion_coefficient=0.0002)
>>> transporter.calc_transport_capacity()
>>> transporter.calc_abrasion_rate()
>>> int(transporter._abrasion[4] * 1e8)
19
calc_implied_depth(grain_diameter=0.01)[source]#

Utility function that calculates and returns water depth implied by slope and grain diameter, using Wickert & Schildgen (2019) equation 8.

The equation is

h = ((rho_s - rho / rho)) (1 + epsilon) tau_c* (D / S)

where the factors on the right are sediment and water density, excess shear-stress factor, critical Shields stress, grain diameter, and slope gradient. Here the prefactor on D/S assumes sediment density of 2650 kg/m3, water density of 1000 kg/m3, shear-stress factor of 0.2, and critical Shields stress of 0.0495, giving a value of 0.09801.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 3), xy_spacing=1000.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[3:] = 10.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid)
>>> depth = transporter.calc_implied_depth()
>>> int(depth[4] * 1000)
98
calc_implied_width(grain_diameter=0.01, time_unit='y')[source]#

Utility function that calculates and returns channel width implied by discharge, slope, and grain diameter, using Wickert & Schildgen (2019) equation 16.

The equation is

b = kb Q S**(7/6) / D**(3/2)

where the dimensional prefactor, which includes sediment and water density, gravitational acceleration, critical Shields stress, and the transport factor epsilon, is

kb = 0.17 g**(-1/2) (((rho_s - rho) / rho) (1 + eps) tau_c*)**(-5/3)

Using g = 9.8 m/s2, rho_s = 2650 (quartz), rho = 1000 kg/m3, eps = 0.2, and tau_c* = 0.0495, kb ~ 2.61 s/m**(1/2). Converting to years, kb = 8.26e-8.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 3), xy_spacing=10000.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[3:] = 100.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid)
>>> width = transporter.calc_implied_width()
>>> int(width[4] * 100)
3833
>>> grid.at_node["surface_water__discharge"] *= 1.0 / (3600 * 24 * 365.25)
>>> width = transporter.calc_implied_width(time_unit="s")
>>> int(width[4] * 100)
3838
calc_sediment_rate_of_change()[source]#

Update the rate of thickness change of coarse sediment at each core node.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 4), xy_spacing=100.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[:] = 0.01 * grid.x_of_node
>>> grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[4] = grid.BC_NODE_IS_FIXED_VALUE
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid)
>>> transporter.calc_sediment_rate_of_change()
>>> np.round(transporter._sediment_outflux[4:7], 3)
array([ 0.   ,  0.038,  0.019])
>>> np.round(transporter._sediment_influx[4:7], 3)
array([ 0.038,  0.019,  0.   ])
>>> np.round(transporter._dzdt[5:7], 8)
array([ -2.93000000e-06,  -2.93000000e-06])
calc_transport_capacity()[source]#

Calculate and return bed-load transport capacity.

Calculation uses Wickert-Schildgen approach, and provides volume per time rate.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 3), xy_spacing=100.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[3:] = 1.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid)
>>> transporter.calc_transport_capacity()
>>> round(transporter._sediment_outflux[4], 4)
0.019
run_one_step_matrix_inversion(dt)[source]#

Advance solution by time interval dt.

WARNING: EXPERIMENTAL AND NOT FULLY TESTED - USE AT OWN RISK!

Notes

Does not update abrasion rate or sediment outflux fields.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 4), xy_spacing=100.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[:] = 0.01 * grid.x_of_node
>>> grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[4] = grid.BC_NODE_IS_FIXED_VALUE
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid, solver="matrix")
>>> transporter.run_one_step == transporter.run_one_step_matrix_inversion
True
>>> transporter.run_one_step(1000.0)
run_one_step_simple_explicit(dt)[source]#

Advance solution by time interval dt.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid = RasterModelGrid((3, 4), xy_spacing=100.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[:] = 0.01 * grid.x_of_node
>>> grid.status_at_node[grid.perimeter_nodes] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[4] = grid.BC_NODE_IS_FIXED_VALUE
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> transporter = GravelRiverTransporter(grid, solver="explicit")
>>> transporter.run_one_step(1000.0)
>>> np.round(elev[4:7], 4)
array([ 0.    ,  0.9971,  1.9971])
make_empty_matrix_and_rhs(grid)[source]#
zero_out_matrix(grid, mat, rcvr, mat_id)[source]#