{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The transport-length hillslope diffuser" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "For more Landlab tutorials, click here: https://landlab.readthedocs.io/en/latest/user_guide/tutorials.html\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This Jupyter notebook illustrates running the transport-length-model hillslope diffusion component in a simple example.\n", "\n", "# The Basics\n", "\n", "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.\n", "\n", "## Theory\n", "\n", "The elevation $z$ of a point of the landscape (such as a grid node) changes according to:\n", "\n", "\\begin{equation}\n", " \\frac{\\partial z}{\\partial t} = -\\epsilon + D + U \\tag{1}\\label{eq:1},\n", "\\end{equation}\n", "\n", "\n", "and we define:\n", "\\begin{equation}\n", " D = \\frac{q_s}{L} \\tag{2}\\label{eq:2},\n", "\\end{equation}\n", "\n", "\n", "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**.\n", "\n", "We specify the erosion rate $\\epsilon$ and the transport length $L$:\n", "\n", "\\begin{equation}\n", " \\epsilon = \\kappa S \\tag{3}\\label{eq:3}\n", "\\end{equation}\n", "\n", "\n", "\n", "\\begin{equation}\n", " L = \\frac{dx}{1-({S}/{S_c})^2} \\tag{4}\\label{eq:4}\n", "\\end{equation}\n", "\n", "where $\\kappa$ [*L/T*] is an erodibility coefficient, $S$ is the local slope [*L/L*] and $S_c$ is the critical slope [*L/L*]. \n", "\n", "Thus, the elevation variation results from the difference between local rates of detachment and deposition. \n", "\n", "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:\n", "- 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\".\n", "- 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.\n", "- 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).\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Contrast with the non-linear diffusion model\n", "\n", "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:\n", "\n", "$ $\n", "\n", "\\begin{equation} \n", " \\frac{\\partial z}{\\partial t} = \\frac{\\partial q_s}{\\partial x} \\tag{5}\\label{eq:5}\n", "\\end{equation}\n", "\n", "$ $\n", "\\begin{equation}\n", " q_s = \\frac{\\kappa' S}{1-({S}/{S_c})^2} \\tag{6}\\label{eq:6}\n", "\\end{equation}\n", "\n", "where $\\kappa'$ [*L$^2$/T*] is a diffusion coefficient.\n", "\n", "This description is thus based on the definition of a flux of transported sediment parallel to the slope:\n", "- when the slope is small, this flux refers to diffusion-like processes such as biogenic soil disturbance, rain splash, or diffuse runoff\n", "- when the slope gets closer to the specified critical slope, the flux increases dramatically, simulating on average the cumulative effect of mass wasting events.\n", "\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Example 1:\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "First, we import what we'll need:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib.pyplot import figure, plot, show, title, xlabel, ylabel\n", "\n", "from landlab import RasterModelGrid\n", "from landlab.components import FlowDirectorSteepest, TransportLengthHillslopeDiffuser\n", "from landlab.plot import imshow_grid\n", "\n", "# to plot figures in the notebook:\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a grid and set boundary conditions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mg = RasterModelGrid(\n", " (20, 20),\n", " xy_spacing=50.0) # raster grid with 20 rows, 20 columns and dx=50m\n", "z = np.random.rand(mg.size('node')) # random noise for initial topography\n", "mg.add_field(\"topographic__elevation\", z, at=\"node\")\n", "\n", "mg.set_closed_boundaries_at_grid_edges(\n", " False, True, False,\n", " True) # N and S boundaries are closed, E and W are open" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the initial and run conditions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "total_t = 2000000.0 # total run time (yr)\n", "dt = 1000.0 # time step (yr)\n", "nt = int(total_t // dt) # number of time steps\n", "uplift_rate = 0.0001 # uplift rate (m/yr)\n", "\n", "kappa = 0.001 # erodibility (m/yr)\n", "Sc = 0.6 # critical slope" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instantiate the components:\n", "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)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fdir = FlowDirectorSteepest(mg)\n", "tl_diff = TransportLengthHillslopeDiffuser(mg,\n", " erodibility=kappa,\n", " slope_crit=Sc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the components for 2 Myr and trace an East-West cross-section of the topography every 100 kyr:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for t in range(nt):\n", " fdir.run_one_step()\n", " tl_diff.run_one_step(dt)\n", " z[mg.core_nodes] += uplift_rate * dt # add the uplift\n", "\n", " # add some output to let us see we aren't hanging:\n", " if t % 100 == 0:\n", " print(t * dt)\n", "\n", " # plot east-west cross-section of topography:\n", " x_plot = range(0, 1000, 50)\n", " z_plot = z[100:120]\n", " figure(\"cross-section\")\n", " plot(x_plot, z_plot)\n", "\n", "figure(\"cross-section\")\n", "title(\"East-West cross section\")\n", "xlabel(\"x (m)\")\n", "ylabel(\"z (m)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And plot final topography:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "figure(\"final topography\")\n", "im = imshow_grid(mg,\n", " \"topographic__elevation\",\n", " grid_units=[\"m\", \"m\"],\n", " var_name=\"Elevation (m)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This behaviour corresponds to the evolution observed using a classical non-linear diffusion model." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Example 2: \n", "\n", "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.\n", "\n", "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.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create grid and topographic elevation field:\n", "mg2 = RasterModelGrid((20, 20), xy_spacing=50.0)\n", "\n", "z = np.zeros(mg2.number_of_nodes)\n", "z[mg2.node_x > 500] = mg2.node_x[mg2.node_x > 500] / 10\n", "mg2.add_field(\"topographic__elevation\", z, at=\"node\")\n", "\n", "# Set boundary conditions:\n", "mg2.set_closed_boundaries_at_grid_edges(False, True, False, True)\n", "\n", "# Show initial topography:\n", "im = imshow_grid(mg2,\n", " \"topographic__elevation\",\n", " grid_units=[\"m\", \"m\"],\n", " var_name=\"Elevation (m)\")\n", "\n", "# Plot an east-west cross-section of the initial topography:\n", "z_plot = z[100:120]\n", "x_plot = range(0, 1000, 50)\n", "figure(2)\n", "plot(x_plot, z_plot)\n", "title(\"East-West cross section\")\n", "xlabel(\"x (m)\")\n", "ylabel(\"z (m)\")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Set the run conditions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "total_t = 1000000.0 # total run time (yr)\n", "dt = 1000.0 # time step (yr)\n", "nt = int(total_t // dt) # number of time steps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instantiate the components:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fdir = FlowDirectorSteepest(mg2)\n", "tl_diff = TransportLengthHillslopeDiffuser(mg2,\n", " erodibility=0.001,\n", " slope_crit=0.6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run for 1 Myr, plotting the cross-section regularly:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for t in range(nt):\n", " fdir.run_one_step()\n", " tl_diff.run_one_step(dt)\n", "\n", " # add some output to let us see we aren't hanging:\n", " if t % 100 == 0:\n", " print(t * dt)\n", " z_plot = z[100:120]\n", " figure(2)\n", " plot(x_plot, z_plot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The material is diffused from the top and along the slope and it accumulates at the bottom, where the topography flattens.\n", "\n", "As a comparison, the following code uses linear diffusion on the same slope:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import Linear diffuser:\n", "from landlab.components import LinearDiffuser\n", "\n", "# Create grid and topographic elevation field:\n", "mg3 = RasterModelGrid((20, 20), xy_spacing=50.0)\n", "z = np.ones(mg3.number_of_nodes)\n", "z[mg.node_x > 500] = mg.node_x[mg.node_x > 500] / 10\n", "mg3.add_field(\"topographic__elevation\", z, at=\"node\")\n", "\n", "# Set boundary conditions:\n", "mg3.set_closed_boundaries_at_grid_edges(False, True, False, True)\n", "\n", "# Instantiate components:\n", "fdir = FlowDirectorSteepest(mg3)\n", "diff = LinearDiffuser(mg3, linear_diffusivity=0.1)\n", "\n", "# Set run conditions:\n", "total_t = 1000000.0\n", "dt = 1000.0\n", "nt = int(total_t // dt)\n", "\n", "# Run for 1 Myr, plotting east-west cross-section regularly:\n", "for t in range(nt):\n", " fdir.run_one_step()\n", " diff.run_one_step(dt)\n", "\n", " # add some output to let us see we aren't hanging:\n", " if t % 100 == 0:\n", " print(t * dt)\n", " z_plot = z[100:120]\n", " figure(2)\n", " plot(x_plot, z_plot)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.1" } }, "nbformat": 4, "nbformat_minor": 2 }