FastscapeEroder: Compute fluvial erosion using stream power theory (“fastscape” algorithm)#

Fastscape stream power erosion.

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

Bases: Component

Fastscape stream power erosion.

This class uses the Braun-Willett Fastscape approach to calculate the amount of erosion at each node in a grid, following a stream power framework. This should allow it to be stable against larger timesteps than an explicit stream power scheme.

Note that although this scheme is nominally implicit, and will reach a numerically-correct solution under topographic steady state regardless of timestep length, the accuracy of transient solutions is not timestep independent (see Braun & Willett 2013, Appendix B for further details). Although the scheme remains significantly more robust and permits longer timesteps than a traditional explicit solver under such conditions, it is still possible to create numerical instability through use of too long a timestep while using this component. The user is cautioned to check their implementation is behaving stably before fully trusting it.

Stream power erosion is implemented as:

\[E = K A ^ m S ^ n - \textit{threshold_sp}\]

if \(K A ^ m S ^ n > \textit{threshold_sp}\), and:

\[E = 0,\]

if \(K A^m S^n <= \textit{threshold_sp}\).

This module assumes you have already run landlab.components.flow_accum.flow_accumulator.FlowAccumulator.run_one_step in the same timestep. It looks for ‘flow__upstream_node_order’, ‘flow__link_to_receiver_node’, ‘drainage_area’, ‘flow__receiver_node’, and ‘topographic__elevation’ at the nodes in the grid. ‘drainage_area’ should be in area upstream, not volume (i.e., set runoff_rate=1.0 when calling FlowAccumulator.run_one_step).

The primary method of this class is run_one_step.

Examples

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

>>> grid = RasterModelGrid((5, 5), xy_spacing=10.0)
>>> z = [
...     [7.0, 7.0, 7.0, 7.0, 7.0],
...     [7.0, 5.0, 3.2, 6.0, 7.0],
...     [7.0, 2.0, 3.0, 5.0, 7.0],
...     [7.0, 1.0, 1.9, 4.0, 7.0],
...     [7.0, 0.0, 7.0, 7.0, 7.0],
... ]
>>> z = grid.add_field("topographic__elevation", z, at="node")
>>> fr = FlowAccumulator(grid, flow_director="D8")
>>> sp = FastscapeEroder(grid, K_sp=1.0)
>>> fr.run_one_step()
>>> sp.run_one_step(dt=1.0)
>>> z
array([ 7.        ,  7.        ,  7.        ,  7.        ,  7.        ,
        7.        ,  2.92996598,  2.02996598,  4.01498299,  7.        ,
        7.        ,  0.85993197,  1.87743897,  3.28268321,  7.        ,
        7.        ,  0.28989795,  0.85403051,  2.42701526,  7.        ,
        7.        ,  0.        ,  7.        ,  7.        ,  7.        ])

>>> grid = RasterModelGrid((3, 7), xy_spacing=1.0)
>>> z = np.array(grid.node_x**2.0)
>>> z = grid.add_field("topographic__elevation", z, at="node")
>>> grid.status_at_node[grid.nodes_at_left_edge] = grid.BC_NODE_IS_FIXED_VALUE
>>> grid.status_at_node[grid.nodes_at_top_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_bottom_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_right_edge] = grid.BC_NODE_IS_CLOSED
>>> fr = FlowAccumulator(grid, flow_director="D8")
>>> sp = FastscapeEroder(grid, K_sp=0.1, m_sp=0.0, n_sp=2.0, threshold_sp=2.0)
>>> fr.run_one_step()
>>> sp.run_one_step(dt=10.0)
>>> z.reshape(grid.shape)[1, :]
array([  0.        ,   1.        ,   4.        ,   8.52493781,
        13.29039716,  18.44367965,  36.        ])

>>> grid = RasterModelGrid((3, 7), xy_spacing=1.0)
>>> z = np.array(grid.node_x**2.0)
>>> z = grid.add_field("topographic__elevation", z, at="node")
>>> grid.status_at_node[grid.nodes_at_left_edge] = grid.BC_NODE_IS_FIXED_VALUE
>>> grid.status_at_node[grid.nodes_at_top_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_bottom_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_right_edge] = grid.BC_NODE_IS_CLOSED
>>> cell_area = 1.0
>>> fr = FlowAccumulator(grid, flow_director="D8", runoff_rate=2.0)
>>> grid.at_node["water__unit_flux_in"]
array([ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,
        2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.])
>>> K_field = grid.ones(at="node")  # K can be a field
>>> sp = FastscapeEroder(
...     grid,
...     K_sp=K_field,
...     m_sp=1.0,
...     n_sp=0.6,
...     threshold_sp=grid.node_x,
...     discharge_field="surface_water__discharge",
... )
>>> fr.run_one_step()
>>> sp.run_one_step(1.0)
>>> z.reshape(grid.shape)[1, :]
array([  0.        ,   0.0647484 ,   0.58634455,   2.67253503,
         8.49212152,  20.92606987,  36.        ])

References

Required Software Citation(s) Specific to this Component

None Listed

Additional References

Braun, J., Willett, S. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. Geomorphology 180-181(C), 170-179. https://dx.doi.org/10.1016/j.geomorph.2012.10.008

Initialize the Fastscape stream power component. Note: a timestep, dt, can no longer be supplied to this component through the input file. It must instead be passed directly to the run method.

Parameters:

  • grid (ModelGrid) – A grid.

  • K_sp (float, array, or field name) – K in the stream power equation (units vary with other parameters).

  • m_sp (float, optional) – m in the stream power equation (power on drainage area).

  • n_sp (float, optional) – n in the stream power equation (power on slope).

  • threshold_sp (float, array, or field name) – Erosion threshold in the stream power equation.

  • discharge_field (float, field name, or array, optional) – Discharge [L^2/T]. The default is to use the grid field ‘drainage_area’. To use custom spatially/temporally varying rainfall, use ‘water__unit_flux_in’ to specify water input to the FlowAccumulator and use “surface_water__discharge” for this keyword argument.

  • 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 K#

Erodibility (units depend on m_sp).

__init__(grid, K_sp=0.001, m_sp=0.5, n_sp=1.0, threshold_sp=0.0, discharge_field='drainage_area', erode_flooded_nodes=True)[source]#

Initialize the Fastscape stream power component. Note: a timestep, dt, can no longer be supplied to this component through the input file. It must instead be passed directly to the run method.

Parameters:

  • grid (ModelGrid) – A grid.

  • K_sp (float, array, or field name) – K in the stream power equation (units vary with other parameters).

  • m_sp (float, optional) – m in the stream power equation (power on drainage area).

  • n_sp (float, optional) – n in the stream power equation (power on slope).

  • threshold_sp (float, array, or field name) – Erosion threshold in the stream power equation.

  • discharge_field (float, field name, or array, optional) – Discharge [L^2/T]. The default is to use the grid field ‘drainage_area’. To use custom spatially/temporally varying rainfall, use ‘water__unit_flux_in’ to specify water input to the FlowAccumulator and use “surface_water__discharge” for this keyword argument.

  • 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.

run_one_step(dt)[source]#

Erode for a single time step.

This method implements the stream power erosion across one time interval, dt, following the Braun-Willett (2013) implicit Fastscape algorithm.

This follows Landlab standardized component design, and supercedes the old driving method erode.

Parameters:

dt (float) – Time-step size

StreamPower: Compute fluvial erosion using stream power theory (also uses “fastscape” algorithm but provides slightly different formulation and options)#

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

Bases: Component

Erode where channels are.

Implemented as:

\[E = K A^m S^n - sp_{crit},\]

and if \(E < 0\), \(E = 0\).

If channel_width_field is declared and True, the module instead implements:

\[E = K A^m S^n / W - sp_{crit}\]

Note that although the Braun-Willett (2013) scheme that underlies this component is nominally implicit, and will reach a numerically-correct solution under topographic steady state regardless of timestep length, the accuracy of transient solutions is not timestep independent (see Braun & Willett 2013, Appendix B for further details). Although the scheme remains significantly more robust and permits longer timesteps than a traditional explicit solver under such conditions, it is still possible to create numerical instability through use of too long a timestep while using this component. The user is cautioned to check their implementation is behaving stably before fully trusting it.

Examples

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

>>> mg = RasterModelGrid((5, 5), xy_spacing=10.0)
>>> z = np.array(
...     [
...         [7.0, 7.0, 7.0, 7.0, 7.0],
...         [7.0, 5.0, 3.2, 6.0, 7.0],
...         [7.0, 2.0, 3.0, 5.0, 7.0],
...         [7.0, 1.0, 1.9, 4.0, 7.0],
...         [7.0, 0.0, 7.0, 7.0, 7.0],
...     ]
... )
>>> z = mg.add_field("topographic__elevation", z, at="node")
>>> fr = FlowAccumulator(mg, flow_director="D8")
>>> sp = StreamPowerEroder(mg, K_sp=1.0)
>>> fr.run_one_step()
>>> sp.run_one_step(dt=1.0)
>>> z
array([ 7.        ,  7.        ,  7.        ,  7.        ,  7.        ,
        7.        ,  2.92996598,  2.02996598,  4.01498299,  7.        ,
        7.        ,  0.85993197,  1.87743897,  3.28268321,  7.        ,
        7.        ,  0.28989795,  0.85403051,  2.42701526,  7.        ,
        7.        ,  0.        ,  7.        ,  7.        ,  7.        ])

>>> mg2 = RasterModelGrid((3, 7))
>>> z = np.array(mg2.node_x**2.0)
>>> z = mg2.add_field("topographic__elevation", z, at="node")
>>> mg2.status_at_node[mg2.nodes_at_left_edge] = mg2.BC_NODE_IS_FIXED_VALUE
>>> mg2.status_at_node[mg2.nodes_at_top_edge] = mg2.BC_NODE_IS_CLOSED
>>> mg2.status_at_node[mg2.nodes_at_bottom_edge] = mg2.BC_NODE_IS_CLOSED
>>> mg2.status_at_node[mg2.nodes_at_right_edge] = mg2.BC_NODE_IS_CLOSED
>>> fr2 = FlowAccumulator(mg2, flow_director="D8")
>>> sp2 = StreamPowerEroder(mg2, K_sp=0.1, m_sp=0.0, n_sp=2.0, threshold_sp=2.0)
>>> fr2.run_one_step()
>>> sp2.run_one_step(dt=10.0)
>>> z.reshape((3, 7))[1, :]
array([  0.        ,   1.        ,   4.        ,   8.52493781,
        13.29039716,  18.44367965,  36.        ])

>>> mg3 = RasterModelGrid((5, 5), xy_spacing=2.0)
>>> z = mg.node_x / 100.0
>>> z = mg3.add_field("topographic__elevation", z, at="node")
>>> mg3.status_at_node[mg3.nodes_at_left_edge] = mg3.BC_NODE_IS_FIXED_VALUE
>>> mg3.status_at_node[mg3.nodes_at_top_edge] = mg3.BC_NODE_IS_CLOSED
>>> mg3.status_at_node[mg3.nodes_at_bottom_edge] = mg3.BC_NODE_IS_CLOSED
>>> mg3.status_at_node[mg3.nodes_at_right_edge] = mg3.BC_NODE_IS_CLOSED
>>> mg3.at_node["water__unit_flux_in"] = mg3.node_y
>>> fr3 = FlowAccumulator(mg3, flow_director="D8")
>>> sp3 = StreamPowerEroder(
...     mg3,
...     K_sp=1.0,
...     sp_type="Unit",
...     a_sp=1.0,
...     b_sp=0.5,
...     c_sp=1.0,
...     discharge_field="surface_water__discharge",
... )
>>> fr3.run_one_step()
>>> sp3.run_one_step(1.0)
>>> z
array([ 0.        ,  0.1       ,  0.2       ,  0.3       ,  0.4       ,
        0.        ,  0.02898979,  0.0859932 ,  0.17463772,  0.4       ,
        0.        ,  0.02240092,  0.06879049,  0.14586033,  0.4       ,
        0.        ,  0.01907436,  0.05960337,  0.12929386,  0.4       ,
        0.        ,  0.1       ,  0.2       ,  0.3       ,  0.4       ])

References

Required Software Citation(s) Specific to this Component

None Listed

Additional References

Braun, J., Willett, S. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. Geomorphology 180-181(C), 170-179. https://dx.doi.org/10.1016/j.geomorph.2012.10.008

Initialize the StreamPowerEroder.

Parameters:

  • grid (ModelGrid) – A grid.

  • K_sp (float, array, or field name) – K in the stream power equation (units vary with other parameters).

  • threshold_sp (positive float, array, or field name, optional) – The threshold stream power, below which no erosion occurs. This threshold is assumed to be in “stream power” units, i.e., if sp_type is ‘Shear_stress’, the value should be tau**a.

  • sp_type ({'set_mn', 'Total', 'Unit', 'Shear_stress'}) –

    Controls how the law is implemented. If ‘set_mn’, use the supplied values of m_sp and n_sp. Else, component will derive values of m and n from supplied values of a_sp, b_sp, and c_sp, following Whipple and Tucker:

    • If 'Total', m = a * c, n = a.

    • If 'Unit', m = a * c *(1 - b), n = a.

    • If 'Shear_stress', m = 2 * a * c * (1 - b) / 3, n = 2 * a / 3.

  • m_sp (float, optional) – m in the stream power equation (power on drainage area). Overridden if a_sp, b_sp, and c_sp are supplied.

  • n_sp (float, optional, ) – n in the stream power equation (power on slope). Overridden if a_sp, b_sp, and c_sp are supplied.

  • a_sp (float, optional) – The power on the SP/shear term to get the erosion rate; the “erosional process” term. Only used if sp_type is not ‘set_mn’.

  • b_sp (float, optional) – The power on discharge to get width; the “hydraulic geometry” term. Only used if sp_type in (‘Unit’, ‘Shear_stress’).

  • c_sp (float, optional) – The power on area to get discharge; the “basin hydology” term. Only used if sp_type is not ‘set_mn’.

  • channel_width_field (None, float, array, or field name, optional) – If not None, component will look for node-centered data describing channel width or if an array, will take the array as the channel widths. It will use the widths to implement incision ~ stream power per unit width. If sp_type is ‘set_mn’, follows the equation given above. If sp_type in (‘Unit’, ‘Shear_stress’), the width value will be implemented directly. W has no effect if sp_type is ‘Total’.

  • discharge_field (float, field name, or array, optional) – Discharge [L^2/T]. The default is to use the grid field ‘drainage_area’. To use custom spatially/temporally varying rainfall, use ‘water__unit_flux_in’ to specify water input to the FlowAccumulator and use “surface_water__discharge” for this keyword argument.

  • 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 K#

Erodibility (units depend on m_sp).

__init__(grid, K_sp=0.001, threshold_sp=0.0, sp_type='set_mn', m_sp=0.5, n_sp=1.0, a_sp=None, b_sp=None, c_sp=None, channel_width_field=1.0, discharge_field='drainage_area', erode_flooded_nodes=True)[source]#

Initialize the StreamPowerEroder.

Parameters:

  • grid (ModelGrid) – A grid.

  • K_sp (float, array, or field name) – K in the stream power equation (units vary with other parameters).

  • threshold_sp (positive float, array, or field name, optional) – The threshold stream power, below which no erosion occurs. This threshold is assumed to be in “stream power” units, i.e., if sp_type is ‘Shear_stress’, the value should be tau**a.

  • sp_type ({'set_mn', 'Total', 'Unit', 'Shear_stress'}) –

    Controls how the law is implemented. If ‘set_mn’, use the supplied values of m_sp and n_sp. Else, component will derive values of m and n from supplied values of a_sp, b_sp, and c_sp, following Whipple and Tucker:

    • If 'Total', m = a * c, n = a.

    • If 'Unit', m = a * c *(1 - b), n = a.

    • If 'Shear_stress', m = 2 * a * c * (1 - b) / 3, n = 2 * a / 3.

  • m_sp (float, optional) – m in the stream power equation (power on drainage area). Overridden if a_sp, b_sp, and c_sp are supplied.

  • n_sp (float, optional, ) – n in the stream power equation (power on slope). Overridden if a_sp, b_sp, and c_sp are supplied.

  • a_sp (float, optional) – The power on the SP/shear term to get the erosion rate; the “erosional process” term. Only used if sp_type is not ‘set_mn’.

  • b_sp (float, optional) – The power on discharge to get width; the “hydraulic geometry” term. Only used if sp_type in (‘Unit’, ‘Shear_stress’).

  • c_sp (float, optional) – The power on area to get discharge; the “basin hydology” term. Only used if sp_type is not ‘set_mn’.

  • channel_width_field (None, float, array, or field name, optional) – If not None, component will look for node-centered data describing channel width or if an array, will take the array as the channel widths. It will use the widths to implement incision ~ stream power per unit width. If sp_type is ‘set_mn’, follows the equation given above. If sp_type in (‘Unit’, ‘Shear_stress’), the width value will be implemented directly. W has no effect if sp_type is ‘Total’.

  • discharge_field (float, field name, or array, optional) – Discharge [L^2/T]. The default is to use the grid field ‘drainage_area’. To use custom spatially/temporally varying rainfall, use ‘water__unit_flux_in’ to specify water input to the FlowAccumulator and use “surface_water__discharge” for this keyword argument.

  • 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.

run_one_step(dt)[source]#

A simple, explicit implementation of a stream power algorithm.

If you are routing across flooded depressions in your flow routing scheme, be sure to set erode_flooded_nodes flag in the instantiation of the component to ensure erosion cannot occur in the lake. Erosion is always zero if the gradient is adverse, but can still procede as usual on the entry into the depression unless this flag is set.

Parameters:

dt (float) – Time-step size

SedDepEroder: Compute fluvial erosion using using “tools and cover” theory#

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

Bases: Component

This module implements sediment flux dependent channel incision following:

E = f(Qs, Qc) * ([a stream power-like term] - [an optional threshold]),

where E is the bed erosion rate, Qs is the volumetric sediment flux into a node, and Qc is the volumetric sediment transport capacity at that node.

This component is under active research and development; proceed with its use at your own risk.

The details of the implementation are a function of the two key arguments, sed_dependency_type and Qc. The former controls the shape of the sediment dependent response function f(Qs, Qc), the latter controls the way in which sediment transport capacities are calculated (primarily, whether a full Meyer-Peter Muller approach is used, or whether simpler stream-power-like equations can be assumed). For Qc, ‘power_law’ broadly follows the assumptions in Gasparini et al. 2006, 2007; ‘MPM’ broadly follows those in Hobley et al., 2011. Note that a convex-up channel can result in many cases assuming MPM, unless parameters b and c are carefully tuned.

If Qc == 'power_law':

E  = K_sp * f(Qs, Qc) * A ** m_sp * S ** n_sp;
Qc = K_t * A ** m_t * S ** n_t

If Qc == 'MPM':

shear_stress = fluid_density * g * depth * S
             = fluid_density * g * (mannings_n/k_w) ** 0.6 * (
               k_Q* A ** c_sp) ** (0.6 * (1. - b_sp)) * S ** 0.7,
               for consistency with MPM

E = K_sp * f(Qs, Qc) * (shear_stress ** a_sp - [threshold_sp])

Qc = 8 * C_MPM * sqrt((sed_density-fluid_density)/fluid_density *
     g * D_char**3) * (shields_stress - threshold_shields)**1.5

shields_stress = shear_stress / (g * (sed_density-fluid_density) *
                 D_char)

If you choose Qc=’MPM’, you may provide thresholds for both channel incision and shields number, or alternatively set either or both of these threshold dynamically. The minimum shear stress can be made equivalent to the Shields number using set_threshold_from_Dchar, for full consistency with the MPM approach (i.e., the threshold becomes a function of the characteristic grain size on the bed). The Shields threshold itself can also be a weak function of slope if slope_sensitive_threshold, following Lamb et al. 2008, taustar_c = 0.15 * S ** 0.25.

The component is able to handle flooded nodes, if created by a lake filler. It assumes the flow paths found in the fields already reflect any lake routing operations, and then requires the optional argument flooded_depths be passed to the run method. A flooded depression acts as a perfect sediment trap, and will be filled sequentially from the inflow points towards the outflow points.

References

Required Software Citation(s) Specific to this Component

None Listed

Additional References

Hobley, D. E. J., Sinclair, H. D., Mudd, S. M., and Cowie, P. A.: Field calibration of sediment flux dependent river incision, J. Geophys. Res., 116, F04017, doi:10.1029/2010JF001935, 2011.

Constructor for the class.

Parameters:

  • grid (a ModelGrid) – A grid.

  • K_sp (float (time unit must be years)) – K in the stream power equation; the prefactor on the erosion equation (units vary with other parameters).

  • g (float (m/s**2)) – Acceleration due to gravity.

  • rock_density (float (Kg m**-3)) – Bulk intact rock density.

  • sediment_density (float (Kg m**-3)) – Typical density of loose sediment on the bed.

  • fluid_density (float (Kg m**-3)) – Density of the fluid.

  • runoff_rate (float, array or field name (m/s)) – The rate of excess overland flow production at each node (i.e., rainfall rate less infiltration).

  • pseudoimplicit_repeats (int) – Number of loops to perform with the pseudoimplicit iterator, seeking a stable solution. Convergence is typically rapid.

  • return_stream_properties (bool) – Whether to perform a few additional calculations in order to set the additional optional output fields, ‘channel__width’, ‘channel__depth’, and ‘channel__discharge’ (default False).

  • sed_dependency_type ({'generalized_humped', 'None', 'linear_decline',) – ‘almost_parabolic’} The shape of the sediment flux function. For definitions, see Hobley et al., 2011. ‘None’ gives a constant value of 1. NB: ‘parabolic’ is currently not supported, due to numerical stability issues at channel heads.

  • Qc ({'power_law', 'MPM'}) – Whether to use simple stream-power-like equations for both sediment transport capacity and erosion rate, or more complex forms based directly on the Meyer-Peter Muller equation and a shear stress based erosion model consistent with MPM (per Hobley et al., 2011).

  • 'generalized_humped'... (If sed_dependency_type ==) –

  • kappa_hump (float) – Shape parameter for sediment flux function. Primarily controls function amplitude (i.e., scales the function to a maximum of 1). Default follows Leh valley values from Hobley et al., 2011.

  • nu_hump (float) – Shape parameter for sediment flux function. Primarily controls rate of rise of the “tools” limb. Default follows Leh valley values from Hobley et al., 2011.

  • phi_hump (float) – Shape parameter for sediment flux function. Primarily controls rate of fall of the “cover” limb. Default follows Leh valley values from Hobley et al., 2011.

  • c_hump (float) – Shape parameter for sediment flux function. Primarily controls degree of function asymmetry. Default follows Leh valley values from Hobley et al., 2011.

  • 'power_law'... (If Qc ==) –

  • m_sp (float) – Power on drainage area in the erosion equation.

  • n_sp (float) – Power on slope in the erosion equation.

  • K_t (float (time unit must be in years)) – Prefactor in the transport capacity equation.

  • m_t (float) – Power on drainage area in the transport capacity equation.

  • n_t (float) – Power on slope in the transport capacity equation.

  • 'MPM'... (if Qc ==) –

  • C_MPM (float) – A prefactor on the MPM relation, allowing tuning to known sediment saturation conditions (leave as 1. in most cases).

  • a_sp (float) – Power on shear stress to give erosion rate.

  • b_sp (float) – Power on drainage area to give channel width.

  • c_sp (float) – Power on drainage area to give discharge.

  • k_w (float (unit variable with b_sp)) – Prefactor on A**b_sp to give channel width.

  • k_Q (float (unit variable with c_sp, but time unit in seconds)) – Prefactor on A**c_sp to give discharge.

  • mannings_n (float) – Manning’s n for the channel.

  • threshold_shear_stress (None or float (Pa)) – The threshold shear stress in the equation for erosion rate. If None, implies that set_threshold_from_Dchar is True, and this parameter will get set from the Dchar value and critical Shields number.

  • Dchar (None, float, array, or field name (m)) – The characteristic grain size on the bed, that controls the relationship between critical Shields number and critical shear stress. If None, implies that set_Dchar_from_threshold is True, and this parameter will get set from the threshold_shear_stress value and critical Shields number.

  • set_threshold_from_Dchar (bool) – If True (default), threshold_shear_stress will be set based on Dchar and threshold_Shields.

  • set_Dchar_from_threshold (bool) – If True, Dchar will be set based on threshold_shear_stress and threshold_Shields. Default is False.

  • threshold_Shields (None or float) – The threshold Shields number. If None, implies that slope_sensitive_threshold is True.

  • slope_sensitive_threshold (bool) – If True, the threshold_Shields will be set according to 0.15 * S ** 0.25, per Lamb et al., 2008 & Hobley et al., 2011.

  • flooded_depths (array or field name (m)) – Depths of flooding at each node, zero where no lake. Note that the component will dynamically update this array as it fills nodes with sediment (…but does NOT update any other related lake fields).

__init__(grid, K_sp=1e-06, g=9.80665, rock_density=2700, sediment_density=2700, fluid_density=1000, runoff_rate=1.0, sed_dependency_type='generalized_humped', kappa_hump=13.683, nu_hump=1.13, phi_hump=4.24, c_hump=0.00181, Qc='power_law', m_sp=0.5, n_sp=1.0, K_t=0.0001, m_t=1.5, n_t=1.0, C_MPM=1.0, a_sp=1.0, b_sp=0.5, c_sp=1.0, k_w=2.5, k_Q=2.5e-07, mannings_n=0.05, threshold_shear_stress=None, Dchar=0.05, set_threshold_from_Dchar=True, set_Dchar_from_threshold=False, threshold_Shields=0.05, slope_sensitive_threshold=False, pseudoimplicit_repeats=5, return_stream_properties=False, flooded_depths=None)[source]#

Constructor for the class.

Parameters:

  • grid (a ModelGrid) – A grid.

  • K_sp (float (time unit must be years)) – K in the stream power equation; the prefactor on the erosion equation (units vary with other parameters).

  • g (float (m/s**2)) – Acceleration due to gravity.

  • rock_density (float (Kg m**-3)) – Bulk intact rock density.

  • sediment_density (float (Kg m**-3)) – Typical density of loose sediment on the bed.

  • fluid_density (float (Kg m**-3)) – Density of the fluid.

  • runoff_rate (float, array or field name (m/s)) – The rate of excess overland flow production at each node (i.e., rainfall rate less infiltration).

  • pseudoimplicit_repeats (int) – Number of loops to perform with the pseudoimplicit iterator, seeking a stable solution. Convergence is typically rapid.

  • return_stream_properties (bool) – Whether to perform a few additional calculations in order to set the additional optional output fields, ‘channel__width’, ‘channel__depth’, and ‘channel__discharge’ (default False).

  • sed_dependency_type ({'generalized_humped', 'None', 'linear_decline',) – ‘almost_parabolic’} The shape of the sediment flux function. For definitions, see Hobley et al., 2011. ‘None’ gives a constant value of 1. NB: ‘parabolic’ is currently not supported, due to numerical stability issues at channel heads.

  • Qc ({'power_law', 'MPM'}) – Whether to use simple stream-power-like equations for both sediment transport capacity and erosion rate, or more complex forms based directly on the Meyer-Peter Muller equation and a shear stress based erosion model consistent with MPM (per Hobley et al., 2011).

  • 'generalized_humped'... (If sed_dependency_type ==) –

  • kappa_hump (float) – Shape parameter for sediment flux function. Primarily controls function amplitude (i.e., scales the function to a maximum of 1). Default follows Leh valley values from Hobley et al., 2011.

  • nu_hump (float) – Shape parameter for sediment flux function. Primarily controls rate of rise of the “tools” limb. Default follows Leh valley values from Hobley et al., 2011.

  • phi_hump (float) – Shape parameter for sediment flux function. Primarily controls rate of fall of the “cover” limb. Default follows Leh valley values from Hobley et al., 2011.

  • c_hump (float) – Shape parameter for sediment flux function. Primarily controls degree of function asymmetry. Default follows Leh valley values from Hobley et al., 2011.

  • 'power_law'... (If Qc ==) –

  • m_sp (float) – Power on drainage area in the erosion equation.

  • n_sp (float) – Power on slope in the erosion equation.

  • K_t (float (time unit must be in years)) – Prefactor in the transport capacity equation.

  • m_t (float) – Power on drainage area in the transport capacity equation.

  • n_t (float) – Power on slope in the transport capacity equation.

  • 'MPM'... (if Qc ==) –

  • C_MPM (float) – A prefactor on the MPM relation, allowing tuning to known sediment saturation conditions (leave as 1. in most cases).

  • a_sp (float) – Power on shear stress to give erosion rate.

  • b_sp (float) – Power on drainage area to give channel width.

  • c_sp (float) – Power on drainage area to give discharge.

  • k_w (float (unit variable with b_sp)) – Prefactor on A**b_sp to give channel width.

  • k_Q (float (unit variable with c_sp, but time unit in seconds)) – Prefactor on A**c_sp to give discharge.

  • mannings_n (float) – Manning’s n for the channel.

  • threshold_shear_stress (None or float (Pa)) – The threshold shear stress in the equation for erosion rate. If None, implies that set_threshold_from_Dchar is True, and this parameter will get set from the Dchar value and critical Shields number.

  • Dchar (None, float, array, or field name (m)) – The characteristic grain size on the bed, that controls the relationship between critical Shields number and critical shear stress. If None, implies that set_Dchar_from_threshold is True, and this parameter will get set from the threshold_shear_stress value and critical Shields number.

  • set_threshold_from_Dchar (bool) – If True (default), threshold_shear_stress will be set based on Dchar and threshold_Shields.

  • set_Dchar_from_threshold (bool) – If True, Dchar will be set based on threshold_shear_stress and threshold_Shields. Default is False.

  • threshold_Shields (None or float) – The threshold Shields number. If None, implies that slope_sensitive_threshold is True.

  • slope_sensitive_threshold (bool) – If True, the threshold_Shields will be set according to 0.15 * S ** 0.25, per Lamb et al., 2008 & Hobley et al., 2011.

  • flooded_depths (array or field name (m)) – Depths of flooding at each node, zero where no lake. Note that the component will dynamically update this array as it fills nodes with sediment (…but does NOT update any other related lake fields).

property characteristic_grainsize#

Return the characteristic grain size used by the component.

Particularly useful if the set_Dchar_from_threshold flag was True at initialization.

Returns:

Dchar

Return type:

float or array

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, SedDepEroder
>>> mg1 = RasterModelGrid((3, 4))
>>> z1 = mg1.add_zeros("node", "topographic__elevation")
>>> fa1 = FlowAccumulator(mg1)
>>> thresh_shields = np.arange(1, mg1.number_of_nodes + 1, dtype=float)
>>> thresh_shields /= 100.0
>>> sde1 = SedDepEroder(
...     mg1,
...     threshold_shear_stress=100.0,
...     Qc="MPM",
...     Dchar=None,
...     set_threshold_from_Dchar=False,
...     set_Dchar_from_threshold=True,
...     threshold_Shields=thresh_shields,
...     g=9.81,
... )
>>> sde1.characteristic_grainsize.reshape(mg1.shape)
array([[ 0.59962823,  0.29981412,  0.19987608,  0.14990706],
       [ 0.11992565,  0.09993804,  0.08566118,  0.07495353],
       [ 0.06662536,  0.05996282,  0.05451166,  0.04996902]])

>>> mg2 = RasterModelGrid((3, 4))
>>> z2 = mg2.add_zeros("node", "topographic__elevation")
>>> fa2 = FlowAccumulator(mg2)
>>> sde2 = SedDepEroder(
...     mg2,
...     threshold_shear_stress=100.0,
...     Qc="MPM",
...     Dchar=None,
...     set_threshold_from_Dchar=False,
...     set_Dchar_from_threshold=True,
...     threshold_Shields=None,
...     slope_sensitive_threshold=True,
...     g=9.81,
... )
>>> S = mg2.at_node["topographic__steepest_slope"]
>>> S[:] = 0.05  # thresh = 100 Pa @ 5pc slope
>>> sde2.characteristic_grainsize.reshape(mg2.shape)
array([[ 0.08453729,  0.08453729,  0.08453729,  0.08453729],
       [ 0.08453729,  0.08453729,  0.08453729,  0.08453729],
       [ 0.08453729,  0.08453729,  0.08453729,  0.08453729]])
get_sed_flux_function(rel_sed_flux)[source]#

Get the sediment flux function.

Parameters:

rel_sed_flux –

get_sed_flux_function_pseudoimplicit(sed_in, trans_cap_vol_out, prefactor_for_volume, prefactor_for_dz)[source]#

Get the pseudoimplicit sediment flux function.

Parameters:

  • sed_in –

  • trans_cap_vol_out –

  • prefactor_for_volume –

  • prefactor_for_dz –

run_one_step(dt)[source]#

Run the component across one timestep increment, dt.

Erosion occurs according to the sediment dependent rules specified during initialization. Method is fully equivalent to the erode method.

Parameters:

dt (float (years, only!)) – Timestep for which to run the component.

StreamPowerSmoothThresholdEroder: Compute fluvial erosion using stream power theory with a numerically smoothed threshold#

stream_power_smooth_threshold.py: Defines the StreamPowerSmoothThresholdEroder, which is a version of the FastscapeEroder (derived from it).

StreamPowerSmoothThresholdEroder uses a mathematically smooth threshold formulation, rather than one with a singularity. The erosion rate is defined as follows:

$omega = K A^m S$

$E = omega - omega_c left[ 1 - exp ( -omega / omega_c ) right]$

Created on Sat Nov 26 08:36:49 2016

@author: gtucker

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

Bases: FastscapeEroder

Stream erosion component with smooth threshold function.

Parameters:

  • grid (ModelGrid) – A grid.

  • K_sp (float, array, or field name) – K in the stream power equation (units vary with other parameters).

  • m_sp (float, optional) – m in the stream power equation (power on drainage area).

  • n_sp (float, optional, ) – n in the stream power equation (power on slope). NOTE: NOT PRESENTLY HONORED BY StreamPowerSmoothThresholdEroder (TODO)

  • threshold_sp (float (TODO: array, or field name)) – The threshold stream power.

  • discharge_field (float, field name, or array, optional) – Discharge [L^2/T]. The default is to use the grid field ‘drainage_area’. To use custom spatially/temporally varying rainfall, use ‘water__unit_flux_in’ to specify water input to the FlowAccumulator and use “surface_water__discharge” for this keyword argument.

  • 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.

Examples

>>> from landlab import RasterModelGrid
>>> rg = RasterModelGrid((3, 4))
>>> rg.set_closed_boundaries_at_grid_edges(False, True, True, True)
>>> z = rg.add_zeros("node", "topographic__elevation")
>>> z[5] = 2.0
>>> z[6] = 1.0
>>> from landlab.components import FlowAccumulator
>>> fr = FlowAccumulator(rg, flow_director="D4")
>>> fr.run_one_step()
>>> from landlab.components import StreamPowerSmoothThresholdEroder
>>> sp = StreamPowerSmoothThresholdEroder(rg, K_sp=1.0)
>>> sp.thresholds
1.0
>>> sp.run_one_step(dt=1.0)
>>> import numpy as np
>>> np.round(z[5:7], 3)
array([ 1.646,  0.667])
>>> z[5] = 2.0
>>> z[6] = 1.0
>>> import numpy as np
>>> q = np.zeros(rg.number_of_nodes) + 0.25
>>> q[6] = 100.0
>>> sp = StreamPowerSmoothThresholdEroder(rg, K_sp=1.0, discharge_field=q)
>>> sp.run_one_step(dt=1.0)
>>> np.round(z[5:7], 3)
array([ 1.754,  0.164])

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

Braun, J., Willett, S. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. Geomorphology 180-181(C), 170-179. https://dx.doi.org/10.1016/j.geomorph.2012.10.008

Initialize StreamPowerSmoothThresholdEroder.

__init__(grid, K_sp=None, m_sp=0.5, n_sp=1.0, threshold_sp=1.0, discharge_field='drainage_area', erode_flooded_nodes=True)[source]#

Initialize StreamPowerSmoothThresholdEroder.

property alpha#

Erosion term divided by link length.

Alpha is given as:

alpha = K A^m dt / L

where K is the erodibility, A is the drainage area, m is the drainage area exponent, dt is the timestep, and L is the link length.

property delta#

Erosion term divided by link length and erosion threshold.

delta is given as:

delta = K A^m dt / (L * omega_c)

where K is the erodibility, A is the drainage area, m is the drainage area exponent, dt is the timestep, L is the link length, and omega_c is the erosion threshold.

property gamma#

Erosion threshold times timestep.

run_one_step(dt, runoff_rate=None)[source]#

Run one forward iteration of duration dt.

Parameters:

  • dt (float) – Time step size

  • runoff_rate ((not used yet)) – (to be added later)

Examples

>>> from landlab import RasterModelGrid
>>> rg = RasterModelGrid((3, 3))
>>> rg.set_closed_boundaries_at_grid_edges(False, True, True, True)
>>> z = rg.add_zeros("node", "topographic__elevation")
>>> z[4] = 1.0
>>> from landlab.components import FlowAccumulator
>>> fr = FlowAccumulator(rg, flow_director="D4")
>>> fr.run_one_step()
>>> from landlab.components import StreamPowerSmoothThresholdEroder
>>> sp = StreamPowerSmoothThresholdEroder(rg, K_sp=1.0)
>>> sp.run_one_step(dt=1.0)
>>> sp.alpha
array([ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.])
>>> sp.gamma
array([ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.])
>>> sp.delta
array([ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.])
property thresholds#

Erosion thresholds.