This Jupyter notebook illustrates running the transport-length-model hillslope diffusion component in a simple example.
This component uses an approach similar to the Davy and Lague (2009) equation for fluvial erosion and transport, and applies it to hillslope diffusion. The formulation and implementation were inspired by Carretier et al. (2016); see this paper and references therein for justification.
The elevation $z$ of a point of the landscape (such as a grid node) changes according to:
\begin{equation} \frac{\partial z}{\partial t} = -\epsilon + D + U \tag{1}\label{eq:1}, \end{equation}and we define: \begin{equation} D = \frac{q_s}{L} \tag{2}\label{eq:2}, \end{equation}
where $\epsilon$ is the local erosion rate [L/T], $D$ the local deposition rate [L/T], $U$ the uplift (or subsidence) rate [L/T], $q_s$ the incoming sediment flux per unit width [L$^2$/T] and $L$ is the transport length.
We specify the erosion rate $\epsilon$ and the transport length $L$:
\begin{equation} \epsilon = \kappa S \tag{3}\label{eq:3} \end{equation}\begin{equation} L = \frac{dx}{1-({S}/{S_c})^2} \tag{4}\label{eq:4} \end{equation}where $\kappa$ [L/T] is an erodibility coefficient, $S$ is the local slope [L/L] and $S_c$ is the critical slope [L/L].
Thus, the elevation variation results from the difference between local rates of detachment and deposition.
The detachment rate is proportional to the local gradient. However, the deposition rate ($q_s/L$) depends on the local slope and the critical slope:
Previous models typically use a "non-linear" diffusion model proposed by different authors (e.g. Andrews and Hanks, 1985; Hanks, 1999; Roering et al., 1999) and supported by $^{10}$Be-derived erosion rates (e.g. Binnie et al., 2007) or experiments (Roering et al., 2001). It is usually presented in the followin form:
$ $
\begin{equation} \frac{\partial z}{\partial t} = \frac{\partial q_s}{\partial x} \tag{5}\label{eq:5} \end{equation}$ $ \begin{equation} q_s = \frac{\kappa' S}{1-({S}/{S_c})^2} \tag{6}\label{eq:6} \end{equation}
where $\kappa'$ [L$^2$/T] is a diffusion coefficient.
This description is thus based on the definition of a flux of transported sediment parallel to the slope:
Despite these conceptual differences, equations ($\ref{eq:3}$) and ($\ref{eq:4}$) predict similar topographic evolution to the 'non-linear' diffusion equations for $\kappa' = \kappa dx$, as shown in the following example.
First, we import what we'll need:
import numpy as np
from matplotlib.pyplot import figure, plot, show, title, xlabel, ylabel
from landlab import RasterModelGrid
from landlab.components import FlowDirectorSteepest, TransportLengthHillslopeDiffuser
from landlab.plot import imshow_grid
# to plot figures in the notebook:
%matplotlib inline
Make a grid and set boundary conditions:
mg = RasterModelGrid(
(20, 20),
xy_spacing=50.0) # raster grid with 20 rows, 20 columns and dx=50m
z = np.random.rand(mg.size('node')) # random noise for initial topography
mg.add_field("topographic__elevation", z, at="node")
mg.set_closed_boundaries_at_grid_edges(
False, True, False,
True) # N and S boundaries are closed, E and W are open
Set the initial and run conditions:
total_t = 2000000.0 # total run time (yr)
dt = 1000.0 # time step (yr)
nt = int(total_t // dt) # number of time steps
uplift_rate = 0.0001 # uplift rate (m/yr)
kappa = 0.001 # erodibility (m/yr)
Sc = 0.6 # critical slope
Instantiate the components: The hillslope diffusion component must be used together with a flow router/director that provides the steepest downstream slope for each node, with a D4 method (creates the field topographic__steepest_slope at nodes).
fdir = FlowDirectorSteepest(mg)
tl_diff = TransportLengthHillslopeDiffuser(mg,
erodibility=kappa,
slope_crit=Sc)
Run the components for 2 Myr and trace an East-West cross-section of the topography every 100 kyr:
for t in range(nt):
fdir.run_one_step()
tl_diff.run_one_step(dt)
z[mg.core_nodes] += uplift_rate * dt # add the uplift
# add some output to let us see we aren't hanging:
if t % 100 == 0:
print(t * dt)
# plot east-west cross-section of topography:
x_plot = range(0, 1000, 50)
z_plot = z[100:120]
figure("cross-section")
plot(x_plot, z_plot)
figure("cross-section")
title("East-West cross section")
xlabel("x (m)")
ylabel("z (m)")
And plot final topography:
figure("final topography")
im = imshow_grid(mg,
"topographic__elevation",
grid_units=["m", "m"],
var_name="Elevation (m)")
This behaviour corresponds to the evolution observed using a classical non-linear diffusion model.
In this example, we show that when the slope is steep ($S \ge S_c$), the transport-length hillsope diffusion simulates mass wasting, with long transport distances.
First, we create a grid: the western half of the grid is flat at 0 m of elevation, the eastern half is a 45-degree slope.
# Create grid and topographic elevation field:
mg2 = RasterModelGrid((20, 20), xy_spacing=50.0)
z = np.zeros(mg2.number_of_nodes)
z[mg2.node_x > 500] = mg2.node_x[mg2.node_x > 500] / 10
mg2.add_field("topographic__elevation", z, at="node")
# Set boundary conditions:
mg2.set_closed_boundaries_at_grid_edges(False, True, False, True)
# Show initial topography:
im = imshow_grid(mg2,
"topographic__elevation",
grid_units=["m", "m"],
var_name="Elevation (m)")
# Plot an east-west cross-section of the initial topography:
z_plot = z[100:120]
x_plot = range(0, 1000, 50)
figure(2)
plot(x_plot, z_plot)
title("East-West cross section")
xlabel("x (m)")
ylabel("z (m)")
Set the run conditions:
total_t = 1000000.0 # total run time (yr)
dt = 1000.0 # time step (yr)
nt = int(total_t // dt) # number of time steps
Instantiate the components:
fdir = FlowDirectorSteepest(mg2)
tl_diff = TransportLengthHillslopeDiffuser(mg2,
erodibility=0.001,
slope_crit=0.6)
Run for 1 Myr, plotting the cross-section regularly:
for t in range(nt):
fdir.run_one_step()
tl_diff.run_one_step(dt)
# add some output to let us see we aren't hanging:
if t % 100 == 0:
print(t * dt)
z_plot = z[100:120]
figure(2)
plot(x_plot, z_plot)
The material is diffused from the top and along the slope and it accumulates at the bottom, where the topography flattens.
As a comparison, the following code uses linear diffusion on the same slope:
# Import Linear diffuser:
from landlab.components import LinearDiffuser
# Create grid and topographic elevation field:
mg3 = RasterModelGrid((20, 20), xy_spacing=50.0)
z = np.ones(mg3.number_of_nodes)
z[mg.node_x > 500] = mg.node_x[mg.node_x > 500] / 10
mg3.add_field("topographic__elevation", z, at="node")
# Set boundary conditions:
mg3.set_closed_boundaries_at_grid_edges(False, True, False, True)
# Instantiate components:
fdir = FlowDirectorSteepest(mg3)
diff = LinearDiffuser(mg3, linear_diffusivity=0.1)
# Set run conditions:
total_t = 1000000.0
dt = 1000.0
nt = int(total_t // dt)
# Run for 1 Myr, plotting east-west cross-section regularly:
for t in range(nt):
fdir.run_one_step()
diff.run_one_step(dt)
# add some output to let us see we aren't hanging:
if t % 100 == 0:
print(t * dt)
z_plot = z[100:120]
figure(2)
plot(x_plot, z_plot)