{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# 3. Slug Test - Falling Head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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",
    "\n",
    "A well partially penetrates a sandy unconfined aquifer that has a saturated depth of 32.57 ft. The top of the screen is located 0.47 ft below the water table and has 13.8 ft in length. The well and casing radii are 5 and 2 inches, respectively. The slug displacement is 1.48 ft. Head change has been recorded at the slug well. This slug test, taken from the AQTESOLV examples (Duffield, 2007), was reported in Batu (1998). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figs/Falling_head.png\" style=\"width:400pt\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt(\"data/falling_head.txt\", skiprows=2)\n",
    "to = data[:, 0] / 60 / 60 / 24  # convert time from seconds to days\n",
    "ho = (10 - data[:, 1]) * 0.3048  # convert drawdown from ft to meters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters and model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rw = 5 * 0.0254  # well radius, m\n",
    "rc = 2 * 0.0254  # well casing radius, m\n",
    "L = 13.8 * 0.3048  # screen length, m\n",
    "b = 32.57 * 0.3048  # aquifer thickness, m\n",
    "zt = 0  # top of aquifer, m\n",
    "zst = -0.47 * 0.3048  # top of screen, m\n",
    "zsb = zst - L  # bottom of screen, m\n",
    "zb = zt - b  # bottom of aquifer, m\n",
    "H0 = 1.48 * 0.3048  # initial displacement in well, m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "convert measured displacement into volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = np.pi * rc**2 * H0\n",
    "print(f\"slug: {Q:.5f} m^3\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create a multi-layer model. For this, we divide the second and third layers into 0.5 m thick layers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# z = np.array([zt, zt - 0.1, zst, zsb, zb])\n",
    "z = np.hstack(\n",
    "    (\n",
    "        zt,\n",
    "        np.linspace(zt - 0.01, zst, 5)[:-1],\n",
    "        np.linspace(zst, zsb, 5)[:-1],\n",
    "        np.linspace(zsb, zb, 10),\n",
    "    )\n",
    ")\n",
    "nlay = len(z) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.Model3D(\n",
    "    kaq=10,\n",
    "    z=z,\n",
    "    Saq=[0.1] + (nlay - 1) * [1e-4],\n",
    "    kzoverkh=1,\n",
    "    tmin=1e-5,\n",
    "    tmax=0.01,\n",
    "    topboundary=\"phreatic\",\n",
    ")\n",
    "w = tft.Well(\n",
    "    ml,\n",
    "    xw=0,\n",
    "    yw=0,\n",
    "    rw=rw,\n",
    "    tsandQ=[(0, -Q)],\n",
    "    layers=range(6, 6 + 5),\n",
    "    rc=rc,\n",
    "    wbstype=\"slug\",\n",
    ")\n",
    "ml.solve()"
   ]
  },
  {
   "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\", layers=list(range(1, nlay)), initial=10, pmin=0)\n",
    "cal.set_parameter(name=\"Saq\", layers=list(range(1, nlay)), initial=1e-4, pmin=0)\n",
    "cal.seriesinwell(name=\"obs\", element=w, t=to, h=ho)\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():.3f} m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tm = np.logspace(np.log10(to[0]), np.log10(to[-1]), 100)\n",
    "hm = w.headinside(tm)\n",
    "plt.semilogx(to, ho / H0, \".\", label=\"obs\")\n",
    "plt.semilogx(tm, hm[0] / H0, label=\"timflow\")\n",
    "plt.xlabel(\"time [d]\")\n",
    "plt.ylabel(\"Normalized head (h/H0)\")\n",
    "plt.title(\"Model results - multi-layer model\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of results\n",
    "Here, the `timflow` performance in analysing slug tests is checked. The solution in `timflow` is compared with the KGS analytical model (Hyder et al. 1994) implemented in AQTESOLV (Duffield, 2007). The parameters of `timflow` and AQTESOLV are similar, even though AQTESOLV only used one layer."
   ]
  },
  {
   "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]\"],\n",
    "    index=[\"timflow\", \"AQTESOLV\"],\n",
    ")\n",
    "\n",
    "t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
    "t.loc[\"AQTESOLV\"] = [0.4211, 5.70e-4, \"-\"]\n",
    "\n",
    "t_formatted = t.style.format(\n",
    "    {\n",
    "        \"k [m/d]\": \"{:.2f}\",\n",
    "        \"Ss [1/m]\": \"{:.2e}\",\n",
    "        \"RMSE [m]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.3f}\",\n",
    "    }\n",
    ")\n",
    "t_formatted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## References\n",
    "\n",
    "* Batu, V. (1998), Aquifer hydraulics: a comprehensive guide to hydrogeologic data analysis, John Wiley & Sons\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. and 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
}