landlab.components.gravel_bedrock_eroder.gravel_bedrock_eroder

Model bedrock incision and gravel transport and abrasion in a network of rivers.

@author: gtucker

class GravelBedrockEroder

Bases: Component

Drainage network evolution of rivers with gravel alluvium overlying bedrock.

Model drainage network evolution for a network of rivers that have a layer of gravel alluvium overlying bedrock.

GravelBedrockEroder is designed to operate together with a flow-routing component such as FlowAccumulator, 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.

As with the GravelRiverTransporter component, 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 alluvium thickness will be (Qin - Qs / (A * (1 - phi)), plus gravel produced by plucking erosion of bedrock (phi is porosity).

Bedrock is eroded by a combination of abrasion and plucking. Abrasion per unit channel length is calculated as the product of volumetric sediment discharge and an abrasion coefficient. Sediment produced by abrasion is assumed to go into wash load that is removed from the model domain. Plucking is calculated using a discharge-slope expression, and a user-defined fraction of plucked material is added to the coarse alluvium.

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

• depth_decay_scale (float (default 1.0)) – Scale for depth decay in bedrock exposure function

• plucking_coefficient (float or (n_core_nodes,) array of float (default 1.0e-4 1/m)) – Rate coefficient for bedrock erosion by plucking

• coarse_fraction_from_plucking (float or (n_core_nodes,) array of float (default 1.0)) – Fraction of plucked material that becomes part of gravel sediment load

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")
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[4] = 300.0
>>> 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()
>>> eroder = GravelBedrockEroder(grid, abrasion_coefficient=0.0005)
>>> rock_elev = grid.at_node["bedrock__elevation"]
>>> for _ in range(200):
...     rock_elev[grid.core_nodes] += 1.0
...     elev[grid.core_nodes] += 1.0
...     fa.run_one_step()
...     eroder.run_one_step(10000.0)
...
>>> int(elev[4] * 100)
2266

Initialize GravelBedrockEroder.

__init__(grid, intermittency_factor=0.01, transport_coefficient=0.041, abrasion_coefficient=0.0, sediment_porosity=0.35, depth_decay_scale=1.0, plucking_coefficient=0.0001, coarse_fraction_from_plucking=1.0)

Initialize GravelBedrockEroder.

static __new__(cls, *args, **kwds)

calc_abrasion_rate()

Update the volume 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 result is stored in the bedload_sediment__rate_of_loss_to_abrasion field.

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
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[3:] = 100.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> eroder = GravelBedrockEroder(grid, abrasion_coefficient=0.0002)
>>> eroder.calc_transport_rate()
>>> eroder.calc_abrasion_rate()
>>> int(eroder._abrasion[4] * 1e8)
19

calc_bedrock_abrasion_rate()

Update the rate of bedrock abrasion.

Note: assumes _abrasion (of sediment) and _rock_exposure_fraction have already been updated. Like _abrasion, the rate is a length per time (equivalent to rate of lowering of the bedrock surface by abrasion). Result is stored in the field bedrock__abrasion_rate.

>>> 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
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[:] = 1.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> eroder = GravelBedrockEroder(grid, abrasion_coefficient=1.0e-4)
>>> eroder.calc_rock_exposure_fraction()
>>> round(eroder._rock_exposure_fraction[6], 4)
0.3679
>>> eroder.calc_transport_rate()
>>> np.round(eroder._sediment_outflux[5:7], 3)
array([0.024, 0.012])
>>> eroder.calc_abrasion_rate()
>>> np.round(eroder._abrasion[5:7], 9)
array([1.2e-08, 6.0e-09])
>>> eroder.calc_bedrock_abrasion_rate()
>>> np.round(eroder._rock_abrasion_rate[5:7], 10)
array([4.4e-09, 2.2e-09])

calc_bedrock_plucking_rate()

Update the rate of bedrock erosion by plucking.

The rate is a volume per area per time [L/T], equivalent to the rate of lowering of the bedrock surface relative to the underlying material as a result of plucking. Result is stored in the field bedrock__plucking_rate.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> grid_res = 100.0
>>> grid = RasterModelGrid((3, 3), xy_spacing=grid_res)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[4] = 1.0
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> eroder = GravelBedrockEroder(grid)
>>> eroder.calc_rock_exposure_fraction()
>>> eroder.calc_bedrock_plucking_rate()
>>> predicted_plucking_rate = 1.0e-6 * 1.0e4 * 0.01 ** (7.0 / 6.0) / grid_res
>>> round(predicted_plucking_rate, 9)  # Kp Q S^(7/6)
4.64e-07
>>> int(round(eroder._pluck_rate[4] * 1e9))
464

calc_implied_depth(grain_diameter=0.01)

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
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[3:] = 100.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> eroder = GravelBedrockEroder(grid)
>>> water_depth = eroder.calc_implied_depth()
>>> int(water_depth[4] * 1000)
98

calc_implied_width(grain_diameter=0.01, time_unit='y')

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
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[3:] = 100.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> eroder = GravelBedrockEroder(grid)
>>> chan_width = eroder.calc_implied_width()
>>> int(chan_width[4] * 100)
3833
>>> grid.at_node["surface_water__discharge"] *= 1.0 / (3600 * 24 * 365.25)
>>> chan_width = eroder.calc_implied_width(time_unit="s")
>>> int(chan_width[4] * 100)
3838

calc_rock_exposure_fraction()

Update the bedrock exposure fraction.

The result is stored in the bedrock__exposure_fraction field.

Examples

>>> 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")
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[4] = 1000.0
>>> sed[5] = 0.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> eroder = GravelBedrockEroder(grid)
>>> eroder.calc_rock_exposure_fraction()
>>> eroder._rock_exposure_fraction[4:6]
array([0., 1.])
>>> sed[4] = 1.0  # exposure frac should be 1/e ~ 0.3679
>>> sed[5] = 2.0  # exposure frac should be 1/e^2 ~ 0.1353
>>> eroder.calc_rock_exposure_fraction()
>>> np.round(eroder._rock_exposure_fraction[4:6], 4)
array([0.3679, 0.1353])

calc_sediment_influx()

Update the volume influx at each node.

Result is stored in the field bedload_sediment__volume_influx.

calc_sediment_rate_of_change()

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

Result is stored in the field sediment__rate_of_change.

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
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[:] = 100.0
>>> 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()
>>> eroder = GravelBedrockEroder(grid)
>>> eroder.calc_transport_rate()
>>> eroder.calc_sediment_influx()
>>> eroder.calc_sediment_rate_of_change()
>>> np.round(eroder._sediment_outflux[4:7], 3)
array([0.   , 0.038, 0.019])
>>> np.round(eroder._sediment_influx[4:7], 3)
array([0.038, 0.019, 0.   ])
>>> np.round(eroder._dHdt[5:7], 8)
array([-2.93e-06, -2.93e-06])

calc_transport_rate()

Calculate and return bed-load transport rate.

Calculation uses Wickert-Schildgen approach, and provides volume per time rate. Transport rate is modulated by available sediment, using the exponential function (1 - exp(-H / Hs)), so that transport rate approaches zero as sediment thickness approaches zero. Rate is a volume per time. The result is stored in the bedload_sediment__volume_outflux field.

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
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[3:] = 100.0
>>> fa = FlowAccumulator(grid)
>>> fa.run_one_step()
>>> eroder = GravelBedrockEroder(grid)
>>> eroder.calc_transport_rate()
>>> round(eroder._sediment_outflux[4], 4)
0.019

cite_as = ''

property coords

Return the coordinates of nodes on grid attached to the component.

property current_time

Current time.

Some components may keep track of the current time. In this case, the current_time attribute is incremented. Otherwise it is set to None.

Return type:

current_time

definitions = (('bedload_sediment__rate_of_loss_to_abrasion', 'Rate of bedload sediment volume loss to abrasion per unit area'), ('bedload_sediment__volume_influx', 'Volumetric incoming streamwise bedload sediment transport rate'), ('bedload_sediment__volume_outflux', 'Volumetric outgoing streamwise bedload sediment transport rate'), ('bedrock__abrasion_rate', 'rate of bedrock lowering by abrasion'), ('bedrock__elevation', 'elevation of the bedrock surface'), ('bedrock__exposure_fraction', 'fractional exposure of bedrock'), ('bedrock__lowering_rate', 'Rate of lowering of bedrock surface'), ('bedrock__plucking_rate', 'rate of bedrock lowering by plucking'), ('flow__link_to_receiver_node', 'ID of link downstream of each node, which carries the discharge'), ('flow__receiver_node', 'Node array of receivers (node that receives flow from current node)'), ('flow__upstream_node_order', 'Node array containing downstream-to-upstream ordered list of node IDs'), ('sediment__rate_of_change', 'Time rate of change of sediment thickness'), ('soil__depth', 'Depth of soil or weathered bedrock'), ('surface_water__discharge', 'Volumetric discharge of surface water'), ('topographic__elevation', 'Land surface topographic elevation'), ('topographic__steepest_slope', 'The steepest *downhill* slope'))

classmethod from_path(grid, path)

Create a component from an input file.

Parameters:

• grid (ModelGrid) – A landlab grid.

• path (str or file_like) – Path to a parameter file, contents of a parameter file, or a file-like object.

Returns:

A newly-created component.

Return type:

Component

property grid

Return the grid attached to the component.

initialize_optional_output_fields()

Create fields for a component based on its optional field outputs, if declared in _optional_var_names.

This method will create new fields (without overwrite) for any fields output by the component as optional. New fields are initialized to zero. New fields are created as arrays of floats, unless the component also contains the specifying property _var_type.

initialize_output_fields(values_per_element=None)

Create fields for a component based on its input and output var names.

This method will create new fields (without overwrite) for any fields output by, but not supplied to, the component. New fields are initialized to zero. Ignores optional fields. New fields are created as arrays of floats, unless the component specifies the variable type.

Parameters:

values_per_element (int (optional)) – On occasion, it is necessary to create a field that is of size (n_grid_elements, values_per_element) instead of the default size (n_grid_elements,). Use this keyword argument to acomplish this task.

input_var_names = ('flow__link_to_receiver_node', 'flow__receiver_node', 'flow__upstream_node_order', 'soil__depth', 'surface_water__discharge', 'topographic__elevation', 'topographic__steepest_slope')

name = 'GravelBedrockEroder'

optional_var_names = ()

output_var_names = ('bedload_sediment__rate_of_loss_to_abrasion', 'bedload_sediment__volume_influx', 'bedload_sediment__volume_outflux', 'bedrock__abrasion_rate', 'bedrock__elevation', 'bedrock__exposure_fraction', 'bedrock__lowering_rate', 'bedrock__plucking_rate', 'sediment__rate_of_change', 'topographic__elevation')

run_one_step(global_dt)

Advance solution by time interval global_dt, subdividing into sub-steps as needed.

property shape

Return the grid shape attached to the component, if defined.

unit_agnostic = True

units = (('bedload_sediment__rate_of_loss_to_abrasion', 'm/y'), ('bedload_sediment__volume_influx', 'm**3/y'), ('bedload_sediment__volume_outflux', 'm**3/y'), ('bedrock__abrasion_rate', 'm/y'), ('bedrock__elevation', 'm'), ('bedrock__exposure_fraction', '-'), ('bedrock__lowering_rate', 'm/y'), ('bedrock__plucking_rate', 'm/y'), ('flow__link_to_receiver_node', '-'), ('flow__receiver_node', '-'), ('flow__upstream_node_order', '-'), ('sediment__rate_of_change', 'm/y'), ('soil__depth', 'm'), ('surface_water__discharge', 'm**3/y'), ('topographic__elevation', 'm'), ('topographic__steepest_slope', '-'))

update_rates()

Update rate of sediment thickness change, and rate of bedrock lowering by abrasion and plucking.

Combined rate of rock lowering relative to underlying material is stored in the field bedrock__lowering_rate.

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
>>> sed = grid.add_zeros("soil__depth", at="node")
>>> sed[:] = 1000.0
>>> 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()
>>> eroder = GravelBedrockEroder(grid)
>>> eroder.run_one_step(1000.0)
>>> np.round(elev[4:7], 4)
array([0.    , 0.9971, 1.9971])

classmethod var_definition(name)

Get a description of a particular field.

Parameters:

name (str) – A field name.

Returns:

A description of each field.

Return type:

tuple of (name, *description*)

classmethod var_help(name)

Print a help message for a particular field.

Parameters:

name (str) – A field name.

classmethod var_loc(name)

Location where a particular variable is defined.

Parameters:

name (str) – A field name.

Returns:

The location (‘node’, ‘link’, etc.) where a variable is defined.

Return type:

str

var_mapping = (('bedload_sediment__rate_of_loss_to_abrasion', 'node'), ('bedload_sediment__volume_influx', 'node'), ('bedload_sediment__volume_outflux', 'node'), ('bedrock__abrasion_rate', 'node'), ('bedrock__elevation', 'node'), ('bedrock__exposure_fraction', 'node'), ('bedrock__lowering_rate', 'node'), ('bedrock__plucking_rate', 'node'), ('flow__link_to_receiver_node', 'node'), ('flow__receiver_node', 'node'), ('flow__upstream_node_order', 'node'), ('sediment__rate_of_change', 'node'), ('soil__depth', 'node'), ('surface_water__discharge', 'node'), ('topographic__elevation', 'node'), ('topographic__steepest_slope', 'node'))

classmethod var_type(name)

Returns the dtype of a field (float, int, bool, str…).

Parameters:

name (str) – A field name.

Returns:

The dtype of the field.

Return type:

dtype

classmethod var_units(name)

Get the units of a particular field.

Parameters:

name (str) – A field name.

Returns:

Units for the given field.

Return type:

str