chapter_2_root_finding/Green_Ampt_infiltration/GreenAmpt_infiltration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Green-Ampt model for infiltraion and runoff\n",
    "\n",
    "by Xiaofeng Liu, Ph.D., P.E.\n",
    "Associate Professor\n",
    "\n",
    "Department of Civil and Environmental Engineering\n",
    "\n",
    "Institute of Computational and Data Sciences\n",
    "\n",
    "Penn State University\n",
    "\n",
    "223B Sackett Building, University Park, PA 16802\n",
    "\n",
    "Web: http://water.engr.psu.edu/liu/\n",
    "\n",
    "(with contribution from Dapeng Feng, who was the TA for CE 497 in Fall 2019)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "This notebook demonstrates the root finding methods on how to calculate the infiltration process based on the Green-Ampt model. During rainfall, water will infiltrate into the soil. As water moves down, a so-called \"wetting front\" marches downward too. If the rainfall is so intense and exceedst the infiltration capacity of soil, water will start to accumulate or \"pond\" on the surface and generate runoff. Based on this simplified phyiscal description, there are several models to capture the processes. The Green-Ampt model is one proposed in 1911 and still used today. \n",
    "\n",
    "Good references for the Green-Ampt model include the following:\n",
    "```\n",
    "- L. W. Mays (2010). Water Resources Engineering, 2nd edition, John Wiley & Sons, Inc.\n",
    "- V. T. Chow, D. R. Maidment, and L. W. Mays (1988). Applied Hydrology, McGraw-Hill, New York\n",
    "```\n",
    "\n",
    "In the following, the Green-Ampt model is briefly introduced. One quantity of interest is the cumulative infiltration $F(t)$ as a function of time $t$. As you will see, we need to solve a nonlinear algebraic equation for $F(t)$. The root of the cumulative infiltration equation can not be sloved directily. Therefore, this is a good case to use the root finding methods. In this example, the secant method in the Scipy package is used. It can be adapted to use other methods covered in class.   \n",
    "\n",
    "## The Green-Ampt method\n",
    "Figure 1 shows the variables defined in the Green-Ampt model and Figure 2 plots the infiltration rate and cumulative infiltration of a typical infiltration process. \n",
    "\n",
    "<img src=\"GA_definition.png\" width=\"400\"/>\n",
    "\n",
    "<h3 align=\"center\">Figure 1. Variables in the Green-Ampt infiltartion model (adapted from Chow et al. (1988)).</h3> \n",
    "\n",
    "<img src=\"GA_curves.png\" width=\"400\"/>\n",
    "\n",
    "<h3 align=\"center\">Figure 2. Rainfall infiltration rate and cumulative infiltration (adapted from Larry Mays (2010)).</h3> \n",
    "\n",
    "Applying Darcy's law to the model in Figure 1, we can get\n",
    "\\begin{equation}\n",
    "f = K \\left(\\frac{\\psi+L}{L}\\right)\n",
    "\\end{equation}\n",
    "where $f$ is the infiltration rate, $K$ is the hydraulic conductivity, and $\\psi$ is the wetting front soil suction head, and $L$ is the length of the wetted zone. This equation neglects the depth of surface ponding $h_0$, which is comparitively small to other terms. The infiltration rate is defined as $f = \\frac{dF}{dt}$. Then, the above equation can be expressed as a differential equation with one unknown variable, i.e., the cumulative infiltration $F$\n",
    "\\begin{equation}\n",
    "\\frac{dF}{dt} = K \\left(\\frac{\\psi+L}{L}\\right)\n",
    "\\end{equation}\n",
    "From mass balance, the infilted water is used to wet the soil between the surface and the wetting front (minus the initial moisture). This means $L=F/\\Delta \\theta$. Plugging this into the previous equation, we get\n",
    "\\begin{equation}\n",
    "\\frac{dF}{dt} = K \\left(\\frac{\\psi \\Delta \\theta+F}{F}\\right)\n",
    "\\end{equation}\n",
    "This is a simple first-order ordinary differential equation (ODE) subject to the initial condition of \n",
    "\\begin{equation}\n",
    "F(t=0) = 0\n",
    "\\end{equation}\n",
    "i.e. there is no infiltration at the beginning. One can use any numerical methods covered in the ODE chapter to solve this problem. Luckily, there is an analytical solution. Integrating the ODE with the use of the initial condition, one can get\n",
    "\\begin{equation}\n",
    "F(t) = K t + \\psi\\Delta\\theta\\ln\\left(1+\\frac{F(t)}{\\psi\\Delta\\theta}\\right)\n",
    "\\end{equation}\n",
    "Unluckily, there is no explicit solution for $F(t)$. To get $F(t)$, we have to solve the above equation using a root finding method. \n",
    "\n",
    "Before we show how to solve the problem using root finding method, other parameters in the problem have to be defined and specified. The effective saturation $S_e$ is defined as\n",
    "\\begin{equation}\n",
    "S_e = \\frac{\\theta-\\theta_r}{\\eta-\\theta_r}\n",
    "\\end{equation}\n",
    "where $\\theta$ is the soil moisture, $\\theta_r$ is the residual moisture content, and $\\eta$ is the porosity. $\\theta_e=\\eta-\\theta_r$ is called the effective porosity. For the initial condition, $\\theta=\\theta_i$, thus $\\theta_i-\\theta_r=S_e\\theta_e$. The change of the moisture $\\Delta\\theta$ can be derived as\n",
    "\\begin{equation}\n",
    "   \\Delta\\theta = \\eta-\\theta_i = \\eta - (S_e\\theta_e+\\theta_r) = (1-S_e) \\theta_e\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Green-Ampt method example\n",
    "\n",
    "This is the modified example similar to that in Mays (2010), which uses the Green-Ampt method to evaluate the infiltration rate and cumulative infiltration depth for a silty clay soil from 0 to 8 hr. The initial effective saturation $S_e$ is 0.15. For a silty clay soil, $\\theta_e=0.45$, $\\psi=30 cm$, and $K=0.04 cm/hr$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import optimize\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Green Ampt model F=K*t+psi*dtheta*ln(1+F/(psi*dtheta))\n",
    "def greenampt(F, t, K=K, psi=psi, dtheta=dtheta):\n",
    "    \n",
    "    y = -F+K*t+psi*dtheta*np.log(1+F/(psi*dtheta))\n",
    "    return y\n",
    "\n",
    "# SECANT method to find the root by calling functions in \n",
    "# the optimize package. Without passing the first order derivative,\n",
    "# the \"newton\" funciton uses the SECANT method.\n",
    "def getroot(t, finit):\n",
    "    fun = lambda F: greenampt(F, t) \n",
    "    Froot = optimize.newton(fun, finit)\n",
    "    return Froot\n",
    "\n",
    "# define parameters\n",
    "K = 0.04        # hydraulic conductivity in cm/hr\n",
    "theta_e = 0.45  # effective porosity: porosit - residual moisture content\n",
    "Se = 0.15       # initial effective saturation\n",
    "psi = 30        # wetting front soil suction head in cm \n",
    "\n",
    "dtheta = (1-Se)*theta_e # calculate delta moisture\n",
    "\n",
    "trange = np.arange(0.0, 8.1, 0.1) # time range in hr\n",
    "F = np.zeros(len(trange)) # cumulative infiltration\n",
    "\n",
    "for i in range(len(trange)):\n",
    "    t = trange[i]\n",
    "    fg=K*t # initial guess of F(t)=Kt\n",
    "    F[i] = getroot(t, fg) # calculate cumulative infiltration\n",
    "\n",
    "# calculate infiltration rate. In the denominator, 1e-9 is \n",
    "# added to avoid the division-by-zero error (because F(t=0)=0).\n",
    "f = K*(psi*dtheta/(F+1e-9)+1) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the simulated infiltration rate and cumulative inflitration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAELCAYAAAAhuwopAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZRcZZk/8O/Te6eTTqezkBAS6O4EQtiydAIooiQdNlGWLIjIUecHSXA5o3MgAY9zBo8zg0FmHLcDCYzjOIqSAAqoo2ZhRCRgOitCJEsn0UDWXtJZeu/n98f73qrblVq7b9W9Vf39nHNPVd26VfXQdOrpd3teUVUQERGlIs/vAIiIKPsweRARUcqYPIiIKGVMHkRElDImDyIiShmTBxERpSyQyUNE6kRkbZznm0Vks4isyGRcRERkSFDXeYjIWlWdF+O5OlVdl+mYiIjICGTLIwkVIlLtdxBERINVtiaPSgBNIrLS70CIiAajAr8D6A9VXQUAItIiIgtU9TnnORFZDGAxAJSVlc2cMmVK6HW7GnehV3sxZdSUyLckIiJr8+bNx1V1dLxrsi552OTQZBNGY+TzNrGsAoDa2lqtr68PPbfk5SV4fufzqF9WH/kyIiKyRORAomsC2W0lIgsA1Npb55wz+2o1gBYRqQMAd6sjkckjJ6OxrRHNbc2exktENNgEsuVhE8JzEefm2dsWAM5Mq5RmXE2qnAQA2NO0B7PGzxp4oEREg1QgWx7pMrlyMgBgd9NunyMhIspugyp51FTWQCDY07TH71CIiLLaoEoeJQUlOK/8PLY8iIgGaFAlD8AMmu9uZPIgIhqIwZc8Kiez24qIaIAGXfKYVDmJ03WJiAZo0CUPzrgiIhq4wZc8Rprkwa4rIqL+G3TJo3pENQTCQXMiogEYdMmjpKAEE4ZPYLcVEdEADLrkAZhBc3ZbERH136BMHpMrJ7PlQUQ0AIM2eTS1NaGprcnvUIiIslLak4eITBORO0RkjohMS/fnJcNdXZeIiFKXlpLsIlIFYAkABdAAoAnACAA1IvIJe36lqu5Px+cn4kzX3d24G7PHz/YjBCKirOZ58hCRuQBUVR9KdJ2IVKvqBq9jSMSZrsuWBxFR/wwoeYhIuaq2RpyuV9UTcZ4HAKjqehEZPpDP7y9O1yUiGpiUkoeIzIk4tRDA/e4TTuKwnhaRZbG6pyKuzSjOuCIi6r+kB8xF5EkA1wOY6TpqErxsJYCKiPe5N8UY04LVdYmI+i+VlscaVV3vPiEiz8W62FoJYK+IVMIMmgtM0nk6pSjTYFLlpNB03crSSr/DISLKKqlM1W20027LXefmJ3jNclW9QVVn2dvrASxKPUzvuWdcERFRalJpeWwAsAmAiIjCtCKqADwe6wWq+nyU05tSijBN3KXZrzzvSp+jISLKLqkkj4VRuq2mJ3pRMoPsfqiprEF5cTle2fcKPnX5p/wOh4goqySdPOzU2nIAdQAqAaxW1a3xXmMH2VsANLpOJxpkz4ii/CLccuEtePHdF7GydyUK8tKyXpKIKCcl/Y1pF/8tBLDXnnpKRFYmWOTXn0H2jJl/8Xw889YzePXAq5hTFdlAIiKiWFIZMK9Q1aWq+k173AlTciSe/gyyQ0TqRGRtnOcX2GuWJRl7VDdOuhGlBaV4/p1oQzNERBRLKn01zcmcE5Em2IH1iPNJDbIDgKquE5Hl0Z4TkQWua6pFpE5V1yX3n9DXkMIhuGnyTfj5X36O7978XeTJoCwyTESUslSSxwgReQDAFvt4Bsx4RmS31VkD627JDLInMAvAs/Z+g42jX8kDMF1XL+x8AW8cfAMfmPCBAYZGRDQ4JP2ntp12ux5mncYiAFtVNdpiv7hTcRMNsiehIuLxSPcDEVksIvUiUn/s2LGEb3bLhbegKL+IXVdERClIqZ9GVbfacY+ldvZVtFIj80TkQRG5N2KswystMLO9YsW4SlVrVbV29OjRCd+svLgc86rn4fmdz0NVvYyTiChnJey2EpHdMCVFtiA80woIj1/0aX04CwNtxdw77W2Dqr7gUcybEG59VAOIObCerDsuvgO/2v0rbD28FTPGzRjo2xER5byELQ9VnWzLqi+0JUacI26pEVU9oapPqerjALaKyH0i8kCURYNnsYPitc7guD231r7vcwCqRaTOPu73eIfj1otuRb7ks+uKiChJkumuGjtgXgezm+A6Vd2Wrs+qra3V+vr6pK6t+1EdDrYexM7P74SIJH4BEVGOEpHNqlob75pUSrJPc92vsjsBpryyzo6bfNO2SALzLT3/4vl4t/Fd7Dy+0+9QiIgCL5UB82rnjqrus9Nxq+NcH5WIXOB6n4HOvPLMbVNug0Cw+u3VfodCRBR4cQfM7WD3IgBLAAwXkYcRbi00wuzXEe/1FwBYCtNFBfvauTBrNQJl3LBxuGnyTfjun76LL131JVSURM4IJiIiR9yWh2vQuxbAUrsvR619/BDMtNl4lsLMhlptj2cBbPYg7rT41zn/iua2ZnzjtW/4HQoRUaClssK8XkQexdmL9OIVRlwbpTDigOpRpdMVY6/A3ZffjW+/+W18ftbnMWH4BL9DIiIKpFTGPB6CaT08B2CV64inwU7PvUNE5tgB9hX9CzUzvn7d19GrvXjk/x7xOxQiosBKJXmsszOl1gOosoPdiarqLgcwCmYPj5n2CMR+HrFcUHEBvjDrC/jh9h/i7aNv+x0OEVEgpVSeREScgoSzbZHEhQleslJVH3KVcf8mzOB7oH3lQ1/BsKJheGj9Q36HQkQUSKkURlxv9/CAqj4E4ASAxxK9rD/7efht5JCRePiah/HLXb/Eqwde9TscIqLASXqFuYiU2zIlyb953709Qvt5qOrkVAPtj1RWmEdq62rD5O9Oxqgho7Dx/21EaWGpx9EREQWTpyvMATztXuCXJKce1vXJ1MMKktLCUqz62CrsOLIDS3+1lBV3iYhcUkkeKxExTTdGSXa3JhF5wnX9NETfkTCQbp58Mx75yCP40fYf4Xt/+p7f4RARBUYq6zxWAtgrIpUAmmC6oGYioiR7FKFxEVXdZqfr7k8xTt989dqvYvOhzfiH3/0Drhh7Ba49/1q/QyIi8l0qLY/ltutpVgpdUFWqui/iXFbV/ciTPPzoth+hekQ1Fq5ZiIOtB/0OiYjId6luQxt5LuZe5ZaIyBN2keAdtgsr5i6AQTW8ZDh+cecvcKbrDG5/9na0tCeqykJElNvSWpLdJpxVAGbbY1WMfc8D7+LRF+OZO57B9sPbcd1/X4cjp474HRIRkW/SXpLdrkp/yB6BKcHeHx+76GN4+a6XsatxF675r2twoOWA3yEREfkibvIQkeF2+9h6ACtEZJPr+C3MwHnka+5I9sP7s5mU326YdAPW3rMWx88cxwd/8EG8c+wdv0MiIsq4VBYJzk1ijMO59j6YRYGrIxcW2tXmdwLQdHdhDWSRYCI7juzA9f9zPbp6u/Dj23+MmybflJbPISLKtGQWCQ5oD3MRuTdWArAbSS2GKYToLqC4CcDzUWZheS6dyQMA9jbtxe3P3o63jr6FL135JXyj7hsoLihO2+cREWXCgJOHiOwBsBd9y4sAPpQa6Y90Jw8AaO9ux4O/exDf2/Q9TBs7DT+d/1NMGTUlrZ9JRJROXiSPuF1VIjI9yIPgmUgejpfffRmfffGzONN1Bv/04X/Cl676ElshRJSVBlzbKtEYR5ATR6Z97KKPYcf9OzCvZh4eWv8QLn3iUrz87susiUVEOSnhVF2nnLotrR55PJHo9YPJucPOxYufeBG/ufs3KMgrwMd/9nHc9JOb8NaRt/wOjYjIU4mm6t6H8LaxT8PMkvqE64i7K6Cd6vukiGyyj6e7Fxvmqhsm3YAdS3fg36//d2w8uBGXP3k5bvvZbdj03ia/QyMi8kSiMY87VPUFe39+ZIkSEYlWu8r9/AMAngJQ7XRxicgcVd0QNyiRBQBaAMxQ1bM2nBKRZgANMFvjLo/1Ppkc84ilqa0J33nzO/jOm99Bc3szrq+5Hss+sAxzquZARBK/ARFRhnmxn0eNaw+PaFlmboLX71PVEzFeG5VNHFDVdQBaRKQuymULVXVmvMQRFJWllXjkI4/gwJcOYEXdCmw7vA11/1OHKd+fgm9t/Baa2s5aZ0lEFHgJ13mIyJMwZUiqcfa03ZmqOjLOa++zr2sEsAXADAAt8RYHisgKAM+q6habOM5qfdgEs0VVG+LFHoSWR6S2rjY8985zeHLzk3j9b6+jpKAEC6YuwN2X3Y266joU5KVSJZ+IyHvJtDwSflOp6lL7ZmdNyxWRuC0PVX1KRKbDjJUsArAmiVXqkSXboyWnSpiNplaq6pKImBbDLE7ExIkTE3xU5pUWluKeK+7BPVfcg+2Ht+PJ+ifxs7d/hh/v+DFGDxmNRZcswl2X3oWrzrsK+Xn5fodLRBTVQFeYp7yveRLvuQLAWlVdZ1se82J1T9lrN6nqc9GeD2LLI5qO7g78Zs9v8Myfn8FL776E9u52jCkbg1sm34Jbp9yKuuo6DCkc4neYRDRIeNLyiHjDyEKGCwHcH+f6aaq6zd6vgunC0gQD5psQbn1UA1gb8Z6LATTZhNGYSvxBVVxQjFun3Ipbp9yK1o5W/GrXr/DSrpfw3M7n8INtP0BJQQmuPf9azKueh3nV83DZOZchT1IpiExE5K1UCiM+CTMDyv2FPc/uKBjrNaHZWq5zMethua5ZBjNGUq2qq+y5tao6T0QqADgZMepsLEe2tDxi6ezpxKsHXsXL776MtQ1rsfP4TgDAmLIx+MgFH8GHJn4I155/LS4dcymTCRF5xtPCiNFKlUSbqmsLIi4CsATAcJiE42gCsDIyoaRLtiePSAdbD2J9w3qsbViL3x/4fWhL3IqSClx93tW4cvyVuOq8qzB7/GyMKB2R4N2IiKLzOnk4i/sanHEOEXlAVR+P85qky7inQ64lDzdVxYETB/CHA3/AqwdexcaDG/HOsXegdlb05MrJmDFuRuiYPnY6Rg6JOTGOiCjE6+TRBDMe4Z6qm1JVXZuAKhMtEvRKLiePaFo7WlH/fj3ePPgm6g/VY8uhLdjfsj/0/LnDzsXl51yOy8ZchsvGXIZLxlyCKaOmcDCeiPrwesB8YZRuq+kJAigH8BWY7isgnHgykjwGm/LicsypmoM5VeF5DU1tTdh6aCu2Ht6Kt46+hR1HdmDDvg3o7OkEAAgE51ecj6mjp+KikRfhwpEXhm7PHXYuV8ETUVRxk4e7lEiM7qdEHesPA3gWdl1GvyKkAaksrcTc6rmYWx1ektPV04XdTbux89hO7Dy+E+8cewfvHHsHr+x7BW3dbaHrhhQOQfWIatSMqMGkykmoHlGNqooqVI2owvnDz0dpYakf/0lEFACJWh6POUUNY6gFMCvO8+tcNa3uUNUXsnHf8lxTmF+IqaOnYuroqX3O92ov3mt9D7sad+Hdxnexp2kP9jbvxe6m3fjt3t+ivbu9z/Vjh47F+cPPx8ThE0O3E4ZPwITyCTiv/DyMLhvNWWBEOSpR8miAmTIbS8IRWBF5VlXvBDBbRKphKvGy2yqA8iTPfPkPn9CnpQKYxHL41GHsa96HfS37sL9lP/a37MeBEwew7fA2vPTuS+jo6ejzmqL8Ipw77FyMHzYe48vHY/yw8Rg3dBzGDRuHcUPH4dxh52Ls0LGoKKlg9xhRlklUVTfuToGp7iRoa12ty8T+5cDgGzD3U6/24tjpYzjYehAHWw/ib61/w99O/A3vnXwP7598H++dfA8HWw/iTNeZs15blF+EsUPH4pyyc3DO0HNwTtk5GFM2JnSMHjLa3JaNxqgho1CUX+TDfyHR4OHpbKt+BuB5+ZJUMHkEi6riZOdJHDp5CIdOHcKhk4dw+NRhHDl9BIdPHQ7dP3r6KI6ePoru3u6o71NeXI5RQ0aFjpGlI80xJHxbWVqJytJKjCwdiRGlIzCsaBhbN0RJ8rw8ST88LSLLVHV/mj+HsoCIoLy4HOXF5bho1EVxr+3VXrS0t+DIqSM4duYYjp0+FrptbGvE8TPHcfzMcRw+dRhvH30bjW2NONV5Kub75Us+RpSOwIiSERhROgIVJRUYURK+HV4yHBUlFagoqcDw4uEYXjK8z+3QoqFMPkQu6U4eKxFRJTeZ8iREeZIXaj1cjIuTek1Hdwea2ppCR2NbIxrPNKK5vRnNbc3mtr0ZTW1NaGlvwf6W/Whua0ZLewu6ersSxuMkvuHFw0P3hxUPQ3mRvS0ux7CiYaH7Q4uGhh4794cWDUVZURknElDWy0Ty2CsizlRdATATZktbIk8VFxSbwfhh41J6naqirbsNJ9pPoKW9BS3tLTjRcQIn2k/0uW3taEVrR2vo/vEzx9HQ3BA6f7rrdNKfWVZYhqFFQ0NHWZF57JwvKyxDWVFZ1NshhUNQVmRv7WPnXGlBKUv5U0akO3ksj7J1baLdB4kySkRCX8CpJh63nt4enO46jZMdJ9Ha0YpTnadwsvMkTnacxMnOkzjVeSp0nOw4idNdp8OPO0/iRPsJvH/yfZzqPIXTnadxuut01AkGiRTnF4f+e4YUDkFpYam5LSgNPQ7dLyhFSUFJ6FzkrfNcSUFJ6HHkuYK8AnbpDUJpTR6RicOe863WFVE65eflh7qzxmO8J+/Zq71o62rD6a7TfRKK+75znO48jbbutrMeu88dPnXYnOsy55z7ibrt4smTvFBSKc4vDt13H8UFxX2uKc4vDp1z7rtvnddEe64ovyjm/aL8IiayDOGep0QBlid5pruqqAwoS9/ndPd2o727HW1dbaGE4ty2d7eb5+zjjp6O0Pm27jZ0dHeErgmd6zHnOro70NbdhlNnToXOOec7ejpCr3UKenqhMK8QxQXhZFKUX9QnuYTORVxTmFd41jXRnivML+xzLt7jyPvRbrO15cbkQUQoyCsIjb9kmqqiu7e7TzJx7kfedvZ09jnX2dPZ5777cVdPV/ga+9j9+jNdZ9DS3tLndV29XaHHnT2d6OrpGlCrLFkFeQVxE4z7Ntq1zrmCvALzWAqinu9zTV5B6Ig8n1TMaf6ZEBHFJSLmiy2/0JfklYiqhpKKk4BSeRx5P+lb+7ru3u4+z3X3dqOrtwunO0/3Oeecd1/jvNZ53Ku9nv1cmDyIiOIQkVCXVbbr1d5QonEnpZ7enj7JZuojUxO+F5MHEdEgkSd5niVCrlQiIqKUMXkQEVHKmDyIiChlTB5ERJQyJg8iIkoZkwcREaWMyYOIiFIWyHUeIrIAQAuAGar6WKrPExFRegWu5WETA1R1HYAWEalL5XkiIkq/wCUPALMANNj7DQBmpPg8ERGlWRC7rSoiHo9M5XkRWQxgsX3YISJ/9jC2dBkF4LjfQSSBcXqLcXonG2IEsifOixJdEMTk0QKgsr/Pq+oqAKsAQETqVbXW2/C8xzi9xTi9lQ1xZkOMQHbFmeiaIHZbbUK4dVENYG2KzxMRUZoFLnmo6nMAqp2BcDswDhFZG+95IiLKnCB2WyHa9FtVnRfv+RhWeRZUejFObzFOb2VDnNkQI5BDcYqqd3sHExHR4BC4bisiIgq+nE0eIrJAROpEZJnfscRjYwz0oL+IVNif5wIRWeF3PLHYn2VdkGN0C3qcItIsIpuzIM4Zzu+n37HEYmPca3+egf6Zur47F8e7LieTRzatQs+SAf9FACrtZAUk+qXyg4jMADDP/jxniEi13zHFY38nAx0jgIWqOlNVl/sdSAJLXBNpgvozrVTVGlWdCeA+ACv9Diga+3vZYP8dNdh/V1HlZPIAV6F7SlVX2fUzgPnCC1zCU9UtqrpcRCpgfvkbEr7IJ/YLLrDxuVQE+MsYQOgPmc0iUq2qjwX1/3vEH4nVQY0TQD2ANTZpVKvqllgX5mrySLRKnfrBfpE0BfgXHwBqYRaSBlmQvzzcKgE0iUgg/0q2auzRJCIr7R8PgSUii50WfBCpagtMq2gNzM81plxNHolWqVP/LFDVJX4HEY/9C68iqP3fIlKXJV2VTouzBabrN5A/T2uvjXMzwqWJgmpe4kv8Y/8/r1PVGtfjqHI1eXAVusdEZIGzviaIY0gissI1FhPkPx6a7GDkApg++kB2qYrIYtcXR6OvwcS3yXW/AgFudQa9VWS5u6oeRZx/RzmZPLJpFbr9B1ob5L/s7M9xhTNTxO94YlgJM8BXB6DCNUYTKHZsZh3MP8ogf5mshmuySVC7WmxcFa44A/n/3aoE0OR3EAmssn841AFYFO/nyUWCRESUspxseRARUXoxeRARUcoCWRgxHtdg7bwsWLxERJSTsqrlkW2riImIclVWtTzsFLIt2bCKmIgol2VV8nCJuYrYvYd5WVnZzClTpmQyLiKirLd58+bjqjo63jVZO1XXlkxYG2/+eW1trdbXJ9yKl4iIAKC3F/jrXyFVVZsT7bWeVS0PW8Z4r124EuRVxEREwXXmDLBrF/CXv/Q9du0C2tqSeousSh4wq4irg76KmIjId6rA0aN9k8POneb2wIHwdSLABRcAF18MzJ0LXHQRsCRxCbus7bZKBrutiCjn9fQA+/aZxOAkBydRtLiGhocMAaZMMcnh4ovN/SlTgMmTgZKSPm8pIrnVbUVENGi1tZluJSdJOMeuXUBnZ/i6c84xyeGuu8IJ4uKLgfHjgTzvVmcweRARBUlrq0kK77wTThDvvGNaF05PkQhQXW2Swo03mlunNTFiREbCZPIgIvLD8eMmKThJwrn//vvha4qKTDdTbS1wzz3hJHHhhWd1NWUakwcRUbqoAkeOAG+/HU4OznH8ePi6sjJg6lSgrs4kh6lTzW1VFVAQzK9pX6MSkWkwmzW1wGxvus3PeIiI+kXVtBjcycFJGM3N4esqKkxiuO22cIKYOhWYMMF0RWWRjCcPEakCsASAAmiA2RxlBIAaEfmEPb9SVfdnOjYiorhUgUOHTGJwkoNz/8SJ8HWVlcAllwCLFpnbqVPNMXZs1iWJWDKaPERkLgBV1YcSXSci1aq6IUOhERGFOWsknMTw9tvAn/98dkti5EiTHD75SZMcLrnEHKNH50ySiCXTLY96VT2R6CJVXS8iwzMREBENck1N4eTgvnWPSTgtiTvvDLckLrkEGDMm55NELBlNHokSh4hc4HRXJZNkiIiSdvp0ODm4j0OHwteUl5ukcPvt4VbEpZeatRODNEnE4tuAuYhcAGApzBgHAAiAuQBm+RQSEeWCri6zcO6tt0xycG4bXDs4lJaawerrrzfJ4dJLTaI47zwmiST5OdtqKYC1MAPmjgqfYiGibKMK/PWvJjk4x5//bEpzdHWZa/Lzw+skPvvZcKKoqjLPUb/5mTzWqup69wkRWeZXMEQUYC0tJjHs2NE3WbS2hq+ZMAG47DLg5pvN7WWXmcRRXOxf3DnMz+TRICIPwEzXdap3LQRwv38hEZGvuruB3btNknAff/1r+JqKCpMYPvWpcJK45BJznjLGz+SxHCZpuDsYa3yKhYgy7dixs5PE228DHR3m+YICU6vpmmtMgrj8cnPLcYlA8DN5rFTVre4TIhJzV0AiylLd3cC77wLbt4ePHTv6znIaO9Ykhy9+0SSIK64wiYNdToHlSfKwazJWAJipqrNEZDrMYsB45UbUlidpUFWn43I+gMe9iImIfNDc3DdJbN/etzVRVGRmOc2bZxLE5ZebY8wYf+OmlHnV8rgPphuqGgBUdauIzEnwmg0ANgEQEVGY7qsqMHkQBZ8qsH8/sG2bObZvN7fuHerGjDEJ4otfNLdOa6Kw0LewyTteJY99qnrCJoFkLYwy22q6R/EQkVc6O01Zjm3bgK1bw8nCqeWUl2dmNV19NXD//SZJTJtmuqIoZ3mVPCpF5FEAjSJSCWAGzGB4vNpUTSLyhKreD4Qq7DbHuZ6I0q211SQGJ0ls3Wq6nZx1E0OGmOTwyU+aBDFtmlk3MWSIv3FTxnmSPFT1KdtquBPAIgCrkyxq+JjrPbbZrq79XsRERAkcO2aSw5Yt4ds9e8LPjxljksMNNwDTp5v7kyZxcR0B8G7A/FFVfRjA1oQXh1Wp6gsR5zhRm8hrzl4Tmzf3TRQHD4avqaoyCeLTnza306cD48ZxSizF5FW3VY1rgLzeNXsqHhGRJ2BKlADAPACbPYqHaHBySnZs2WKShZMwjh41z4uYQesPfzicJKZPz9i+15Q7vOq2WuTcF5ELRGQxgBZVfTrOa54XkQaYri4AWBW57oOI4lA1s5ucJFFfbxJFY6N5Pj/frLz+6EeBGTPMcfnlwNCh/sZNOcGrbqtpMAPkSwHUwUzBXZnodTZZMGEQJaJqupnq68OJor4+nCgKCszA9W23ATNnmuPyy4GSEn/jppzlVbfVBpjup1XxdgkUkTuijHPEunYOdxKkQevIEWDTJpMgnFun6yk/36zCdhJFba15zERBGeTZIkFVfT7RRar6gojcB7OHx+rIsRERKYfpxtJoXV4iUgHTsgGAWaq6fOChE/mspSWcJJzDGczOyzMrsm++2SSKWbNMi6K01N+YadDzaszjrMQhIvdGSwB2Wu9wAEtEpAaAe6RuE4DnVXVfjI9aZN9jlYjMEpHFqrrKg/8EosxobzfrJ/70J3Ns2mQ2LnJMmgR86EMmSdTWmsFsjlFQAA0oeYjIbgAzAWwBsBemxIi71EjUAXO7xew3U/28iERRjSTGVYh809trCgK++WY4WWzfbgoFAmYq7OzZZnqskyw464myxICSh6pOBgARWRilQm7aSo2ISDWAJlVtiPLcYgCLAWDixInpCoHobEePmkTx5pvAG2+YVoWzWdGwYSZBPPCASRizZpnS4kRZSlRTKUcVDCKyTFUfS3RdbW2t1tfXZyIkGmw6Okwr4o03wsc+29uan2/GJa680hyzZ5u1FXl5/sZMlCQR2ayqtfGu8WyqrlN+XUSqYLqUNNXZUiJSnmiBoYgscBKHiNSp6rr+xk2UtIMHgY0bzfHGG2Y9hVNm/LzzTJL43OeAq64y6ylY64lynFezraoBbAMAO9i9T0TuRfzCiIhStj3uNrQiUgdghYg8bE9xthV5r7PTlPB4/fVwwnBmP5WUmFlPX/yiSRRXXsnuJxqU+j2TlycAAA/rSURBVJ087IypRQCWABju+kIHgCYkGMwWkSdhFhY2uk7H3YbWtjK4VS156+hRkyBef90cmzaFWxUTJwIf/CDwgQ+YkuNXXGE2NCIa5PqdPOyMqacAPCUicyP35kjCmij7eXAbWkovVeAvfwH++MfwsXu3ea6w0LQqPvc5kzCuvho491x/4yUKKK/WeUQmgWkAKhOMeTRyG1pKu44OMz7x2mvm+OMfwyU9Ro0yLYp77zXJYuZMrtImSpJXA+blAL4CYLhzCma9R7zkwW1oyXutrabr6bXXgD/8waytaG83z114IXDrrSZRfPCD5jFLjhP1i1cD5g8DeBZAJcx4RzK4DS0N3PHjJkm8+qo5tm0zi/Py883q7PvvB665xhxjxvgdLVHO8Cp5rHMWCTrFD6PMpOpDVdfbFksdTNJZzZLslNChQ8Dvf28Sxe9/b/bWBkx301VXAV/9qkkUV1/Nsh5EaeRV8oCIPKuqdwKYbVeA1yBOt5WIzIWZmrvXnnpKRFayki718d57Jkn83/+ZW6cO1NChJkl86lNmY6OZM4HiYl9DJRpMvBwwX2/vP2Qr565I8LIKVV3qPiEi872Ih7LY4cPAK6+YZPHKK+GZUMOHm4KB990HfOQjZj/tAs/+9iGiFKXlX5+tnHsvYhRGtJqTPEe5rLHRtCg2bDDHzp3mfHk5cO21wNKlJllccYUZxyCiQBjIIsE9CHc5nfU04lTVtUaIyAMwFXkBYAbMokF2W+WyU6fMAPf69ebYvt2svSgrM8nis58FrrvODHYzWRAF1kBaHkviLQxMNHPK7mE+HWaFOhBl0SDlgO5us2J77Vpg3TpTF6qry6zS/sAHgK99DZg711SZLSz0O1oiStJAVpgn+qKP1Spxv8dWmH3PAcTeQIqyiCqwZw/wu9+ZhPHKK2bthYhpTXz5yyZZXHMNiwcSZTGvFglGTsutgGlR3BDl2sgNpEJPIXFXFwVRS4vpgvrd78yxf785f8EFwJ13AvPmma6oUaP8jJKIPOTVgPljAH4GkwAAYBbCYxl9+LWBFHmot9eU/PjNb8zxxhtAT4/Z8GjOHODBB4HrrwdqariCmyhHeZU8FkbuO27XccQUbUEgFwkGWGOjaVX87/+ahHHsmDk/cyawfDlw441mkR7HLYgGBa/WeUQmjnKYLqiYvNpAitJE1cyE+tWvgF//2rQuenuBkSOBG24AbrrJtC5Y8oNoUPJqzKMJtsihPdWMBPt5oJ8bSFEanTljxi5++UuTNN57z5yvrTVlP26+2dznFFqiQc/LbqvIIocXRF4UYwMpJ+E0InHCIa8dOgS8/DLw0ksmcbS3m7GLG24APvpR08I45xy/oySigPEqeTSIyKOuxwJgLszAeUi8DaTsYPkIj+KhWFTNKu5f/AJ48UVTshwwM6MWLwZuucXUiuJueUQUh1fJYwmAtehbjr0iwWvqbcKJvI7dVl7r7TVJ4uc/N4dTL2r2bOCf/xn4+MeBSy/lzCgiSppXyWNtlG6rZQle8xCA1UhtDxBKVne3KQPy/PMmYbz/vikkeN11ZqHexz8OjB/vd5RElKW87LZ6AEADTH0qwJRbvz/Oa1LeA4QS6O421WjXrDEJ49gxoLTUTKO94w4zhjGCPYNENHBeJY/lMEnD3e9Rk+hFqe4BQlH09JiqtKtXm1bG8eOmyODHPgbMn28GvMvK/I6SiHKMV8ljZZTV4s/Fe0GMPUAe8yie3KZq1l387GcmaRw+HE4YixaZlkZpqd9RElEO8yp5qIhMA9Cgqq323HwAj8d6gYiUu66Fqj7lUSy5a+dO4Cc/Mcf+/WbnvI9+FLjrLrMGg4UGiShDvEoeG2AXCYqIIlzkMGbyAPC0iCxT1f0exZCbjhwBnnkG+PGPTT2pvDxTaPBrXwNuu81smkRElGHpXCSYqMjhSkRM002mJLuI1AFYrqrz+hVpNujoMAv3/vu/TS2pnh6zsvs//sNUqR071u8IiWiQ83IP80iJ9vNYCWCviDhTdQWmVHvc5KGq60Rkeb8CDbrt24H//E/TymhuNlNpH3wQ+PSngSlT/I6OiCgk4/t5uCxX1ecj3iduJd6c1NoK/PSnwNNPA/X1ZmX3HXeY7VjnzmUdKSIKpIzv5+GITBz23IC3oRWRxQAWA8DEiRMH+nbps3Ur8OSTZvD79GngssuAb38buPtuU7mWiCjAfNvPI10l2VV1FYBVAFBbW6sDeS/PtbcDzz4LPPEE8OabZjrtXXcBS5aYPbxZHoSIskSeF2/Sn/08YJJF6PW21VEd5/rs9f77wD/+IzBxIvCZzwAnTphWxnvvmTGO2bOZOIgoq2R0P48YJdkdTdFeE+U9FgCoFZEFqhp3IaLvtm4FHn/cLOTr6TEVa//+781WrUwWRJTF0jZVN5p4JdmTZRNGcJOGKvDb35qksX692RvjC18wR03Cii1ERFmh391WIvKgiGyyZdX32nPDRWSOXW0eV7TEYXcSzE49PWY8Y9o0U09q505gxQrgb38DvvUtJg4iyikDaXmsg6syLhBqWWwQkSr3gLhDRPbAJBoB4KxEh+t+FRKs8wic7m4z1fZf/gV4912zHuO//gv45Ce5oRIR5ayBJI8qVX0h2hOqui9GefUl8bqqkliVHhw9PSZpPPIIsHevmWq7erVZo8G1GUSU4zyZbRXDWTsJJhrjiKzMG0iqwC9/CUyfDtxzjxnT+PnPgW3bgIULmTiIaFAYSPKYneD5qJ38dhovRGRalOOJAcSTfhs3AtdcY0qft7WZkuibN5sChXnpzMNERMEykG6rZ0XkWQD3qupJ56RNDk8BeDTyBXbPjhkwOww+DbPveUobSPni/feB5ctNzalx48zK8L/7O6Cw0O/IiIh80e/koapbReQpAFtFZATMFrQ1ABoBLI0cLLcaVdXZmvbRKLWtEi0szKyODlPJ9utfB7q6gK98BXj4YWDoUL8jIyLy1YDWeajqOgCT7Jf+DABbIlebR6gRkQvsHh7RSofMRVBmW23caFoXf/kLcOutwL/9G6fbEhFZXpVk3wcgXtJwrvumiDxp9yuvFpEl6DttN2FJ9rQ7c8aUEvnWt4AJE4Bf/9qs2yAiohCvVpgnTVWXAmZabpR9z/0tyf7aa6YU+p49wNKlZpEfd+ojIjqLb1OEYkzL3ZTxQAAz/faxx4APf9is31i/3lS+ZeIgIooq4y0PtygLCRfCzMTKnNZW09p44QVgwQLgBz8wazeIiCgm35KHiDwJoAVmdpYjsyPSb79tVoTv3WsGxL/8ZVa7JSJKgp8tjzWRK85FJHPVcv/wB+Dmm4GyMmDDBuDaazP20URE2c7PZdGNdlW5e2BhfkY++bXXzAyq8ePNCnEmDiKilPjZ8tgAu4GUiLir6j6e1k99/fVw4njlFbNinIiIUuJn8jhrA6m0V9V94w3gxhtNwmDiICLqNz+n6karsDszbR+4dy9www3AOeeYxHHuuWn7KCKiXJfxloeI7IZJEltw9sZQ6dsM6jvfMZVw160zXVZERNRvfqwwnwwAIrIwygrz9HRbnT4N/PCHZr+N889Py0cQEQ0mgVphnrbNoH7yE7MY8POfT8vbExENNrm/g5Eq8P3vA1dcAVx9td/REBHlBF/Lk2TExo3Ajh3AypVcPU5E5JHcb3l8//umwOHdd/sdCRFRzsjt5NHdDaxZA3zmM6YMCREReSLruq1EZAFMQcUZqvpY3IuPHTPbx96f2UK9RES5LqtaHjZxONvftohIXdwXHD8OzJ0LTJmSifCIiAaNrEoeAGYBaLD3G2D2TY+tsxP43OfSHRMR0aCTbd1WFRGPR0ZeICKLASy2Dztk/vw/pz2qgRsF4LjfQSSBcXqLcXonG2IEsifOixJdkG3JowVAZbwLVHUVgFUAICL1qlqbicAGgnF6i3F6KxvizIYYgeyKM9E12dZttQnh1kc1gLU+xkJENGhlVfJQ1ecAVDsD5XbgnIiIMizbuq2QcHpuX6vSFoi3GKe3GKe3siHObIgRyKE4RVUzEQgREeWQrOq2IiKiYMjZ5CEiC0SkTkSW+R1LPDbGQA/8i0iF/XkuEJEVfscTi/1Z1gU5RregxykizSKyOQvinOH8fvodSyw2xr325xnon6nru3NxvOtyMnmkvBLdR1ky6L8IQKWdsIBEv1R+EJEZAObZn+cMEan2O6Z47O9koGMEsFBVZ6rqcr8DSWCJazJNUH+mlapao6ozAdwHYKXfAUVjfy8b7L+jBvvvKqqcTB5IdSU6xaWqq+z6GcB84QUu4anqFlVdLiIVML/8DQlf5BP7BRfY+FwqAvxlDCD0h8xmEalW1ceC+v894o/E6qDGCaAewBqbNKpVdUusC3M1eSRciU6ps18kTQH+xQeAWpjFpEEW5C8Pt0oATSISyL+SrRp7NInISvvHQ2CJyGKnBR9EqtoC0ypaA/NzjSlXk0fClejULwtUdYnfQcRj/8KrCGr/t4jUZUlXpdPibIHp+g3kz9Paa+PcjHBpoqCa53cA8dj/z+tUtcb1OKpcTR5cie4xEVngrLEJ4hiSiKxwjcUE+Y+HJjsYuQCmjz6QXaoistj1xdHoazDxbXLdr0CAW51BbxVZ7q6qRxHn31FOJo9sWolu/4HWBvkvO/tzXOHMFPE7nhhWwgzw1QGocI3RBIodm1kH848yyF8mq+GabBLUrhYbV4UrzkD+f7cqATT5HUQCq+wfDnUAFsX7eXKRIBERpSwnWx5ERJReTB5ERJQyJg8iIkoZkwcREaWMyYOIiFLG5EEUwU5VXGNvV9iCdgtEZFkmi1g607ftupDmJF9TEcR1OJR7sm4zKKJMUNWFQGiNywxXUcgWe7vXWYWbDraMxSoby7pk9pS217aICGytp2wogUJZii0PorPF+9J1npuZrg8f6Epkuwgx0GVkKPsxeRBFiFeRwLYCZsDUUQp1Kdn9GmbY4nzVzr4i7qq0TskP2/0VL0Eshlnh3Yd9/zp3oUL7fnvt+TVZUgKDcgCTB1GKbO2fBnt/HUwZa+d8C0wByXUwNdWWAKH9Rmps99cqAPE2A6qxhf6ifXafPWrs+zXAlL1Y7nrd3qDWzaLcwDEPIm843VmNrvvuOkZ3Amh0faHXxnmvWMXo3J/R5/oo+y40wRQFjbkfA9FAMHkQeSCipRCt1VABYK3rSz7lMZNYrRHELrYX2AqzlP3YbUWUGWvg2tEyA11K2VDBlbIYkwdRDHadxRKY8v6LncFoZ4tOe865v8wOjs8DsNBeuwRAnZ02GxqrSGIdxib3wHe8z7DvVRtlX/maeFuIEg0US7ITBYxNEHUD2ZtCRFao6nIPwyLqgy0PooAZ6OI+2xoJ8r7jlAOYPIgCSFVX9afMiNPdxdXllG7stiIiopSx5UFERClj8iAiopQxeRARUcqYPIiIKGVMHkRElDImDyIiStn/BztUANw1uDxTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rc('text', usetex=True)  #allow the use of Latex for math expressions and equations\n",
    "plt.rc('font', family='serif') #specify the default font family to be \"serif\"\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "# infiltration rate\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(trange, f, color='green')\n",
    "plt.xlabel('Time (hr)',fontsize=12)\n",
    "plt.ylabel('Infiltration \\n rate (cm/hr)',fontsize=12)\n",
    "plt.xlim([0,8])\n",
    "plt.ylim([0,1.5])\n",
    "\n",
    "# cumulative infiltration\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(trange, F, color='red')\n",
    "plt.xlabel('Time (hr)',fontsize=12)\n",
    "plt.ylabel('Cumulative \\n infiltration (cm)',fontsize=12)\n",
    "plt.xlim([0,8])\n",
    "plt.ylim([0,3])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "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": 2
}