# The transport-length hillslope diffuser

# The basics:

This component uses an approach similar to Davy and Lague (2009)'s equation for fluvial erosion and transport, and applies it to hillslope diffusion.
Formulation and implementation were inspired by Carretier et al (2016), see this paper and references therein for justification.

## Theory

The elevation z of a point of the landscape (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:
- when $S \ll S_c$, most of the sediment entering a node is deposited there, this is the pure diffusion case. In this case, the sediment flux $q_s$ does not include sediment eroded from above and is thus "local".
- when $S \approx S_c$, *L* becomes infinity and there is no redeposition on the node, the sediments are transferred further downstream.  This behaviour corresponds to mass wasting, grains can travel a long distance before being deposited. In that case, the flux $q_s$ is "non-local" as it incorporates sediments that have both been detached  locally and transited from upslope.
- for an intermediate $S$, there is a prgogressive transition between pure creep and "balistic" transport of the material. This is consistent with experiments (Roering et al., 2001; Gabet and Mendoza, 2012).

## Contrast with the non-linear diffusion model

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:
- when the slope is small, this flux refers to diffusion processes such as aoil creep, rain splash or diffuse runoff
- when the slope gets closer to critical slope, the flux increases dramatically, simulating on average the cumulative effect of mass wasting events.


Despite these conceptual differences, Eq ($\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.


# Example 1:

First, we import what we'll need:


```python
%matplotlib inline

import numpy as np
from matplotlib.pyplot import figure, show, plot, xlabel, ylabel, title
import pymt.models
```

    [33;01m➡ models: Avulsion, Plume, Sedflux3D, Subside, FrostNumber, Ku, ExponentialWeatherer, Flexure, FlowAccumulator, FlowDirectorD8, FlowDirectorDINF, FlowDirectorSteepest, FlowRouter, LinearDiffuser, OverlandFlow, SoilMoisture, StreamPowerEroder, TransportLengthHillslopeDiffuser, Vegetation, Hydrotrend, Child, Cem, Waves[39;49;00m


Set the initial and run conditions:


```python
total_t = 2000000.  # total run time (yr)
dt = 1000.  # 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).


```python
fdir = pymt.models.FlowDirectorSteepest()
tl_diff = pymt.models.TransportLengthHillslopeDiffuser()
```


```python
config_file, config_dir = fdir.setup(
    grid_row_spacing=10.,
    grid_column_spacing=10.,
    grid_rows=100,
    grid_columns=100,
    clock_start=0.0,
    clock_stop=total_t,
    clock_step=dt,
)
fdir.initialize(config_file, config_dir)
```


```python
config_file, config_dir = tl_diff.setup(
    grid_row_spacing=10.,
    grid_column_spacing=10.,
    grid_rows=100,
    grid_columns=100,

    clock_start=0.0,
    clock_stop=total_t,
    clock_step=dt,

    erodibility=kappa,
    slope_crit=Sc,
)
tl_diff.initialize(config_file, config_dir)
```

Set the boundary conditions. The **FlowDirector** component uses a variable called *boundary_condition_flag* to set its boundary conditions. A value of 1, means the boundary is open and sediment is free to leave the grid. A value of 4 means the nodes are closed and so there is no flux through them. The **TransportLengthHillslopeDiffuser** uses these boundary conditions as input so we'll set them both here.


```python
status = fdir.get_value("boundary_condition_flag").reshape((100, 100))
status[:, (0, -1)] = 1 # E and W boundaries are open 
status[(0, -1), :] = 4 # N and S boundaries are closed
```


```python
fdir.set_value("boundary_condition_flag", status)
tl_diff.set_value("boundary_condition_flag", status)
```

Start with an initial surface that's just random noise.


```python
z = np.random.rand(100 * 100)
```


```python
fdir.set_value("topographic__elevation", z)
```

Get the input values for **TransportLengthHillslopeDiffuser** from the flow director.


```python
tl_diff.set_value("topographic__elevation", z)
tl_diff.set_value("flow__receiver_node", fdir.get_value("flow__receiver_node"))
tl_diff.set_value("topographic__steepest_slope", fdir.get_value("topographic__steepest_slope"))
```

Run the components for 2 Myr and trace an East-West cross-section of the topography every 100 kyr:


```python
for t in range(nt - 1):
    fdir.update()
    
    tl_diff.set_value("topographic__elevation", z)
    tl_diff.set_value(
        "flow__receiver_node", fdir.get_value("flow__receiver_node")
    )
    tl_diff.set_value(
        "topographic__steepest_slope",
        fdir.get_value("topographic__steepest_slope"),
    )

    tl_diff.update()
    
    z = tl_diff.get_value("topographic__elevation").reshape((100, 100))
    z[1:-1, 1:-1] += uplift_rate * dt  # add the uplift
    
    fdir.set_value("topographic__elevation", z)
    
    # 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, 10)
        z_plot = z[1, :]
        figure('cross-section')
        plot(x_plot, z_plot)
```

And plot final topography:


```python
tl_diff.quick_plot("topographic__elevation")
```

# Example 2

In this example, we show that when the slope is steep ($S \geq 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 0m of elevation, the eastern half is a 45-degree slope.


```python
fdir = pymt.models.FlowDirectorSteepest()
tl_diff = pymt.models.TransportLengthHillslopeDiffuser()
```


```python
total_t = 1000000.  # total run time (yr)
dt = 1000.  # time step (yr)
nt = int(total_t // dt)  # number of time steps

kappa = 0.001  # erodibility (m / yr)
Sc = 0.6  # critical slope
```


```python
grid_params = {
    "grid_row_spacing": 10.,
    "grid_column_spacing": 10.,
    "grid_rows": 100,
    "grid_columns": 100,
}
clock_params = {
    "clock_start": 0.0,
    "clock_stop": total_t,
    "clock_step": dt,
}
```


```python
config_file, config_dir = fdir.setup(**grid_params, **clock_params)
fdir.initialize(config_file, config_dir)
```


```python
config_file, config_dir = tl_diff.setup(
    **grid_params,
    **clock_params,
    erodibility=kappa,
    slope_crit=Sc,
)
tl_diff.initialize(config_file, config_dir)
```

As before, set the boundary conditions for both components.


```python
status = fdir.get_value("boundary_condition_flag").reshape((100, 100))
status[:, (0, -1)] = 1 # E and W boundaries are open 
status[(0, -1), :] = 4 # N and S boundaries are closed
```


```python
fdir.set_value("boundary_condition_flag", status)
tl_diff.set_value("boundary_condition_flag", status)
```

In this example, we'll use a different initial surface: a dipping plane.


```python
grid = fdir.var_grid("topographic__elevation")
n_vals = fdir.grid_size(grid)
x, y = fdir.grid[0].node_x, fdir.grid[0].node_y
```


```python
z = np.zeros(n_vals)
z[x > 500] = x[x < 490] / 10.0
```

To make sure we've set things up correctly.


```python
fdir.set_value("topographic__elevation", z)
fdir.quick_plot("topographic__elevation")
```

Now time step through the model, plotting things along the way.


```python
for t in range(1000):
    fdir.update()

    tl_diff.set_value("topographic__elevation", fdir.get_value("topographic__elevation"))
    tl_diff.set_value("flow__receiver_node", fdir.get_value("flow__receiver_node"))
    tl_diff.set_value("topographic__steepest_slope", fdir.get_value("topographic__steepest_slope"))

    tl_diff.update()
    
    fdir.set_value("topographic__elevation", tl_diff.get_value("topographic__elevation"))

    # add some output to let us see we aren't hanging:
    if t % 100 == 0:
        print(t * dt)
        
        z = tl_diff.get_value("topographic__elevation").reshape((100, 100))
        
        # plot east-west cross-section of topography:
        x_plot = range(0, 1000, 10)
        z_plot = z[1, :]
        figure('cross-section')
        plot(x_plot, z_plot)
```


```python
fdir.quick_plot("topographic__elevation")
```

The material is diffused from the top and along the slope and it accumulates at the bottom, where the topography flattens.

# Example 3

As a comparison, the following code uses linear diffusion on the same slope. Instead of using the **TransportLengthHillslopeDiffuser** component, we'll swap in the **LinearDiffuser** component. Everything else will be pretty much the same.


```python
fdir = pymt.models.FlowDirectorSteepest()
diff = pymt.models.LinearDiffuser()
```

Setup and initialize the models.


```python
config_file, config_dir = fdir.setup(**grid_params, **clock_params)
fdir.initialize(config_file, config_dir)
```


```python
config_file, config_dir = diff.setup(
    **grid_params,
    **clock_params,
    linear_diffusivity=0.1,
)
diff.initialize(config_file, config_dir)
```

Set boundary conditions.


```python
status = fdir.get_value("boundary_condition_flag").reshape((100, 100))
status[:, (0, -1)] = 1 # E and W boundaries are open 
status[(0, -1), :] = 4 # N and S boundaries are closed
```


```python
fdir.set_value("boundary_condition_flag", status)
diff.set_value("boundary_condition_flag", status)
```

Set the initial topography.


```python
grid = fdir.var_grid("topographic__elevation")
n_vals = fdir.grid_node_count(grid)
x, y = fdir.grid[0].node_x, fdir.grid[0].node_y
```


```python
z = np.zeros(n_vals)
z[x > 500] = x[x < 490] / 10.0
```


```python
fdir.set_value("topographic__elevation", z)
fdir.quick_plot("topographic__elevation")
```

Run the model!


```python
for t in range(1000):
    fdir.update()

    diff.set_value("topographic__elevation", fdir.get_value("topographic__elevation"))

    diff.update()
    
    fdir.set_value("topographic__elevation", diff.get_value("topographic__elevation"))

    # add some output to let us see we aren't hanging:
    if t % 100 == 0:
        print(t * dt)
        
        z = diff.get_value("topographic__elevation").reshape((100, 100))
        
        # plot east-west cross-section of topography:
        x_plot = range(0, 1000, 10)
        z_plot = z[1, :]
        figure('cross-section')
        plot(x_plot, z_plot)
```


```python
fdir.quick_plot("topographic__elevation")
```


```python

```

129b8f49bcaa8f0d436d29618cc654ce21a73c76
ipynb
Jupyter Notebook
nb/transport_length_hillslope_diffuser.ipynb
mcflugen/pymt-live
dcfc430a4953ec33bcf3c6efce0fe4b35d3358a9
MIT
nb/transport_length_hillslope_diffuser.ipynb
mcflugen/pymt-live
dcfc430a4953ec33bcf3c6efce0fe4b35d3358a9
MIT
2019-05-22T20:47:43.000Z
2019-12-30T14:07:12.000Z
nb/transport_length_hillslope_diffuser.ipynb
csdms/pymt-live
dcfc430a4953ec33bcf3c6efce0fe4b35d3358a9
MIT
Qwen/Qwen-72B
1. YES
2. YES