Stream Power with Alluvium Conservation and Entrainment#

Space: Stream Power with Alluvium Conservation and Entrainment#

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

Bases: _GeneralizedErosionDeposition

Stream Power with Alluvium Conservation and Entrainment (SPACE)

See the publication:

Shobe, C. M., Tucker, G. E., and Barnhart, K. R.: The SPACE 1.0 model: a Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution, Geosci. Model Dev., 10, 4577-4604, https://doi.org/10.5194/gmd-10-4577-2017, 2017.

Unlike other some other fluvial erosion componets in Landlab, in this component (and ErosionDeposition) 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.

Note: If timesteps are large enough that Es*dt (sediment erosion) exceeds sediment thickness H, the ‘adaptive’ solver is necessary to subdivide timesteps. Compare Es and H arrays to determine whether timesteps are appropriate or too large for the ‘basic’ solver.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import (
...     FlowAccumulator,
...     DepressionFinderAndRouter,
...     Space,
...     FastscapeEroder,
... )
>>> np.random.seed(seed=5000)

Define grid and initial topography:

  • 5x5 grid with base level in the lower left corner

  • All other boundary nodes closed

  • Initial topography is plane tilted up to the upper right with noise

>>> mg = RasterModelGrid((5, 5), xy_spacing=10.0)
>>> _ = mg.add_zeros("topographic__elevation", at="node")
>>> mg.at_node["topographic__elevation"] += (
...     mg.node_y / 10.0 + mg.node_x / 10.0 + np.random.rand(len(mg.node_y)) / 10.0
... )
>>> 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
>>> space_dt = 100.0

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 _ in range(100):
...     fr.run_one_step()
...     df.map_depressions()
...     fsc.run_one_step(dt=fsc_dt)
...     mg.at_node["topographic__elevation"][0] -= 0.001  # Uplift
...

Add some soil to the drainage network:

>>> _ = mg.add_zeros("soil__depth", at="node", dtype=float)
>>> mg.at_node["soil__depth"] += 0.5
>>> mg.at_node["topographic__elevation"] += mg.at_node["soil__depth"]

Instantiate the Space component:

>>> ha = Space(
...     mg,
...     K_sed=0.00001,
...     K_br=0.00000000001,
...     F_f=0.5,
...     phi=0.1,
...     H_star=1.0,
...     v_s=0.001,
...     m_sp=0.5,
...     n_sp=1.0,
...     sp_crit_sed=0,
...     sp_crit_br=0,
... )

Now run the Space component for 2000 short timesteps:

>>> for _ in range(2000):  # Space component loop
...     fr.run_one_step()
...     df.map_depressions()
...     ha.run_one_step(dt=space_dt)
...     mg.at_node["bedrock__elevation"][0] -= 2e-6 * space_dt
...

Now we test to see if soil depth and topography are right:

>>> np.around(mg.at_node["soil__depth"], decimals=3)
array([ 0.5  ,  0.5  ,  0.5  ,  0.5  ,  0.5  ,  0.5  ,  0.495,  0.492,
        0.491,  0.5  ,  0.5  ,  0.492,  0.492,  0.49 ,  0.5  ,  0.5  ,
        0.491,  0.49 ,  0.484,  0.5  ,  0.5  ,  0.5  ,  0.5  ,  0.5  ,
        0.5  ])

>>> np.around(mg.at_node["topographic__elevation"], decimals=3)
array([ 0.423,  1.536,  2.573,  3.511,  4.561,  1.582,  0.424,  0.428,
        0.438,  5.51 ,  2.54 ,  0.428,  0.428,  0.438,  6.526,  3.559,
        0.438,  0.438,  0.45 ,  7.553,  4.559,  5.541,  6.57 ,  7.504,
        8.51 ])

References

Required Software Citation(s) Specific to this Component

Shobe, C., Tucker, G., Barnhart, K. (2017). The SPACE 1.0 model: a Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution. Geoscientific Model Development 10(12), 4577 - 4604. https://dx.doi.org/10.5194/gmd-10-4577-2017

Additional References

None Listed

Initialize the Space model.

Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • K_sed (float, field name, or array) – Erodibility for sediment (units vary).

  • K_br (float, field name, or array) – Erodibility for bedrock (units vary).

  • F_f (float) – Fraction of permanently suspendable fines in bedrock [-].

  • phi (float) – Sediment porosity [-].

  • H_star (float) – Sediment thickness required for full entrainment [L].

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

  • m_sp (float) – Drainage area exponent (units vary)

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

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

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

  • 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’: subdivides global time step as needed to prevent slopes from reversing and alluvium from going negative.

property Er#

Bedrock erosion term.

property Es#

Sediment erosion term.

property H#

Sediment thickness.

property K_br#

Erodibility of bedrock(units depend on m_sp).

property K_sed#

Erodibility of sediment(units depend on m_sp).

__init__(grid, K_sed=0.002, K_br=0.002, F_f=0.0, phi=0.3, H_star=0.1, v_s=1.0, m_sp=0.5, n_sp=1.0, sp_crit_sed=0.0, sp_crit_br=0.0, discharge_field='surface_water__discharge', solver='basic', dt_min=0.001)[source]#

Initialize the Space model.

Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • K_sed (float, field name, or array) – Erodibility for sediment (units vary).

  • K_br (float, field name, or array) – Erodibility for bedrock (units vary).

  • F_f (float) – Fraction of permanently suspendable fines in bedrock [-].

  • phi (float) – Sediment porosity [-].

  • H_star (float) – Sediment thickness required for full entrainment [L].

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

  • m_sp (float) – Drainage area exponent (units vary)

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

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

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

  • 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’: subdivides global time step as needed to prevent slopes from reversing and alluvium from going negative.

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]#

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

Parameters:

dt (float) – Model timestep [T]

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> import numpy as np

>>> rg = RasterModelGrid((3, 4))
>>> z = rg.add_zeros("topographic__elevation", at="node")
>>> z[:] = 0.1 * rg.x_of_node
>>> H = rg.add_zeros("soil__depth", at="node")
>>> H += 0.1
>>> br = rg.add_zeros("bedrock__elevation", at="node")
>>> br[:] = z - H

>>> fa = FlowAccumulator(rg, flow_director="FlowDirectorSteepest")
>>> fa.run_one_step()
>>> sp = Space(
...     rg,
...     K_sed=1.0,
...     K_br=0.1,
...     F_f=0.5,
...     phi=0.0,
...     H_star=1.0,
...     v_s=1.0,
...     m_sp=0.5,
...     n_sp=1.0,
...     sp_crit_sed=0,
...     sp_crit_br=0,
...     solver="adaptive",
... )
>>> sp.run_one_step(dt=10.0)

>>> np.round(sp.Es[5:7], 4)
array([ 0.0029,  0.0074])
>>> np.round(sp.Er[5:7], 4)
array([ 0.0032,  0.0085])
>>> np.round(H[5:7], 3)
array([ 0.088,  0.078])

SpaceLargeScaleEroder: SPACE large-scale eroder#

Grid-based simulation of lateral erosion by channels in a drainage network.

Benjamin Campforts

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

Bases: Component

Stream Power with Alluvium Conservation and Entrainment (SPACE) large scale eroder

The SPACE_large_Scale_eroder is based on the SPACE component and is designed to be more robust against large time steps and coded in such a way that mass conservation is explicitly conserved during calculation.

See the publication: Shobe, C. M., Tucker, G. E., and Barnhart, K. R.: The SPACE 1.0 model: a Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution, Geosci. Model Dev., 10, 4577-4604, https://doi.org/10.5194/gmd-10-4577-2017, 2017.

Unlike other some other fluvial erosion componets in Landlab, in this component (and ErosionDeposition) 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.

Note: In the current version, we do not provide an adaptive time stepper. This will be addded in future versions of this component.

For more explanation and examples, check out the correponding notebook of this component

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import PriorityFloodFlowRouter, SpaceLargeScaleEroder
>>> import matplotlib.pyplot as plt  # For plotting results; optional
>>> from landlab import imshow_grid  # For plotting results; optional

>>> num_rows = 20
>>> num_columns = 20
>>> node_spacing = 100.0
>>> mg = RasterModelGrid((num_rows, num_columns), xy_spacing=node_spacing)
>>> node_next_to_outlet = num_columns + 1
>>> np.random.seed(seed=5000)
>>> _ = mg.add_zeros("topographic__elevation", at="node")
>>> _ = mg.add_zeros("soil__depth", at="node")
>>> mg.at_node["soil__depth"][mg.core_nodes] = 2.0
>>> _ = mg.add_zeros("bedrock__elevation", at="node")
>>> mg.at_node["bedrock__elevation"] += (
...     mg.node_y / 10.0 + mg.node_x / 10.0 + np.random.rand(len(mg.node_y)) / 10.0
... )
>>> mg.at_node["bedrock__elevation"][:] = mg.at_node["topographic__elevation"]
>>> mg.at_node["topographic__elevation"][:] += mg.at_node["soil__depth"]
>>> 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
... )

>>> fr = PriorityFloodFlowRouter(mg, flow_metric="D8", suppress_out=True)
>>> sp = SpaceLargeScaleEroder(
...     mg,
...     K_sed=0.01,
...     K_br=0.001,
...     F_f=0.0,
...     phi=0.0,
...     H_star=1.0,
...     v_s=5.0,
...     m_sp=0.5,
...     n_sp=1.0,
...     sp_crit_sed=0,
...     sp_crit_br=0,
... )
>>> timestep = 10.0
>>> elapsed_time = 0.0
>>> count = 0
>>> run_time = 1e4
>>> sed_flux = np.zeros(int(run_time // timestep))
>>> while elapsed_time < run_time:
...     fr.run_one_step()
...     _ = sp.run_one_step(dt=timestep)
...     sed_flux[count] = mg.at_node["sediment__flux"][node_next_to_outlet]
...     elapsed_time += timestep
...     count += 1
...

Plot the results.

>>> fig = plt.figure()
>>> plot = plt.subplot()
>>> _ = imshow_grid(
...     mg,
...     "topographic__elevation",
...     plot_name="Sediment flux",
...     var_name="Sediment flux",
...     var_units=r"m$^3$/yr",
...     grid_units=("m", "m"),
...     cmap="terrain",
... )
>>> _ = plt.figure()
>>> _ = imshow_grid(
...     mg,
...     "sediment__flux",
...     plot_name="Sediment flux",
...     var_name="Sediment flux",
...     var_units=r"m$^3$/yr",
...     grid_units=("m", "m"),
...     cmap="terrain",
... )
>>> fig = plt.figure()
>>> sedfluxplot = plt.subplot()
>>> _ = sedfluxplot.plot(
...     np.arange(len(sed_flux)) * timestep, sed_flux, color="k", linewidth=1.0
... )
>>> _ = sedfluxplot.set_xlabel("Time [yr]")
>>> _ = sedfluxplot.set_ylabel(r"Sediment flux [m$^3$/yr]")

References

Required Software Citation(s) Specific to this Component

Shobe, C., Tucker, G., Barnhart, K. (2017). The SPACE 1.0 model: a Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution. Geoscientific Model Development 10(12), 4577 - 4604. https://dx.doi.org/10.5194/gmd-10-4577-2017

Additional References

None Listed

Initialize the SpaceLargeScaleEroder model.

Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • K_sed (float, array of float, or str, optional) – Erodibility for sediment (units vary) as either a number or a field name.

  • K_br (float, array of float, or str, optional) – Erodibility for bedrock (units vary) as either a number or a field name.

  • F_f (float, optional) – Fraction of permanently suspendable fines in bedrock [-].

  • phi (float, optional) – Sediment porosity [-].

  • H_star (float, optional) – Sediment thickness required for full entrainment [L].

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

  • v_s_lake (float, optional) – Effective settling velocity in lakes for chosen grain size metric [L/T].

  • m_sp (float, optional) – Drainage area exponent (units vary).

  • n_sp (float, optional) – Slope exponent (units vary).

  • sp_crit_sed (float, array of float, or str, optional) – Critical stream power to erode sediment [E/(TL^2)].

  • sp_crit_br (float, array of float, or str, optional) – Critical stream power to erode rock [E/(TL^2)]

  • discharge_field (float, array of float, or str, optional) – 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.

  • erode_flooded_nodes (bool, optional) – Whether erosion occurs in flooded nodes identified by a depression/lake mapper (e.g., DepressionFinderAndRouter). When set to false, the field flood_status_code must be present on the grid (this is created by the DepressionFinderAndRouter). Default True.

property Er#

Bedrock erosion term.

property Es#

Sediment erosion term.

property K_br#

Erodibility of bedrock(units depend on m_sp).

property K_sed#

Erodibility of sediment(units depend on m_sp).

__init__(grid, K_sed=0.02, K_br=0.02, F_f=0.0, phi=0.3, H_star=0.1, v_s=1.0, v_s_lake=None, m_sp=0.5, n_sp=1.0, sp_crit_sed=0.0, sp_crit_br=0.0, discharge_field='surface_water__discharge', erode_flooded_nodes=False, thickness_lim=100)[source]#

Initialize the SpaceLargeScaleEroder model.

Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • K_sed (float, array of float, or str, optional) – Erodibility for sediment (units vary) as either a number or a field name.

  • K_br (float, array of float, or str, optional) – Erodibility for bedrock (units vary) as either a number or a field name.

  • F_f (float, optional) – Fraction of permanently suspendable fines in bedrock [-].

  • phi (float, optional) – Sediment porosity [-].

  • H_star (float, optional) – Sediment thickness required for full entrainment [L].

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

  • v_s_lake (float, optional) – Effective settling velocity in lakes for chosen grain size metric [L/T].

  • m_sp (float, optional) – Drainage area exponent (units vary).

  • n_sp (float, optional) – Slope exponent (units vary).

  • sp_crit_sed (float, array of float, or str, optional) – Critical stream power to erode sediment [E/(TL^2)].

  • sp_crit_br (float, array of float, or str, optional) – Critical stream power to erode rock [E/(TL^2)]

  • discharge_field (float, array of float, or str, optional) – 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.

  • erode_flooded_nodes (bool, optional) – Whether erosion occurs in flooded nodes identified by a depression/lake mapper (e.g., DepressionFinderAndRouter). When set to false, the field flood_status_code must be present on the grid (this is created by the DepressionFinderAndRouter). Default True.

property drainage_area_exp#

Drainage area exponent (units vary).

property fraction_fines#

Fraction of permanently suspendable fines in bedrock [-].

run_one_step(dt)[source]#
run_one_step_basic(dt=10)[source]#
property sediment_influx#

Volumetric sediment influx to each node.

property sediment_porosity#

Sediment porosity [-].

property settling_velocity#

Effective settling velocity for chosen grain size metric [L/T].

property slope_exp#

Slope exponent (units vary).