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.8.6 | packaged by conda-forge | (default, Oct  7 2020, 18:42:56) \n",
      "[Clang 10.0.1 ]\n",
      "numpy version: 1.18.5\n",
      "matplotlib version: 3.2.2\n",
      "flopy version: 3.3.3\n"
     ]
    }
   ],
   "source": [
    "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's coordinate information 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 0x1296806d0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3hU1dn38e8KBBAQTQAfKZAEJQIhCUkINDQcVBQ5STWX1GI4CFJQ1NIqVdQK4kupVmupSmt5UEBBUEGRekQQxPigQASRk4oaBBTMCRQIh5D1/rGTcAoQyGT27Jnf57rmcvbMZK97hvHOytpr3ctYaxEREe8KczsAERGpGiVyERGPUyIXEfE4JXIREY9TIhcR8biabjTaqFEjGxMT40bTIiKelZ2dnWetbXzi464k8piYGFavXu1G0yIinmWM2VrR4xpaERHxOCVyERGPUyIXEfE4V8bIK3L48GG2b9/OgQMH3A5FAkidOnVo1qwZ4eHhbociErACJpFv376d888/n5iYGIwxbocjAcBaS35+Ptu3b6dFixZuhyMSsHw2tGKMqWGMWWOMeeNcfv7AgQM0bNhQSVzKGWNo2LCh/koTOQNfjpGPBjZV5QRK4nIifSdEzswnidwY0wzoA0zzxflExPsOHoTiYrejCA2+6pFPBu4BSk71AmPMCGPMamPM6tzcXB81G7hiYmLIy8s7p9dMmjSpusIS8YviYrj2Whgxwu1IQkOVE7kxpi/wo7U2+3Svs9ZOtdamWmtTGzc+aYWpHEOJXLxo3z4YOBBefhlq1oQuXaBrV7ejCg2+6JGnA/2MMTnAXOBKY8wsH5zXr3JycmjdujXDhw8nPj6ezMxMFi9eTHp6OrGxsaxcuRKAgoICrrvuOhITE0lLS2PdunUA5Ofn06NHD5KTkxk5ciTH7rw0a9YsOnbsSFJSEiNHjuTIkSOnjGPs2LEUFRWRlJREZmYmAE888QTx8fHEx8czefLk4+IdMmQIiYmJ3HDDDezfv7+6Ph6RM6pVC7Zuhe3bneMHH4Sbb3bu//ija2GFBOPLrd6MMZcDY6y1fU/3utTUVHtirZVNmzbRpk0bACb8dwMbv//JZ3EBxP2iAeOvbXvK53NycmjZsiVr1qyhbdu2dOjQgXbt2vHss8+ycOFCpk+fzoIFC7jzzjtp1KgR48eP5/333+euu+5i7dq1/P73v6dRo0aMGzeON998k759+5Kbm0tubi733HMPr776KuHh4YwaNYq0tDQGDx5cXnOmUaNGx8VSv3599u7dC0B2djY333wzH3/8MdZafvnLXzJr1iwiIiJo0aIFWVlZpKenM2zYMOLi4hgzZoxPP7dAcOx3QwJLcTFMmQK/+x3Uresc1zxhUvPmzZCWBo8/DsOHuxNnsDDGZFtrU098XCs7j9GiRQsSEhIICwujbdu2dO/eHWMMCQkJ5OTkAJCVlcWgQYMAuPLKK8nPz2fPnj0sX76cgQMHAtCnTx8iIiIAWLJkCdnZ2XTo0IGkpCSWLFnCN998U+mYsrKyuP7666lXrx7169cnIyODDz/8EIDmzZuTnp4OwMCBA8nKyvLVRyFSKatXwx//CC+95ByfmMQBWrZ0euZXX+3X0EKKTxcEWWuXAcuqep7T9ZyrU+3atcvvh4WFlR+HhYVRXHr5vaK/YMqmyFU0Vc5ay5AhQ/jrX/96TjGd7i+mE9vTVD3xly1bnASdlgaffgpJSad+bc2aUDoiiLXw1Vdw2WX+iTNUqEd+lrp27crs2bMBWLZsGY0aNaJBgwbHPf72229TWFgIQPfu3Zk3bx4/lg4SFhQUsHVrhZUoy4WHh3P48OHy9hYsWMD+/fvZt28fr732Gl26dAHgu+++Y8WKFQDMmTOHzp07+/4Ni5xg8mSIj4cvvnCOT5fETzRlCiQmwuefV09soSpgluh7xUMPPcTQoUNJTEykbt26zJw5E4Dx48czYMAAUlJS6NatG1FRUQDExcUxceJEevToQUlJCeHh4UyZMoXo6OhTtjFixAgSExNJSUlh9uzZ3HzzzXTs2BGA4cOHk5ycTE5ODm3atGHmzJmMHDmS2NhYbrvttur/ACRkWQvGwE03waFDcOmlZ3+OAQNg715o684f3UHLpxc7K+tMFzvlzHJycujbty/r1693O5Rqp++G+yZNgjVrnKmFvhrBy8uD/Hxo1co35wsFutgpIuesdm2oU8fpifvKgAHQpw+UjiJKFWhoxaNiYmJCojcu7nnxRYiKgs6d4a67fNcTL/OPfzg9clUorjolchE5yYEDMG4cdOjgJPLqmBAVH3/0/pIlkJwMkZG+bycUaGhFRMqtWQNHjjjDKO+/Dy+8UP1t5uXBr38NY8dWf1vBSolcRABnSmBqKjz1lHMcFVXxAh9fa9QI/vtfZ+WnnBslcpEQV1ZqNj7emec9dKj/Y7jiCmjQ4OiSf5W/PTtK5MeoX7++X9t7+umnadmyJcaYU5a8nTFjBnfccUelz7l7927+9a9/+SpECXJvv+2sstyxwxkHv/VWuOAC9+J56y244w4nLqk8JXI/sdZSUnJ8ufb09HQWL1582sVBZ0uJXM7GpZdCbKwzLh4I+vWDTz5xaplL5SmRV2Dv3r10796dlJQUEhISeP311wF48MEH+ec//1n+ugceeIAnn3wSgMcee4wOHTqQmJjI+PHjAcpXX44aNYqUlBS2bdt2XDvJycnExMScMZ5t27bRs2dPWrVqxYQJE8ofr6i87dixY/n6669JSkriT3/6U5U+BwlO//d/zgIfcHrj777rjIcHitJFzHzxxdHxejm9wJx++PZY2OnjYgwXJ0CvRyr10jp16vDaa6/RoEED8vLySEtLo1+/ftxyyy1kZGQwevRoSkpKmDt3LitXrmTRokV89dVXrFy5Emst/fr1Y/ny5URFRfHFF18wffr0KvWSV65cyfr166lbty4dOnSgT58+GGOYPn06n3zySXl5227duvHII4+wfv161q5de87tSXCbNw8WLHCGMBo0cDuaU/vXv2DuXKckQMOGbkcT2AIzkbvMWsv999/P8uXLCQsLY8eOHezatYuYmBgaNmzImjVr2LVrF8nJyTRs2JBFixaxaNEikpOTAadH/9VXXxEVFUV0dDRpaWlViufqq6+mYek3OSMjg6ysLIwx5eVtyx7/8MMP6devX9XevASlnBxnBWVsrNMbHz8+sJM4wN//DnffrSReGYGZyCvZc64us2fPJjc3l+zsbMLDw4mJieHAgQOAU7RqxowZ7Ny5k2HDhgFO4r/vvvsYOXLkcefJyckpT7RVUVG5Wjdq5Ig3FRfDVVdB8+awdKkzR7xOHbejOrOaNY8O+Tz1FERHO2PocjKNkVdgz549XHTRRYSHh7N06dLjys5ef/31vPPOO6xatYprrrkGgGuuuYbnnnuufFefHTt2lJet9YX33nuPgoICioqKWLBgAenp6acsb3v++efz888/+6xt8a79+52KhTVrwrPPwvTpbkd0bg4fhlmznGEWqVhg9shdlpmZybXXXktqaipJSUm0bt26/LlatWpxxRVXcOGFF1KjRg0AevTowaZNm+jUqRPgTGOcNWtW+fOn8uSTT/K3v/2NnTt3kpiYSO/evZk2bdpJr+vcuTODBg1iy5Yt3HTTTaSmOsXPKipvC85smPj4eHr16sVjjz1W9Q9EPOe77+Dyy51l9jffDN26uR3RuQsPh0WLwAd/3AYtlbE9SyUlJaSkpPDKK68QGxvrdjghwSvfjUBy5AgMGwYjR8KvfuV2NL6zZ4+zYGniRIiLczsa/1MZWx/YuHEjLVu2pHv37kriEnByciAzE376CWrUgJkzgyuJAxQUQHa2s6GzHKWhlbMQFxd3Vhsni/jTjh3OisjPP4fSPbmDTosWzvzysou1ZbsWhTr1yEU87Kef4M03nfvp6bB1a/Am8TJlSXzZMujSxemlhzolchEPe+ABuOEGKJskdf757sbjT8XFUFTk1E4PdUrkIh5TXOxc9AOYMMGpG37RRe7G5IarroJVq+AXv3CGWEJ5yzglchEPsRZ694Ybb3TuR0ZC6azXkBRWmsF+/3tnD9BAKf7lb0rkx/B3GdvMzExatWpFfHw8w4YN43AFXQqVsZVjGQO//S0MGqSLfMdq2dK5EBoWohktRN+2/1VUxjYzM5PNmzfz+eefU1RUVOFioLOlRB589u6FIUPgvfec42HDnGmGctTo0fDYY84vt59+cv5aCSVK5BXwVxnb3r17Y4zBGEPHjh3Zvn17hfGojG1oCwtz9tLcuNHtSAJfbi60bw+PuFuuye8Cch75oysfZXOBb2f8t45szb0d763Ua/1dxvbw4cO88MILx/2SOJbK2Iae4mKnPsqwYVC3LqxeDbVquR1V4GvYEHr29HZJgnMRkIncbf4uYztq1Ci6du1Kly5dKnxeZWxDz9KlzrZrF17oXNhUEq+csLDjN6PYsQOaNnUvHn8JyERe2Z5zdfFnGdsJEyaQm5vLf/7zn1O+RmVsQ8e2bU652auvdnbyCeUZKVX11luQkeGsdr3iCrejqV4aI6+Av8rYTps2jXfffZc5c+YQdprL7SpjGxoeecTZyb7sUoqSeNV06QK3335067hgVuUeuTGmDrAcqF16vnnW2vFVPa+b/FXG9tZbbyU6Orr85zIyMhg3btxJr1MZ2+BWVi/kxhudsfEmTdyOKDicf76zyxDAoUPw7bfQqpW7MVWXKpexNc7f/fWstXuNMeFAFjDaWvvxqX5GZWzlbHjlu3Eu/vIX2L4d/v1vtyMJbiNHwquvwpdfQkSE29Gcu2orY2sde0sPw0tvQTmAqzK24mv79jnznouL3Y4kuN1zj9M793ISPx2fXOw0xtQAsoGWwBRr7ScVvGYEMAIgqmwjPo9RGVupqpISZ3f49HRITnY2SPDCasTIyEgKCwv90lZERAQFPi5peOmlzg1gwwZnmuLFF/u0CVf55CtkrT1irU0CmgEdjTHxFbxmqrU21Vqb2rhxY180K+I5P/3kDKfMnOkceyGJAxQWFmKt9cutOn9hHDzo1KoZOrTamnCFT6cfWmt3G2OWAT2B9b48t4hXWetMgevVy5kX/sknzhRD8b/atZ1fotHRbkfiW1XuDxhjGhtjLiy9fx5wFaCNmERKLVwIffpAaaUHoqJU8MpNl1/uFNgCeOklp5fudb74w64JsNQYsw5YBbxnrX3DB+cV8SxrYedO5/6118LLL4MW3QaW7GynkuRp1uJ5hi9mrayz1iZbaxOttfHW2od9EZjXxcTEkJeXd06vOdtyujNmzOD7778/q5+R6jV2LKSmwu7dzjh4//7eGQ8PFe3bw+LFzqIhr9NXKwgokQeOsmUZ/fvDnXeG1tZrXtS9O9So4ez7OWWKd8vfKpGXysnJoXXr1gwfPpz4+HgyMzNZvHgx6enpxMbGsnLlSgAKCgq47rrrSExMJC0tjXXr1gGQn59Pjx49SE5OZuTIkcfVQpk1axYdO3YkKSmJkSNHcqQS25jcfffdpKSk0L17d3JzcwFYu3YtaWlpJCYmcv3111NYWMi8efNYvXo1mZmZJCUlUVRUVA2fjpxJcTEMHuxsvQZOb/zee50kIYHvf/8X7rrLWTDkRQFZNGvnpEkc3OTb66W127Tm4vvvP+1rtmzZwiuvvMLUqVPp0KEDL774IllZWSxcuJBJkyaxYMECxo8fT3JyMgsWLOD9999n8ODBrF27lgkTJtC5c2fGjRvHm2++ydSpUwFnVeJLL73ERx99RHh4OKNGjWL27NkMHjz4lHHs27ePlJQU/v73v/Pwww8zYcIEnn76aQYPHsxTTz1Ft27dGDduHBMmTGDy5Mk8/fTTPP744+VL98X/ataE8HDnv+I9f/qTMy3Rq0v49bU7RosWLUhISACgbdu2dO/eHWMMCQkJ5OTkAJCVlcX8+fMBuPLKK8nPz2fPnj0sX76cV199FYA+ffoQUbqEbMmSJWRnZ9OhQwcAioqKuOgMO+WGhYVx4403AjBw4EAyMjLYs2cPu3fvpltpoeUhQ4bQv39/334AclZ+/tnZxX7MGGcmyrRpmo3iVWFhUPq/Pu++C/n5cNNN7sZ0NgIykZ+p51xdateuXX4/LCys/DgsLIzi0jXUFdWmKSsze2K52bLXDxkyhL/+9a/nHFdF5xX35eU5c5KTkpwNIPTP5H3Wwj/+4STyG2/0ztCYxsjPUteuXZk9ezYAy5Yto1GjRjRo0OC4x99+++3y1Wndu3dn3rx55WVtCwoKjiuLW5GSkhLmzZsHwIsvvkjnzp254IILiIiI4MMPPwTghRdeKO+dq3St/+zfD3PnOvdbtIBvvnGSuAQHY+CVV5z9Ub2SxCFAe+SB7KGHHmLo0KEkJiZSt25dZpautR4/fjwDBgwgJSWFbt26ldeTiYuLY+LEifTo0YOSkhLCw8OZMmUK0adZWlavXj02bNhA+/btueCCC3jppZcAmDlzJrfeeiv79+/nkksuYfr06YBTzvbWW2/lvPPOY8WKFZx33nnV/CmErilTnAJMiYkQF+fU7JDgUjbT6MgR+OMfnYvYgX75qcplbM+Fl8vYiv+5/d04eBB27XLGwQ8ehJUrnU0LQok/d6UKlB2wcnPhl790tty75x63o3GcqoyteuQiZ9Cvn7NK89NPnVodoZbEQ1XjxrBmDVxwgduRnJnGyEUqcPjw0cUhd98NkyZ5a8xUfKMsiW/Z4uyjWlZ2IdAokYucIC/P2S9z2jTnuEcPp+iVhK68PPjqq8BN5BpaETlBZCS0bOn8aS0CkJbmrPqsVcs5LikJrNo5ARSKiHu+/hp+85ujRa7mzoXrrnM7KgkkZUn8mWecv9ACqfytErkIUFgIS5fCem2HImdQq5ZzKylxO5KjlMgrYfjw4WzcuNHtMMpNmjTptM/37t2b3bt3+yWWzZs3k5SURHJyMl9//bVf2vSV776DWbOc+6mpsHUrdO7sbkzif5GRkRhjKn275RbDwoWGunUr/zPH3iIjI33+HjSP3IPq16/P3r17T3q8bM/DMD8N3h05coTHHnuMoqIiJpSV/asG1fXdGDHCWcWXk+ONKWZuCuZ55F5q71TzyNUjL1VWxnbIkCEkJiZyww03sH//fgAuv/xyyn7xLFq0iE6dOpGSkkL//v3Zu3cvq1evJikpiaSkJBISEspro1RUdrbsfH/84x/p2rUrbdq0YdWqVWRkZBAbG8uf//zn8pgqKn87duxYioqKSEpKIjMzk5ycHNq0acOoUaNISUlh27Ztx21Y8fzzz5OYmEi7du0YNGjQSe/7gw8+KI89OTmZn3/+mWXLltG3b9/y19xxxx3MmDEDcDbDePjhh+ncuTMvvfQSkydPZtq0aVxxxRUAXHfddbRv3562bduWV4AEeOedd0hJSaFdu3Z0794dcKo8Dhs2jA4dOpCcnMzrZXuhVaMffoAdO5z7jz7qzA1XEhfP89fO2Mfe2rdvb0+0cePG4467dbN2+nTn/qFDzvELLzjH+/Y5x3PnOse7dzvH8+c7x7m5zvHChc7xDz+c1NxJvv32WwvYrKwsa621Q4cOtY899lhpLN3sqlWrbG5uru3SpYvdu3evtdbaRx55xE6YMOG484wZM8aOGTPGWmttQkKCXbZsmbXW2gcffNCOHj26/Hz33HOPtdbayZMn2yZNmtjvv//eHjhwwDZt2tTm5eXZjRs32r59+9pDhw5Za6297bbb7MyZM6211tarV++4uI0xdsWKFeWPRUdH29zcXLt+/Xp72WWX2dzcXGuttfn5+Se97759+5a/559//tkePnzYLl261Pbp06f8NbfffrudXvqPER0dbR999NHy58aPH1/+OR3bxv79+23btm1tXl6e/fHHH22zZs3sN998c9xr7rvvPvtC6T9qYWGhjY2NLf9sj3Xid+NcHTxobVSUtdde65PThRQnVQRfW15rD1htK8ipmn54jObNm5Oeng445WOffPJJxowZU/78xx9/zMaNG8tfc+jQITp16lT+/Msvv8ynn37KokWLzlh2tl/pBo4JCQm0bduWJk2aAHDJJZewbds2srKyKl3+Njo6mrS0tJMef//997nhhhto1KgRQIVjc+np6dx1111kZmaSkZFBs2bNzvg5lZXYrciTTz7Ja6+9BsC2bdv46quvyM3NpWvXrrQo3fG2LI5FixaxcOFCHn/8cQAOHDjAd9995/NhlH37oF495wLVk0+CRvAk2ARsIl+27Oj98PDjj+vWPf74gguOP27U6Pjjiy+uXJsnlos98dhay9VXX82cOXNO+tkNGzYwfvx4li9fTo1KLAE8tkTuieVzi4uLz6r8bb169Sp83Fp7xhK4Y8eOpU+fPrz11lukpaWxePFiatasSckxl+QPHDhQqfaWLVvG4sWLWbFiBXXr1uXyyy/nwIEDp4zDWsv8+fNpVY3V/Nevd7bzmj7d2Tjg17+utqZEXKMx8mN89913rFixAoA5c+bQ+YQpDGlpaXz00Uds2bIFgP379/Pll1+yZ88efvvb3/L888/TuHQVyenKzlbG6crfhoeHc/jw4Uqd4+WXXyY/P7/8HCf6+uuvSUhI4N577yU1NZXNmzcTHR3Nxo0bOXjwIHv27GHJkiWVinnPnj1ERERQt25dNm/ezMcffwxAp06d+OCDD/j222+Pi+Oaa67hqaeeKr/ws2bNmkq1czZiY+Gqq6B5c5+fWiRgBGyP3A1t2rRh5syZjBw5ktjYWG677bbjnm/cuDEzZsxgwIABHCxdDTBx4kRWrFjB1q1b+d3vflf+2rVr156y7GxlnK787YgRI0hMTCQlJYW//OUvpzxH27ZteeCBB+jWrRs1atQgOTm5/KJlmcmTJ7N06VJq1KhBXFwcvXr1onbt2vzmN78hMTGR2NhYkpOTKxVzz549eeaZZ0hMTKRVq1blwz2NGzdm6tSpZGRkUFJSwkUXXcR7773Hgw8+yB/+8AcSExOx1hITE8Mbb7xR6c/oVJYsgSeegFdfdYpclZaJFwlamn5YKicnh759+7JeK0ICztl+N954w9n4+K234DRl3+UsaPphYLSn6YcS1D74AEr336BvX/jsMyVxCR0aWikVExOj3rhHWQsTJzr7LPbv79RK0W72EkoC6utemVkWElpO9yfoJ59Aq1Zw4YXOUvv69QOrIp2IvwTM175OnTrk5+cHxBZPEhisteTn51OnTp2Tntu5E7p2hYcfdo7/53+cueIioShgeuTNmjVj+/bt5Obmuh2KBJA6deoct0gpL89ZJ3Dxxc6YeGllAJGQFjCzVqrThP9uYOP3P/mtPakeOz6L4ONpsVwxZgOR0fvcDiekfPDBB2e1DsIrbVWlvbhfNGD8tW3P+uc0a0VCUtl3vnHLn2jxq1zqNQygiv4iAaDKPXJjTHPgeeBioASYaq395+l+xt89cgksxcXORsbGODvy/PADXHaZ89g33zjL6vv2dS5czp8PM2bAwoXO68UdmkceGO1VZ4+8GLjbWtsGSANuN8bE+eC84ieHD8ORI879oiLYtAnKyp3n5jpJtHSVP99+66yaLNuEdu1auO22o6Vhly51dhvfvt05nj8fYmKOHv/nP07tnF27nOPnn4e4ONizxzl+9VWnHkppBWG++MKZSrhPIykip1Tli53W2h+AH0rv/2yM2QQ0BQJnS523x8LOz92Oolrs2nMhc1d25fJWn9Mu6ltyf2rAuAUDGfyr9+nUcjPf5Tfmlumjua/Py1zZZh0bdzSn9+QJPDN4Cj0Tsvnk61ak/eUJ3vzDeHonrmb1l23p+sjfeO/uB7iq7VrWbWzHrx+fxAf33kPXVhv44vP23P2Ph/lVwV1cfOkX/LCuPfNn38WoX9xP02ZbKdnYjv05Azn84uPQeBf/82Uc3ZpeQ/i86XDBbjr8cCn/7/pUzpv/X6i7n2v2NmXurZdQd97HEH6Ym8IiuWJcJHXmfgM1Sri7UU1q9yuGV9z+pEPb0iF1YXqfoGurSu1dnAC9HvF9QOfAp7NWjDExQDLwSQXPjQBGAERFRfmy2ZBWdLgWD742iEdumEG7qG85dCSc+dnpdL1sA51abibMlLD/UG2KjzgVGS+ou59urT6nYX3n4m9Uwx95+LoXaHnR9wC0abKNOSMfJb6ZU6CrQ4svWT1uNK0udrrUV7b5jN1P96d+HaciYq/EbH78Z2Z5PN3jPqN73Gflx50v20jny47+Tk+J/pqU6KNbwrVqsoNWTXaUH/8iooBfRBwt7lU7vNg3H5RIEPPZrBVjTH3gA+Av1tpXT/dajZFXTVGRMwSRWZo/d+92FsNoNaNUF42RB0Z71TprxRgTDswHZp8piUvV/fvfMHAgrFvnHF94oZK4SCirciI3zpr6Z4FN1tonqh6SVKSk5OgFwjvvhA8/hMREd2MSkcDgix55OjAIuNIYs7b01tsH55VjDBsGl18OBw44sz5O2PNCREKYL2atZAGa4VvNBg6EtDRnowQRkWNpZWeAOnQIbr8d/vUv5/iqq+DWW7UoRkROpktkASo8HHJynAuZIiKno0QeQEpKnJWPv/kNNGwIr7+u2SgicmYaWgkgX34Jo0fDc885x0riIlIZShUBYMsWaNkSWreGVas0rVBEzo565C6bNctJ4CtXOsft2umCpoicHSVyl5St0O3XDx56yEngwS4yMhJjjN9ukZGRbr9lEb9QInfBM89Az55O6dgGDeDPfw6N+eGFhYVYa/12KywsdPsti/iFErkL6tRxbmU1t0VEqkKJ3E/eeAPeesu5P2QILFgA55/vbkwiEhw0a8UPjhyB8eMhIgJ69dLFTBHxLSXyarRhA1x6qTOMsnAhNGqkJC4ivqehlWqybRu0bw8TJzrHTZuGxgVNEfE/9ch9rLjYWZHZvLlT8Oraa92OSESCnXrkPvTRRxAb6+z8Dk4N8caN3Y1J/CPY58hHRET47b1FRET49b0FA/XIfSgmBi65xO0oxA1lc+T9xfj5YktBQcGZXySuUY+8itatcxb0WOuMgy9ZAq1auR2ViIQSJfIq+u9/4dlnYedOtyMRkVClRH4Odu6Ezz937o8dC+vXQ5Mm7sYkIsJ+I84AAAifSURBVKFLY+RnyVpnJsqBA/DZZ1CjhrMJhIiIW5TIK+nAAahVC8LCYMoUqF/fuS8i4jalokrIz4fUVHjiCee4Y0eIi3M3JhGRMkrklRAZCV26aOceEQlMSuSnsHMnDBwIu3Y59VH+/W/o0cPtqERETqZEfgp5eU7Z2dWr3Y5EROT0lMiPUVQEr7/u3I+Ph61boU8fd2MSETkTJfJj/O1vkJHh7GoP2vjB1/xZr0M1OySUhPz0w5IS2L3buaA5Zgx07QotW7odVXBSvQ6R6hHyPfKBA51de4qLoV49uOIKtyMSETk7Id8jz8iA3FxnhaaIiBeFXCI/eNAZQuncGW68EW64we2IRESqxidDK8aY54wxPxpj1vvifNWpRg3IzoaNG92ORETEN3w1Rj4D6Omjc/mctU6p2X37nG3Yli6FCRPcjkpExDd8MrRirV1ujInxxbmqw+gXZvDU7wbzv5/Ops2177sdjgShFmNbMPSdoUHbXjA718+ydWRr7u14bzVEdPb8NkZujBkBjACIioryV7MANG29i95/n0Tj1l/7tV0REX8wvtpnsLRH/oa1Nv5Mr01NTbWrtfZdgogxxu97dvqzvWDmpX87Y0y2tTb1xMdDfh65iIjXKZGLiHicr6YfzgFWAK2MMduNMbf44rwiInJmvpq1MsAX5xERkbOnoRUREY9TIhcR8TglchERj1MiFxHxOCVyERGPUyIXEfE4JXIREY9TIhcR8TglchERj1MiFxHxOCVyERGPUyIXEfE4JXIREY9TIhcR8Ti/7dkplRMZGUlhYaFf2oqIiKCgoMAvbQW7iIgIjDF+bU+kjBJ5gCksLPTb/oH+TDzBTr8QxU0aWhER8TglchERj1MiFxHxOCVyERGPUyIXEfE4JXIREY9TIhcR8TglchERj1MiFxHxOCVyERGPUyIXEfE4JXIREY9TIhcR8TglchERj/NJIjfG9DTGfGGM2WKMGeuLc4qISOVUOZEbY2oAU4BeQBwwwBgTV9XziohI5fhiY4mOwBZr7TcAxpi5wK+BjT44t0/snDSJg5s2ux1GpcxoHsXWQYODri2RQHWu/x/UbtOai++/vxoiOnu+GFppCmw75nh76WPHMcaMMMasNsaszs3N9UGzIiICYKq6rZgxpj9wjbV2eOnxIKCjtfbOU/1MamqqXb16dZXaDVbGGL9u9eavtkQClb//P6hKe8aYbGtt6omP+6JHvh1ofsxxM+B7H5xXREQqwReJfBUQa4xpYYypBfwWWOiD84qISCVU+WKntbbYGHMH8C5QA3jOWruhypGJiEil+GLWCtbat4C3fHEuERE5O1rZKSLicUrkIiIe55OhFRERr4qIiMAY49f2fE2JXERCWkFBgdshVJnnhlYiIyMxxvjtFhkZ6fZbFhE5Lc/1yAsLC/2+CktEJJB5rkcuIiLHUyIXEfE4JXIREY9TIhcR8TglchERj1MiFxHxOCVyERGPUyIXEfE4zy0ICnb+rPtQHTUfRMT/lMgDTDDUfRAR/9LQioiIxymRi4h4nBK5iIjHKZGLiHicErmIiMcpkYuIeJwSuYiIxymRi4h4nBK5iIjHKZGLiHicErmIiMcpkYuIeJwSuYiIxymRi4h4nBK5iIjHVSmRG2P6G2M2GGNKjDGpvgpKREQqr6o98vVABrDcB7GIiMg5qNIOQdbaTYDftiYTEZGT+W2rN2PMCGAEQFRU1Dmfx597Wpa1JyISyM6YyI0xi4GLK3jqAWvt65VtyFo7FZgKkJqaaisd4Qm0p6WIyPHOmMittVf5IxARETk3mn4oIuJxVZ1+eL0xZjvQCXjTGPOub8ISEZHKquqsldeA13wUi4iInAMNrYiIeJwSuYiIxymRi4h4nBK5iIjHGWvPeW3OuTdqTC6w1c/NNgLy/NymvwTze4Pgfn96b97lxvuLttY2PvFBVxK5G4wxq621QVmhMZjfGwT3+9N7865Aen8aWhER8TglchERjwulRD7V7QCqUTC/Nwju96f35l0B8/5CZoxcRCRYhVKPXEQkKCmRi4h4XNAncmNMT2PMF8aYLcaYsW7H40vGmOeMMT8aY9a7HYuvGWOaG2OWGmM2lW7wPdrtmHzJGFPHGLPSGPNZ6fub4HZMvmaMqWGMWWOMecPtWHzJGJNjjPncGLPWGLPa7XggyMfIjTE1gC+Bq4HtwCpggLV2o6uB+YgxpiuwF3jeWhvvdjy+ZIxpAjSx1n5qjDkfyAauC6J/OwPUs9buNcaEA1nAaGvtxy6H5jPGmLuAVKCBtbav2/H4ijEmB0i11gbMYqdg75F3BLZYa7+x1h4C5gK/djkmn7HWLgeCcu87a+0P1tpPS+//DGwCmroble9Yx97Sw/DSW9D0qowxzYA+wDS3YwkFwZ7ImwLbjjneThAlg1BhjIkBkoFP3I3Et0qHHtYCPwLvWWuD6f1NBu4BStwOpBpYYJExJrt0U3nXBXsiNxU8FjS9nlBgjKkPzAf+YK39ye14fMlae8RamwQ0AzoaY4JieMwY0xf40Vqb7XYs1STdWpsC9AJuLx3idFWwJ/LtQPNjjpsB37sUi5yl0rHj+cBsa+2rbsdTXay1u4FlQE+XQ/GVdKBf6VjyXOBKY8wsd0PyHWvt96X//RFnh7SO7kYU/Il8FRBrjGlhjKkF/BZY6HJMUgmlFwOfBTZZa59wOx5fM8Y0NsZcWHr/POAqYLO7UfmGtfY+a20za20Mzv9z71trB7oclk8YY+qVXnzHGFMP6AG4PmssqBO5tbYYuAN4F+di2cvW2g3uRuU7xpg5wAqglTFmuzHmFrdj8qF0YBBOb25t6a2320H5UBNgqTFmHU6H4z1rbVBN0wtS/wNkGWM+A1YCb1pr33E5puCefigiEgqCukcuIhIKlMhFRDxOiVxExOOUyEVEPE6JXETE45TIRUQ8TolcRMTj/j95Enawh0xQ9QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}