{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear diffusion exercise with Landlab\n", "\n", "This notebook is adapted from *Landscape Evolution Modeling with CHILD* by Gregory Tucker and Stephen Lancaster. This notebook was created by Nicole Gasparini at Tulane University." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "For tutorials on learning Landlab, click here: https://github.com/landlab/landlab/wiki/Tutorials\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**What is this notebook?**\n", "\n", "This notebook illustrates the evolution of landforms dominated by processes that result in linear diffusion of sediment. In other words, the downhill flow of soil is proportional to the (downhill) gradient of the land surface multiplied by a transport coefficient.\n", "\n", "The notebook first illustrates a simple example of a diffusing hillslope. We then provide a number of exercises for students to do on their own. This set of exercises is recomended for students in a quantitative geomorphology class, who have been introduced to the linear diffusion equation in class. \n", "\n", "**Application of linear diffusion transport law:**\n", "\n", "For relatively gentle, soil-mantled slopes, there is reasonably strong support for a transport law of the form:\n", "\\begin{equation}\n", "q_s = -D \\nabla z\n", "\\end{equation}\n", "where ${q}_s$ is the transport rate with dimensions of L$^2$T$^{-1}$; $D$ is a transport coefficient with dimensions of L$^2$T$^{-1}$; and $z$ is elevation. $\\nabla z$ is the gradient in elevation. If distance is increasing downslope, $\\nabla z$ is negative downslope, hence the negative in front of $D$. \n", " \n", "Changes in elevation, or erosion, are calculated from conservation of mass:\n", "\\begin{equation}\n", "\\frac{dz}{dt} = U-\\nabla q_s\n", "\\end{equation}\n", "where $U$ is the rock uplift rate, with dimensions LT$^{-1}$.\n", "\n", "**How will we explore this with Landlab?**\n", "\n", "We will use the Landlab component *LinearDiffuser*, which implements the equations above, to explore how hillslopes evolve when linear diffusion describes hillslope sediment transport. We will explore both steady state, here defined as erosion rate equal to rock uplift rate, and also how a landscape gets to steady state.\n", "\n", "The first example illustrates how to set-up the model and evolve a hillslope to steady state, along with how to plot some variables of interest. We assume that you have knowledge of how to derive the steady-state form of a uniformly uplifting, steady-state, diffusive hillslope. For more information on hillslope sediment transport laws, this paper is a great overview:\n", "\n", "Roering, Joshua J. (2008) \"How well can hillslope evolution models “explain” topography? Simulating soil transport and production with high-resolution topographic data.\" Geological Society of America Bulletin.\n", "\n", "Based on the first example, you are asked to first think about what will happen as you change a parameter, and then you explore this numerically by changing the code.\n", "\n", "Start at the top by reading each block of text and sequentially running each code block (shift - enter OR got to the _Cell_ pulldown menu at the top and choose _Run Cells_). \n", "\n", "Remember that you can always go to the _Kernel_ pulldown menu at the top and choose _Restart & Clear Output_ or _Restart & Run All_ if you change things and want to start afresh. If you just change one code block and rerun only that code block, only the parts of the code in that code block will be updated. (E.g. if you change parameters but don't reset the code blocks that initialize run time or topography, then these values will not be reset.) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Now on to the code example**\n", "\n", "Import statements. You should not need to edit this." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# below is to make plots show up in the notebook\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 1\n", "\n", "import numpy as np\n", "from matplotlib.pyplot import figure, legend, plot, show, title, xlabel, ylabel, ylim\n", "\n", "from landlab.plot.imshow import imshow_grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will create a grid with 41 rows and 5 columns, and dx is 5 m (a long, narrow, hillslope). The initial elevation is 0 at all nodes.\n", "\n", "We set-up boundary conditions so that material can leave the hillslope at the two short ends." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 2\n", "\n", "# setup grid\n", "from landlab import RasterModelGrid\n", "\n", "mg = RasterModelGrid((41, 5), 5.0)\n", "z_vals = mg.add_zeros(\"topographic__elevation\", at=\"node\")\n", "\n", "# initialize some values for plotting\n", "ycoord_rast = mg.node_vector_to_raster(mg.node_y)\n", "ys_grid = ycoord_rast[:, 2]\n", "\n", "# set boundary condition.\n", "mg.set_closed_boundaries_at_grid_edges(True, False, True, False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we import and initialize the *LinearDiffuser* component. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 3\n", "\n", "from landlab.components import LinearDiffuser\n", "\n", "D = 0.01 # initial value of 0.01 m^2/yr\n", "lin_diffuse = LinearDiffuser(mg, linear_diffusivity=D)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now initialize a few more parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 4\n", "\n", "# Uniform rate of rock uplift\n", "uplift_rate = 0.0001 # meters/year, originally set to 0.0001\n", "\n", "# Total time in years that the model will run for.\n", "runtime = 1000000 # years, originally set to 1,000,000\n", "\n", "# Stability criteria for timestep dt. Coefficient can be changed\n", "# depending on our tolerance for stability vs tolerance for run time.\n", "dt = 0.5 * mg.dx * mg.dx / D\n", "\n", "# nt is number of time steps\n", "nt = int(runtime // dt)\n", "\n", "# Below is to keep track of time for labeling plots\n", "time_counter = 0\n", "\n", "# length of uplift over a single time step, meters\n", "uplift_per_step = uplift_rate * dt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we figure out the analytical solution for the elevation of the steady-state profile." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 5\n", "\n", "ys = np.arange(mg.number_of_node_rows * mg.dx - mg.dx)\n", "\n", "# location of divide or ridge crest -> middle of grid\n", "# based on boundary conds.\n", "divide_loc = (mg.number_of_node_rows * mg.dx - mg.dx) / 2\n", "\n", "# half-width of the ridge\n", "half_width = (mg.number_of_node_rows * mg.dx - mg.dx) / 2\n", "\n", "# analytical solution for elevation under linear diffusion at steady state\n", "zs = (uplift_rate / (2 * D)) * (np.power(half_width, 2) - np.power(ys - divide_loc, 2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we evolve the landscape, let's look at the initial topography. (This is just verifying that it is flat with zero elevation.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 6\n", "\n", "figure(1)\n", "imshow_grid(mg, \"topographic__elevation\")\n", "title(\"initial topography\")\n", "figure(2)\n", "elev_rast = mg.node_vector_to_raster(mg.at_node[\"topographic__elevation\"])\n", "plot(ys_grid, elev_rast[:, 2], \"r-\", label=\"model\")\n", "plot(ys, zs, \"k--\", label=\"analytical solution\")\n", "ylim((-5, 50)) # may want to change upper limit if D changes\n", "xlabel(\"horizontal distance (m)\")\n", "ylabel(\"vertical distance (m)\")\n", "legend(loc=\"lower center\")\n", "title(\"initial topographic cross section\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to evolve the landscape and compare it to the steady state solution.\n", "\n", "Below is the time loop that does all the calculations. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 7\n", "\n", "for i in range(nt):\n", " mg[\"node\"][\"topographic__elevation\"][mg.core_nodes] += uplift_per_step\n", " lin_diffuse.run_one_step(dt)\n", " time_counter += dt\n", "\n", " # All landscape evolution is the first two lines of loop.\n", " # Below is simply for plotting the topography halfway through the run\n", " if i == int(nt // 2):\n", " figure(1)\n", " imshow_grid(mg, \"topographic__elevation\")\n", " title(\"topography at time %s, with D = %s\" % (time_counter, D))\n", " figure(2)\n", " elev_rast = mg.node_vector_to_raster(mg.at_node[\"topographic__elevation\"])\n", " plot(ys_grid, elev_rast[:, 2], \"k-\", label=\"model\")\n", " plot(ys, zs, \"g--\", label=\"analytical solution - SS\")\n", " plot(ys, zs * 0.75, \"b--\", label=\"75% of analytical solution\")\n", " plot(ys, zs * 0.5, \"r--\", label=\"50% of analytical solution\")\n", " xlabel(\"horizontal distance (m)\")\n", " ylabel(\"vertical distance (m)\")\n", " legend(loc=\"lower center\")\n", " title(\"topographic__elevation at time %s, with D = %s\" % (time_counter, D))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we plot the final cross-section." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 8\n", "\n", "elev_rast = mg.node_vector_to_raster(mg.at_node[\"topographic__elevation\"])\n", "plot(ys_grid, elev_rast[:, 2], \"k-\", label=\"model\")\n", "plot(ys, zs, \"g--\", label=\"analytical solution - SS\")\n", "plot(ys, zs * 0.75, \"b--\", label=\"75% of analytical solution\")\n", "plot(ys, zs * 0.5, \"r--\", label=\"50% of analytical solution\")\n", "xlabel(\"horizontal distance (m)\")\n", "ylabel(\"vertical distance (m)\")\n", "legend(loc=\"lower center\")\n", "title(\"topographic cross section at time %s, with D = %s\" % (time_counter, D))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we plot the steepest slope in the downward direction across the landscape.\n", "\n", "(To calculate the steepest slope at a location, we need to route flow across the landscape.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Code Block 9\n", "\n", "from landlab.components import FlowAccumulator\n", "\n", "fr = FlowAccumulator(mg) # intializing flow routing\n", "fr.run_one_step()\n", "plot(\n", " mg.node_y[mg.core_nodes],\n", " mg.at_node[\"topographic__steepest_slope\"][mg.core_nodes],\n", " \"k-\",\n", ")\n", "xlabel(\"horizontal distance (m)\")\n", "ylabel(\"topographic slope (m/m)\")\n", "title(\"slope of the hillslope at time %s, with D = %s\" % (time_counter, D))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Has the landscape reached steady state yet? How do you know?\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Answer: Not quite, but it is getting close. Go back and rerun Code Blocks 7, 8 and 9 (time loop and plotting). (Remember you can rerun a cell with shift-return, or from the cell pull-down menu.) Has it reached steady state yet? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**What to do and hand in:**\n", "1. In the example illustrated here ($D$ = 0.01 m$^2$yr$^{-1}$ and $U$ = 0.0001 m yr$^{-1}$). Restart everything, and use the model to determine how long it takes for the landscape to go from a flat to reach 50%, 75% and 100% of its steady-state morphology. Does the landscape approach steady state linearly in time? (You can run the time loop (Code Block 7) multiple times without running other code blocks again to continually evolve the landscape. You will initially want to rerun all the code blocks and change the value of **run_time** (Code Block 4). Determining the correct value of **run_time** to use will take some iteration.)\n", "2. What do you think will happen when you increase $D$ (Code Block 3) by a factor of 10? Will the time to steady state differ? If yes, how? Will the topography be different? If yes, how and why? What does it mean physically, about processes, if $D$ increases? Answer these questions before running any code. \n", "3. Now set $D$ = 0.1 m$^2$yr$^{-1}$ and rerun landscape evolution from an initial flat. Illustrate the final steady state topography and record the time to steady state. Discuss how the landscape differs from the results in question 1. Discuss how the results are similar to or different from your intuition. It is OK if your intuition was wrong! \n", "4. What do you think will happen when you increase **uplift_rate** (Code Block 4) by a factor of 10? Will the time to steady state differ? If yes, how? Will the topography be different? If yes, how and why? Answer these questions first, and then rerun the code with **uplift_rate** = 0.001 m yr$^{-1}$. (Make sure you change $D$ - Code Block 3 - back to the original value of 0.01 m$^2$yr$^{-1}$ and restart from a flat surface.) Illustrate the final steady state topography. Discuss how these results differ from the results in question 1 and how the results match (or do not) your intuition. It is OK if your intuition was wrong.\n", "\n", "You should hand in a typed document that answers the above questions with supporting plots. Plots should be embedded in the text, or, if they all fall at the end, they need to be clearly labeled, e.g. each plot has a figure number and plots are referred to by figure number in the text.\n", "\n", "Other questions you can explore.\n", "\n", "1. What happens to time to steady state as you increase the length of your hillslope? \n", "2. Does grid resolution affect your answers? If so, how?\n" ] } ], "metadata": { "anaconda-cloud": {}, "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.0" } }, "nbformat": 4, "nbformat_minor": 1 }