diff --git a/lessons/landlab/landlab/fault-scarp.ipynb b/lessons/landlab/landlab/fault-scarp.ipynb
deleted file mode 100644
index 43c4257..0000000
--- a/lessons/landlab/landlab/fault-scarp.ipynb
+++ /dev/null
@@ -1,1334 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- ""
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Introduction to Landlab: Creating a simple 2D scarp diffusion model"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "
\n",
- "For more Landlab tutorials, click here: https://landlab.readthedocs.io/en/latest/user_guide/tutorials.html\n",
- "
\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "This tutorial illustrates how you can use Landlab to construct a simple two-dimensional numerical model on a regular (raster) grid, using a simple forward-time, centered-space numerical scheme. The example is the erosional degradation of an earthquake fault scarp, and which evolves over time in response to the gradual downhill motion of soil. Here we use a simple \"geomorphic diffusion\" model for landform evolution, in which the downhill flow of soil is assumed to be proportional to the (downhill) gradient of the land surface multiplied by a transport coefficient.\n",
- "\n",
- "We start by importing the [numpy](https://numpy.org) and [matplotlib](https://matplotlib.org) libraries:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "import matplotlib.pyplot as plt\n",
- "import numpy as np"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Part 1: 1D version using numpy\n",
- "\n",
- "This example uses a finite-volume numerical solution to the 2D diffusion equation. The 2D diffusion equation in this case is derived as follows. Continuity of mass states that:\n",
- "\n",
- "$\\frac{\\partial z}{\\partial t} = -\\nabla \\cdot \\mathbf{q}_s$,\n",
- "\n",
- "where $z$ is elevation, $t$ is time, the vector $\\mathbf{q}_s$ is the volumetric soil transport rate per unit width, and $\\nabla$ is the divergence operator (here in two dimensions). (Note that we have omitted a porosity factor here; its effect will be subsumed in the transport coefficient). The sediment flux vector depends on the slope gradient:\n",
- "\n",
- "$\\mathbf{q}_s = -D \\nabla z$,\n",
- "\n",
- "where $D$ is a transport-rate coefficient---sometimes called *hillslope diffusivity*---with dimensions of length squared per time. Combining the two, and assuming $D$ is uniform, we have a classical 2D diffusion equation:\n",
- "\n",
- "$\\frac{\\partial z}{\\partial t} = -\\nabla^2 z$.\n",
- "\n",
- "In this first example, we will create a our 1D domain in $x$ and $z$, and set a value for $D$.\n",
- "\n",
- "This means that the equation we solve will be in 1D. \n",
- "\n",
- "$\\frac{d z}{d t} = \\frac{d q_s}{dx}$,\n",
- "\n",
- "where \n",
- "\n",
- "$q_s = -D \\frac{d z}{dx}$\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "dx = 1\n",
- "x = np.arange(0, 100, dx, dtype=float)\n",
- "z = np.zeros(x.shape, dtype=float)\n",
- "D = 0.01"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Next we must create our fault by uplifting some of the domain. We will increment all elements of `z` in which `x>50`."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "z[x > 50] += 100"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Finally, we will diffuse our fault for 1,000 years.\n",
- "\n",
- "We will use a timestep with a [Courant–Friedrichs–Lewy condition](https://en.wikipedia.org/wiki/Courant–Friedrichs–Lewy_condition) of $C_{cfl}=0.2$. This will keep our solution numerically stable. \n",
- "\n",
- "$C_{cfl} = \\frac{\\Delta t D}{\\Delta x^2} = 0.2$"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "dt = 0.2 * dx * dx / D\n",
- "total_time = 1e3\n",
- "nts = int(total_time / dt)\n",
- "z_orig = z.copy()\n",
- "for i in range(nts):\n",
- " qs = -D * np.diff(z) / dx\n",
- " dzdt = -np.diff(qs) / dx\n",
- " z[1:-1] += dzdt * dt\n",
- "\n",
- "plt.plot(x, z_orig, label=\"Original Profile\")\n",
- "plt.plot(x, z, label=\"Diffused Profile\")\n",
- "plt.legend()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The prior example is pretty simple. If this was all you needed to do, you wouldn't need Landlab. \n",
- "\n",
- "But what if you wanted...\n",
- "\n",
- "... to use the same diffusion model in 2D instead of 1D.\n",
- "\n",
- "... to use an irregular grid (in 1 or 2D). \n",
- "\n",
- "... wanted to combine the diffusion model with a more complex model. \n",
- "\n",
- "... have a more complex model you want to use over and over again with different boundary conditions.\n",
- "\n",
- "These are the sorts of problems that Landlab was designed to solve. \n",
- "\n",
- "In the next two sections we will introduce some of the core capabilities of Landlab. \n",
- "\n",
- "In Part 2 we will use the RasterModelGrid, fields, and a numerical utility for calculating flux divergence. \n",
- "\n",
- "In Part 3 we will use the HexagonalModelGrid. \n",
- "\n",
- "In Part 4 we will use the LinearDiffuser component. "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Part 2: 2D version using Landlab's Model Grids\n",
- "\n",
- "The Landlab model grids are data structures that represent the model domain (the variable `x` in our prior example). Here we will use `RasterModelGrid` which creates a grid with regularly spaced square grid elements. The RasterModelGrid knows how the elements are connected and how far apart they are.\n",
- "\n",
- "Lets start by creating a RasterModelGrid class. First we need to import it. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from landlab import RasterModelGrid"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "### (a) Explore the RasterModelGrid\n",
- "\n",
- "Before we make a RasterModelGrid for our fault example, lets explore the Landlab model grid. \n",
- "\n",
- "Landlab considers the grid as a \"dual\" graph. Two sets of points, lines and polygons that represent 2D space. \n",
- "\n",
- "The first graph considers points called \"nodes\" that are connected by lines called \"links\". The area that surrounds each node is called a \"cell\".\n",
- "\n",
- "First, the nodes"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from landlab.plot.graph import plot_graph\n",
- "\n",
- "grid = RasterModelGrid((4, 5), xy_spacing=(3, 4))\n",
- "plot_graph(grid, at=\"node\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "You can see that the nodes are points and they are numbered with unique IDs from lower left to upper right. \n",
- "\n",
- "Next the links"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "plot_graph(grid, at=\"link\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which are lines that connect the nodes and each have a unique ID number. \n",
- "\n",
- "And finally, the cells"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "plot_graph(grid, at=\"cell\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which are polygons centered around the nodes. \n",
- "\n",
- "Landlab is a \"dual\" graph because it also keeps track of a second set of points, lines, and polygons (\"corners\", \"faces\", and \"patches\"). We will not focus on them further."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Exercises for section 2a\n",
- "\n",
- "(2a.1) Create an instance of a `RasterModelGrid` with 5 rows and 7 columns, with a spacing between nodes of 10 units. Plot the node layout, and identify the ID number of the center-most node."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2a.1 here)\n",
- "\n",
- "rmg = RasterModelGrid((5, 7), 10.0)\n",
- "plot_graph(rmg, at=\"node\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2a.2) Find the ID of the cell that contains this node."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2a.2 here)\n",
- "plot_graph(rmg, at=\"cell\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2a.3) Find the ID of the horizontal link that connects to the last node on the right in the middle column."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2a.3 here)\n",
- "plot_graph(rmg, at=\"link\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### (b) Use the RasterModelGrid for 2D diffusion \n",
- "\n",
- "Lets continue by making a new grid that is bigger. We will use this for our next fault diffusion example.\n",
- "\n",
- "The syntax in the next line says: create a new *RasterModelGrid* object called **mg**, with 25 rows, 40 columns, and a grid spacing of 10 m."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mg = RasterModelGrid((25, 40), 10.0)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Note the use of object-oriented programming here. `RasterModelGrid` is a class; `mg` is a particular instance of that class, and it contains all the data necessary to fully describe the topology and geometry of this particular grid.\n",
- "\n",
- "Next we'll add a *data field* to the grid, to represent the elevation values at grid nodes. The \"dot\" syntax below indicates that we are calling a function (or *method*) that belongs to the *RasterModelGrid* class, and will act on data contained in **mg**. The arguments indicate that we want the data elements attached to grid nodes (rather than links, for example), and that we want to name this data field `topographic__elevation`. The `add_zeros` method returns the newly created NumPy array."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "z = mg.add_zeros(\"topographic__elevation\", at=\"node\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The above line of code creates space in memory to store 1,000 floating-point values, which will represent the elevation of the land surface at each of our 1,000 grid nodes."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Let's plot the positions of all the grid nodes. The nodes' *(x,y)* positions are stored in the arrays `mg.x_of_node` and `mg.y_of_node`, respectively."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "plt.plot(mg.x_of_node, mg.y_of_node, \".\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "If we bothered to count, we'd see that there are indeed 1,000 grid nodes, and a corresponding number of `z` values:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "len(z)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now for some tectonics. Let's say there's a fault trace that angles roughly east-northeast. We can describe the trace with the equation for a line. One trick here: by using `mg.x_of_node`, in the line of code below, we are calculating a *y* (i.e., north-south) position of the fault trace for each grid node---meaning that this is the *y* coordinate of the trace at the *x* coordinate of a given node."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "fault_trace_y = 50.0 + 0.25 * mg.x_of_node"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Here comes the earthquake. For all the nodes north of the fault (i.e., those with a *y* coordinate greater than the corresponding *y* coordinate of the fault trace), we'll add elevation equal to 10 meters plus a centimeter for every meter east along the grid (just to make it interesting):"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "z[mg.y_of_node > fault_trace_y] += (\n",
- " 10.0 + 0.01 * mg.x_of_node[mg.y_of_node > fault_trace_y]\n",
- ")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(A little bit of Python under the hood: the statement `mg.y_of_node > fault_trace_y` creates a 1000-element long boolean array; placing this within the index brackets will select only those array entries that correspond to `True` in the boolean array)\n",
- "\n",
- "Let's look at our newly created initial topography using Landlab's *imshow_node_grid* plotting function (which we first need to import)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from landlab.plot.imshow import imshow_grid\n",
- "\n",
- "imshow_grid(mg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "To finish getting set up, we will define two parameters: the transport (\"diffusivity\") coefficient, `D`, and the time-step size, `dt`. (The latter is set using the Courant condition for a forward-time, centered-space finite-difference solution; you can find the explanation in most textbooks on numerical methods)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "D = 0.01 # m2/yr transport coefficient\n",
- "dt = 0.2 * mg.dx * mg.dx / D\n",
- "dt"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Boundary conditions: for this example, we'll assume that the east and west sides are closed to flow of sediment, but that the north and south sides are open. (The order of the function arguments is east, north, west, south)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mg.set_closed_boundaries_at_grid_edges(True, False, True, False)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "*A note on boundaries:* with a Landlab raster grid, all the perimeter nodes are boundary nodes. In this example, there are 24 + 24 + 39 + 39 = 126 boundary nodes. The previous line of code set those on the east and west edges to be **closed boundaries**, while those on the north and south are **open boundaries** (the default). All the remaining nodes are known as **core** nodes. In this example, there are 1000 - 126 = 874 core nodes:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "len(mg.core_nodes)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "One more thing before we run the time loop: we'll create an array to contain soil flux. In the function call below, the first argument tells Landlab that we want one value for each grid link, while the second argument provides a name for this data *field*:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "qs = mg.add_zeros(\"sediment_flux\", at=\"link\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "And now for some landform evolution. We will loop through 25 iterations, representing 50,000 years. On each pass through the loop, we do the following:\n",
- "\n",
- "1. Calculate, and store in the array `g`, the gradient between each neighboring pair of nodes. These calculations are done on **links**. The gradient value is a positive number when the gradient is \"uphill\" in the direction of the link, and negative when the gradient is \"downhill\" in the direction of the link. On a raster grid, link directions are always in the direction of increasing $x$ (\"horizontal\" links) or increasing $y$ (\"vertical\" links).\n",
- "\n",
- "2. Calculate, and store in the array `qs`, the sediment flux between each adjacent pair of nodes by multiplying their gradient by the transport coefficient. We will only do this for the **active links** (those not connected to a closed boundary, and not connecting two boundary nodes of any type); others will remain as zero.\n",
- "\n",
- "3. Calculate the resulting net flux at each node (positive=net outflux, negative=net influx). The negative of this array is the rate of change of elevation at each (core) node, so store it in a node array called `dzdt'.\n",
- "\n",
- "4. Update the elevations for the new time step."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "for i in range(25):\n",
- " g = mg.calc_grad_at_link(z)\n",
- " qs[mg.active_links] = -D * g[mg.active_links]\n",
- " dzdt = -mg.calc_flux_div_at_node(qs)\n",
- " z[mg.core_nodes] += dzdt[mg.core_nodes] * dt"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Let's look at how our fault scarp has evolved."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "imshow_grid(mg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Notice that we have just created and run a 2D model of fault-scarp creation and diffusion with fewer than two dozen lines of code. How long would this have taken to write in C or Fortran?\n",
- "\n",
- "While it was very very easy to write in 1D, writing this in 2D would mean we would have needed to keep track of the adjacency of the different parts of the grid. This is the primary problem that the Landlab grids are meant to solve. \n",
- "\n",
- "Think about how difficult this would be to hand code if the grid were irregular or hexagonal. In order to conserve mass and implement the differential equation you would need to know how nodes were conected, how long the links were, and how big each cell was.\n",
- "\n",
- "We do such an example after the next section. "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Exercises for section 2b\n",
- "\n",
- "(2b .1) Create an instance of a `RasterModelGrid` called `mygrid`, with 16 rows and 25 columns, with a spacing between nodes of 5 meters. Use the `plot` function in the `matplotlib` library to make a plot that shows the position of each node marked with a dot (hint: see the plt.plot() example above)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.1 here)\n",
- "mygrid = RasterModelGrid((16, 25), xy_spacing=5.0)\n",
- "plt.plot(mygrid.x_of_node, mygrid.y_of_node, \".\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2b.2) Query the grid variables `number_of_nodes` and `number_of_core_nodes` to find out how many nodes are in your grid, and how many of them are core nodes."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.2 here)\n",
- "print(mygrid.number_of_nodes)\n",
- "print(mygrid.number_of_core_nodes)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2b.3) Add a new field to your grid, called `temperature` and attached to nodes. Have the initial values be all zero."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.3 here)\n",
- "temp = mygrid.add_zeros(\"temperature\", at=\"node\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2b.4) Change the temperature of nodes in the top (north) half of the grid to be 10 degrees C. Use the `imshow_grid` function to display a shaded image of the elevation field."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.4 here)\n",
- "temp[mygrid.y_of_node >= 40.0] = 10.0\n",
- "imshow_grid(mygrid, \"temperature\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2b.5) Use the grid function `set_closed_boundaries_at_grid_edges` to assign closed boundaries to the right and left sides of the grid."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.5 here)\n",
- "mygrid.set_closed_boundaries_at_grid_edges(True, False, True, False)\n",
- "imshow_grid(mygrid, \"temperature\", color_for_closed=\"c\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2b.6) Create a new field of zeros called `heat_flux` and attached to links. Using the `number_of_links` grid variable, verify that your new field array has the correct number of items. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.6 here)\n",
- "Q = mygrid.add_zeros(\"heat_flux\", at=\"link\")\n",
- "print(mygrid.number_of_links)\n",
- "print(len(Q))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2b.7) Use the `calc_grad_at_link` grid function to calculate the temperature gradients at all the links in the grid. Given the node spacing and the temperatures you assigned to the top versus bottom grid nodes, what do you expect the maximum temperature gradient to be? Print the values in the gradient array to verify that this is indeed the maximum temperature gradient."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.7 here)\n",
- "print(\"Expected max gradient is 2 C/m\")\n",
- "temp_grad = mygrid.calc_grad_at_link(temp)\n",
- "print(temp_grad)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2b.8) Back to hillslopes: Reset the values in the elevation field of the grid `mg` to zero. Then copy and paste the time loop above (i.e., the block in Section 2b that starts with `for i in range(25):`) below. Modify the last line to add uplift of the hillslope material at a rate `uplift_rate` = 0.0001 m/yr (hint: the amount of uplift in each iteration should be the uplift rate times the time-step duration). Then run the block and plot the resulting topography. Try experimenting with different uplift rates and different values of `D`."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2b.8 here)\n",
- "z[:] = 0.0\n",
- "uplift_rate = 0.0001\n",
- "for i in range(25):\n",
- " g = mg.calc_grad_at_link(z)\n",
- " qs[mg.active_links] = -D * g[mg.active_links]\n",
- " dzdt = -mg.calc_flux_div_at_node(qs)\n",
- " z[mg.core_nodes] += (dzdt[mg.core_nodes] + uplift_rate) * dt\n",
- "imshow_grid(mg, z)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### (c) What's going on under the hood?\n",
- "\n",
- "This example uses a finite-volume numerical solution to the 2D diffusion equation. The 2D diffusion equation in this case is derived as follows. Continuity of mass states that:\n",
- "\n",
- "$\\frac{\\partial z}{\\partial t} = -\\nabla \\cdot \\mathbf{q}_s$,\n",
- "\n",
- "where $z$ is elevation, $t$ is time, the vector $\\mathbf{q}_s$ is the volumetric soil transport rate per unit width, and $\\nabla$ is the divergence operator (here in two dimensions). (Note that we have omitted a porosity factor here; its effect will be subsumed in the transport coefficient). The sediment flux vector depends on the slope gradient:\n",
- "\n",
- "$\\mathbf{q}_s = -D \\nabla z$,\n",
- "\n",
- "where $D$ is a transport-rate coefficient---sometimes called *hillslope diffusivity*---with dimensions of length squared per time. Combining the two, and assuming $D$ is uniform, we have a classical 2D diffusion equation:\n",
- "\n",
- "$\\frac{\\partial z}{\\partial t} = -\\nabla^2 z$.\n",
- "\n",
- "For the numerical solution, we discretize $z$ at a series of *nodes* on a grid. The example in this notebook uses a Landlab *RasterModelGrid*, in which every interior node sits inside a cell of width $\\Delta x$, but we could alternatively have used any grid type that provides nodes, links, and cells.\n",
- "\n",
- "The gradient and sediment flux vectors will be calculated at the *links* that connect each pair of adjacent nodes. These links correspond to the mid-points of the cell faces, and the values that we assign to links represent the gradients and fluxes, respectively, along the faces of the cells.\n",
- "\n",
- "The flux divergence, $\\nabla \\mathbf{q}_s$, will be calculated by summing, for every cell, the total volume inflows and outflows at each cell face, and dividing the resulting sum by the cell area. Note that for a regular, rectilinear grid, as we use in this example, this finite-volume method is equivalent to a finite-difference method.\n",
- "\n",
- "To advance the solution in time, we will use a simple explicit, forward-difference method. This solution scheme for a given node $i$ can be written:\n",
- "\n",
- "$\\frac{z_i^{t+1} - z_i^t}{\\Delta t} = -\\frac{1}{A_i} \\sum\\limits_{j=1}^{N_i} \\delta (l_{ij}) q_s (l_{ij}) \\lambda(l_{ij})$.\n",
- "\n",
- "Here the superscripts refer to time steps, $\\Delta t$ is time-step size, $q_s(l_{ij})$ is the sediment flux per width associated with the link that crosses the $j$-th face of the cell at node $i$, $\\lambda(l_{ij})$ is the width of the cell face associated with that link ($=\\Delta x$ for a regular uniform grid), and $N_i$ is the number of active links that connect to node $i$. The variable $\\delta(l_{ij})$ contains either +1 or -1: it is +1 if link $l_{ij}$ is oriented away from the node (in which case positive flux would represent material leaving its cell), or -1 if instead the link \"points\" into the cell (in which case positive flux means material is entering).\n",
- "\n",
- "To get the fluxes, we first calculate the *gradient*, $G$, at each link, $k$:\n",
- "\n",
- "$G(k) = \\frac{z(H_k) - z(T_k)}{L_k}$.\n",
- "\n",
- "Here $H_k$ refers the *head node* associated with link $k$, $T_k$ is the *tail node* associated with link $k$. Each link has a direction: from the tail node to the head node. The length of link $k$ is $L_k$ (equal to $\\Delta x$ is a regular uniform grid). What the above equation says is that the gradient in $z$ associated with each link is simply the difference in $z$ value between its two endpoint nodes, divided by the distance between them. The gradient is positive when the value at the head node (the \"tip\" of the link) is greater than the value at the tail node, and vice versa.\n",
- "\n",
- "The calculation of gradients in $z$ at the links is accomplished with the `calc_grad_at_link` function. The sediment fluxes are then calculated by multiplying the link gradients by $-D$. Once the fluxes at links have been established, the `calc_flux_div_at_node` function performs the summation of fluxes."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Exercises for section 2c\n",
- "\n",
- "(2c.1) Make a 3x3 `RasterModelGrid` called `tinygrid`, with a cell spacing of 2 m. Use the `plot_graph` function to display the nodes and their ID numbers."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2c.1 here)\n",
- "tinygrid = RasterModelGrid((3, 3), 2.0)\n",
- "plot_graph(tinygrid, at=\"node\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2c.2) Give your `tinygrid` a node field called `height` and set the height of the center-most node to 0.5. Use `imshow_grid` to display the height field."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2c.2 here)\n",
- "ht = tinygrid.add_zeros(\"height\", at=\"node\")\n",
- "ht[4] = 0.5\n",
- "imshow_grid(tinygrid, ht)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2c.3) The grid should have 12 links (extra credit: verify this with `plot_graph`). When you compute gradients, which of these links will have non-zero gradients? What will the absolute value(s) of these gradients be? Which (if any) will have positive gradients and which negative? To codify your answers, make a 12-element numpy array that contains your predicted gradient value for each link."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2c.3 here)\n",
- "plot_graph(tinygrid, at=\"link\")\n",
- "pred_grad = np.array([0, 0, 0, 0.25, 0, 0.25, -0.25, 0, -0.25, 0, 0, 0])\n",
- "print(pred_grad)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2c.4) Test your prediction by running the `calc_grad_at_link` function on your tiny grid. Print the resulting array and compare it with your predictions."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2c.4 here)\n",
- "grad = tinygrid.calc_grad_at_link(ht)\n",
- "print(grad)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2c.5) Suppose the flux of soil per unit cell width is defined as -0.01 times the height gradient. What would the flux be at the those links that have non-zero gradients? Test your prediction by creating and printing a new array whose values are equal to -0.01 times the link-gradient values."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2c.5 here)\n",
- "flux = -0.01 * grad\n",
- "print(flux)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(2c.6) Consider the net soil accumulation or loss rate around the center-most node in your tiny grid (which is the only one that has a cell). The *divergence* of soil flux can be represented numerically as the sum of the total volumetric soil flux across each of the cell's four faces. What is the flux across each face? (Hint: multiply by face width) What do they add up to? Test your prediction by running the grid function `calc_flux_div_at_node` (hint: pass your unit flux array as the argument). What are the units of the divergence values returned by the `calc_flux_div_at_node` function?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 2c.6 here)\n",
- "print(\"predicted div is 0 m/yr\")\n",
- "dqsdx = tinygrid.calc_flux_div_at_node(flux)\n",
- "print(dqsdx)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Part 3: Hexagonal grid\n",
- "\n",
- "Next we will use an non-raster Landlab grid.\n",
- "\n",
- "We start by making a random set of points with x values between 0 and 400 and y values of 0 and 250. We then add zeros to our grid at a field called \"topographic__elevation\" and plot the node locations. \n",
- "\n",
- "Note that the syntax here is exactly the same as in the RasterModelGrid example (once the grid has been created)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from landlab import HexModelGrid\n",
- "\n",
- "mg = HexModelGrid((25, 40), 10, node_layout=\"rect\")\n",
- "z = mg.add_zeros(\"topographic__elevation\", at=\"node\")\n",
- "plt.plot(mg.x_of_node, mg.y_of_node, \".\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Next we create our fault trace and uplift the hanging wall. \n",
- "\n",
- "We can plot just like we did with the RasterModelGrid. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "fault_trace_y = 50.0 + 0.25 * mg.x_of_node\n",
- "z[mg.y_of_node > fault_trace_y] += (\n",
- " 10.0 + 0.01 * mg.x_of_node[mg.y_of_node > fault_trace_y]\n",
- ")\n",
- "imshow_grid(mg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "And we can use the same code as before to create a diffusion model!\n",
- "\n",
- "Landlab supports multiple grid types. You can read more about them [here](https://landlab.readthedocs.io/en/latest/reference/grid/index.html)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "qs = mg.add_zeros(\"sediment_flux\", at=\"link\")\n",
- "for i in range(25):\n",
- " g = mg.calc_grad_at_link(z)\n",
- " qs[mg.active_links] = -D * g[mg.active_links]\n",
- " dzdt = -mg.calc_flux_div_at_node(qs)\n",
- " z[mg.core_nodes] += dzdt[mg.core_nodes] * dt\n",
- "imshow_grid(mg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### Exercises for section 3\n",
- "\n",
- "(3.1-6) Repeat the exercises from section 2c, but this time using a hexagonal tiny grid called `tinyhex`. Your grid should have 7 nodes: one core node and 6 perimeter nodes. (Hints: use `node_layout = 'hex'`, and make a grid with 3 rows and 2 base-row columns.)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 3.1 here)\n",
- "tinyhex = HexModelGrid((3, 2), 2.0)\n",
- "plot_graph(tinyhex, at=\"node\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 3.2 here)\n",
- "hexht = tinyhex.add_zeros(\"height\", at=\"node\")\n",
- "hexht[3] = 0.5\n",
- "imshow_grid(tinyhex, hexht)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 3.3 here)\n",
- "plot_graph(tinyhex, at=\"link\")\n",
- "pred_grad = np.array([0, 0, 0.25, 0.25, 0, 0.25, -0.25, 0, -0.25, -0.25, 0, 0])\n",
- "print(pred_grad)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 3.4 here)\n",
- "hexgrad = tinyhex.calc_grad_at_link(hexht)\n",
- "print(hexgrad)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 3.5 here)\n",
- "hexflux = -0.01 * hexgrad\n",
- "print(hexflux)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 3.6 here)\n",
- "print(tinyhex.length_of_face)\n",
- "print(tinyhex.area_of_cell)\n",
- "total_outflux = 6 * 0.0025 * tinyhex.length_of_face[0]\n",
- "divergence = total_outflux / tinyhex.area_of_cell[0]\n",
- "print(total_outflux)\n",
- "print(divergence)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Part 4: Landlab Components\n",
- "\n",
- "Finally we will use a Landlab component, called the LinearDiffuser [link to its documentation](https://landlab.readthedocs.io/en/latest/reference/components/diffusion.html).\n",
- "\n",
- "Landlab was designed to have many of the utilities like `calc_grad_at_link`, and `calc_flux_divergence_at_node` to help you make your own models. Sometimes, however, you may use such a model over and over and over. Then it is nice to be able to put it in its own python class with a standard interface. \n",
- "\n",
- "This is what a Landlab Component is. \n",
- "\n",
- "There is a whole [tutorial on components](../component_tutorial/component_tutorial.ipynb) and a [page on the User Guide](https://landlab.readthedocs.io/en/latest/user_guide/components.html). For now we will just show you what the prior example looks like if we use the LinearDiffuser. \n",
- "\n",
- "First we import it, set up the grid, and uplift our fault block. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from landlab.components import LinearDiffuser\n",
- "\n",
- "mg = HexModelGrid((25, 40), 10, node_layout=\"rect\")\n",
- "z = mg.add_zeros(\"topographic__elevation\", at=\"node\")\n",
- "fault_trace_y = 50.0 + 0.25 * mg.x_of_node\n",
- "z[mg.y_of_node > fault_trace_y] += (\n",
- " 10.0 + 0.01 * mg.x_of_node[mg.y_of_node > fault_trace_y]\n",
- ")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Next we instantiate a LinearDiffuser. We have to tell the component what value to use for the diffusivity. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "ld = LinearDiffuser(mg, linear_diffusivity=D)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Finally we run the component forward in time and plot. Like many Landlab components, the LinearDiffuser has a method called \"run_one_step\" that takes one input, the timestep dt. Calling this method runs the LinearDiffuser forward in time by an increment dt. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "for i in range(25):\n",
- " ld.run_one_step(dt)\n",
- "imshow_grid(mg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Exercises for section 4\n",
- "\n",
- "(4.1) Repeat the steps above that instantiate and run a `LinearDiffuser` component, but this time give it a `RasterModelGrid`. Use `imshow_grid` to display the topography below."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 4.1 here)\n",
- "rmg = RasterModelGrid((25, 40), 10)\n",
- "z = rmg.add_zeros(\"topographic__elevation\", at=\"node\")\n",
- "fault_trace_y = 50.0 + 0.25 * rmg.x_of_node\n",
- "z[rmg.y_of_node > fault_trace_y] += (\n",
- " 10.0 + 0.01 * rmg.x_of_node[rmg.y_of_node > fault_trace_y]\n",
- ")\n",
- "ld = LinearDiffuser(rmg, linear_diffusivity=D)\n",
- "for i in range(25):\n",
- " ld.run_one_step(dt)\n",
- "imshow_grid(rmg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(4.2) Using either a raster or hex grid (your choice) with a `topographic__elevation` field that is initially all zeros, write a modified version of the loop that adds uplift to the core nodes each iteration, at a rate of 0.0001 m/yr. Run the model for enough time to accumulate 10 meters of uplift. Plot the terrain to verify that the land surface height never gets higher than 10 m. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 4.2 here)\n",
- "rmg = RasterModelGrid((40, 40), 10) # while we're at it, make it a bit bigger\n",
- "z = rmg.add_zeros(\"topographic__elevation\", at=\"node\")\n",
- "ld = LinearDiffuser(rmg, linear_diffusivity=D)\n",
- "for i in range(50):\n",
- " ld.run_one_step(dt)\n",
- " z[rmg.core_nodes] += dt * 0.0001\n",
- "imshow_grid(rmg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(4.3) Now run the same model long enough that it reaches (or gets very close to) a dynamic equilibrium between uplift and erosion. What shape does the hillslope have? "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to 4.3 here)\n",
- "z[:] = 0.0\n",
- "uplift_rate = 0.0001\n",
- "for i in range(4000):\n",
- " ld.run_one_step(dt)\n",
- " z[rmg.core_nodes] += dt * uplift_rate\n",
- "imshow_grid(rmg, \"topographic__elevation\")\n",
- "plt.figure()\n",
- "plt.plot(rmg.x_of_node, z, \".\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "(BONUS CHALLENGE QUESTION) Derive an analytical solution for the cross-sectional shape of your steady-state hillslope. Plot this solution next to the actual model's cross-section."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "#### SOLUTION (derivation)\n",
- "\n",
- "##### Derivation of the original governing equation\n",
- "\n",
- "(Note: you could just start with the governing equation and go from there, but we include this here for completeness).\n",
- "\n",
- "Consider a topographic profile across a hillslope. The horizontal coordinate along the profile is $x$, measured from the left side of the profile (i.e., the base of the hill on the left side, where $x=0$). The horizontal coordinate perpendicular to the profile is $y$. Assume that at any time, the hillslope is perfectly symmetrical in the $y$ direction, and that there is no flow of soil in this direction.\n",
- "\n",
- "Now consider a vertical column of soil somewhere along the profile. The left side of the column is at position $x$, and the right side is at position $x+\\Delta x$, with $\\Delta x$ being the width of the column in the $x$ direction. The width of the column in the $y$ direction is $W$. The height of the column, $z$, is also the height of the land surface at that location. Height is measured relative to the height of the base of the slope (in other words, $z(0) = 0$).\n",
- "\n",
- "The total mass of soil inside the column, and above the slope base, is equal to the volume of soil material times its density times the fraction of space that it fills, which is 1 - porosity. Denoting soil particle density by $\\rho$ and porosity by $\\phi$, the soil mass in a column of height $h$ is\n",
- "\n",
- "$m = (1-\\phi ) \\rho \\Delta x W z$.\n",
- "\n",
- "Conservation of mass dictates that the rate of change of mass equals the rate of mass inflow minus the rate of mass outflow. Assume that mass enters or leaves only by (1) soil creep, and (2) uplift of the hillslope material relative to the elevation of the hillslope base. The rate of the latter, in terms of length per time, will be denoted $U$. The rate of soil creep at a particular position $x$, in terms of bulk volume (including pores) per time per width, will be denoted $q_s(x)$. With this definition in mind, mass conservation dictates that:\n",
- "\n",
- "$\\frac{\\partial (1-\\phi ) \\rho \\Delta x W z}{\\partial t} = \\rho (1-\\phi ) \\Delta x W U + \\rho (1-\\phi ) q_s(x) - \\rho (1-\\phi ) q_s(x+\\Delta x)$.\n",
- "\n",
- "Assume that porosity and density are steady and uniform. Then,\n",
- "\n",
- "$\\frac{\\partial z}{\\partial t} = U + \\frac{q_s(x) - q_s(x+\\Delta x)}{\\Delta x}$.\n",
- "\n",
- "Factoring out -1 from the right-most term, and taking the limit as $\\Delta x\\rightarrow 0$, we get a differential equation that expresses conservation of mass for this situation:\n",
- "\n",
- "$\\frac{\\partial z}{\\partial t} = U - \\frac{\\partial q_s}{\\partial x}$.\n",
- "\n",
- "Next, substitute the soil-creep rate law\n",
- "\n",
- "$q_s = -D \\frac{\\partial z}{\\partial x}$,\n",
- "\n",
- "to obtain\n",
- "\n",
- "$\\frac{\\partial z}{\\partial t} = U + D \\frac{\\partial^2 z}{\\partial x^2}$.\n",
- "\n",
- "##### Steady state\n",
- "\n",
- "Steady means $dz/dt = 0$. If we go back to the mass conservation law a few steps ago and apply steady state, we find\n",
- "\n",
- "$\\frac{dq_s}{dx} = U$.\n",
- "\n",
- "If you think of a hillslope that slopes down to the right, you can think of this as indicating that for every step you take to the right, you get another increment of incoming soil via uplift relative to baselevel. (Turns out it works the same way for a slope that angles down the left, but that's less obvious in the above math)\n",
- "\n",
- "Integrate to get:\n",
- "\n",
- "$q_s = Ux + C_1$, where $C_1$ is a constant of integration.\n",
- "\n",
- "To evaluate the integration constant, let's assume the crest of the hill is right in the middle of the profile, at $x=L/2$, with $L$ being the total length of the profile. Net downslope soil flux will be zero at the crest (where the slope is zero), so for this location:\n",
- "\n",
- "$q_s = 0 = UL/2 + C_1$, \n",
- "\n",
- "and therefore,\n",
- "\n",
- "$C_1 = -UL/2$, \n",
- "\n",
- "and\n",
- "\n",
- "$q_s = U (x - L/2)$.\n",
- "\n",
- "Now substitute the creep law for $q_s$ and divide both sides by $-D$:\n",
- "\n",
- "$\\frac{dz}{dx} = \\frac{U}{D} (L/2 - x)$.\n",
- "\n",
- "Integrate:\n",
- "\n",
- "$z = \\frac{U}{D} (Lx/2 - x^2/2) + C_2$.\n",
- "\n",
- "To evaluate $C_2$, recall that $z(0)=0$ (and also $z(L)=0$), so $C_2=0$. Hence, here's our analytical solution, which describes a parabola:\n",
- "\n",
- "$\\boxed{z = \\frac{U}{2D} (Lx - x^2)}$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# (enter your solution to the bonus challenge question here)\n",
- "L = 390.0 # hillslope length, m\n",
- "x_analytic = np.arange(0.0, L)\n",
- "z_analytic = 0.5 * (uplift_rate / D) * (L * x_analytic - x_analytic * x_analytic)\n",
- "plt.plot(rmg.x_of_node, z, \".\")\n",
- "plt.plot(x_analytic, z_analytic, \"r\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Hey, hang on a minute, that's not a very good fit! What's going on? \n",
- "\n",
- "Turns out our 2D hillslope isn't as tall as the idealized 1D profile because of the boundary conditions: with soil free to flow east and west as well as north and south, the crest ends up lower than it would be if it were perfectly symmetrical in one direction.\n",
- "\n",
- "So let's try re-running the numerical model, but this time with the north and south boundaries closed so that the hill shape becomes uniform in the $y$ direction:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "rmg = RasterModelGrid((40, 40), 10)\n",
- "z = rmg.add_zeros(\"topographic__elevation\", at=\"node\")\n",
- "rmg.set_closed_boundaries_at_grid_edges(False, True, False, True) # closed on N and S\n",
- "ld = LinearDiffuser(rmg, linear_diffusivity=D)\n",
- "for i in range(4000):\n",
- " ld.run_one_step(dt)\n",
- " z[rmg.core_nodes] += dt * uplift_rate\n",
- "imshow_grid(rmg, \"topographic__elevation\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "plt.plot(rmg.x_of_node, z, \".\")\n",
- "plt.plot(x_analytic, z_analytic, \"r\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "That's more like it!"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Congratulations on making it to the end of this tutorial!\n",
- "\n",
- "**Click here for more** Landlab tutorials"
- ]
- }
- ],
- "metadata": {
- "anaconda-cloud": {},
- "kernelspec": {
- "display_name": "Python 3 (ipykernel)",
- "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.11.0"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 1
-}