{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# 2. Confined Aquifer Test - Gridley"
   ]
  },
  {
   "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)  # default figure size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction and Conceptual Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "This example is a pumping test conducted in 1953 in Gridley, Illinois, US. It was reported by Walton (1962). The aquifer is an 18 ft thick sand and gravel layer under confined conditions. The pumping well fully penetrates the formation, and pumping was conducted for 8 hours at a rate of 220 gallons per minute. The effect of pumping was observed at observation well 1, located 824 ft away from the well.\n",
    "\n",
    "The time-drawdown data for the observation well were obtained from the AQTESOLV documentation (Duffield, 2007), while data from the pumping well were obtained from the original paper from Walton (1962). Following AQTESOLV documentation (Duffield, 2007), radii of 0.5 ft were used for both the pumping and observation wells. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figs/gridley.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": [
    "# Loading Observation well (Well 1)\n",
    "data1 = np.loadtxt(\"data/gridley_well_1.txt\")\n",
    "to1 = data1[:, 0]\n",
    "ho1 = data1[:, 1]\n",
    "\n",
    "# Loading Pumping Well data (Well 3)\n",
    "data2 = np.loadtxt(\"data/gridley_well_3.txt\")\n",
    "to2 = data2[:, 0]\n",
    "ho2 = data2[:, 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 = 18 * 0.3048  # aquifer thickness in m (18 ft = 5.48 m)\n",
    "Q = (\n",
    "    (220 * 0.00378541) * 60 * 24\n",
    ")  # constant discharge in m^3/d (220 gallons/minute = 1199.22 m^3/d)\n",
    "r = (\n",
    "    824 * 0.3048\n",
    ")  # distance between observation well to test well in m (824 ft = 251.16 m)\n",
    "rw = 0.5 * 0.3048  # screen radius of test well in m (0.5 ft = 0.15 m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Wellbore storage is added to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=10, z=[0, -H], Saq=0.001, tmin=0.001, tmax=1, topboundary=\"conf\")\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=rw, tsandQ=[(0, Q)], rc=rw, layers=0)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Estimate aquifer parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The aquifer parameters are estimated using both datasets simultaneously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unknown parameters: kaq, Saq\n",
    "cal = tft.Calibrate(ml)  # create Calibrate object\n",
    "cal.set_parameter(\n",
    "    name=\"kaq\", initial=10, layers=0\n",
    ")  # setting the parameters for calibration\n",
    "cal.set_parameter(name=\"Saq\", initial=1e-4, pmin=1e-7, layers=0)\n",
    "cal.series(name=\"obs1\", x=r, y=0, t=to1, h=ho1, layer=0)  # setting the observation data\n",
    "cal.seriesinwell(name=\"obs2\", element=w, t=to2, h=ho2)  # setting the observation data\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": [
    "hm1 = ml.head(r, 0, to1)\n",
    "hm2 = w.headinside(to2)\n",
    "plt.semilogx(to1, hm1[0], \"C0\", label=\"timflow piezometer\")\n",
    "plt.semilogx(to1, ho1, \"C0.\", label=\"obs piezometer\")\n",
    "plt.semilogx(to2, hm2[0], \"C1\", label=\"timflow well\")\n",
    "plt.semilogx(to2, ho2, \"C1.\", label=\"obs well\")\n",
    "plt.title(\"Model Results\")\n",
    "plt.xlabel(\"time [d]\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "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` with two datasets simultaneously was evaluated by comparison with AQTESOLV (Duffield, 2007), and MLU (Carlson and Randall, 2012). Results from `timflow` with added wellbore storage and MLU are identical, while those from AQTESOLV show small deviations."
   ]
  },
  {
   "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]\", \"RMSE [m]\"], index=[\"timflow\", \"AQTESOLV\", \"MLU\"]\n",
    ")\n",
    "\n",
    "t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
    "t.loc[\"MLU\"] = [38.094, 1.193e-06, 0.259]\n",
    "t.loc[\"AQTESOLV\"] = [37.803, 1.356e-06, 0.270]\n",
    "\n",
    "t_formatted = t.style.format(\n",
    "    {\"k [m/d]\": \"{:.2f}\", \"Ss [1/m]\": \"{:.2e}\", \"RMSE [m]\": \"{:.2f}\"}\n",
    ")\n",
    "t_formatted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "* 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",
    "* Walton, W.C. (1962), Selected analytical methods for well and aquifer evaluation, Illinois, department of Registration & Education.bulletin 49."
   ]
  }
 ],
 "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
}