{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Slug Test - Pratt County"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This slug test, taken from AQTESOLV examples (Duffield, 2007), was conducted in Pratt County Monitoring Site, US, and reported by Butler (1998). \n",
    "\n",
    "A partially penetrating well is screened in unconsolidated alluvial deposits, consisting of sand and gravel interbedded by clay. The total thickness of the aquifer is 47.87 m. The screen is located at 16.77 m depth and has a screen length 1.52 m the well radius is 0.125 m, and the casing radius 0.064 m. The slug displacement is 0.671 m. Head change has been recorded at the slug well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figs/Pratt_County.png\" style=\"width:400pt\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt(\"data/slug.txt\", skiprows=1)\n",
    "to = data[:, 0] / 60 / 60 / 24  # convert time in seconds to days\n",
    "ho = data[:, 1]  # m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters and model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rw = 0.125  # well radius in m\n",
    "rc = 0.064  # well casing radius in m\n",
    "L = 1.52  # screen length in m\n",
    "b = -47.87  # aquifer thickness in m\n",
    "zt = -16.77  # depth to top of screen in m\n",
    "H0 = 0.671  # initial displacement in the well in m\n",
    "zb = zt - L  # bottom of screen in m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The input for the slug test in `timflow` is the added or removed volume. Therefore, we must first convert our measured displacement into volume."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = np.pi * rc**2 * H0\n",
    "print(\"slug:\", round(Q, 5), \"m^3\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We conceptualize the aquifer as a three-layer model, one layer above the screen, one layer at the screen top and bottom and another layer just below it.\n",
    "\n",
    "The setting of the slug well is slightly different from the pumping well. We detail the differences below:\n",
    "* the ```tsandQ``` argument in the ```Well``` object has a different meaning. Instead of meaning the time of start or shutdown and the pumping rate of the pumping well, it means the time a volume is instantaneously added or removed from the well. In our case, we defined it as ```[(0, -Q)]``` where ```0``` is the moment in time when we added the slug and ```Q``` is the added volume. A negative sign means a volume is added. Otherwise, it would mean an extracted volume.\n",
    "* the ```wbstype``` argument is set to ```' slug'```, so `timflow` knows the ```tsandQ``` argument means time and instant volumes, instead of pumping rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.Model3D(kaq=10, z=[0, zt, zb, b], Saq=1e-4, kzoverkh=1, tmin=1e-6, tmax=0.01)\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, -Q)], layers=1, wbstype=\"slug\")\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate aquifer parameters "
   ]
  },
  {
   "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\", layers=[0, 1, 2], initial=10)\n",
    "cal.set_parameter(name=\"Saq\", layers=[0, 1, 2], initial=1e-4)\n",
    "cal.seriesinwell(name=\"obs\", element=w, t=to, h=ho)\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": [
    "tm = np.logspace(np.log10(to[0]), np.log10(to[-1]), 100)\n",
    "hm = ml.head(0, 0, tm, layers=1)\n",
    "plt.semilogx(to, ho / H0, \".\", label=\"obs\")\n",
    "plt.semilogx(tm, hm[-1] / H0, \"r\", label=\"timflow\")\n",
    "plt.ylim([0, 1])\n",
    "plt.xlabel(\"time [d]\")\n",
    "plt.ylabel(\"h/H0\")\n",
    "plt.title(\"Model Results - Three layers\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we reproduce the work of Yang (2020) to check the `timflow` performance in analysing slug tests. We compare the solution in `timflow` with the KGS analytical model (Hyder et al. 1994) implemented in AQTESOLV (Duffield, 2007). Both models show similarly low RMSE values. However, the estimated hydraulic conductivity and specific storage parameters differ substantially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = pd.DataFrame(\n",
    "    columns=[\"k [m/d]\", \"Ss [1/m]\", \"RMSE [m]\"],\n",
    "    index=[\"timflow\", \"AQTESOLV\"],\n",
    ")\n",
    "\n",
    "t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values.tolist(), cal.rmse())\n",
    "t.loc[\"AQTESOLV\"] = [4.034, 3.834e-04, 0.002976]\n",
    "\n",
    "t_formatted = t.style.format(\n",
    "    {\"k [m/d]\": \"{:.2f}\", \"Ss [1/m]\": \"{:.2e}\", \"RMSE [m]\": \"{:.3f}\"}\n",
    ")\n",
    "t_formatted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## References\n",
    "* Butler Jr., J.J. (1998), The Design, Performance, and Analysis of Slug Tests, Lewis Publishers, Boca Raton, Florida, 252p.\n",
    "* Duffield, G.M. (2007), AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n",
    "* Hyder, Z., Butler Jr, J.J., McElwee, C.D., Liu, W. (1994), Slug tests in partially penetrating wells, Water Resources Research 30, 2945–2957."
   ]
  }
 ],
 "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
}