Note

This page was generated from a jupyter notebook.

Component Overview: DepthDependentTaylorDiffuser

Introduction and background

This tutorial introduces the DepthDependentTaylorDiffuser component, which we’ll refer to here as “DDTD”. The DDTD component models the process of downslope soil creep and its role in modifying topography. It combines the mathematics behind two other components: DepthDependentDiffuser and TaylorNonLinearDiffuser. The component is described (as one element in the terrainBento package) in Barnhart et al. (2019), which is the appropriate paper to cite for it.

Theory

Consider a topographic surface in which the elevation at any time \(t\) and horizontal position \((x,y)\) is \(\eta (x,y,t)\). The thickness of the mobile soil layer is \(H(x,y,t)\). Let \(\mathbf{q}_s\) be a 2D vector that represents the rate of soil volume flow per unit slope width (with dimensions of length squared per time; we’ll assume that \(\mathbf{q}_s\) represents a “bulk” flux that includes pore spaces between soil grains). In the absence of any “local” input sources (such as weathering of rock) or output (such as removal by wash erosion), conservation of mass dictates that:

\[\frac{\partial \eta}{\partial t} = -\nabla \cdot \mathbf{q}_s\]

The DDTD component represents the soil flux as:

\[\mathbf{q}_s = K H_* \mathbf{S} [1 - \exp ( - H / H_*)] [ 1 + (S/S_c)^2 + (S/S_c)^4 + ... + (S/S_c)^2(n-1) ]\]

where \(\mathbf{S} = -\nabla \eta\) is the downslope topographic gradient, and \(S\) is its magnitude. Parameter \(H_*\) is a depth scale that determines how rapidly transport rate decays as the soil thins. Parameter \(K\) is a transport coefficient with dimensions of velocity. The effective diffusion-like coefficient is \(D=KH_*\). This is the effective diffusivity when the soil is much thicker than \(H_*\).

The above can be written slightly more compactly:

\[\mathbf{q}_s = K H_* \mathbf{S} [1 - \exp ( - H / H_*)] \left[ 1 + \sum_{i=1}^N \left(\frac{S}{S_c}\right)^{2i} \right]\]

where \(i\) is the number of additional terms desired. If \(i=0\), the expression is the same as the depth-dependent, slope-linear transport function implemented by the DepthDependentDiffuser component and described, for example, by Johnstone and Hilley (2015).

The use of a truncated Taylor series is meant to approximate the Andrews-Bucknam transport function (e.g., Roering et al., 1999) while avoiding that equation’s blow-up at \(S=S_c\); the idea of using a truncated Taylor series comes from Ganti et al. (2012).

Numerical implementation

The component uses an explicit finite-volume solution method. Soil flux values are calculated from the gradient values on the active links, using the grid method calc_grad_at_link. Flux divergence is then calculated using the grid method calc_flux_div_at_node. The calculation updates soil thickness, bedrock elevation (using the user-provided values of the soil_production__rate field), and total elevation as the sum of the two.

An optional dynamic timestep capability will check the local Courant condition (which can vary in time and space when nonlinear terms are included) and sub-divide the user-specified time step as needed to ensure stability.

Examples

Needed imports

Like all Landlab components, DDTD requires a grid object on which to operate, so for this example we’ll import RasterModelGrid as well as the component itself.

import matplotlib.pyplot as plt
import numpy as np
from tqdm.notebook import trange

from landlab import RasterModelGrid
from landlab.components import DepthDependentTaylorDiffuser

Example 1: equilibrium hillslope profile with linear diffusion

For the first example, we’ll use a long and skinny grid to effectively create a 1D domain. We’ll test the ability of DDTD to reduce to a simple linear, depth-independent diffusive model when \(i=0\) and \(H \gg H_*\). We’ll impose (relative) rock uplift by raising the interior of the domain at a specified rate \(U\) relative to the fixed boundary nodes on either side. The expectation is that:

\[\frac{d\eta}{dx} = -\frac{U}{D}x\]

where \(x\) is distance from the ridge top (because the ridge top will form in the middle of the domain, \(x<0\) on the left and \(x>0\) on the right). Integrating this, we get

\[\eta = -\frac{U}{2D} x^2 + C\]

We can evaluate the integration constant by noting that \(\eta = 0\) at \(x = \pm L\), where \(L\) is the distance from base to crest. Therefore,

\[\boxed{\eta = \frac{U}{2D} \left( L^2 - x^2 \right)}\]

We’ll test this using a hill that is 100 m long (51 nodes, two of which are fixed boundaries, with 2 m spacing between them; 50 m from base to crest on each side), a soil layer that is much thicker than the characteristic decay depth \(H^*\), a transport coefficient of 0.01 m\(^2\)/y, and an uplift rate of 0.0001 m/y. With these parameters, the predicted ridge height (at \(x=0\)) is calculated below.

# define parameters
L = 50.0  # distance from base to ridgeline, m
dx = 2.0  # node spacing, m
Hstar = 0.1  # characteristic transport depth, m
V0 = 0.1  # transport velocity coefficient, m/y
U = 0.0001  # uplift rate, m/y
H = 100.0  # initial soil thickness, m
num_steps = 20000  # number of time steps

# time step size (calculate using Courant condition for linear diffusion)
D = V0 * Hstar  # effective (maximum) diffusivity
dt = 0.1 * dx * dx / D

# prediction
predicted_crest_height = 0.5 * (U / D) * L * L
print("Crest height should be " + str(predicted_crest_height))

# create grid
grid = RasterModelGrid((3, 51), xy_spacing=dx)
grid.set_closed_boundaries_at_grid_edges(False, True, False, True)

# create fields
elev = grid.add_zeros("topographic__elevation", at="node")  # this is eta
rock = grid.add_zeros("bedrock__elevation", at="node")  # this is eta - H
rock[:] = -H
soil = grid.add_zeros("soil__depth", at="node")  # this is H
soil_production_rate = grid.add_zeros("soil_production__rate", at="node")

# instantiate component
ddtd = DepthDependentTaylorDiffuser(
    grid, soil_transport_velocity=V0, soil_transport_decay_depth=Hstar, nterms=1
)

# run the model in a time loop with uplift applied
for _ in trange(num_steps):
    elev[grid.core_nodes] += U * dt
    rock[grid.core_nodes] += U * dt
    ddtd.run_one_step(dt)

midrow = np.arange(51, 102, dtype=int)
plt.plot(grid.x_of_node[midrow], elev[midrow])
plt.xlabel("Distance (m)")
plt.ylabel("Elevation (m)")

print(np.amax(elev))

Example 2: linear, depth-dependent diffusion

In this example we add a rule for soil production that will limit the soil thickness and hence reduce the transport efficiency. The rate of soil production from bedrock will be:

\[P = P_0 \exp ( H / H_0 )\]

where \(P_0\) is the maximum production rate and \(H_0\) is a characteristic decay depth. In our example, we’ll set \(P_0\) to twice the uplift rate. At equilibrium, the actual production rate \(P = U\), which means that the equilibrium soil thickness can be found from:

\[P = U = 2 U \exp ( -H / H_0 )\]

or

\[H = -H_0 \ln 1/2\]

The effective diffusion coefficient is therefore

\[D_{eff} = D [1 - \exp ( -H / H_* )] = D [1 - \exp ( H_0 \ln 1/2 / H_* )]\]

For the sake of example, we’ll assume \(H_0 = H_*\), so

\[D_{eff} = D / 2\]

and therefore our hill crest should be twice as high as in the prior case.

# new parameter: maximum soil production rate
P0 = 2 * U  # m/yr

# create grid
grid = RasterModelGrid((3, 51), xy_spacing=dx)
grid.set_closed_boundaries_at_grid_edges(False, True, False, True)

# create fields
elev = grid.add_zeros("topographic__elevation", at="node")  # this is eta
rock = grid.add_zeros("bedrock__elevation", at="node")  # this is eta - H
soil = grid.add_zeros("soil__depth", at="node")  # this is H
soil_production_rate = grid.add_zeros("soil_production__rate", at="node")

# instantiate component
ddtd = DepthDependentTaylorDiffuser(
    grid, soil_transport_velocity=V0, soil_transport_decay_depth=Hstar, nterms=1
)

# run the model in a time loop with uplift applied
for _ in trange(num_steps):
    elev[grid.core_nodes] += U * dt
    rock[grid.core_nodes] += U * dt
    soil_production_rate[grid.core_nodes] = P0 * np.exp(-soil[grid.core_nodes] / Hstar)
    ddtd.run_one_step(dt)

plt.plot(grid.x_of_node[midrow], elev[midrow])
plt.xlabel("Distance (m)")
plt.ylabel("Elevation (m)")

Here we haven’t quite reached equilibrium yet, but we can see that the hilltop crest is approaching our expected height of 25 m: twice as high as it would be if the soil flux were not limited by soil thickness.

Example 3: Nonlinear behavior

When we include nonlinear terms in the transport law, we expect to see slopes that become more planar in character. We’ll test this by setting a critical slope value \(S_c = 0.6\) (about 31\(^\circ\)), and using a higher uplift rate. We’ll have two terms, one linear and one cubic. We will also invoke the dynamic_dt option, which allows the component to subdivide each “global” timestep if needed for numerical stability: a useful thing to do because now our Courant condition varies according to slope gradient.

U = 0.0005  # uplift rate, m/yr
Sc = 0.6  # critical slope gradient, m/m
H = 1000.0  # plenty of soil
num_steps = 2000  # number of time steps

# create grid
grid = RasterModelGrid((3, 51), xy_spacing=dx)
grid.set_closed_boundaries_at_grid_edges(False, True, False, True)

# create fields
elev = grid.add_zeros("topographic__elevation", at="node")  # this is eta
rock = grid.add_zeros("bedrock__elevation", at="node")  # this is eta - H
rock[:] = -H
soil = grid.add_zeros("soil__depth", at="node")  # this is H
soil_production_rate = grid.add_zeros("soil_production__rate", at="node")

# instantiate component
ddtd = DepthDependentTaylorDiffuser(
    grid,
    soil_transport_velocity=V0,
    soil_transport_decay_depth=Hstar,
    slope_crit=Sc,
    dynamic_dt=True,
    nterms=2,
)

# run the model in a time loop with uplift applied
for _ in trange(num_steps):
    elev[grid.core_nodes] += U * dt
    rock[grid.core_nodes] += U * dt
    ddtd.run_one_step(dt)

plt.plot(grid.x_of_node[midrow], elev[midrow])
plt.xlabel("Distance (m)")
plt.ylabel("Elevation (m)")

The resulting hill is taller (due to the higher uplift rate) and no longer has uniform convexity.

How do we know whether it has reached equilibrium? One way is to inspect the soil flux: it should increase linearly with \(x\), and be zero at the crest. The values at the base of the slope should equal slope length times uplift rate, or 50 m x 0.0005 m/yr = 0.025 m\(^2\)/yr.

active_link_midpts = (
    grid.x_of_node[grid.node_at_link_tail[grid.active_links]] + 0.5 * dx
)
plt.plot(active_link_midpts, grid.at_link["soil__flux"][grid.active_links])
plt.grid(True)
plt.xlabel("Distance (m)")
plt.ylabel("Soil flux (m2/yr)")

So it appears as if we are not quite there, but pretty close.

Example 4: Nonlinear, depth dependent, and 2D

In the final example we’ll use a proper 2D domain, with both a soil-depth dependence and a nonlinear term in the flux law.

U = 0.0002  # uplift rate, m/yr

# create grid
grid = RasterModelGrid((21, 31), xy_spacing=dx)

# create fields
elev = grid.add_zeros("topographic__elevation", at="node")  # this is eta
rock = grid.add_zeros("bedrock__elevation", at="node")  # this is eta - H
soil = grid.add_zeros("soil__depth", at="node")  # this is H
soil_production_rate = grid.add_zeros("soil_production__rate", at="node")

# instantiate component
ddtd = DepthDependentTaylorDiffuser(
    grid,
    soil_transport_velocity=V0,
    soil_transport_decay_depth=Hstar,
    slope_crit=Sc,
    dynamic_dt=True,
    nterms=2,
)

# run the model in a time loop with uplift applied
for _ in trange(num_steps):
    elev[grid.core_nodes] += U * dt
    rock[grid.core_nodes] += U * dt
    soil_production_rate[grid.core_nodes] = P0 * np.exp(-soil[grid.core_nodes] / Hstar)
    ddtd.run_one_step(dt)

grid.imshow(elev)

References

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.

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.

Roering, J. J., Kirchner, J. W., & Dietrich, W. E. (1999). Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology. Water Resources Research, 35(3), 853-870.

---

Generated by nbsphinx from a Jupyter notebook.