{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 1: Well in a Multi-Aquifer System"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uniform flow to a well in a multi-aquifer system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will simulate steady flow to an extraction well in the middle aquifer of a three-aquifer system. \n",
    "Aquifer properties are given in the table shown below. The well is located at $(x,y)=(0,0)$, the discharge is $Q=10,000$ m$^3$/d and the radius is 0.2 m. There is a uniform flow from West to East with a gradient of 0.002. The head is fixed to 20 m at a distance of 10,000 m downstream of the well. \n",
    "\n",
    "### Aquifer properties\n",
    "|Layer        |$k$ (m/d)|$z_b$ (m)|$z_t$|$c$ (days)|\n",
    "|-------------|--------:|--------:|----:|---------:|\n",
    "|Aquifer 0    |    10   |   -20   |  0  |    -     |\n",
    "|Leaky Layer 1|     -   |   -40   | -20 |  4000    |   \n",
    "|Aquifer 1    |    20   |   -80   | -40 |    -     |\n",
    "|Leaky Layer 2|     -   |   -90   | -80 |  10000   | \n",
    "|Aquifer 2    |     5   |   -140  | -90 |    -     ||\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start the model by importing the regular packages and `timflow.steady`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import timflow.steady as tfs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a multi-aquifer model with the `ModelMaq` class. We specify the hydraulic conductivitity for each aquifer (`kaq`), followed by a vector `z` with the top and bottom of each aquifer from the top down (i.e., 6 values, as there are 3 aquifers). \n",
    "The model instance is stored in the variable `ml`. Next, three elements are added to the model: a well (`Well`), a reference point with a fixed head (`Constant`), and uniform flow (`Uflow`).\n",
    "The well is screened in the middle aquifer layer (which is number 1). Since it is screened in only 1 layer, it suffices to provide one value for the keyword `layers`. The well is stored in the variable `w` for later use. \n",
    "The reference point is added by specifying the location, the head, and the layer. Note that only one reference point can be added to the model and in only one layer (hence, the keyword is `layer` and not `layers`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(\n",
    "    kaq=[10, 20, 5],  # hydraulic conductivity, m/d\n",
    "    z=[0, -20, -40, -80, -90, -140],  # tops and bottoms of aquifers, m\n",
    "    c=[4000, 10000],  # resistance of leaky layers, d\n",
    "    npor=0.3,  # porosity of the aquifers, one value so the same for all aquifers, -\n",
    ")\n",
    "w = tfs.Well(\n",
    "    model=ml,  # model to which element is added\n",
    "    xw=0,  # x-location of well, m\n",
    "    yw=0,  # y-location of well, m\n",
    "    Qw=10_000,  # discharge of well, positive for extraction, m^3/d\n",
    "    rw=0.2,  # well radius, m\n",
    "    layers=1,  # layer numbere where well is screened (may also be a list)\n",
    "    label=\"well 1\",\n",
    ")\n",
    "tfs.Constant(\n",
    "    model=ml,  # model to which element is added\n",
    "    xr=10_000,  # x-location of fixed head, m\n",
    "    yr=0,  # y-location of fixed head, m\n",
    "    hr=20,  # fixed head, m\n",
    "    layer=0,  # layer number where head is fixed\n",
    ")\n",
    "tfs.Uflow(\n",
    "    model=ml,  # model to which element is added\n",
    "    slope=0.002,  # head drop / distance in direction of flow\n",
    "    angle=0,  # direction of uniform flow, straight East is zero degrees\n",
    ")\n",
    "ml.solve()  # solve the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.nunknowns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transmissivity of each aquifer is computed and stored by `timflow`. Each `timflow` model stores all aquifer properties as `aq`. The transmissivities are stored in the variable `T`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.aq.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that the transmissivity of the middle aquifer is 800 m$^2$/d (as is indeed the product of the hydraulic conductivity and thickness of aquifer layer 1, see the values in the table). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A confined aquifer system with three aquifer layers and two leaky layers has two leakage factors. The heads in all three aquifers are approximately equal at a distance of three times the largest leakage factor away from the well. The leakage factors are computed by `timflow` and may be obtained as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"leakage factors of model (m)\")\n",
    "print(ml.aq.lab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `timflow` actually returns three leakage factors (as the system consists of\n",
    "three aquifers), but since the top of the aquifer system is impermeable, the first\n",
    "leakage factor is equal to 0 (it will be a non-zero value for a semi-confined system,\n",
    "where the top aquifer is covered by a leaky layer with a fixed head above it)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One contour plot of the heads in all three aquifers is created. \n",
    "Contours are drawn inside a window with lower-left hand corner $(x,y)=(−3000,−3000)$ and upper right-hand corner $(x,y)=(3000,3000)$. Notice that the heads in the three aquifers are almost equal at three times the largest leakage factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(\n",
    "    win=[-3000, 3000, -3000, 3000],  # window to contour [xmin, xmax, ymin, ymax]\n",
    "    ngr=50,  # number of points where to compute the head\n",
    "    layers=[0, 1, 2],  # layers to contour\n",
    "    levels=10,  # draw 10 contours in each layer\n",
    "    labels=True,  # add labels along the contours\n",
    "    decimals=1,  # print labels with 1 decimal place\n",
    "    legend=True,  # add a legend\n",
    "    figsize=(5, 5),  # specify a figure size of 5 by 5 inches\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The head in the well is computed as before. Note that the drawdown is quite large, as the head at the well is 40 m in absence of the well.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"The head at the well is: {w.headinside()} m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The head at the well may also be computed with the `ml.head` function, which returns the head at the well in all three layers. (The careful reader may note that the head in layer 1 is slightly different than the value obtained with `w.headinside`, because the head is computed slightly different, but the difference is less than 1 mm.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.head(0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, six pathlines are computed and plotted in both a plan view and a cross-section. Three pathlines are started from from $(x,y)=(-2000,-1000)$, but at three different levels $z=-120$, $z=-60$, and $z=-10$ (i.e., one from each aquifer). The other three pathlines are started from the same elevations, but from $(x,y)=(0, 1000)$. A plot of the aquifer must be created before pathlines can be visualized. In the code cell below, the plot is created with the `ml.plot` function, which only plots the elements to the screen (in this case only the well, as there are no elements to visualize). In addition, the `orientation` is set to `both`, which means that both a plan view and a vertical cross-section are plotted.\n",
    "Note that the pathlines are automatically colored: they are blue (`C0`) in aquifer 0, orange (`C1`) in aquifer 1 and green (`C2`) in aquifer 2. Note that the pathline that starts in the top aquifer at $(-2000, -1000)$ leaks to the middle aquifer (when it changes from blue to orange) and ends at the well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "win = [-3000, 3000, -3000, 3000]  # window that is plotted [xmin, xmax, ymin, ymax]\n",
    "axes = ml.plots.topview_and_xsection(\n",
    "    win=win,\n",
    "    figsize=(5, 5),  # specify a figure size of 5 by 5 inches\n",
    ")\n",
    "ml.plots.tracelines(\n",
    "    win=win,\n",
    "    xstart=-2000 * np.ones(3),  # x-locations of starting points, m\n",
    "    ystart=-1000 * np.ones(3),  # y-locations of starting points, m\n",
    "    zstart=[-120, -60, -10],  # z-locations of starting points, m\n",
    "    hstepmax=50,  # maximum horizontal step size, m\n",
    "    orientation=\"both\",  # plot both in plan and cross-sectional view\n",
    "    ax=axes,  # use the axes created before\n",
    ")\n",
    "ml.plots.tracelines(\n",
    "    win=win,\n",
    "    xstart=0 * np.ones(3),\n",
    "    ystart=1000 * np.ones(3),\n",
    "    zstart=[-120, -50, -10],\n",
    "    hstepmax=50,\n",
    "    orientation=\"both\",\n",
    "    ax=axes,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}