{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A system with wells, rivers, and recharge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a system of three aquifers. The aquifer parameters are presented in Table 1. Note that an average thickness is specified for the top unconfined aquifer. A river with three branches cuts through the upper aquifer. The river is modeled with a string of 7 head-specified line-sinks and each branch is modeled with strings of 5 head-specified line-sinks. The heads are specified at the ends of the line-sinks and are shown in Figure 1. \n",
    "\n",
    "Three wells are present. Well 1 is screened in aquifer 0 and has a discharge of 1000 m$^3$/d, well 2 is screened in aquifer 2 and has a discharge of 5000 m$^3$/d, and well 3 is screened in aquifers 1 and 2 and has a total discharge of 5000 m$^3$/d. A constant recharge through the upper boundary of aquifer 0 is simulated by one large circular infiltration area that covers the entire model area; the recharge rate is 0.2 mm/day. A head of 175 m is specified in layer 0 at the upper righthand corner of the model domain. A layout of all analytic elements, except the boundary of the infiltration area, is shown in Figure 1. \n",
    "\n",
    "**Table 1: Aquifer data for Exercise 2**\n",
    "|Layer        | $k$ (m/d) | $z_b$ (m) | $z_t$ | $c$ (days) | $n$ (-) | $n_{ll}$ (-) |\n",
    "|------------:|----------:|----------:|------:|-----------:|--------:|----------:|\n",
    "|Aquifer 0    |   2       |   140     | 165   |      -     |  0.3    |    -      | \n",
    "|Leaky Layer 1|   -       |   120     | 140   |    2000    |   -     |   0.2     |    \n",
    "|Aquifer 1    |   6       |   80      | 120   |      -     |  0.25   |    -      |  \n",
    "|Leaky Layer 2|   -       |   60      | 80    |    20000   |   -     |   0.25    |  \n",
    "|Aquifer 2    |   4       |   0       | 60    |      -     |  0.3    |    -      ||\n",
    "\n",
    "**Figure 1: Layout of elements for Exercise 2. Heads at centers of line-sinks are indicated.**\n",
    "\n",
    "![Layout of elements for Exercise 2. Heads at centers of line-sinks are indicated.](figs/timml_notebook2_layout.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import timflow.steady as tfs\n",
    "\n",
    "figsize = (6, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create basic model elements\n",
    "ml = tfs.ModelMaq(kaq=[2, 6, 4], z=[165, 140, 120, 80, 60, 0], c=[2000, 20000], npor=0.3)\n",
    "rf = tfs.Constant(ml, xr=20000, yr=20000, hr=175, layer=0)\n",
    "p = tfs.CircAreaSink(ml, xc=10000, yc=10000, R=15000, N=0.0002, layer=0)\n",
    "w1 = tfs.Well(ml, xw=10000, yw=8000, Qw=1000, rw=0.3, layers=0, label=\"well 1\")\n",
    "w2 = tfs.Well(ml, xw=12100, yw=10700, Qw=5000, rw=0.3, layers=2, label=\"well 2\")\n",
    "w3 = tfs.Well(ml, xw=10000, yw=4600, Qw=5000, rw=0.3, layers=[1, 2], label=\"maq well\")\n",
    "#\n",
    "xy1 = [\n",
    "    (833, 14261),\n",
    "    (3229, 14843),\n",
    "    (6094, 15885),\n",
    "    (8385, 15677),\n",
    "    (10781, 14895),\n",
    "    (12753, 14976),\n",
    "]\n",
    "hls1 = [176, 166]\n",
    "xy2 = [\n",
    "    (356, 6976),\n",
    "    (4043, 7153),\n",
    "    (6176, 8400),\n",
    "    (9286, 9820),\n",
    "    (12266, 9686),\n",
    "    (15066, 9466),\n",
    "]\n",
    "hls2 = [174, 162]\n",
    "xy3 = [\n",
    "    (1376, 1910),\n",
    "    (4176, 2043),\n",
    "    (6800, 1553),\n",
    "    (9953, 2086),\n",
    "    (14043, 2043),\n",
    "    (17600, 976),\n",
    "]\n",
    "hls3 = [170, 156]\n",
    "xy4 = [\n",
    "    (9510, 19466),\n",
    "    (12620, 17376),\n",
    "    (12753, 14976),\n",
    "    (13020, 12176),\n",
    "    (15066, 9466),\n",
    "    (16443, 7910),\n",
    "    (17510, 5286),\n",
    "    (17600, 976),\n",
    "]\n",
    "hls4 = [170, np.nan, 166, np.nan, 162, np.nan, np.nan, 156]\n",
    "\n",
    "ls1 = tfs.RiverString(ml, xy=xy1, hls=hls1, layers=0)\n",
    "ls2 = tfs.RiverString(ml, xy=xy2, hls=hls2, layers=0)\n",
    "ls3 = tfs.RiverString(ml, xy=xy3, hls=hls3, layers=0)\n",
    "ls4 = tfs.RiverString(ml, xy=xy4, hls=hls4, layers=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions:\n",
    "### Exercise 2a\n",
    "Solve the model and create a contour plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.solve()\n",
    "ml.plots.contour(\n",
    "    win=[0, 20000, 0, 20000],\n",
    "    ngr=50,\n",
    "    layers=[0, 1, 2],\n",
    "    levels=10,\n",
    "    color=[\"C0\", \"C1\", \"C2\"],\n",
    "    legend=True,\n",
    "    figsize=figsize,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2b\n",
    "What are the heads at the three wells?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The head at well 1 is:\", w1.headinside())\n",
    "print(\"The head at well 2 is:\", w2.headinside())\n",
    "print(\"The head at well 3 is:\", w3.headinside())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2c\n",
    "Create a contour plot including a cross-section.\n",
    "Create 50-year capture zones for all three wells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axes = ml.plots.topview_and_xsection(win=[0, 20000, 0, 20000], figsize=figsize)\n",
    "ml.plots.plotcapzone(\n",
    "    well=w1,\n",
    "    hstepmax=50,\n",
    "    nt=10,\n",
    "    zstart=150,\n",
    "    tmax=250 * 365.25,\n",
    "    orientation=\"both\",\n",
    "    ax=axes,\n",
    ")\n",
    "ml.plots.plotcapzone(\n",
    "    well=w2, hstepmax=50, nt=10, zstart=30, tmax=250 * 365.25, orientation=\"both\", ax=axes\n",
    ")\n",
    "ml.plots.plotcapzone(\n",
    "    well=w3, hstepmax=50, nt=10, zstart=30, tmax=250 * 365.25, orientation=\"both\", ax=axes\n",
    ")\n",
    "ml.plots.plotcapzone(\n",
    "    well=w3,\n",
    "    hstepmax=50,\n",
    "    nt=10,\n",
    "    zstart=100,\n",
    "    tmax=250 * 365.25,\n",
    "    orientation=\"both\",\n",
    "    ax=axes,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axes = ml.plots.topview_and_xsection(\n",
    "    win=[0, 20000, 0, 20000], topfigfrac=0.7, figsize=figsize\n",
    ")\n",
    "ml.plots.plotcapzone(\n",
    "    well=w1, hstepmax=50, nt=10, zstart=150, tmax=50 * 365.25, orientation=\"both\", ax=axes\n",
    ")\n",
    "ml.plots.plotcapzone(\n",
    "    well=w2, hstepmax=50, nt=10, zstart=30, tmax=50 * 365.25, orientation=\"both\", ax=axes\n",
    ")\n",
    "ml.plots.plotcapzone(\n",
    "    well=w3, hstepmax=50, nt=10, zstart=30, tmax=50 * 365.25, orientation=\"both\", ax=axes\n",
    ")\n",
    "ml.plots.plotcapzone(\n",
    "    well=w3, hstepmax=50, nt=10, zstart=100, tmax=50 * 365.25, orientation=\"both\", ax=axes\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}