{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Well near a straight river"
   ]
  },
  {
   "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 well in the middle aquifer of a three aquifer system located at $(x,y)=(0,0)$. The well starts pumping at time $t=0$ at a discharge of $Q=1000$ m$^3$/d. Aquifer properties are the shown in table 3 (same as exercise 2). A stream runs North-South along the line $x=50$. The head along the stream is fixed. \n",
    "\n",
    "**Table 3 - Aquifer properties for exercise 3.**\n",
    "|      Layer    | $k$ (m/d) | $c$ (d) |  $S$  | $S_s$ | $z_t$ (m) | $z_b$ (m)|\n",
    "|---------------| ---------:| -------:| -----:| -----:| ---------:| --------:|\n",
    "|Aquifer 0      |      1    |         |   0.1 |       |   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",
    "\n",
    "## Exercise 3a\n",
    "Model a 1000 m long section of the stream using 12 linesinks with $y$-endpoints at [-500,-300,-200,-100,-50,0,50,100,200,300,500]. Create a cross-section of the head along $y=0$ from $x=-200$ to $x=200$ in all 3 layers."
   ]
  },
  {
   "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=[0.1, 1e-4, 1e-4],\n",
    "    Sll=[0, 0],\n",
    "    phreatictop=True,\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, label=\"well 1\")\n",
    "yls = [-500, -300, -200, -100, -50, 0, 50, 100, 200, 300, 500]\n",
    "xls = 50 * np.ones(len(yls))\n",
    "ls1 = tft.RiverString(\n",
    "    ml, list(zip(xls, yls, strict=False)), tsandh=\"fixed\", layers=0, label=\"river\"\n",
    ")\n",
    "ml.solve()\n",
    "ml.plots.head_along_line(\n",
    "    x1=-200, x2=200, npoints=100, t=100, layers=[0, 1, 2], sstart=-200, figsize=(8, 3)\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3b\n",
    "Compute the discharge of the stream section (the stream depletion) as a function of time from $t=0.1$ till $t=1000$ days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.logspace(-1, 3, 100)\n",
    "Q = ls1.discharge(t)\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.semilogx(t, Q[0])\n",
    "plt.ylabel(\"Q [m$^3$/d]\")\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3c\n",
    "Make a contour plot of each layer after 100 days of pumping. Use 20 grid points in each direction (this may take a little time)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16, 4))\n",
    "for ilay in [0, 1, 2]:\n",
    "    ax = plt.subplot(1, 3, ilay + 1)\n",
    "    ml.plots.contour(\n",
    "        win=[-200, 200, -200, 200],\n",
    "        ngr=[20, 20],\n",
    "        t=100,\n",
    "        layers=ilay,\n",
    "        levels=10,\n",
    "        color=\"C0\",\n",
    "        labels=True,\n",
    "        decimals=2,\n",
    "        ax=ax,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3d\n",
    "The discharge of the well is $Q=1000$ m$^3$/d for 100 days every summer. Compute the stream depletion for a five year period."
   ]
  },
  {
   "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=[0.1, 1e-4, 1e-4],\n",
    "    Sll=[0, 0],\n",
    "    phreatictop=True,\n",
    "    tmin=0.1,\n",
    "    tmax=2000,\n",
    ")\n",
    "tsandQ = [\n",
    "    (0, 1000),\n",
    "    (100, 0),\n",
    "    (365, 1000),\n",
    "    (465, 0),\n",
    "    (730, 1000),\n",
    "    (830, 0),\n",
    "    (1095, 1000),\n",
    "    (1195, 0),\n",
    "    (1460, 1000),\n",
    "    (1560, 0),\n",
    "]\n",
    "w = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=tsandQ, layers=1, label=\"well 1\")\n",
    "yls = [-500, -300, -200, -100, -50, 0, 50, 100, 200, 300, 500]\n",
    "xls = 50 * np.ones(len(yls))\n",
    "ls1 = tft.RiverString(\n",
    "    ml, list(zip(xls, yls, strict=False)), tsandh=\"fixed\", layers=0, label=\"river\"\n",
    ")\n",
    "ml.solve()\n",
    "\n",
    "t = np.linspace(0.1, 2000, 200)\n",
    "Q = ls1.discharge(t)\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.plot(t, Q[0])\n",
    "plt.ylabel(\"Q [m$^3$/d]\")\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.grid()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}