stream_power_channels_steepness_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"http://landlab.github.io\"><img style=\"float: left\" src=\"../../../landlab_header.png\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stream power, channel steepness, and relief \n",
    "This notebook was created by Nicole Gasparini at Tulane University."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "<small>For tutorials on learning Landlab, click here: <a href=\"https://github.com/landlab/landlab/wiki/Tutorials\">https://github.com/landlab/landlab/wiki/Tutorials</a></small>\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What is this notebook?**\n",
    "\n",
    "This notebook illustrates the evolution of detachment-limited channels in an actively uplifting landscape. The landscape evolves according to the equation:\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{d z}{d t} = -K_\\text{sp} A^{m_{sp}} S^{n_{sp}} + U\n",
    "\\end{equation}\n",
    "Here, $K_{sp}$ is the erodibility coefficient on fluvial incision, which is thought to be positively correlated with climate wetness, or storminess (this is hard to quantify) and to be negatively correlated with rock strength (again, rock strength is hard to quantify). $m_{sp}$ and $n_{sp}$ are positive exponents, usually thought to have a ratio, $m_{sp}/n_{sp} \\approx 0.5$. $A$ is drainage area and $S$ is the slope of steepest descent ($-\\frac{dz}{dx}$) where $x$ is horizontal distance (positive in the downslope direction) and $z$ is elevation. (If slope is negative there is no fluvial erosion.) $U$ is an externally-applied rock uplift rate field.\n",
    "\n",
    "The fluvial erosion term is also known as the stream power equation. Before using this notebook you should be familiar with this equation from class lectures and reading. \n",
    "\n",
    "For a great overview of the stream power equation, see: \n",
    "\n",
    "- Whipple and Tucker, 1999, Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs, Journal of Geophysical Research.\n",
    "\n",
    "For some great illustrations of modeling with the sream power equation, see:\n",
    "\n",
    "- Tucker and Whipple, 2002, Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodel comparison, Journal of Geophysical Research.\n",
    "\n",
    "Helpful background on landscape sensitivity to rock uplift rates and patterns can be found here:\n",
    "\n",
    "- Kirby and Whipple, 2012, Expression of active tectonics in erosional landscapes, Journal of Structural Geology.\n",
    "\n",
    "**What will you do?**\n",
    "\n",
    "In this exercise you will modify the code to get a better understanding of how rock uplift rates and patterns and the erodibility coefficient control fluvial channel form.\n",
    "\n",
    "Start at the top by reading each block of text and sequentially running each code block (shift - enter OR got to the _Run_ pulldown menu at the top and choose _Run Selected Cells_). \n",
    "\n",
    "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 the cahnged parameters will apply on the current topography.) \n",
    "\n",
    "**STUDENTS: Questions to answer before starting this assignment.**\n",
    "\n",
    "Answer these questions before running the notebook.\n",
    "\n",
    "1. What do you think will happen to total relief (defined as the maximum minus the minimum elevation, here area is fixed) and channel slope at steady state if $K_{sp}$ is uniformly increased?\n",
    "2. What do you think will happen to total relief and channel slope at steady state if $U$ is uniformly increased?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now on to the code...**\n",
    "\n",
    "First we have to import the parts of Python and Landlab that are needed to run this code. You should not have to change this first code block.\n",
    "\n",
    "**NOTE** If you are running this on the EarthscapeHub, the kernal (in the upper right hand corner of this notebook) should be set to CSDMS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code block 1\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from landlab import RasterModelGrid, imshow_grid\n",
    "from landlab.components import (\n",
    "    ChannelProfiler,\n",
    "    ChiFinder,\n",
    "    FlowAccumulator,\n",
    "    SteepnessFinder,\n",
    "    StreamPowerEroder,\n",
    ")\n",
    "from landlab.io import write_esri_ascii"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a grid and set boundary conditions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 2\n",
    "\n",
    "number_of_rows = 100  # number of raster cells in vertical direction (y)\n",
    "number_of_columns = 200  # number of raster cells in horizontal direction (x)\n",
    "dxy = 100  # side length of a raster model cell, or resolution [m]\n",
    "\n",
    "# Below is a raster (square cells) grid, with equal width and height\n",
    "mg1 = RasterModelGrid((number_of_rows, number_of_columns), dxy)\n",
    "\n",
    "# Set boundary conditions - only the south side of the grid is open.\n",
    "# Boolean parameters are sent to function in order of\n",
    "# east, north, west, south.\n",
    "mg1.set_closed_boundaries_at_grid_edges(True, True, True, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we make the initial grid of elevation of zeros with a very small amount of noise to make a more pleasing network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 3\n",
    "\n",
    "np.random.seed(77)  # seed set so our figures are reproducible\n",
    "mg1_noise = (np.random.rand(mg1.number_of_nodes) / 1000.0\n",
    "             )  # intial noise on elevation gri\n",
    "\n",
    "# set up the elevation on the grid\n",
    "z1 = mg1.add_zeros(\"topographic__elevation\", at=\"node\")\n",
    "z1 += mg1_noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set parameters related to time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 4\n",
    "\n",
    "tmax = 1e6  # time for the model to run [yr] (Original value was 2E6 yr)\n",
    "dt = 5000  # time step [yr] (Original value was 5000 yr)\n",
    "total_time = 0  # amount of time the landscape has evolved [yr]\n",
    "# total_time will increase as you keep running the code.\n",
    "\n",
    "t = np.arange(0, tmax, dt)  # each of the time steps that the code will run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set parameters for incision and intializing all of the process components that do the work. We also initialize tools for quantifying the landscape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 5\n",
    "\n",
    "# Original K_sp value is 1e-5\n",
    "K_sp = 1.0e-5  # units vary depending on m_sp and n_sp\n",
    "m_sp = 0.5  # exponent on drainage area in stream power equation\n",
    "n_sp = 1.0  # exponent on slope in stream power equation\n",
    "\n",
    "frr = FlowAccumulator(mg1, flow_director='FlowDirectorD8')  # intializing flow routing\n",
    "spr = StreamPowerEroder(mg1, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp,\n",
    "                        threshold_sp=0.0)  # initializing stream power incision\n",
    "\n",
    "theta = m_sp / n_sp\n",
    "# initialize the component that will calculate channel steepness\n",
    "sf = SteepnessFinder(mg1, reference_concavity=theta, min_drainage_area=1000.0)\n",
    "# initialize the component that will calculate the chi index\n",
    "cf = ChiFinder(mg1,\n",
    "               min_drainage_area=1000.0,\n",
    "               reference_concavity=theta,\n",
    "               use_true_dx=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize rock uplift rate. This will need to be changed later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 6\n",
    "\n",
    "#  uplift_rate [m/yr] (Original value is 0.0001 m/yr)\n",
    "uplift_rate = np.ones(mg1.number_of_nodes) * 0.0001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the code loop. \n",
    "\n",
    "You can rerun Code Block 7 many times, and as long as you don't reset the elevation field (Code Block 3), it will take the already evolved landscape and evolve it even more. If you want to change parameters in other code blocks (e.g. Code Block 5 or 6), you can do that too, and as long as you don't reset the elevation field (Code Block 3) the new parameters will apply on the already evolved topography. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Code Block 7\n",
    "\n",
    "for ti in t:\n",
    "    z1[mg1.\n",
    "       core_nodes] += uplift_rate[mg1.core_nodes] * dt  # uplift the landscape\n",
    "    frr.run_one_step()  # route flow\n",
    "    spr.run_one_step(dt)  # fluvial incision\n",
    "    total_time += dt  # update time keeper\n",
    "    print(total_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the topography."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 8\n",
    "\n",
    "imshow_grid(mg1,\n",
    "            \"topographic__elevation\",\n",
    "            grid_units=(\"m\", \"m\"),\n",
    "            var_name=\"Elevation (m)\")\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]*1000} mm/yr; $time$={total_time/1000} kyr\"\n",
    "plt.title(title_text)\n",
    "\n",
    "max_elev = np.max(z1)\n",
    "print(\"Maximum elevation is \", np.max(z1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the slope and area data at each point on the landscape (in log-log space). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Code Block 9\n",
    "\n",
    "plt.loglog(\n",
    "    mg1.at_node[\"drainage_area\"][np.where(mg1.node_y[mg1.core_nodes]>200)],\n",
    "    mg1.at_node[\"topographic__steepest_slope\"][np.where(mg1.node_y[mg1.core_nodes]>200)],\n",
    "    \"b.\",\n",
    ")\n",
    "plt.ylabel(\"Topographic slope\")\n",
    "plt.xlabel(\"Drainage area (m^2)\")\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]*1000} mm/yr; $time$={total_time/1000} kyr\"\n",
    "\n",
    "plt.title(title_text)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is slightly easier to interpret slope-area data when we look at a single channel or channels, rather than the entire landscape. Below we plot the profile and slope-area data for the three largest channels on the landscape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 10\n",
    "\n",
    "# profile the largest channels, set initially to find the mainstem channel in the three biggest watersheds\n",
    "# you can change the number of watersheds, or choose to plot all the channel segments in the watershed that\n",
    "# have drainage area below the threshold (here we have set the threshold to the area of a grid cell).\n",
    "prf = ChannelProfiler(mg1,\n",
    "                      number_of_watersheds=3,\n",
    "                      main_channel_only=True,\n",
    "                      minimum_channel_threshold=dxy**2)\n",
    "prf.run_one_step()\n",
    "\n",
    "# plot the elevation as a function of distance upstream\n",
    "plt.figure(1)\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]*1000} mm/yr; $time$={total_time/1000} kyr\"\n",
    "prf.plot_profiles(xlabel='distance upstream (m)',\n",
    "                  ylabel='elevation (m)',\n",
    "                  title=title_text)\n",
    "ax=plt.gca()\n",
    "ax.invert_xaxis()\n",
    "\n",
    "# plot the location of the channels in map view\n",
    "plt.figure(2)\n",
    "prf.plot_profiles_in_map_view()\n",
    "\n",
    "# slope-area data in just the profiled channels\n",
    "plt.figure(3)\n",
    "for i, outlet_id in enumerate(prf.data_structure):\n",
    "    for j, segment_id in enumerate(prf.data_structure[outlet_id]):\n",
    "        if j == 0:\n",
    "            label = \"channel {i}\".format(i=i + 1)\n",
    "        else:\n",
    "            label = '_nolegend_'\n",
    "        segment = prf.data_structure[outlet_id][segment_id]\n",
    "        profile_ids = segment[\"ids\"]\n",
    "        color = segment[\"color\"]\n",
    "        plt.loglog(\n",
    "            mg1.at_node[\"drainage_area\"][profile_ids],\n",
    "            mg1.at_node[\"topographic__steepest_slope\"][profile_ids],\n",
    "            '.',\n",
    "            color=color,\n",
    "            label=label,\n",
    "        )\n",
    "\n",
    "plt.legend(loc=\"lower left\")\n",
    "plt.xlabel(\"drainage area (m^2)\")\n",
    "plt.ylabel(\"channel slope [m/m]\")\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]*1000} mm/yr; $time$={total_time/1000} kyr\"\n",
    "plt.title(title_text)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The chi index is a useful way to quantitatively interpret fluvial channels. Below we plot the chi index in the three largest channels and also a chi map across the entire landscape. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 11\n",
    "\n",
    "# calculate the chi index\n",
    "cf.calculate_chi()\n",
    "\n",
    "# # chi-elevation plots in the profiled channels\n",
    "plt.figure(4)\n",
    "\n",
    "for i, outlet_id in enumerate(prf.data_structure):\n",
    "    for j, segment_id in enumerate(prf.data_structure[outlet_id]):\n",
    "        if j == 0:\n",
    "            label = \"channel {i}\".format(i=i + 1)\n",
    "        else:\n",
    "            label = '_nolegend_'\n",
    "        segment = prf.data_structure[outlet_id][segment_id]\n",
    "        profile_ids = segment[\"ids\"]\n",
    "        color = segment[\"color\"]\n",
    "        plt.plot(\n",
    "            mg1.at_node[\"channel__chi_index\"][profile_ids],\n",
    "            mg1.at_node[\"topographic__elevation\"][profile_ids],\n",
    "            color=color,\n",
    "            label=label,\n",
    "        )\n",
    "        \n",
    "ax=plt.gca()\n",
    "ax.invert_xaxis()\n",
    "plt.xlabel(\"chi index (m)\")\n",
    "plt.ylabel(\"elevation (m)\")\n",
    "plt.legend(loc=\"upper right\")\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]} m/yr; $time$={total_time/1000} kyr; concavity={theta}\"\n",
    "plt.title(title_text)\n",
    "\n",
    "# chi map\n",
    "plt.figure(5)\n",
    "imshow_grid(\n",
    "    mg1,\n",
    "    \"channel__chi_index\",\n",
    "    grid_units=(\"m\", \"m\"),\n",
    "    var_name=\"Chi index (m)\",\n",
    "    cmap=\"jet\",\n",
    ")\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]*1000} mm/yr; $time$={total_time/1000} kyr; concavity={theta}\"\n",
    "plt.title(title_text)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The channel steepness index is another useful index to quantify fluvial channels. Below we plot the steepness index in the same three largest channels, and also plot steepness index across the grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 12\n",
    "\n",
    "# calculate channel steepness\n",
    "sf.calculate_steepnesses()\n",
    "\n",
    "# plots of steepnes vs. distance upstream in the profiled channels\n",
    "plt.figure(6)\n",
    "\n",
    "for i, outlet_id in enumerate(prf.data_structure):\n",
    "    for j, segment_id in enumerate(prf.data_structure[outlet_id]):\n",
    "        if j == 0:\n",
    "            label = \"channel {i}\".format(i=i + 1)\n",
    "        else:\n",
    "            label = '_nolegend_'\n",
    "        segment = prf.data_structure[outlet_id][segment_id]\n",
    "        profile_ids = segment[\"ids\"]\n",
    "        distance_upstream = segment[\"distances\"]\n",
    "        color = segment[\"color\"]\n",
    "        plt.plot(\n",
    "            distance_upstream,\n",
    "            mg1.at_node[\"channel__steepness_index\"][profile_ids],\n",
    "            'x',\n",
    "            color=color,\n",
    "            label=label,\n",
    "        )\n",
    "\n",
    "ax=plt.gca()\n",
    "ax.invert_xaxis()\n",
    "plt.xlabel(\"distance upstream (m)\")\n",
    "plt.ylabel(\"steepness index\")\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.title(\n",
    "    f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]*1000} mm/yr; $time$={total_time/1000} kyr; concavity={theta}\"\n",
    ")\n",
    "\n",
    "# channel steepness map\n",
    "plt.figure(7)\n",
    "imshow_grid(\n",
    "    mg1,\n",
    "    \"channel__steepness_index\",\n",
    "    grid_units=(\"m\", \"m\"),\n",
    "    var_name=\"Steepness index \",\n",
    "    cmap=\"jet\",\n",
    ")\n",
    "#title_text = (\"$K_{sp}$=\" + str(K_sp) + \"; $time$=\" + str(total_time) +\n",
    "#              \"yr; $dx$=\" + str(dxy) + \"m\" + \"; concavity=\" + str(theta))\n",
    "plt.title(\n",
    "    f\"$K_{{sp}}$={K_sp}; $U$={uplift_rate[50]*1000} mm/yr; $time$={total_time/1000} kyr; concavity={theta}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Code Block 13\n",
    "\n",
    "# from landlab.utils import window_statistic\n",
    "\n",
    "# radius = 500\n",
    "\n",
    "# relief = window_statistic.calculate_window_statistic(mg1,\"topographic__elevation\", np.ptp,search_radius=radius, calc_on_closed_nodes=False,)\n",
    "\n",
    "# plt.figure(8)\n",
    "# imshow_grid(mg1, relief/radius)\n",
    "# plt.title(\"map of local relief with window radius = %s m\" % (radius))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have a grid that you want to export, uncomment and edit the appropriate lines below and run the code block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 14\n",
    "\n",
    "## Below has the name of the file that data will be written to.\n",
    "## You need to change the name of the file every time that you want\n",
    "## to write data, otherwise you will get an error.\n",
    "## This will write to the directory that you are running the code in.\n",
    "# write_file_name = 'data_file.txt'\n",
    "## Below is writing elevation data in the ESRI ascii format so that it can\n",
    "## easily be read into Arc GIS or back into Landlab.\n",
    "# write_esri_ascii(write_file_name, mg1, 'topographic__elevation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After running every code block once, has the landscape reached steady state? Answer: NO! How do you know? After you think about this, you are ready to complete this project.\n",
    "\n",
    "Answer the following questions using the code above and below. All answers and supporting figures (produced using the code) should be embedded in a file (e.g. ppt, txt, link to a google doc, Jupyter Notebook, ...) that you hand in. Code Blocks 8-12 produce different figures that you may find useful. You can use any or all of these different figures to help you with the questions below.  \n",
    "\n",
    "Anything with a question mark should be answered in the file that you hand in. Make sure your write in full sentences and proofread the document that you hand in.\n",
    "\n",
    "1. **Steady state with low rock uplift rate.** Using the parameters provided in the initial notebook, run the landscape to steady state. (Note that you can keep running the main evolution loop - Code Block 7 - and the different plotting blocks without running the code blocks above them. You may also want to change $tmax$ in Code Block 4.) How did you decide that the landscape reached steady state? [These landscapes may not reach a perfect steady state. Close is fine.] Include appropriate plots. \n",
    "\n",
    "2. **Steady-state landscape with higher rock uplift rate.** Now increase rock uplift uniformly by a factor of 2 to 0.0002 m/yr (Code Block 6). Rerun the entire notebook, and continue to run the main evolution loop - Code Block 7 - until the landscape reaches steady state. Provide a plot that illustrates that the landscape is in steady state. What aspects of the landscape have changed in comparison with the base landscape from question 1?\n",
    "\n",
    "3. **Increase erodibility.** Now increase $K_{sp}$ to 2E-5 (Code Block 5). Make sure rock uplift rate is set to the original value of 0.0001 m/yr (Code Block 6). Rerun the entire notebook, and continue to run the main evolution loop - Code Block 7 - until the landscape reaches steady state. Compare this landscape with the one produced in the first question. What, if anything, has changed? \n",
    "\n",
    "4. **General Reflection** In 4 sentences or less, describe how the steady-state landscape form changes with rock uplift rate? What metrics could you measure in the landscape that would illustrate higher rock uplift rate? In 4 sentences or less, describe how the steady-state landscape form changes with erodibility? What metrics could you measure in different landscapes that would illustrate they have different erodibilities? What might lead to different erodibilities in different landscapes?\n",
    "\n",
    "5. **Final Reflection** Was your initial insight into how parameters would affect the landscape correct? Discuss in 6 sentences or less."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "CSDMS",
   "language": "python",
   "name": "csdms-2023"
  },
  "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}