{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# 2. Slug Test - Falling Head"
   ]
  },
  {
   "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",
    "\n",
    "This slug test, taken from the AQTESOLV examples (Duffield, 2007), was reported in Batu (1998). \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."
   ]
  },
  {
   "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 in m (5 inch = 0.127 m)\n",
    "rc = 2 * 0.0254  # well casing radius in m (2 inch = 0.0508 m)\n",
    "L = 13.8 * 0.3048  # screen length in m (13.8 ft = 4.2 m)\n",
    "b = -32.57 * 0.3048  # aquifer thickness in m (-32.57 ft = -9.92 m)\n",
    "zt = -0.47 * 0.3048  # depth to top of the screen in m (-0.47 ft = -0.14 m)\n",
    "H0 = 1.48 * 0.3048  # initial displacement in the well in m (1.48 ft = 0.4511 m)\n",
    "zb = zt - L  # bottom of the screen in 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": [
    "# The thickness of each layer is set to be 0.5 m\n",
    "z0 = np.arange(zt, zb, -0.5)\n",
    "z1 = np.arange(zb, b, -0.5)\n",
    "zlay = np.append(z0, z1)\n",
    "zlay = np.append(zlay, b)\n",
    "zlay = np.insert(zlay, 0, 0)\n",
    "nlay = len(zlay) - 1  # number of layers\n",
    "Saq = 1e-4 * np.ones(nlay)\n",
    "Saq[0] = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.Model3D(\n",
    "    kaq=10, z=zlay, Saq=Saq, kzoverkh=1, tmin=1e-5, tmax=0.01, phreatictop=True\n",
    ")\n",
    "w = tft.Well(\n",
    "    ml,\n",
    "    xw=0,\n",
    "    yw=0,\n",
    "    rw=rw,\n",
    "    tsandQ=[(0, -Q)],\n",
    "    layers=[1, 2, 3, 4, 5, 6, 7, 8],\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(nlay)), initial=10, pmin=0)\n",
    "cal.set_parameter(name=\"Saq\", layers=list(range(nlay)), initial=1e-4, pmin=0)\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 = 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, 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). AQTESOLV parameters are quite different from the set parameters in `timflow`. Furthermore, AQTESOLV also has a better RMSE performance."
   ]
  },
  {
   "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-multi\", \"AQTESOLV\"],\n",
    ")\n",
    "\n",
    "t.loc[\"timflow-multi\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
    "t.loc[\"AQTESOLV\"] = [2.616, 7.894e-5, 0.001197]\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",
    "\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
}