examples/Notebooks/flopy3_get_transmissivities_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Demonstration of `flopy.utils.get_transmissivities` method\n",
    "for computing open interval transmissivities (for weighted averages of heads or fluxes)\n",
    "In practice this method might be used to:  \n",
    "\n",
    "* compute vertically-averaged head target values representative of observation wells of varying open intervals (including variability in saturated thickness due to the position of the water table)\n",
    "* apportion boundary fluxes (e.g. from an analytic element model) among model layers based on transmissivity.\n",
    "* any other analysis where a distribution of transmissivity by layer is needed for a specified vertical interval of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.6.4 | packaged by conda-forge | (default, Dec 23 2017, 16:54:01) \n",
      "[GCC 4.2.1 Compatible Apple LLVM 6.1.0 (clang-602.0.53)]\n",
      "numpy version: 1.14.0\n",
      "matplotlib version: 2.1.2\n",
      "flopy version: 3.2.9\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# run installed version of flopy or add local path\n",
    "try:\n",
    "    import flopy\n",
    "except:\n",
    "    fpth = os.path.abspath(os.path.join('..', '..'))\n",
    "    sys.path.append(fpth)\n",
    "    import flopy\n",
    "\n",
    "print(sys.version)\n",
    "print('numpy version: {}'.format(np.__version__))\n",
    "print('matplotlib version: {}'.format(mpl.__version__))\n",
    "print('flopy version: {}'.format(flopy.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make up some open interval tops and bottoms and some heads\n",
    "* (these could be lists of observation well screen tops and bottoms)\n",
    "* the heads array contains the simulated head in each model layer,\n",
    "  at the location of each observation well (for example, what you would get back from HYDMOD if you had an entry for each layer at the location of each head target).\n",
    "* make up a model grid with uniform horizontal k of 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[2., 2., 2., 2., 2., 2.],\n",
       "        [2., 2., 2., 2., 2., 2.],\n",
       "        [2., 2., 2., 2., 2., 2.],\n",
       "        [2., 2., 2., 2., 2., 2.],\n",
       "        [2., 2., 2., 2., 2., 2.],\n",
       "        [2., 2., 2., 2., 2., 2.]],\n",
       "\n",
       "       [[1., 1., 1., 1., 1., 1.],\n",
       "        [1., 1., 1., 1., 1., 1.],\n",
       "        [1., 1., 1., 1., 1., 1.],\n",
       "        [1., 1., 1., 1., 1., 1.],\n",
       "        [1., 1., 1., 1., 1., 1.],\n",
       "        [1., 1., 1., 1., 1., 1.]],\n",
       "\n",
       "       [[0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0.]]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sctop = [-.25, .5, 1.7, 1.5, 3., 2.5] # screen tops\n",
    "scbot = [-1., -.5, 1.2, 0.5, 1.5, -.2] # screen bottoms\n",
    "# head in each layer, for 6 head target locations\n",
    "heads = np.array([[1., 2.0, 2.05, 3., 4., 2.5],\n",
    "                  [1.1, 2.1, 2.2, 2., 3.5, 3.],\n",
    "                  [1.2, 2.3, 2.4, 0.6, 3.4, 3.2]\n",
    "        ])\n",
    "nl, nr = heads.shape\n",
    "nc = nr\n",
    "botm = np.ones((nl, nr, nc), dtype=float)\n",
    "top = np.ones((nr, nc), dtype=float) * 2.1\n",
    "hk = np.ones((nl, nr, nc), dtype=float) * 2.\n",
    "for i in range(nl):\n",
    "    botm[nl-i-1, :, :] = i\n",
    "botm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a flopy modflow model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = flopy.modflow.Modflow('junk', version='mfnwt', model_ws='data')\n",
    "dis = flopy.modflow.ModflowDis(m, nlay=nl, nrow=nr, ncol=nc, botm=botm, top=top)\n",
    "upw = flopy.modflow.ModflowUpw(m, hk=hk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get transmissivities along the diagonal cells\n",
    "* alternatively, if a model `SpatialReference` has been set up, the real-world x and y coordinates could be supplied with the `x` and `y` arguments\n",
    "* if `sctop` and `scbot` arguments are given, the transmissivites are computed for the open intervals only\n",
    "  (cells that are partially within the open interval have reduced thickness, cells outside of the open interval have transmissivities of 0). If no `sctop` or `scbot` arguments are supplied, trasmissivites reflect the full saturated thickness in each column of cells (see plot below, which shows different open intervals relative to the model layering)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0. , 0. , 0. , 0. , 0.2, 0.2],\n",
       "       [0. , 0. , 1. , 1. , 1. , 2. ],\n",
       "       [2. , 1. , 0. , 0.2, 0. , 2. ]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r, c = np.arange(nr), np.arange(nc)\n",
    "T = flopy.utils.get_transmissivities(heads, m, r=r, c=c, sctop=sctop, scbot=scbot)\n",
    "np.round(T, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2., 2., 2., 2., 2., 2.],\n",
       "       [1., 1., 1., 1., 1., 1.],\n",
       "       [0., 0., 0., 0., 0., 0.]], dtype=float32)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.dis.botm.array[:, r, c]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the model top and layer bottoms (colors)\n",
    "open intervals are shown as boxes\n",
    "* well 0 has zero transmissivities for each layer, as it is below the model bottom\n",
    "* well 1 has T values of 0 for layers 1 and 2, and 1 for layer 3 (K=2 x 0.5 thickness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x12022f668>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XtcVVX+//HXQlFTs0Dtm6MilqQiIiA6OHipLPOWUzxyyvCS5mhZjTPllNWk2c9xamoap7Jp/FpqaVppmdPVNM3oa6mkmdeywtRSAdHyLrJ+f2zAS6ggm7MP+7yfj8d5dDZs9vqc0/HDYu21PstYaxEREf8I8zoAERFxlxK7iIjPKLGLiPiMEruIiM8osYuI+IwSu4iIzyixi4j4jBK7iIjPKLGLiPhMVS8arVevno2OjvaiaRGRSiszMzPHWlv/bOd5ktijo6NZuXKlF02LiFRaxpgtpTlPQzEiIj6jxC4i4jNK7CIiPuPJGHtJjh49yrZt2zh06JDXoUgQqVGjBo0aNSI8PNzrUEQqjaBJ7Nu2beP8888nOjoaY4zX4UgQsNaSm5vLtm3baNq0qdfhiFQarg3FGGOqGGNWGWPeOpefP3ToEHXr1lVSl2LGGOrWrau/4kTKyM0x9pHAhvJcQEldTqXPhEjZuZLYjTGNgF7AFDeuJyKV3+HDkJ/vdRShya0e+0TgXqDgdCcYY4YZY1YaY1ZmZ2e71Gzwio6OJicn55zOmTBhQkWFJRIQ+flw7bUwbJjXkYSmcid2Y0xvYJe1NvNM51lrJ1trk621yfXrn3VFbEhTYpfKaP9+6N8fXn0VqlaFTp2gc2evowpNbvTYU4E+xpgsYDZwpTFmhgvXDaisrCxatGjB0KFDiYuLIz09nYULF5KamkpMTAzLly8HYPfu3Vx33XXEx8eTkpLCmjVrAMjNzaVbt24kJiYyfPhwrLXF154xYwbt27cnISGB4cOHc+zYsdPGMXr0aA4ePEhCQgLp6ekAPPnkk8TFxREXF8fEiRNPinfQoEHEx8dzww03cODAgYp6e0TOqlo12LIFtm1zjh96CG65xXm+a5dnYYUkc2ICKvfFjLkcGGWt7X2m85KTk+2ptWI2bNhAy5YtARj333Ws/+En1+ICiP1VHcZe2+q038/KyqJZs2asWrWKVq1a0a5dO9q0acPzzz/P/PnzmTp1KvPmzeOuu+6iXr16jB07lg8//JC7776b1atX84c//IF69eoxZswY3n77bXr37k12djbZ2dnce++9vP7664SHhzNixAhSUlIYOHBgcc2cevXqnRRL7dq12bdvHwCZmZnccsstfPrpp1hr+fWvf82MGTOIiIigadOmZGRkkJqaypAhQ4iNjWXUqFGuvm/B4MTPhgSX/HyYNAl+/3uoWdM5rnrKJOqNGyElBZ54AoYO9SZOvzDGZFprk892nlaenqBp06a0bt2asLAwWrVqRdeuXTHG0Lp1a7KysgDIyMhgwIABAFx55ZXk5uayd+9eli5dSv/+/QHo1asXERERACxatIjMzEzatWtHQkICixYt4ttvvy11TBkZGVx//fXUqlWL2rVrk5aWxscffwxA48aNSU1NBaB///5kZGS49VaIlMrKlfCnP8ErrzjHpyZ1gGbNnJ771VcHNLSQ5uoCJWvtEmBJea9zpp51RapevXrx87CwsOLjsLAw8gtv75f0F07RlLySpuZZaxk0aBB/+9vfzimmM/1FdWp7mhoogbJ5s5OwU1Lg888hIeH051atCoUjiFgLX38Nl10WmDhDlXrsZdS5c2dmzpwJwJIlS6hXrx516tQ56evvvvsueXl5AHTt2pU5c+awq3CQcffu3WzZcubKm+Hh4Rw9erS4vXnz5nHgwAH279/PG2+8QadOnQD4/vvvWbZsGQCzZs2iY8eO7r9gkVNMnAhxcbBpk3N8pqR+qkmTID4evvyyYmITR9CUFKgsHn74YQYPHkx8fDw1a9Zk+vTpAIwdO5Z+/fqRlJREly5diIqKAiA2Npbx48fTrVs3CgoKCA8PZ9KkSTRp0uS0bQwbNoz4+HiSkpKYOXMmt9xyC+3btwdg6NChJCYmkpWVRcuWLZk+fTrDhw8nJiaG22+/veLfAAlZ1oIxcPPNcOQIXHpp2a/Rrx/s2wetvPmjPGS4evO0tM5281TOLisri969e7N27VqvQ6lw+mx4b8IEWLXKmcro1ohfTg7k5kLz5u5cLxTo5qmIuKZ6dahRw+mpu6VfP+jVCwpHHcVFGoqppKKjo0Oity7eeflliIqCjh3h7rvd66kX+ec/nR67KjK7T4ldRH7h0CEYMwbatXMSe0VMuIqLO/580SJITITISPfbCUUaihGRYqtWwbFjzrDLhx/CSy9VfJs5OfDb38Lo0RXfVqhQYhcRwJmCmJwMTz/tHEdFlbzgyG316sF//+usTBV3KLGLhLii0rpxcc4888GDAx/DFVdAnTrHSxSo3G/5KLGfoHbt2gFt75lnnqFZs2YYY05b4nfatGnceeedpb7mnj17ePbZZ90KUXzu3XedVaDbtzvj6LfdBhdc4F0877wDd97pxCXnTok9QKy1FBScXK4+NTWVhQsXnnGxUlkpsUtZXHopxMQ44+rBoE8f+Owzp5a7nDsl9hLs27ePrl27kpSUROvWrXnzzTcBeOihh/jXv/5VfN6DDz7IU089BcDjjz9Ou3btiI+PZ+zYsQDFq0NHjBhBUlISW7duPamdxMREoqOjzxrP1q1b6d69O82bN2fcuHHFXy+pnO/o0aP55ptvSEhI4M9//nO53gfxp//7P2fBETi99fffd8bTg0XhIms2bTo+3i9lE5zTHd8dDTtcLiZxcWvo8WipTq1RowZvvPEGderUIScnh5SUFPr06cOtt95KWloaI0eOpKCggNmzZ7N8+XIWLFjA119/zfLly7HW0qdPH5YuXUpUVBSbNm1i6tSp5epFL1++nLVr11KzZk3atWtHr169MMYwdepUPvvss+Jyvl26dOHRRx9l7dq1rF69+pzbE3+bMwfmzXOGPOrU8Tqa03v2WZg92ylhULeu19FULsGZ2D1mreWBBx5g6dKlhIWFsX37dnbu3El0dDR169Zl1apV7Ny5k8TEROrWrcuCBQtYsGABiYmJgNPj//rrr4mKiqJJkyakpKSUK56rr76auoWf7LS0NDIyMjDGFJfzLfr6xx9/TJ8+fcr34sWXsrKcFZ4xMU5vfezY4E7qAP/4B9xzj5L6uQjOxF7KnnVFmTlzJtnZ2WRmZhIeHk50dDSHDh0CnCJc06ZNY8eOHQwZMgRwfhHcf//9DB8+/KTrZGVlFSfe8iipPK8XNX6kcsrPh6uugsaNYfFiZ456jRpeR3V2VaseHyJ6+mlo0sQZg5ez0xh7Cfbu3ctFF11EeHg4ixcvPqnM7vXXX897773HihUruOaaawC45ppreOGFF4p3Pdq+fXtxmV43fPDBB+zevZuDBw8yb948UlNTT1vO9/zzz+fnn392rW2pvA4ccCoyVq0Kzz8PU6d6HdG5OXoUZsxwhmWkdIKzx+6x9PR0rr32WpKTk0lISKBFixbF36tWrRpXXHEFF154IVWqVAGgW7dubNiwgQ4dOgDOtMkZM2YUf/90nnrqKf7+97+zY8cO4uPj6dmzJ1OmTPnFeR07dmTAgAFs3ryZm2++meRkp7hbSeV8wZltExcXR48ePXj88cfL/4ZIpfP993D55U5ZgFtugS5dvI7o3IWHw4IF4MIfvyFDZXvLqKCggKSkJF577TViYmK8DickVJbPRjA5dgyGDIHhw+E3v/E6Gvfs3essoBo/HmJjvY4m8FS2twKsX7+eZs2a0bVrVyV1CTpZWZCeDj/9BFWqwPTp/krqALt3Q2ams0G2nJ6GYsogNja2TBtRiwTS9u3Ois0vv4TCPc59p2lTZ3570c3fol2d5GTqsYtUYj/9BG+/7TxPTYUtW/yb1IsUJfUlS6BTJ6cXLydTYhepxB58EG64AYomYZ1/vrfxBFJ+Phw86NSOl5MpsYtUMvn5zk1EgHHjnLrpF13kbUxeuOoqWLECfvUrZ0hGW+wdp8QuUolYCz17wo03Os8jI6Fwlm1ICivMYH/4g7OHarAUM/OaEvsJAl22Nz09nebNmxMXF8eQIUM4WkKXQ2V75UTGwE03wYABuml4ombNnBurYcpogBJ7wJRUtjc9PZ2NGzfy5ZdfcvDgwRIXJ5WVErv/7NsHgwbBBx84x0OGONMa5biRI+Hxx51fdj/95Pw1E8qU2EsQqLK9PXv2xBiDMYb27duzbdu2EuNR2d7QFhbm7EW6fr3XkQS/7Gxo2xYe9bbclOeCch77Y8sfY+Nud1cgtIhswX3t7yvVuYEu23v06FFeeumlk35pnEhle0NPfr5T32XIEKhZE1auhGrVvI4q+NWtC927V+4SCm4IysTutUCX7R0xYgSdO3emU6dOJX5fZXtDz+LFzjZ1F17o3ChVUi+dsLCTN+fYvh0aNvQuHq8EZWIvbc+6ogSybO+4cePIzs7mP//5z2nPUdne0LF1q1Ne9+qrnZ2OQnnGS3m98w6kpTmrca+4wutoAqvcY+zGmBrGmOXGmC+MMeuMMePO/lPBLVBle6dMmcL777/PrFmzCDvD7XyV7Q0Njz4KcXFOcgcl9fLq1AnuuOP4VnuhxI0e+2HgSmvtPmNMOJBhjHnXWvupC9f2RKDK9t522200adKk+OfS0tIYM2bML85T2V5/K6p3cuONzth6gwZeR+QP55/v7MIEcOQIfPcdNG/ubUyB4mrZXmNMTSADuN1a+9npzlPZXimLyvLZOBd//Sts2wb//rfXkfjb8OHw+uvw1VcQEeF1NOcuoGV7jTFVjDGrgV3AB2dK6pWZyvaK2/bvd+Zd5+d7HYm/3Xuv03uvzEm9LFy5eWqtPQYkGGMuBN4wxsRZa9eeeI4xZhgwDCCqaCPDSkZle6W8Cgrg2WedCoyJic6GEZVhtWRkZCR5eXkBaSsiIoLdLpdsvPRS5wGwbp0zLfLii11tIqi4+pGy1u4BlgDdS/jeZGttsrU2uX79+m42K1Jp/PSTM/wyfbpzXBmSOkBeXh7W2oA8KvIXyOHDTq2dwYMrrImgUO4euzGmPnDUWrvHGHMecBXwWLkjE/EJa50pdz16OPPSP/vMmdIogVe9uvNLtUkTryOpWG70FxoAi40xa4AVOGPsb7lwXRFfmD8fevWCwsoUREWpgJeXLr/cKRgG8MorTi/eb8qd2K21a6y1idbaeGttnLX2ETcCE6nMrIUdO5zn114Lr74KWhQcXDIznUqZZ1gbWGlVkhG+yic6OpqcnJxzOqes5YOnTZvGDz/8UKafkYo1ejQkJ8OePc44et++lWc8PVS0bQsLFzqLmPxGHzUfUGIPHkXLQvr2hbvuCq2t6iqjrl2hShVn39RJk/xT7leJvVBWVhYtWrRg6NChxMXFkZ6ezsKFC0lNTSUmJobly5cDsHv3bq677jri4+NJSUlhzZo1AOTm5tKtWzcSExMZPnz4SbVcZsyYQfv27UlISGD48OEcK8U2L/fccw9JSUl07dqV7OxsAFavXk1KSgrx8fFcf/315OXlMWfOHFauXEl6ejoJCQkcPHiwAt4dOZv8fBg40NmqDpze+n33OUlDgt///i/cfbezgMkPgrII2I4JEzi8wd2yvdVbtuDiBx444zmbN2/mtddeY/LkybRr146XX36ZjIwM5s+fz4QJE5g3bx5jx44lMTGRefPm8eGHHzJw4EBWr17NuHHj6NixI2PGjOHtt99m8uTJgLNq8pVXXuGTTz4hPDycESNGMHPmTAYOHHjaOPbv309SUhL/+Mc/eOSRRxg3bhzPPPMMAwcO5Omnn6ZLly6MGTOGcePGMXHiRJ555hmeeOKJ4lIDEnhVq0J4uPNfqXz+/GdnGqRfSg7oY3iCpk2b0rp1awBatWpF165dMcbQunVrsrKyAMjIyGDu3LkAXHnlleTm5rJ3716WLl3K66+/DkCvXr2IKFzitmjRIjIzM2nXrh0ABw8e5KKz7DwcFhbGjTfeCED//v1JS0tj79697Nmzhy6FhaYHDRpE37593X0DpEx+/hkefBBGjXJmukyZotkulVVYGBT+0+f99yE3F26+2duYyiMoE/vZetYVpXr16sXPw8LCio/DwsLIL1zzXVJtnaKyuqeW1y06f9CgQfztb38757hKuq54LyfHmROdkOBsiKH/TZWftfDPfzqJ/cYbK+9QmsbYy6hz587MnDkTgCVLllCvXj3q1Klz0tfffffd4tVzXbt2Zc6cOcVlfHfv3n1SGeCSFBQUMGfOHABefvllOnbsyAUXXEBERAQff/wxAC+99FJx712legPnwAGYPdt53rQpfPutk9TFH4yB115z9petrEkdgrTHHswefvhhBg8eTHx8PDVr1mR64drwsWPH0q9fP5KSkujSpUtxPZzY2FjGjx9Pt27dKCgoIDw8nEmTJtHkDEvfatWqxbp162jbti0XXHABr7zyCgDTp0/ntttu48CBA1xyySVMnToVcMr33nbbbZx33nksW7aM8847r4LfhdA1aZJTUCo+HmJjnZoj4i9FM5mOHYM//cm5KV7Zbl+5Wra3tCpz2V4JPK8/G4cPw86dzjj64cOwfLmziUMoCeSuXcGyQ1h2Nvz6184Whffe63U0jtKW7VWPXeQs+vRxVpF+/rlTayTUknqoql8fVq2CCy7wOpKy0xi7SAmOHj2+WOWee2DChMo95irnpiipb97s7ENbVCYi2Cmxi5wiJ8fZb3TKFOe4WzeniJeErpwc+PrrypPYNRQjcorISGjWzPlTXAQgJcVZlVqtmnNcUBDctX+CODSRwPnmG/jd744X7Zo9G667zuuoJJgUJfXnnnP+ggvmcr9K7CJAXh4sXgxr1579XAlt1ao5j4ICryM5PSX2Uhg6dCjr16/3OoxiEyZMOOP3e/bsyZ49ewISy8aNG0lISCAxMZFvvvkmIG265fvvYcYM53lyMmzZAh07ehuTBF5kZCTGmFI/br3VMH++oWbN0v/MiY/IyMgKf02ax14J1a5dm3379v3i60V7RoYFaPDv2LFjPP744xw8eJBxRWUNK0BFfTaGDXNWGWZlVc4pbYHk53nslam90s5jV4+9UFHZ3kGDBhEfH88NN9zAgQMHALj88ssp+kW0YMECOnToQFJSEn379mXfvn2sXLmShIQEEhISaN26dXFtl5LK7BZd709/+hOdO3emZcuWrFixgrS0NGJiYvjLX/5SHFNJ5X5Hjx7NwYMHSUhIID09naysLFq2bMmIESNISkpi69atJ23g8eKLLxIfH0+bNm0YMGDAL173Rx99VBx7YmIiP//8M0uWLKF3797F59x5551MmzYNcDYHeeSRR+jYsSOvvPIKEydOZMqUKVxxxRUAXHfddbRt25ZWrVoVV7gEeO+990hKSqJNmzZ07doVcKpYDhkyhHbt2pGYmMibRXvHVaAff4Tt253njz3mzE1XUhffCdTO4yc+2rZta0+1fv36k467dLF26lTn+ZEjzvFLLznH+/c7x7NnO8d79jjHc+c6x9nZzvH8+c7xjz/+orlf+O677yxgMzIyrLXWDh482D7++OOFsXSxK1assNnZ2bZTp05237591lprH330UTtu3LiTrjNq1Cg7atQoa621rVu3tkuWLLHWWvvQQw/ZkSNHFl/v3nvvtdZaO3HiRNugQQP7ww8/2EOHDtmGDRvanJwcu379etu7d2975MgRa621t99+u50+fbq11tpatWqdFLcxxi5btqz4a02aNLHZ2dl27dq19rLLLrPZ2dnWWmtzc3N/8bp79+5d/Jp//vlne/ToUbt48WLbq1ev4nPuuOMOO7Xwf0aTJk3sY489Vvy9sWPHFr9PJ7Zx4MAB26pVK5uTk2N37dplGzVqZL/99tuTzrn//vvtS4X/U/Py8mxMTEzxe3uiUz8b5+rwYWujoqy99lpXLhdSnFThv7YqW3vASluKHKvpjido3LgxqampgFMu96mnnmLUqFHF3//0009Zv3598TlHjhyhQ4cOxd9/9dVX+fzzz1mwYMFZy+z2KdwAs3Xr1rRq1YoGDRoAcMkll7B161YyMjJKXe63SZMmpKSk/OLrH374ITfccAP16tUDKHFsLzU1lbvvvpv09HTS0tJo1KjRWd+nopLCJXnqqad44403ANi6dStff/012dnZdO7cmaaFOwgXxbFgwQLmz5/PE088AcChQ4f4/vvvXR922b8fatVybng99RRoxE/8LmgT+5Ilx5+Hh598XLPmyccXXHDycb16Jx9ffHHp2jy1PO6px9Zarr76ambNmvWLn123bh1jx45l6dKlVCnFEsUTSwKfWi44Pz+/TOV+a9WqVeLXrbVnLfk7evRoevXqxTvvvENKSgoLFy6katWqFJxwy//QoUOlam/JkiUsXLiQZcuWUbNmTS6//HIOHTp02jistcydO5fmFbi7wdq1zvZnU6c6Gyn89rcV1pRI0NAY+wm+//57li1bBsCsWbPoeMoUiZSUFD755BM2b94MwIEDB/jqq6/Yu3cvN910Ey+++CL1C1e1nKnMbmmcqdxveHg4R48eLdU1Xn31VXJzc4uvcapvvvmG1q1bc99995GcnMzGjRtp0qQJ69ev5/Dhw+zdu5dFixaVKua9e/cSERFBzZo12bhxI59++ikAHTp04KOPPuK77747KY5rrrmGp59+uvhG0qpVq0rVTlnExMBVV0Hjxq5fWiRoBW2P3QstW7Zk+vTpDB8+nJiYGG6//faTvl+/fn2mTZtGv379OFy4OmH8+PEsW7aMLVu28Pvf/7743NWrV5+2zG5pnKnc77Bhw4iPjycpKYm//vWvp71Gq1atePDBB+nSpQtVqlQhMTGx+CZokYkTJ7J48WKqVKlCbGwsPXr0oHr16vzud78jPj6emJgYEhMTSxVz9+7dee6554iPj6d58+bFw0P169dn8uTJpKWlUVBQwEUXXcQHH3zAQw89xB//+Efi4+Ox1hIdHc1bb71V6vfodBYtgiefhNdfd4p2FZbJFwkZmu5YKCsri969e7NWK1SCTlk/G2+95Wwk/c47cIay91IGmu4YHO1puqOElI8+gsL9SOjdG774QkldQpeGYgpFR0ert15JWQvjxzv7VPbt69R6qapPtoSwoPr4l2YWh4SWM/3J+tln0Lw5XHihUxqgdu3grrgnEihB88+gRo0a5ObmBsWWWBIcrLXk5uZSo0aNX3xvxw7o3BkeecQ5/p//ceaqi0gQ9dgbNWrEtm3byM7O9joUCSI1atQ4adFUTo6zTuHii50x9cJKBiJygqCZFVORxv13Het/+Clg7UnF2P5FBJ9OieGKUeuIbLLf63BCykcffVSmdRiVpa3ytBf7qzqMvbZVmX9Os2JEOL73aP1mP9H0N9nUqhvEOxyIBIFy99iNMY2BF4GLgQJgsrX2X2f6mUD32CW45Oc7G0Mb4+xY9OOPcNllzte+/dYpA9C7t3MjdO5cmDYN5s93zhdvaB57cLQXyB57PnCPtbYlkALcYYyJdeG6EiBHj8KxY87zgwdhwwYoKveene0k1cKqBHz3nbOqs2hT39Wr4fbbj5fCXbzY2c192zbneO5ciI4+fvyf/zi1f3budI5ffBFiY2HvXuf49dedei6FFZPZtMmZurhfIy8ipVbum6fW2h+BHwuf/2yM2QA0BIJny6F3R8OOL72OokLs3Hshs5d35vLmX9Im6juyf6rDmHn9GfibD+nQbCPf59bn1qkjub/Xq1zZcg3rtzem58RxPDdwEt1bZ/LZN81J+euTvP3HsfSMX8nKr1rR+dG/88E9D3JVq9WsWd+G3z4xgY/uu5fOzdex6cu23PPPR/jN7ru5+NJN/LimLXNn3s2IXz1Aw0ZbKFjfhgNZ/Tn68hNQfyf/81UsXRpeQ/icqXDBHtr9eCn/7/pkzpv7X6h5gGv2NWT2bZdQc86nEH6Um8MiuWJMJDVmfwtVCrinXlWq98mH17x+p0Pb4kE1YWov37VVrvYubg09HnU/IBe4OivGGBMNJAKflfC9YcAwgKioKDebDWkHj1bjoTcG8OgN02gT9R1HjoUzNzOVzpeto0OzjYSZAg4cqU7+Mafi5AU1D9Cl+ZfUre3cTI6qu4tHrnuJZhf9AEDLBluZNfwx4ho5BcfaNf2KlWNG0vxip8t9Zcsv2PNMX2rXcCo+9ojPZNe/0ovj6Rr7BV1jvyg+7njZejpedvx3fFKTb0hqcnwLveYNttO8wfbi419F7OZXEceLlVUPz3fnjRIJIa7NijHG1AY+Av5qrX39TOdqjL18Dh50hizSC/Ppnj3O4hyttpSKojH24GgvoLNijDHhwFxg5tmSupTfv/8N/fvDmjXO8YUXKqmLyHHlTuzGqQHwPLDBWvtk+UOSkhQUHL/heNdd8PHHEB/vbUwiEpzc6LGnAgOAK40xqwsfPV24rpxgyBC4/HI4dMiZVXLKHiAiIsXcmBWTAWiGcQXr3x9SUpyNI0REzkQrT4PUkSNwxx3w7LPO8VVXwW23aZGOiJydbrkFqfBwyMpyboyKiJSFEnsQKShwVmb+7ndQty68+aZmu4hI2WkoJoh89RWMHAkvvOAcK6mLyLlQ6ggCmzdDs2bQogWsWKFpjCJSPuqxe2zGDCehL1/uHLdpoxukIlI+SuweKVpR3KcPPPywk9D9LjIyEmNMwB6RkZFev2QRTyixe+C556B7d6dUbp068Je/hMb89Ly8PKy1AXvk5eV5/ZJFPKHE7oEaNZxHUc1xERE3KbEHyFtvwTvvOM8HDYJ58+D8872NSUT8SbNiAuDYMRg7FiIioEcP3RwVkYqlxF6B1q2DSy91hl3mz4d69ZTURaTiaSimgmzdCm3bwvjxznHDhqFxg1REvKceu8vy850Vo40bOwW8rr3W64hEJNSox+6iTz6BmBjYtMk5HjIE6tf3NiYJDL/P0Y+IiAjYa4uIiAjoa/Mj9dhdFB0Nl1zidRTihaI5+oFiAnyzZvfu3Wc/SYKGeuzltGaNs8DIWmccfdEiaN7c66gyYYpsAAAIyUlEQVREJJQpsZfTf/8Lzz8PO3Z4HYmIiEOJ/Rzs2AFffuk8Hz0a1q6FBg28jUlEpIjG2MvIWmemy6FD8MUXUKWKsymGiEiwUGIvpUOHoFo1CAuDSZOgdm3nuYhIsFFqKoXcXEhOhiefdI7bt4fYWG9jEhE5HSX2UoiMhE6dtLORiFQOSuynsWMH9O8PO3c69V3+/W/o1s3rqEREzk6J/TRycpwyuytXeh2JiEjZKLGf4OBBePNN53lcHGzZAr16eRuTiEhZKbGf4O9/h7Q02LzZOdZGGO4KZL0R1RyRUBby0x0LCmDPHucG6ahR0LkzNGvmdVT+pHojIoER8j32/v2dXY3y86FWLbjiCq8jEhEpn5DvsaelQXa2s4JURMQPXEnsxpgXgN7ALmttnBvXrCiHDztDLh07wo03wg03eB2RiIi73BqKmQZ0d+laFapKFcjMhPXrvY5ERKRiuNJjt9YuNcZEu3GtimAtvPAC3HSTM46+eLH2HxUR/wqJMfaRL03j6d8P5H8/n0nLaz/0OhzxoaajmzL4vcG+bc/PzvW9bBHZgvva31cBEZVfwBK7MWYYMAwgKioqUM0C0LDFTnr+YwL1W3wT0HZFRLxg3NqnsXAo5q3S3DxNTk62K7VWX3zEGBPwPU8D2Z6fVab/d8aYTGtt8tnOC/l57CIifuNKYjfGzAKWAc2NMduMMbe6cV0RESk7t2bF9HPjOiIiUn4aihER8RkldhERn1FiFxHxGSV2ERGfUWIXEfEZJXYREZ9RYhcR8RkldhERn1FiFxHxGSV2ERGfUWIXEfEZJXYREZ9RYhcR8RkldhERn1FiFxHxmZDYzLoyiYyMJC8vLyBtRUREsHv37oC05XcREREYYwLansjpKLEHmby8vIDtvxjIROR3+gUpwURDMSIiPqPELiLiM0rsIiI+o8QuIuIzSuwiIj6jxC4i4jNK7CIiPqPELiLiM0rsIiI+o8QuIuIzSuwiIj6jxC4i4jNK7CIiPqPELiLiM64kdmNMd2PMJmPMZmPMaDeuKSIi56bcid0YUwWYBPQAYoF+xpjY8l5XRETOjRsbbbQHNltrvwUwxswGfgusd+HartgxYQKHN2z0OoxSmdY4ii0DBvquLZFgda7/Dqq3bMHFDzxQARGVnxtDMQ2BrSccbyv82kmMMcOMMSuNMSuzs7NdaFZEREpiyrsNmzGmL3CNtXZo4fEAoL219q7T/UxycrJduXJludr1K2NMQLfGC1RbIsEq0P8OytOeMSbTWpt8tvPc6LFvAxqfcNwI+MGF64qIyDlwI7GvAGKMMU2NMdWAm4D5LlxXRETOQblvnlpr840xdwLvA1WAF6y168odmYiInBM3ZsVgrX0HeMeNa4mISPlo5amIiM8osYuI+IwrQzEiIpVVREQExpiAtlfRlNhFJKTt3r3b6xBcV+mGYiIjIzHGBOwRGRnp9UsWESmTStdjz8vLC/gqMRGRyqTS9dhFROTMlNhFRHxGiV1ExGeU2EVEfEaJXUTEZ5TYRUR8RoldRMRnlNhFRHym0i1Q8rtA1q0IRM0KEQk8JfYg48e6FSISWBqKERHxGSV2ERGfUWIXEfEZJXYREZ9RYhcR8RkldhERn1FiFxHxGSV2ERGfUWIXEfEZJXYREZ9RYhcR8RkldhERn1FiFxHxGSV2ERGfKVdiN8b0NcasM8YUGGOS3QpKRETOXXl77GuBNGCpC7GIiIgLyrXRhrV2AxCwHX9EROTsNMYuIuIzZ+2xG2MWAheX8K0HrbVvlrYhY8wwYBhAVFRUqQM8VSD3BC1qT0SkMjlrYrfWXuVGQ9baycBkgOTkZHuu19GeoCIiZ6ahGBERnynvdMfrjTHbgA7A28aY990JS0REzlV5Z8W8AbzhUiwiIuICDcWIiPiMEruIiM8osYuI+IwSu4iIzyixi4j4jLH2nNcKnXujxmQDWwLcbD0gJ8BtBoqfXxv4+/XptVVeXry+Jtba+mc7yZPE7gVjzEprrS9LC/v5tYG/X59eW+UVzK9PQzEiIj6jxC4i4jOhlNgnex1ABfLzawN/vz69tsoraF9fyIyxi4iEilDqsYuIhISQSOzGmO7GmE3GmM3GmNFex+MWY8wLxphdxpi1XsfiNmNMY2PMYmPMhsIN00d6HZObjDE1jDHLjTFfFL6+cV7H5DZjTBVjzCpjzFtex+ImY0yWMeZLY8xqY8xKr+Mpie+HYowxVYCvgKuBbcAKoJ+1dr2ngbnAGNMZ2Ae8aK2N8zoeNxljGgANrLWfG2POBzKB6/zw/w3AONuA1bLW7jPGhAMZwEhr7aceh+YaY8zdQDJQx1rb2+t43GKMyQKSrbVBO0c/FHrs7YHN1tpvrbVHgNnAbz2OyRXW2qWAL7eUstb+aK39vPD5z8AGoKG3UbnHOvYVHoYXPnzTyzLGNAJ6AVO8jiUUhUJibwhsPeF4Gz5KEKHAGBMNJAKfeRuJuwqHKlYDu4APrLV+en0TgXuBAq8DqQAWWGCMySzcyznohEJiL2nna9/0jPzOGFMbmAv80Vr7k9fxuMlae8xamwA0AtobY3wxnGaM6Q3sstZmeh1LBUm11iYBPYA7CodEg0ooJPZtQOMTjhsBP3gUi5RB4djzXGCmtfZ1r+OpKNbaPcASoLvHobglFehTOBY9G7jSGDPD25DcY639ofC/u3B2kGvvbUS/FAqJfQUQY4xpaoypBtwEzPc4JjmLwpuLzwMbrLVPeh2P24wx9Y0xFxY+Pw+4CtjobVTusNbeb61tZK2Nxvn39qG1tr/HYbnCGFOr8GY+xphaQDcg6Gal+T6xW2vzgTuB93FuwL1qrV3nbVTuMMbMApYBzY0x24wxt3odk4tSgQE4vb3VhY+eXgflogbAYmPMGpzOxwfWWl9NC/Sp/wEyjDFfAMuBt62173kc0y/4frqjiEio8X2PXUQk1Cixi4j4jBK7iIjPKLGLiPiMEruIiM8osYuI+IwSu4iIzyixi4j4zP8HQx+ElyyOO70AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11e1af6a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.plot(m.dis.top.array[r, c], label='model top')\n",
    "for i, l in enumerate(m.dis.botm.array[:, r, c]):\n",
    "    label = 'layer {} bot'.format(i+1)\n",
    "    if i == m.nlay -1:\n",
    "        label = 'model bot'\n",
    "    plt.plot(l, label=label)\n",
    "plt.plot(heads[0], label='piezometric surface', color='b', linestyle=':')\n",
    "for iw in range(len(sctop)):\n",
    "    ax.fill_between([iw-.25, iw+.25], scbot[iw], sctop[iw], \n",
    "                    facecolor='None', edgecolor='k')\n",
    "ax.legend(loc=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### example of transmissivites without `sctop` and `scbot`\n",
    "* well zero has T=0 in layer 1 because it is dry; T=0.2 in layer 2 because the sat. thickness there is only 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0. , 0. , 0.1, 0.2, 0.2, 0.2],\n",
       "       [0.2, 2. , 2. , 2. , 2. , 2. ],\n",
       "       [2. , 2. , 2. , 1.2, 2. , 2. ]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = flopy.utils.get_transmissivities(heads, m, r=r, c=c)\n",
    "np.round(T, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}