DepthDependentTaylorDiffuser: Model depth dependent non-linear soil creep after Ganti et al. (2012) and Johnstone and Hilley (2014)#

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

Bases: Component

This component implements a depth-dependent Taylor series diffusion rule, combining concepts of Ganti et al. (2012) and Johnstone and Hilley (2014).

Hillslope sediment flux uses a Taylor Series expansion of the Andrews- Bucknam formulation of nonlinear hillslope flux derived following following Ganti et al., 2012 with a depth dependent component inspired Johnstone and Hilley (2014). The flux \(q_s\) is given as:

\[q_s = - K H_* \nabla \eta ( 1 + (S/S_c)^2 + (S/S_c)^4 + .. + (S/S_c)^2(n-1) ) (1 - exp( - H / H_*)\]

where \(K\) is a transport velocity coefficient, \(\eta\) is land surface elevation, \(S\) is the slope gradient (defined as positive downward), \(S_c\) is the critical slope, \(n\) is the number of terms, \(H\) is the soil depth on links, and \(H_*\) is the soil transport decay depth.

The default behavior uses two terms to produce a slope dependence as described by Equation 6 of Ganti et al. (2012).

This component will ignore soil thickness located at non-core nodes.

Examples

First lets make a simple example with flat topography.

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import ExponentialWeatherer
>>> from landlab.components import DepthDependentTaylorDiffuser
>>> mg = RasterModelGrid((5, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg)
>>> expweath.calc_soil_prod_rate()
>>> np.allclose(mg.at_node['soil_production__rate'][mg.core_nodes], 1.)
True
>>> DDdiff.run_one_step(2.)
>>> np.allclose(mg.at_node['topographic__elevation'][mg.core_nodes], 0.)
True
>>> np.allclose(mg.at_node['bedrock__elevation'][mg.core_nodes], -2.)
True
>>> np.allclose(mg.at_node['soil__depth'][mg.core_nodes], 2.)
True

Now a more complicated example with a slope.

>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> z += mg.node_x.copy()
>>> BRz += mg.node_x / 2.
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg)
>>> expweath.calc_soil_prod_rate()
>>> np.allclose(
...     mg.at_node['soil_production__rate'][mg.core_nodes],
...     np.array([ 0.60653066, 0.36787944, 0.22313016]))
True
>>> DDdiff.run_one_step(0.1)
>>> np.allclose(
...     mg.at_node['topographic__elevation'][mg.core_nodes],
...     np.array([ 1.04773024, 2.02894986, 3.01755898]))
True
>>> np.allclose(mg.at_node['bedrock__elevation'][mg.core_nodes],
...     np.array([ 0.43934693, 0.96321206, 1.47768698]))
True
>>> np.allclose(mg.at_node['soil__depth'], z - BRz)
True

The DepthDependentTaylorDiffuser makes and moves soil at a rate proportional to slope, this means that there is a characteristic time scale for soil transport and an associated stability criteria for the timestep. The maximum characteristic time scale, \(De_{max}\), is given as a function of the hillslope diffustivity, \(D\), the maximum slope, \(S_{max}\), and the critical slope \(S_c\).

\[De_{max} = D \left( 1 + \left( \frac{S_{max}{S_c}\right )^2 + \left( \frac{S_{max}{S_c}\right )^4 + \dots + \left( \frac{S_{max}{S_c}\right )^{( 2 * ( n - 1 ))} \right)\]

The maximum stable time step is given by

\[dtmax = courant_factor * dx * dx / Demax\]

Where the courant factor is a user defined scale (default is 0.2)

The DepthDependentTaylorDiffuser has a boolean flag that permits a user to be warned if timesteps are too large for the slopes in the model grid (if_unstable = ‘warn’) and a boolean flag that turns on dynamic timestepping (dynamic_dt = False).

>>> DDdiff = DepthDependentTaylorDiffuser(mg, if_unstable='warn')
>>> DDdiff.run_one_step(2.)
Topographic slopes are high enough such that the Courant condition is
exceeded AND you have not selected dynamic timestepping with
dynamic_dt=True. This may lead to infinite and/or nan values for slope,
elevation, and soil depth. Consider using a smaller time step or dynamic
timestepping. The Courant condition recommends a timestep of
0.0953407607307 or smaller.

Alternatively you can specify if_unstable=’raise’, and a Runtime Error will be raised if this condition is not met.

Next, lets do an example with dynamic timestepping.

>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')

We’ll use a steep slope and very little soil.

>>> z += mg.node_x.copy()**2
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)

Lets try to move the soil with a large timestep. Without dynamic time steps, this gives a warning that we’ve exceeded the dynamic timestep size and should use a smaller timestep. We could either use the smaller timestep, or specify that we want to use a dynamic timestep.

>>> DDdiff = DepthDependentTaylorDiffuser(mg, if_unstable='warn', dynamic_dt=False)
>>> expweath.calc_soil_prod_rate()
>>> DDdiff.run_one_step(10)
Topographic slopes are high enough such that the Courant condition is
exceeded AND you have not selected dynamic timestepping with
dynamic_dt=True. This may lead to infinite and/or nan values for slope,
elevation, and soil depth. Consider using a smaller time step or dynamic
timestepping. The Courant condition recommends a timestep of
0.004 or smaller.

Now, we’ll re-build the grid and do the same example with dynamic timesteps.

>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> z += mg.node_x.copy()**2
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg, if_unstable='warn', dynamic_dt=True)
>>> expweath.calc_soil_prod_rate()
>>> DDdiff.run_one_step(10)
>>> np.any(np.isnan(z))
False

Now, we’ll test that changing the transport decay depth behaves as expected.

>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> z += mg.node_x.copy()**0.5
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg, soil_transport_decay_depth = 0.1)
>>> DDdiff.run_one_step(1)
>>> soil_decay_depth_point1 = mg.at_node['topographic__elevation'][mg.core_nodes]
>>> z[:] = 0
>>> z += mg.node_x.copy()**0.5
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> DDdiff = DepthDependentTaylorDiffuser(mg, soil_transport_decay_depth = 1.0)
>>> DDdiff.run_one_step(1)
>>> soil_decay_depth_1 = mg.at_node['topographic__elevation'][mg.core_nodes]
>>> np.greater(soil_decay_depth_1[1], soil_decay_depth_point1[1])
False

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

Ganti, V., Passalacqua, P., Foufoula-Georgiou, E. (2012). A sub-grid scale closure for nonlinear hillslope sediment transport models Journal of Geophysical Research: Earth Surface 117(F2). https://dx.doi.org/10.1029/2011jf002181

Johnstone, S., Hilley, G. (2015). Lithologic control on the form of soil-mantled hillslopes Geology 43(1), 83-86. https://doi.org/10.1130/G36052.1

Initialize the DepthDependentTaylorDiffuser.

Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • linear_diffusivity (float, optional, DEPRECATED) – Hillslope diffusivity / decay depth, m/yr Default = 1.0

  • slope_crit (float, optional) – Critical gradient parameter, m/m Default = 1.0

  • soil_transport_decay_depth (float, optional) – characteristic transport soil depth, m Default = 1.0

  • nterms (int, optional. default = 2) – number of terms in the Taylor expansion. Two terms (default) gives the behavior described in Ganti et al. (2012).

  • dynamic_dt (bool, optional, default = False) – Whether internal timestepping is used.

  • if_unstable (str, optional, default = "pass") – What to do if unstable (options are “pass”, “raise”, “warn”)

  • courant_factor (float, optional, default = 0.2) – Courant factor for timestep calculation.

  • soil_transport_velocity (float, optional, default = 1.0) – Velocity parameter for soil transport, m/yr. Diffusivity is the product of this parameter and soil_transport_decay_depth.

__init__(grid, linear_diffusivity=None, slope_crit=1.0, soil_transport_decay_depth=1.0, nterms=2, dynamic_dt=False, if_unstable='pass', courant_factor=0.2, soil_transport_velocity=1.0)[source]#

Initialize the DepthDependentTaylorDiffuser.

Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • linear_diffusivity (float, optional, DEPRECATED) – Hillslope diffusivity / decay depth, m/yr Default = 1.0

  • slope_crit (float, optional) – Critical gradient parameter, m/m Default = 1.0

  • soil_transport_decay_depth (float, optional) – characteristic transport soil depth, m Default = 1.0

  • nterms (int, optional. default = 2) – number of terms in the Taylor expansion. Two terms (default) gives the behavior described in Ganti et al. (2012).

  • dynamic_dt (bool, optional, default = False) – Whether internal timestepping is used.

  • if_unstable (str, optional, default = "pass") – What to do if unstable (options are “pass”, “raise”, “warn”)

  • courant_factor (float, optional, default = 0.2) – Courant factor for timestep calculation.

  • soil_transport_velocity (float, optional, default = 1.0) – Velocity parameter for soil transport, m/yr. Diffusivity is the product of this parameter and soil_transport_decay_depth.

run_one_step(dt)[source]#
Parameters:

dt (float (time)) – The imposed timestep.

soilflux(dt)[source]#

Calculate soil flux for a time period ‘dt’.

Parameters:

dt (float (time)) – The imposed timestep.