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": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:28.605674Z",
     "iopub.status.busy": "2022-12-15T13:10:28.605116Z",
     "iopub.status.idle": "2022-12-15T13:10:29.916592Z",
     "shell.execute_reply": "2022-12-15T13:10:29.914338Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.7.15 (default, Oct 18 2022, 09:36:49) \n",
      "[GCC 11.2.0]\n",
      "numpy version: 1.21.6\n",
      "matplotlib version: 3.5.3\n",
      "flopy version: 3.3.6\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "from tempfile import TemporaryDirectory\n",
    "\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": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:29.967462Z",
     "iopub.status.busy": "2022-12-15T13:10:29.966349Z",
     "iopub.status.idle": "2022-12-15T13:10:29.983926Z",
     "shell.execute_reply": "2022-12-15T13:10:29.982929Z"
    }
   },
   "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 = [-0.25, 0.5, 1.7, 1.5, 3.0, 2.5]  # screen tops\n",
    "scbot = [-1.0, -0.5, 1.2, 0.5, 1.5, -0.2]  # screen bottoms\n",
    "# head in each layer, for 6 head target locations\n",
    "heads = np.array(\n",
    "    [\n",
    "        [1.0, 2.0, 2.05, 3.0, 4.0, 2.5],\n",
    "        [1.1, 2.1, 2.2, 2.0, 3.5, 3.0],\n",
    "        [1.2, 2.3, 2.4, 0.6, 3.4, 3.2],\n",
    "    ]\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.0\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": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:29.988455Z",
     "iopub.status.busy": "2022-12-15T13:10:29.987862Z",
     "iopub.status.idle": "2022-12-15T13:10:29.998337Z",
     "shell.execute_reply": "2022-12-15T13:10:29.997267Z"
    }
   },
   "outputs": [],
   "source": [
    "# temporary directory\n",
    "temp_dir = TemporaryDirectory()\n",
    "model_ws = temp_dir.name\n",
    "\n",
    "m = flopy.modflow.Modflow(\"junk\", version=\"mfnwt\", model_ws=model_ws)\n",
    "dis = flopy.modflow.ModflowDis(\n",
    "    m, nlay=nl, nrow=nr, ncol=nc, botm=botm, top=top\n",
    ")\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": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:30.003092Z",
     "iopub.status.busy": "2022-12-15T13:10:30.002401Z",
     "iopub.status.idle": "2022-12-15T13:10:30.011354Z",
     "shell.execute_reply": "2022-12-15T13:10:30.010323Z"
    }
   },
   "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(\n",
    "    heads, m, r=r, c=c, sctop=sctop, scbot=scbot\n",
    ")\n",
    "np.round(T, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:30.015909Z",
     "iopub.status.busy": "2022-12-15T13:10:30.015400Z",
     "iopub.status.idle": "2022-12-15T13:10:30.022217Z",
     "shell.execute_reply": "2022-12-15T13:10:30.021163Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:30.027438Z",
     "iopub.status.busy": "2022-12-15T13:10:30.026892Z",
     "iopub.status.idle": "2022-12-15T13:10:30.306849Z",
     "shell.execute_reply": "2022-12-15T13:10:30.305777Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc2f47248d0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNUUlEQVR4nO3deVxU9f7H8dcAsigDCqJC4A6GlvtGllLatc2bdbNNU0ztd7tqejVTs9TUm9TNtEW93haX1OqWWZZpi4Wa5YZhWriGaaZhLiCorOf3x5FBcgNlPDPD+/l4zEPmzGHOm3F0PnxXm2EYBiIiIiIW8LI6gIiIiFRcKkRERETEMipERERExDIqRERERMQyKkRERETEMipERERExDIqRERERMQyKkRERETEMj5WB7iQwsJCfvvtN+x2Ozabzeo4IiIiUgqGYXD8+HEiIiLw8rpwm4dLFyK//fYbUVFRVscQERGRS7Bv3z4iIyMveI5LFyJ2ux0wf5CgoCCL04iIiEhpZGZmEhUV5fgcvxCXLkSKumOCgoJUiIiIiLiZ0gyr0GBVERERsYwKEREREbGMChERERGxjEuPESkNwzDIz8+noKDA6igiTuHt7Y2Pj4+msIuIR3LrQiQ3N5cDBw5w4sQJq6OIOFXlypUJDw/H19fX6igiIuXKbQuRwsJC0tLS8Pb2JiIiAl9fX/3GKB7HMAxyc3M5dOgQaWlpREdHX3RxIBERd+K2hUhubi6FhYVERUVRuXJlq+OIOE1AQACVKlXil19+ITc3F39/f6sjiYiUG7f/1Uq/HUpFoPe5iHgq/e8mIiIilrlihUhiYiI2m42hQ4deqUuKiIiIi7sihciGDRuYNWsWTZs2vRKXkzPEx8eXqfibM2cOVatWdVoeERGRMzm9EMnKyqJnz5689tprVKtWzdmXEyfbs2cPNpuNlJQUq6OIiIgHcHohMnDgQG6//Xa6dOly0XNzcnLIzMwscRMREff29dcwYwbk51udRFyRUwuRd955h02bNjF58uRSnT958mSCg4Mdt6ioqDJdzzAMTuTmW3IzDKPUOePj4xk8eDBDhw6lWrVq1KxZk9dee43s7Gz69u2L3W6nYcOGLFu2rMT3rVy5krZt2+Ln50d4eDijRo0i/4x/2dnZ2fTu3ZvAwEDCw8OZMmXKWdfOycnh8ccf56qrrqJKlSq0a9eOpKSkUmevV68eAC1atMBmsxEfHw+Y67pMmDCByMhI/Pz8aN68OcuXL3d8X1FLyjvvvMN1112Hv78/11xzDStXriz1tUXEPe3aBQMHwogRVicRV+S0dUT27dvHkCFD+OKLL0q97sHo0aMZNmyY435mZmaZipGTeQU0HvtZmbOWh58mdKWyb+lfzrlz5/LEE0+wfv163n33XR599FEWL17MXXfdxZNPPsnUqVN56KGH2Lt3L5UrV2b//v3cdtttJCQkMG/ePLZt28aAAQPw9/dn/PjxAIwYMYKVK1fy0UcfUaNGDZ588kk2bdpE8+bNHdcdNGgQP/30E++88w4REREsXryYW265hS1bthAdHX3R3OvXr6dt27Z8+eWXNGnSxLHS50svvcSUKVOYNWsWLVq04M033+Svf/0rP/74Y4nnHTFiBNOmTaNx48a8+OKLdOvWjbS0NEJDQ0v92omIe2nWzPyzWjUwDNDak3Imm1GWX+XL4MMPP+Suu+7C29vbcaygoACbzYaXlxc5OTklHjuXzMxMgoODycjIICgoqMRjp06dIi0tjXr16jkKnRO5+W5RiMTHx1NQUMDq1asB83UJDg7m7rvvZt68eQAcPHiQ8PBwvvvuO9q3b8+YMWNYtGgRqampjhVkZ8yYwciRI8nIyODEiROEhoYyf/58evToAcCRI0eIjIzkkUceYdq0aezdu5f69euzd+9eIiIiHHm6dOlC27ZtefbZZ5kzZw5Dhw7l2LFj58y+Z88e6tWrx/fff1+iwLnqqqsYOHAgTz75pONY27ZtadOmDdOnT3d8X2JiIiNHjgQgPz+fevXqMXjwYJ544onSvdAV1Lne7yKuau9emDQJXn4Zit6u27bB1Vdbm0uunAt9fv+Z01pEOnfuzJYtW0oc69u3L1dffTUjR468aBFyKQIqefPThK7l/rylvXZZnDmDyNvbm9DQUK699lrHsZo1awKQnp4OQGpqKnFxcSWWse/QoQNZWVn8+uuvHD16lNzcXNq1a+d4PCQkhEaNGjnub9myhYKCAmJiYkpkycnJuawWiczMTH777Tc6dOhQ4niHDh3YvHlziWNxcXGOr318fGjdujWpqamXfG0RcS0FBdC1q1l4BATASy+Zx88sQgoLISPDbCERcVohYrfbueaaa0ocq1KlCqGhoWcdLy82m61M3SNWqlSpUon7NputxLGigqOwsLDcrpmVlYW3tzfJyclnFYKBgYHldh0Rqbi8vc3i48kn4YyedofMTOjVC9LTISmpuMVEKi6trOomYmNj+e6770oMil2zZg12u53IyEgaNGhApUqVWLdunePxo0ePsmPHDsf9Fi1aUFBQQHp6Og0bNixxq1WrVqlyFI0JKSgocBwLCgoiIiKCNWvWlDh3zZo1NG7cuMSxtWvXOr7Oz88nOTmZ2NjYUl1bRFxTfj7s21d8/y9/gfXroU6ds8/94w/45htISYGNG69YRHFhV7T5oCyzM6Skf/zjH0ybNo3BgwczaNAgtm/fzrhx4xg2bBheXl4EBgbSr18/RowYQWhoKDVq1GDMmDEl9iiJiYmhZ8+e9O7dmylTptCiRQsOHTrEihUraNq0KbfffvtFc9SoUYOAgACWL19OZGQk/v7+BAcHM2LECMaNG0eDBg1o3rw5s2fPJiUlhQULFpT4/unTpxMdHU1sbCxTp07l6NGjPPzww+X+eonIlXH0KNx7L+zebRYf1aubx8+3PVL9+rBoEQQGQps2Vy6nuC736McQrrrqKj799FNGjBhBs2bNCAkJoV+/fjz11FOOc/7973+TlZVFt27dsNvtDB8+nIyMjBLPM3v2bCZNmsTw4cPZv38/1atXp3379txxxx2lyuHj48PLL7/MhAkTGDt2LDfccANJSUk89thjZGRkMHz4cNLT02ncuDFLliw5ayZOYmIiiYmJpKSk0LBhQ5YsWUL1ov+5RMTtFBbCzz/D77/Dli1w440X/57SnCMVh9NmzZSHss6aEdd1vtk2Ujp6v4sr+/FHs3umaJpuWezdCw8/DLNmQYMG5Z9NrFGWWTMaIyIiIqVmGPDcc/DZGSslNGlyaUUIwKBBsGIFDBhQPvnE/ahrRkRESu2NN2DUKAgOhu3b4fRKA5fsP/8xp/zOmFE++cT9qBCRK6Ju3bplWgZfRFzTQw/B/Plw332XX4QARETA0qWX/zzivlSIiIjIBaWlweltpvDzg6++Ov+smMuVnAwHDkApx8+LB9AYEREROa958yA21uxCKeKsImTjRrj+erj/fti61TnXENejQkRERM5r/37IyTEHlDq7d7V5c7MQiY+HMm6+Lm5MXTMiInJeo0aZi5D16OH8XXN9fOCDD6ByZXOpeKkY1CIiIiIO27fDwIHmTBYwi4/77nNed8yf2e0li5A/7Z0qHkiFiAXi4+MZOnSo1THKzZ49e7DZbKSkpFgdRUQuw8mT5qqnM2bA5MlWp4F//QuaNjWnDIvnUiEiJZw6dYqEhASuvfZafHx86N69+xW7tqcVaCLuJiAApk41x2m4wgJjRZuP//STtTnEuTRGpILKzc117KR7poKCAgICAnjsscdYtGiRBclE5Eo6dcrcuC483Lx/333meJAr1RVzIWPGQLt25m6+4rlc4K0mb731Fq1bt8Zut1OrVi0efPBB0tPTATAMg4YNG/LCCy+U+J6UlBRsNhu7du0C4NixY/Tv35+wsDCCgoK46aab2Lx5s+P88ePH07x5c15//fUL7ldSpUoVZs6cyYABA6hVq1aZfo5t27Zx3XXX4e/vzzXXXMPKlStLPL5y5Uratm2Ln58f4eHhjBo1ivz8fAASEhJYuXIlL730EjabDZvNxp49e8p0fREpmwMHzBkqt9wCWVnFx12hCAEzx5lFiGFAXp51ecQ5XOTtVk4MA3Kzrbldxry2vLw8Jk6cyObNm/nwww/Zs2cPCQkJANhsNh5++GFmz55d4ntmz55Nx44dadiwIQA9evQgPT2dZcuWkZycTMuWLencuTNHjhxxfM+uXbtYtGgRH3zwgVPGc4wYMYLhw4fz/fffExcXR7du3Th8+DAA+/fv57bbbqNNmzZs3ryZmTNn8sYbbzBp0iQAXnrpJeLi4hgwYAAHDhzgwIEDRGn+nohTFRTAnj2wb585SNWVnToFffuaG+RpkWbP4lldM3kn4NkIa6795G/gW+WSvvXhhx92fF2/fn1efvll2rRpQ1ZWFoGBgSQkJDB27FjWr19P27ZtycvLY+HChY5Wkm+++Yb169eTnp6On58fAC+88AIffvgh77//Po888ghgdsfMmzePsLCwy/xhz23QoEH87W9/A2DmzJksX76cN954gyeeeIIZM2YQFRXFq6++is1m4+qrr+a3335j5MiRjB07luDgYHx9falcuXKZW2JE5NJERsLHH0O1anD6dxqXlZJiLi1vGDB0KLRqZXUiKS+e1SLippKTk+nWrRu1a9fGbrfTqVMnAPbu3QtAREQEt99+O2+++SYAH3/8MTk5OfTo0QOAzZs3k5WVRWhoKIGBgY5bWloau3fvdlynTp06TitCAOLi4hxf+/j40Lp1a1JTUwFITU0lLi4O2xkLEXTo0IGsrCx+/fVXp2USkWIFBea4i3Xrio+1aeP6RQhA+/bm6q6ffaYixNN4VotIpcpmy4RV174E2dnZdO3ala5du7JgwQLCwsLYu3cvXbt2JTc313Fe//79eeihh5g6dSqzZ8/mvvvuo3Jl85pZWVmEh4eTlJR01vNXrVrV8XWVKpfWYiMinuHf/4Znn4XZs82uGLvd6kRl07+/1QnEGTyrELHZLrl7xCrbtm3j8OHDJCYmOsZEbNy48azzbrvtNsdA0uXLl7Nq1SrHYy1btuTgwYP4+PhQt27dKxX9LGvXrqVjx44A5Ofnk5yczKBBgwCIjY1l0aJFGIbhaBVZs2YNdrudyMhIAHx9fSkoWkVJRMrdwIHmyqX//Kf7FSF/9scfMHo0TJkCQUFWp5HLoa4Zi9WuXRtfX19eeeUVfv75Z5YsWcLEiRPPOs/b25uEhARGjx5NdHR0iW6QLl26EBcXR/fu3fn888/Zs2cP3377LWPGjDlnUXMxP/30EykpKRw5coSMjAxSUlJKNbh1+vTpLF68mG3btjFw4ECOHj3qGP/yj3/8g3379jF48GC2bdvGRx99xLhx4xg2bBhep4fo161bl3Xr1rFnzx7++OMPCosWERCRS5aWVvy13Q5r18IDD1iXpzwYBtx1F7z+Ovzf/1mdRi6b4cIyMjIMwMjIyDjrsZMnTxo//fSTcfLkSQuSXZ5OnToZQ4YMcdxfuHChUbduXcPPz8+Ii4szlixZYgDG999/X+L7du/ebQDG888/f9ZzZmZmGoMHDzYiIiKMSpUqGVFRUUbPnj2NvXv3GoZhGOPGjTOaNWtWqnx16tQxgLNu55OWlmYAxsKFC422bdsavr6+RuPGjY2vvvqqxHlJSUlGmzZtDF9fX6NWrVrGyJEjjby8PMfj27dvN9q3b28EBAQYgJGWllaqvBWBO7/fxTovvWQYPj6G8d57Vicpf+vXG0aTJoaxdavVSeRcLvT5/Wc2w3DdiVCZmZkEBweTkZFB0J/a3k6dOkVaWtoF18TwNKtXr6Zz587s27ePmjVrWh1HrqCK+H6Xyzd8OLz4Ijz6qLlsu6cpLHSdNU+kpAt9fv+ZZ40R8VA5OTkcOnSI8ePH06NHDxUhIlIqzz1nzja55x6rkzjHmUXIjh3mWiNNm1qXRy6Nakk38Pbbb1OnTh2OHTvG888/b3UcEXFRmzbB448XL/jl42Mu137GrHmPtHYttG0Lf/0rHDpkdRopK7WIuIGEhATHSqsiIudy9Ki5c25mprkuyN//bnWiK6dRI6hRA8LCijfKE/ehQkRExANUqwaTJ5srpbr7rJiyqlYNvvwSataE04tLixtR14yIiJs6ftxcT6PIo4/C0qUQHGxdJqvUrl2yCFEXjftQISIi4oZ+/hni4sz1NIoWYbbZNIsEzKXg69aF1autTiKlobesiIgbysuDX3+F3bvhl1+sTuM6DAO+/hpOnID//c/qNFIaGiMiIuKGGjUyx4PUrw9XXWV1Gtdhs8Gbb0KXLtqbxl2oRURExA3k5pp7xPz0U/GxG25QEXIuVarAgAGeP23ZU6gQsUB8fDxDhw61Oka5SUpKwmazcezYMaujiHisJ5+EadPMMSF5eVancR8FBTBihLm4m7gmFSJSQlJSEnfeeSfh4eFUqVKF5s2bs2DBgity7bp16zJt2rQrci0RdzNyJDRrBlOnQqVKVqdxH0uXwgsvmIXczp1Wp5Fz0RiRCio3NxdfX9+zjn/77bc0bdqUkSNHUrNmTT755BN69+5NcHAwd9xxhwVJRSqutDSoV8/8OizMXDlVs2LK5q9/haFD4brrIDra6jRyLnpLu4C33nqL1q1bY7fbqVWrFg8++CDp6ekAGIZBw4YNeeGFF0p8T0pKCjabjV27dgFw7Ngx+vfvT1hYGEFBQdx0001s3rzZcf748eNp3rw5r7/++gU3TnvyySeZOHEi1113HQ0aNGDIkCHccsstfPDBBxf9OdasWUPTpk3x9/enffv2bN26tcTjixYtokmTJvj5+VG3bl2mTJnieCw+Pp5ffvmFf/7zn9hsNmzq3JUKzDDgmWcgJgZWrCg+riLk0kydai51L67Jo97WhmFwIu+EJbfL2cQ4Ly+PiRMnsnnzZj788EP27NnjWNLdZrPx8MMPM3v27BLfM3v2bDp27EjDhg0B6NGjB+np6Sxbtozk5GRatmxJ586dOXLkiON7du3axaJFi/jggw9ISUkpdb6MjAxCQkIuet6IESOYMmUKGzZsICwsjG7dupF3ujM7OTmZe++9l/vvv58tW7Ywfvx4nn76aebMmQPABx98QGRkJBMmTODAgQMcOHCg1PlEPNHPP0N+PqxaZXUSz5KVBYmJ5tgRcQ0e1TVzMv8k7Ra2s+Ta6x5cR+VKlS/pex9++GHH1/Xr1+fll1+mTZs2ZGVlERgYSEJCAmPHjmX9+vW0bduWvLw8Fi5c6Ggl+eabb1i/fj3p6en4nV5a8IUXXuDDDz/k/fff55FHHgHM7ph58+YRFhZW6mz/+9//2LBhA7NmzbrouePGjePmm28GYO7cuURGRrJ48WLuvfdeXnzxRTp37szTTz8NQExMDD/99BP//ve/SUhIICQkBG9vb0erkEhFZrPBrFnmwNTu3a1O4zkKC6FzZ1i/3tyT59lnrU4k4GEtIu4qOTmZbt26Ubt2bex2O506dQJg7969AERERHD77bfz5ptvAvDxxx+Tk5NDj9NtjZs3byYrK4vQ0FACAwMdt7S0NHbv3u24Tp06dcpUhHz99df07duX1157jSZNmlz0/Li4OMfXISEhNGrUiNTUVABSU1Pp0KFDifM7dOjAzp07KdCvJiKsXg3jxhXf9/dXEVLevLzM8SLh4dCtm9VppIhHtYgE+ASw7sF1ll37UmRnZ9O1a1e6du3KggULCAsLY+/evXTt2pXconWbgf79+/PQQw8xdepUZs+ezX333UflymYLTFZWFuHh4SQlJZ31/FWrVnV8XaVKlVLnWrlyJd26dWPq1Kn07t37kn42ESmdffvMBbhyc+Haa+Gee6xO5LkeeADuuAPsdquTSBGPKkRsNtsld49YZdu2bRw+fJjExESioqIA2Lhx41nn3XbbbVSpUoWZM2eyfPlyVp3RcdyyZUsOHjyIj48PdevWvexMSUlJ3HHHHTz33HOObp3SWLt2LbVr1wbg6NGj7Nixg9jYWABiY2NZs2ZNifPXrFlDTEwM3t7eAPj6+qp1RCqkqCh46in48Ue47Tar03i+M4uQ/fvBx8fcuVesoa4Zi9WuXRtfX19eeeUVfv75Z5YsWcLEiRPPOs/b25uEhARGjx5NdHR0iW6QLl26EBcXR/fu3fn888/Zs2cP3377LWPGjDlnUXMhX3/9NbfffjuPPfYYf/vb3zh48CAHDx4sMej1fCZMmMCKFSvYunUrCQkJVK9ene6n25aHDx/OihUrmDhxIjt27GDu3Lm8+uqrPP74447vr1u3LqtWrWL//v38ceaWoiIe6PBhc5xCkaeegrffhsru9buUW9uwAVq3hr/9rXjjQLnyVIhYLCwsjDlz5vDee+/RuHFjEhMTz5qqW6Rfv37k5ubSt2/fEsdtNhuffvopHTt2pG/fvsTExHD//ffzyy+/ULOMZf7cuXM5ceIEkydPJjw83HG7++67L/q9iYmJDBkyhFatWnHw4EE+/vhjx1olLVu25H//+x/vvPMO11xzDWPHjmXChAmO2UFgFjJ79uyhQYMGZRrLIuJutm6FNm2gZ8/i2Rs2m5Ykv9KCgszN8TIyzMJQrGEzLmfeqZNlZmYSHBxMRkYGQUFBJR47deoUaWlpF1wTw9OsXr2azp07s2/fvjIXGOLeKuL73ZNt2AAdO0JEhDlINSLC6kQV1/r10LgxBAZancSzXOjz+888aoyIp8rJyeHQoUOMHz+eHj16qAgRcXNt2sCSJdCyJYSGWp2mYmvbtuT93Fw4x6LT4kTqmnEDb7/9NnXq1OHYsWM8//zzVscRkTI6cQIefRROz8gH4OabVYS4mvfeM1ez/eUXq5NULCpE3EBCQgIFBQUkJydzlfb8FnE7jz4K//mPucy463aGV2wFBeYOvb/8Yi4JL1eOumZERJxswgRITjZ3gdWAVNfk7Q2LF8Mbb8DpBaDlClEhIiLiBGfunFunDvzwgzatc3VRUTB+vNUpKh79sxARKUf5+eYy4o0bw5nL+KgIcS+GYXbVfPyx1Uk8n/5piIiUI5sNdu2CU6fg22+tTiOXau5cGDXKXOtFm4E7l7pmRETKkbc3LFwI33yj5drdWc+esGCBuepqeLjVaTybChERkcu0dCmkpkLRjgVBQSpC3F2lSvDZZ+pSuxL0Enu4+Ph4hg4dWurz58yZU2LH3j/bs2cPNpuNlJSUy84m4gm2bDG3lB8xAr76yuo0Up7OLEJyc2HRIuuyeDIVImKJshZIIq7q2mth0CD4+9/h+uutTiPOkJcHnTvDPffA/PlWp/E86poRESmjAwegWjUo2vZn2jQ14XuySpWgUyez9Uur4ZY//dOxQHx8PIMHD2bo0KFUq1aNmjVr8tprr5GdnU3fvn2x2+00bNiQZcuWlfi+lStX0rZtW/z8/AgPD2fUqFHk5+c7Hs/OzqZ3794EBgYSHh7OlClTzrp2Tk4Ojz/+OFdddRVVqlShXbt2JCUllfln2LZtG9dddx3+/v5cc801rFy5stRZExISWLlyJS+99BI2mw2bzcaePXvKnEHECuvWQatWZgtI0SqpKkI834QJ5lowt95qdRLP41H/fAzDoPDECUtuZd3EeO7cuVSvXp3169czePBgHn30UXr06MF1113Hpk2b+Mtf/sJDDz3EiRMnANi/fz+33XYbbdq0YfPmzcycOZM33niDSZMmOZ5zxIgRrFy5ko8++ojPP/+cpKQkNm3aVOK6gwYN4rvvvuOdd97hhx9+oEePHtxyyy3s3LmzTPlHjBjB8OHD+f7774mLi6Nbt24cPr2P9sWyvvTSS8TFxTFgwAAOHDjAgQMHiIqKKtP1Raxy/Dikp5srpWZkWJ1GrhQvL6hdu/j+0aOQk2NdHk9iM8r6CVoGM2fOZObMmY7fdps0acLYsWO5tZQl5YW2ET7XtuiFJ06wvWWrcv0ZSqvRpmS8Klcu1bnx8fEUFBSwevVqAAoKCggODubuu+9m3rx5ABw8eJDw8HC+++472rdvz5gxY1i0aBGpqanYTq8RPWPGDEaOHElGRgYnTpwgNDSU+fPn06NHDwCOHDlCZGQkjzzyCNOmTWPv3r3Ur1+fvXv3EnHGvuNdunShbdu2PPvss8yZM4ehQ4dy7Nixc2bfs2cP9erVIzExkZEjRwKQn59PvXr1GDx4ME888cRFs3p5eREfH0/z5s2ZNm1amV/riuhc73exziefmE31drvVScQKP/4Id94J8fHw2mtatv9cLvT5/WdOHSMSGRlJYmIi0dHRGIbB3LlzufPOO/n+++9p0qSJMy/t8po2ber42tvbm9DQUK699lrHsZo1awKQnp4OQGpqKnFxcY4PdoAOHTqQlZXFr7/+ytGjR8nNzaVdu3aOx0NCQmjUqJHj/pYtWygoKCAmJqZElpycHELL2PEZFxfn+NrHx4fWrVuTmppaqqy1z/y1QsTFHTtmzoiZPBmqVzeP3XGHpZHEYr/+ai7hX1AAhw8Xvy/k0ji1EOnWrVuJ+//617+YOXMma9eudUohYgsIoNGm5HJ/3tJeuywqVapU8vttthLHij7ECwsLLz/caVlZWXh7e5OcnIy3t3eJxwIDA8vtOiKepGdP+PRT2L/f/FOka1d45x248UYVIeXhis2aKSgo4L333iM7O7vEb9NnysnJIeeMTrfMzMwyXcNms2ErZfeIu4mNjWXRokUYhuEoUtasWYPdbicyMpKQkBAqVarEunXrHC0OR48eZceOHXTq1AmAFi1aUFBQQHp6OjfccMNl5Vm7di0dO3YEzK6Z5ORkBg0aVKqsAL6+vhQUFFxWBpEr4fnnza3h//Uvq5NcWTt37uT48eNWx7ggu91OdHS0Jdc+3QPuYBjqorlUTi9EtmzZQlxcHKdOnSIwMJDFixfTuHHjc547efJknnnmGWdHckv/+Mc/mDZtGoMHD2bQoEFs376dcePGMWzYMLy8vAgMDKRfv36MGDGC0NBQatSowZgxY/A6Yzh/TEwMPXv2pHfv3kyZMoUWLVpw6NAhVqxYQdOmTbn99ttLnWf69OlER0cTGxvL1KlTOXr0KA8//HCpsgLUrVuXdevWsWfPHgIDAwkJCSmRVcQq+fnmKqlFPaVNmlS8nXN37tx5Vheuq9qxY4dlxUiRL780C9WPPwY1Lped0wuRRo0akZKSQkZGBu+//z59+vRh5cqV5yxGRo8ezbBhwxz3MzMzNZvitKuuuopPP/2UESNG0KxZM0JCQujXrx9PPfWU45x///vfZGVl0a1bN+x2O8OHDyfjT8P6Z8+ezaRJkxg+fDj79++nevXqtG/fnjvK2OmdmJhIYmIiKSkpNGzYkCVLllD9dBtlabI+/vjj9OnTh8aNG3Py5EnS0tKoW7fupb9AIuXg6FFzb5HkZHPn3KLPt4pUhACOlpD58+cTGxtrcZpzS01NpVevXpa32pw6BQkJZtddYiKcMZFRSsmps2bOpUuXLjRo0IBZs2Zd9NyyzpoR8VR6v18ZeXlw003w/ffw3nsVd82ITZs20apVK5KTk2nZsqXVcc7JlTJ+9x28+Sa88krxIncVncvMmjmXwsLCEuNARESsVPSrmM1mrqD5v/+ZLSPn6UEWOUtcnHmTS+PUQmT06NHceuut1K5dm+PHj7Nw4UKSkpL47LPPnHlZEZFSyc42V0ht2xYGDzaPhYdr23e5PG++ab6nrrnG6iTuwamFSHp6Or179+bAgQMEBwfTtGlTPvvsM26++WZnXlZEpFTee8/cxGzRIrjvPqhRw+pE4u5mzICBA6F+fdi0CYKDrU7k+pxaiLzxxhvOfHoRkcvSp485MLVHDxUhUj7uvRemTIG+feEiQyPkNO2+KyIVRkEBvPEGPPww+PiY40JeecXqVOJJqlc3p3tXqWJ1EvdRwSaliUhFdu+98H//B08+aXUS8WRnFiEFBWarm5yfChERqTAeeMBccMpFZ6SKhzlxwtwcr0MHWLfO6jSuS4WIiHgswzCn4ha55x74+We4/37rMknF4e9vTgm32cyN8uTcNEZERDzSiRPm1NyNG2H9+uKlt8PCrM0lFYeXF8ybB7t3Q/PmVqdxXWoRcVF169Zl2rRpVse4oubMmUPVqlWtjnGWEydO8Le//Y2goCBsNhvHjh2zOpKUwvHjsGIFbN8OX39tdRqpqOz2kkXIqVPFi+iJSS0iLmrDhg1U8YBh1wkJCRw7dowPP/zwoufed9993Hbbbc4PVUZz585l9erVfPvtt1SvXp1gLQzgFmrWhPffh5wciI+3Oo2IuYtz9+7mtPGhQ61O4zpUiLiosArWfpyXl0dAQAABAQFWR3HIzc3F19eX3bt3ExsbyzVaJtGlFRTA+PHwl7/ADTeYx7TstriSpUshJQV+/x3699dOvUXUNWOB+Ph4Bg0axKBBgwgODqZ69eo8/fTTnLn/4J+7Zo4dO0b//v0JCwsjKCiIm266ic2bN5c432aznXUrsmXLFm666SYCAgIIDQ3lkUceISsry/F4QkIC3bt359lnn6VmzZpUrVqVCRMmkJ+fz4gRIwgJCSEyMpLZs2eX+Fn27dvHvffeS9WqVQkJCeHOO+9kz549AIwfP565c+fy0UcfOfIkJSWxZ88ebDYb7777Lp06dcLf358FCxacs2vm448/pk2bNvj7+1O9enXuuuuu876umzdv5sYbb8RutxMUFESrVq3YuHGjI0vzP3XSTps2rcSOv0Wvwb/+9S8iIiJo1KgR8fHxTJkyhVWrVmGz2Yg//av1W2+9RevWrbHb7dSqVYsHH3yQ9PT0Es//448/cscddxAUFITdbueGG25g9+7djsdff/11YmNj8ff35+qrr2bGjBnn/dnk4qZMMXc+vfdeyMy0Oo3I2R59FCZOhLVrVYScySNbRLKzzT8rVzZHKwPk5po7a/r4gJ/f2ecGBBRv9Z2XZ57v7V1yJ8XznVupUtkzzp07l379+rF+/Xo2btzII488Qu3atRkwYMA5z+/RowcBAQEsW7aM4OBgZs2aRefOndmxYwchISFs2LCBgoICAAoKCrjnnnuodDpYdnY2Xbt2JS4ujg0bNpCenk7//v0ZNGgQc+bMcVzjq6++IjIyklWrVrFmzRr69evHt99+S8eOHVm3bh3vvvsu//d//8fNN99MZGQkeXl5juddvXo1Pj4+TJo0iVtuuYUffviBxx9/nNTUVDIzMx0FTEhICL/99hsAo0aNYsqUKbRo0QJ/f/+z9iBaunQpd911F2PGjGHevHnk5uby6aefnvc17dmzJy1atGDmzJl4e3uTkpLieA1Ka8WKFQQFBfHFF18AEB4ezqhRo9i6dSsffPABvr6+gNmCM3HiRBo1akR6ejrDhg0jISHBkW///v107NiR+Ph4vvrqK4KCglizZg35+fkALFiwgLFjx/Lqq6/SokULvv/+ewYMGECVKlXo06dPmTKLaeBA+OADGDJEK1qKa7LZ4KmnrE7hggwXlpGRYQBGRkbGWY+dPHnS+Omnn4yTJ0+e9Zg5FMgw0tOLj02aZB7r37/kuZUrm8fT0oqPTZ1qHnvwwZLnVq9uHt+6tfjYf/9b9p+rU6dORmxsrFFYWOg4NnLkSCM2NtZxv06dOsbUqVMNwzCM1atXG0FBQcapU6dKPE+DBg2MWbNmnfX8jz32mFGnTh0j/fQL8N///teoVq2akZWV5Thn6dKlhpeXl3Hw4EHDMAyjT58+Rp06dYyCggLHOY0aNTJuuOEGx/38/HyjSpUqxttvv20YhmG89dZbRqNGjUr8HDk5OUZAQIDx2WefOZ73zjvvLJEvLS3NAIxp06aVOD579mwjODjYcT8uLs7o2bPnWT/f+djtdmPOnDnnfGzcuHFGs2bNShybOnWqUadOHcf9Pn36GDVr1jRycnJKnDdkyBCjU6dOF7z2hg0bDMA4fvy4YRiGMXr0aKNevXpGbm7uOc9v0KCBsXDhwhLHJk6caMTFxZ3z/Au93yuqwkLDWLOm5LEz3r5yGZKTkw3ASE5OtjrKeblDxotZv94w/v53z3zfXujz+8/UNWOR9u3bl+g6iYuLY+fOnY5WjTNt3ryZrKwsQkNDCQwMdNzS0tJKNPUD/Pe//+WNN95gyZIljnEmqampNGvWrMTg1w4dOlBYWMj27dsdx5o0aYKXV/FbombNmlx77bWO+97e3oSGhjq6IDZv3syuXbuw2+2OTCEhIZw6deqsXOfSunXrCz6ekpJC586dL/o8RYYNG0b//v3p0qULiYmJpcrwZ9dee62j1eNCkpOT6datG7Vr18Zut9OpUycA9u7d68h+ww03nLNFJjs7m927d9OvX78Sf5+TJk26pMwVkWFAv37mQlHvvlt83Ev/o4mbyMiAm2+G//xH2wx4ZNdM0dCHypWLj40YYY5S9vnTT1zUrX/mGMmBA2HAALNr5kynhz6UODchoRwCX0RWVhbh4eEkJSWd9diZYyq+/vprBg8ezNtvv03Tpk3LfJ0/f2jabLZzHissLHTkatWqFQsWLDjruUoz2PZis4LKOnB1/PjxPPjggyxdupRly5Yxbtw43nnnHe666y68vLxKjMEBs3ulrJmguKura9euLFiwgLCwMPbu3UvXrl3Jzc29aPaisTmvvfYa7dq1K/GY95/fdHJONps5K8bLCw4csDqNSNkFB8NLL5k7P/fta3Uaa3lkIXKuzxJfX/NWmnMrVTr3uI/znXsp1v1pvd+1a9cSHR19zg+ili1bcvDgQXx8fEoMrjzTrl27uOeee3jyySe5++67SzwWGxvLnDlzyM7OdnzQrlmzBi8vLxo1anRpP8DpXO+++y41atQg6Dyd8r6+vuds5SmNpk2bsmLFCvqW4V9pTEwMMTEx/POf/+SBBx5g9uzZ3HXXXYSFhXHw4EEMw3C0RKWkpFxSrm3btnH48GESExOJiooCcAyKPTP73LlzycvLO6uYq1mzJhEREfz888/07NnzkjJUVIWFxa0eEyfC3XdDmzbWZhK5VH36QO/exWMZKyo1ZFpk7969DBs2jO3bt/P222/zyiuvMGTIkHOe26VLF+Li4ujevTuff/45e/bs4dtvv2XMmDFs3LiRkydP0q1bN1q0aMEjjzzCwYMHHTcwB3H6+/vTp08ftm7d6mg5eeihh6hZs+Yl/ww9e/akevXq3HnnnaxevZq0tDSSkpJ47LHH+PX0esZ169blhx9+YPv27fzxxx/nbIU4n3HjxvH2228zbtw4UlNT2bJlC88999w5zz158iSDBg0iKSmJX375hTVr1rBhwwZiY2MBc6bSoUOHeP7559m9ezfTp09n2bJll/Rz165dG19fX1555RV+/vlnlixZwsSJE0ucM2jQIDIzM7n//vvZuHEjO3fu5K233nJ0hT3zzDNMnjyZl19+mR07drBlyxZmz57Niy++eEmZPF1BAYwdC/fdV7wYlI+PihBxf2cWIZ98An/8YV0Wq6gQsUjv3r05efIkbdu2ZeDAgQwZMoRHHnnknOfabDY+/fRTOnbsSN++fYmJieH+++/nl19+oWbNmvz+++9s27aNFStWEBERQXh4uOMGULlyZT777DOOHDlCmzZtuOeee+jcuTOvvvrqZf0MlStXZtWqVdSuXZu7776b2NhY+vXrx6lTpxwtJAMGDKBRo0a0bt2asLAw1qxZU+rnj4+P57333mPJkiU0b96cm266ifXr15/zXG9vbw4fPkzv3r2JiYnh3nvv5dZbb+WZZ54BzFahGTNmMH36dJo1a8b69et5/PHHL+nnDgsLY86cObz33ns0btyYxMREXnjhhRLnhIaG8tVXX5GVlUWnTp1o1aoVr732mqN1pH///rz++uvMnj2ba6+9lk6dOjFnzhzq1at3SZk8XWoqJCaaC5R9+aXVaUTK3/Tp0K2bOf28DL+veQSb8eeOcxeSmZlJcHAwGRkZZzX9nzp1irS0NOrVq4f/mXNs3UB8fDzNmzevcEu4y6Vz5/d7eXn9dXM6fa9eVifxfJs2baJVq1YkJyfT0kW3KnaHjGXx44/Qvj08/LC5Js6fxzO6mwt9fv+Zm/+oIuKJDAPmzjVnFVx1lXmsf39rM4k4U5Mm8NNPcHrYWYWirhkRcTnPPmvOJOjRw1xcUKQiOLMIMQw4vfajx1MhYoGkpCR1y4hcwP33Q2go/PWv7t9ELVJWp06ZM2patoTT4/49mgoREXEJZ27V06AB7NoFo0ZpkTKpePLzzc3x/vgDyjC+323pn7iIWKqwEJ55BurXhy1bio//af9DkQojMBA++gg+/9ycsu7p3L4QceFJPyLlxpPf54YB335rbiq5ZInVaURcQ716cNNNxfc9+L8A9y1EitZjOHHihMVJRJyv6H1e1t2E3YG3NyxcCG+/DWPGWJ1GxPUcOgSdO8OfNij3GG47DMzb25uqVas6NmCrXLlyiU3kRDyBYRicOHGC9PR0qlat6jF70bzxBhw9CkVryoWGmgNUReRsU6bA11/D7t2wc+e5tytxZ25biADUqlULwFGMiHiqqlWrOt7v7m7VKnNNEC8vuPFGaNXK6kQizrNz506OHz9+Wc/RvbuNrVtr07fv72zdeqqckhWz2+1ER0eX+/OWllsXIjabjfDwcGrUqFGmPUxE3EmlSpU8piUE4IYbzEKkXj1o0cLqNCLOs3PnTmJiYsrt+ZYuLbenOsuOHTssK0bcuhAp4u3t7VH/UYt4mtWroW1b8PMzN/n673+146h4vqKWkPnz5zs24HQ1qamp9OrV67JbbS6HRxQiIuK6XnzRHAvy97/DjBnmMRUhUpHExsZ6xH44zuK2s2ZExD00bmz+WVho3kREzqQWEREpd/n5xUuz33ILfP89NGtmbSYRcU1qERGRcjV7NjRvbk7PLaIiRETOR4WIiJSb7GyYMAF+/LF4PIiIyIWoa0ZEyk2VKrBoEXz6KYwebXUaEXEHKkRE5LJ8/rk5C+bmm837LVuaNxGR0lAhIiKX7NNP4Y47oFo1c0Bq7dpWJxIRd6NCREQu2U03mUu0t2gBNWpYnUZE3JEKEREpk337ICrK/NrfH5KSzLEhIiKXQrNmRKTU5s2DmBh4883iYypCRORyqBARkVLbtw9OnYJPPgHDsDqNiHgCdc2ISKmNHm3umnv//dovRkTKh1pEROS8vvwSHnoICgrM+15e8OCD5p8iIuVBLSIick6HD8Ndd0FWFsTFwT/+YXUiEfFEKkRE5JxCQ81l2pOSoG9fq9OIiKdSISIiDj/8AIGBUL++ef+hh8ybiIizqKdXRABzldT27eFvf4OTJ61OIyIVhQoREQGgaVOzNaRmTcjJsTqNiFQU6poRqcBycsDPz/w6MhK+/dacnuvtbW0uEak41CIiUkF99RU0bAgbNhQfa9hQRYiIXFkqREQqqJkz4ddfYdIkq5OISEWmrhmRCuqNNyA6Gp5+2uokIlKRqRApB4ZhcDKvwOoYIhe0dSusSrLxj0HmJjE+/vDUeDCAE7mWRhMXcyq/EFslP07lF3IiN9/qOOdUkTMGVPLG5kF7LNgMw3W3rsrMzCQ4OJiMjAyCgoKsjnNeJ3LzaTz2M6tjiJxXfqY/v73eCSPPhxr3rCegwSGrI4nIJfppQlcq+5ZPO8KmTZto1aoVycnJtGzZslyeE8r2+a0WEZEKwCfoFIFN95F3OBDf8GNWxxERcVAhUg4CKnnz04SuVscQKeG3/VC1GlSubN7Pe9rcrM7b+yZrg4nLS0lJ4frrr+ebb76hefPmVsc5p4qcMaCSZ01tUyFSDmw2W7k1k4mcT34+HDlirv0RFVV8fONGSEuDZs0gJsY8tm8fxLWFW26BefPAZgN8LYktbsjfxwsjLwd/Hy+X/b9NGT2HU6fvTp48mTZt2mC326lRowbdu3dn+/btzrykiCXOHGmVnQ0pKebtTEuXwksvQWpq8bGffzb3cvnzzraPPgq1a8OCBcXHUlPNVU9btSp57osvwr33mku0FykshD/+gC1bIDPzcn4yERHncmohsnLlSgYOHMjatWv54osvyMvL4y9/+QvZ2dnOvKzIWU6eND+Yz1y6/NAhWL4cVq0qee5rr8GoUSULhk2boEuXs3ehvfVWqFQJFi0qPvb999CiBfToUfLcGTNg6FBYu7b42PHjMH8+LF5c8twjR8xWjcOHi49VqWL++efh5Y0bQ8eOUKtW8bGQEHjuOXOl1ODgs14OERGX4dS2ouXLl5e4P2fOHGrUqEFycjIdO3Z05qWvLMOAvBNWp/AIBQWQus1G/XqGY2zD7t02Vq72pmYNg9tvK54mPfaZSvyy18aYkXnExJifzp997s0/R/jSvFkhC+cVVx3t2gawZasXX3x6ki6dCwFY94033e72p3WrAjZ8e8px7uw3/flurTftWp0itoF5vayjXqxYEcDVjQoht3hHOKPAj/x8H05k5sDp6XlBATZq1fKnRpgBucXP27FDJexVvKgTkQe5ZobIGvBCog/VquH4foAJT9kYMdRG3TqFcHpqbb2roPDU6W6WM6bbPvWEeYPi43Y/eGJoyWMipeWVf5LKlcw/yXXNXxwrdMZKlU//R+AZruj03V27dhEdHc2WLVu45pprzno8JyeHnDN+Zc3MzCQqKsrlp++Smw3PRlidwu0dPlGNzvOWsPn3pmwYEE/riO8BWPBDD3otfp0u9b/mi4e6O86/ZsZ3/HioMSt6d+OmemazxpLtt3LnO+/Q7qoNrO3fxXFu+9e/ZN3+Niy5/z66NTIL5HW/tuIfn75I47BtvHXX/znOnbb2UX45FkWf5m/TvNYWAA5lh/LFzzdSvfIR/tLgK8e5B47XpNDwIrTyEfx9tFOciFwBT/4GvlXK5akq1PTdwsJChg4dSocOHc5ZhIA5puSZZ565UpHExVQLOEaDkD1s/r0pBYXFvYZ1qu7j9ujltKj1Q4nzh7T7D8dOBVO/2h7Hseui1vFV7zuoXvlwiXO/7nMHvt65eHsVOo61i0wm+ZFOZ+UY2n7mWcfCqhzmwWvfP+t4uP33Uv98IiJytivWIvLoo4+ybNkyvvnmGyIjI895jtu2iKhr5rIYRnErY1YWHDtmIzLSZdfZE/F4KSkpdLj+eta4+NTYCpuxHLtmKkyLyKBBg/jkk09YtWrVeYsQAD8/P/yK9iR3JzZbuTWTVSR5eTBkiDmwsmjjtcAQ8yYi1in0CeBEnvmnq/7fpoyew6mFiGEYDB48mMWLF5OUlES9evWceTlxM19+ae4Aa7NBr15w9dVWJxIRkSvNqYXIwIEDWbhwIR999BF2u52DBw8CEBwcTEBAgDMvLW7g1lth3Dho2VJFiIhIReXUdURmzpxJRkYG8fHxhIeHO27vvvuuMy8rLuyLL+DEGcNpxo+Hv/7VsjgiImIxp3fNiBSZNg3++U+47z54+22PmgYvIiKXyKktIiJnatkSfHwgLMxcuExERES78IhTnTk1t2NHc+8TjQcREZEiahERp/nmG2jTBk6PUQZUhIiISEkqRMQpCgth4EBIToann7Y6jYiIuCoVIuIUXl7w3nvmbrXTplmdRkREXJUKESk3R47A6tXF92Ni4M03i7evFxER+TMVIlIu9u2Dtm3htttg61ar04iIiLvQrBkpF7VqQe3a5rRcLR8jIiKlpUJELllRwWGzQaVK5piQwkJznRAREZHSUNeMXJKTJ82N6l56qfhYaKiKEBERKRsVInJJ3n8fFi6EUaPgt9+sTiMiIu5KXTNySXr1gpQU6NYNIiKsTiMiIu5KhYiU2qefws03m+NBbDaYMsXqRJ5l586dHD9+3OoYF2S324mOjrY6hoh4EBUiUipPPw2TJpmrpb76qtVpPM/OnTuJiYmxOkap7NixQ8WIiJQbFSJSKq1bm6ulhoSU3MhOykdRS8j8+fOJjY21OM25paam0qtXL5dvtRER96JCRM7rzILjzjvhxx+1aZ2zxcbG0rJlS6tjiIhcMZo1I+e0bBnExUFGRvExFSEiIlLeVIjIWXJy4NFHYd06SEy0Oo2IiHgydc3IWfz84IMP4I03YPx4q9OIiIgnUyEiABw4ULxxHUDLluZNRETEmVSICNu2QefOkJsLGzZA3bpWJxIRkYpChYhQu7a5e+7Jk+buuSIiIleKCpEK6sypuZUrw8cfQ2AgBAVZm0tERCoWzZqpgDIyzD1i5s8vPhYRoSJERESuPBUiFdCbb8LSpTB4cMl1QkRERK40dc1UQEOGwPbt0L8/BAdbnUakfGnzwPKTmppqdYTzcuVsUjYqRCoAw4BPPoHbbzf3i/Hygv/8x+pUIuVPmweWD7vdDkCvXr0sTnJxRVnFfakQqQAGDoSZM2HcOC1QJp5NmweWj+joaHbs2OHSGcF9WpbkwlSIVABt2sCsWeqGkYpDmwdePn3Ay5WiQsRDnTk9t29fcwM7bVonIiKuRrNmPNC778JNN5kLlBVRESIiIq5IhYiHycgwx4QkJZnjQkRERFyZumY8THAwvPceLF9uTtMVERFxZSpEPMDPP0N2Nlx7rXn/xhvNm4iIiKtT14yb27jRnBVzxx2Qnm51GhERkbJRi4iba9AAQkOhalXtnCsiIu5HhYgbOnNqbrVq8MUXUKMGBARYm0tERKSs1DXjZg4dMqfmfvpp8bE6dVSEiIiIe1KLiJuZNs2cmrt7N+zcCX5+VicSERG5dCpE3My4cXDgADzxhIoQERFxf+qacXGFhbBkiTkuBMDXF958UyulioiIZ1Ah4sIMAx54AO68E2bMsDqNiIhI+VMh4sJsNnONEF9fqFLF6jQiIiLlT2NEXNCZ03OHD4e//hViYqzNJCIi4gxqEXEx//2vuUpqfr5532ZTESIiIp5LhYgL2b8fhg0z1whZuNDqNCIiIs6nrhkXctVV8NZbkJoKDz1kdRqxQmpqqtURzsuVs4mI+1IhYrGtW81VURs0MO/fdZd5k4rFbrcD0KtXL4uTXFxRVhGR8qBCxEJff20ORK1dG777DoKCrE4kVomOjmbHjh0cP37c6igXZLfbiY6OtjqGiHgQFSIWatTILD5q1SoenCoVlz7gRaQiUiFyhZ05NTciAlatMltEKlWyNpeIiIgVNGvmCtq3D667DtasKT7WoIGKEBERqbhUiFxBkybB2rXw97+be8iIiIhUdOqauYJefBFOnDALEi+VgCIiImoRcab8fFi6tPh+lSrmOiF16liXSURExJWoEHGSvDy49VZzufa337Y6jYiIiGtyaiGyatUqunXrRkREBDabjQ8//NCZl3MplSpBy5ZmK4i/v9VpREREXJNTC5Hs7GyaNWvG9OnTnXkZl2IYxV8/+yykpGilVBERkfNx6mDVW2+9lVtvvdWZl3AJhmFwIu8kL77gw5YfvJg9L7d4rZA6cCLP2nwiFUVOYQ42Xxs5hTmcyDthdZxzcoeMUj6c9Xcd4BOArehDxgO41KyZnJwccnJyHPczMzMtTFN6J/NP0nzK/eweuwijwIf1EY9hb7ba6lgiFVKT/zbhkR8fgR+tTnJ+7pBRyocz/q7XPbiOypUql98TWsylCpHJkyfzzDPPWB3jkvhHpBHeeyJGQSUVISIiIqXkUoXI6NGjGTZsmON+ZmYmUVFRFiYqnQCfANY9uA4eLDoy3Mo4IhXW5s2b6dChA2vWrKFZs2ZWxzknd8go5cNZf9cBPgHl9lyuwKUKET8/P/z8/KyOUWY2m82jmslE3JWflx9GroGfl5/L/pt0h4xSPvR3XTpaR0REREQs49QWkaysLHbt2uW4n5aWRkpKCiEhIdSuXduZlxYRERE34NRCZOPGjdx4442O+0XjP/r06cOcOXOceWkRERFxA04tROLj4zHOXOFLRERE5AwaIyIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIilvGxOoC4v507d3L8+HGrY1yQ3W4nOjra6hhyhaSmplod4bxcOZuIFVSIyGXZuXMnMTExVscolR07dqgY8XB2ux2AXr16WZzk4oqyilR0KkTkshS1hMyfP5/Y2FiL05xbamoqvXr1cvlWG7l80dHR7Nixw+X/rtVCJ1JMhYiUi9jYWFq2bGl1DBF9wIu4GQ1WFREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLqBARERERy6gQEREREcuoEBERERHLXJFCZPr06dStWxd/f3/atWvH+vXrr8RlRURExMU5vRB59913GTZsGOPGjWPTpk00a9aMrl27kp6e7uxLi4iIiIvzcfYFXnzxRQYMGEDfvn0B+M9//sPSpUt58803GTVqlLMvf0UYhoFx8qTVMaxx6hQBNhucOkXhiRNWpzk3d8goIp7HSf/32AICsNls5fZ8VnNqIZKbm0tycjKjR492HPPy8qJLly589913Z52fk5NDTk6O435mZqYz45Ub4+RJtrdsZXUMSwQAyTGN4OF+bLc6zHm4Q0YR8TzO+r+n0aZkbJUrl+MzWsupXTN//PEHBQUF1KxZs8TxmjVrcvDgwbPOnzx5MsHBwY5bVFSUM+OJiIiIxZzeNVMWo0ePZtiwYY77mZmZblGM2AICaLQp2eoYlkhJSeH666/nm2++oXnz5lbHOSd3yCginsdZ//fYAgLK7blcgVMLkerVq+Pt7c3vv/9e4vjvv/9OrVq1zjrfz88PPz8/Z0ZyCpvN5lHNZGXi789JwwB/f7xc9TVwh4wi4nn0f0+pOLVrxtfXl1atWrFixQrHscLCQlasWEFcXJwzLy0iIiJuwOldM8OGDaNPnz60bt2atm3bMm3aNLKzsx2zaERERKTicnohct9993Ho0CHGjh3LwYMHad68OcuXLz9rAKuIiIhUPFdksOqgQYMYNGjQlbiUiIiIuBHtNSMiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZxqb1mREREPE1qaqrVEc7LFbKpEBEREXECu90OQK9evSxOcnFFWa2gQkRERMQJoqOj2bFjB8ePH7c6ygXZ7Xaio6Mtu36FLkR27typN4iIiDiN/v++uApbiOzcuZOYmBirY5TKjh079GYWERGPVGELkaKWkPnz5xMbG2txmnNLTU2lV69eLt9qIyIicqkqbCFSJDY2lpYtW1odQ0REpELSOiIiIiJiGRUiIiIiYhkVIiIiImIZFSIiIiJiGRUiIiIiYhkVIiIiImIZFSIiIiJiGRUiIiIiYhkVIiIiImIZFSIiIiJiGRUiIiIiYhkVIiIiImIZFSIiIiJimQq/+66Uj9TUVKsjnJcrZxMRqehUiMhlsdvtAPTq1cviJBdXlFVERFyHChG5LNHR0ezYsYPjx49bHeWC7HY70dHRVscQEZE/USEil00f8CIicqk0WFVEREQso0JERERELKNCRERERCyjQkREREQso0JERERELKNCRERERCyjQkREREQso0JERERELKNCRERERCyjQkREREQso0JERERELKNCRERERCyjQkREREQso0JERERELKNCRERERCyjQkREREQso0JERERELKNCRERERCyjQkREREQso0JERERELKNCRERERCyjQkREREQso0JERERELKNCRERERCyjQkREREQso0JERERELOO0QuRf//oX1113HZUrV6Zq1arOuoyIiIi4MacVIrm5ufTo0YNHH33UWZcQERERN+fjrCd+5plnAJgzZ46zLiEiIiJuzmmFyKXIyckhJyfHcT8zM9Pp10xNTXX6NS6VK2cTEREpDy5ViEyePNnRkuJsdrsdgF69el2R612OoqwiIiKepkyFyKhRo3juuecueE5qaipXX331JYUZPXo0w4YNc9zPzMwkKirqkp7rYqKjo9mxYwfHjx93yvOXF7vdTnR0tNUxREREnKJMhcjw4cNJSEi44Dn169e/5DB+fn74+fld8veXlT7gRURErFWmQiQsLIywsDBnZREREZEKxmljRPbu3cuRI0fYu3cvBQUFpKSkANCwYUMCAwOddVkRERFxI04rRMaOHcvcuXMd91u0aAHA119/TXx8vLMuKyIiIm7EZhiGYXWI88nMzCQ4OJiMjAyCgoKsjiMiIiKlUJbPb+01IyIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIillEhIiIiIpZRISIiIiKWUSEiIiIilnHaEu/loWjR18zMTIuTiIiISGkVfW6XZvF2ly5Ejh8/DkBUVJTFSURERKSsjh8/TnBw8AXPcem9ZgoLC/ntt9+w2+3YbDar41xQZmYmUVFR7Nu3T/viXAa9juVDr2P50OtYPvQ6lg93eh0Nw+D48eNERETg5XXhUSAu3SLi5eVFZGSk1THKJCgoyOXfIO5Ar2P50OtYPvQ6lg+9juXDXV7Hi7WEFNFgVREREbGMChERERGxjAqRcuLn58e4cePw8/OzOopb0+tYPvQ6lg+9juVDr2P58NTX0aUHq4qIiIhnU4uIiIiIWEaFiIiIiFhGhYiIiIhYRoWIiIiIWEaFSDmZPn06devWxd/fn3bt2rF+/XqrI7mVVatW0a1bNyIiIrDZbHz44YdWR3JLkydPpk2bNtjtdmrUqEH37t3Zvn271bHczsyZM2natKlj4ai4uDiWLVtmdSy3l5iYiM1mY+jQoVZHcSvjx4/HZrOVuF199dVWxyo3KkTKwbvvvsuwYcMYN24cmzZtolmzZnTt2pX09HSro7mN7OxsmjVrxvTp062O4tZWrlzJwIEDWbt2LV988QV5eXn85S9/ITs72+pobiUyMpLExESSk5PZuHEjN910E3feeSc//vij1dHc1oYNG5g1axZNmza1OopbatKkCQcOHHDcvvnmG6sjlRtN3y0H7dq1o02bNrz66quAuUdOVFQUgwcPZtSoURancz82m43FixfTvXt3q6O4vUOHDlGjRg1WrlxJx44drY7j1kJCQvj3v/9Nv379rI7idrKysmjZsiUzZsxg0qRJNG/enGnTplkdy22MHz+eDz/8kJSUFKujOIVaRC5Tbm4uycnJdOnSxXHMy8uLLl268N1331mYTAQyMjIA80NULk1BQQHvvPMO2dnZxMXFWR3HLQ0cOJDbb7+9xP+TUjY7d+4kIiKC+vXr07NnT/bu3Wt1pHLj0pveuYM//viDgoICatasWeJ4zZo12bZtm0WpRMyWuaFDh9KhQweuueYaq+O4nS1bthAXF8epU6cIDAxk8eLFNG7c2OpYbuedd95h06ZNbNiwweoobqtdu3bMmTOHRo0aceDAAZ555hluuOEGtm7dit1utzreZVMhIuKhBg4cyNatWz2qL/lKatSoESkpKWRkZPD+++/Tp08fVq5cqWKkDPbt28eQIUP44osv8Pf3tzqO27r11lsdXzdt2pR27dpRp04d/ve//3lEV6EKkctUvXp1vL29+f3330sc//3336lVq5ZFqaSiGzRoEJ988gmrVq0iMjLS6jhuydfXl4YNGwLQqlUrNmzYwEsvvcSsWbMsTuY+kpOTSU9Pp2XLlo5jBQUFrFq1ildffZWcnBy8vb0tTOieqlatSkxMDLt27bI6SrnQGJHL5OvrS6tWrVixYoXjWGFhIStWrFB/slxxhmEwaNAgFi9ezFdffUW9evWsjuQxCgsLycnJsTqGW+ncuTNbtmwhJSXFcWvdujU9e/YkJSVFRcglysrKYvfu3YSHh1sdpVyoRaQcDBs2jD59+tC6dWvatm3LtGnTyM7Opm/fvlZHcxtZWVklqvu0tDRSUlIICQmhdu3aFiZzLwMHDmThwoV89NFH2O12Dh48CEBwcDABAQEWp3Mfo0eP5tZbb6V27docP36chQsXkpSUxGeffWZ1NLdit9vPGp9UpUoVQkNDNW6pDB5//HG6detGnTp1+O233xg3bhze3t488MADVkcrFypEysF9993HoUOHGDt2LAcPHqR58+YsX778rAGscn4bN27kxhtvdNwfNmwYAH369GHOnDkWpXI/M2fOBCA+Pr7E8dmzZ5OQkHDlA7mp9PR0evfuzYEDBwgODqZp06Z89tln3HzzzVZHkwro119/5YEHHuDw4cOEhYVx/fXXs3btWsLCwqyOVi60joiIiIhYRmNERERExDIqRERERMQyKkRERETEMipERERExDIqRERERMQyKkRERETEMipERERExDIqRERERMQyKkRERETEMipERERExDIqRERERMQyKkRERETEMv8PtnlxO4ssKqUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "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(\n",
    "        [iw - 0.25, iw + 0.25],\n",
    "        scbot[iw],\n",
    "        sctop[iw],\n",
    "        facecolor=\"None\",\n",
    "        edgecolor=\"k\",\n",
    "    )\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": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:30.311182Z",
     "iopub.status.busy": "2022-12-15T13:10:30.310814Z",
     "iopub.status.idle": "2022-12-15T13:10:30.319823Z",
     "shell.execute_reply": "2022-12-15T13:10:30.318956Z"
    }
   },
   "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": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:10:30.323896Z",
     "iopub.status.busy": "2022-12-15T13:10:30.323492Z",
     "iopub.status.idle": "2022-12-15T13:10:30.330230Z",
     "shell.execute_reply": "2022-12-15T13:10:30.329197Z"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    # ignore PermissionError on Windows\n",
    "    temp_dir.cleanup()\n",
    "except:\n",
    "    pass"
   ]
  }
 ],
 "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.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}