{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pathline tracing"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pathlines are computed through numerical integration of the velocity vector. Pathlines are computed with the `timtrace` function. The `timtrace` function takes as input arguments the starting locations of the pathline (`xstart`, `ystart`, `zstart`), a list with the starting time and end time of the pathline (`tstartend`), a time offset (`tstartoffset`), and a maximum time step (`tstep`). The pathline starts at the starting time plus the time offset (`tstartoffset`). When the starting time is at the time of a change in boundary condition (e.g., a well starts pumping), the time offset can not be smaller than `tmin`. The `timtrace` function has several keyword arguments to affect the numerical integration procedure. The `hstepmax` keyword is used here, which is the maximum horizontal step (in the length units consistently used throughout the model). \n",
    "\n",
    "The `timtrace` function returns a dictionary with three entries:\n",
    "\n",
    "- `'xyzt'` : 2D array with four columns: x, y, z, t along pathline\n",
    "- `'message'` : list with termination text messages of each section of the pathline\n",
    "- `'status'` : numerical indication of the result. Negative is likely undesirable.\n",
    "    \n",
    "The `'status'` can be one of the following:\n",
    "- -2 : reached maximum number of steps before reaching maximum time\n",
    "- -1 : starting $z$ value not inside aquifer\n",
    "- +1 : reached maximum time\n",
    "- +2 : reached element   \n",
    "\n",
    "The `message` is automatically printed to the screen unless the keyword `silent` is set to `True`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1, a pumping well in a confined aquifer\n",
    "Consider a well that starts pumping at $t=0$ in a confined aquifer. The well is located at the origin of the coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "Q = 100  # discharge of well, m^3/d\n",
    "k = 10  # hydraulic conductivity, m/d\n",
    "H = 10  # thickness of aquifer, m\n",
    "Ss = 1e-4  # specific storage, m^(-1)\n",
    "npor = 0.3  # porosity, -\n",
    "xw = 0  # x-location of well\n",
    "yw = 0  # y-location of well\n",
    "rw = 0.3  # radius of well, m\n",
    "tmin = 0.001  # first time of simulation after change in bc, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=[k], z=[H, 0], Saq=[Ss], poraq=npor, tmin=tmin, tmax=1000, M=10)\n",
    "w = tft.Well(ml, xw=0, yw=0, tsandQ=[(0, Q)], rw=rw)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A pathline is started at $(x,y,z)=(10, 10, 0.5H)$ and time $t=0$ and is followed for 10 days. \n",
    "The output of the `timtrace` function is stored in the `trace` dictionary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = tft.timtrace(\n",
    "    ml,\n",
    "    xstart=10,\n",
    "    ystart=10,\n",
    "    zstart=0.5 * H,\n",
    "    tstartend=[0, 10],\n",
    "    tstartoffset=tmin,\n",
    "    tstep=1,\n",
    "    hstepmax=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entries of the `trace` dictionary are printed below. Note that the pathline starts at $t=0.001$ d, which is the specified `tstartoffset`, that steps are taken with a length of 1 day, and that the pathline is terminated when the maximum time of 10 days is reached. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"xyzt array of pathline:\")\n",
    "print(trace[\"xyzt\"])\n",
    "print(\"trace message:\", trace[\"message\"])\n",
    "print(\"trace status:\", trace[\"status\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the pathline is computed for a max of 20 days (i.e., the last entry in `tstartend` equals 20), the well is reached before $t=20$ d is reached"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = tft.timtrace(\n",
    "    ml,\n",
    "    xstart=10,\n",
    "    ystart=10,\n",
    "    zstart=0.5 * H,\n",
    "    tstartend=[0, 20],\n",
    "    tstartoffset=tmin,\n",
    "    tstep=1,\n",
    "    hstepmax=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"last two entries in xyzt:\")\n",
    "print(trace[\"xyzt\"][-2:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the horizontal stepsize is larger than the specified `hstepmax` value, the timestep is reduced. For example, specifying `hstepmax=0.5` and following the pathline for 4 days gives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = tft.timtrace(\n",
    "    ml,\n",
    "    xstart=10,\n",
    "    ystart=10,\n",
    "    zstart=0.5 * H,\n",
    "    tstartend=[0, 1],\n",
    "    tstartoffset=tmin,\n",
    "    tstep=1,\n",
    "    hstepmax=0.1,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyzt = trace[\"xyzt\"]\n",
    "print(\"xyzt:\")\n",
    "print(trace[\"xyzt\"])\n",
    "x0, y0, z0, t0 = xyzt[1]\n",
    "x1, y1, z1, t1 = xyzt[2]\n",
    "print(f\"length of first step: {np.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2):.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other keyword arguments of `timtrace` are `nstepmax=100`, `silent=False`, and `correctionstep=True`. The `nstepmax` is the maximum number of steps. Numerical integration stops when the maximum number of steps is reached and the returned `status` is `-2`. For most practical cases, this is an undesirable result and the maximum number of steps should be increased. Setting `silent=True` prevents the message from being printed to the screen, which may be useful when a lot of pathlines are computed, for example in a loop. When the `correctionstep` is set to `True` (default), the numerical integration scheme is a predictor-correcter scheme. When `correctionstep` is set to `False`, the integration scheme is forward integration through time, which is less accurate but quicker. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2, injection and recovery of a pumping well in a confined aquifer\n",
    "Consider a well in a confined aquifer. The well starts pumping with discharge $Q$ at time $t=0$ and starts injection with discharge $Q$ at time $t=10$ d. The well is located at the origin of the coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=[k], z=[H, 0], Saq=[Ss], poraq=npor, tmin=tmin, tmax=1000, M=10)\n",
    "w = tft.Well(ml, xw=0, yw=0, tsandQ=[(0, Q), (10, -Q)], rw=rw)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A pathline is started at $(x,y,z)=(10, 10, 5)$ and followed for 20 days. The `tstartend` list is now `[0, 10, 20]`. This means that the pathline consists of two sections. The first section starts at $t=0$ d (or really at `tmin`) and is followed up till $t=10$ d. The second section starts from the endpoint at $t=10$ d of the first section at $t=10$ d (or really at 10+`tmin` days) and is followed up till $t=20$ d. The `nstepmax` keyword is the maximum number of steps for each section of the pathline. There are two messages printed to the screen, one for each section. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = tft.timtrace(\n",
    "    ml,\n",
    "    xstart=10,\n",
    "    ystart=10,\n",
    "    zstart=0.5 * H,\n",
    "    tstartend=[0, 10, 20],\n",
    "    tstartoffset=tmin,\n",
    "    tstep=1,\n",
    "    hstepmax=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $x,y,z,t$ array is printed below. Note that the first section ends at $t=10$ d. \n",
    "After 20 d, the pathline ends at (almost) the same location as it started. The small discrepancy is caused by numerical error of the integration scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyzt = trace[\"xyzt\"]\n",
    "print(\"xyzt:\")\n",
    "print(trace[\"xyzt\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3, a well in a semi-confined aquifer\n",
    "Consider an injection well in a semi-confined aquifer. The well starts injecting with discharge $Q$ at time $t=0$. The well is located at the origin of the coordinate system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "k = 20  # hydraulic conductivity aquifer, m/d\n",
    "H = 10  # thickness of aquifers, m\n",
    "Hstar = 2  # thickness of leaky layer, m\n",
    "c = 100  # resistance of leaky layer, d\n",
    "Ss = 1e-4  # specific storage of both aquifers, m^(-1)\n",
    "npor = 0.3  # porosity of both aquifers, -\n",
    "Q = 1000  # discharge of well in aquifer 1, m^3/d\n",
    "xw = 0  # x-location of well\n",
    "yw = 0  # y-location of well\n",
    "rw = 0.3  # radius of well, m\n",
    "tmin = 0.001  # first time of simulation after change in bc, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=[k],\n",
    "    z=[H + Hstar, H, 0],\n",
    "    c=[c],\n",
    "    Saq=Ss,\n",
    "    poraq=npor,\n",
    "    topboundary=\"semi\",\n",
    "    tmin=tmin,\n",
    "    tmax=1000,\n",
    "    M=10,\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, tsandQ=[(0, -Q)], layers=0, rw=0.3)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nine pathlines are started at the well screen at different elevations in the aquifer and followed for a max of 1000 days. Five of the nine pathlines end at the top of the semi-confining layer within the 1000 days (the five pathlines that start at the highest elevations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(111, aspect=10)\n",
    "zstart = np.arange(1.0, 10, 1)\n",
    "for zs in zstart:\n",
    "    trace = tft.timtrace(\n",
    "        ml,\n",
    "        xstart=rw,\n",
    "        ystart=0,\n",
    "        zstart=zs,\n",
    "        tstartend=[0, 1000],\n",
    "        tstartoffset=tmin,\n",
    "        tstep=100,\n",
    "        nstepmax=100,\n",
    "        hstepmax=2,\n",
    "        silent=True,\n",
    "    )\n",
    "    xyzt = trace[\"xyzt\"]\n",
    "    plt.plot(xyzt[:, 0], xyzt[:, 2])\n",
    "for y in [0, H, H + Hstar]:\n",
    "    plt.axhline(y, color=\"k\")\n",
    "plt.xlabel(\"$x$ (m)\")\n",
    "plt.ylabel(\"$z$ (m)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 4. A well in the bottom aquifer of a two-aquifer system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a pumping well in the bottom aquifer of a two-aquifer system; the aquifers are separated by a leaky layer. The well starts pumping with discharge $Q$ at time $t=0$. The well is located at the origin of the coordinate system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "k0 = 20  # hydraulic conductivity aquifer 0, m/d\n",
    "k1 = 40  # hydraulic conductivity aquifer 1, m/d\n",
    "H = 10  # thickness of both aquifers, m\n",
    "Hstar = 2  # thickness of leaky layer, m\n",
    "c = 100  # resistance of leaky layer, d\n",
    "Ss = 1e-4  # specific storage of both aquifers, m^(-1)\n",
    "npor = 0.3  # porosity of both aquifers, -\n",
    "Q = 100  # discharge of well in aquifer 1, m^3/d\n",
    "xw = 0  # x-location of well\n",
    "yw = 0  # y-location of well\n",
    "rw = 0.3  # radius of well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=[k0, k1],\n",
    "    z=[H + Hstar + H, H + Hstar, H, 0],\n",
    "    c=[c],\n",
    "    Saq=[Ss, Ss],\n",
    "    poraq=npor,\n",
    "    porll=npor,\n",
    "    tmin=0.01,\n",
    "    tmax=10000,\n",
    "    M=10,\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, tsandQ=[(0, Q)], layers=1, rw=0.3)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Five pathlines are started in the top aquifer an five in the bottom aquifer at a distance of 200 m from the well and followed for a maximum of 10,000 days. The pathlines in the bottom aquifer all reach the well within that time, but the pathlines in the top aquifer don't. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zstart = np.hstack((np.arange(1, 10, 2), np.arange(13, 22, 2)))\n",
    "plt.subplot(111, aspect=5)\n",
    "for zs in zstart:\n",
    "    trace = tft.timtrace(\n",
    "        ml,\n",
    "        xstart=200,\n",
    "        ystart=0,\n",
    "        zstart=zs,\n",
    "        tstartend=[0, 10000],\n",
    "        tstartoffset=0.01,\n",
    "        tstep=100,\n",
    "        nstepmax=100,\n",
    "        hstepmax=10,\n",
    "        silent=True,\n",
    "    )\n",
    "    xyzt = trace[\"xyzt\"]\n",
    "    plt.plot(xyzt[:, 0], xyzt[:, 2])\n",
    "for y in [0, H, H + Hstar, H + Hstar + H]:\n",
    "    plt.axhline(y, color=\"k\")\n",
    "plt.xlabel(\"$x$ (m)\")\n",
    "plt.ylabel(\"$z$ (m) - VE:5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5. A partially penetrating well in a four-layer system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a partially penetrating injection well. The aquifer is simulated by 4 layers and the well is screened in the second layer. The well starts pumping with discharge $Q$ at time $t=0$. The well is located at the origin of the coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 10  # hydraulic conductivity of aquifer, m/d\n",
    "H = 5  # thickness of each layer, m\n",
    "Q = 100  # injection rate of well, m^3/d\n",
    "Ss = 1e-4  # specific storage of both aquifers, m^(-1)\n",
    "npor = 0.3  # porosity of both aquifers, -\n",
    "tmin = 0.01  # minimum time, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.Model3D(\n",
    "    kaq=k, z=np.arange(4, -1, -1) * H, Saq=Ss, poraq=npor, tmin=tmin, tmax=1000\n",
    ")\n",
    "w = tft.Well(ml, xw=0, yw=0, tsandQ=[(0, -Q)], layers=1, rw=0.1)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start 11 pathlines near the well and trace for 100 days. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zstart = np.linspace(10.01, 14.99, 11)\n",
    "for zs in zstart:\n",
    "    trace = tft.timtrace(\n",
    "        ml,\n",
    "        xstart=1,\n",
    "        ystart=0,\n",
    "        zstart=zs,\n",
    "        tstartend=[0, 100],\n",
    "        tstartoffset=tmin,\n",
    "        tstep=5,\n",
    "        nstepmax=200,\n",
    "        hstepmax=2,\n",
    "        silent=True,\n",
    "    )\n",
    "    xyzt = trace[\"xyzt\"]\n",
    "    plt.plot(xyzt[:, 0], xyzt[:, 2])\n",
    "for z in np.arange(4, -1, -1) * H:\n",
    "    plt.axhline(z, color=\"k\")\n",
    "plt.axis(\"scaled\")\n",
    "plt.xlabel(\"$x$ (m)\")\n",
    "plt.ylabel(\"$z$ (m)\")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}