This module attempts to “componentify” GT’s Fastscape stream power erosion. Created DEJH, March 2014.
FastscapeEroder
(*args, **kwds)[source]¶Bases: landlab.core.model_component.Component
This class uses the BraunWillett 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 numericallycorrect 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 * (rainfall_intensity*A)**m * S**n  threshold_sp,
if K * A**m * S**n > threshold_sp, and:
E = 0,
if K * A**m * S**n <= threshold_sp.
This module assumes you have already run
landlab.components.flow_routing.route_flow_dn.FlowRouter.route_flow()
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
FlowRouter.route_flow).
The primary method of this class is run_one_step()
.
Construction:
FastscapeEroder(grid, K_sp=None, m_sp=0.5, n_sp=1., threshold_sp=0.,
rainfall_intensity=1.)
Parameters:  grid : ModelGrid
K_sp : float, array, or field name
m_sp : float, optional
n_sp : float, optional, ~ 0.5<n_sp<4.
threshold_sp : float, array, or field name
rainfall_intensity : float; optional


Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab import CLOSED_BOUNDARY, FIXED_VALUE_BOUNDARY
>>> from landlab.components import FlowRouter
>>> mg = RasterModelGrid((5, 5), 10.)
>>> z = np.array([7., 7., 7., 7., 7.,
... 7., 5., 3.2, 6., 7.,
... 7., 2., 3., 5., 7.,
... 7., 1., 1.9, 4., 7.,
... 7., 0., 7., 7., 7.])
>>> z = mg.add_field('node', 'topographic__elevation', z)
>>> fr = FlowRouter(mg)
>>> sp = FastscapeEroder(mg, K_sp=1.)
>>> fr.run_one_step()
>>> sp.run_one_step(dt=1.)
>>> 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), 1.)
>>> z = np.array(mg2.node_x**2.)
>>> z = mg2.add_field('node', 'topographic__elevation', z)
>>> mg2.status_at_node[mg2.nodes_at_left_edge] = FIXED_VALUE_BOUNDARY
>>> mg2.status_at_node[mg2.nodes_at_top_edge] = CLOSED_BOUNDARY
>>> mg2.status_at_node[mg2.nodes_at_bottom_edge] = CLOSED_BOUNDARY
>>> mg2.status_at_node[mg2.nodes_at_right_edge] = CLOSED_BOUNDARY
>>> fr2 = FlowRouter(mg2)
>>> sp2 = FastscapeEroder(mg2, K_sp=0.1, m_sp=0., n_sp=2.,
... threshold_sp=2.)
>>> fr2.run_one_step()
>>> sp2.run_one_step(dt=10.)
>>> z.reshape((3, 7))[1, :]
array([ 0. , 1. , 4. , 8.52493781,
13.29039716, 18.44367965, 36. ])
>>> mg3 = RasterModelGrid((3, 7), 1.)
>>> z = np.array(mg3.node_x**2.)
>>> z = mg3.add_field('node', 'topographic__elevation', z)
>>> mg3.status_at_node[mg3.nodes_at_left_edge] = FIXED_VALUE_BOUNDARY
>>> mg3.status_at_node[mg3.nodes_at_top_edge] = CLOSED_BOUNDARY
>>> mg3.status_at_node[mg3.nodes_at_bottom_edge] = CLOSED_BOUNDARY
>>> mg3.status_at_node[mg3.nodes_at_right_edge] = CLOSED_BOUNDARY
>>> fr3 = FlowRouter(mg3)
>>> K_field = mg3.ones('node') # K can be a field
>>> sp3 = FastscapeEroder(mg3, K_sp=K_field, m_sp=1., n_sp=0.6,
... threshold_sp=mg3.node_x,
... rainfall_intensity=2.)
>>> fr3.run_one_step()
>>> sp3.run_one_step(1.)
>>> z.reshape((3, 7))[1, :]
array([ 0. , 0.0647484 , 0.58634455, 2.67253503,
8.49212152, 20.92606987, 36. ])
>>> previous_z = z.copy()
>>> sp3.run_one_step(1., rainfall_intensity_if_used=0.)
>>> np.allclose(z, previous_z)
True
erode
(grid_in, dt=None, K_if_used=None, flooded_nodes=None, rainfall_intensity_if_used=None)[source]¶This method implements the stream power erosion, following the Braun Willett (2013) implicit Fastscape algorithm. This should allow it to be stable against larger timesteps than an explicit stream power scheme.
This driving method for this component is now superceded by the new,
standardized wrapper run_one_step()
, but is retained for
back compatibility.
Set ‘K_if_used’ as a field name or nnodeslong array if you set K_sp as ‘array’ during initialization.
It returns the grid, in which it will have modified the value of value_field, as specified in component initialization.
Parameters:  grid_in : a grid
dt : float
K_if_used : array (optional)
flooded_nodes : ndarray of int (optional)
rainfall_intensity_if_used : float or None (optional)


Returns:  grid

run_one_step
(dt, flooded_nodes=None, rainfall_intensity_if_used=None, **kwds)[source]¶This method implements the stream power erosion across one time interval, dt, following the BraunWillett (2013) implicit Fastscape algorithm.
This follows Landlab standardized component design, and supercedes the
old driving method erode()
.
Parameters:  dt : float
flooded_nodes : ndarray of int (optional)
rainfall_intensity_if_used : float or None (optional)


StreamPowerEroder
(*args, **kwds)[source]¶Bases: landlab.core.model_component.Component
Erode where channels are.
Implemented as:
and if \(E < 0\), \(E = 0\).
If use_W
is declared and True
, the module instead implements:
DEJH Sept 2013, major modifications Sept 14 and May 16. This component now wraps Fastscapestyle functionality under the hood.
Note that although the BraunWillett (2013) scheme that underlies this component is nominally implicit, and will reach a numericallycorrect 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.
NB: If you want spatially or temporally variable runoff, pass the runoff values at each pixel to the flow router using the input argument use_Q.
Construction:
StreamPowerEroder(grid, K_sp=None, threshold_sp=0., sp_type='set_mn',
m_sp=0.5, n_sp=1., a_sp=None, b_sp=None, c_sp=None,
use_W=None, use_Q=None)
Parameters:  grid : ModelGrid
K_sp : float, array, or field name
threshold_sp : positive float, optional
sp_type : {‘set_mn’, ‘Total’, ‘Unit’, ‘Shear_stress’}
m_sp : float, optional
n_sp : float, optional, ~ 0.5<n_sp<4.
a_sp : float, optional
b_sp : float, optional
c_sp : float, optional
use_W : None, array, or field name, optional
use_Q : None, array, or field name, optional


Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab import CLOSED_BOUNDARY, FIXED_VALUE_BOUNDARY
>>> from landlab.components import FlowRouter
>>> from landlab.components import StreamPowerEroder
>>> mg = RasterModelGrid((5, 5), 10.)
>>> z = np.array([7., 7., 7., 7., 7.,
... 7., 5., 3.2, 6., 7.,
... 7., 2., 3., 5., 7.,
... 7., 1., 1.9, 4., 7.,
... 7., 0., 7., 7., 7.])
>>> z = mg.add_field('node', 'topographic__elevation', z)
>>> fr = FlowRouter(mg)
>>> sp = StreamPowerEroder(mg, K_sp=1.)
>>> fr.run_one_step()
>>> sp.run_one_step(dt=1.)
>>> 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), 1.)
>>> z = np.array(mg2.node_x**2.)
>>> z = mg2.add_field('node', 'topographic__elevation', z)
>>> mg2.status_at_node[mg2.nodes_at_left_edge] = FIXED_VALUE_BOUNDARY
>>> mg2.status_at_node[mg2.nodes_at_top_edge] = CLOSED_BOUNDARY
>>> mg2.status_at_node[mg2.nodes_at_bottom_edge] = CLOSED_BOUNDARY
>>> mg2.status_at_node[mg2.nodes_at_right_edge] = CLOSED_BOUNDARY
>>> fr2 = FlowRouter(mg2)
>>> sp2 = StreamPowerEroder(mg2, K_sp=0.1, m_sp=0., n_sp=2.,
... threshold_sp=2.)
>>> fr2.run_one_step()
>>> sp2.run_one_step(dt=10.)
>>> z.reshape((3, 7))[1, :]
array([ 0. , 1. , 4. , 8.52493781,
13.29039716, 18.44367965, 36. ])
>>> mg3 = RasterModelGrid((5, 5), 2.)
>>> z = mg.node_x/100.
>>> z = mg3.add_field('node', 'topographic__elevation', z)
>>> mg3.status_at_node[mg3.nodes_at_left_edge] = FIXED_VALUE_BOUNDARY
>>> mg3.status_at_node[mg3.nodes_at_top_edge] = CLOSED_BOUNDARY
>>> mg3.status_at_node[mg3.nodes_at_bottom_edge] = CLOSED_BOUNDARY
>>> mg3.status_at_node[mg3.nodes_at_right_edge] = CLOSED_BOUNDARY
>>> mg3.at_node['water__unit_flux_in'] = mg3.node_y
>>> fr3 = FlowRouter(mg3)
>>> Q = mg3.at_node['surface_water__discharge']
>>> sp3 = StreamPowerEroder(mg3, K_sp=1., sp_type='Unit', a_sp=1.,
... b_sp=0.5, c_sp=1., use_Q=Q)
>>> fr3.run_one_step()
>>> sp3.run_one_step(1.)
>>> 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 ])
erode
(grid, dt, node_elevs='topographic__elevation', node_drainage_areas='drainage_area', flow_receiver='flow__receiver_node', node_order_upstream='flow__upstream_node_order', slopes_at_nodes='topographic__steepest_slope', link_node_mapping='flow__link_to_receiver_node', link_slopes=None, slopes_from_elevs=None, W_if_used=None, Q_if_used=None, K_if_used=None, flooded_nodes=None)[source]¶Note
deprecated
This run method is now DEPRECATED. Use the fully standardized
method run_one_step()
instead.
A simple, explicit implementation of a stream power algorithm.
Parameters:  grid : RasterModelGrid
dt : float
node_elevs : str or ndarray, optional
node_drainage_areas: str or ndarray, optional
flow_receiver, node_order_upstream : str or ndarray, optional
slopes_from_elevs : str, optional
W_if_used, Q_if_used : str or ndarray, optional


Returns:  tuple

run_one_step
(dt, flooded_nodes=None, **kwds)[source]¶A simple, explicit implementation of a stream power algorithm.
This component now looks exclusively for the field ‘topographic__steepest_slope’ at each node to determine the local slope (previoiusly it was possible to map values from links explicitly within the component, but this functionality is now deprecated).
If you are routing across flooded depressions in your flow routing scheme, be sure to set flooded_nodes with a boolean array or array of IDs 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 flooded_nodes is set.
Parameters:  dt : float
flooded_nodes : ndarray of int (optional)


SedDepEroder
(grid, K_sp=1e06, g=9.81, 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.5e07, 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, **kwds)[source]¶Bases: landlab.core.model_component.Component
This module implements sediment flux dependent channel incision following:
E = f(Qs, Qc) * ([a stream powerlike 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 MeyerPeter Muller approach is used, or whether simpler streampowerlike 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 convexup 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_densityfluid_density)/fluid_density *
g * D_char**3) * (shields_stress  threshold_shields)**1.5
shields_stress = shear_stress / (g * (sed_densityfluid_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.
Construction:
SedDepEroder(grid, K_sp=1.e6, g=9.81, rock_density=2700,
sediment_density=2700, fluid_density=1000,
runoff_rate=1.,
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.,
K_t=1.e4, m_t=1.5, n_t=1., C_MPM=1., a_sp=1.,
b_sp=0.5, c_sp=1., k_w=2.5, k_Q=2.5e7,
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)
Parameters:  grid : a ModelGrid
K_sp : float (time unit must be years)
g : float (m/s**2)
rock_density : float (Kg m**3)
sediment_density : float (Kg m**3)
fluid_density : float (Kg m**3)
runoff_rate : float, array or field name (m/s)
pseudoimplicit_repeats : int
return_stream_properties : bool
sed_dependency_type : {‘generalized_humped’, ‘None’, ‘linear_decline’,
Qc : {‘power_law’, ‘MPM’}
If ``sed_dependency_type == ‘generalized_humped’``... kappa_hump : float
nu_hump : float
phi_hump : float
c_hump : float
If ``Qc == ‘power_law’``... m_sp : float
n_sp : float
K_t : float (time unit must be in years)
m_t : float
n_t : float
if ``Qc == ‘MPM’``... C_MPM : float
a_sp : float
b_sp : float
c_sp : float
k_w : float (unit variable with b_sp)
k_Q : float (unit variable with c_sp, but time unit in seconds)
mannings_n : float
threshold_shear_stress : None or float (Pa)
Dchar :None, float, array, or field name (m)
set_threshold_from_Dchar : bool
set_Dchar_from_threshold : bool
threshold_Shields : None or float
slope_sensitive_threshold : bool


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 : float or array 

Examples
>>> from landlab import RasterModelGrid
>>> mg1 = RasterModelGrid((3,4), 1.)
>>> thresh_shields = np.arange(1, mg1.number_of_nodes+1, dtype=float)
>>> thresh_shields /= 100.
>>> sde1 = SedDepEroder(mg1, threshold_shear_stress=100., Qc='MPM',
... Dchar=None, set_threshold_from_Dchar=False,
... set_Dchar_from_threshold=True,
... threshold_Shields=thresh_shields)
>>> sde1.characteristic_grainsize
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), 1.)
>>> sde2 = SedDepEroder(mg2, threshold_shear_stress=100., Qc='MPM',
... Dchar=None, set_threshold_from_Dchar=False,
... set_Dchar_from_threshold=True,
... threshold_Shields=None,
... slope_sensitive_threshold=True)
>>> S = mg2.add_ones('node', 'topographic__steepest_slope')
>>> S *= 0.05 # thresh = 100 Pa @ 5pc slope
>>> sde2.characteristic_grainsize
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])
erode
(dt, flooded_depths=None, **kwds)[source]¶Erode and deposit on the channel bed for a duration of dt.
Erosion occurs according to the sediment dependent rules specified during initialization.
Parameters:  dt : float (years, only!)
flooded_depths : array or field name (m)


get_sed_flux_function_pseudoimplicit
(sed_in, trans_cap_vol_out, prefactor_for_volume, prefactor_for_dz)[source]¶run_one_step
(dt, flooded_depths=None, **kwds)[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!)
flooded_depths : array or field name (m)

