ErosionDeposition: Fluvial erosion in the style of Davy and Lague (2009)#

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

Bases: _GeneralizedErosionDeposition

Erosion-Deposition model in the style of Davy and Lague (2009). It uses a mass balance approach across the total sediment mass both in the bed and in transport coupled with explicit representation of the sediment transport lengthscale (the “xi-q” model) to derive a range of erosional and depositional responses in river channels.

This implementation is close to the Davy & Lague scheme, with a few deviations:

  • A fraction of the eroded sediment is permitted to enter the wash load, and lost to the mass balance (F_f).

  • Here an incision threshold \(\omega\) is permitted, where it was not by Davy & Lague. It is implemented with an exponentially smoothed form to prevent discontinuities in the parameter space. See the StreamPowerSmoothThresholdEroder for more documentation.

  • This component uses an “effective” settling velocity, v_s, as one of its inputs. This parameter is simply equal to Davy & Lague’s d_star * V dimensionless number.

Erosion of the bed follows a stream power formulation, i.e.,

Note that the transition between transport-limited and detachment-limited behavior is controlled by the dimensionless ratio (v_s/r) where r is the runoff ratio (Q=Ar). r can be changed in the flow accumulation component but is not changed within ErosionDeposition. Because the runoff ratio r is not changed within the ErosionDeposition component, v_s becomes the parameter that fundamentally controls response style. Very small v_s will lead to a detachment-limited response style, very large v_s will lead to a transport-limited response style. v_s == 1 means equal contributions from transport and erosion, and a hybrid response as described by Davy & Lague.

Unlike other some other fluvial erosion componets in Landlab, in this component (and SPACE) no erosion occurs in depressions or in areas with adverse slopes. There is no ability to pass a keyword argument erode_flooded_nodes.

If a depressions are handled (as indicated by the presence of the field “flood_status_code” at nodes), then deposition occurs throughout the depression and sediment is passed out of the depression. Where pits are encountered, then all sediment is deposited at that node only.

A note about sediment porosity: Prior to Landlab v2.0 this component took a porositiy keyworkd argument phi. For an explaination of why it no longer does (including a mathematical derivation), see Pull Request 1186. If phi is passed to this component a value error will be raised.

Component written by C. Shobe, K. Barnhart, and G. Tucker.

References

Required Software Citation(s) Specific to this Component

Barnhart, K., Glade, R., Shobe, C., Tucker, G. (2019). Terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution. Geoscientific Model Development 12(4), 1267–1297. https://dx.doi.org/10.5194/gmd-12-1267-2019

Additional References

Davy, P., Lague, D. (2009). Fluvial erosion/transport equation of landscape evolution models revisited Journal of Geophysical Research 114(F3), F03007. https://dx.doi.org/10.1029/2008jf001146

Initialize the ErosionDeposition model.

Parameters:
  • grid (ModelGrid) – Landlab ModelGrid object

  • K (float, field name, or array) – Erodibility for substrate (units vary).

  • v_s (float) – Effective settling velocity for chosen grain size metric [L/T].

  • m_sp (float) – Discharge exponent (units vary)

  • n_sp (float) – Slope exponent (units vary)

  • sp_crit (float, field name, or array) – Critical stream power to erode substrate [E/(TL^2)]

  • F_f (float) – Fraction of eroded material that turns into “fines” that do not contribute to (coarse) sediment load. Defaults to zero.

  • discharge_field (float, field name, or array) – Discharge [L^2/T]. The default is to use the grid field ‘surface_water__discharge’, which is simply drainage area multiplied by the default rainfall rate (1 m/yr). To use custom spatially/temporally varying rainfall, use ‘water__unit_flux_in’ to specify water input to the FlowAccumulator.

  • solver (string) –

    Solver to use. Options at present include:
    1. ’basic’ (default): explicit forward-time extrapolation. Simple but will become unstable if time step is too large.

    2. ’adaptive’: adaptive time-step solver that estimates a stable step size based on the shortest time to “flattening” among all upstream-downstream node pairs.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import DepressionFinderAndRouter
>>> from landlab.components import ErosionDeposition
>>> from landlab.components import FastscapeEroder
>>> np.random.seed(seed=5000)
Define grid and initial topography:

-5x5 grid with baselevel in the lower left corner -all other boundary nodes closed -Initial topography is plane tilted up to the upper right + noise

>>> nr = 5
>>> nc = 5
>>> dx = 10
>>> mg = RasterModelGrid((nr, nc), xy_spacing=10.0)
>>> _ = mg.add_zeros("node", "topographic__elevation")
>>> mg.at_node["topographic__elevation"] += (
...     mg.node_y / 10 + mg.node_x / 10 + np.random.rand(len(mg.node_y)) / 10
... )
>>> mg.set_closed_boundaries_at_grid_edges(
...     bottom_is_closed=True,
...     left_is_closed=True,
...     right_is_closed=True,
...     top_is_closed=True,
... )
>>> mg.set_watershed_boundary_condition_outlet_id(
...     0, mg.at_node["topographic__elevation"], -9999.0
... )
>>> fsc_dt = 100.0
>>> ed_dt = 1.0

Check initial topography

>>> mg.at_node["topographic__elevation"].reshape(mg.shape)
array([[ 0.02290479,  1.03606698,  2.0727653 ,  3.01126678,  4.06077707],
       [ 1.08157495,  2.09812694,  3.00637448,  4.07999597,  5.00969486],
       [ 2.04008677,  3.06621577,  4.09655859,  5.04809001,  6.02641123],
       [ 3.05874171,  4.00585786,  5.0595697 ,  6.04425233,  7.05334077],
       [ 4.05922478,  5.0409473 ,  6.07035008,  7.0038935 ,  8.01034357]])

Instantiate Fastscape eroder, flow router, and depression finder

>>> fr = FlowAccumulator(mg, flow_director="D8")
>>> df = DepressionFinderAndRouter(mg)
>>> fsc = FastscapeEroder(mg, K_sp=0.001, m_sp=0.5, n_sp=1)

Burn in an initial drainage network using the Fastscape eroder:

>>> for x in range(100):
...     fr.run_one_step()
...     df.map_depressions()
...     flooded = np.where(df.flood_status == 3)[0]
...     fsc.run_one_step(dt=fsc_dt)
...     mg.at_node["topographic__elevation"][0] -= 0.001  # uplift
...

Instantiate the E/D component:

>>> ed = ErosionDeposition(
...     mg, K=0.00001, v_s=0.001, m_sp=0.5, n_sp=1.0, sp_crit=0
... )

Now run the E/D component for 2000 short timesteps:

>>> for x in range(2000):  # E/D component loop
...     fr.run_one_step()
...     df.map_depressions()
...     ed.run_one_step(dt=ed_dt)
...     mg.at_node["topographic__elevation"][0] -= 2e-4 * ed_dt
...

Now we test to see if topography is right:

>>> np.around(mg.at_node["topographic__elevation"], decimals=3).reshape(
...     mg.shape
... )
array([[-0.477,  1.036,  2.073,  3.011,  4.061],
       [ 1.082, -0.08 , -0.065, -0.054,  5.01 ],
       [ 2.04 , -0.065, -0.065, -0.053,  6.026],
       [ 3.059, -0.054, -0.053, -0.035,  7.053],
       [ 4.059,  5.041,  6.07 ,  7.004,  8.01 ]])
property K#

Erodibility (units depend on m_sp).

__init__(grid, K=0.002, v_s=1.0, m_sp=0.5, n_sp=1.0, sp_crit=0.0, F_f=0.0, discharge_field='surface_water__discharge', solver='basic', dt_min=0.001, **kwds)[source]#

Initialize the ErosionDeposition model.

Parameters:
  • grid (ModelGrid) – Landlab ModelGrid object

  • K (float, field name, or array) – Erodibility for substrate (units vary).

  • v_s (float) – Effective settling velocity for chosen grain size metric [L/T].

  • m_sp (float) – Discharge exponent (units vary)

  • n_sp (float) – Slope exponent (units vary)

  • sp_crit (float, field name, or array) – Critical stream power to erode substrate [E/(TL^2)]

  • F_f (float) – Fraction of eroded material that turns into “fines” that do not contribute to (coarse) sediment load. Defaults to zero.

  • discharge_field (float, field name, or array) – Discharge [L^2/T]. The default is to use the grid field ‘surface_water__discharge’, which is simply drainage area multiplied by the default rainfall rate (1 m/yr). To use custom spatially/temporally varying rainfall, use ‘water__unit_flux_in’ to specify water input to the FlowAccumulator.

  • solver (string) –

    Solver to use. Options at present include:
    1. ’basic’ (default): explicit forward-time extrapolation. Simple but will become unstable if time step is too large.

    2. ’adaptive’: adaptive time-step solver that estimates a stable step size based on the shortest time to “flattening” among all upstream-downstream node pairs.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import DepressionFinderAndRouter
>>> from landlab.components import ErosionDeposition
>>> from landlab.components import FastscapeEroder
>>> np.random.seed(seed=5000)
Define grid and initial topography:

-5x5 grid with baselevel in the lower left corner -all other boundary nodes closed -Initial topography is plane tilted up to the upper right + noise

>>> nr = 5
>>> nc = 5
>>> dx = 10
>>> mg = RasterModelGrid((nr, nc), xy_spacing=10.0)
>>> _ = mg.add_zeros("node", "topographic__elevation")
>>> mg.at_node["topographic__elevation"] += (
...     mg.node_y / 10 + mg.node_x / 10 + np.random.rand(len(mg.node_y)) / 10
... )
>>> mg.set_closed_boundaries_at_grid_edges(
...     bottom_is_closed=True,
...     left_is_closed=True,
...     right_is_closed=True,
...     top_is_closed=True,
... )
>>> mg.set_watershed_boundary_condition_outlet_id(
...     0, mg.at_node["topographic__elevation"], -9999.0
... )
>>> fsc_dt = 100.0
>>> ed_dt = 1.0

Check initial topography

>>> mg.at_node["topographic__elevation"].reshape(mg.shape)
array([[ 0.02290479,  1.03606698,  2.0727653 ,  3.01126678,  4.06077707],
       [ 1.08157495,  2.09812694,  3.00637448,  4.07999597,  5.00969486],
       [ 2.04008677,  3.06621577,  4.09655859,  5.04809001,  6.02641123],
       [ 3.05874171,  4.00585786,  5.0595697 ,  6.04425233,  7.05334077],
       [ 4.05922478,  5.0409473 ,  6.07035008,  7.0038935 ,  8.01034357]])

Instantiate Fastscape eroder, flow router, and depression finder

>>> fr = FlowAccumulator(mg, flow_director="D8")
>>> df = DepressionFinderAndRouter(mg)
>>> fsc = FastscapeEroder(mg, K_sp=0.001, m_sp=0.5, n_sp=1)

Burn in an initial drainage network using the Fastscape eroder:

>>> for x in range(100):
...     fr.run_one_step()
...     df.map_depressions()
...     flooded = np.where(df.flood_status == 3)[0]
...     fsc.run_one_step(dt=fsc_dt)
...     mg.at_node["topographic__elevation"][0] -= 0.001  # uplift
...

Instantiate the E/D component:

>>> ed = ErosionDeposition(
...     mg, K=0.00001, v_s=0.001, m_sp=0.5, n_sp=1.0, sp_crit=0
... )

Now run the E/D component for 2000 short timesteps:

>>> for x in range(2000):  # E/D component loop
...     fr.run_one_step()
...     df.map_depressions()
...     ed.run_one_step(dt=ed_dt)
...     mg.at_node["topographic__elevation"][0] -= 2e-4 * ed_dt
...

Now we test to see if topography is right:

>>> np.around(mg.at_node["topographic__elevation"], decimals=3).reshape(
...     mg.shape
... )
array([[-0.477,  1.036,  2.073,  3.011,  4.061],
       [ 1.082, -0.08 , -0.065, -0.054,  5.01 ],
       [ 2.04 , -0.065, -0.065, -0.053,  6.026],
       [ 3.059, -0.054, -0.053, -0.035,  7.053],
       [ 4.059,  5.041,  6.07 ,  7.004,  8.01 ]])
run_one_step_basic(dt=1.0)[source]#

Calculate change in rock and alluvium thickness for a time period ‘dt’.

Parameters:

dt (float) – Model timestep [T]

run_with_adaptive_time_step_solver(dt=1.0)[source]#

CHILD-like solver that adjusts time steps to prevent slope flattening.

Parameters:

dt (float) – Model timestep [T]