<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <a href="http
Step1:  Make a grid and set boundary conditions
Step2:  Set the initial and run conditions
Step3:  Instantiate the components
Step4:  Run the components for 2 Myr and trace an East-West cross-section of the topography every 100 kyr
Step5:  And plot final topography
Step6:  This behaviour corresponds to the evolution observed using a classical non-linear diffusion model.
Example 2
Step7:  Set the run conditions
Step8:  Instantiate the components
Step9:  Run for 1 Myr, plotting the cross-section regularly
Step10:  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
<ASSISTANT_TASK:>
Python Code:
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
Explanation: <a href="http://landlab.github.io"><img style="float: left" src="../../landlab_header.png"></a>
The transport-length hillslope diffuser
<hr>
<small>For more Landlab tutorials, click here: <a href="https://landlab.readthedocs.io/en/latest/user_guide/tutorials.html">https://landlab.readthedocs.io/en/latest/user_guide/tutorials.html</a></small>
<hr>
This Jupyter notebook illustrates running the transport-length-model hillslope diffusion component in a simple example.
The Basics
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.
Theory
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:
- 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-like processes such as biogenic soil disturbance, rain splash, or diffuse runoff
- when the slope gets closer to the specified critical slope, the flux increases dramatically, simulating on average the cumulative effect of mass wasting events.
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.
Example 1:
First, we import what we'll need:
End of explanation
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
Explanation: Make a grid and set boundary conditions:
End of explanation
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
Explanation: Set the initial and run conditions:
End of explanation
fdir = FlowDirectorSteepest(mg)
tl_diff = TransportLengthHillslopeDiffuser(mg,
                                           erodibility=kappa,
                                           slope_crit=Sc)
Explanation: 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).
End of explanation
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)")
Explanation: Run the components for 2 Myr and trace an East-West cross-section of the topography every 100 kyr:
End of explanation
figure("final topography")
im = imshow_grid(mg,
                 "topographic__elevation",
                 grid_units=["m", "m"],
                 var_name="Elevation (m)")
Explanation: And plot final topography:
End of explanation
# 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)")
Explanation: This behaviour corresponds to the evolution observed using a classical non-linear diffusion model.
Example 2:
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.
End of explanation
total_t = 1000000.0  # total run time (yr)
dt = 1000.0  # time step (yr)
nt = int(total_t // dt)  # number of time steps
Explanation: Set the run conditions:
End of explanation
fdir = FlowDirectorSteepest(mg2)
tl_diff = TransportLengthHillslopeDiffuser(mg2,
                                           erodibility=0.001,
                                           slope_crit=0.6)
Explanation: Instantiate the components:
End of explanation
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)
Explanation: Run for 1 Myr, plotting the cross-section regularly:
End of explanation
# 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)
Explanation: 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:
End of explanation
<END_TASK>