{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Well in multi-layer system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.transient as tft\n",
    "\n",
    "plt.rcParams[\"font.size\"] = 8.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a three-aquifer system. Aquifer properties are given in Table 1. All aquifers have elastic storage. A well is located at $(x,y)=(0,0)$ and is screened in layer 1. The well starts pumping at time $t=0$ with a discharge $Q=1000$ m$^3$/d. The radius of the well is 0.2 m.\n",
    "\n",
    "**Table 1 - Aquifer properties for exercise 1**\n",
    "|Layer          | $k$ (m/d) | $c$ (d) | $S_s$ (m$^{-1}$) | $z_t$ (m) | $z_b$ (m)|\n",
    "|---------------| ---------:| -------:| -----:| ---------:| --------:|\n",
    "|Aquifer 0      |      1    |    -    |0.0001 |   25      |        20|\n",
    "|Leaky layer 1  |      -    |  1000   |0      |   20      |        18|\n",
    "|Aquifer 1      |     20    |    -    |0.0001 |   18      |        10|\n",
    "|Leaky layer 2  |      -    |  2000   |0      |   10      |         8|\n",
    "|Aquifer 2      |      2    |    -    |0.0001 |    8      |         0|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1a\n",
    "\n",
    "Compute the head as a function of time at $(x,y)=(50,0)$. Make a plot of the head vs. time from $t=0.1$ till $t=1000$ days using a linear scaling on both axis. Do the same using a logarithmic time axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=[1, 20, 2],\n",
    "    z=[25, 20, 18, 10, 8, 0],\n",
    "    c=[1000, 2000],\n",
    "    Saq=[1e-4, 1e-4, 1e-4],\n",
    "    Sll=[0, 0],\n",
    "    phreatictop=False,\n",
    "    tmin=0.1,\n",
    "    tmax=1000,\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000)], layers=1)\n",
    "ml.solve()\n",
    "\n",
    "t = np.logspace(-1, 3, 100)\n",
    "h = ml.head(50, 0, t)\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.subplot(121)\n",
    "plt.plot(t, h[0], label=\"layer 0\")\n",
    "plt.plot(t, h[1], label=\"layer 1\")\n",
    "plt.plot(t, h[2], label=\"layer 2\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.grid()\n",
    "plt.subplot(122)\n",
    "plt.semilogx(t, h[0], label=\"layer 0\")\n",
    "plt.semilogx(t, h[1], label=\"layer 1\")\n",
    "plt.semilogx(t, h[2], label=\"layer 2\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1b\n",
    "Create a plot of the head vs. distance from the well after 10 days of pumping. Plot the head in all three layers up to a distance of 1000 m from the well. Make the same plot after 1000 days of pumping. Is there much difference between 100 and 1000 days of pumping?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.head_along_line(\n",
    "    x1=-1000, x2=1000, y1=0, y2=0, npoints=100, t=10, layers=[0, 1, 2], figsize=(8, 3)\n",
    ")\n",
    "plt.xlabel(\"distance (m)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "\n",
    "ml.plots.head_along_line(\n",
    "    x1=-1000, x2=1000, y1=0, y2=0, npoints=100, t=1000, layers=[0, 1, 2], figsize=(8, 3)\n",
    ")\n",
    "plt.xlabel(\"distance (m)\")\n",
    "plt.ylabel(\"head (m)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1c\n",
    "The well is turned off after 100 days. Compute the head as a function of time at $(x,y)=(50,0)$. Make a plot of the head vs. time from $t=0.1$ till $t=1000$ days using a logarithmic time axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=[1, 20, 2],\n",
    "    z=[25, 20, 18, 10, 8, 0],\n",
    "    c=[1000, 2000],\n",
    "    Saq=[1e-4, 1e-4, 1e-4],\n",
    "    Sll=[0, 0],\n",
    "    phreatictop=False,\n",
    "    tmin=0.1,\n",
    "    tmax=1000,\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000), (100, 0)], layers=1)\n",
    "ml.solve()\n",
    "t = np.logspace(-1, 3, 100)\n",
    "h = ml.head(50, 0, t)\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.semilogx(t, h[0], label=\"layer 0\")\n",
    "plt.semilogx(t, h[1], label=\"layer 1\")\n",
    "plt.semilogx(t, h[2], label=\"layer 2\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1d\n",
    "Compute the head inside the well as a function of time from $t=0.1$ till $t=1000$ days using a logarithmic time axis. On the same graph, plot the head inside the well vs. time when the entry resistance of the well is 0.1 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "h = w.headinside(t)\n",
    "plt.semilogx(t, h[0], label=\"res=0\")  # head from previous solution\n",
    "w.res = 0.1\n",
    "ml.solve()\n",
    "h = w.headinside(t)\n",
    "plt.semilogx(t, h[0], label=\"res=0.1\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1e\n",
    "Conside again the case of a well without skin effect. \n",
    "Compute the head inside the well as a function of time from $t=1$ min till $t=1$ day using a logarithmic time axis. On the same graph, plot the head inside the well vs. time when the wellbore storage is taken into account."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmin = 1.0 / 24 / 60  # 1 minute\n",
    "ml = tft.ModelMaq(\n",
    "    kaq=[1, 20, 2],\n",
    "    z=[25, 20, 18, 10, 8, 0],\n",
    "    c=[1000, 2000],\n",
    "    Saq=[1e-4, 1e-4, 1e-4],\n",
    "    Sll=[0, 0],\n",
    "    phreatictop=False,\n",
    "    tmin=1e-4,\n",
    "    tmax=1,\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000)], layers=1)\n",
    "ml.solve()\n",
    "t = np.logspace(np.log10(tmin), 0, 100)\n",
    "h = w.headinside(t)\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.semilogx(t, h[0], label=\"no wellbore storage\")  # head from previous solution\n",
    "w.rc = 0.2\n",
    "ml.solve()\n",
    "h = w.headinside(t)\n",
    "plt.semilogx(t, h[0], label=\"wellbore storage\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.xticks([tmin, 10 * tmin, 60 * tmin, 1], [\"1 min\", \"10 min\", \"1 hr\", \"1 day\"])\n",
    "plt.grid()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}