Note

This page was generated from a jupyter notebook.

Component Overview: TaylorNonLinearDiffuser

Introduction and background

This tutorial introduces the TaylorNonLinearDiffuser component, which we’ll refer to here as “TNLD”. The TNLD component models the process of downslope soil creep and its role in modifying topography. Inspired by Ganti et al. (2012), it uses a slope-dependent flux law with a user-specified number of terms in a Taylor expansion described in that paper. 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)\). 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 TNLD component represents the soil flux as:

\[\mathbf{q}_s = D \mathbf{S} [ 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 \(D\) is a diffusion-like coefficient with dimensions of length squared per time.

The above can be written slightly more compactly:

\[\mathbf{q}_s = D \mathbf{S} \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 reduces to plain old linear diffusion (see LinearDiffuser).

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 code then updates the elevation field.

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, TNLD 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 landlab import RasterModelGrid
from landlab.components import TaylorNonLinearDiffuser

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 TNLD to reduce to a simple linear, depth-independent diffusive model when \(i=0\). 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
D = 0.01  # diffusion-like coefficient, m2/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)
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 elevation field
elev = grid.add_zeros("topographic__elevation", at="node")  # this is eta

# instantiate component
tnld = TaylorNonLinearDiffuser(grid, linear_diffusivity=D, nterms=1)

# run the model in a time loop with uplift applied
for i in range(num_steps):
    elev[grid.core_nodes] += U * dt
    tnld.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: 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
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

# instantiate component
tnld = TaylorNonLinearDiffuser(
    grid, linear_diffusivity=D, slope_crit=Sc, dynamic_dt=True, nterms=2
)

# run the model in a time loop with uplift applied
for i in range(num_steps):
    elev[grid.core_nodes] += U * dt
    tnld.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.

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.

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.