From 83a6833b4945536373323ee1cdb54a9a7613b32f Mon Sep 17 00:00:00 2001 From: DavidLitwin <45630470+DavidLitwin@users.noreply.github.com> Date: Sat, 23 Nov 2019 22:28:03 -0700 Subject: [PATCH] Add tutorial for gdp (#80) --- groundwater/groundwater_flow.ipynb | 424 +++++++++++++++++++++++++++++ 1 file changed, 424 insertions(+) create mode 100644 groundwater/groundwater_flow.ipynb diff --git a/groundwater/groundwater_flow.ipynb b/groundwater/groundwater_flow.ipynb new file mode 100644 index 0000000..1685f6d --- /dev/null +++ b/groundwater/groundwater_flow.ipynb @@ -0,0 +1,424 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Modeling groundwater flow in a conceptual catchment\n", + "\n", + "
\n", + " For instructions on how to run an interactive iPython notebook, click here: https://github.com/landlab/tutorials/blob/release/README.md
\n", + "For more Landlab tutorials, click here: https://github.com/landlab/landlab/wiki/Tutorials\n", + "
\n", + "\n", + "This tutorial demonstrates how the GroundwaterDupuitPercolator can be used to model groundwater flow and seepage (groundwater return flow). It is recommended to read the documentation for the component before starting this tutorial to be familiar with the mechanics of the model.\n", + "\n", + "In this tutorial you will:\n", + "* Create a raster grid on which to run the model\n", + "* Simulate constant recharge and check that the component conserves mass\n", + "* Simulate recharge from storm events, check conservation of mass, and look at the outflow hydrograph" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A bit of magic so that we can plot within this notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from __future__ import print_function\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Import libraries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from landlab import RasterModelGrid, FIXED_VALUE_BOUNDARY, CLOSED_BOUNDARY, imshow_grid\n", + "from landlab.components import GroundwaterDupuitPercolator, FlowAccumulator\n", + "from landlab.components.uniform_precip import PrecipitationDistribution" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create a RasterModelGrid\n", + "\n", + "Here you will make the grid on which we will run the model. You will create three fields: topographic elevation, aquifer base elevation, and initial water table elevation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "boundaries = {'top': 'closed','bottom': 'closed','right':'closed','left':'closed'}\n", + "grid = RasterModelGrid((51, 51), spacing=10.0,bc=boundaries)\n", + "grid.status_at_node[1] = FIXED_VALUE_BOUNDARY \n", + "\n", + "elev = grid.add_zeros('node', 'topographic__elevation')\n", + "elev[:] = (0.001*grid.x_of_node**2 + 0.001*grid.y_of_node**2)+2\n", + "\n", + "base = grid.add_zeros('node', 'aquifer_base__elevation')\n", + "base[:] = elev - 2\n", + "\n", + "wt = grid.add_zeros('node', 'water_table__elevation')\n", + "wt[:] = elev" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "imshow_grid(grid,'topographic__elevation')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The grid is square with dimensions 500x500m. The surface elevation and aquifer base have the same concave parabolic shape, with thickness 2m between them. The aquifer is initially fully saturated (water table at the surface). Water is only allowed to exit the domain through a single node in the the lower left corner. All other boundaries are closed. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulate constant groundwater recharge\n", + "\n", + "Now initialize the model components. In addition to the grid, the GroundwaterDupuitPercolator takes four optional arguments: hydraulic conductivity, porosity, recharge rate, and a regularization factor that smooths the transition between subsurface and surface flow as the water table approaches the ground surface. The greater the value, the smoother the transition.\n", + "\n", + "You will also initialize a FlowAccumulator in order to use an included method to calculate the surface water discharge out of the domain. The runoff rate used by the FlowAccumulator is the surface water specific discharge from the groundwater model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "K = 0.01 # hydraulic conductivity, (m/s)\n", + "R = 1e-7 # recharge rate, (m/s)\n", + "n = 0.2 # porosity, (-)\n", + "gdp = GroundwaterDupuitPercolator(grid, hydraulic_conductivity=K, porosity = n, recharge_rate=R,regularization_f=0.01)\n", + "fa = FlowAccumulator(grid, surface='topographic__elevation',\n", + " flow_director='FlowDirectorSteepest', runoff_rate='surface_water__specific_discharge')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, run the model forward in time, and track the fluxes leaving the domain." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "N = 1000\n", + "dt = 1E2\n", + "\n", + "recharge_flux = np.zeros(N)\n", + "gw_flux = np.zeros(N)\n", + "sw_flux = np.zeros(N)\n", + "storage = np.zeros(N)\n", + "\n", + "for i in range(N):\n", + " gdp.run_one_step(dt)\n", + " \n", + " fa.run_one_step()\n", + " \n", + " recharge_flux[i] = gdp.calc_recharge_flux_in()\n", + " gw_flux[i] = gdp.calc_gw_flux_out()\n", + " sw_flux[i] = gdp.calc_sw_flux_out()\n", + " storage[i] = gdp.calc_total_storage()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now visualize some results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "imshow_grid(grid,(wt-base)/(elev-base),cmap='Blues')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above shows how saturated the aquifer is. Note that it is most saturated at the lowest area of the domain, nearest the outlet.\n", + "\n", + "Now look at the mass balance by ploting cumulative fluxes. The cumulative recharge in should be equal to cumulative fluxes out (groundwater and surface water) plus the change in storage from the initial condition." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t = np.arange(0,N*dt,dt)\n", + "\n", + "plt.figure(figsize=(8,6))\n", + "plt.plot(t/3600,np.cumsum(gw_flux)*dt+np.cumsum(sw_flux)*dt+storage-storage[0],\n", + " 'b-',linewidth=3, alpha=0.5,label='Total Fluxes + Storage' )\n", + "plt.plot(t/3600,np.cumsum(recharge_flux)*dt,'k:',label='recharge flux')\n", + "plt.plot(t/3600,np.cumsum(gw_flux)*dt,'b:',label='groundwater flux')\n", + "plt.plot(t/3600,np.cumsum(sw_flux)*dt,'g:',label='surface water flux')\n", + "plt.plot(t/3600,storage-storage[0], 'r:', label='storage')\n", + "plt.ylabel('Cumulative Volume $[m^3]$')\n", + "plt.xlabel('Time [h]')\n", + "plt.legend(frameon=False)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The thick blue line (cumulative fluxes plus storage) matches the black cumulative recharge flux line, which indicates that the model has conserved mass. Because the initial domain was fully saturated, the primary feature that shows up in this mass balance is the loss of that initial water. It will be easier to see what is going on here in the second example. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulate time-varying recharge\n", + "\n", + "Lastly, simulate time-varying recharge, look at the mass balance, and the outflow hydrograph. This will use the same grid and groundwater model instance as above, taking the final condition of the previous model run as the new initial condition here. This time the adaptive timestep solver will be used to make sure the model remains stable.\n", + "\n", + "First, we need a distribution of recharge events. We will use landlab's precipitation distribution tool to create a lists paired recharge events and intensities." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#generate storm timeseries\n", + "T = 10*24*3600 #sec\n", + "Tr = 1*3600 #sec\n", + "Td = 24*3600 #sec\n", + "dt = 1e3 #sec\n", + "p = 1e-3 #m\n", + "\n", + "precip = PrecipitationDistribution(mean_storm_duration=Tr, mean_interstorm_duration=Td, \n", + " mean_storm_depth=p, total_t=T, delta_t=dt)\n", + "durations = []\n", + "intensities = []\n", + "precip.seed_generator(seedval=1)\n", + "for (interval_duration, rainfall_rate_in_interval) in (\n", + " precip.yield_storm_interstorm_duration_intensity(subdivide_interstorms=True)\n", + "):\n", + " durations.append(interval_duration)\n", + " intensities.append(rainfall_rate_in_interval)\n", + "N = len(durations) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next run the model forward with the run_with_adaptive_time_step_solver. This method is the same as run_one_step, except that it subdivides the provided timestep (event or inter-event duration in this case) in order to meet a Courant-type stability criterion. The argument courant_coefficient indicates how large the maximum allowed timestep is relative to the Courant limit. Values close to 0.1 are recommended for best results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "recharge_flux = np.zeros(N)\n", + "gw_flux = np.zeros(N)\n", + "sw_flux = np.zeros(N)\n", + "storage = np.zeros(N)\n", + "num_substeps = np.zeros(N)\n", + "\n", + "for i in range(N):\n", + " gdp.recharge = intensities[i]*np.ones_like(gdp.recharge)\n", + " gdp.run_with_adaptive_time_step_solver(durations[i], courant_coefficient=0.2)\n", + "\n", + " num_substeps[i] = gdp.number_of_substeps\n", + " \n", + " fa.run_one_step()\n", + "\n", + " recharge_flux[i] = gdp.calc_recharge_flux_in()\n", + " gw_flux[i] = gdp.calc_gw_flux_out()\n", + " sw_flux[i] = gdp.calc_sw_flux_out()\n", + " storage[i] = gdp.calc_total_storage()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Again, visualize the mass balance:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t = np.cumsum(durations)\n", + "\n", + "plt.figure()\n", + "plt.plot(t/3600,np.cumsum(gw_flux*durations)+np.cumsum(sw_flux*durations)+storage-storage[0],\n", + " 'b-',linewidth=3, alpha=0.5,label='Total Fluxes + Storage' )\n", + "plt.plot(t/3600,np.cumsum(recharge_flux*durations)-recharge_flux[0]*durations[0],'k:',label='recharge flux')\n", + "plt.plot(t/3600,np.cumsum(gw_flux*durations),'b:',label='groundwater flux')\n", + "plt.plot(t/3600,np.cumsum(sw_flux*durations),'g:',label='surface water flux')\n", + "plt.plot(t/3600,storage-storage[0], 'r:', label='storage')\n", + "plt.ylabel('Cumulative Volume $[m^3]$')\n", + "plt.xlabel('Time [h]')\n", + "plt.legend(frameon=False)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize numer of substeps that the model took for stability:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(num_substeps,'.')\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel('Numer of Substeps')\n", + "plt.yticks([1,5,10,15,20])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "max(num_substeps)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The method has subdivided the timestep up to 18 times in order to meet the stability criterion. This is dependent on a number of factors, including the Courant coefficient, the hydraulic conductivity, and hydraulic gradient.\n", + "\n", + "Now look at the timeseries of recharge in and groundwater and surface water leaving the domain at the open node:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(figsize=(8,6))\n", + "ax.plot(t/(3600*24),sw_flux, label='Surface water flux')\n", + "ax.plot(t/(3600*24),gw_flux, label='Groundwater flux')\n", + "ax.set_ylim((0,0.04))\n", + "ax.set_ylabel('Flux out $[m^3/s]$')\n", + "ax.set_xlabel('Time [d]')\n", + "ax.legend(frameon=False,loc=7)\n", + "ax1 = ax.twinx()\n", + "ax1.plot(t/(3600*24),recharge_flux,'0.6')\n", + "ax1.set_ylim((1.2,0))\n", + "ax1.set_ylabel('Recharge flux in $[m^3/s]$')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The relationship between maximum flux that can be passed through the subsurface and the occurrence of groundwater seepage is clear from this figure." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Click here for more Landlab tutorials" + ] + }, + { + "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.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}