{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1974f9bb-529b-4788-9fe9-cbfe24e19ab5",
   "metadata": {},
   "source": [
    "## Well near a partially penetrating river with leaky bed - comparison with exact solution of Hunt (1999)\n",
    "Consider flow to a well near a river with a leaky bed. The well is located at $(x,y)=(d, 0)$ and starts pumping at a discharge $Q$ at time $t=0$. The river is infinitely long and runs along $x=0$. The head in the river is fixed (i.e., the drawdown in the river is zero). The width of the river is $b$ and the resistance of the streambed is $c$. The width of the river is neglected in the model but taken into account when computing the riverbed conductance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4a26698-6ac6-4dba-ab29-21c00ed4e6a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import quad\n",
    "from scipy.special import exp1\n",
    "\n",
    "import timflow.transient as tft\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (6, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3211ac6a-9508-4b42-917f-b50df7d7687d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "k = 10  # hydraulic conductivity, m/d\n",
    "H = 10  # thickness of aquifer, m\n",
    "S = 0.1  # phreatic storage, -\n",
    "c = 5  # streambed resistance, d\n",
    "b = 5  # stream width, m\n",
    "d = 100  # distance of well from river, m\n",
    "Q = 500  # well discharge, m^3/d\n",
    "# computed parameters\n",
    "T = k * H\n",
    "C = b / c  # riverbed conductance, = lambda in Hunt (1999)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa2e2d5d-0906-4e63-b420-c8020bd5bb24",
   "metadata": {},
   "source": [
    "The exact solution is given in Hunt (1999), Equation 30. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1e40642-bfb8-4fec-a485-0d0595937c40",
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(theta, x, y, t):\n",
    "    rv = np.exp(-theta) * exp1(\n",
    "        ((d + np.abs(x) + 2 * T * theta / C) ** 2 + y**2) / (4 * T * t / S)\n",
    "    )\n",
    "    return rv\n",
    "\n",
    "\n",
    "def hunt(x, y, t):\n",
    "    part1 = exp1(((d - x) ** 2 + y**2) / (4 * T * t / S))\n",
    "    part2 = quad(func, 0, np.inf, args=(x, y, t))[0]\n",
    "    return -Q / (4 * np.pi * T) * (part1 - part2)\n",
    "\n",
    "\n",
    "huntvec = np.vectorize(hunt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4340897f-86e6-43fe-ab0d-29c154826c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in [1, 20, 200, 500]:\n",
    "    nx = 100\n",
    "    x = np.linspace(-d, 2 * d, nx)\n",
    "    y = np.zeros(nx)\n",
    "    h = huntvec(x, y, t)\n",
    "    plt.plot(x, h, label=f\"time={t} d\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.title(\"head along x=0\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48e74205-50ae-42a0-8ec2-4de8a0eca660",
   "metadata": {},
   "source": [
    "### Timflow model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4d7069e-2efc-4213-b224-1dbb3acd808a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=k, z=[0, -10], Saq=S, phreatictop=True, tmin=1, tmax=1000)\n",
    "w = tft.Well(ml, d, 0, tsandQ=[(0, Q)])\n",
    "yls = np.hstack(\n",
    "    (\n",
    "        np.linspace(-1000, -200, 10),\n",
    "        np.linspace(-200, 200, 10)[1:],\n",
    "        np.linspace(200, 1000, 10)[1:],\n",
    "    )\n",
    ")\n",
    "xls = np.zeros(len(yls))\n",
    "xy = np.array([xls, yls]).T\n",
    "riv = tft.RiverString(ml, xy=xy, tsandh=\"fixed\", res=c, wh=b)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a73a0f51-a358-4d72-8ffd-c68eb6f5a996",
   "metadata": {},
   "source": [
    "### Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6097091-3af2-4b25-9af6-a89017347d6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-100, 200, 101)\n",
    "y = np.zeros(101)\n",
    "t = 20\n",
    "hhunt = huntvec(x, y, t)\n",
    "plt.plot(x, hhunt, label=\"hunt\")\n",
    "htim = ml.headalongline(x, y, t)\n",
    "plt.plot(x, htim[0, 0], \"--\", label=\"timflow\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.title(f\"t={t} d\")\n",
    "plt.axvline(0, color=\"k\", ls=\"--\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c300372-ad31-4f81-b214-f4aa2285f0dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.meshgrid(np.linspace(-100, 200, 50), np.linspace(0, 100, 50))\n",
    "t = 20\n",
    "h = huntvec(x, y, t)\n",
    "plt.subplot(111, aspect=1)\n",
    "cs = plt.contour(x, y, h, levels=np.arange(-4, 0, 0.2), colors=\"C0\")\n",
    "plt.clabel(cs, fmt=\"%1.1f\")\n",
    "x, y = np.linspace(-100, 200, 50), np.linspace(-100, 0, 50)\n",
    "h = ml.headgrid(x, y, t).squeeze()\n",
    "cs = plt.contour(x, y, h, levels=np.arange(-4, 0, 0.2), colors=\"C1\")\n",
    "plt.clabel(cs, fmt=\"%1.1f\")\n",
    "plt.axvline(0, color=\"k\", ls=\"--\")\n",
    "plt.text(-80, 80, \"hunt\")\n",
    "plt.text(-80, -80, \"timflow\")\n",
    "plt.title(f\"t={t} d\")\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"y (m)\")\n",
    "plt.plot(d, 0, \"k.\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d6db098-6a6b-4c50-89df-596ee48f38c9",
   "metadata": {},
   "source": [
    "### Reference\n",
    "Hunt, B., 1999. Unsteady stream depletion from ground water pumping. Groundwater, 37(1), pp.98-102, [https://doi.org/10.1111/j.1745-6584.1999.tb00962.x](https://doi.org/10.1111/j.1745-6584.1999.tb00962.x)"
   ]
  }
 ],
 "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": 5
}