{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Horizontal wells"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A horizontal well is located in a 20 m thick aquifer; the hydraulic conductivity is $k = 10$ m/d and the vertical\n",
    "anisotropy factor is 0.1. The horizontal well is placed 5 m above the bottom of the\n",
    "aquifer. The well has a discharge of 10,000 m$^3$/d and radius of $r=0.2$ m. The well\n",
    "is 200 m long and runs from $(x, y) = (−100, 0)$ to $(x, y) = (100, 0)$. A long\n",
    "straight river with a head of 40 m runs to the right of the horizontal well along the\n",
    "line $x = 200$. The head is fixed to 42 m at $(x, y) = (−1000, 0)$.\n",
    "\n",
    "Three-dimensional flow to the horizontal well is modeled by dividing the aquifer up in\n",
    "11 layers; the elevations are: `[20, 15, 10, 8, 6, 5.5, 5.2, 4.8, 4.4, 4, 2, 0]`. At\n",
    "the depth of the well, the layer thickness is equal to the diameter of the well, and it\n",
    "increases in the layers above and below the well. A steady timflow model is created with the\n",
    "Model3D command. The horizontal well is located in layer 6 and is modeled with the\n",
    "`Ditch` element. Initially, the entry resistance of the well is set to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "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": [
    "z = [20, 15, 10, 8, 6, 5.5, 5.2, 4.8, 4.4, 4, 2, 0]\n",
    "ml = tfs.Model3D(kaq=10, z=z, kzoverkh=0.1)\n",
    "ls1 = tfs.Ditch(ml, x1=-100, y1=0, x2=100, y2=0, Qls=10000, order=5, layers=6)\n",
    "ls2 = tfs.RiverString(\n",
    "    ml,\n",
    "    [(200, -1000), (200, -200), (200, 0), (200, 200), (200, 1000)],\n",
    "    hls=40,\n",
    "    order=5,\n",
    "    layers=0,\n",
    ")\n",
    "rf = tfs.Constant(ml, xr=-1000, yr=0, hr=42, layer=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions:\n",
    "### Exercise 4a\n",
    "Solve the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4b\n",
    "Create contour plots of layers 0 and 6 and note the difference between the layers. Also,\n",
    "compute the head at $(x, y) = (0, 0.2)$ (on the edge of the well) and notice that there is a very large head\n",
    "difference between the top of the aquifer and the well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(\n",
    "    win=[-150, 150, -150, 150], ngr=[50, 100], layers=[0, 6], figsize=figsize\n",
    ")\n",
    "print(\"The head at the top and in layer 6 are:\")\n",
    "print(ml.head(0, 0.2, [0, 6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4c\n",
    "\n",
    "Draw a number of pathlines from different elevations using the `tracelines` command. First make a plot with a cross section below it.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axes = ml.plots.topview_and_xsection(win=[-1000, 1000, -1000, 1000], figsize=figsize)\n",
    "ml.plots.tracelines(\n",
    "    xstart=[-500, -500, -500],\n",
    "    ystart=[-500, -500, -500],\n",
    "    zstart=[5, 9, 15],\n",
    "    hstepmax=20,\n",
    "    tmax=10 * 365.25,\n",
    "    orientation=\"both\",\n",
    "    color=\"C0\",\n",
    "    ax=axes,\n",
    ")\n",
    "ml.plots.tracelines(\n",
    "    xstart=[250, 250, 250],\n",
    "    ystart=[50, 50, 50],\n",
    "    zstart=[5, 9, 15],\n",
    "    hstepmax=20,\n",
    "    tmax=10 * 365.25,\n",
    "    orientation=\"both\",\n",
    "    color=\"C1\",\n",
    "    ax=axes,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4d\n",
    "Make a contour plot of the heads in a vertical cross-section using the `vcontour` command. Use a cross-section along the well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.vcontour(win=[-200, 300, 0, 0], n=50, levels=20, figsize=figsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4e\n",
    "Change the entry resistance of the horizontal well to 0.01 days and set the width to 0.4 m and resolve the model. Notice the difference in the head inside the horizontal well with the `headinside` function of the horizontal well. Use a "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"head inside w/o resistance:\")\n",
    "print(ls1.headinside())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.Model3D(kaq=10, z=z, kzoverkh=0.1)\n",
    "ls = tfs.Ditch(\n",
    "    ml, x1=-100, y1=0, x2=100, y2=0, Qls=10000, order=5, layers=6, wh=0.4, res=0.01\n",
    ")\n",
    "tfs.RiverString(\n",
    "    ml,\n",
    "    [(200, -1000), (200, -200), (200, 0), (200, 200), (200, 1000)],\n",
    "    hls=40,\n",
    "    order=5,\n",
    "    layers=0,\n",
    ")\n",
    "rf = tfs.Constant(ml, xr=-1000, yr=0, hr=42, layer=0)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"head inside horizontal well:\", ls.headinside())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.vcontour(win=[-200, 300, 0, 0], n=50, levels=20, vinterp=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slanted well\n",
    "A slanted well may be simulated by putting different sections in different layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = [20, 15, 10, 8, 6, 5.5, 5.2, 4.8, 4.4, 4, 2, 0]\n",
    "ml = tfs.Model3D(kaq=10, z=np.linspace(20, 0, 21), kzoverkh=0.1)\n",
    "rf = tfs.Constant(ml, 0, 1000, 20)\n",
    "x = np.linspace(-200, 200, 21)\n",
    "y = np.zeros(21)\n",
    "ls = tfs.RiverString(\n",
    "    ml, xy=list(zip(x, y, strict=False)), hls=10, layers=np.arange(0, 20, 1), order=0\n",
    ")\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.vcontour(win=[-400, 400, 0, 0], n=200, levels=20, vinterp=False, figsize=(16, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the head is not exactly constant along the line-sink segments when using order 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = ml.headalongline(np.linspace(-201, -170, 100), 0)\n",
    "plt.plot(np.linspace(-201, -170, 100), h[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Qtot = np.sum(ls.discharge())\n",
    "print(\"Discharge of slanted well when modeled with fixed head hls=10:\", Qtot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = [20, 15, 10, 8, 6, 5.5, 5.2, 4.8, 4.4, 4, 2, 0]\n",
    "ml = tfs.Model3D(kaq=10, z=np.linspace(20, 0, 21), kzoverkh=0.1)\n",
    "rf = tfs.Constant(ml, 0, 1000, 20)\n",
    "x = np.linspace(-200, 200, 21)\n",
    "y = np.zeros(21)\n",
    "ls = tfs.DitchString(\n",
    "    ml, xy=list(zip(x, y, strict=False)), Qls=Qtot, layers=np.arange(0, 20, 1), order=0\n",
    ")\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.vcontour(win=[-400, 400, 0, 0], n=200, levels=20, vinterp=False, figsize=(16, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Head in slanted well when modeled with fixed discharge:\")\n",
    "[print(lspart.headinside()) for lspart in ls.lslist];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lake with a horizontal well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timflow.steady as tfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[1, 2], z=[10, 5, 4, 0], c=20)\n",
    "xy = [(-5, 0), (0, 0), (5, 0), (5, 8), (-5, 8)]\n",
    "p1 = tfs.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy,\n",
    "    kaq=[0.2, 8],\n",
    "    z=[11, 10, 5, 4, 0],\n",
    "    c=[2, 20],\n",
    "    topboundary=\"semi\",\n",
    "    hstar=1.0,\n",
    "    order=3,\n",
    "    ndeg=1,\n",
    ")\n",
    "w = tfs.Well(ml, xw=0, yw=-10, Qw=100, layers=1)\n",
    "rf = tfs.Constant(ml, xr=0, yr=-100, hr=2)\n",
    "ml.solve()\n",
    "\n",
    "ml.plots.contour([-10, 10, -15, 10], 50, [0, 1], 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[1, 2], z=[10, 5, 4, 0], c=20)\n",
    "xy = [(-5, 0), (0, 0), (5, 0), (5, 8), (-5, 8)]\n",
    "p1 = tfs.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy,\n",
    "    kaq=[0.2, 8],\n",
    "    z=[11, 10, 5, 4, 0],\n",
    "    c=[2, 20],\n",
    "    topboundary=\"semi\",\n",
    "    hstar=1.0,\n",
    "    order=5,\n",
    "    ndeg=2,\n",
    ")\n",
    "ls1 = tfs.DitchString(ml, [(0, -4), (0, 0), (0, 4)], 100, order=3, layers=[1])\n",
    "rf = tfs.Constant(ml, xr=0, yr=-100, hr=2)\n",
    "ml.solve()\n",
    "\n",
    "ml.plots.contour([-10, 10, -5, 10], 50, [0], 30)\n",
    "ml.plots.contour([-10, 10, -5, 10], 50, [1], 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ls in ls1.lslist:\n",
    "    for i in range(ls.ncp):\n",
    "        print(ml.head(ls.xc[i], ls.yc[i], ls.layers))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[1, 2], z=[10, 5, 4, 0], c=20)\n",
    "xy = [(-5, 0), (0, 0), (5, 0), (5, 8), (-5, 8)]\n",
    "p1 = tfs.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy,\n",
    "    kaq=[1, 2],\n",
    "    z=[11, 10, 5, 4, 0],\n",
    "    c=[2, 2],\n",
    "    topboundary=\"semi\",\n",
    "    hstar=1.0,\n",
    "    order=5,\n",
    "    ndeg=2,\n",
    ")\n",
    "ls1 = tfs.DitchString(ml, [(0, -4), (0, 0), (0, 4)], 100, order=3, layers=[1])\n",
    "rf = tfs.Constant(ml, xr=0, yr=-100, hr=2)\n",
    "ml.solve()\n",
    "\n",
    "ml.plots.contour([-10, 10, -5, 10], 51, [0], 30)\n",
    "ml.plots.contour([-10, 10, -5, 10], 51, [1], 30)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}