import numpy as np
from scipy.integrate import quad
from landlab.components.erosion_deposition.generalized_erosion_deposition import (
DEFAULT_MINIMUM_TIME_STEP,
)
from landlab.components.erosion_deposition.generalized_erosion_deposition import (
_GeneralizedErosionDeposition,
)
from landlab.utils.return_array import return_array_at_node
from .ext.calc_qs import calculate_qs_in
ROOT2 = np.sqrt(2.0) # syntactic sugar for precalculated square root of 2
TIME_STEP_FACTOR = 0.5 # factor used in simple subdivision solver
[docs]
class Space(_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 <https://www.geosci-model-dev.net/10/4577/2017/>`_, 2017.
Unlike other some other fluvial erosion componets in Landlab, in this
component (and :py:class:`~landlab.components.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
""" # noqa: B950
_name = "Space"
_unit_agnostic = True
_info = {
"flow__link_to_receiver_node": {
"dtype": int,
"intent": "in",
"optional": False,
"units": "-",
"mapping": "node",
"doc": "ID of link downstream of each node, which carries the discharge",
},
"flow__receiver_node": {
"dtype": int,
"intent": "in",
"optional": False,
"units": "-",
"mapping": "node",
"doc": "Node array of receivers (node that receives flow from current node)",
},
"flow__upstream_node_order": {
"dtype": int,
"intent": "in",
"optional": False,
"units": "-",
"mapping": "node",
"doc": "Node array containing downstream-to-upstream ordered list of node IDs",
},
"sediment__influx": {
"dtype": float,
"intent": "out",
"optional": False,
"units": "m3/s",
"mapping": "node",
"doc": "Sediment flux (volume per unit time of sediment entering each node)",
},
"sediment__outflux": {
"dtype": float,
"intent": "out",
"optional": False,
"units": "m3/s",
"mapping": "node",
"doc": "Sediment flux (volume per unit time of sediment leaving each node)",
},
"soil__depth": {
"dtype": float,
"intent": "inout",
"optional": False,
"units": "m",
"mapping": "node",
"doc": "Depth of soil or weathered bedrock",
},
"surface_water__discharge": {
"dtype": float,
"intent": "in",
"optional": False,
"units": "m**3/s",
"mapping": "node",
"doc": "Volumetric discharge of surface water",
},
"topographic__elevation": {
"dtype": float,
"intent": "inout",
"optional": False,
"units": "m",
"mapping": "node",
"doc": "Land surface topographic elevation",
},
"topographic__steepest_slope": {
"dtype": float,
"intent": "in",
"optional": False,
"units": "-",
"mapping": "node",
"doc": "The steepest *downhill* slope",
},
}
_cite_as = """
@Article{gmd-10-4577-2017,
AUTHOR = {Shobe, C. M. and Tucker, G. E. and Barnhart, K. R.},
TITLE = {The SPACE~1.0 model: a~Landlab component for 2-D calculation
of sediment transport, bedrock erosion, and landscape evolution},
JOURNAL = {Geoscientific Model Development},
VOLUME = {10},
YEAR = {2017},
NUMBER = {12},
PAGES = {4577--4604},
URL = {https://www.geosci-model-dev.net/10/4577/2017/},
DOI = {10.5194/gmd-10-4577-2017}
}"""
[docs]
def __init__(
self,
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=DEFAULT_MINIMUM_TIME_STEP,
):
"""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.
"""
if grid.at_node["flow__receiver_node"].size != grid.size("node"):
raise NotImplementedError(
"A route-to-multiple flow director has been "
"run on this grid. The landlab development team has not "
"verified that SPACE is compatible with "
"route-to-multiple methods. Please open a GitHub Issue "
"to start this process."
)
super().__init__(
grid,
m_sp=m_sp,
n_sp=n_sp,
F_f=F_f,
v_s=v_s,
dt_min=dt_min,
discharge_field=discharge_field,
)
if phi >= 1.0:
raise ValueError("Porosity must be < 1.0")
if phi < 0.0:
raise ValueError("Porosity must be > 0.0")
self._phi = float(phi)
self._porosity_factor = 1.0 / (1.0 - self._phi)
# space specific inits
self._H_star = H_star
self._soil__depth = grid.at_node["soil__depth"]
if "bedrock__elevation" in grid.at_node:
self._bedrock__elevation = grid.at_node["bedrock__elevation"]
else:
self._bedrock__elevation = grid.add_zeros(
"bedrock__elevation", at="node", dtype=float
)
self._bedrock__elevation[:] = (
self._topographic__elevation - self._soil__depth
)
self._sed_erosion_term = np.zeros(grid.number_of_nodes)
self._br_erosion_term = np.zeros(grid.number_of_nodes)
self._Es = np.zeros(grid.number_of_nodes)
self._Er = np.zeros(grid.number_of_nodes)
# K's and critical values can be floats, grid fields, or arrays
# use setters defined below
self.K_sed = K_sed
self.K_br = K_br
self._sp_crit_sed = return_array_at_node(grid, sp_crit_sed)
self._sp_crit_br = return_array_at_node(grid, sp_crit_br)
# Handle option for solver
if solver == "basic":
self.run_one_step = self.run_one_step_basic
elif solver == "adaptive":
self.run_one_step = self.run_with_adaptive_time_step_solver
self._time_to_flat = np.zeros(grid.number_of_nodes)
self._time_to_zero_alluv = np.zeros(grid.number_of_nodes)
self._dzdt = np.zeros(grid.number_of_nodes)
else:
raise ValueError(
"Parameter 'solver' must be one of: " + "'basic', 'adaptive'"
)
@property
def K_br(self):
"""Erodibility of bedrock(units depend on m_sp)."""
return self._K_br
@K_br.setter
def K_br(self, new_val):
self._K_br = return_array_at_node(self._grid, new_val)
@property
def K_sed(self):
"""Erodibility of sediment(units depend on m_sp)."""
return self._K_sed
@K_sed.setter
def K_sed(self, new_val):
self._K_sed = return_array_at_node(self._grid, new_val)
def _calc_erosion_rates(self):
"""Calculate erosion rates."""
# if sp_crits are zero, then this colapses to correct all the time.
if self._n_sp == 1.0:
S_to_the_n = self._slope
else:
S_to_the_n = np.power(self._slope, self._n_sp)
omega_sed = self._K_sed * self._Q_to_the_m * S_to_the_n
omega_br = self._K_br * self._Q_to_the_m * S_to_the_n
omega_sed_over_sp_crit = np.divide(
omega_sed,
self._sp_crit_sed,
out=np.zeros_like(omega_sed),
where=self._sp_crit_sed != 0,
)
omega_br_over_sp_crit = np.divide(
omega_br,
self._sp_crit_br,
out=np.zeros_like(omega_br),
where=self._sp_crit_br != 0,
)
self._sed_erosion_term = omega_sed - self._sp_crit_sed * (
1.0 - np.exp(-omega_sed_over_sp_crit)
) / (
1 - self._phi
) # convert from a volume to a mass flux.
self._br_erosion_term = omega_br - self._sp_crit_br * (
1.0 - np.exp(-omega_br_over_sp_crit)
)
H_over_Hstar = self._soil__depth / self._H_star
self._Es = self._sed_erosion_term * (1.0 - np.exp(-H_over_Hstar))
self._Er = self._br_erosion_term * np.exp(-H_over_Hstar)
@property
def Es(self):
"""Sediment erosion term."""
return self._Es
@property
def Er(self):
"""Bedrock erosion term."""
return self._Er
@property
def H(self):
"""Sediment thickness."""
return self._H
def _calc_qs_in_and_depo_rate(self):
# Choose a method for calculating erosion:
self._calc_hydrology()
self._calc_erosion_rates()
is_flooded_core_node = self._get_flooded_core_nodes()
self._Es[is_flooded_core_node] = 0.0
self._Er[is_flooded_core_node] = 0.0
self._sed_erosion_term[is_flooded_core_node] = 0.0
self._br_erosion_term[is_flooded_core_node] = 0.0
self.sediment_influx[:] = 0
self._depo_rate[:] = 0.0
# iterate top to bottom through the stack, calculate qs
# cythonized version of calculating qs_in
calculate_qs_in(
np.flipud(self._stack),
self._flow_receivers,
self._cell_area_at_node,
self._q,
self._qs,
self.sediment_influx,
self._Es,
self._Er,
self._v_s,
self._F_f,
)
self._depo_rate[self._q > 0] = self._qs[self._q > 0] * (
self._v_s / self._q[self._q > 0]
)
if not self._depressions_are_handled(): # all sed dropped here
self._depo_rate[is_flooded_core_node] = (
self.sediment_influx[is_flooded_core_node]
/ self._cell_area_at_node[is_flooded_core_node]
)
return is_flooded_core_node
[docs]
def run_one_step_basic(self, dt=1.0):
"""Calculate change in rock and alluvium thickness for a time period
'dt'.
Parameters
----------
dt : float
Model timestep [T]
"""
self._calc_qs_in_and_depo_rate()
cores = self._grid.core_nodes
H0 = self._soil__depth.copy()
# now, the analytical solution to soil thickness in time:
# need to distinguish D=kqS from all other cases to save from blowup!
# distinguish cases:
# there is a numerical blow up as well, which is not discussed in the
# paper. when H>>H* the analytical solution (Eq 34) has a term that is
# exp (H/H*) This can becom infinite. We will consider the case of
# H/H*>50 as distinct.
H_over_H_star = self._soil__depth / self._H_star
too_thick = H_over_H_star > 100
H_over_H_star[too_thick] = 100
no_entrainment = self._sed_erosion_term <= 0.0 # this will include pits.
blowup = (
(self._depo_rate == self._sed_erosion_term)
& (self._sed_erosion_term > 0.0)
& (~no_entrainment)
)
full_equation = (~blowup) & (~no_entrainment)
# First do blowup.
# this is Space paper Eq 34. This example has the same issues with
# very thick sed. If sed is very thick AND ero=depo, there is no change.
self._soil__depth[blowup * (~too_thick)] = self._H_star * np.log(
((self._sed_erosion_term[blowup]) / self._H_star) * dt
+ np.exp(H_over_H_star[blowup])
)
# Second, do no entrainment of sediment. This is equation 35.
self._soil__depth[no_entrainment] += (
self._depo_rate[no_entrainment] / (1 - self._phi)
) * dt
# Treat the case of very thick sediments and ero != depo.
self._soil__depth[full_equation * too_thick] += (
(
self._depo_rate[full_equation * (full_equation * too_thick)]
/ (1 - self._phi)
)
- (self._sed_erosion_term[full_equation * too_thick] / (1 - self._phi))
) * dt
# Finally do the full equation (Eq 32) when not too thick.
self._soil__depth[full_equation * (~too_thick)] = self._H_star * np.log(
(
1
/ (
(self._depo_rate[full_equation * (~too_thick)] / (1 - self._phi))
/ (
self._sed_erosion_term[full_equation * (~too_thick)]
/ (1 - self._phi)
)
- 1
)
)
* (
np.exp(
(
self._depo_rate[full_equation * (~too_thick)] / (1 - self._phi)
- (
self._sed_erosion_term[full_equation * (~too_thick)]
/ (1 - self._phi)
)
)
* (dt / self._H_star)
)
* (
(
(
self._depo_rate[full_equation * (~too_thick)]
/ (1 - self._phi)
/ (
self._sed_erosion_term[full_equation * (~too_thick)]
/ (1 - self._phi)
)
)
- 1
)
* np.exp(H_over_H_star[full_equation * (~too_thick)])
+ 1
)
- 1
)
)
# Equation 12 gives dRdt, and Equation 36 gives R. However, we don't
# include dH/dt within timestep in integrating for R.
# This matters when we are starting with very little soil and increasing.
# to do this right I think we need a numerical integral for
# R(t). Unfortunately we have three cases for H(t). These are now
# implemented in the function _dRdt. This loop and the solver can be
# cythonized. But not today.
dR = self._grid.zeros(at="node")
for idx in self.grid.core_nodes:
args = (
self._br_erosion_term[idx],
1.0 / self._H_star,
self._depo_rate[idx] / (1.0 - self._phi),
self._sed_erosion_term[idx] / (1.0 - self._phi),
H0[idx],
)
dR[idx] = quad(_dRdt, 0, dt, args)[0]
self._bedrock__elevation += dR
# finally, determine topography by summing bedrock and soil
self._topographic__elevation[cores] = (
self._bedrock__elevation[cores] + self._soil__depth[cores]
)
[docs]
def run_with_adaptive_time_step_solver(self, dt=1.0):
"""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])
"""
# Initialize remaining_time, which records how much of the global time
# step we have yet to use up.
remaining_time = dt
z = self._grid.at_node["topographic__elevation"]
br = self._grid.at_node["bedrock__elevation"]
H = self._grid.at_node["soil__depth"]
r = self._flow_receivers
cores = self._grid.core_nodes
first_iteration = True
is_flooded_core_node = self._get_flooded_core_nodes()
# Outer WHILE loop: keep going until time is used up
while remaining_time > 0.0:
# Update all the flow-link slopes.
#
# For the first iteration, we assume this has already been done
# outside the component (e.g., by flow router), but we need to do
# it ourselves on subsequent iterations.
if not first_iteration:
# update the link slopes
self._update_flow_link_slopes()
# update where nodes are flooded. This shouuldn't happen because
# of the dynamic timestepper, but just in case, we update here.
is_flooded_core_node[self._slope < 0] = True
else:
first_iteration = False
is_flooded_core_node = (
self._calc_qs_in_and_depo_rate()
) # THIS IS THE SPEED BOTTLENECK
# Now look at upstream-downstream node pairs, and recording the
# time it would take for each pair to flatten. Take the minimum.
self._dzdt[cores] = self._depo_rate[cores] * self._porosity_factor - (
self._Es[cores] + self._Er[cores]
)
rocdif = self._dzdt - self._dzdt[r]
zdif = z - z[r]
self._time_to_flat[:] = remaining_time
converging = np.where(rocdif < 0.0)[0]
self._time_to_flat[converging] = -(
TIME_STEP_FACTOR * zdif[converging] / rocdif[converging]
)
self._time_to_flat[np.where(zdif <= 0.0)[0]] = remaining_time
# From this, find the maximum stable time step with regard to slope
# evolution.
dt_max1 = np.amin(self._time_to_flat)
# Next we consider time to exhaust regolith
self._time_to_zero_alluv[:] = remaining_time
# poof deposition by phi
dHdt = self._porosity_factor * (self._depo_rate - self._Es)
decreasing_H = np.where(dHdt < 0.0)[0]
self._time_to_zero_alluv[decreasing_H] = -(
TIME_STEP_FACTOR * H[decreasing_H] / dHdt[decreasing_H]
)
# Now find the smallest time that would lead to near-empty alluv
dt_max2 = np.amin(self._time_to_zero_alluv)
# Take the smaller of the limits
dt_max = max(self._dt_min, min(dt_max1, dt_max2))
# Now a vector operation: apply dzdt and dhdt to all nodes
br[cores] -= self._Er[cores] * dt_max
H[cores] += dHdt[cores] * dt_max
z[cores] = br[cores] + H[cores]
# Update remaining time and continue
remaining_time -= dt_max
def _dRdt(t, a, b, c, d, H0):
"""
dRdt = a * exp(-b * H (t))
a = Kr q S*n
b = 1/H*
c = VQs/(Q * (1-phi))
d = Ks q Sn
if d <= 0:
H (t) = (1/b) ln [ d*b*t + exp(H0/H*) ]
"""
# truncate H/H*
if b * H0 > 100:
too_thick = True
else:
too_thick = False
bH0 = min(b * H0, 100)
# set cases for H depending on values of constants.
# deposition only case.
if d <= 0:
H = H0 + (c * t)
else:
# sediment deposition = sediment entrainment.
if c == d:
if too_thick:
H = H0
else:
H = (1 / b) * np.log((d * b * t) + np.exp(bH0))
# full equation 32
else:
if too_thick:
H = H0 + (c - d) * t
else:
term1 = 1.0 / ((c / d) - 1)
term2 = np.exp((c - d) * t * b)
term3 = (c / d - 1) * np.exp(bH0) + 1
interior = term1 * (term2 * (term3) - 1)
H = (1 / b) * np.log(interior)
# calculate dRdt
dRdt = -a * np.exp(-b * H)
return dRdt