gw_modelling_python_session_02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0f4ee7a8",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Groundwater Modelling in Python - Session 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c262a0a4-6132-443a-a472-243be6908bbb",
   "metadata": {},
   "source": [
    "<img src=\"figs/part_of_cover_bakker_post.png\" width=\"800px\"> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e62b54a2-dcb3-4512-8fc9-85d871e7e2a9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Before starting we need to import NumPy and Matplotlib and set some defaults\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# and set some parameters to make the figures look good\n",
    "plt.rcParams[\"figure.autolayout\"] = True # same as tight_layout after every plot\n",
    "plt.rcParams[\"figure.figsize\"] = (9, 3.5) # set default figure size\n",
    "plt.rcParams[\"contour.negative_linestyle\"] = 'solid' # set default line style\n",
    "plt.rc('font', size=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1361e7de-488b-40bb-aeef-85739f75ff4f",
   "metadata": {},
   "source": [
    "Sections and Exercises with a star are additional material that is not covered in detail during the workshop."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11572768-6b42-4a3b-9576-fdaa695eb239",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "This session deals with one-dimensional transient flow, which means that the head and flow are a function of $x$ and time $t$, respectively. As for steady flow, the starting point of the derivation of the differential equation is the volume balance for a small piece of the aquifer\n",
    "\\begin{equation}\\label{continuity}\n",
    "\\text{Volume in} - \\text{Volume out} = \\text{Increase in volume}\n",
    "\\end{equation}\n",
    "Recall from the previous session that the inflow consists of horizontal flow from the left, $Q_x(x, t)$ [L$^2$/T], and recharge at the top, $N$ [L/T] and that outflow consists of horizontal flow at the right, $Q_x(x+\\Delta x, t)$ (see Figure). \n",
    "\n",
    "<IMG src=\"figs/merged_nbs5_4_0.png\"  width=400>>\n",
    "\n",
    "Unlike for steady flow, the change in storage within the volume is no longer zero. The increase in storage equals the increase in head multiplied by the storage coefficient $S$, so\n",
    "\\begin{equation}\n",
    "Q_x(x,t)\\Delta t + N\\Delta x \\Delta t - Q_x(x+\\Delta x,t)\\Delta t =\n",
    "S\\Delta x[h(x, t + \\Delta t) - h(x, t)]\n",
    "\\end{equation}\n",
    "Division by $\\Delta x \\Delta t$ and rearrangement of terms gives\n",
    "\\begin{equation}\n",
    "\\frac{Q_x(x + \\Delta x, t) - Q_x(x, t)}{\\Delta x}\n",
    "= -S \\frac{h(x, t + \\Delta t) - h(x, t)}{\\Delta t} + N\n",
    "\\end{equation}\n",
    "In the limit for $\\Delta x \\to 0$ and $\\Delta t \\to 0$ the volume balance turns into the differential equation for one-dimensional horizontal transient flow\n",
    "\\begin{equation}\\label{basicdeq1dtr}\n",
    "\\frac{\\partial Q_x}{\\partial x} = - S\\frac{\\partial h}{\\partial t} + N\n",
    "\\end{equation}\n",
    "\n",
    "Substitution of $Q_x = -T\\tfrac{\\partial h}{\\partial x}$ gives\n",
    "\\begin{equation} \\label{deqtransient}\n",
    "\\frac{\\partial^2 h}{\\partial x^2} = \\frac{S}{T} \\frac{\\partial h}{\\partial t} - \\frac{N}{T}\n",
    "\\end{equation}\n",
    "The term $T/S$ is also referred to as the aquifer diffusivity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01f72e17-bb62-45ba-9593-d9938a671224",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Solution 1. Riverbank storage"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f615022c-93fc-4722-aa82-d88ac6ae5754",
   "metadata": {},
   "source": [
    "## Problem definition\n",
    "The figure below shows an aquifer bounded on the left by a fully penetrating canal or river that is in direct hydraulic contact with the aquifer. Towards the right, the aquifer has an infinite extent in the $x$ direction. Initially, the head in the aquifer is equal to the canal level $h_0$ everywhere. At time $t=t_0$, the stage in the canal is raised by an amount $\\Delta h$. It is further assumed that the head change does not extend indefintively inland, so $h|_{x\\to\\infty, t\\ge 0} = h_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e92ce962-fc0f-4d5a-861b-f7f959a1be4d",
   "metadata": {},
   "source": [
    "<img src=\"figs/fig5.1.png\" width=400>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d27d0c7-9343-4c5b-848d-85fb870f3d0f",
   "metadata": {
    "tags": []
   },
   "source": [
    "The solution for the head under these conditions is\n",
    "\\begin{equation}\\label{hedelman}\n",
    "h(x, t) = \\Delta h \\, \\text{erfc}(u) + h_0 {\\hskip 2em} t>t_0\n",
    "\\end{equation}\n",
    "where \n",
    "\\begin{equation}\\label{uedelman}\n",
    "u = \\sqrt{\\frac{Sx^2}{4T(t-t_0)}}\n",
    "\\end{equation}\n",
    "and erfc is the complimentary error function\n",
    "\\begin{equation}\n",
    "\\text{erfc}(u) = \\int_u^\\infty \\frac{2}{\\sqrt{\\pi}} \\text{e}^{-\\tau^2}\\text{d}\\tau\n",
    "\\end{equation}\n",
    "The solution for the discharge vector is\n",
    "\\begin{equation}\n",
    "Q_x= T\\Delta h \\frac{2u}{x\\sqrt{\\pi}} \\text{e}^{-u^2} {\\hskip 2em} t>t_0\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29c7f40d-2c35-41c1-8b15-43282049f19e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# parameters\n",
    "T = 100 # transmissivity, m^2/d\n",
    "S = 0.2 # storage coefficient, -\n",
    "delh = 2 # change in river level, m\n",
    "t0 = 0 # time of change in river level, d\n",
    "h0 = 0 # starting head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1039f70f-b172-4e1b-9c3e-964e26664efc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# solution\n",
    "from scipy.special import erfc\n",
    "\n",
    "# Function for the head\n",
    "def head(x, t, T, S, delh=1, t0=0):\n",
    "    u = np.sqrt(S * x ** 2 / (4 * T * (t - t0)))\n",
    "    return delh * erfc(u)\n",
    "\n",
    "# Function for the discharge vector\n",
    "def disvec(x, t, T, S, delh, t0=0):\n",
    "    u = np.sqrt(S * x ** 2 / (4 * T * (t - t0)))\n",
    "    return T * delh * 2 * u / (x * np.sqrt(np.pi)) * np.exp(-u ** 2) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3715969b-a1fa-4df3-8de0-55fde7419f69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "x = np.linspace(1e-12, 200, 100)\n",
    "plt.subplot(121)\n",
    "for t in [1, 10, 100]:\n",
    "    h = head(x, t, T, S, delh, t0) + h0\n",
    "    plt.plot(x, h, label=f'time={t} d')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$ (m)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.legend()\n",
    "plt.subplot(122)\n",
    "for t in [1, 10, 100]:\n",
    "    Qx = disvec(x, t, T, S, delh, t0)\n",
    "    plt.plot(x, Qx, label=f'time={t} d')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$ (m)')\n",
    "plt.ylabel('$Q_x$ (m$^2$/d)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "535346f4-718f-43f9-8b7d-672b24a07e66",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "tags": [
     "hide_output"
    ]
   },
   "outputs": [],
   "source": [
    "# basic plot head and Qx vs t\n",
    "t = np.linspace(1e-12, 100, 100)\n",
    "plt.subplot(121)\n",
    "for x in [50, 100, 200]:\n",
    "    h = head(x, t, T, S, delh, t0)\n",
    "    plt.plot(t, h, label=f'distance {x} m')\n",
    "plt.grid()\n",
    "plt.xlabel('$t$ (d)')\n",
    "plt.ylabel('head (m)')    \n",
    "plt.legend()\n",
    "plt.subplot(122)\n",
    "for x in [50, 100, 200]:\n",
    "    Qx = disvec(x, t, T, S, delh, t0)\n",
    "    plt.plot(t, Qx, label=f'distance {x} m')\n",
    "plt.grid()\n",
    "plt.xlabel('$t$ (d)')\n",
    "plt.ylabel('$Q_x$ (m$^2$/d)');    \n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d619d0c-eb0e-4f14-b585-24c7f54f43cb",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "In this flow problem, the head is only a function of the dimensionless parameter $u$, also referred to as the similarity variable (which combines the independent variables). Different combinations of the independent variables $x$, $t$, $T$ and $S$ can give the same head. For example, the head change at $x=100$ m and $t=10$ d is equal to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7294540-9ed1-4d12-aa5c-a2f9c9594faf",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(f'head at x = 100 m, t = 10 d for S = 0.2:')\n",
    "print(f'h = {head(100, 10, T, 0.2, delh, t0):.6f} m')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd56e024-bb45-4944-b62c-057255e30727",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "If the storage coefficient is not $S=0.2$ but $S=0.002$ and $T = 1$ m$^2$/d instead of $T = 100$ m$^2$/d, the same head is obtained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2abc933b-e107-4db0-962d-8b8a2b5ff0ff",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(f'head at x = 100 m, t = 10 d for T = 1 m2/d and S = 0.002:')\n",
    "print(f'h = {head(100, 10, 1, 0.002, delh, t0):.6f} m')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d991de43-3210-43e6-af0e-4b506ba8954c",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "### Exercise 1.1\n",
    "At what distance from the canal is the head change at $t=10$ d and $S=0.002$ equal to the head change at $x=10$ m and $t=10$ d when $S=0.2$? Use $T = 100$ m$^2$/d."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb866b2a-6923-48cc-9be6-75513ce9a2bc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "bbf9de5b-b755-4a69-a6a8-5097f13652af",
   "metadata": {},
   "source": [
    "### Water balance\n",
    "\n",
    "The total inflow $Q_{in}$ from the canal into the aquifer over a period $\\Delta t$  may be computed through integration of the discharge vector\n",
    "\\begin{equation}\n",
    "Q_{in} = \\int_{t_0}^{t_0+\\Delta t} Q_x|_{x=0,t} = 2 \\Delta h \\sqrt{\\frac{ST \\Delta t}{\\pi}}\n",
    "\\end{equation}\n",
    "The total inflow must be equal to the total increase in storage. \n",
    "The total increase in storage may be computed through numerical integration of the transient head in the aquifer at $t=\\Delta t$ and multiplication by the storage coefficient. The integation may be carried out, for example, by using the `quad` method of `scipy.integrate`. The `quad` function returns a number of items of which the first one is the computed integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86c2fb4a-2ded-49f7-bffc-7fd5a372e06e",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "t = 10 # time period for which water balance is checked\n",
    "Qin = 2 * delh * np.sqrt(S * T * t / (np.pi)) # from equation\n",
    "from scipy.integrate import quad\n",
    "stored = S * quad(head, 1e-12, np.infty, args=(t, T, S, delh, 0))[0]\n",
    "print(f'total inflow from canal  : {Qin:.6f} m^3')\n",
    "print(f'total increase in storage: {stored:.6f} m^3')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b32c8361-2d95-42fa-b66b-379e63469d53",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "### Exercise 1.2\n",
    "Plot the increase in groundwater storage as a function of time between $t = 1$ and $t = 100$ d. Why does the volume keep increasing?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90971888-7611-4260-a086-3fca1e1d5221",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a0fd26a7-3386-42e2-a045-3bf8a52458a2",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Effect of a flood wave$^*$\n",
    "The governing differential equation is linear, which means that superposition of the solutions can be applied to simulate the effect of multiple river level stages, for example the passage of a flood wave. For a single flood wave the boundary condition at the river is\n",
    "\\begin{equation}\n",
    "h|_{x=0}=h_0 {\\hskip 2em} t \\le t_0\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "h|_{x=0}=h_0 + \\Delta h {\\hskip 2em} t_0\\le t \\le t_1\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "h|_{x=0}=h_0  {\\hskip 2em} t \\ge t_1\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6aebfc2-ace4-4e4b-8522-d7fa662bdb88",
   "metadata": {},
   "outputs": [],
   "source": [
    "h0 = 0 # m\n",
    "delh = 2 # m\n",
    "t0 = 0\n",
    "t1 = 4 # days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "422a4228-2ffb-4ba1-b49f-4cc26e732387",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5, 3))\n",
    "plt.plot([-1, t0, t0, t1, t1, 3 * t1], [h0, h0, h0 + delh, h0 + delh, h0, h0])\n",
    "plt.xlabel('time (days)')\n",
    "plt.ylabel('river stage (m)')\n",
    "plt.xticks(np.arange(0, 12.1, 4))\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ad615f2-fcf7-407f-ae92-cd384d814b35",
   "metadata": {},
   "source": [
    "The solution for the head is obtained from superposition as\n",
    "\\begin{equation}\n",
    "h(x, t) = h_0 + \\Delta h \\text{erfc}(u_0) {\\hskip 2em} t_0\\le t \\le t_1\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "h(x, t) = h_0 + \\Delta h \\text{erfc}(u_0) - \\Delta h \\text{erfc}(u_1) {\\hskip 2em} t \\ge t_1\n",
    "\\end{equation}\n",
    "where\n",
    "\\begin{equation}\n",
    "u_0 = \\sqrt{\\frac{Sx^2}{4T(t-t_0)}} {\\hskip 2em} u_1 = \\sqrt{\\frac{Sx^2}{4T(t-t_1)}}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eee19835-bce5-46c6-8511-c983f2594542",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# basic plot head and Qx vs x\n",
    "x = np.linspace(1e-12, 200, 100)\n",
    "plt.subplot(121)\n",
    "for t in [t1 + 1e-12, 2 * t1, 3 * t1]:\n",
    "    h = head(x, t, T, S, delh, t0) - head(x, t, T, S, delh, t1) + h0\n",
    "    plt.plot(x, h, label=f'time={t:.0f} d')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$ (m)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.legend()\n",
    "plt.subplot(122)\n",
    "t = 2 * t1\n",
    "Qx = disvec(x, t, T, S, delh, t0) - disvec(x, t, T, S, delh, t1)\n",
    "plt.plot(x, Qx, 'C1', label=f'time={t} d')\n",
    "t = 3 * t1\n",
    "Qx = disvec(x, t, T, S, delh, t0) - disvec(x, t, T, S, delh, t1)\n",
    "plt.plot(x, Qx, 'C2', label=f'time={t} d')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$ (m)')\n",
    "plt.ylabel('$Q_x$ (m$^2$/d)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8507f094-568e-487a-b792-43c08136dc78",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "tags": []
   },
   "source": [
    "### Exercise 1.3$^*$\n",
    "Model the effect of a second flood wave that has $\\Delta h = 1$ m and passes between 8 and 10 days. Plot the $h$ and $Q_x$ as a function of $x$ for the same times as in the above figure."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14b3c74c-0a0d-4954-9962-e2f569c32d28",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Solution 2. Transient recharge"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43ddadd7-0086-44f9-a095-e4e08c803ff5",
   "metadata": {},
   "source": [
    "## Problem definition\n",
    "This problem is the transient version of Solution 1 of Session 1. An aquifer is bounded on the left and right sides by two long parallel rivers that fully penetrate the aquifer; for simplicity, the heads in the two rivers are equal to 0. The rivers are a distance $L$ apart and in direct hydraulic contact (i.e., no entry resistance). Initially, the head in the aquifer is equal to the head in the rivers everywhere\n",
    "\\begin{equation}\n",
    "h\\vert_{x=-L/2, t} = h\\vert_{x=L/2, t} = h\\vert_{x, t=0}\n",
    "\\end{equation}\n",
    "At time $t=0$, it starts to rain and the groundwater recharge is equal to $N$. \n",
    "The transmissivity of the aquifer is approximated as constant and equal to $T$. The storage coefficient is equal to $S$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ff1708-7cd8-4946-86e4-cd4c611a08c8",
   "metadata": {},
   "source": [
    "<img src=\"figs/merged_nbs5_67_0.png\" width=400>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cae5445c-9b4c-48dd-8464-612062441719",
   "metadata": {
    "tags": []
   },
   "source": [
    "The solution is somewhat complicated, but the Python script is only a few lines. \n",
    "\\begin{equation}\\label{KvdL}\n",
    "h = \\frac{-N}{2T}\\left(x^2 - \\tfrac{1}{4}L^2\\right) -\n",
    "\\frac{4NL^2}{\\pi^3T} \\sum_{n=0}^\\infty\n",
    "\\frac{(-1)^n}{(2n + 1)^3} \\cos\\left[\\frac{(2n+1)\\pi x}{L}\\right]\n",
    "\\exp\\left[-\\frac{(2n+1)^2\\pi^2 Tt}{SL^2}\\right]\n",
    "\\end{equation}\n",
    "Note that the first term of the solution is the solution for steady flow (but the origin of the coordinate system is different than for Solution 1 of Session 1):\n",
    "\\begin{equation}\n",
    "h_\\text{steady} = \\frac{-N}{2T}\\left(x^2 - \\tfrac{1}{4}L^2\\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ecfec5f-1a40-4edb-a855-dd223676442c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# parameters\n",
    "L = 1000 # aquifer length, m\n",
    "S = 0.1 # storage coefficient, -\n",
    "T = 200 # transmissivity, m^2/d\n",
    "N = 0.001 # recharge rate, m/d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "185af0f4-8947-4113-bfe4-8e150e46e6ca",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# solution\n",
    "def hsteady(x, T=T, L=L, N=N):\n",
    "    return -N / (2 * T) * (x ** 2 - L ** 2 / 4)\n",
    "\n",
    "def head(x, t, T=T, S=S, L=L, N=N, nterms=10):\n",
    "    h = 0\n",
    "    for n in range(nterms):\n",
    "        h += (-1)**n / (2 * n + 1)**3 * \\\n",
    "        np.cos(((2 * n + 1) * np.pi * x) / L) * \\\n",
    "        np.exp(-((2 * n + 1)**2 * np.pi**2 * T * t) / (S * L ** 2))\n",
    "    h = hsteady(x, T, L, N) - 4 * N * L**2 / (np.pi ** 3 * T) * h\n",
    "    return h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bdfc298-8a27-4a7d-aeeb-9bb76a63213c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# plot\n",
    "x = np.linspace(-L / 2, L / 2, 100)\n",
    "for t in [50, 100, 150]:\n",
    "    plt.plot(x, head(x, t, nterms=10), label=f't={t} d')\n",
    "plt.plot(x, head(x, 0, nterms=0), 'k', label='steady'); # 0 terms gives steady solution\n",
    "plt.grid()\n",
    "plt.xlabel('x (m)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40ab06fe",
   "metadata": {},
   "source": [
    "Next, the head in an observation well at $x=200$ is plotted as a function of time for $t$ varying from 0 to 300 days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da9f20ce-a98c-4d36-b17c-52cc8fbc0d5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0, 300, 100)\n",
    "h = head(200, t)\n",
    "plt.plot(t, h)\n",
    "plt.xlabel('time (d)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49e0b459-6b15-409b-931d-48d051a11813",
   "metadata": {},
   "source": [
    "The head in an observation well may be approximated as\n",
    "\\begin{equation}\\label{KvdLapprox}\n",
    "h \\approx \\frac{-N}{2T}\\left(x^2 - \\tfrac{1}{4}L^2\\right)(1 - \\text{e}^{-t/\\tau})\n",
    "\\end{equation}\n",
    "where $\\tau$ is the characteristic time of the system, defined as\n",
    "\\begin{equation}\\label{aKvdL}\n",
    "\\tau = \\frac{SL^2}{\\pi^2T}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f1759b-7c5f-4406-ac75-d0266608cb0d",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Response time\n",
    "95% of the final head is reached at $t=3\\tau= 3SL^2/(\\pi^2 T)$, which is also called the response time (or the memory) of the system. For this case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36f95270-d8ae-44de-b435-26f1a28816f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = S * L ** 2 / (np.pi ** 2 * T)\n",
    "t_memory = 3 * tau\n",
    "print(f'Memory of the system: {t_memory:.0f} d')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77404a67-baa3-4236-baa0-91c6daf93e42",
   "metadata": {},
   "source": [
    "This means that if it rains for 1 day, it takes approximaty 152 days before most of the water has left the aquifer and has flown into the river for this example."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc8035b9",
   "metadata": {},
   "source": [
    "### Exercise 2.1\n",
    "Use both the exact solution and the approximate solution to plot the head in an observation well at $x=200$ vs. time for $t$ varying from 0 to 300 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0845658-8652-423e-b9bc-e84c37c128f9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "92c9ca5b-7d2f-441e-b22d-37dc5db81edc",
   "metadata": {},
   "source": [
    "## Block response\n",
    "The head response due to a recharge event of intensity $N$ from $t_0$ to till $t_1$ is called the block response and is obtained by superposition as\n",
    "\\begin{equation}\\label{KvdLapprox2}\n",
    "h = h_\\text{steady}(1 - \\text{e}^{-(t-t_0)/\\tau}) - \n",
    "h_\\text{steady}(1 - \\text{e}^{-(t-t_1)/\\tau}) \\qquad t\\ge t_1\n",
    "\\end{equation}\n",
    "which may be simplified to\n",
    "\\begin{equation}\\label{KvdLapprox3}\n",
    "h = h_\\text{steady}(\\text{e}^{-(t - t_1)/\\tau} - \\text{e}^{-(t - t_0)/\\tau})\n",
    "\\qquad t\\ge t_1\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81535718-faa8-4235-b535-764263ca4e52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "t0 = 0 # start of recharge event, d\n",
    "t1 = 1 # end of recharge event, d\n",
    "N = 0.01 # recharge, m/d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc13d1fa-c827-4455-9a9d-15fa39dafa8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(t1, 300, 100)\n",
    "h = hsteady(200, T, L, N) * (np.exp(-(t - t1) / tau) - np.exp(-(t - t0) / tau))\n",
    "plt.plot(t, h)\n",
    "plt.xlabel('time (d)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "677e7dce-ce6a-426b-99ef-4bfc3ff1a98d",
   "metadata": {},
   "source": [
    "# Solution 3. Time series analysis of measured head series"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35d89155-d420-4339-9e59-051e83dac567",
   "metadata": {},
   "source": [
    "The head response due to a series of recharge events is computed through superposition. \n",
    "The unit block response, the block response to a recharge event with unit intensity, is written as a function of two parameters: $A$ and $a$ (compare to the block response defined above). \n",
    "\\begin{equation}\\label{KvdLapprox3}\n",
    "h = A(\\text{e}^{-(t - t_1)/a} - \\text{e}^{-(t - t_0)/a})\n",
    "\\qquad t\\ge t_1\n",
    "\\end{equation}\n",
    "A Python function is written for the unit block response response. The function takes as input arguments the parameters `A` and `a`, an array of times `t`, and the start of the recharge event `tstart` and the length of the recharge event `delt`. The function returns an array of head responses for all times in the array `t`, where the head is 0 for times before `tstart`.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9f8c5a4-a337-4819-b80c-962a07495dc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def unitblockresponse(A, a, t, tstart, delt):\n",
    "    h = np.zeros(len(t))\n",
    "    h[t >= tstart + delt] = A * (np.exp(-(t[t >= tstart + delt] - (tstart + delt)) / a) - \n",
    "                                 np.exp(-(t[t >= tstart + delt] - tstart) / a))\n",
    "    return h"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbad22f0-f1a9-4d62-8517-62bfd51254f4",
   "metadata": {},
   "source": [
    "The head response is plotted for a recharge of $N=0.001$ m/d starting at $t=0$ and ending at $t=30$ days. The head is computed every 30 days up to 180 days. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efd272fa-33b4-4408-bd03-40276181834f",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 200 # d\n",
    "a = 40 # d\n",
    "tstart = 0\n",
    "delt = 30\n",
    "t = np.arange(0, 200, 30)\n",
    "N = 0.001 # m/d\n",
    "h = N * unitblockresponse(A, a, t, tstart, delt)\n",
    "plt.plot(t, h, marker='.')\n",
    "plt.xticks(np.arange(0, 200, delt))\n",
    "plt.grid()\n",
    "plt.title('head response to 1 mm/d for 30 days');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "345303cf-18a3-4fba-a6a0-a4c2b0f224a2",
   "metadata": {},
   "source": [
    "As an example, the head response for 5 months of recharge is computed below. First, the head response is plotted for each month separately. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "211db485-e31d-4fd6-9401-d403529b7108",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(0, 301, 30)\n",
    "N = 0.001 * np.array([0.4, 1.4, 0.8, 0.2, 1.0]) # 5 months of average daily recharge in mm/d\n",
    "plt.figure(figsize=(9, 7))\n",
    "for i in range(5):\n",
    "    plt.subplot(5, 1, i + 1)\n",
    "    plt.plot(t, N[i] * unitblockresponse(A, a, t, i * delt, delt))\n",
    "    plt.xticks(np.arange(0, 301, delt), 11 * [''])\n",
    "    plt.ylim(0, 0.15)\n",
    "    if i == 2: plt.ylabel('head (m)')\n",
    "    plt.grid()\n",
    "plt.xticks(np.arange(0, 301, delt), np.arange(0, 301, delt))\n",
    "plt.xlabel('days');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "609f825a-eee9-41b0-b4e3-bba23ae383a1",
   "metadata": {},
   "source": [
    "The head response is obtained through superposition by simply adding the five reponses together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5bc12e9-dc3c-422f-bf72-4bbdeee0c89e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "h = np.zeros(len(t))\n",
    "for i in range(5):\n",
    "    h += N[i] * unitblockresponse(A, a, t, i * delt, delt)\n",
    "plt.plot(t, h)\n",
    "plt.xticks(np.arange(0, 301, delt), np.arange(0, 301, delt))\n",
    "plt.xlabel('days')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78becf0b-afc4-49ed-9ef7-f31eb2ef40bd",
   "metadata": {},
   "source": [
    "## Four years of monthly recharge\n",
    "Consider four years of monthly recharge (a month is defined as 365 days divided by 12 here). The recharge is plotted and the head in the aquifer is computed and plotted at the end of every month for the chosen values of `A` and `a`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84c46d90-9846-44dc-8e75-d044adf59d32",
   "metadata": {},
   "outputs": [],
   "source": [
    "recharge = 0.001 * np.array( \n",
    "[ 1.2,  1.5,  0.8, -1.2,  0.2, -2.9, -0.6,  2.8,  1.2,  0.9,  3.8,\n",
    "  1.8,  1.9,  1.8, -0.6, -1.3, -2.5,  0.6,  0.9,  1.8, -0.3,  0.8,\n",
    " -0.2,  3.8,  2.8,  0.3, -0.5,  0.3, -1.1,  1.5,  1. , -1.5,  0.1,\n",
    "  1.5,  1. ,  4. ,  1.3,  1.3,  0.3, -1.1,  0.6, -0.8, -1.2, -1.6,\n",
    "  1.3,  1.4,  2.4,  2.6]) # monthly recharge, m/d\n",
    "#\n",
    "tmax = 4 * 365 # maximum time, d\n",
    "delt = 365 / 12 # time step, d\n",
    "t = np.arange(0, tmax + 1, delt) # time, d\n",
    "#\n",
    "A = 200 # d\n",
    "a = 20 # d\n",
    "h = np.zeros(48 + 1)\n",
    "for i in range(48):\n",
    "    h += recharge[i] * unitblockresponse(A, a, t, tstart=i * delt, delt=delt)\n",
    "#\n",
    "plt.subplot(211)\n",
    "plt.step(t, 1000 * np.hstack((recharge, recharge[-1])), where='post')\n",
    "plt.xticks(np.linspace(0, 4 * 365, 5), 5 * [''])\n",
    "plt.ylabel('recharge (mm/d)')\n",
    "plt.grid()\n",
    "plt.subplot(212)\n",
    "plt.plot(t, h)\n",
    "plt.xticks(np.linspace(0, 4 * 365, 5), np.linspace(0, 4, 5))\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d760216-7aa3-4b9d-be53-15ddb57a7fea",
   "metadata": {},
   "source": [
    "Head observations are available at the end of every month."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebb729bf-6230-41e8-8820-413223d57941",
   "metadata": {},
   "outputs": [],
   "source": [
    "hobs = np.array(\n",
    "    [ 0.  ,  0.33,  0.6 ,  0.56, -0.01,  0.04, -0.76, -0.61,  0.42,\n",
    "      0.58,  0.58,  1.36,  1.27,  1.23,  1.19,  0.51, -0.08, -0.74,\n",
    "     -0.27,  0.1 ,  0.55,  0.24,  0.35,  0.14,  1.11,  1.41,  0.88,\n",
    "      0.35,  0.27, -0.15,  0.31,  0.46, -0.14, -0.07,  0.37,  0.49,\n",
    "      1.37,  1.14,  0.99,  0.64,  0.06,  0.18, -0.11, -0.39, -0.66,\n",
    "      -0.03,  0.38,  0.87,  1.21])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c43591be-bcf7-4534-980b-7a18246e13bd",
   "metadata": {},
   "source": [
    "The model is plotted versus the observations. For these values of `A` and `a`, the fit isn't very good. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06538ae4-07a5-46e2-bbbb-947b4534b95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 200 # d\n",
    "a = 20 # d\n",
    "h = np.zeros(48 + 1)\n",
    "for i in range(48):\n",
    "    h += recharge[i] * unitblockresponse(A, a, t, tstart=i * delt, delt=delt)\n",
    "#\n",
    "plt.plot(t, hobs, 'k.', label='observed')\n",
    "plt.plot(t, h, label='model')\n",
    "plt.xticks(np.linspace(0, 4 * 365, 5), np.linspace(0, 4, 5))\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid()\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "380acbf1-e6b2-4bec-88db-ce56a1fce70b",
   "metadata": {},
   "source": [
    "### Exercise 3.1\n",
    "Find the optimal parameters of `A` and `a` by varying them manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "684858fa-ab9f-4c6a-8dba-30382765772e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7e1b4ccd-7e0e-4cb3-b961-43f7e5480161",
   "metadata": {},
   "source": [
    "## Minimize sum of squared errors\n",
    "The optimal values of `A` and `a` are obtained by fitting the modeled heads $h_\\text{model}$ to the observed heads $h_\\text{obs}$. The objective function $f_\\text{obj}$ is the sum of squared errors\n",
    "\\begin{equation}\n",
    "f_\\text{obj} = \\sum(h_\\text{obs} - h_\\text{model})^2\n",
    "\\end{equation}\n",
    "The optimal values of `A` and `a` are obtained by minimizing the  objective function $f_\\text{obj}$  using the `fmin` function of `scipy.optimize`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ece5d739-0746-4fb4-a0f2-f71c82af2d41",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fobj(p, return_heads=False):\n",
    "    A, a = p\n",
    "    h = np.zeros(48 + 1)\n",
    "    for i in range(48):\n",
    "        h += recharge[i] * unitblockresponse(A, a, t, tstart=i * delt, delt=delt)\n",
    "    if return_heads: return h\n",
    "    return np.sum((h - hobs) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d289a5b0-ffdd-43e4-bc07-e3f269c1d731",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fmin\n",
    "A, a = fmin(fobj, [200, 20])\n",
    "print(f'Estimated parameters A: {A:.2f} d, a: {a:.2f} d')\n",
    "hm = fobj([A, a]) # compute modeled heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d7d3f2a-991a-4871-8312-cf99c156b379",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = fobj([A, a], return_heads=True) # compute heads with optimal parameters\n",
    "plt.plot(t, hobs, 'k.', label='observed')\n",
    "plt.plot(t, h, label='model')\n",
    "plt.xticks(np.linspace(0, 4 * 365, 5), np.linspace(0, 4, 5))\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid()\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "532bff5f-5b62-473e-870f-bbcafbd7fc10",
   "metadata": {},
   "source": [
    "The above procedure to fit a time series of head observations is referred to as time series analysis. More information on time series analysis can be found in [Collenteur et al., 2019. Pastas: open source software for the analysis of groundwater time series. Groundwater, 57(6), pp.877-885](https://ngwa.onlinelibrary.wiley.com/doi/full/10.1111/gwat.12925). Or check out the [Pastas](https://github.com/pastas/pastas) website. An example analysis of a 10 year time series is shown below, including the model fit and the contribution of rainfall and evaporation to the head variation (from the paper by Collenteur et al., referenced above). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e0b724c-a59a-4674-86cf-07172e5a9763",
   "metadata": {},
   "source": [
    "<img src=\"figs/fig_pastas.png\" width=800>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5da01629-91c9-4a8a-8cce-6addcd5690b7",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Solution 4. Pumping test analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f51439-8aa1-452f-9344-07bd31f132fe",
   "metadata": {},
   "source": [
    "### Confined aquifer (Theis)\n",
    "\n",
    "The Theis solution is a famous formula in groundwater hydrology to compute the drawdown around a well in a confined aquifer. The aquifer is approximated as fully-confined and the transmissivity $T$ and storage coefficient $S$ are approximated as constant. The well extraction rate is $Q$.\n",
    "\n",
    "<img src=\"figs/figwell.png\"  width=500>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db23802-f17a-42de-b610-6fc6d0ef5bd0",
   "metadata": {
    "tags": []
   },
   "source": [
    "The solution for a well that starts pumping at $t=t_0$ (Theis solution) is\n",
    "\\begin{equation}\\label{headtheis}\n",
    "h = - \\frac{Q}{4\\pi T}\\text{E}_1\\left( u\\right) {\\hskip 2em} t\\ge t_0\n",
    "\\end{equation}\n",
    "where\n",
    "\\begin{equation}\n",
    "u =  \\frac{S r^2}{4T(t-t_0)}\n",
    "\\end{equation}\n",
    "and E$_1$ is the exponential integral\n",
    "\\begin{equation}\n",
    "\\text{E}_1(u) = \\int_u^{\\infty} \\frac{\\exp(-s)}{s}\\text{d} s \n",
    "\\end{equation}\n",
    "\n",
    "Because the pumped water can only come from storage in the fully confined aquifer, the Theis solution does not reach steady state. It is only suitable to analyze the early stages of a pumping test, before other sources of water (i.e. from an overlying aquifer) start to contribute significantly to the well flow. \n",
    "\n",
    "The following example shows how the Theis solution is implemented in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c957ba30-7aa5-4029-a255-cf5d1007062c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# parameters\n",
    "T = 200 # transmissivity of aquifer, m^2/d\n",
    "S = 0.0005 # storage coefficient of aquifer, -\n",
    "Q = 800 # discharge of well, m^3/d\n",
    "rw = 0.3 # radius of well, m\n",
    "t0 = 0 # start of pumping, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5355c2b9-afbc-43c7-84b2-8874ffbc0bae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solution\n",
    "from scipy.special import exp1\n",
    "\n",
    "# Function for equation 15\n",
    "def theis(r, t, T, S, Q, t0=0):\n",
    "    u = S * r ** 2 / (4 * T * (t - t0))\n",
    "    return -Q / (4 * np.pi * T) * exp1(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd232751-9ea2-48c3-b0fa-0e2b9fbd1b5e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "r = np.linspace(rw, 100, 100)\n",
    "for i, t in enumerate([1, 2, 10]):\n",
    "    h = theis(r, t, T, S, Q)\n",
    "    plt.plot(r, h, 'C' + str(i), label=f't={t} d')\n",
    "    plt.plot(-r, h, 'C' + str(i))\n",
    "plt.xlabel('radial distance (m)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "519d891a-cac1-4833-a04a-f768b01438fa",
   "metadata": {},
   "source": [
    "### Semi-confined aquifer (Hantush)\n",
    "\n",
    "For flow in a semi-confined aquifer, the pumpted groundwater initially comes from aquifer storage but as the heads drop around the well, water from the overlying aquifer (or surface water) will start to contribute.\n",
    "\n",
    "<IMG src=\"figs/merged_nbs6_81_0.png\" width=400>\n",
    "\n",
    "The solution to this problem (for the case that $h^*=0$) is referred to as the Hantush function\n",
    "\\begin{equation}\\label{hhantush}\n",
    "h = -\\frac{Q}{4\\pi T}\\int_u^\\infty \\frac{1}{\\tau}\\exp\\left(-\\tau -\\frac{r^2}{4\\lambda^2 \\tau}\\right) \\text{d}\\tau\n",
    "\\end{equation}\n",
    "where $\\lambda=\\sqrt{Tc}$ is the leakage factor, and $u$ is the same as in the Theis solution. For the case that $h^* \\ne 0$, it can simply be added to this solution. Much has been written in the literature about how to solve the integral in the Hantush equation but in Python it can simply be calculated using the `quad` function of the `scipy.integrate` package, which we have seen before in Solution 1 of this session.\n",
    "    \n",
    "The Hantush and Theis heads are very similar for early times when most of the water pumped by the well comes from storage. The Hantush and Theis heads deviate once a significant amount of water pumped by the well starts to come from leakage through the semi-confining layer. The Theis solution never reaches steady state, but the Hantush solution does. To estimate the time it takes to reach steady state, the Hantush solution can be compared to the solution for steady flow to a well in a semi-confined aquifer, which is\n",
    "\\begin{equation}\\label{besselwell}\n",
    "h = h^* -\\frac{Q}{2\\pi T} \\text{K}_0(r/\\lambda)\n",
    "\\end{equation}\n",
    "in which K$_0$ is the modified Bessel function of the second kind and order 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbd8ca15-58c1-4e04-ba80-1f4e6c10f751",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "T = 200 # transmissivity of aquifer, m^2/d\n",
    "S = 0.0005 # storage coefficient of aquifer, -\n",
    "c = 1000 # resistance of leaky layer, d\n",
    "Q = 800 # discharge of well, m^3/d\n",
    "rw = 0.3 # radius of well, m\n",
    "lab = np.sqrt(c * T) # leakage factor, m\n",
    "print(f'leakage factor: {lab:.2f} m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ed52f04-6544-422e-aa4b-c0256f78dea8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solution\n",
    "from scipy.special import exp1, k0\n",
    "from scipy.integrate import quad\n",
    "\n",
    "def integrand(tau, r, T, lab):\n",
    "    return 1 / tau * np.exp(-tau - r ** 2 / (4 * lab ** 2 * tau))\n",
    "\n",
    "def hantush(r, t, T, S, c, Q):\n",
    "    lab = np.sqrt(T * c)\n",
    "    u = S * r ** 2 / (4 * T * t)\n",
    "    F = quad(integrand, u, np.inf, args=(r, T, lab))[0]\n",
    "    return -Q / (4 * np.pi * T) * F\n",
    "\n",
    "hantushvec = np.vectorize(hantush) # vectorize hantush function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0e77660-4d05-49ca-8da6-558239c11df7",
   "metadata": {
    "tags": [
     "hide_output"
    ]
   },
   "outputs": [],
   "source": [
    "# basic plot for one value of r and one value of t\n",
    "plt.subplot(121)\n",
    "r = np.linspace(rw, 2 * lab, 100)\n",
    "t = 6 / 24\n",
    "plt.plot(r, hantushvec(r, t, T, S, c, Q), label='Hantush')\n",
    "plt.plot(r, -Q / (2 * np.pi * T) * k0(r / lab), 'C0--', label='Steady semi-confined')\n",
    "plt.plot(r, theis(r, t, T, S, Q), label='Theis') # theis func from above\n",
    "plt.title('Head after 6 hours')\n",
    "plt.grid()\n",
    "plt.xlabel('$x$ (m)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.legend()\n",
    "plt.subplot(122)\n",
    "t = np.linspace(0.01, 2, 100)\n",
    "r = lab # evaluate head at r=lambda\n",
    "plt.plot(t, hantushvec(r, t, T, S, c, Q), label='Hantush')\n",
    "plt.axhline(-Q / (2 * np.pi * T) * k0(lab / lab), color='C0', ls='--', label='Steady semi-confined')\n",
    "plt.plot(t, theis(r, t, T, S, Q), '--', label='Theis')\n",
    "plt.title(f'Head at {lab:0.1f} m')\n",
    "plt.grid()\n",
    "plt.xlabel('$t$ (d)')\n",
    "plt.ylabel('head (m)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08036050-4239-4828-afe8-6b5be81133df",
   "metadata": {},
   "source": [
    "### Pumping test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cc0a9a6-b6de-427f-a6d7-bcf5d056edbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data at Oude Korendijk taken from Kruseman and de Ridder (1990)\n",
    "robs = 30 # distance to of observation well, m\n",
    "Q = 788 # discharge of well, m^3/d\n",
    "time = np.array(\n",
    "    [0.1 ,  0.25,  0.5 ,  0.7 ,  1.  ,  1.4 ,  1.9 ,  2.33,  2.8 ,\n",
    "     3.36,  4.  ,  5.35,  6.8 ,  8.3 ,  8.7 , 10.  , 13.1 , 18.  ,\n",
    "     27. ,  33. ,  41. ,  48. ,  59. ,  80. , 95.  , 139. , 181. , \n",
    "     245.,  300.,  360.,  480.,  600.,  728., 830.]) # in minutes\n",
    "drawdown = np.array(\n",
    "    [0.04 , 0.08 , 0.13 , 0.18 , 0.23 , 0.28 , 0.33 , 0.36 , 0.39 ,\n",
    "     0.42 , 0.45 , 0.5  , 0.54 , 0.57 , 0.58 , 0.6  , 0.64 , 0.68 ,\n",
    "     0.742, 0.753, 0.779, 0.793, 0.819, 0.855, 0.873, 0.915, 0.935,\n",
    "     0.966, 0.99 , 1.007, 1.05 , 1.053, 1.072, 1.088]) # in meters\n",
    "tobs = time / 24 / 60 # convert observation time to days\n",
    "hobs = -drawdown # convert drawdown to heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6466fc6-fb99-491e-b79a-933b0de0c413",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(121)\n",
    "plt.plot(tobs, hobs, 'C1.')\n",
    "plt.xlabel('time (d)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.title(f'measured head at r={robs} m')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0e493a5-d8d1-476f-9e2f-cd9680792ba7",
   "metadata": {},
   "source": [
    "Find optimal values for transmissivity $T$ and storage coefficient $S$ using least squares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70d38800-5523-4284-aba3-766ed01bebb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define objective function to minimize\n",
    "def fobj(p, ho, to, ro, Q):\n",
    "    T, S = np.exp(p)\n",
    "    hm = theis(ro, to, T, S, Q)\n",
    "    return np.sum((ho - hm) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a0ad72c-55ce-4225-aad0-2668718d33db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# least squares solution\n",
    "from scipy.optimize import fmin\n",
    "logp = fmin(fobj, np.log(np.array([100, 1e-4])), args=(hobs, tobs, robs, Q), disp=True)\n",
    "T, S = np.exp(logp)\n",
    "print(f'Estimated parameters T: {T:.1f} m^2/d, S: {S:.2e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bdd0a7f-12ef-46d7-9cae-8fe990ad29a6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# plot of best fit\n",
    "hm = theis(robs, tobs, T, S, Q)\n",
    "plt.subplot(121)\n",
    "plt.plot(tobs, hobs, 'C1.')\n",
    "plt.plot(tobs, hm, 'C0')\n",
    "plt.xlabel('time (d)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid()\n",
    "plt.subplot(122)\n",
    "plt.semilogx(tobs, hobs, 'C1.')\n",
    "plt.semilogx(tobs, hm, 'C0')\n",
    "plt.xlabel('time (d)')\n",
    "plt.ylabel('head (m)')\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bec16640-103a-4097-9b3e-1ac8959630bb",
   "metadata": {},
   "source": [
    "## Exercise 4.1\n",
    "\n",
    "Fit the pump test data using the Hantush formula. You can do that by trial and error (varying the values of the parameters $T$, $S$ and $c$ by hand), or by formally minimizing the sum of the squared errors (note that you may see some warnings and that it could take a while before the calculation is finished)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6127a4bd-ed31-4bb9-83af-7860800856c7",
   "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.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}