{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Leaky Aquifer Test - Dalem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "source": [
    "### Introduction and Conceptual Model\n",
    "This example is a pumping test from Dalem (Kruseman et al., 1970), the Netherlands. The hydrogeological cross-section is composed of the following elements: \n",
    "* an initial 8 m deep aquitard layer,\n",
    "* followed by an aquifer from 8 m to 45 m depth,\n",
    "* the layer underlying the aquifer is considered an aquiclude.\n",
    "\n",
    "The pumping well is placed at 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 has been previously corrected for the tide effect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "source": [
    "<img src=\"./figs/Dalem_model1.png\" style=\"width:400pt\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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]  # days\n",
    "h1 = data1[:, 1]  # m\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": {},
   "source": [
    "### Parameters and model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": {},
   "source": [
    "In many situations, we cannot ignore the leakage potential of overlying and underlying formations to aquifers, and we cannot conceptualize them as confined. `timflow` is capable of modelling and adjusting parameters to leaky, semi-confined aquifers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unknown parameters: kaq, Saq, c\n",
    "ml = tft.ModelMaq(\n",
    "    kaq=10, z=[0, zt, zb], c=500, Saq=0.001, topboundary=\"semi\", tmin=0.01, tmax=1\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, tsandQ=[(0, Q), (0.34, 0)])\n",
    "ml.solve(silent=\"True\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate aquifer parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(cal.parameters)\n",
    "print(\"RMSE:\", cal.rmse())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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(\"drawdown(m)\")\n",
    "plt.title(\"Model Results\")\n",
    "plt.legend(bbox_to_anchor=(1.05, 1))\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The performance of `timflow` was evaluated by comparison with other numerical models and with parameter values reported in Kruseman and de Ridder (1970). The published values were obtained by graphical matching to the Hantush family of type curves (Hantush, 1955). In addition to the `timflow` results and the literature values, results from AQTESOLV (Duffield, 2007) and MLU (Carlson & Randall, 2012), as reported by Yang (2020), are included for comparison.\n",
    "\n",
    "Overall, all models yield similar estimates of the aquifer hydraulic conductivity. 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": {},
   "outputs": [],
   "source": [
    "t = pd.DataFrame(\n",
    "    columns=[\"k [m/d]\", \"Ss [1/m]\", \"c [d]\", \"RMSE [m]\"],\n",
    "    index=[\"timflow\", \"AQTESOLV\", \"MLU\", \"Hantush\"],\n",
    ")\n",
    "\n",
    "t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
    "t.loc[\"AQTESOLV\"] = [49.286, 4.559e-05, 745.156, 0.007245]\n",
    "t.loc[\"MLU\"] = [45.186, 3.941e-05, 769.200, 0.005941]\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]\": \"{:.2f}\"}\n",
    ")\n",
    "t_formatted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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",
    "* Carlson, F. and Randall, J. (2012), MLU: a Windows application for the analysis of aquifer tests and the design of well fields in layered systems, Ground Water 50(4):504–510.\n",
    "* Duffield, G.M., (2007), AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\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)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}