example_1d.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples of using GRLP\n",
    "\n",
    "GRLP is a model that solves for the long-profile evolution of a fluvial channel using a built-in semi-implicit solver. \n",
    "In particular, it simulates the evolution of a gravel-bed river that shapes its bed by eroding, transporting, and depositing sediments. \n",
    "The model takes a few basic inputs, and since the numerical solver has been set up implicitly, the model can take long timesteps while remaining stable. \n",
    "As its basic inputs, the model takes a sediment and water input, a spatial domain, an uplift or subsidence rate (or equivalently, other local source or sink of material), and a floodplain width, plus a slew of appropriate constants. \n",
    "The model can also replace some of these with a set of power-law scaling relationships between downstream distance and (1) valley width, (2) drainage area, and (3) discharge.\n",
    "If this approach is taken, it allows an analytical solution, assuming that these variables *actually* follow a power-law downstream. \n",
    "However, for the numerical solution, these can be an arbitrary functions up and downstream, instead of power laws. \n",
    "\n",
    "[Wickert & Schildgen (2019, ESurf)](https://esurf.copernicus.org/articles/7/17/2019/) describe the theory behind this model and the equations that are implemented numerically here.\n",
    "\n",
    "First, we will present an example where we compare the numerical solution with power-law functions for important variables to the corresponding analytical solution, which will demonstrate some basic outcomes and methods of interacting with the model. \n",
    "The very first thing to do is load the relevant libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import numerical and plotting libraries\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from copy import deepcopy\n",
    "\n",
    "# Import the GRLP module\n",
    "import grlp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%autoreload"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GRLP is an object-oriented model, so next, we instantiate `grlp` model object,('`lp`') which creates the data structure, and associates methods with `lp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lp = grlp.LongProfile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To provide `lp` with the necessary starting values, we set some constants and create the domain, using some handy built-in methods that set a suite of reasonable constants for Earth rivers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the x domain\n",
    "lp.set_x(dx=1000, nx=90, x0=10000)\n",
    "\n",
    "# S0 is the upstream-end slope that determines the sediment input to the catchment\n",
    "S0 = 1.5E-2\n",
    "\n",
    "# Set up a starting set of channel-bed elevations (z) on a uniform slope (S0)\n",
    "lp.set_z(S0=-S0, z1=0)\n",
    "\n",
    "# Intermittency: What fraction of the total time is the river experiencing a\n",
    "# geomorphically-effective flood? This assumes a binary on--off state, common \n",
    "# for gravel-bed rivers with floodplains (see Blom et al., 2017:\n",
    "# https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017JF004213).\n",
    "# For an evaluation of the intermittency in general and its variability across\n",
    "# a range of fluvial systems, see Hayden et al. (2021):\n",
    "# https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL092598\n",
    "lp.set_intermittency(0.80)\n",
    "\n",
    "# Set up the number of iterations in semi-implicit solver; defaults to 3.\n",
    "lp.set_niter()\n",
    "\n",
    "# Utility functions to assign reasonable constants \n",
    "# (defined by Wickert & Schildgen, 2019)\n",
    "lp.basic_constants()\n",
    "lp.bedload_lumped_constants()\n",
    "lp.set_hydrologic_constants()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As noted in the comments, the \"intermittency\" is a factor that simplifies a river channel's behavior into times of \"shaping its channel\" and times of performing no morphological work.\n",
    "* For a discussion of the intermittency and its appropriateness as an approximation to the full hydrograph for the purposes of fluvial geomorphology, see [Blom et al., 2017, GRL](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017JF004213)\n",
    "* For a catalog of intermittency values across a wide range of fluvial systems and a broader discussion of its significance, see [Hayden et al., 2021, JGR - Earth Surface](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL092598) \n",
    "\n",
    "Pursuant to our assumption that valley width, discharge, and drainage area vary as power-law functions downstream, we will use handy setting functions that provide these values for the whole `x` domain.\n",
    "* Distance downstream is known to scale with drainage area through Hack's Law (countless publications)\n",
    "* The drainage-area--discharge relationship is well established, see:\n",
    "    * 50 streams in PA and NJ, USA: [Aron and Miller, 1978](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1752-1688.1978.tb02169.x)\n",
    "    * All USGS gages, including PR: [O'Connor and Costa, 2004](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2003WR002247)\n",
    "    * Theoretical approach: [Sólyom and Tucker, 2004](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2003JF000032)\n",
    "* The valley width -- drainage-area scaling is described in detail by [Langston & Temme (2019, GRL)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL085164), and they note a common exponent of *~0.4* in the power-law relation. Other relevant references include:\n",
    "    * [Tomkin et al. (2003)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2001JB000862): exponent 0.76\n",
    "    * [May et al. (2013)](https://pubs.geoscienceworld.org/gsa/geology/article/41/4/503/131211/Controls-on-valley-width-in-mountainous-landscapes): exponent 0.574\n",
    "    * [Schanz & Montgomery (2016)](https://www-sciencedirect-com.ezp3.lib.umn.edu/science/article/pii/S0169555X16300150): exponents 0.22, 0.34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up transfer functions between drainage area (A), discharge (Q), and\n",
    "# valley width (B).\n",
    "\n",
    "# drainage area: A = k_xA * x**P_xA\n",
    "lp.set_A(k_xA=1., P_xA = 7/4)\n",
    "\n",
    "# discharge: Q = k_xQ * x**P_xQ\n",
    "lp.set_Q(k_xQ = 1.43e-5, P_xQ = 7/4*0.7)\n",
    "\n",
    "# Valley width: B = k_xB * x**P_xB\n",
    "lp.set_B(k_xB = 4, P_xB = 0.4)\n",
    "\n",
    "# Set the uplift rate [m/s]; positive upwards\n",
    "lp.set_uplift_rate(0)\n",
    "\n",
    "# Set the base level \n",
    "# This is currently redundant, but will be useful later. \n",
    "# the set_z function already assumes that z_bl starts at 0.\n",
    "lp.set_z_bl(0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the last part, we will set the input sediment discharge, which arrives in the upstream (in our case left-most) model node. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input sediment discharge: this is set based on your defined S0, above.\n",
    "# (this ficticious boundary-condition slope is the transport slope for the\n",
    "#  amount of sediment being supplied)\n",
    "\n",
    "Qs0 = lp.k_Qs * lp.Q[0] * (S0)**(7/6.)\n",
    "\n",
    "lp.set_Qs_input_upstream(Qs0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what the initial model domain looks like, just making a plot of the long-profile. You'll notice that it's just a flat line, because of `lp.set_z(S0=-S0, z1 = 0)` above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(lp.x/1000., lp.z, '0.6', linewidth=6, label='Numerical')\n",
    "plt.xlim(lp.x.min()/1000, lp.x.max()/1000)\n",
    "plt.ylim(0, lp.z.max())\n",
    "plt.xlabel('Downstream distance [km]', fontsize=26)\n",
    "plt.ylabel('Elevation [m]', fontsize=26)\n",
    "plt.tick_params(axis='both', which='major', labelsize=16)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can evolve the river profile, and calculate the analytical and numerical solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numerical: (number of time steps, length of time step [s])\n",
    "# note that this is a very long time step (31 Ma)\n",
    "lp.evolve_threshold_width_river(1, 1E15)\n",
    "# Analytical, no time steps, just solve for equilibrium profile.\n",
    "zanalytical = lp.analytical_threshold_width() # suppress printing returned data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the analytical solution, and the numerical solution together. They match!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(lp.x/1000., lp.z, '0.6', linewidth=6, label='Numerical')\n",
    "plt.plot(lp.x/1000., lp.zanalytical, 'k', linewidth=2, label='Analytical')\n",
    "plt.xlim(lp.x.min()/1000, lp.x.max()/1000)\n",
    "plt.ylim(0, lp.z.max())\n",
    "plt.xlabel('Downstream distance [km]', fontsize=26)\n",
    "plt.ylabel('Elevation [m]', fontsize=26)\n",
    "plt.tick_params(axis='both', which='major', labelsize=16)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what is the geometry of this channel? The model object `lp` has a lot of attributes and methods. You can query them (like any other python object), using `dir()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(lp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can see that the object has an attribute `B` for valley width, and `b` for river width. If we query `B`..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lp.B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but if we query `b` we get..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(lp.b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "because this information about channels (width, depth, etc.) is computed after the fact, and it requires knowledge of the grain size and sediment discharge everywhere, which either cancels out of our is calculated within the equations for long-profile evolution.\n",
    "The `LongProfile` class contains methods to calculate these, but because we follow the near-threshold channel-geometry closure of [Parker (1978)](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/selfformed-straight-rivers-with-equilibrium-banks-and-mobile-bed-part-2-the-gravel-river/3AD9322C1939528ED73D409654E35E22) (Gravel on banks exactly at the threshold of motion at a bankfull, or equivalently in this framework, geomorphologically effective flow), computing the channel width and flow depth (= channel depth) also requires us to know the grain size.\n",
    "Let's say that this input grain size is 10 cm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lp.compute_Q_s()\n",
    "lp.D = 0.1\n",
    "lp.compute_channel_width()\n",
    "lp.compute_flow_depth()\n",
    "print(lp.b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plot these things too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 1, figsize=(12,10))\n",
    "axes[0].plot(lp.x/1000., lp.b, 'xkcd:blue', linewidth=4)\n",
    "axes[1].plot(lp.x/1000., lp.B, 'xkcd:red', linewidth=4)\n",
    "axes[2].plot(lp.x/1000., lp.h, 'xkcd:orange', linewidth=4)\n",
    "axes[2].set_xlabel('Downstream distance [km]', fontsize=16)\n",
    "axes[0].set_ylabel('Channel Width [m]', fontsize=16)\n",
    "axes[1].set_ylabel('Valley Width [m]', fontsize=16)\n",
    "axes[2].set_ylabel('Flow Depth [m]', fontsize=16)\n",
    "plt.tick_params(axis='both', which='major', labelsize=16)\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the basic usage of the model, which is for a given long-profile of a single stream, to calculate the equilibrium profile. \n",
    "Now, there are plenty of situations where the analytical profile would not be expected to match the numerical profile, because there are several simplifying assumptions applied to the analytical solution to render the math tractable.\n",
    "1. Uplift is assumed to match river incision exactly for the analytical solution.\n",
    "2. Discharge, Drainage Area, and Valley Width must be specified functions of downstream distance. However, the model does not specify _why_ that should be, only that it is. So if the downstream function of these variables is anything but a power law, the solutions should be different. \n",
    "\n",
    "So, for fun, let's violate #1 and see how it goes.\n",
    "First, we will make a copy of the existing object using `deepcopy`, then we will change the boundary conditions, and run each of them forward. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lp_uplift = deepcopy(lp)\n",
    "lp_uplift.set_uplift_rate(1e-10) # in [m/s]\n",
    "\n",
    "lp_uplift.evolve_threshold_width_river(1, 1E15) # evolve for the same amount of time as lp\n",
    "_ = lp_uplift.analytical_threshold_width()\n",
    "\n",
    "lp.evolve_threshold_width_river(1, 1E15) # this is unnecessary, but it's not expensive\n",
    "_ = lp.analytical_threshold_width()\n",
    "\n",
    "lp.compute_Q_s()\n",
    "lp.compute_channel_width()\n",
    "lp.compute_flow_depth()\n",
    "\n",
    "lp_uplift.compute_Q_s()\n",
    "lp_uplift.compute_channel_width()\n",
    "lp_uplift.compute_flow_depth()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(12,10))\n",
    "axes[0].plot(lp.x/1000., lp.z, 'xkcd:green', linewidth=8, label = 'Numerical', alpha = 0.3)\n",
    "axes[0].plot(lp.x/1000., lp.zanalytical, 'xkcd:green', linewidth=2, label = 'Analytical', alpha = 1)\n",
    "axes[1].plot(lp.x/1000., lp_uplift.z, 'xkcd:purple', linewidth=8, label = 'Numerical', alpha = 0.3)\n",
    "axes[1].plot(lp.x/1000., lp_uplift.zanalytical, 'xkcd:purple', linewidth=2, label = 'Analytical', alpha = 1)\n",
    "axes[1].set_xlabel('Downstream distance [km]', fontsize=16)\n",
    "axes[0].set_title('Uplift = Incision', fontsize = 20)\n",
    "axes[1].set_title('Uplift ≠ Incision', fontsize = 20)\n",
    "for a in axes:\n",
    "    a.set_ylabel('Elevation [m]', fontsize=16)\n",
    "    a.legend()\n",
    "    a.tick_params(axis='both', which='major', labelsize=10)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Violating #2 does other interesting stuff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Screw with transfer functions between drainage area (A), discharge (Q), and valley width (B).\n",
    "# Instead of power laws, these are now either constants, or a linear function \n",
    "# (drainage area has to increase downstream --- that's just geometry)\n",
    "\n",
    "lp_constantQBA = deepcopy(lp)\n",
    "\n",
    "lp_constantQBA.set_Q(Q = np.repeat(lp.Q[0], len(lp.x)))\n",
    "lp_constantQBA.set_B(B = np.repeat(lp.B[0], len(lp.x)))\n",
    "lp_constantQBA.set_A(A = lp.x * np.repeat(lp.A[0], len(lp.x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lp_constantQBA.evolve_threshold_width_river(1, 1E15) # evolve for the same amount of time as lp\n",
    "_ = lp_constantQBA.analytical_threshold_width()\n",
    "\n",
    "lp_constantQBA.compute_Q_s()\n",
    "lp_constantQBA.compute_channel_width()\n",
    "lp_constantQBA.compute_flow_depth()\n",
    "\n",
    "lp.evolve_threshold_width_river(1, 1E15) # this is unnecessary, but it's not expensive\n",
    "_ = lp.analytical_threshold_width()\n",
    "\n",
    "lp.compute_Q_s()\n",
    "lp.compute_channel_width()\n",
    "lp.compute_flow_depth()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(12,10))\n",
    "axes[0].plot(lp.x/1000., lp.z, 'xkcd:gold', linewidth=8, label = 'Numerical', alpha = 0.3)\n",
    "axes[0].plot(lp.x/1000., lp.zanalytical, 'xkcd:gold', linewidth=2, label = 'Analytical', alpha = 1)\n",
    "axes[1].plot(lp.x/1000., lp_constantQBA.z, 'xkcd:dark grey', linewidth=8, label = 'Numerical', alpha = 0.3)\n",
    "axes[1].plot(lp.x/1000., lp_constantQBA.zanalytical, 'xkcd:dark grey', linewidth=2, label = 'Analytical', alpha = 1)\n",
    "axes[1].set_xlabel('Downstream distance [km]', fontsize=16)\n",
    "axes[0].set_title('QBA vary as power-law', fontsize = 20)\n",
    "axes[1].set_title('QBA are all constants', fontsize = 20)\n",
    "for a in axes:\n",
    "    a.set_ylabel('Elevation [m]', fontsize=16)\n",
    "    a.legend()\n",
    "    a.tick_params(axis='both', which='major', labelsize=10)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The river can also evolve through time! Till now, this tutorial uses _really_ big timesteps (300 My), but the model can also resolve shorter time. Here's a quick example, where an sudden increase in uplift rate is applied to the upstream end of our original channel profile `lp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lp_increase = deepcopy(lp)\n",
    "fig = plt.figure(figsize=(12,5))\n",
    "ax1 = fig.add_subplot(1,1,1)\n",
    "plt.xlabel('Downstream distance [km]', fontsize=14, fontweight='bold')\n",
    "plt.ylabel('Elevation [m]', fontsize=14, fontweight='bold')\n",
    "plt.tight_layout()\n",
    "\n",
    "# Initial condition\n",
    "ax1.plot(lp.x/1000., lp_increase.z, color='.5', linewidth=3, label = 'Start')\n",
    "\n",
    "# Transient\n",
    "U = 1E-3\n",
    "lp_increase.set_uplift_rate(U/3.15E7)\n",
    "for i in range(5):\n",
    "    lp_increase.evolve_threshold_width_river(1, 1E12) # five timesteps of 31 ka\n",
    "    ax1.plot(lp_increase.x/1000., lp_increase.z, color='.5', linewidth=1)\n",
    "\n",
    "# New equilibrium\n",
    "lp_increase.evolve_threshold_width_river(1, 1E14) # final state 3 Ma later\n",
    "ax1.plot(lp.x/1000., lp_increase.z, color='0', linewidth=3, label = 'Finish')\n",
    "\n",
    "ax1.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}