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.) TrueNow 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) TrueThe 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)) FalseNow, 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]) FalseReferences
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.