stream_power_channels_notebook_part2.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": [
    "# Quantifying river channel evolution with Landlab\n",
    "These exercises are based on a project orginally designed by Kelin Whipple at Arizona State University. 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 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 _Cell_ pulldown menu at the top and choose _Run 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 these values will not be reset.) \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?\n",
    "3. How do you think a steady-state landscape with a uniform low rock uplift rate will respond if rock uplift is uniformly increased (relative to a steady base level)? How will channel slopes change through time?"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code block 1\n",
    "\n",
    "import copy\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 = 50  # number of raster cells in vertical direction (y)\n",
    "number_of_columns = 100  # number of raster cells in horizontal direction (x)\n",
    "dxy = 200  # 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(35)  # seed set so our figures are reproducible\n",
    "mg1_noise = (\n",
    "    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 = 5e5  # time for the model to run [yr] (Original value was 5E5 yr)\n",
    "dt = 1000  # time step [yr] (Original value was 100 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(\n",
    "    mg1, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp, threshold_sp=0.0\n",
    ")  # 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(\n",
    "    mg1, min_drainage_area=1000.0, reference_concavity=theta, use_true_dx=True\n",
    ")"
   ]
  },
  {
   "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",
    "Note that 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.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(\n",
    "    mg1, \"topographic__elevation\", grid_units=(\"m\", \"m\"), var_name=\"Elevation (m)\"\n",
    ")\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m\"\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). We will only plot the core nodes because the boundary nodes have slopes that are influenced by the boundary conditions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Code Block 9\n",
    "\n",
    "plt.loglog(\n",
    "    mg1.at_node[\"drainage_area\"][mg1.core_nodes],\n",
    "    mg1.at_node[\"topographic__steepest_slope\"][mg1.core_nodes],\n",
    "    \"b.\",\n",
    ")\n",
    "plt.ylabel(\"Topographic slope\")\n",
    "plt.xlabel(\"Drainage area (m^2)\")\n",
    "title_text = f\"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m\"\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, 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(\n",
    "    mg1,\n",
    "    number_of_watersheds=3,\n",
    "    main_channel_only=True,\n",
    "    minimum_channel_threshold=dxy**2,\n",
    ")\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}; $time$={total_time} yr; $dx$={dxy} m\"\n",
    "prf.plot_profiles(\n",
    "    xlabel=\"distance upstream (m)\", ylabel=\"elevation (m)\", title=title_text\n",
    ")\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}; $time$={total_time} yr; $dx$={dxy} m\"\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",
    "plt.xlabel(\"chi index (m)\")\n",
    "plt.ylabel(\"elevation (m)\")\n",
    "plt.legend(loc=\"lower right\")\n",
    "title_text = (\n",
    "    f\"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m; concavity={theta}\"\n",
    ")\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 = (\n",
    "    f\"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m; concavity={theta}\"\n",
    ")\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",
    "plt.xlabel(\"distance upstream (m)\")\n",
    "plt.ylabel(\"steepness index\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.title(f\"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m; concavity={theta}\")\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 = (\n",
    "    \"$K_{sp}$=\"\n",
    "    + str(K_sp)\n",
    "    + \"; $time$=\"\n",
    "    + str(total_time)\n",
    "    + \"yr; $dx$=\"\n",
    "    + str(dxy)\n",
    "    + \"m\"\n",
    "    + \"; concavity=\"\n",
    "    + str(theta)\n",
    ")\n",
    "plt.title(f\"$K_{{sp}}$={K_sp}; $time$={total_time} yr; $dx$={dxy} m; concavity={theta}\")"
   ]
  },
  {
   "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 13\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": "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 should be typed, and supporting figures (produced using the code) should be embedded in one document that you hand in. Code Blocks 8-12 and 18-21 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. (Download or screenshoot the figures.) \n",
    "\n",
    "Anything with a question mark should be answered in the document 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 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 know that the landscape reached steady state? Note the approximate time that it took to reach steady state for your own reference. (This will be usefull for later questions.) Include appropriate plots. (If you want to analyze these landscapes outside of Landlab or save for later, make sure you save the elevation data to a text file (Code Block 13).)\n",
    "\n",
    "** NOTE, For the rest of the questions you should use Code Blocks 14 - 21. These will allow you to use the steady-state landscape created for question 1 - referred to here as the 'base landscape' - as the initial condition. Start by editing what you need to in Code Blocks 14 - 16. Run these each once, sequentially. You can run Code Block 17, the time loop, as many times as you need to, along with Code Blocks 18-21, which produce plots.**\n",
    "\n",
    "2. **Transient landscape responding to an increase in rock uplift. ** Use the base landscape and increase rock uplift uniformly by a factor of 4 to 0.0004 m/yr. Make sure you update the rock uplift rate (Code Block 16) and ensure that $tmax$ is 1e5 yrs and $dt$ is 500 yrs (Code Block 15). Run this until the maximum elevation in the grid is ~ 170 m and observe how the landscape gets to this elevation, i.e. plot intermediate steps. What patterns do you see in the supporting plots that illustrate this type of transient? Which patterns, if any, are diagnostic of a landscape response to uniform increase in rock uplift rate? (You may need to answer this after completing all of the questions.)\n",
    "\n",
    "3. ** Steady-state landscape with increased rock uplift. ** Now run the landscape from question 2 until it reaches steady state. (I.e. run the time loop, Code Block 17, a bunch of times. You can increase $tmax$ and $dt$ to make this run faster.) 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",
    "4. ** Increase erodibility. ** Start again from the base landscape, but this time increase $K_{sp}$ to 2E-5 (Code Block 14). Make sure rock uplift rate is set to the original value of 0.0001 m/yr (Code Block 16). Set $tmax$ to 1e5 yrs (Code Block 15). Run for 1e5 yrs and save the plots that you think are diagnostic. Run for another 1e5 yrs and save plots again. Now run for 5e5 yrs and save plots again. Quantitatively describe how the landscape evolves in response to the increase in erodibility and provide supporting plots. What could cause a uniform increase in erodibility?\n",
    "\n",
    "5. ** Spatially varible uplift - discrete, massive earthquake. ** Start again from the base landscape, and make sure that $K_{sp}$ = 1E-5 (Code Block 14). Now add a seismic event to this steady state landscape - a fault that runs horizontally across the landscape at y = 4000 m, and instantaneously uplifts half the landscape by 10 meters (Code Block 16). In this case, we will keep background uplift uniform at 0.0001 m/yr. Set $tmax$ to 1e5 yrs and $dt$ to 500 yrs (Code Block 15) before evolving the landscape after the fault. Now run the time loop four times and look at the different plots after each loop. How does the landscape respond to this fault? What patterns do you see in the supporting plots that illustrate this type of transient? Which patterns, if any, are diagnostic of a channel response to an earthquake? (You may need to answer this after completing all of the questions.)\n",
    "\n",
    "6. ** Spatially Varible Rock Uplift - discrete fault with two different uplift rates. ** Start again from the base landscape, and make sure that $K_{sp}$ = 1E-5 (Code Block 14). Now we will add a fault (at y = 4000 m) to this landscape. In this case the uplift rate on the footwall is higher (0.0004 m/yr) than on the hanging wall (uplift rate = 0.0001 m/yr). (Edit Code Block 16.) Set $tmax$ to 1e5 yrs and $dt$ to 500 yrs (Code Block 15). Now run the time loop four separate times and look at the different plots after each loop. How does the landscape respond to this fault? What patterns do you see in the supporting plots that illustrate this type of transient? Which patterns, if any, are diagnostic of a channel response to a this type of gradient in rock uplift rates? (You may need to answer this after completing all of the questions.)\n",
    "\n",
    "7. ** Spatially Varible Rock Uplift - gradient in uplift across the range. ** Start again from the base landscape, and make sure that $K_{sp}$ = 1E-5 (Code Block 14). Now we will add a linear gradient in uplift rate across the entire range (edit Code Block 16). The maximum uplift rate will be 0.0004 m/yr at the core of the range, and 0.0001 m/yr at the front of the range. Set $tmax$ to 1e5 yrs (Code Block 4) and $dt$ to 500 yrs before you start running the time loop for the fault before you start running the time loop with the rock uplift gradient. Now run the time loop four separate times and look at the different plots after each loop. How does the landscape respond to this gradient in uplift rate? What patterns do you see in the supporting plots that illustrate this type of transient? Which patterns, if any, are diagnostic of a channel response to this type of gradient in rock uplift rates? (You may need to answer this after completing all of the questions.)\n",
    "\n",
    "8. ** 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": [
    "# Code Block 14\n",
    "\n",
    "number_of_rows = 50  # number of raster cells in vertical direction (y)\n",
    "number_of_columns = 100  # number of raster cells in horizontal direction (x)\n",
    "dxy2 = 200  # 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",
    "mg2 = RasterModelGrid((number_of_rows, number_of_columns), dxy2)\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",
    "mg2.set_closed_boundaries_at_grid_edges(True, True, True, False)\n",
    "\n",
    "z2 = copy.copy(z1)  # initialize the elevations with the steady state\n",
    "# topography produced for question 1\n",
    "z2 = mg2.add_field(\"topographic__elevation\", z2, at=\"node\")\n",
    "\n",
    "# K_sp value for base landscape is 1e-5\n",
    "K_sp2 = 1e-5  # units vary depending on m_sp and n_sp\n",
    "m_sp2 = 0.5  # exponent on drainage area in stream power equation\n",
    "n_sp2 = 1.0  # exponent on slope in stream power equation\n",
    "\n",
    "frr2 = FlowAccumulator(mg2, flow_director=\"FlowDirectorD8\")  # intializing flow routing\n",
    "spr2 = StreamPowerEroder(\n",
    "    mg2, K_sp=K_sp2, m_sp=m_sp2, n_sp=n_sp2, threshold_sp=0.0\n",
    ")  # initializing stream power incision\n",
    "\n",
    "theta2 = m_sp2 / n_sp2\n",
    "# initialize the component that will calculate channel steepness\n",
    "sf2 = SteepnessFinder(mg2, reference_concavity=theta2, min_drainage_area=1000.0)\n",
    "# initialize the component that will calculate the chi index\n",
    "cf2 = ChiFinder(\n",
    "    mg2, min_drainage_area=1000.0, reference_concavity=theta2, use_true_dx=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 15\n",
    "\n",
    "tmax = 1e5  # time for the model to run [yr] (Original value was 5E5 yr)\n",
    "dt = 500  # time step [yr] (Original value was 500 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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 16\n",
    "\n",
    "#  uplift_rate [m/yr] (value was 0.0001 m/yr for base landscape)\n",
    "uplift_rate = np.ones(mg2.number_of_nodes) * 0.0001\n",
    "\n",
    "## If you want to add a one-time event that uplifts only part of the\n",
    "## landscape, uncomment the 3 lines below\n",
    "# fault_location = 4000  # [m]\n",
    "# uplift_amount = 10 # [m]\n",
    "# z2[np.nonzero(mg2.node_y>fault_location)] += uplift_amount\n",
    "\n",
    "## IMPORTANT! To use the below fault generator, comment the one-time\n",
    "## uplift event above if it isn't already commented out.\n",
    "## Code below creates a fault horizontally across the grid.\n",
    "## Uplift rates are greater where y values  > fault location.\n",
    "## To use, uncomment the 5 code lines below and edit to your values\n",
    "# fault_location = 4000  # [m]\n",
    "# low_uplift_rate = 0.0001 # [m/yr]\n",
    "# high_uplift_rate = 0.0004 # [m/yr]\n",
    "# uplift_rate[np.nonzero(mg2.node_y<fault_location)] = low_uplift_rate\n",
    "# uplift_rate[np.nonzero(mg2.node_y>fault_location)] = high_uplift_rate\n",
    "\n",
    "## IMPORTANT! To use below rock uplift gradient, comment the two\n",
    "## uplift options above if they aren't already commented out.\n",
    "## If you want a linear gradient in uplift rate\n",
    "## (increasing uplift into the range),\n",
    "## uncomment the 4 code lines below and edit to your values.\n",
    "# low_uplift_rate = 0.0001 # [m/yr]\n",
    "# high_uplift_rate = 0.0004 # [m/yr]\n",
    "## below is uplift gradient per node row index, NOT row value in meters\n",
    "# uplift_rate_gradient = (high_uplift_rate - low_uplift_rate)/(number_of_rows-3)\n",
    "# uplift_rate = low_uplift_rate + ((mg2.node_y / dxy)-1) * uplift_rate_gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 17\n",
    "\n",
    "for ti in t:\n",
    "    z2[mg2.core_nodes] += uplift_rate[mg2.core_nodes] * dt  # uplift the landscape\n",
    "    frr2.run_one_step()  # route flow\n",
    "    spr2.run_one_step(dt)  # fluvial incision\n",
    "    total_time += dt  # update time keeper\n",
    "    print(total_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 18\n",
    "# Plot topography\n",
    "\n",
    "plt.figure(8)\n",
    "imshow_grid(\n",
    "    mg2, \"topographic__elevation\", grid_units=(\"m\", \"m\"), var_name=\"Elevation (m)\"\n",
    ")\n",
    "plt.title(f\"$K_{{sp}}$={K_sp2}; $time$={total_time} yr; $dx$={dxy2} m\")\n",
    "\n",
    "max_elev = np.max(z2)\n",
    "print(\"Maximum elevation is \", np.max(z2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 19\n",
    "# Plot Channel Profiles and slope-area data along the channels\n",
    "\n",
    "prf2 = ChannelProfiler(\n",
    "    mg2,\n",
    "    number_of_watersheds=3,\n",
    "    main_channel_only=True,\n",
    "    minimum_channel_threshold=dxy**2,\n",
    ")\n",
    "prf2.run_one_step()\n",
    "\n",
    "# plot the elevation as a function of distance upstream\n",
    "plt.figure(9)\n",
    "title_text = f\"$K_{{sp}}$={K_sp2}; $time$={total_time} yr; $dx$={dxy} m\"\n",
    "prf2.plot_profiles(\n",
    "    xlabel=\"distance upstream (m)\", ylabel=\"elevation (m)\", title=title_text\n",
    ")\n",
    "\n",
    "# plot the location of the channels in map view\n",
    "plt.figure(10)\n",
    "prf2.plot_profiles_in_map_view()\n",
    "\n",
    "# slope-area data in just the profiled channels\n",
    "plt.figure(11)\n",
    "for i, outlet_id in enumerate(prf2.data_structure):\n",
    "    for j, segment_id in enumerate(prf2.data_structure[outlet_id]):\n",
    "        if j == 0:\n",
    "            label = \"channel {i}\".format(i=i + 1)\n",
    "        else:\n",
    "            label = \"_nolegend_\"\n",
    "        segment = prf2.data_structure[outlet_id][segment_id]\n",
    "        profile_ids = segment[\"ids\"]\n",
    "        color = segment[\"color\"]\n",
    "        plt.loglog(\n",
    "            mg2.at_node[\"drainage_area\"][profile_ids],\n",
    "            mg2.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_sp2}; $time$={total_time} yr; $dx$={dxy2} m\"\n",
    "plt.title(title_text)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 20\n",
    "# Chi Plots\n",
    "\n",
    "# calculate the chi index\n",
    "cf2.calculate_chi()\n",
    "\n",
    "# chi-elevation plots in the profiled channels\n",
    "plt.figure(12)\n",
    "for i, outlet_id in enumerate(prf2.data_structure):\n",
    "    for j, segment_id in enumerate(prf2.data_structure[outlet_id]):\n",
    "        if j == 0:\n",
    "            label = \"channel {i}\".format(i=i + 1)\n",
    "        else:\n",
    "            label = \"_nolegend_\"\n",
    "        segment = prf2.data_structure[outlet_id][segment_id]\n",
    "        profile_ids = segment[\"ids\"]\n",
    "        color = segment[\"color\"]\n",
    "        plt.plot(\n",
    "            mg2.at_node[\"channel__chi_index\"][profile_ids],\n",
    "            mg2.at_node[\"topographic__elevation\"][profile_ids],\n",
    "            color=color,\n",
    "            label=label,\n",
    "        )\n",
    "\n",
    "plt.xlabel(\"chi index (m)\")\n",
    "plt.ylabel(\"elevation (m)\")\n",
    "plt.legend(loc=\"lower right\")\n",
    "title_text = (\n",
    "    f\"$K_{{sp}}$={K_sp2}; $time$={total_time} yr; $dx$={dxy2} m; concavity={theta2}\"\n",
    ")\n",
    "plt.title(title_text)\n",
    "\n",
    "# chi map\n",
    "plt.figure(13)\n",
    "imshow_grid(\n",
    "    mg2,\n",
    "    \"channel__chi_index\",\n",
    "    grid_units=(\"m\", \"m\"),\n",
    "    var_name=\"Chi index (m)\",\n",
    "    cmap=\"jet\",\n",
    ")\n",
    "plt.title(\n",
    "    f\"$K_{{sp}}$={K_sp2}; $time$={total_time} yr; $dx$={dxy2} m; concavity={theta2}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code Block 21\n",
    "# Plot channel steepness along profiles and across the landscape\n",
    "\n",
    "# calculate channel steepness\n",
    "sf2.calculate_steepnesses()\n",
    "\n",
    "# plots of steepnes vs. distance upstream in the profiled channels\n",
    "plt.figure(14)\n",
    "\n",
    "for i, outlet_id in enumerate(prf2.data_structure):\n",
    "    for j, segment_id in enumerate(prf2.data_structure[outlet_id]):\n",
    "        if j == 0:\n",
    "            label = \"channel {i}\".format(i=i + 1)\n",
    "        else:\n",
    "            label = \"_nolegend_\"\n",
    "        segment = prf2.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",
    "            mg2.at_node[\"channel__steepness_index\"][profile_ids],\n",
    "            \"x\",\n",
    "            color=color,\n",
    "            label=label,\n",
    "        )\n",
    "\n",
    "plt.xlabel(\"distance upstream (m)\")\n",
    "plt.ylabel(\"steepness index\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.title(\n",
    "    f\"$K_{{sp}}$={K_sp2}; $time$={total_time} yr; $dx$={dxy2} m; concavity={theta2}\"\n",
    ")\n",
    "\n",
    "# channel steepness map\n",
    "plt.figure(15)\n",
    "imshow_grid(\n",
    "    mg2,\n",
    "    \"channel__steepness_index\",\n",
    "    grid_units=(\"m\", \"m\"),\n",
    "    var_name=\"Steepness index \",\n",
    "    cmap=\"jet\",\n",
    ")\n",
    "plt.title(\n",
    "    f\"$K_{{sp}}$={K_sp2}; $time$={total_time} yr; $dx$={dxy2} m; concavity={theta2}\"\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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}