{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# 1. Leaky Aquifer Test - Dalem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import timflow.transient as tft\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = [5, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Introduction and Conceptual Model\n",
    "This example is a pumping test from Dalem, the Netherlands (Kruseman and de Ridder, 1970). An aquifer of 37 m thickness is overlain by a semi-confining layer of 8 m thick. Water is pumped from the aquifer and drawdown is recorded at four different piezometers, 30, 60, 90 and 120 m away from the well. The pumping lasted 8 hours in total at a rate of 761 m$^3$/d. There is a river 1500 m away from the well. The tide affects both river and well levels. Data was corrected for the tide effect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "source": [
    "<img src=\"./figs/Dalem.png\" style=\"width:400pt\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# data of observation well 30 m away from pumping well\n",
    "data1 = np.loadtxt(\"data/dalem_p30.txt\", skiprows=1)\n",
    "t1 = data1[:, 0]\n",
    "h1 = data1[:, 1]\n",
    "\n",
    "# data of observation well 60 m away from pumping well\n",
    "data2 = np.loadtxt(\"data/dalem_p60.txt\", skiprows=1)\n",
    "t2 = data2[:, 0]\n",
    "h2 = data2[:, 1]\n",
    "\n",
    "# data of observation well 90 m away from pumping well\n",
    "data3 = np.loadtxt(\"data/dalem_p90.txt\", skiprows=1)\n",
    "t3 = data3[:, 0]\n",
    "h3 = data3[:, 1]\n",
    "\n",
    "# data of observation well 120 m away from pumping well\n",
    "data4 = np.loadtxt(\"data/dalem_p120.txt\", skiprows=1)\n",
    "t4 = data4[:, 0]\n",
    "h4 = data4[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Parameters and model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# known parameters\n",
    "H = 37  # aquifer thickness in m\n",
    "zt = -8  # top boundary of aquifer in m\n",
    "zb = zt - H  # bottom boundary of the aquifer in m\n",
    "Q = 761  # constant pumping rate in m^3/d\n",
    "t = 8 / 24  # total pumping time in days (8 hours = 0.34 days)\n",
    "\n",
    "r1 = 30  # distance from pumping well to observation well 1 in m\n",
    "r2 = 60  # distance from pumping well to observation well 2 in m\n",
    "r3 = 90  # distance from pumping well to observation well 3 in m\n",
    "r4 = 120  # distance from pumping well to observation well 4 in m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The `timflow` model consists of a single semi-confined aquifer and a pumping well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# unkonwn parameters: kaq, Saq, c\n",
    "ml = tft.ModelMaq(\n",
    "    kaq=10, z=[0, zt, zb], c=500, Saq=1e-4, topboundary=\"semi\", tmin=0.01, tmax=1\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, tsandQ=[(0, Q)])\n",
    "ml.solve(silent=\"True\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Estimate aquifer parameters\n",
    "The hydraulic condivity and specific storage coefficient of the aquifer are estimated and the resistance of the leaky semi-confining layer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "cal = tft.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq\", initial=10, layers=0)\n",
    "cal.set_parameter(name=\"Saq\", initial=1e-4, layers=0)\n",
    "cal.set_parameter(name=\"c\", initial=500, pmin=0, layers=0)\n",
    "\n",
    "cal.series(name=\"obs1\", x=30, y=0, t=t1, h=h1, layer=0)\n",
    "cal.series(name=\"obs2\", x=60, y=0, t=t2, h=h2, layer=0)\n",
    "cal.series(name=\"obs3\", x=90, y=0, t=t3, h=h3, layer=0)\n",
    "cal.series(name=\"obs4\", x=120, y=0, t=t4, h=h4, layer=0)\n",
    "cal.fit(report=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "display(cal.parameters.loc[:, [\"optimal\"]])\n",
    "print(f\"RMSE: {cal.rmse():.5f} m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "hm_1 = ml.head(r1, 0, t1)\n",
    "plt.semilogx(t1, h1, \".\", label=\"obs at 30 m\")\n",
    "plt.semilogx(t1, hm_1[0], label=\"timflow at 30 m\")\n",
    "\n",
    "hm_2 = ml.head(r2, 0, t2)\n",
    "plt.semilogx(t2, h2, \".\", label=\"obs at 60 m\")\n",
    "plt.semilogx(t2, hm_2[0], label=\"timflow at 60 m\")\n",
    "\n",
    "hm_3 = ml.head(r3, 0, t3)\n",
    "plt.semilogx(t3, h3, \".\", label=\"obs at 90 m\")\n",
    "plt.semilogx(t3, hm_3[0], label=\"timflow at 90 m\")\n",
    "\n",
    "hm_4 = ml.head(r4, 0, t4)\n",
    "plt.semilogx(t4, h4, \".\", label=\"obs at 120 m\")\n",
    "plt.semilogx(t4, hm_4[0], label=\"timflow at 120 m\")\n",
    "\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head change (m)\")\n",
    "plt.title(\"Model Results\")\n",
    "plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Estimate aquifer parameters using weights\n",
    "The head change is smaller in the observation wells that are farther away from the pumping well. It may possibly be expected that the difference between the modeled and measured values is smaller when the total drawdown is smaller. This may be taken into account by giving the measurements weights. In the calibration below, the weight for each observation well is equal to one divided by the maximum absolute head change in that overservation well. It is up to the modeler to decide whether such a weighing is desirable. Note that the root mean squared error (computed without weights) is a bit larger than for the case without the weights, but the increase is small. The fitted values of $k$ and $S_s$ are similar to the case without weights, but the resistance $c$ is quite a bit larger. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "cal2 = tft.Calibrate(ml)\n",
    "cal2.set_parameter(name=\"kaq\", initial=10, layers=0)\n",
    "cal2.set_parameter(name=\"Saq\", initial=1e-4, layers=0)\n",
    "cal2.set_parameter(name=\"c\", initial=500, pmin=0, layers=0)\n",
    "\n",
    "cal2.series(name=\"obs1\", x=30, y=0, t=t1, h=h1, layer=0, weights=1 / np.max(np.abs(h1)))\n",
    "cal2.series(name=\"obs2\", x=60, y=0, t=t2, h=h2, layer=0, weights=1 / np.max(np.abs(h2)))\n",
    "cal2.series(name=\"obs3\", x=90, y=0, t=t3, h=h3, layer=0, weights=1 / np.max(np.abs(h3)))\n",
    "cal2.series(name=\"obs4\", x=120, y=0, t=t4, h=h4, layer=0, weights=1 / np.max(np.abs(h4)))\n",
    "cal2.fit(report=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "display(cal2.parameters[[\"optimal\"]])\n",
    "print(f\"RMSE: {cal2.rmse(weighted=False):.3f} m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Comparison of results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The performance of `timflow` is compared with other models and with parameter values reported in Kruseman and de Ridder (1970). The published values were obtained by graphical matching to a Hantush type curve (Hantush, 1955). In addition to the `timflow` results and the literature values, results from AQTESOLV (Duffield, 2007) and MLU (Hemker & Post, 2014) are included for comparison. Both AQTESOLV and MLU were applied to minimize the sum of the squared log-transformed drawdowns. AQTESOLV only fitted the model to the observation well at r = 90 m, whereas the other models were calibrated using data from multiple observation wells.\n",
    "\n",
    "Overall, all models yield similar estimates of the hydraulic conductivity and specific storage coefficient. The fitted value of the resistance differs. However, the `timflow` results differ substantially from the MLU and AQTESOLV solutions with respect to aquitard storage and hydraulic resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "t = pd.DataFrame(\n",
    "    columns=[\"k [m/d]\", \"Ss [1/m]\", \"c [d]\", \"RMSE [m]\"],\n",
    "    index=[\"timflow\", \"timflow weights\", \"AQTESOLV\", \"MLU\", \"Hantush\"],\n",
    ")\n",
    "\n",
    "t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
    "t.loc[\"timflow weights\"] = np.append(\n",
    "    cal2.parameters[\"optimal\"].values, cal2.rmse(weighted=False)\n",
    ")\n",
    "t.loc[\"AQTESOLV\"] = [41.5575, 4.4625e-05, 327.920, 0.03818]\n",
    "t.loc[\"MLU\"] = [45.297, 4.865e-05, 328, 0.042426]\n",
    "t.loc[\"Hantush\"] = [45.332, 4.762e-5, 331.141, 0.005917]\n",
    "\n",
    "t_formatted = t.style.format(\n",
    "    {\"k [m/d]\": \"{:.2f}\", \"Ss [1/m]\": \"{:.2e}\", \"c [d]\": \"{:.0f}\", \"RMSE [m]\": \"{:.4f}\"}\n",
    ")\n",
    "t_formatted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## References\n",
    "\n",
    "* Bakker, M. (2013), Semi-analytic modeling of transient multi-layer flow with TTim, Hydrogeol J 21, 935–943, https://doi.org/10.1007/s10040-013-0975-2 \n",
    "* Duffield, G.M., (2007), AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n",
    "* Hemker, K. en Post V. (2014) MLU for Windows: well flow modeling in multilayer aquifer systems; MLU User's guide. https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmicrofem.com%2Fdownload%2Fmlu-user.pdf&data=05%7C02%7CMark.Bakker%40tudelft.nl%7Cad7f16364d2d4fd55dbf08de73832eaa%7C096e524d692940308cd38ab42de0887b%7C0%7C0%7C639075204580287861%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=OBoe8seXZUfoat89Dfr4g6lF%2Bn1FdtXqtp%2F18BMXCn0%3D&reserved=0\n",
    "* Kruseman, G.P., De Ridder, N.A. and Verweij, J.M., (1970), Analysis and evaluationof pumping test data, volume 11, International institute for land reclamation and improvement The Netherlands.\n",
    "* Newville, M., Stensitzki, T., Allen, D.B. and Ingargiola, A. (2014), LMFIT: Non Linear Least-Squares Minimization and Curve Fitting for Python, https://dx.doi.org/10.5281/zenodo.11813, https://lmfit.github.io/lmfit-py/intro.html (last access: August,2021)."
   ]
  },
  {
   "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.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}