{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Connected Wells\n",
    "\n",
    "WellString elements can be used to model a series of connected wells. There are several types of\n",
    "WellString:\n",
    "\n",
    "- `WellString`: connected wells with equal head inside the wells and a user-specified total discharge\n",
    "\n",
    "This notebook shows how to model the `WellString` element with two wells\n",
    "and compares the results to models with two individual wells, which should give the same answer for this simple case with only two wells pumping in an infinite aquifer. Examples are shown for both\n",
    "elements in single- and multi-layer models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow\n",
    "import timflow.transient as tft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timflow.show_versions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WellString"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A series of wells that have equal head (pressure) inside the well and pump with a total specified discharge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single-layer model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters\n",
    "kh = 10  # m/day\n",
    "Ss = 1e-4  # 1/m\n",
    "\n",
    "ctop = 1000.0  # resistance top leaky layer in days\n",
    "\n",
    "ztop = 0.0  # surface elevation\n",
    "zbot = -20.0  # bottom elevation of the model\n",
    "\n",
    "z = [1.0, ztop, zbot]\n",
    "kaq = np.array([kh])\n",
    "Ss = np.array([Ss])\n",
    "c = np.array([ctop])\n",
    "\n",
    "Qtot = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reference model with 2 wells with equal discharge of `Qtot / 2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlref = tft.ModelMaq(kaq=kaq, Saq=Ss, c=c, z=z, topboundary=\"semi\", tmin=0.1, tmax=10)\n",
    "w1 = tft.Well(mlref, 0, -10, tsandQ=(0, Qtot / 2), rw=0.1)\n",
    "w2 = tft.Well(mlref, 0, 10, tsandQ=(0, Qtot / 2), rw=0.1)\n",
    "mlref.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model with a WellString that has total discharge equal to the sum of the two discharge-specified wells. This should give the same solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=kaq, Saq=Ss, c=c, z=z, topboundary=\"semi\", tmin=0.1, tmax=10)\n",
    "ws = tft.WellString(ml, xy=[(0, -10), (0, 10)], tsandQ=(0, Qtot), rw=0.1)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contour the heads of the first reference model and the WellString model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = 0.2\n",
    "levels = 10\n",
    "ax = mlref.plots.topview(win=(-50, 50, -50, 50))\n",
    "mlref.plots.contour(\n",
    "    win=(-50, 50, -50, 50),\n",
    "    ngr=51,\n",
    "    t=time,\n",
    "    levels=levels,\n",
    "    decimals=2,\n",
    "    layers=[0],\n",
    "    ax=ax,\n",
    ")\n",
    "ml.plots.contour(\n",
    "    win=(-50, 50, -50, 50),\n",
    "    ngr=51,\n",
    "    t=time,\n",
    "    levels=levels,\n",
    "    decimals=2,\n",
    "    layers=[0],\n",
    "    color=\"C1\",\n",
    "    ax=ax,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare heads along y between the two models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linspace(-50, 50, 101)\n",
    "x = np.zeros_like(y)\n",
    "href = mlref.headalongline(x, y, t=time)\n",
    "h0 = ml.headalongline(x, y, t=time)\n",
    "\n",
    "plt.figure(figsize=(10, 2))\n",
    "plt.plot(y, href[0, 0], label=\"Reference Wells\")\n",
    "plt.plot(y, h0[0, 0], \"k.\", ms=3, label=\"WellString model\")\n",
    "plt.xlabel(\"y [m]\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.legend(loc=(0, 1), frameon=False, ncol=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check computed discharges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 discharge wells, Q specified not computed\n",
    "w1.discharge(time), w2.discharge(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WellString\n",
    "ws.discharge_per_well(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare heads inside the wells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 discharge wells\n",
    "w1.headinside(time), w2.headinside(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WellString\n",
    "ws.headinside(time), ws.wlist[0].headinside(time), ws.wlist[1].headinside(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multilayer model\n",
    "\n",
    "The following example compares the WellString element to the reference models with 2\n",
    "individual Wells in an aquifersystem consisting of 3 layers. All wells are screened in the bottom two aquifers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters\n",
    "kh = [5, 10, 20]  # m/day\n",
    "Ss = 1e-4  # 1/m\n",
    "\n",
    "c = [1000.0, 100.0, 10.0]  # resistance leaky layers in days\n",
    "\n",
    "ztop = 0.0  # surface elevation\n",
    "zbot = -50.0  # bottom elevation of the model\n",
    "\n",
    "z = [1.0, ztop, -10, -15, -25, -26, zbot]\n",
    "kaq = np.array(kh)\n",
    "c = np.array(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reference model with discharge wells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlref = tft.ModelMaq(kaq=kaq, Saq=Ss, c=c, z=z, topboundary=\"semi\", tmin=0.1, tmax=10)\n",
    "w1 = tft.Well(mlref, 0, -10, tsandQ=[(0, Qtot / 2)], rw=0.1, layers=[1, 2])\n",
    "w2 = tft.Well(mlref, 0, 10, tsandQ=[(0, Qtot / 2)], rw=0.1, layers=[1, 2])\n",
    "mlref.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model with WellString"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=kaq, Saq=Ss, c=c, z=z, topboundary=\"semi\", tmin=0.1, tmax=10)\n",
    "ws = tft.WellString(ml, xy=[(0, -10), (0, 10)], tsandQ=[(0, Qtot)], rw=0.1, layers=[1, 2])\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare head contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = 0.2\n",
    "ilay = 1\n",
    "levels = 10\n",
    "ax = mlref.plots.topview(win=(-50, 50, -50, 50))\n",
    "mlref.plots.contour(\n",
    "    win=(-50, 50, -50, 50),\n",
    "    ngr=51,\n",
    "    t=time,\n",
    "    levels=levels,\n",
    "    decimals=2,\n",
    "    layers=[ilay],\n",
    "    ax=ax,\n",
    ")\n",
    "ml.plots.contour(\n",
    "    win=(-50, 50, -50, 50),\n",
    "    ngr=51,\n",
    "    t=time,\n",
    "    levels=levels,\n",
    "    decimals=2,\n",
    "    layers=[ilay],\n",
    "    color=\"C1\",\n",
    "    ax=ax,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare drawdowns along $y$ in layer 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linspace(-50, 50, 101)\n",
    "x = np.zeros_like(y)\n",
    "href = mlref.headalongline(x, y, t=time)\n",
    "h = ml.headalongline(x, y, t=time)\n",
    "\n",
    "ilay = 1\n",
    "plt.figure(figsize=(10, 2))\n",
    "plt.plot(y, href[ilay, 0], label=\"Reference Well\")\n",
    "plt.plot(y, h[ilay, 0], \"k.\", ms=3, label=\"WellString model\")\n",
    "plt.xlabel(\"y [m]\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.legend(loc=(0, 1), frameon=False, ncol=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare discharges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 discharge wells, total Q specified\n",
    "w1.discharge(time), w2.discharge(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WellString, transposed to compare to above\n",
    "ws.discharge_per_well(time).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 discharge wells\n",
    "w1.headinside(time), w2.headinside(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws.headinside(time), ws.wlist[0].headinside(time), ws.wlist[1].headinside(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-layer with skin effect "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlref = tft.ModelMaq(kaq=kaq, Saq=Ss, c=c, z=z, topboundary=\"semi\", tmin=0.1, tmax=10)\n",
    "w1 = tft.Well(mlref, 0, -10, tsandQ=[(0, Qtot / 2)], rw=0.1, res=0.1, layers=[1, 2])\n",
    "w2 = tft.Well(mlref, 0, 10, tsandQ=[(0, Qtot / 2)], rw=0.1, res=0.1, layers=[1, 2])\n",
    "mlref.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=kaq, Saq=Ss, c=c, z=z, topboundary=\"semi\", tmin=0.1, tmax=10)\n",
    "ws = tft.WellString(\n",
    "    ml, xy=[(0, -10), (0, 10)], tsandQ=[(0, Qtot)], rw=0.1, res=0.1, layers=[1, 2]\n",
    ")\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare discharges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 discharge wells, total Q specified\n",
    "w1.discharge(time), w2.discharge(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WellString, transposed to compare to above\n",
    "ws.discharge_per_well(time).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 discharge wells\n",
    "w1.headinside(time), w2.headinside(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws.headinside(time), ws.wlist[0].headinside(time), ws.wlist[1].headinside(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Well string with 3 wells\n",
    "In general, the discharge of each individual well in a well string will differ depending on the setting. Here is a simple example of a string with three wells. The solution will show that the discharge of the cell in the middle is less than the other two cells, but the total discharge of all three wells is the specified `Qtot` and the heads in the three wells are equal at all times. The aquifer is changed to confined. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=kaq, Saq=Ss, c=[800, 400], z=z[1:], topboundary=\"conf\", tmin=0.1, tmax=10\n",
    ")\n",
    "ws = tft.WellString(\n",
    "    ml, xy=[(0, -10), (0, 0), (0, 10)], tsandQ=[(0, Qtot)], rw=0.1, res=0.1, layers=[1, 2]\n",
    ")\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discharges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = 0.2\n",
    "# WellString, transposed to compare to above\n",
    "ws.discharge_per_well(time).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Heads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws.headinside(time), ws.wlist[0].headinside(time), ws.wlist[1].headinside(time)"
   ]
  }
 ],
 "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
}