effective_groundwater_model_calibration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Effective Groundwater Model Calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import flopy\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conceptual Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/conceptual_model_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/conceptual_model_2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](img/conceptual_model_3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Steady-State model setup \n",
    "\n",
    "We are creating a square model with:\n",
    "\n",
    "* 2 layers\n",
    "* 18 rows\n",
    "* 18 columns\n",
    "* GHB-Boundary (right side) \n",
    "* RIV-Boundary (left side)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set workspace for steady-state model\n",
    "\n",
    "* and make shure that directory exists  \n",
    "* save the starting path\n",
    "* change to the working directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "workspace = os.path.join('steady-state')\n",
    "\n",
    "if not os.path.exists(workspace):\n",
    "    os.makedirs(workspace)\n",
    "    \n",
    "cwdpth = os.getcwd()\n",
    "\n",
    "os.chdir(workspace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the MODFLOW model \n",
    "Store it (in this case in the variable `ml`, but you can call it whatever you want).  \n",
    "The modelname will be the name given to all MODFLOW files (input and output).  \n",
    "The exe_name should be the full path to your MODFLOW executable.  The version is either 'mf2k' for MODFLOW2000 or 'mf2005'for MODFLOW2005."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "name = 'steady-state'\n",
    "\n",
    "ml = flopy.modflow.Modflow(\n",
    "    modelname=name, \n",
    "    exe_name='mf2005',\n",
    "    version='mf2005'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretization\n",
    "\n",
    "We are defining all kind of discretization parameters.  \n",
    "\n",
    "As we have a confining bed between top- and bottom-layer we have to specify  \n",
    "the botm-parameter with three layer values `[50, 40, -10]` and set the   \n",
    "laycbd-parameter with a flag (important: the flag must be 1) on the first layer `[1, 0]`\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "nlay = 2\n",
    "nrow = ncol = 18\n",
    "L = W = 18000\n",
    "delr = L/nrow \n",
    "delc = W/ncol\n",
    "top = 100\n",
    "botm = [50, 40, -10]\n",
    "laycbd = [1, 0]\n",
    "itmuni = 1  # seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dis = flopy.modflow.ModflowDis(\n",
    "    ml, \n",
    "    nlay=nlay, \n",
    "    nrow=nrow, \n",
    "    ncol=ncol, \n",
    "    delr=delr,\n",
    "    delc=delc, \n",
    "    top=top, \n",
    "    botm=botm, \n",
    "    laycbd=laycbd,\n",
    "    itmuni=itmuni\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boundaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BAS-Package \n",
    "\n",
    "Specify boundary conditions and starting heads. \n",
    "\n",
    "We can set one value for all cells.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ibound = 1\n",
    "strt = 200\n",
    "\n",
    "bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=strt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RIV-Package \n",
    "\n",
    "The river boundary is located in the left side of the model and takes place in the upper layer.   \n",
    "\n",
    "The river width is 10m and depth is 10m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.axes._subplots.AxesSubplot at 0x7fdd5bbef9e8>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEICAYAAAB4TcDdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGxpJREFUeJzt3Xu4XVV57/Hvz4QERUIS0BSCkCixNdiqASFWy1G0ScihjX1KbaiagKl4QWu9VKHaEwQ8j1gV4YhYFI4EKSHGS3ksNqSIcrQmIQgCQTDhKjEQIBdULBh4zx/j3TCz3HtnZ42dvdcKv8/z7GfNNcaYc44x91rvnrc9X0UEZmbtetZwd8DMupuDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQeRAZIUkg4Z4nV+UdI/DeU6d0eS3izpqjbnPU3SVwepH9+T9LeDsaxOstsEEUl3S/qNpF9Jul/SVyQ9t1H/FUlnStpT0hZJR/eyjLMlLR3anj+17hMk/aBZFhHvjIgzhrAPkzJYjhyqdQ6FiLg0ImYM9nIljZK0ND97Iem1g72OXUHSSyUtk/SQpOobxXabIJL+LCKeC7wceAVwamuDiPhv4HJgXrNc0gjgeODiIehn1+q2ADME/f0B8Bbg/l28nrb0Mf7fAkuABYOxjt0tiAAQEfcDyyjBpDcXA38p6TmNspmU7fGdHS1f0j6SFkl6UNI9kj4m6VmN+rdL+qmkX0q6VdK0LD9F0h2N8r/I8pcAXwRelXtSW7L8K5LObFnuOkmbJF0h6YBGXUh6p6S1uad1niT10f8jJK2W9IikByR9Nquuzdct2Y9X5R7SD3Mv7WHgtFzG23KMm/Ov2sFZrmy7MZd/s6SXZt3sHPcvJa2X9KE++tezzs9L2irpNkmvb9n+F0rakMs5M/8INOd9qr+te3mS/ljSdbns6yT9caNusqTvZx+XA/v19TmIiMcj4nMR8QPgib7a9THGF0n6rqSHc4/gUkljs+4fJH29pf25ks5pZ/y99Pv2iLgQWLMzfe5TROwWP8DdwBty+kDgZuCcRv1XgDMb738GvKXx/jLgc/0sP4BDcnoR8G/A3sCkXNaCrPsrYD3wSkDAIcDBjboDKMHqr4FfA/tn3QnAD1rW+VSfgaOBh4BpwGjg/wDXtvTv28BY4CDgQWBWH2P5EfDWnH4uMD2nJ+VyRjbangBsA94LjASeDcwB1gEvybKPAf+V7WcC12c/lG16xrgB+JOcHgdM66N/Pet8P7BHbqutwPis/ybwL8BewPOBVcA7+unvU9sWGA9sBt6a9cfn+30b2+azuY2PAn4JfHUAn7/7gNfuoM33gL/N6UOAP831PI8SwD+XdfvnZ2Nsvh8JbAQOa2f8/fTnECCqv3vD/eUfrB9KEPlV/tIDuLrnl9D6hcz3HwOuyukxwKPAK/pZfuRGHwE8Dkxt1L0D+F5OLwPeN8A+3wjMafzy+wsiFwKfatQ9l7JbOqnRv9c06pcAp/Sx3muBjwP7tZRPovcgcm9Lu++QQTPfPyu338GUYPczYDrwrJb57s1tNWYH2+UE4BeAGmWrKF/8CcBjzS8HJRBc009/n9q2uYxVLfU/yjYH5Rdwr0bdv7ILgkgvdW8EbmjZxm/P6WOBW3N6p8ffT38GJYjsboczb4yIvYHXAn9AP7uiwCXA6/KQ4Djgjoi4YQDr2I/y1/GeRtk9wMScfgFwR28zSpon6cY83NgCvHQHfWw6oLnOiPgV8HBjvbD9cfmjlEDTmwXAi4Hbcnf+2B2s++ct7w8GzmmMYxNlr2NiRHwX+DxwHrBR0gWSxuR8fwnMBu7JQ4ZX9bPO9ZGf9HQPZRscTNn+Gxrr/xfKX+S++tu03XZsLHti1m2OiF+31A06SRMkLc7DkUeAr7L9Z+FiyrkW8vWSnK4d/6Db3YIIABHxfcpf8U/30+Ye4P9RfkFvZeAnVB+i7AEc3Cg7iHIIA+UX+KLWmfKcwZeA91B2nccCt1C+fFD2APrzi+Y6Je0F7NtY74BFxNqIOJ7ywTsLWJrL66sPreU/p+w+j238PDsi/iuXf25EHAZMpQSrf8jy6yJiTq73W5S9pb5MbDmncxBlG/yc8pd4v8a6x0TEof30t2m77dhY9nrK4da43BbNul3hf1P6+YcRMYbyOWyO91vAH+X5pGOBS7O8dvyDbrcMIulzwJ9Kelk/bS6mfKlfzdO/pH5FxBOUD/8nJO2dweEDlL8kAF8GPiTpsDzJeEi26fmSPggg6UTKnkiPB4ADJY3qY9WXASdKermk0ZQP4cqIuHsg/W6S9BZJz4uIJ4EtWfxk9u1J4IU7WMQXgVMlHZrL20fSX+X0KyUdKWkPynH9fwNPqlwOfbOkfSLit8Ajua6+PB/4O0l75LJfAlwZERuAq4DPSBoj6Vl5kvJ/DHD4VwIvlvQ3kkZK+mtKsPt2/mFZDXw8+/sa4M/6W5ik0ZL2zLejVG4h6PWEdou9KYffWyVNJANtjyhXEZdSDqdWRcS9WV47/p6T33sCo/L9nvmZastuG0Qi4kHKCdD/1U+zr1NOtF2dv5yBei/lC3In5RLfvwIX5Xq/Bnwiy35J+YsyPiJuBT5DOf5+APhD4IeNZX6Xcrb8fkkP9TKe/wT+Kfu8gbK3M3cn+tw0C1gj6VfAOcDciPhNRDyaff9h7ipP723miPgmZQ9mce6K3wIck9VjKHtcmymHAg8D/5x1bwXuznneCby5nz6uBKZQ9vw+ARwXEQ9n3TzKF+DWXM9SysnIHcplHAt8MPv2YeDYiOjZ5n8DHEk5RFtI+Qz153bgN5TDoWU53bqn05uPU06SbwX+HfhGL20upnxOLmkpb3v86eDsZ8/Vmd/kONqi7Q87zYafpBMoJyBfM9x9GU6SDgJuA34vIh4Z7v70ZbfdEzHrZir3HX0AWNzJAQQGEEQkXaRy49AtjbKXS1qRVxpWSzoiy5U3xayTdJPyJqusm69yI9RaSfMb5Yep3JC0LucdyPGk2W4rT+w+QrmPZOEwd2fHBnAt+SjKsdstjbKrgGNyejZP3yMxm3J9W5T7BFZm+XjK+YPxlJuM7gTGxdPX/6fnPN/pWa5//OOf7vjZ4Z5IRFxLOcm0XTHlBBrAPpTLZlDuZFwUxQpgrKT9KXcxLo+ITRGxGVgOzMq6MRGxIiKCchLrjQMJfmbWGdr956S/B5ZJ+jTlkKjnfw8msv2NLvdlWX/l9/VS3itJJwEnAbAHhw34Ni3gsAMOG3hjM+P6669/KCKet6N27QaRdwHvj4ivS3oT5ZbsN7S5rAGLiAuACwB0gIJ3DHze1QtX76Jeme2eJA3obt12r87M5+nr2l8Djsjp9ZTbvnscmGX9lR/YS7mZdYl2g8gvgJ475I4G1ub0FcC8vEozHdga5SauZcAMSeMkjQNmAMuy7hFJ0/OqzDzKf8eaWZfY4eGMpMso/9C2n6T7KJec3k75B6yRlNuaT8rmV1Ku0Kyj/APYiQARsUnSGcB12e70iOg5Wftuyv+5PJtydWaHz/Mws87RtXes7uw5kVjYneM0Gy6Sro+Iw3fUznesmlkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVqWtvDNZ/l5Jt0laI+lTjfJTM4fM7ZJmNspnZdk6Sac0yidLWpnll/eTi9bMOtBA9kS+Qsnd+hRJr6Okh3hZlGzkn87yqZT8sIfmPF+QNELSCOA8Sr7WqcDx2RZKTtezI+IQSl7RBbWDMrOh027emXcBn4yIx7LNxiyfQ0n791hE3EV5TOIR+bMuIu6MiMeBxcCcfK7q0ZSExFASGDvvjFkXafecyIuBP8nDkO9LemWW72zemX2BLRGxraXczLpEu3lnRlJSYk4HXgkskfTCQetVH7ZLXrXPrl6bmQ1Eu0HkPuAbmfpylaQngf3oO78MfZQ/TEm1OTL3RvrNO/M7yavMbNi1ezjzLeB1AJJeDIwCHqLknZkrabSkycAUSsLu64ApeSVmFOXk6xUZhK4Bjsvlzsd5Z8y6Srt5Zy4CLsrLvo8D8zMgrJG0BLgV2AacHBFP5HLeQ0liNQK4KCLW5Co+AiyWdCZwAyUlp5l1CeedMbNeOe+MmQ0JBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqtJ28Kus+KCkk7ZfvJencTER1k6RpjbbzJa3Nn/mN8sMk3ZzznJtpJMysS7SVvApA0guAGcC9jeJjKM9VnUJ5Kvv52XY85bGKR1Jy0CyUNC7nOR94e2O+31mXmXWudpNXAZwNfBhoPndwDrAoihWUJ7nvD8wElkfEpojYDCwHZmXdmIhYkc9oXYSTV5l1lbbOiUiaA6yPiJ+0VO1s8qqJOd1a3td6T5K0WtJqHm2n52Y22HY674yk5wD/SDmUGVLOO2PWedrZE3kRMBn4iaS7KQmnfizp9+g7eVV/5Qf2Um5mXWKng0hE3BwRz4+ISRExiXIIMi0i7qckr5qXV2mmA1sjYgMl38wMSePyhOoMYFnWPSJpel6VmYeTV5l1lYFc4r0M+BHw+5Luk7Sgn+ZXAncC64AvAe8GiIhNwBmUTHjXAadnGdnmyznPHcB32huKmQ0HJ68ys145eZWZDQkHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWpa28M5L+WdJtmVvmm5LGNupOzRwyt0ua2SiflWXrJJ3SKJ8saWWWXy5p1GAO0Mx2rXbzziwHXhoRfwT8DDgVQNJUYC5waM7zBUkjJI0AzqPkpZkKHJ9tAc4Czo6IQ4DNQH9PTjOzDtNW3pmIuCoituXbFTz9sOU5wOKIeCwi7qI88vCI/FkXEXdGxOPAYmBOPlf1aGBpzn8xzjtj1lUG45zI23j6uag7m3dmX2BLIyA574xZl6kKIpI+CmwDLh2c7vQvIi6IiMMj4nCeMxRrNLMd2enkVT0knQAcC7w+nn7ac1/5Zeij/GFKqs2RuTfivDNmXabdNJqzKHl4/zwimgcWVwBzJY2WNJmSoHsVJU3ElLwSM4py8vWKDD7XAMfl/PNx3hmzrtJu3pnPA3sDyyXdKOmLABGxBlgC3Ar8B3ByRDyRexnvoSSx+imwJNsCfAT4gKR1lHMkFw7qCM1sl3LeGTPrlfPOmNmQcBAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlXazTszXtJySWvzdVyWS9K5mUPmJknTGvPMz/ZrJc1vlB8m6eac59x8AryZdYl2886cAlwdEVOAq/M9lLwyU/LnJOB8KEEHWAgcSUkfsbAn8GSbtzfma12XmXWwtvLOUPLLXJzTzVwxc4BFUaygPIR5f2AmsDwiNkXEZkryq1lZNyYiVuTzVhfhvDNmXaXdp71PiIgNOX0/MCGndzbvzMScbi3vlaSTKHs4sE+bPTezQVV9YjX3IIbkAabOO2PWedoNIg/koQj5ujHL+8o701/5gb2Um1mXaDeIXEHJEQPb54q5ApiXV2mmA1vzsGcZMEPSuDyhOgNYlnWPSJqeV2Xm4bwzZl1lh+dEMu/Ma4H9JN1HucrySWBJ5qC5B3hTNr8SmE1J5P0ocCJARGySdAYliRXA6RHRc7L23ZQrQM+m5PTtyetrZl3AeWfMrFfOO2NmQ8JBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVqQoikt4vaY2kWyRdJmlPSZMlrcw8MpdLGpVtR+f7dVk/qbGcU7P8dkkz64ZkZkOp7SAiaSLwd8DhEfFSYAQwFzgLODsiDgE2AwtylgXA5iw/O9shaWrOdygl58wXJI1ot19mNrRqD2dGAs+WNBJ4DrABOBpYmvWtOWl6ctUsBV6fz1WdAyyOiMci4i7KoxWPqOyXmQ2RtoNIRKwHPg3cSwkeW4HrgS0RsS2bNfPIPJV7Juu3AvvSd06a3yHpJEmrJa3m0XZ7bmaDqeZwZhxlL2IycACwF7s4Babzzph1nprDmTcAd0XEgxHxW+AbwKspqTN7niLfzCPzVO6ZrN8HeJi+c9KYWReoCSL3AtMlPSfPbbweuBW4Bjgu27TmpOnJVXMc8N3MnncFMDev3kymJPVeVdEvMxtC7ebiJSJWSloK/BjYBtwAXAD8O7BY0plZdmHOciFwiaR1lAThc3M5ayQtoQSgbcDJEfFEu/0ys6HlvDNm1ivnnTGzIeEgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyq1OadGStpqaTbJP1U0qskjZe0XNLafB2XbSXp3Mwvc5OkaY3lzM/2ayXN73uNZtZpavdEzgH+IyL+AHgZ8FPgFODqiJgCXJ3vAY6hPPpwCnAScD6ApPHAQuBISqqIhT2Bx8w6X83T3vcBjiIffxgRj0fEFrbPL9Oad2ZRFCsoD3TeH5gJLI+ITRGxGVjOLn5qvJkNnpo9kcnAg8D/lXSDpC9L2guYEBEbss39wISc7iu/jPPOmHWxmiAyEpgGnB8RrwB+zdOHLgDk09wH7eGmzjtj1nlqgsh9wH0RsTLfL6UElQfyMIV83Zj1feWXcd4Zsy5Wk0bzfuDnkn4/i3ryzjTzy7TmnZmXV2mmA1vzsGcZMEPSuDyhOiPLzKwLtJ13Jr0XuFTSKOBO4ERKYFoiaQFwD/CmbHslMJuSsPvRbEtEbJJ0BnBdtjs9IjZV9svMhojzzphZr5x3xsyGhIOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6tSHUQkjcgHNX8730+WtDLzy1yeDyxC0uh8vy7rJzWWcWqW3y5pZm2fzGzoDMaeyPso+WZ6nAWcHRGHAJuBBVm+ANic5WdnOyRNBeYCh1JSRXxB0ohB6JeZDYHaDHgHAv8T+HK+F3A05aHN8Lt5Z3ry0SwFXp/t5wCLI+KxiLiL8vjEI2r6ZWZDp3ZP5HPAh4En8/2+wJaI2Jbvmzlknsovk/Vbs73zzph1sZoMeMcCGyPi+kHsT7+cd8as89Q87f3VwJ9Lmg3sCYyh5OYdK2lk7m00c8j05Je5T9JIYB/gYZx3xqyr1eSdOTUiDoyISZQTo9+NiDcD1wDHZbPWvDM9+WiOy/aR5XPz6s1kSsLvVe32y8yGVm3emd58BFgs6UzgBjLhd75eImkdsIkSeIiINZKWUBJfbQNOjogndkG/zGwXcN4ZM+uV886Y2ZBwEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVWqesfoCSddIulXSGknvy/LxkpZLWpuv47Jcks7N/DI3SZrWWNb8bL9W0vy+1mlmnadmT2Qb8MGImApMB07OHDKnAFdHxBTg6nwPcAzl0YdTgJOA86EEHWAhcCQlVcTCnsBjZp2v5hmrGyLixzn9S0oCq4lsn1+mNe/MoihWUB7ovD8wE1geEZsiYjOwnJLEysy6wKCcE8mUmK8AVgITImJDVt0PTMjpvvLLDDjvjJl1nsHIxftc4OvA30fEI826fJr7oD3c1MmrzDpPbRrNPSgB5NKI+EYWP5CHKeTrxizvK7/MgPPOOHmVWeepuTojShqIn0bEZxtVzfwyrXln5uVVmunA1jzsWQbMkDQuT6jOyDIz6wK1GfDeCtws6cYs+0fgk8ASSQuAe4A3Zd2VwGxKwu5HgRMBImKTpDOA67Ld6RGxqaJfZjaEnHfGzHrlvDNmNiQcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6vSMUFE0ixJt2dyq1N2PIeZdYKOCCKSRgDnURJcTQWOz0RYZtbhOiKIUDLfrYuIOyPicWAxJdmVmXW4mgc1D6beElgd2dpI0kmUFJwAj3Eatwx0BTpNVR3sAPsBDw13J4bYM23MnTbegwfSqFOCyIBExAXABQCSVg/kIbK7i2faeOGZN+ZuHW+nHM4MOIGVmXWWTgki1wFTJE2WNAqYS0l2ZWYdriMOZyJim6T3UDLfjQAuiog1O5jtgl3fs47yTBsvPPPG3JXj7drkVWbWGTrlcMbMupSDiJlV6bogsjvdHi/pbkk3S7pR0uosGy9puaS1+TouyyXp3Bz3TZKmNZYzP9uvlTR/uMbTG0kXSdoo6ZZG2aCNUdJhuQ3X5bzDekNQH+M9TdL6/D3fKGl2o+7U7PvtkmY2ynv9nOfFh5VZfnleiBheEdE1P5STrncALwRGAT8Bpg53vyrGczewX0vZp4BTcvoU4Kycng18BxAwHViZ5eOBO/N1XE6PG+6xNcZzFDANuGVXjBFYlW2V8x7TgeM9DfhQL22n5md4NDA5P9sj+vucA0uAuTn9ReBdw/077rY9kWfC7fFzgItz+mLgjY3yRVGsAMZK2h+YCSyPiE0RsRlYDswa6k73JSKuBTa1FA/KGLNuTESsiPKtWtRY1rDoY7x9mQMsjojHIuIuYB3lM97r5zz3so4Glub8zW03bLotiPR2e/zEYerLYAjgKknX5y39ABMiYkNO3w9MyOm+xt6N22Swxjgxp1vLO9F78hDtop7DN3Z+vPsCWyJiW0v5sOq2ILK7eU1ETKP89/LJko5qVuZf1936GvwzYYzA+cCLgJcDG4DPDG93Ble3BZHd6vb4iFifrxuBb1J2Yx/I3XTydWM272vs3bhNBmuM63O6tbyjRMQDEfFERDwJfInye4adH+/DlEO8kS3lw6rbgshuc3u8pL0k7d0zDcwAbqGMp+fqw3zg33L6CmBeXsGYDmzNQ4JlwAxJ43I3eUaWdbJBGWPWPSJpep4vmNdYVsfoCZjpL+Cp/z6/ApgrabSkycAUyoniXj/nudd2DXBczt/cdsNnuM/stnH2ezbwM8rZ648Od38qxvFCyln3nwBresZCOe69GlgL/CcwPstFeXDTHcDNwOGNZb2NclJuHXDicI+tZZyXUXbhf0s5hl8wmGMEDqd8Ke8APk/ehd1h470kx3MTJXDs32j/0ez77TSuLPX1Oc/PzarcDl8DRg/379i3vZtZlW47nDGzDuMgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKr8f+I5dlYt81I8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "stage = 100\n",
    "cond = 1000\n",
    "rbot = 90\n",
    "\n",
    "rivSp = {0: []}\n",
    "\n",
    "for r in range(0, 18):\n",
    "    rivSp[0].append([0, r, 0, stage, cond, rbot])\n",
    "    \n",
    "riv = flopy.modflow.ModflowRiv(ml, stress_period_data=rivSp)    \n",
    "\n",
    "riv.plot(mflay=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GHB-Package \n",
    "\n",
    "The general head boundary is located in the right side of the model and takes place in both layers.   \n",
    "\n",
    "So we'll add two for loops, iterating through layers and rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.axes._subplots.AxesSubplot at 0x7fdd59b83208>,\n",
       " <matplotlib.axes._subplots.AxesSubplot at 0x7fdd59b29668>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEICAYAAAB4TcDdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGydJREFUeJzt3X20XVV57/Hvj8QEQV4SgpRXiRJKg1oNCKl6FYGbhOgwegelcVgTgcpVoLW9tAjaWyjQCkpFMxQ6UBEClBAjtdhqMSLWMa4mEOQ1KCZAkJdAgLyAUMDAc/+Yz4GVzTkn5+x5cs7e5PcZY4+99lxzrTXnfnnOWnOtsx5FBGZm7dpmpBtgZt3NQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDSBsknSnpimHe5j6Sfitp1HBu99VI0nJJh7W5bEjabwjacJikB2vX0wm6PohImi1pqaSnJa3J6RMlKedfKumclmX2zS/D6Ead5/NH+pSkmyW9dyT602jjKklH9ryOiN9ExOsi4oVhbMMr3rtXg4g4MCJ+MtTrlXSMpJ9JekbSkK9/S5F0tqQ7JG2UdOZgl+/qICLpFOArwBeB3wN2Az4JvAsYM8jVfSEiXgfsCFwEXOO/+v3rCcLdYhjauxb4MnDuFt5OW/rp/0rgVOA/2llv1wYRSTsBZwEnRsSiiHgqilsi4qMR8Vw7641yCe+/AOMpQWkgbflg7iKvl/QTSX/QmLe3pGskPSbpCUlfzfI3Sfpxlj0u6UpJO+e8y4F9gO/l3tGpvew97SHpWklrJa2U9InGNs+UtFDS/NyzWi7p4D7aLkkX5F7ck/kX6c2STgA+Cpyabfhe1l8l6TOSbgeeljQ62/Kd7ON9kv6isf5DJC3LdT8q6UtZvq2kK7L/6yXdJKnX9zu3ebqkuyStk/QtSds25n9A0q25np9JemvLsq3tfWkvT9JYSV+W9HA+vixpbGP5v5G0Oucd19/3ICJ+FBELgYf7q9dHH0+TdE9+XndJ+nCWj8nP+C2Nuq/PvZ1d2+l/L+2+LCJ+ADw12Hb3rKArH8AMYCMwejP1LgXOaSnbF4ieZZt1gFGUvZl7gVF9rPNM4Iqc3h94GvifwGsoEX0lZU9oFHAbcAGwPbAt8O5cbr9cZiywK/BT4MuNbawCjuynzT8FLsx1vg14DDi80b5ngZnZhs8DS/roy3TgZmBnQMAfALv3896tAm4F9gZeS/lDdDPwd9nnN+Z7Nz3r/xz4WE6/Dpia0/8b+B6wXbbxIGDHPtq4Crgztzke+H+Nz+vtwBrg0FzP3Kw/trf2tr63lD9ES4DX5+fwM+DsxnfsUeDN+fn9S34G+23mO/dnwE82U+cw4MHG6z8G9sj3808o36mez+FC4LxG3U8D32u3//206QrgzEH/Fkc6GLT7AP4UeKSl7GfAeuC/gfc0fgjPZnnP40leGUSebSz7LPDRfrZ9Ji8Hkf8LLGzM2wZ4KL8kf0T5cfcb6HK5DwG3tPxweg0i+YV4AdihMf/zwKWN9v2oMW8y8N99bPdw4NfAVGCblnmX0nsQOa7x+lDgNy11Tge+ldM/Bf4emNBS57j8vN46gPdmFfDJxuuZwD05fRH5o2/Mvxt4b2/tbX1vgXuAmY1504FVOX0JcG5j3v5soSDSy/xbgVnN9xhQvl4GHNNu//vZZltBpGsPZ4AngAnN3bOIeGdE7Jzzmn07PyJ27nkAb+WVzs952wEHA1+UdNQA2rEHcH+jDS8CDwB7Un7s90fExtaFJO0maYGkhyQ9SfkAJwxgez3bXBsRzd3P+3ObPR5pTD8DbNvHruyPga8CXwPWSLpY0o6b2f4Djek3AHvkrvR6SeuBz/LyoeDxlB/fr/KQ5QNZfjlwHbAgDxW+IOk1A9zm/ZT3oGf7p7Rsf+/G/NZlW23y+bWse49etrtFSJrTOCRZT9n7mQAQEUspn+Fhkg6g7MVem4vW9r9aNweRnwPPAbOGcqVR3EnZZX7/ABZ5mPJBAmWMgfIhPkT58Pbp7ccL/CPlr9pbImJHyp6Vmk3ZzDbHS9qhUbZPbnPQImJeRBxE2WPZH/ibzbShWf4AcF8zSEfEDhExM9e9IiI+QjlcOA9YJGn7iPhdRPx9REwG3gl8AJjTTzP3bkzvw8vjDg8A/9Cy/e0i4qo+2ttqk8+vZd2re9nukJP0BuDrwMnALvnH7E42/T5cRvmOfAxYFBHPZnlt/6t1bRCJiPWU3eQLJR0taQdJ20h6G+X4tW0Z7d8NLB9A9YXA+yUdkX9JT6EEt58BN1K+iOdK2j4HE9+Vy+0A/BbYIGlPXv7h9niUMr7wChHxQK7/87nOt1L+4g/62hVJ75B0aLb9acqh3Iuba0PDjcBTOXj3WkmjcmD2Hbn+P5W0a+6hrc9lXpT0PklvUTkD9iTwu8Z2e3OSpL0kjQc+B1yd5V8HPpl9UL7P728JsP25CvhbSbtKmkAZ2+l5HxcCH5c0WdJ2wBn9rSj7vi3lkHOb/Gz627vqsT3lh/5YrudYyp5I0xXAhymBZH6jvLb/SHpNtnsbYHS2e+BnJgd7/NNpD8oZhBspu3uPAUuBE4AxOf9SBjaw+jzlR/005fjzH2kZI2gsfyY5JpKvPwzcBWwA/gs4sDFvH+C7lEOsx4F5WX4gZUDyt5Tj31PYdKBtVrZjPfDXvbR5L+DfKacV72HTMYPW9m2ybEtfjgBuz3Y8DlwJvC7nTcq2rQe+Gy8fYx/Zso49KD/GR4B1lIHKnjGHKygDf7+lBOUPZflHKMfuT1OC1bze2tfY5un5Hq+n/FXerjF/BnBTzlsNfJscL+qjvS+VUQam5+Vyq3N620bd07JfD1PGcfocEwE+nvObj0v7qHtYy+f9D/lZPg58ifI9+rOWZX6UbVdL+aD630tbLu2l3R8f6G+wZ6DGrGNJWkX5Qf1opNsykiRdAjwcEX870m1p6qqLhcy2VpL2Bf4X5ZRuR9nsmIikS1QuRLqzUfY2SUtyNHmZpEOyXJLmqVz8dLukKY1l5kpakY+5jfKDVC5wWpnLCjN7iaSzKQOtX4yI+0a6Pa8wgDGH9wBTgDsbZT8EjsrpmeQ58Zz+AWVUeSqwNMvHUy5AGg+My+lxOe/GrKtc9qjhHFPxww8/6h6b3ROJiJ9SBnw2Kab8jwnATrx8SmwWMD+KJcDOknanXMCzOCLWRsQ6YDEwI+ftGBFLIiIoo84fGkjwM7PO0O6YyF8C10k6n3JI9M4s35NNL2x5MMv6K3+wl/Jeqfw/xwkA22+//UEHHHBAm8032/rcPOgFbn48InbdXLV2g8ingL+KiO9IOgb4JnDkZpapFhEXAxcDHHzwwbFs2bItvUmzV41BDzZKA7pCt92LzeYC1+T0t4FDcvohNr3Cb68s6698r17KzaxLtBtEHgbem9OHAyty+lpgTp6lmQpsiIjVlP+RmCZpnKRxwDTgupz3pKSpeVZmDvBv7XbGzIbfZg9nJF1Fubpugsrt3M4APgF8Jf8n5FlynAL4PuUMzUrKFaTHAkTE2jxNdVPWOysiegZrT6RcMfdaytmZH1T3ysyGTddeseoxEbPBaWNM5OaI6PVmVk1d+w94ZtYZHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWZW28s5k+Z9L+pWk5ZK+0Cg/PXPI3C1peqN8RpatlHRao3yipKVZfrWkMUPVOTPb8gayJ3IpJdfnSyS9j5Ie4g8j4kDg/CyfDMym5JmdQUm2PSqTA38NOIqSef4jWRdKpvgLImI/Sh7X42s7ZWbDp928M58Czo2I57LOmiyfBSyIiOeiZOpaSbmJ8yHAyoi4NyKeBxYAs/K+qocDi3L5y3DeGbOu0u6YyP7A/8jDkP+S9I4sH2zemV2A9RGxsaXczLpEu3lnRlNSYk4F3gEslPTGIWtVH5rJq/bZZ58tvTkzG4B290QeBK7JdJk3Ai8CExh83pknKKk2R7eU9yoiLo6IgyPi4F133WxiLjMbBu0Gke8C7wOQtD8wBnickndmtqSxkiYCkygJu28CJuWZmDGUwddrM//uDcDRud65OO+MWVdpN+/MJcAledr3eWBuBoTlkhYCdwEbgZMi4oVcz8mUJFajgEsiYnlu4jPAAknnALdQUnKaWZdw3hmzrYTzzphZR3IQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrErbyaty3imSQtKEfC1J8zIR1e2SpjTqzpW0Ih9zG+UHSbojl5mXaSTMrEu0lbwKQNLewDTgN43ioyj3VZ1EuSv7RVl3POW2iodSctCcIWlcLnMR8InGcq/Ylpl1rnaTVwFcAJwKNO+vOAuYn3eBX0K5k/vuwHRgcUSsjYh1wGJgRs7bMSKW5D1a5+PkVWZdpa0xEUmzgIci4raWWYNNXrVnTreW97XdEyQtk7Tssccea6fpZjbEBh1EJG0HfBb4u6FvTv+cd8as87SzJ/ImYCJwm6RVlIRTv5D0eww+edVDOd1abmZdYtBBJCLuiIjXR8S+EbEv5RBkSkQ8QkleNSfP0kwFNkTEakq+mWmSxuWA6jTgupz3pKSpeVZmDk5eZdZVBnKK9yrg58DvS3pQ0vH9VP8+cC+wEvg6cCJARKwFzqZkwrsJOCvLyDrfyGXuAX7QXlfMbCQ4eZXZVsLJq8ysIzmImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqbeWdkfRFSb/K3DL/KmnnxrzTM4fM3ZKmN8pnZNlKSac1yidKWprlV0saM5QdNLMtq928M4uBN0fEW4FfA6cDSJoMzAYOzGUulDRK0ijga5S8NJOBj2RdgPOACyJiP2Ad0N+d08ysw7SVdyYifhgRG/PlEl6+2fIsYEFEPBcR91FueXhIPlZGxL0R8TywAJiV91U9HFiUy1+G886YdZWhGBM5jpfvizrYvDO7AOsbAcl5Z8y6TFUQkfQ5YCNw5dA0p3/OO2PWeUa3u6CkjwMfAI6Il+/23Fd+Gfoof4KSanN07o0474xZl2k3jeYMSh7eD0bEM41Z1wKzJY2VNJGSoPtGSpqISXkmZgxl8PXaDD43AEfn8nNx3hmzrtJu3pmvAjsAiyXdKumfASJiObAQuAv4T+CkiHgh9zJOpiSx+iWwMOsCfAb4P5JWUsZIvjmkPTSzLcp5Z8y2Es47Y2YdyUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVpN+/MeEmLJa3I53FZLknzMofM7ZKmNJaZm/VXSJrbKD9I0h25zLy8A7yZdYl2886cBlwfEZOA6/M1lLwyk/JxAnARlKADnAEcSkkfcUZP4Mk6n2gs17otM+tgbeWdoeSXuSynm7liZgHzo1hCuQnz7sB0YHFErI2IdZTkVzNy3o4RsSTvtzof550x6yrtjonsFhGrc/oRYLecHmzemT1zurW8V847Y9Z5qgdWcw9iWG7U6rwzZp2n3SDyaB6KkM9rsryvvDP9le/VS7mZdYl2g8i1lBwxsGmumGuBOXmWZiqwIQ97rgOmSRqXA6rTgOty3pOSpuZZmTk474xZV9lsBrzMO3MYMEHSg5SzLOcCCzMHzf3AMVn9+8BMSiLvZ4BjASJiraSzKUmsAM6KiJ7B2hMpZ4BeS8np25PX18y6gPPOmG0lnHfGzDqSg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq1IVRCT9laTlku6UdJWkbSVNlLQ088hcLWlM1h2br1fm/H0b6zk9y++WNL2uS2Y2nNoOIpL2BP4CODgi3gyMAmYD5wEXRMR+wDrg+FzkeGBdll+Q9ZA0OZc7kJJz5kJJo9ptl5kNr9rDmdHAayWNBrYDVgOHA4tyfmtOmp5cNYuAI/K+qrOABRHxXETcR7m14iGV7TKzYdJ2EImIh4Dzgd9QgscG4GZgfURszGrNPDIv5Z7J+RuAXeg7J80rOO+MWeepOZwZR9mLmAjsAWzPFk6B6bwzZp2n5nDmSOC+iHgsIn4HXAO8i5I6s+cu8s08Mi/lnsn5OwFP0HdOGjPrAjVB5DfAVEnb5djGEcBdwA3A0VmnNSdNT66ao4EfZ/a8a4HZefZmIiWp940V7TKzYbTZvDN9iYilkhYBvwA2ArcAFwP/ASyQdE6WfTMX+SZwuaSVlAThs3M9yyUtpASgjcBJEfFCu+0ys+HlvDNmWwnnnTGzjuQgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyq1Oad2VnSIkm/kvRLSX8kabykxZJW5PO4rCtJ8zK/zO2SpjTWMzfrr5A0t+8tmlmnqd0T+QrwnxFxAPCHwC+B04DrI2IScH2+BjiKcuvDScAJwEUAksYDZwCHUlJFnNETeMys89Xc7X0n4D3k7Q8j4vmIWM+m+WVa887Mj2IJ5YbOuwPTgcURsTYi1gGL2cJ3jTezoVOzJzIReAz4lqRbJH1D0vbAbhGxOus8AuyW033ll3HeGbMuVhNERgNTgIsi4u3A07x86AJA3s19yG7i6rwzZp2nJog8CDwYEUvz9SJKUHk0D1PI5zU5v6/8Ms47Y9bFatJoPgI8IOn3s6gn70wzv0xr3pk5eZZmKrAhD3uuA6ZJGpcDqtOyzMy6QNt5Z9KfA1dKGgPcCxxLCUwLJR0P3A8ck3W/D8ykJOx+JusSEWslnQ3clPXOioi1le0ys2HivDNmWwnnnTGzjuQgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqVAcRSaPyRs3/nq8nSlqa+WWuzhsWIWlsvl6Z8/dtrOP0LL9b0vTaNpnZ8BmKPZFPU/LN9DgPuCAi9gPWAcdn+fHAuiy/IOshaTIwGziQkiriQkmjhqBdZjYMajPg7QW8H/hGvhZwOOWmzfDKvDM9+WgWAUdk/VnAgoh4LiLuo9w+8ZCadpnZ8KndE/kycCrwYr7eBVgfERvzdTOHzEv5ZXL+hqzvvDNmXawmA94HgDURcfMQtqdfzjtj1nlq7vb+LuCDkmYC2wI7UnLz7ixpdO5tNHPI9OSXeVDSaGAn4Amcd8asq9XknTk9IvaKiH0pA6M/joiPAjcAR2e11rwzPflojs76keWz8+zNRErC7xvbbZeZDa/avDO9+QywQNI5wC1kwu98vlzSSmAtJfAQEcslLaQkvtoInBQRL2yBdpnZFuC8M2ZbCeedMbOO5CBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKrU3GN1b0k3SLpL0nJJn87y8ZIWS1qRz+OyXJLmZX6Z2yVNaaxrbtZfIWluX9s0s85TsyeyETglIiYDU4GTMofMacD1ETEJuD5fAxxFufXhJOAE4CIoQQc4AziUkirijJ7AY2adr+Yeq6sj4hc5/RQlgdWebJpfpjXvzPwollBu6Lw7MB1YHBFrI2IdsJiSxMrMusCQjIlkSsy3A0uB3SJidc56BNgtp/vKLzPgvDNm1nmGIhfv64DvAH8ZEU825+Xd3IfsJq5OXmXWeWrTaL6GEkCujIhrsvjRPEwhn9dkeV/5ZQacd8bJq8w6T83ZGVHSQPwyIr7UmNXML9Oad2ZOnqWZCmzIw57rgGmSxuWA6rQsM7MuUJsB72PAHZJuzbLPAucCCyUdD9wPHJPzvg/MpCTsfgY4FiAi1ko6G7gp650VEWsr2mVmw8h5Z8y2Es47Y2YdyUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKh0TRCTNkHR3Jrc6bfNLmFkn6IggImkU8DVKgqvJwEcyEZaZdbiOCCKUzHcrI+LeiHgeWEBJdmVmHa7mRs1DqbcEVoe2VpJ0AiUFJ8Bzku4chrZ1ignA4yPdiGG2tfW50/r7hoFU6pQgMiARcTFwMYCkZQO5ieyrxdbWX9j6+tyt/e2Uw5kBJ7Ays87SKUHkJmCSpImSxgCzKcmuzKzDdcThTERslHQyJfPdKOCSiFi+mcUu3vIt6yhbW39h6+tzV/a3a5NXmVln6JTDGTPrUg4iZlal64LIq+nyeEmrJN0h6VZJy7JsvKTFklbk87gsl6R52e/bJU1prGdu1l8hae5I9ac3ki6RtKZ5Tc9Q9lHSQfkersxlB51ydij10d8zJT2Un/OtkmY25p2ebb9b0vRGea/f8zz5sDTLr84TESMrIrrmQRl0vQd4IzAGuA2YPNLtqujPKmBCS9kXgNNy+jTgvJyeCfyAkpd5KrA0y8cD9+bzuJweN9J9a/TnPcAU4M4t0UfgxqyrXPaoDuzvmcBf91J3cn6HxwIT87s9qr/vObAQmJ3T/wx8aqQ/427bE9kaLo+fBVyW05cBH2qUz49iCbCzpN2B6cDiiFgbEeuAxcCM4W50XyLip8DaluIh6WPO2zEilkT5Vc1vrGtE9NHfvswCFkTEcxFxH7CS8h3v9Xuee1mHA4ty+eZ7N2K6LYj0dnn8niPUlqEQwA8l3ZyX9APsFhGrc/oRYLec7qvv3fieDFUf98zp1vJOdHIeol3Sc/jG4Pu7C7A+Ija2lI+obgsirzbvjogplP9ePknSe5oz86/rq/oc/NbQR+Ai4E3A24DVwD+NbHOGVrcFkVfV5fER8VA+rwH+lbIb+2juppPPa7J6X33vxvdkqPr4UE63lneUiHg0Il6IiBeBr1M+Zxh8f5+gHOKNbikfUd0WRF41l8dL2l7SDj3TwDTgTkp/es4+zAX+LaevBebkGYypwIY8JLgOmCZpXO4mT8uyTjYkfcx5T0qamuMFcxrr6hg9ATN9mPI5Q+nvbEljJU0EJlEGinv9nude2w3A0bl8870bOSM9stvG6PdM4NeU0evPjXR7KvrxRsqo+23A8p6+UI57rwdWAD8Cxme5KDduuge4Azi4sa7jKINyK4FjR7pvLf28irIL/zvKMfzxQ9lH4GDKj/Ie4KvkVdgd1t/Lsz+3UwLH7o36n8u2303jzFJf3/P83tyY78O3gbEj/Rn7snczq9JthzNm1mEcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFmV/w/OpxgZnd8/sgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEICAYAAAB4TcDdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGy5JREFUeJzt3X20XVV57/Hvj8QECCAJYMqrRAmlAS2GCKlaReAmITqM3svlhmFNBCpXhdb20qugvYWCraBc0QyFDlSEqJeQRmppi8WIWMaoJrwoAkExAQIEA0HygoKAgef+MZ9jVjbn5Jyz58k5e5PfZ4w99tpzzbXWnHuv/Zy15tpnPYoIzMzatdNIN8DMupuDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQeRNkg6X9LXh3mbB0n6taRRw7ndlyNJKyQd2+ayIemQIWjDsZLW1K6nE3R9EJE0V9JySU9LWpfTH5aknH+VpE+2LHNw7gyjG3Wezy/pryTdIeltI9GfRhtXSzqh53VEPBwRu0XEC8PYhpe8dy8HEXF4RHx/qNcr6WRJP5D0jKQhX//2IOlVkq6R9AtJmyT9p6RjBrOOrg4iks4GPg98Bvg9YCLwQeDNwJhBru7TEbEbsAdwOXCd/+pvW08Q7hbD0N71wOeAi7bzdtrSR/93A24DjgImAFcD/yZptwGvOCK68gG8Enga+G/91LsK+GRL2cFAAKN7qwPsmvP362Od5wNfb7x+F7AC2Ah8H/iDxrwDgeuAJ4AngS9k+WuB72XZL4FvAHvmvK8BLwK/AX4NfLSXNu8HXE/ZcVcBH2hp32JgIfCrbNu0Pvoi4FJgHfAUcDdwBHAG8Fvg+WzDv2T91cDHgLuA54DR2ZZvZh8fBP68sf6jgdtz3Y8Dn83ynYGvZ/83UnbkiX20cTVwLnAvsAH4KrBzY/47gTtzPT8AXt+ybGt7VwMn5PyxlC/+L/LxOWBsY/n/DazNeaflZ3BIP/vcnwLf76fOscCaxutzgPvz87oXeE+Wj8nP+HWNuq8CngH2aaf/A/huPQUcNeDv4kgHg3YfwCxgc39vCoMMIsAoytHMA8CoPtZ5PhlEgEMpwey/AK+gfOFX5Yc/CvgJ5Us6Lr84b8nlDsllxgL7ALcAn2v58E/YRptvAS7LdR5J+QIf12jfs8DsbMOngGV99GUmcAewJyWg/AGw7zbeu9W5wx4I7EI5mr0D+Jvs82vyvZuZ9X8IvC+ndwOm5/T/BP6FErBHUf4S7tFHG1cD9+Q2JwD/2fi83kAJgMfkeuZn/bG9tbf1vQUuAJZRvpj7UL6EFzb2sccpQXUc8P/YfkHkv1OC8U7A/6DsUz2fw2XAxY26H2FLUB90//tp15G577xywN/FkQ4G7T6APwEeayn7ASUa/wZ4a+OL8GyW9zye4qVB5NnGss8C793Gts9nSxD5P8DixrydgEdzJ/kjypd7INH/3cCPW744vQaR3CFeAHZvzP8UcFWjfd9tzJsC/KaP7R4H/ByYDuzUMu8qeg8ipzVeHwM83FLnXOCrOX0L8LfA3i11TqPlr+Y23pvVwAcbr2cD9+f05eSXvjH/PuBtvbW39b2l/PWf3Zg3E1id01cCFzXmHcp2CiK9zL8TmNN8jwHl69uBk9vt/za2uQflSPTcwXwXu3lM5Elg7+Z5XkS8KSL2zHnNvl0SEXv2PIDX97K+S3LersA04DOSThxAO/YDHmq04UXgEWB/ypf9oYjY3LqQpImSFkl6VNJTlEP7vQewvZ5tro+IXzXKHspt9nisMf0MsHNv58QR8T3gC8AXgXWSrpC0Rz/bf6Qx/WpgP0kbex7AxynjUwCnU758P5N0m6R3ZvnXgBuBRTmo92lJrxjgNh+ivAc92z+7ZfsHNua3Lttqq8+vZd379bLd7ULSPEl3NvpwBLk/RMRyymd4rKTDKEex1+eitf3v2f4ulCPDZRHxqcG0vZuDyA8p53hzhnKlUdxDOWR+xwAW+QXlgwQgrwodSDkaeQQ4qI8Brb+n/FV7XUTsQTmyUrMp/WxzgqTdG2UH5TYHLSIWRMRRlCOWQynjANtqQ7P8EeDBZpCOiN0jYnaue2VEnEI5XbgYWCJpXET8NiL+NiKmAG+inNfP20YzD2xMH0R5D3q2/3ct2981Iq7po72ttvr8Wta9tpftDjlJrwa+BJwF7JV/zO5h6/3haso+8j5gSUQ8m+W1/UfSWOBbwBrKaeagdG0QiYiNlMPkyySdJGl3STtJOpJy/tq2jPZvoQxI9mcx8A5Jx+df0rMpwe0HwK2UHfEiSeMk7Szpzbnc7pQBy02S9mfLF7fH45TxhZeIiEdy/Z/Kdb6e8hd/0L9dkfRGScdk25+mnMq92F8bGm4FfiXpY5J2kTRK0hGS3pjr/xNJ++QR2sZc5kVJb5f0urwC9hRlEPfF3jcBwJmSDpA0AfgEcG2Wfwn4YPZB+T6/oyXAbss1wF9L2kfS3pSxnZ73cTHwfklTJO0KnLetFWXfd6accu6Un822jq56jKN80Z/I9ZxKORJp+jrwHkogWdgor+p/tm8J5TR+fn5OgzOYc59OfADvpezIz+SHsJxyZWFMzr+KgQ2s9lyFeJpy/vn3tIwRNJY/n62vzryHMqK+CfgP4PDGvIMoUb7nKsyCLD+cMiD5a8r579lsPdA2J9uxEfirXtp8APCvlJH7+9l6zKC1fVst29KX4ykj979my1Wi3XLeZLaM+n8rtpxjn9Cyjv0oX8bHKFdPlrFlzOHrlIG/X1OC8ruz/BTKufvTlGC1oLf2NbbZc3VmI+Wv8q6N+bMoV3c2UoL2P5LjRX2093dllIHpBbnc2pxuXvk5J/vV79UZ4P05v/m4qo+6x7Z83n+Xn+Uvgc9S9qM/bVnmu9l2tZQPqv8ty74t2/lMfkY9jz8e6HewZ6DGrGNJWk35Qn13pNsykiRdCfwiIv56pNvS1FU/FjLbUUk6GPivlEu6HaXfMRFJV6r8nPyeRtmRkpblaPLtko7OcklaIGmVpLskTW0sM1/SynzMb5QfJenuXGZBDkyaWZJ0IWWg9TMR8eBIt+clBjDm8FZgKnBPo+w7wIk5PZu8Jp7T36aMKk8Hlmf5BMoPkCYA43N6fM67Nesqlz1xuMdV/PDDj/Yf/R6JRMQtlAGfrYopP0yB8vPznktic4CFUSwD9pS0L+UHPEsjYn1EbACWArNy3h4RsSwigjLq/O6BBD8z6wztjon8BXCjpEsop0RvyvL92fqHLWuybFvla3op75WkMyhXXhg3btxRhx12WJvNN9vx3DHoBe74ZUTs01+1doPIh4C/jIhvSjoZ+ApwQj/LVIuIK4ArAKZNmxa333779t6k2cvGoAcbpQH9QrfdH5vNp/xnKpRr0kfn9KNs/Qu/A7JsW+UH9FJuZl2i3SDyC8qPVKD8A9fKnL4emJdXaaYDmyJiLeV/JGZIGi9pPDADuDHnPSVpel6VmQf8c7udMbPh1+/pjKRrKL+u21vldm7nAR8APp//E/IsOU4B3EC5QrOK8gu4UwEiYn1eprot610QET2DtR+m/GJ0F8rVmW9X98rMhk3X/mLVYyJmg9PGmMgdETGtv2pd+w94ZtYZHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWZW28s5k+Z9J+pmkFZI+3Sg/N3PI3CdpZqN8VpatknROo3ySpOVZfq2kMUPVOTPb/gZyJHIVJdfn70h6OyU9xB9GxOHAJVk+BZhLyTM7i5Jse1Qmbf4icCIl8/wpWRdKpvhLI+IQSh7X02s7ZWbDp928Mx8CLoqI57LOuiyfAyyKiOeiZOpaRbmJ89HAqoh4ICKeBxYBc/K+qsdRspJDSdTsvDNmXaTdMZFDgT/O05D/kPTGLB9s3pm9gI0Rsbml3My6RLt5Z0ZTUmJOB94ILJb0miFrVR+ayasOOuig7b05MxuAdo9E1gDXZbrMW4EXgb0ZfN6ZJympNke3lPcqIq6IiGkRMW2fffpNzGVmw6DdIPIt4O0Akg4FxgC/pOSdmStprKRJwGRKwu7bgMl5JWYMZfD1+sy/ezNwUq53Ps47Y9ZV2s07cyVwZV72fR6YnwFhhaTFwL3AZuDMiHgh13MWJYnVKODKiFiRm/gYsEjSJ4EfU1JymlmXcN4Zsx2E886YWUdyEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKxK28mrct7ZkkLS3vlakhZkIqq7JE1t1J0vaWU+5jfKj5J0dy6zINNImFmXaCt5FYCkA4EZwMON4hMp91WdTLkr++VZdwLltorHUHLQnCdpfC5zOfCBxnIv2ZaZda52k1cBXAp8FGjeX3EOsDDvAr+Mcif3fYGZwNKIWB8RG4ClwKyct0dELMt7tC7EyavMukpbYyKS5gCPRsRPWmYNNnnV/jndWt7Xds+QdLuk25944ol2mm5mQ2zQQUTSrsDHgb8Z+uZsm/POmHWedo5EXgtMAn4iaTUl4dSPJP0eg09e9WhOt5abWZcYdBCJiLsj4lURcXBEHEw5BZkaEY9RklfNy6s004FNEbGWkm9mhqTxOaA6A7gx5z0laXpelZmHk1eZdZWBXOK9Bvgh8PuS1kg6fRvVbwAeAFYBXwI+DBAR64ELKZnwbgMuyDKyzpdzmfuBb7fXFTMbCU5eZbaDcPIqM+tIDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrEpbeWckfUbSzzK3zD9J2rMx79zMIXOfpJmN8llZtkrSOY3ySZKWZ/m1ksYMZQfNbPtqN+/MUuCIiHg98HPgXABJU4C5wOG5zGWSRkkaBXyRkpdmCnBK1gW4GLg0Ig4BNgDbunOamXWYtvLORMR3ImJzvlzGlpstzwEWRcRzEfEg5ZaHR+djVUQ8EBHPA4uAOXlf1eOAJbn81TjvjFlXGYoxkdPYcl/Uwead2QvY2AhIzjtj1mWqgoikTwCbgW8MTXO2zXlnzDrP6HYXlPR+4J3A8bHlbs995Zehj/InKak2R+fRiPPOmHWZdtNozqLk4X1XRDzTmHU9MFfSWEmTKAm6b6WkiZicV2LGUAZfr8/gczNwUi4/H+edMesq7ead+QKwO7BU0p2S/gEgIlYAi4F7gX8HzoyIF/Io4yxKEqufAouzLsDHgP8laRVljOQrQ9pDM9uunHfGbAfhvDNm1pEcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyoOImZWxUHEzKo4iJhZlXbzzkyQtFTSynwen+WStCBzyNwlaWpjmflZf6Wk+Y3yoyTdncssyDvAm1mXaDfvzDnATRExGbgpX0PJKzM5H2cAl0MJOsB5wDGU9BHn9QSerPOBxnKt2zKzDtZW3hlKfpmrc7qZK2YOsDCKZZSbMO8LzASWRsT6iNhASX41K+ftERHL8n6rC3HeGbOu0u6YyMSIWJvTjwETc3qweWf2z+nW8l4574xZ56keWM0jiGG5Uavzzph1nnaDyON5KkI+r8vyvvLObKv8gF7KzaxLtBtErqfkiIGtc8VcD8zLqzTTgU152nMjMEPS+BxQnQHcmPOekjQ9r8rMw3lnzLpKvxnwMu/MscDektZQrrJcBCzOHDQPASdn9RuA2ZRE3s8ApwJExHpJF1KSWAFcEBE9g7UfplwB2oWS07cnr6+ZdQHnnTHbQTjvjJl1JAcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlalKohI+ktJKyTdI+kaSTtLmiRpeeaRuVbSmKw7Nl+vyvkHN9ZzbpbfJ2lmXZfMbDi1HUQk7Q/8OTAtIo4ARgFzgYuBSyPiEGADcHoucjqwIcsvzXpImpLLHU7JOXOZpFHttsvMhlft6cxoYBdJo4FdgbXAccCSnN+ak6YnV80S4Pi8r+ocYFFEPBcRD1JurXh0ZbvMbJi0HUQi4lHgEuBhSvDYBNwBbIyIzVmtmUfmd7lncv4mYC/6zknzEs47Y9Z5ak5nxlOOIiYB+wHj2M4pMJ13xqzz1JzOnAA8GBFPRMRvgeuAN1NSZ/bcRb6ZR+Z3uWdy/iuBJ+k7J42ZdYGaIPIwMF3Srjm2cTxwL3AzcFLWac1J05Or5iTge5k973pgbl69mURJ6n1rRbvMbBj1m3emLxGxXNIS4EfAZuDHwBXAvwGLJH0yy76Si3wF+JqkVZQE4XNzPSskLaYEoM3AmRHxQrvtMrPh5bwzZjsI550xs47kIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqtTmndlT0hJJP5P0U0l/JGmCpKWSVubz+KwrSQsyv8xdkqY21jM/66+UNL/vLZpZp6k9Evk88O8RcRjwh8BPgXOAmyJiMnBTvgY4kXLrw8nAGcDlAJImAOcBx1BSRZzXE3jMrPPV3O39lcBbydsfRsTzEbGRrfPLtOadWRjFMsoNnfcFZgJLI2J9RGwAlrKd7xpvZkOn5khkEvAE8FVJP5b0ZUnjgIkRsTbrPAZMzOm+8ss474xZF6sJIqOBqcDlEfEG4Gm2nLoAkHdzH7KbuDrvjFnnqQkia4A1EbE8Xy+hBJXH8zSFfF6X8/vKL+O8M2ZdrCaN5mPAI5J+P4t68s4088u05p2Zl1dppgOb8rTnRmCGpPE5oDojy8ysC7Sddyb9GfANSWOAB4BTKYFpsaTTgYeAk7PuDcBsSsLuZ7IuEbFe0oXAbVnvgohYX9kuMxsmzjtjtoNw3hkz60gOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOr4iBiZlUcRMysSnUQkTQqb9T8r/l6kqTlmV/m2rxhEZLG5utVOf/gxjrOzfL7JM2sbZOZDZ+hOBL5CCXfTI+LgUsj4hBgA3B6lp8ObMjyS7MekqYAc4HDKakiLpM0agjaZWbDoDYD3gHAO4Av52sBx1Fu2gwvzTvTk49mCXB81p8DLIqI5yLiQcrtE4+uaZeZDZ/aI5HPAR8FXszXewEbI2Jzvm7mkPldfpmcvynrO++MWReryYD3TmBdRNwxhO3ZJuedMes8NXd7fzPwLkmzgZ2BPSi5efeUNDqPNpo5ZHryy6yRNBp4JfAkzjtj1tVq8s6cGxEHRMTBlIHR70XEe4GbgZOyWmvemZ58NCdl/cjyuXn1ZhIl4fet7bbLzIZXbd6Z3nwMWCTpk8CPyYTf+fw1SauA9ZTAQ0SskLSYkvhqM3BmRLywHdplZtuB886Y7SCcd8bMOpKDiJlVcRAxsyoOImZWxUHEzKo4iJhZFQcRM6viIGJmVRxEzKyKg4iZVXEQMbMqDiJmVsVBxMyqOIiYWRUHETOrUnOP1QMl3SzpXkkrJH0kyydIWippZT6Pz3JJWpD5Ze6SNLWxrvlZf6Wk+X1t08w6T82RyGbg7IiYAkwHzswcMucAN0XEZOCmfA1wIuXWh5OBM4DLoQQd4DzgGEqqiPN6Ao+Zdb6ae6yujYgf5fSvKAms9mfr/DKteWcWRrGMckPnfYGZwNKIWB8RG4CllCRWZtYFhmRMJFNivgFYDkyMiLU56zFgYk73lV9mwHlnzKzzDEUu3t2AbwJ/ERFPNefl3dyH7CauTl5l1nlq02i+ghJAvhER12Xx43maQj6vy/K+8ssMOO+Mk1eZdZ6aqzOipIH4aUR8tjGrmV+mNe/MvLxKMx3YlKc9NwIzJI3PAdUZWWZmXaA2A977gLsl3ZllHwcuAhZLOh14CDg5590AzKYk7H4GOBUgItZLuhC4LetdEBHrK9plZsPIeWfMdhDOO2NmHclBxMyqOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKg4iZlbFQcTMqjiImFkVBxEzq+IgYmZVHETMrIqDiJlVcRAxsyodE0QkzZJ0Xya3Oqf/JcysE3REEJE0CvgiJcHVFOCUTIRlZh2uI4IIJfPdqoh4ICKeBxZRkl2ZWYeruVHzUOotgdUxrZUknUFJwQnwnKR7hqFtnWJv4Jcj3YhhtqP1udP6++qBVOqUIDIgEXEFcAWApNsHchPZl4sdrb+w4/W5W/vbKaczA05gZWadpVOCyG3AZEmTJI0B5lKSXZlZh+uI05mI2CzpLErmu1HAlRGxop/Frtj+LesoO1p/Ycfrc1f2t2uTV5lZZ+iU0xkz61IOImZWpeuCyMvp5/GSVku6W9Kdkm7PsgmSlkpamc/js1ySFmS/75I0tbGe+Vl/paT5I9Wf3ki6UtK65m96hrKPko7K93BVLjvolLNDqY/+ni/p0fyc75Q0uzHv3Gz7fZJmNsp73c/z4sPyLL82L0SMrIjomgdl0PV+4DXAGOAnwJSRbldFf1YDe7eUfRo4J6fPAS7O6dnAtyl5macDy7N8AvBAPo/P6fEj3bdGf94KTAXu2R59BG7NusplT+zA/p4P/FUvdafkPjwWmJT79qht7efAYmBuTv8D8KGR/oy77UhkR/h5/Bzg6py+Gnh3o3xhFMuAPSXtC8wElkbE+ojYACwFZg13o/sSEbcA61uKh6SPOW+PiFgW5Vu1sLGuEdFHf/syB1gUEc9FxIPAKso+3ut+nkdZxwFLcvnmezdiui2I9Pbz+P1HqC1DIYDvSLojf9IPMDEi1ub0Y8DEnO6r7934ngxVH/fP6dbyTnRWnqJd2XP6xuD7uxewMSI2t5SPqG4LIi83b4mIqZT/Xj5T0lubM/Ov68v6GvyO0EfgcuC1wJHAWuD/jmxzhla3BZGX1c/jI+LRfF4H/BPlMPbxPEwnn9dl9b763o3vyVD18dGcbi3vKBHxeES8EBEvAl+ifM4w+P4+STnFG91SPqK6LYi8bH4eL2mcpN17poEZwD2U/vRcfZgP/HNOXw/MyysY04FNeUpwIzBD0vg8TJ6RZZ1sSPqY856SND3HC+Y11tUxegJmeg/lc4bS37mSxkqaBEymDBT3up/nUdvNwEm5fPO9GzkjPbLbxuj3bODnlNHrT4x0eyr68RrKqPtPgBU9faGc994ErAS+C0zIclFu3HQ/cDcwrbGu0yiDcquAU0e6by39vIZyCP9byjn86UPZR2Aa5Ut5P/AF8lfYHdbfr2V/7qIEjn0b9T+Rbb+PxpWlvvbz3G9uzffhH4GxI/0Z+2fvZlal205nzKzDOIiYWRUHETOr4iBiZlUcRMysioOImVVxEDGzKv8f0MmOAxJ9qwMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "stage = 350\n",
    "cond = 1E-7\n",
    "\n",
    "ghbSp = {0: []}\n",
    "\n",
    "for l in range(0, 2):\n",
    "    for r in range(0, 18):\n",
    "        ghbSp[0].append([l, r, 17, stage, cond])\n",
    "\n",
    "ghb = flopy.modflow.ModflowGhb(ml, stress_period_data=ghbSp)\n",
    "\n",
    "ghb.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RCH-Package \n",
    "\n",
    "The recharge to the top layer is separated in two zones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.axes._subplots.AxesSubplot at 0x7fdd59b01ba8>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATsAAAEICAYAAAAgMlPEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAH9dJREFUeJzt3XuYHVWZ7/Hvj4QEkVtCGA6XQaJEncAYhQBRGeSiCTCO0fOAJ+qYiGhGjXg543DRcwYeLueAl0FzFOZEyQlBh4ARlVEwBARxRpMQFAgBMSFcBAOBJFwECXT6PX+s1VDd9O69e1+yd3f9Ps9TT1evWlV7VXfnzVq1qupVRGBmNtxt1+4GmJltCw52ZlYKDnZmVgoOdmZWCg52ZlYKDnZmVgoOdlY6kq6TNKvOfR+Q9M5mt8laz8FuG8v/WP4s6U+SHpW0QNJOhe0LJL2Qt/csdxS2j5J0tqQ1kp7Nx5svaf+8/WZJH+vzmUdJeriOtp4t6bv1n21niojjI+KyZh9X0tGSbpL0lKQHmn18a4yDXXv8XUTsBLwZeAtwZp/tX46InQrLpMK2xcB7gA8CuwKTgNuAY7dBu3tRMmT+hrZBe58F5gP/1MLPsDoNmT/U4SgiHgWWkIJeVXn49C5gekTcGhFdEfFURHwrIi6ttx2STpf0iKRnJN0r6VhJxwFfBP5bsXeZe47nS/pP4DngtZJ2lXSppPX5OOdJGpHrHyDpF7m384SkK3O5JF0kaYOkpyWtknRQhfbdLOl/S1qR6/5Y0tjC9imSfiXpSUl3SDqqz7592/tS71fSdpL+h6QHc1sWStq1sP+H87aNkr400M8xIlZExOXAunp+D9ZaDnZtJGlf4HhgbY27vBNYERF/aGIb3gB8Gjg0InYGpgEPRMTPgP8FXNlP7/LDwGxgZ+BBYAHQBRxA6qlOBXqG0ucC1wNjgH2B/5PLpwJHAq8n9VDfD2wcoKkzgY8Ce+XPmpvbvw/wU+A8YCzwBeAHkvYYoL1FH8nL0cBrgZ2Ab+ZjTwQuyfvvDeyez8GGIAe79viRpGeAPwAbgLP6bP9C7qX0LD3Xl3YH1tdw/LnF/YGfDFB3KzAamChp+4h4ICLuq3L8BRGxOiK6SAHmBOBzEfFsRGwALgJm5LovAq8B9o6I5yPiPwrlOwNvBBQR90TEQOd2eUTcFRHPAv8TeH/uPf49cG1EXBsR3RGxFFiZ2/SK9kbEi32O+yHgXyJiXUT8iXRJYYakkcCJwE8i4paI2JI/t7vKz8Y6lINde7w396KOIv1jH9dn+1cjYrfC0jNzuJHUs6nmM8X9gXdXqhgRa4HPAWcDGyQtkrR3leMXe5avAbYH1heC6/8F/iJvPw0QsELSakkfzZ/7c1IP6lv5c+dJ2qXGz3wwf+a4/Pkn9QnuR9D75zRQT3hvevf2HgRGAnvmbS/tmwPtQL1P62AOdm0UEb8gDQG/WuMuNwCH5eFvM9vxbxFxBClwBHBhz6ZKuxTW/wBsAcYVAuwuEXFgPvajEfHxiNgb+AfgYkkH5G1zI+IQYCJpODvQhf2/LKzvR+oZPpE///I+/zm8OiIuqNDevv6Yz7t47C7gMVIv+qXPlbQjqXdtQ5CDXft9HXiXpEnVKkbEDcBS4IeSDpE0UtLOkj7R02MaLElvkHSMpNHA88CfeXmo9hiw/0AzmHnoeT3wNUm75Av+r5P0jnz8kwrBeTMp8HRLOlTS4ZK2J81iPs/AQ8S/lzQxB5xzgMURsRX4LvB3kqZJGiFpB6VbbWr9D+EK4POSxivdAtRznbKLNPP9bklHSBqVP7fizyKf+w6kXqdyW0bV2A5rMQe7NouIx4GFwD8Xik9T7/vsnihsOxG4FrgSeAq4C5hM6vXVYzRwAamX9Chp+NlzK8z389eNkn4zwDFmAqOAu0kBbTEvDyMPBZZL+hNwDfDZiFgH7AJ8O9d/kDQ8/MoAn3E5qRf8KLAD8BmAPFkznTRz/Dipp/dP1P63PT8f+xbgflLQPTUfezUwB/g3Ui9vMzDQ/YpHkv6zuJbUQ/wz6T8C6wDyyzut00m6GfhuRHyn3W2xocs9OzMrharBTulRpA2S7iqUvVnSMkm3S1op6bBcLklzJa2VdKekgwv7zFJ6xGmNCs8l5mtPq/I+cyWp2SdpZlZ1GCvpSOBPwMKIOCiXXQ9cFBHXSToBOC0ijsrrp5LucToc+EZEHJ7vdl9JurYUpMebDomIzZJWkK6/LCdd65gbEde14mTNrLyq9uwi4hZgU99i0gVmSHe//zGvTycFxYiIZcBukvYi3ZW/NCI2RcRm0ozicXnbLhGxLFLUXQi8t+GzMjPrY2Sd+30OWCLpq6SA+bZcvg+9b+B8OJcNVP5wP+X9kjSb9NgPr95Rh7zxAM/qDxW/v3PHdjfBBukZNj8REXtUrzk01BvsPgl8PiJ+IOn9wKWk5zZbKiLmAfMAJk/aIVYs2a/VH2lNMm3vqrcRWoe5IRb3fY54SKt3NnYWcHVe/z5wWF5/hN53uu+bywYq37efcjOzpqo32P0ReEdePwZYk9evAWbmWdkpwFP5DvslwFRJYySNIb3xYkne9rTSK3pEujn1x/WejJlZJVWHsZKuID2wPk7pbbdnAR8HvpHfDPE8+ToaaTb1BNIri54DTgaIiE2SzgVuzfXOiYieSY9Pke6MfxVwXV7MzJqqarCLiA9U2HRIP3WD9HhNf8eZT3o0p2/5SqDflzaamTWLn6Aws1JwsDOzUnCwM7NScLAzs1JwsDOzUnCwM7NScLAzs1JwsDOzUnCwM7NScLAzs1JwsDOzUnCwM7NScLAzs1JwsDOzUnCwM7NSqCtvbC4/VdLvJK2W9OVC+Zk5B+y9kqYVyo/LZWslnVEoHy9peS6/UpKz6JhZ09XSs1sAHFcskHQ0KW3ipIg4EPhqLp8IzAAOzPtcLGmEpBHAt4DjgYnAB3JdgAtJOWgPADYDpzR6UmZmfdWbN/aTwAURsSXX2ZDLpwOLImJLRNxPej37YXlZGxHrIuIFYBEwPeedOAZYnPe/DOeNNbMWqPea3euBv8nDz19IOjSXDzZv7O7AkxHR1afczKyp6s0bOxIYC0wBDgWukvTaprWqgmKS7P32qbfpZlZG9fbsHgaujmQF0A2MY/B5YzcCu+UsZcXyfkXEvIiYHBGT99h9RJ1NN7MyqjfY/Qg4GkDS64FRwBOkvLEzJI2WNB6YAKwgpVCckGdeR5EmMa7J2chuAk7Mx52F88aaWQvUmzd2PjA/347yAjArB67Vkq4C7ga6gDkRsTUf59OkZNkjgPkRsTp/xOnAIknnAb8FLm3i+ZmZAaAUo4aeyZN2iBVL9mt3M6xG0/ae1O4m2CDdEItvi4jJ7W5Hs/gJCjMrBQc7MysFBzszKwUHOzMrBQc7MysFBzszKwUHOzMrBQc7MysFBzszKwUHOzMrBQc7MysFBzszKwUHOzMrBQc7MysFBzszKwUncjAzAFbHSp5gPaMYzVs1teHj5XzSf0vqVC0FPhttfIFm3Umy87Z/lBSSxuXvJWluTnh9p6SDC3VnSVqTl1mF8kMkrcr7zM3pFc1sG9ub1/AWjmjKsSS9DXg78CbgIFJirnc05eB1qitJNoCkvwSmAg8Vio8n5Z2YQMoCdkmuO5b0OvfDSTlkz5I0Ju9zCfDxwn6v+Cwza70x2oPtGdWrTNLrJP1M0m2SfinpjTUeLoAdSPlpRgPbA481tcGDVG+SbICLgNNIJ9VjOrAwZx1bRsocthcwDVgaEZsiYjOpS3tc3rZLRCzL3duFOEm2WSeZB5waEYcAXwAurmWniPg1KZnW+rwsiYh7WtbKGtR1zU7SdOCRiLijz6hzsEmy98nrfcsrfa7zxpptO9sBbwO+X/h3PhpA0n8Fzulnn0ciYpqkA4C/IqVHBVgq6W8i4pctbnNFg44YknYEvkgawm5TETGP9D8NkyftMDQzBZkNLU9GxJv7FkbE1cDVA+z3PmBZRPwJQNJ1wFuBtgW7em49eR0wHrhD0gOkyP0bSf+FwSfJfoSXI3+x3Mzarxu4X9JJ8NIEZK1p4h4C3iFppKTtSZMTbR3GDjrYRcSqiPiLiNg/IvYnDT0PjohHSUmyZ+YfyhTgqYhYT8oXO1XSmDwxMZU0hl8PPC1pSp6FnYmTZJu1xapYzq3cxHM8wy/jpwDjgA8Bp0i6A1hNui5fi8XAfcAq4A7gjoj49xY0u2Z1JcmOiEqJrK8FTgDWAs8BJwNExCZJ5wK35nrnRETPpMenSDO+rwKuy4uZbWN/rcN7fX9DLH4iIu6njjskImIr8A9NalpTVA12EfGBKtv3L6wHMKdCvfnA/H7KV5LuwzEzaxlPaZpZy0w7+tWxcdPWXmW33bllSURs8/tpHezMrGUe39TFf/5s715lO+79wLh2tMXBzsxaJghejK3VK24DDnZm1jIBvEh3u5sBONiZWQsFsCUc7MxsmIsIXmjfW516cbAzs5YJxIt0xlvbHOzMrGW6geejM16I7mBnZi0TwIsOdmY23AXixRjR7mYADnZm1kLdiOdj+3Y3A3CwM7MWihAvuGdnZsNduqnYwc7MhrlAPN/tYayZDXNpgqIzwkxdeWMlfUXS73Ju2B9K2q2w7cycA/ZeSdMK5cflsrWSziiUj5e0PJdfKal3LjczG7K6Sdfsiku71Js3dilwUES8Cfg9cCaApInADODAvM/FkkZIGgF8i5RXdiLwgVwX4ELgoog4ANgMnNLQGZlZx4gQW2L7Xku71JU3NiKuj4iu/O0yXk6aMx1YFBFb8uuc15KSYh8GrI2IdRHxArAImJ7zThxDel89wGU4b6zZsBHACzGy19Iuzbi1+aO8nDdisHljdyelauvqU94vSbMlrZS08vGNnfGOLDOrrBuxpXv7Xks1/V0667N9V0n/LukOSaslnVxLWxoKdpK+BHQB32vkOLWKiHkRMTkiJu+xe2dMZ5tZZT1PUBSXGixg4CQ/c4C7I2ISKRnY12q51l93n1LSR4B3A8fmRDtQOT8sFco3ArtJGpl7d84bazaMRAz+cbGIuEXS/gNVAXbOl8F2Il1m6xqgPlBnz07SccBpwHsi4rnCpmuAGZJGSxoPTABWkFIoTsgzr6NIkxjX5CB5E3Bi3n8WzhtrNmz03GdXXEhpWVcWltmDPOw3gb8C/kjKS/vZiOpvCK0rbyxp9nU0sDQFV5ZFxCciYrWkq4C7SZF2Ts4fiaRPk5JljwDmR8Tq/BGnA4sknQf8FqiUk9bMhpj01pNX9OyeiIjJDRx2GnA7aXLzdaQ49MuIeHqgnerNG1sxIEXE+cD5/ZRfS0qi3bd8HWm21syGmUB0Nf/eupOBC/LIcK2k+4E3kkaRFXXGi6bMbFjqDrFl68heSxM8BBwLIGlP4A3Aumo7dcZzHGY2LNXz8s4Kl862B4iIfwXOBRZIWgUIOD0inqh2XAc7M2uhwQ9jK1w6K27/IzB1sC1xsDOzlukZxnaCzmiFmQ1LAXQ5B4WZDXshuro742knBzsza5lu4AUHOzMb7gLR1e1hrJkNcxHu2ZlZCbhnZ2alsdWzsWY23EWIF7d6GGtmw1wAWz2MNbNhL2BrqN2tABzszKyFAtG1tTN6dvXmjR0raamkNfnrmFwuSXNzDtg7JR1c2GdWrr9G0qxC+SGSVuV95uZXLZvZMBAB3d3b9Vrapd68sWcAN0bEBODG/D2kvLAT8jIbuARScCS9puVw0os6z+oJkLnOxwv7DZRow8yGmK3d6rVUUy27WK5zlKTbc3axX9TSjrryxpLyw16W14u5XqcDCyNZRkqmsxfpNcpLI2JTRGwmJdk+Lm/bJSKW5beOLsR5Y82GjUBs7d6u11KDBQzQ6ZG0G3AxKQfOgcBJtRy03j7lnhGxPq8/CuyZ1webN3afvN63vF/OG2s2xAREt3otVXfpv4NV9EHg6oh4KNffUEtTGh5A5x5ZVK3YBM4bazb0dHer10Lj2cVeD4yRdLOk2yTNrGWnemdjH5O0V0Ssz0PRnshaKW/sI6TXLBfLb87l+/ZT38yGgQjofuVsbKPZxUYCh5DyULwK+LWkZRHx+4F2qrdndw0pxyv0zvV6DTAzz8pOAZ7Kw90lwFRJY/LExFRgSd72tKQpeRZ2Js4bazaMaNDD2Bo8TIofz+bcE7cAk6rtVMutJ1cAvwbeIOlhSacAFwDvkrQGeGf+HlKqxHXAWuDbwKcAImITKUnGrXk5J5eR63wn73MfcF1Np2tmna+Oa3Y1+DFwhKSRknYk3eVxT7Wd6s0bCzmVWZ+6AcypcJz5wPx+ylcCB1Vrh5kNUYMMcNWyi0XEPZJ+BtxJej/odyKi4m0qPfwEhZm1TjDoYFctu1iu8xXgK4M5roOdmbVUbO2Mh6Ic7MysdQLU3e5GJA52ZtZCGvQwtlUc7MysdQLwMNbMysDDWDMrBW2Th0mrc7Azs9bxMNbMysLDWDMb9uRbT8ysLORhrJkNe+7ZmVlpONiZWRl0Ss+uMxI6mtnwlIexxaWaWrKL5XqHSuqSdGItTWko2En6fE5ldpekKyTtIGm8pOU5D+yVkkbluqPz92vz9v0Lxzkzl98raVojbTKzziFAW3svNVhAlZSqkkYAFwLX19qWuoOdpH2AzwCTI+IgYAQwIzfgoog4ANgMnJJ3OQXYnMsvyvWQNDHvdyDpBC/OJ2JmQ10dPbsasosBnAr8gJfz31TV6DB2JPAqSSOBHYH1wDHA4ry9b07Znlyzi4Fjc96J6cCiiNgSEfeTXs9+WIPtMrMO0U+wayi7WO5ovQ+4ZDD71T1BERGPSPoq8BDwZ1J38jbgyYjoytWKeWBfyh0bEV2SngJ2z+XLCoeumDs2/1BmA+y3j+dWzDpe9Dt0bTS72NeB0yOiO/WXalN3xMhZwqYD44Enge9TZZzdqIiYB8wDmDxphw55vNjMBtKC2djJwKIc6MYBJ0jqiogfDbRTI92jdwL3R8TjAJKuBt4O7CZpZO7dFfPA9uSUfTgPe3cFNlI516yZDXVB0++zi4jxPeuSFgA/qRbooLFrdg8BUyTtmK+9HQvcDdwE9EwF980p25Nr9kTg5zkb2TXAjDxbOx6YAKxooF1m1iEEbNfde6m6Tz/pWyV9QtInGmlLI9fslktaDPwG6AJ+Sxpi/pTUxTwvl12ad7kUuFzSWtJMy4x8nNWSriIFyi5gTkTUNkFtZp1vkD27WrKLFep+pNa6DV3lj4izSDkdi9bRz2xqRDwPnFThOOcD5zfSFjPrQH421szKYrsOGac52JlZ67RggqJeDnZm1jLCw1gzK4OA7bZ2xi2xDnZm1lLu2ZnZ8OfZWDMrg/SKJw9jzWy4c8/OzMpiu67qdbYFBzsza50AdXsYa2bDnO+zM7NyiOiYCQpnFzOz1ok0G1tcqqmWXUzShyTdKWmVpF9JmlRLUxzszKylBptwh+rZxe4H3hERfw2cS357eTUexppZ68Tg77OLiFuKqVb72f6rwrfLSG83r6rRvLG7SVos6XeS7pH0VkljJS2VtCZ/HZPrStLcnB/2TkkHF44zK9dfI2lW5U80s6Gk56biPsPYhrKL9XEKcF0tFRvt2X0D+FlEnJiTYe8IfBG4MSIukHQGcAZwOnA86ZXrE4DDSWnQDpc0lvQC0MmkF8LcJumaiNjcYNvMrN36v/Wk0exiAEg6mhTsjqilfiNJsncFjiS/dj0iXoiIJ+mdH7Zv3tiFkSwjJebZC5gGLI2ITTnALaXFWcrMbFsJtLW719IMkt4EfAeYHhEba9mnkWHseOBx4P9J+q2k70h6NbBnRKzPdR4F9szrL+WNzXryw1YqfwVJs3u6vo9v7JDXn5pZZXXMxlYjaT/gauDDEfH7WvdrJNiNBA4GLomItwDPkoasL8nZw5p2k01EzIuIyRExeY/dRzTrsGbWQuqOXkvV+tWzi/0zsDtwsaTbJa2spR2NXLN7GHg4Ipbn7xeTgt1jkvaKiPV5mLohb6+UH/YR4Kg+5Tc30C4z6xCKQF2DG7pWyy4WER8DPjbYttTds4uIR4E/SHpDLurJG1vMD9s3b+zMPCs7BXgqD3eXAFMljckzt1NzmZkNB93dvZc2aXQ29lTge3kmdh1wMimAXiXpFOBB4P257rXACcBa4Llcl4jYJOlc4NZc75yI2NRgu8ysEwSD7tm1SqN5Y28n3TLS17H91A1gToXjzAfmN9IWM+tAEW3tzRX5CQoza6lOeRGAg52ZtU4EdHXGbWIOdmbWWh7GmtmwFwFb3bMzs+Eu8DDWzMrAs7FmVgaBh7FmVgIRRFdn5FJ0sDOz1mrSa50a5RwUZtY6PbOxxaWKGhLuVHzr+UAc7MysdSKIF7t6LTVYwMAv8C2+9Xw26a3nVTnYmVnLBBBbt/Zaqu4TcQsw0MtAKr31fEC+ZmdmLfNMbFqy9MUrxvUp3qHPCzfnRURN6RCzSm83X99/9cTBzsxaJiI6Jp+Mh7FmNtRUeuv5gBoOdpJG5IQ7P8nfj5e0PM+UXJlf7Imk0fn7tXn7/oVjnJnL75U0rdE2mdmwVumt5wNqRs/us8A9he8vBC6KiAOAzaS8juSvm3P5RbkekiYCM4ADSTMwF0tyNh2zkqoh4c61pDejrwW+DXyqluM2dM1O0r7A3wLnA/9dkoBjgA/mKpcBZ5OmhqfndUjJeb6Z608HFkXEFuB+SWuBw0gna2YlU0PCnYpvPR9Ioz27rwOnAT23SO8OPBkRPTfTFHPAvjSDkrc/les7b6yZtVzdwU7Su4ENEXFbE9szIOeNNbN6NTKMfTvwHkknADsAuwDfIN3gNzL33oqzJD0zKA9LGgnsCmykzpkVM7PBaCRv7JkRsW9E7E+aYPh5RHwIuAk4MVfrmze2J5/sibl+5PIZebZ2POkRkBX1tsvMrD+tuKn4dGCRpPOA3wKX5vJLgcvzBMQmUoAkIlZLuoqUYLsLmBMRviBnZk3VlGAXETcDN+f1daTZ1L51ngdOqrD/+aQZXTOzlvATFGZWCg52ZlYKDnZmVgoOdmZWCg52ZlYKDnZmVgoOdmZWCg52ZlYKDnZmVgoOdmZWCg52ZlYKDnZmVgoOdmZWCg52ZlYKDnZmVgqN5KD4S0k3Sbpb0mpJn83lYyUtlbQmfx2TyyVpbs4Pe6ekgwvHmpXrr5E0q9JnmpnVq5GeXRfwjxExEZgCzMk5YM8AboyICcCN+XuA40mvXJ8AzCalV0TSWOAs4HDSSz/P6gmQZmbN0kgOivUR8Zu8/gwpUfY+pDywl+VqlwHvzevTgYWRLCMl5tkLmAYsjYhNEbEZWEpKlm1m1jRNuWYnaX/gLcByYM+IWJ83PQrsmdcr5YetOW+smVm9Gg52knYCfgB8LiKeLm7L2cOi0c8ofJaTZJtZXRoKdpK2JwW670XE1bn4sTw8JX/dkMsr5YetOW+sk2SbWb0amY0VKT3iPRHxL4VNxfywffPGzsyzslOAp/JwdwkwVdKYPDExNZeZmTVNI6kU3w58GFgl6fZc9kXgAuAqSacADwLvz9uuBU4A1gLPAScDRMQmSecCt+Z650TEpgbaZWb2CnUHu4j4D0AVNh/bT/0A5lQ41nxgfr1tMTOrxk9QmFkpONiZWSk42JlZKTjYmVkpONiZWSk42JlZKTjYmVkpONiZWSk42JlZKTjYmVkpONiZWSk42JlZKTjYmVkpONiZWSk42JlZKTjYmVkpdEywk3ScpHtzEu0zqu9hZla7jgh2kkYA3yIl0p4IfCAn3DYza4qOCHbAYcDaiFgXES8Ai0hJtc3MmqKRhDvN1F+i7MP7VpI0G5idv90yYq81d22DtnWKccAT7W5E/dbUs9MQP+dB67TzfU27G9BMnRLsahIR84B5AJJWRsTkNjdpmynb+UL5zrls57utdcowtuZE2WZm9eiUYHcrMEHSeEmjgBmkpNpmZk3REcPYiOiS9GlgCTACmB8Rq6vsNq/1LesoZTtfKN85l+18tyml3NVmZsNbpwxjzcxaysHOzEphyAW74fRYmaQHJK2SdLuklblsrKSlktbkr2NyuSTNzed9p6SDC8eZleuvkTSrXefTH0nzJW2QdFehrGnnKOmQ/DNcm/fVtj3D3iqc79mSHsm/59slnVDYdmZu+72SphXK+/07z5N4y3P5lXlCz2oREUNmIU1e3Ae8FhgF3AFMbHe7GjifB4Bxfcq+DJyR188ALszrJwDXAQKmAMtz+VhgXf46Jq+Pafe5Fc7nSOBg4K5WnCOwItdV3vf4Djzfs4Ev9FN3Yv4bHg2Mz3/bIwb6OweuAmbk9X8FPtnu3/FQWYZaz64Mj5VNBy7L65cB7y2UL4xkGbCbpL2AacDSiNgUEZuBpcBx27rRlUTELcCmPsVNOce8bZeIWBbpX//CwrHaosL5VjIdWBQRWyLifmAt6W+837/z3Gs9Blic9y/+7KyKoRbs+nusbJ82taUZArhe0m35UTiAPSNifV5/FNgzr1c696H4M2nWOe6T1/uWd6JP56H5/J5hO4M/392BJyOiq0+51WCoBbvh5oiIOJj0tpc5ko4sbsy9lWF9b1AZzhG4BHgd8GZgPfC19jannIZasBtWj5VFxCP56wbgh6Thy2N5eEb+uiFXr3TuQ/Fn0qxzfCSv9y3vKBHxWERsjYhu4Nuk3zMM/nw3kob2I/uUWw2GWrAbNo+VSXq1pJ171oGpwF2k8+mZbZwF/DivXwPMzDOWU4Cn8lBwCTBV0pg8PJqayzpZU84xb3ta0pR8PWtm4VgdoyewZ+8j/Z4hne8MSaMljQcmkCZc+v07z73gm4AT8/7Fn51V0+4ZksEupBm735Nmq77U7vY0cB6vJc2y3QGs7jkX0nWZG0nvRLoBGJvLRXrB6X3AKmBy4VgfJV3cXguc3O5z63OeV5CGbi+SrjGd0sxzBCaTgsd9wDfJTwV12Plens/nTlKA26tQ/0u57fdSmEmu9Hee/25W5J/D94HR7f4dD5XFj4uZWSkMtWGsmVldHOzMrBQc7MysFBzszKwUHOzMrBQc7MysFBzszKwU/j8WxC6eaRguQQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rchSp = {0: [[] for _ in range(nrow)]}\n",
    "\n",
    "rchLeft = 2E-8\n",
    "rchRight = 1E-8\n",
    "\n",
    "for r in range(0, nrow):\n",
    "    for c in range(0, ncol):\n",
    "        if c < ncol/2:\n",
    "            rchSp[0][r].append(rchLeft)\n",
    "            continue\n",
    "        \n",
    "        rchSp[0][r].append(rchRight)\n",
    "\n",
    "rch = flopy.modflow.ModflowRch(ml, rech=rchSp)\n",
    "\n",
    "rch.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LPF-Package\n",
    "The aquifer properties (really only the hydraulic conductivity) are defined with the LPF package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "hk = [3E-4, 4E-5]\n",
    "vkcb = 1E-7\n",
    "layvka = 1\n",
    "\n",
    "lpf = flopy.modflow.ModflowLpf(ml, hk=hk, vkcb=vkcb, layvka=layvka)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCG-Package and OC-Package\n",
    "Finally, we need to specify the solver we want to use (PCG with default values), and the output control (using the default values).  \n",
    "\n",
    "Then we are ready to write all MODFLOW input files and run MODFLOW."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcg = flopy.modflow.ModflowPcg(ml)\n",
    "oc = flopy.modflow.ModflowOc(ml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Head observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "observations = [\n",
    "    ['hd01', 0, 2, 0, 1.018040E+02],\n",
    "    ['hd02', 0, 3, 3, 1.281170E+02],\n",
    "    ['hd03', 0, 9, 8, 1.566780E+02],\n",
    "    ['hd04', 0, 12, 3, 1.248930E+02],\n",
    "    ['hd05', 0, 13, 5, 1.409610E+02],\n",
    "    ['hd06', 1, 3, 3, 1.265370E+02],\n",
    "    ['hd07', 1, 9, 0, 1.011120E+02],\n",
    "    ['hd08', 1, 9, 8, 1.581350E+02],\n",
    "    ['hd09', 1, 9, 17, 1.763740E+02],\n",
    "    ['hd10', 1, 17, 5, 1.420200E+02],\n",
    "]\n",
    "\n",
    "\n",
    "obs = []\n",
    "\n",
    "for o in observations:\n",
    "    ho = flopy.modflow.HeadObservation(\n",
    "        ml, \n",
    "        layer=o[1], \n",
    "        row=o[2], \n",
    "        column=o[3], \n",
    "        time_series_data=[[0,o[4]]]\n",
    "    )\n",
    "    \n",
    "    obs.append(ho)\n",
    "\n",
    "hob = flopy.modflow.ModflowHob(\n",
    "    ml, \n",
    "    iuhobsv=1051, \n",
    "    hobdry=-9999.,\n",
    "    obs_data=obs\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write input and run the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FloPy is using the following executable to run the model: /usr/local/bin/mf2005\n",
      "\n",
      "                                  MODFLOW-2005     \n",
      "    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL\n",
      "                             Version 1.11.00 8/8/2013                        \n",
      "\n",
      " Using NAME file: steady-state.nam \n",
      " Run start date and time (yyyy/mm/dd hh:mm:ss): 2018/05/29 15:40:22\n",
      "\n",
      " Solving:  Stress period:     1    Time step:     1    Ground-Water Flow Eqn.\n",
      " Run end date and time (yyyy/mm/dd hh:mm:ss): 2018/05/29 15:40:22\n",
      " Elapsed run time:  0.011 Seconds\n",
      "\n",
      "  Normal termination of simulation\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(True, [])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ml.write_input()\n",
    "ml.run_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read the heads-file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(cwdpth)\n",
    "hds = flopy.utils.HeadFile(os.path.join(workspace, name + '.hds'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD8CAYAAABQOZBmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VeX9wPHPN3tvCFNBpgxRQJYbFAS1jlILVUFFbR21WuuqP+uqs7a27klFqwKKWlRUEBxVmbJk7xEkIXsRMm6e3x/nSbjh3puEnCQk+H2/XveVc5/znOecc2/yzRnPeb5ijEEppRoq6EhvgFKqddMgopRyRYOIUsoVDSJKKVc0iCilXNEgopRypc4gIiJTRWSfiKzxKjtRRBaJyEoRWSYiQ2y5iMjTIrJFRFaLyECvZSaLyGb7muxVPkhEfrTLPC0i0tg7qZRqOvU5EnkdOPeQsieAB4wxJwJ/se8BxgI97Os64AUAEUkC7gOGAkOA+0Qk0S7zAnCt13KHrksp1YLVGUSMMd8AOYcWA3F2Oh74yU5fCLxhHIuABBFpD4wB5hljcowxucA84Fw7L84Ys8g4vd7eAC5yvVdKqWYT0sDlbgE+F5EncQLRCFveEdjtVS/NltVWnuan3C8RuQ7nCIfo6OhBvXv39lsvbaMT0zr16hBwB9K2Zzp1urYJWAdgd5oTPzt3Sqq13s70XACObZdYaz2AHZlO3S5t6q5bZVuus8xxifVfBmB7YTYAXWOTD2s5b2n79wHQKaptg9vwJ6vU+Z5SwgN/T41tf/l2AKJCuzbbOptMhbMvhDTNvvzwww9Zxpja/0BoeBC5HrjVGDNLRC4FXgPObmBb9WaMeRl4GWDw4MFm2bJlfuvddtZ9APz9ywcCtnXHZS8C8MRbv6t1nbfe9hYAT/39slrr/e7RmQC8ePeltdYDuPLldwF4/bpf1Vm3ysRZMwB455e/rvcyAL9Z8CYAb4+84rCW83bbimcB+PtJNzW4DX9e3XovANd0e6hR263N8r3O5zCw/ZvNts6mUpl9OQBByf9pkvZFZGd96jX07sxk4H07/S7OdQ6APUBnr3qdbFlt5Z38lCulWomGBpGfgDPs9Ehgs52eDUyyd2mGAfnGmL3A58BoEUm0F1RHA5/beQUiMszelZkE/LehO6OUan51ns6IyDvAmUCKiKTh3GW5FviXiIQAB7DXKYA5wDhgC7AfuArAGJMjIg8BS229B40xVRdrb8C5AxQJfGpfSqlWos4gYoyZGGDWID91DXBjgHamAlP9lC8D+tW1HUqplkl7rCqlXNEgopRyRYOIUsoVDSJKKVc0iCilXNEgopRyRYOIUsoVDSJKKVc0iCilXNEgopRyRYOIUsoVDSJKKVc0iCilXNEgopRyRYOIUsqVBuWdseW/F5ENIrJWRJ7wKr/b5pDZKCJjvMrPtWVbROQur/KuIrLYls8QkbDG2jmlVNNrUN4ZETkLJz3EAGNMX+BJW94HmAD0tcs8LyLBIhIMPIeTl6YPMNHWBXgceMoY0x3IBaa43SmlVPNpaN6Z64HHjDGlts4+W34hMN0YU2qM2Y4zTOIQ+9pijNlmjCkDpgMX2nFVRwLv2eWnoXlnlGpVGnpNpCdwmj0N+VpETrblh5t3JhnIM8ZUHFKulGolGpp3JgRIAoYBJwMzReS4RtuqALyTVx1zzDFNvTqlVD009EgkDXjfpstcAlQCKRx+3plsnFSbIYeU+2WMedkYM9gYM7hNmzoTcymlmkFDg8iHwFkAItITCAOycPLOTBCRcBHpipOgewlOqoge9k5MGM7F19l2dPgvgfG23clo3hmlWpWG5p2ZCky1t33LgMk2IKwVkZnAOqACuNEY47Ht3ISTxCoYmGqMWWtXcScwXUT+CqzAScmplGol3OSduTxA/YeBh/2Uz8FJbnVo+TYOpuFUSrUy2mNVKeWKBhGllCsaRJRSrmgQUUq5okFEKeWKBhGllCsaRJRSrmgQUUq5okFEKeWKBhGllCsaRJRSrmgQUUq5okFEKeWKBhGllCsaRJRSrmgQUUq50uDkVXbebSJiRCTFvhcRedomolotIgO96k4Wkc32NdmrfJCI/GiXedqmkVBKtRINSl4FICKdgdHALq/isTjjqvbAGZX9BVs3CWdYxaE4o5jdJyKJdpkXgGu9lvNZl1Kq5Wpo8iqAp4A7AONVdiHwhh0FfhHOSO7tgTHAPGNMjjEmF5gHnGvnxRljFtkxWt9Ak1cp1ao06JqIiFwI7DHGrDpk1uEmr+popw8tD7Te60RkmYgsy8zMbMimK6Ua2WEHERGJAv4M/KXxN6d2mndGqZanIUci3YCuwCoR2YGTcGq5iLTj8JNX7bHTh5YrpVqJww4ixpgfjTFtjTFdjDFdcE5BBhpj0nGSV02yd2mGAfnGmL04+WZGi0iivaA6GvjczisQkWH2rswkNHmVUq1KfW7xvgMsBHqJSJqITKml+hxgG7AFeAW4AcAYkwM8hJMJbynwoC3D1nnVLrMV+LRhu6KUOhLcJK+qmt/Fa9oANwaoNxUnc96h5cuAfnVth1KqZdIeq0opVzSIKKVc0SCilHJFg4hSyhUNIkopVzSIKKVc0SCilHJFg4hSyhUNIkopVzSIKKVc0SCilHJFg4hSyhUNIkopVzSIKKVc0SCilHKlQXlnRORvIrLB5pb5QEQSvObdbXPIbBSRMV7l59qyLSJyl1d5VxFZbMtniEhYY+6gUqppNTTvzDygnzHmBGATcDeAiPQBJgB97TLPi0iwiAQDz+HkpekDTLR1AR4HnjLGdAdygdpGTlNKtTANyjtjjJlrjKmwbxdxcLDlC4HpxphSY8x2nCEPh9jXFmPMNmNMGTAduNCOqzoSeM8uPw3NO6NUq9IY10Su5uC4qIebdyYZyPMKSJp3RqlWxlUQEZF7gArgrcbZnNpp3hmlWp46B2oORESuBM4HRtkBmiFwfhkClGfjpNoMsUcjmndGqVamoWk0z8XJw/sLY8x+r1mzgQkiEi4iXXESdC/BSRPRw96JCcO5+DrbBp8vgfF2+clo3hmlWpWG5p15FogF5onIShF5EcAYsxaYCawDPgNuNMZ47FHGTThJrNYDM21dgDuBP4rIFpxrJK816h4qpZpUQ/POBPxDN8Y8DDzsp3wOTnKrQ8u34dy9UUq1QtpjVSnligYRpZQrGkSUUq5oEFFKuaJBRCnligYRpZQrGkSUUq5oEFFKuaJBRCnligYRpZQrGkSUUq5oEFFKuaJBRCnlihwcT6h1GTx4sFm2bNmR3gyljloi8oMxZnBd9fRIRCnlSkPzziSJyDwR2Wx/JtpyEZGnbQ6Z1SIy0GuZybb+ZhGZ7FU+SER+tMs8bUeAV0q1Eg3NO3MXMN8Y0wOYb9+Dk1emh31dB7wATtAB7gOG4gxAdF9V4LF1rvVa7tB1KaVasAblncHJLzPNTnvnirkQeMM4FuEMwtweGAPMM8bkGGNycZJfnWvnxRljFtnxVt9A884o1ao09JpIqjFmr51OB1Lt9OHmnelopw8t90vzzijV8ri+sGqPIJrlFo/mnVGq5WloEMmwpyLYn/tseaC8M7WVd/JTrpRqJRoaRGbj5IiBmrliZgOT7F2aYUC+Pe35HBgtIon2gupo4HM7r0BEhtm7MpNohLwzj1z2T7/la75dz9+nPA/A47e+7b/Osu08dfe71e//+shsv/V+XLObv/394OD1//fiJ37rrdy0h4de+9yn/PbpPgPfA/DDjj3cO2uuT/kfPvPf/tKf0rjzC9/2q9yy8EO/5csyd3PXko8DLuftkXVv+i1fk7eNv2+YXq82/Jmx6ym/5TuK1/P+7uca3G5t8g784Le8tCKTAxXpTbLOpmLK/O+L8WRiPM23L3WmjLB5Z84EUkQkDecuy2PATJuDZidwqa0+BxiHk8h7P3AVgDEmR0QewkliBfCgMabqYu0NOHeAInFy+lbl9W2wVV+tI33HPg69WxwRE8GiT5YDsHrJVjLScuCQOpFR4Sz5cv3BtlbtIj09z7deRBiLFm2pfr98Qxo/Zeb7rDMqPJTvVm332cal29LYk5uPcEj9sFC+3uBbf/Ge3aQV+NaPDg1jwfZtPvWrl9u3k7TiPN/1hISy4KctAZaqaVXeFtJLcg79CIgIDmdR9lr/C9XD9qI15Jbt8ykPC4pgQ2HTdCRcnn457aIvoHfKIwTJwV//Mk8Wm3MfZWC7N5pkvU3B5FyGibgAiX8U8doXKjMxhY8hSc2zLw3NOwMwyk9dA9wYoJ2pwFQ/5cuAfnVtx+EozC7k3l885vMHDZCYGg9AQe5+7v/t6+CnV0pCSkz1dEFhCffcO8vnDwggMTG6ejq/6AC3/etDnz9WgKS4KJ+yvP0HuHHaf/1uY3KMn/olB7jmow/9bS4pUb71q5crK+Ha/830u10pEdF+lvBVWF7MvWte9bvuxLDYerXhz35PIW/ueAR/X0JMSEKD261NTGhPQoMTWZ3xW/q1fZqQIOcziA0/nnJPdpOss8mE9ISgREzubyHhacTui4T2wVRmNd9mNNuamlFCajyvrP5HrXUSU2J54ZM/1tlWYkI0r70ypc56SXFRvPPXyXXWq5IcE8WHt0yqf/2oKD67rP7tVy8XEc2n51532Mt5SwiL5ZWT73DVhj8xIQnc3NP/qWeTEaFH0t3sLpjGivQr6NfmaSJDO1FRWUzr68AtBMX9GVM8DZNzhRNIQjphmnlfjsogMun+X9dZ5/Kbz6lXW5MnnVqvetddPByArWlZdGwbT0RYaK31bzp7eL3arXLLsBGHVb96ub6nN2g5b5O6NE3/v1GpdX9Pjc856ukcN5mIkI6sSJ9ETFhv9pdvpWPcZUdge9xw9kWiJ0NwR0zOFZjQ3lCxDYlqvn05KoNIbKL/w/TdG/ew/cddnD5+ODNeXEBS2zgGn96rRp0V321m+bebmHLneQAsXrKVE/p3pmPHxBr1SkrKmD5jEVdd6fyRXnCac0Z24xPv0SEljisvGMrpJ3ULuI1Ltu3m1F5daBMbTU7RfqYvWlVj/g2HBJmLe/fhnTWrmbt1C3uLCgkRoWtiIpf1P5FhnToTyPjjBgScd/P3H/D0iIsDzq/SOaot1yx5jKzSfM5pdzIXdDiFF7d+SKUxXNvtArrFBOzaU6vSygMUV+QTHRJPVulPvL/7WTJKd5Mc1p5LOt1Au8guDWq3Nt0Sbq2ebhN1NokRwykoXUV4cCrRYYG/r5ZIYg8eSUvE2RA2HMpXQXAqEtJ8+9Lajt/q5bmbfS69ACBBQbz5gHPnpaLcwyuPfsy8WTUv4A0Y3o3FCw5eWN29O6dGAHnt398AEBkZxncLN/usIyI8lGduH88PG3Zz17MfkZFT6HdbNuzNpE2sE+wSoyP57/J1RIWHERUexntL1/jUv2PeZ6QXFfK7QSczolNnzj6uOxP7DeD5pYuZtmp5wM/ijsUf8fXerXgqK33m7Sis3zWA5za/z+97jufD0x7lxIQe3L36Jca0G8r4zmfyr03v1t1AAIuzPiU6xLlG9dGeVzmtzUXc2/dNxnW4iv/ueanB7dYmLvwEtuT8jUV7xvLNriF8n3Ymm7IfZG/R+5R78ptknU3GlNZ4K0HRSPgIQDAHPmu2zTgqj0Tyswr408j7fcpFYNd6p4NsbEIUj7x+LQ/87nWyMvKZeINznTgoKIjQsODqZYKCal70mzdvDVOuco4+gvxcFDXGEB0Zxq0Tz2TL7kwenjqXoX2P5bKxNZ+o9r6gKiJER4Rz5WmDAPh45QafdtdlZfKPMeMAGNqpM5fMeJtbho1gSMdOjH1rGpMHDPRZBuCklI68tGEhdyz5iDEde3HBsX05uc0xfusGUolhQEJ3AE5p05/ntrzPGW1PBODVbfW7TeyXCJWmkiAJosRTyPHxTl73LtHHU1Z5oOHt1mJN5i0kRY5gYLs3CQtOAZw7M3uLPmBN5q2c1M7/P6CWyBQ8hET4Oy0PwhQ9g0Q0z2NoR2UQiUmI5rq/XeH3zoe3uMRoHn3jtzzxx7e5/7f/ZvT4k9mwchftj0murpOaGsf7Hy7j1BE9mb9gHTEx4bw9fSEREaE17s5cdPtrGGPYl1vEhX96tcZ6Fq/d6RNEOibGMe3b5ZzTtzvvLv2Rfh1TqU1YcDDbcnM4LjGJVRnpRIQ4X11IUBDBQYEPKCd2G8jEbgPJKCnkk13reGTlF2QeKGZc5+PJL6/fH2pYUAjLczcxMLEn3+xbCcCCjOXEhEQSIsF1LB1Y37ihvJ/2LGe2HU/vuJP5NnM2/eNHsLVoNQlhbRvcbm1KK/ZybHzNC81hwSkcG38tewvfa5J1NpnKXCpzrvAzQ6Bia7NtxlE5KNGfxz3MI3PuqXX5Z+/7gJseOHg94KuPVrD0642ktIvn0t+eRXRsBADZ2UU8+/wX7N6dzamn9uTiiwbz+uv/o9IYJl9xCklJzu3g/KISRIQr/vIf3njA96JWfExkjfdZhcU8/snXbMnIpn+ndtx+3unERoQD8L+NOzitV5ca9Ren7ea2eZ8SGhRMpTE8O+4C+rdNJadkP68sX8adp9T/Auquolxm71zL3LQNzB5zTZ31txSm8dSmmezZn0mvuGO4rdcEpm77hKyyfK497gJ6xR3ekY235TkLWJLzOdmlGXhMBfGhyfSNH8ZpbS4iPDiy7gYO04r0q0mMGEb7mIsJD3EenSityCS9+ANySxZzYrvXGn2dTaVy3ylI4kv47acASGhfV+3Xd1CiozKIHEk5Bfv99gtpLLklJSRGNv4f189FuaeAnfkvk1UynzJPDmAIC04mJXIUXRJ+S0hQw/u9NLfKnCkEJTVd0PvZB5Gd63azcPYysvY4HWNTOiYx/MKTOfb4g4/q7NycweIF68jKcC6opaTGM3RkH47tUfPUYtHirSxcuJnMbOciaZvkWIYP78GwoTWvgH+3ahvfrNhGZm6RUy8xhjMGdmPECV39buOWjGy+XL+VfQXFALSNi2Zkn250a5vst/7m7Gy+2L6FjCKn/dSYGM45rjvdk/zXr14uP5P5P20mfb+z/e2iYjm7Y0+6x6XUupy3xdnrWJi1huwy57NKDotneEo/hib3qXcb/uw7sJv1BUspKHcu8saFJnN83Mm0jQh8x0k1j591EHnn0Q/4+t3vOfPXp9C2s/MHtm93Nl/N+I4zfjWCiXdfzIwXF/DNnNWccd4A2rR3ekdm7s3j609Wcfq4E/j170YC8K9n5pKens/oc/rRpo3zXyozs5C589bQoUMCv7/RubD1tzfn81NmAeNO6UNqklMvI6eQOd+to2PbeP50+cga2/jKV0v4bPUmxg7oRbt4p356fiGfrtrIuSf05Nozh9So//zSxczZvJHzevaiQ2wcAHsLC/l40wbG9ejFDScP9ftZvLDuO+bsXs/5x/ShfZRdbn8BH+9ax7jOx3N9n1Pq/Kyf2TSL9APZnNPuZNqE28+qNI956UvpEJnCjT0uqbMNf77eN4sf877nhIRTiA91Alp+eRar876jf8IIzmj7ywa1W5es/V+TtX8+pR6ny314cFvaRI0iOeqMJllfUzIVW+DAfExlBgASlAoRo5CQ7q7brm8QwRjTKl+DBg0ygUzqfqMpKy3zKa8orzCTetxkjDHmqpGPmbLScr91rh71WPX7yya9EHA9l096sXr64ttfDVjvkjte8ykb88RrprS8wqe8vMJjzv3bVJ/yM19/1ZRW+Knv8ZgzXw+87jM/fi7gcmd9/FzA5bxNWvjXwPMWBZ5XlyfXX2/KPX6+p8oK8/f1NzS43dpsyHrArEy/1qQXfWLySn4weSU/mPSiT8zK9GvNxqyHmmSdTaWy8AXjyfyFqSx8yVTun+28Cl+yZYF/b+sLWGbq8bd4VPYTCQ4NqT6N8ZaZlk1wiHM3ISQkmOyMAp86Wen5BAcfvOMQERHKunW+oxOsX/8TEREHe6VGhIfy45affOqt2brXb+/V0OBg9hUU+ZRnFBQS4uduS0hQUPVpjLf0Iv/1q9cjQWSU+PZVSS8pIETq9/VHBIexLn+HT/n6gp1EBIXVqw1/giWYggrf76mgPJugem7b4cop+R8DUl8mNXoc8REDiY8YSGr0OAakvkx2yTdNss6mYkreQ5LfRWKuQyIvcF4x1yHJszAls5ptO47KW7w3/PMq7hr9EB26tyOlo3M6k/1TDns27+X3z10LwG//7xf8+cpX6HBsMintnA5P2RkF7NmRxY33H7xrc8efzuNfT8+lqPgAbVKc046srCKio8O58/bzquv9ZcoYnnhjPoX7S2lrT2cyc4uIiQzjvmvG+GzjXeefyTWvzeKY5ITq05mMgiJ2ZuVy74U+zzbyl9PPYtKH73FsfALtYpw7QhnFRezIy+PBM33rV7l34Ggmff02x8Yk0j7SnjaVFLKjKJcHB9WvH8Gfek/g6U2zKK4oISXc+ayyyvKJDo7k9t6/qVcb/pzXYQr/3vYAyWHtiAt1vqeCihyyS/fyi47unvcJJFgiyT+wkviIE2uU55euIlha2wXrEPBkQMgh14886UDDb70frqPymgg4p2kblmypcWG195DuNfqOGGPYuGo32fbCanJqPL0GdPbbvyQnt5isLOc/ekpKLEkButZn5xfXuLCaHB/4SVljDD/uTifDHpGkxsXQv3O7gP1bjDGsykgn3R6RtIuJYUBq4Po1lsv5qfqIJDUylgFJHepc7lC5ZYVkleYBkBKe4OoJXu9tSyvZTEG58z3FhSbRKbLHYW9bfRWWrmVj9gNUVBYQHtIOgFJPBiESS6/kB4gNP75J1tsUTOm3mIL7IfhYCLY3AzwZ4NmJxN2PhNfvua9AftYXVpWqS5kni9IK52JkeEhqde/V1sYYA+WrodIOQhTUDkJPaJQg3CzJq0TkVhFZKyJrROQdEYkQka4istjmkZkhImG2brh9v8XO7+LVzt22fKOI+B77H6Z7zn+kzjp/uaZ+3Zuvu/7f9ap3xX3/qVe9Kte/7n/EsUCunv3+YdWvMuWbho88VuX6ZU+6bsOfadv/2iTt1kdYcAqx4X2JDe/bagMIOI9MSNgAJGKM8wob0GRHcQHV5+qrvxfOqOzbgUj7fiZwpf05wZa9CFxvp28AXrTTE4AZdroPsAoIB7oCW4HgutZf292Z7L05dV55zt5XUGedprSvoOjw6hcdXv3q5fYXNmi55lBQVvf31NgW77mo2dfZVDw51zRp+zTT3ZkQIFKcsdmigL3ASKDqIYRDc9JU5ap5Dxhlx1W9EJhujCk1xmzHGVqxZieJw5TU7uBTtwXZhRRk+96dSGpz8Hy+ILeYgtziWtvMySli06Z0Nm1KJyfH9y5Jlay8YtbvyGD9jgyy8gK3WfUEL0BecQl5xSW1rr9N9MH6uSUl5JbUXr96uciDo7Tllu4nt3R/vZbzJ6e0gE2Fu9lUuJucUt87W4crNvTg97S/opD9Ff6feG5MQzp80OTraC4SV/cRd3No8N0ZY8weEXkS2AWUAHOBH4A8Y0yFreadR6Y694wxpkJE8oFkW77Iq+mAuWdE5DqczHocc0zg5zX2bs/gtbvfYtWXa4lOcP74inKLOHFkP6557HLadWlL+u4c/v3kHFYt2kpMnHNVvjC/hBOHd+Pq28eR2ikJgE2b0vnXM3MpLi6t7myWlVVIZFQYt9w8hp49nItz63dk8Lc35lNUUlbd2WxfbiFREWHcOWkUvbvU7AWblpPPU599y5Jtu6ufmSk4UMrQ4zrzx7Gn0jExvkb93fn5PPH9NyxK201seIStf4DhnTtz5ymn0ymuZv3q5YryeGL1Ahbt20lcqF2u/ADD2h7LnQNG0im67mEINxXu5plNsyj2lFR3NssqzScqOJybe46nR2zDepfmlGUwd++bbCtaQ0Sw8z2VeIrpFtOPMe0nkdhED+GB87xMqcdeEwlOrX6OpjWR4IPbbCpznbKgxEDVm059Dlf8vYBEYAHQBggFPgQuB7Z41ekMrLHTa4BOXvO2AinAs8DlXuWvAePrWn9tpzM3Db3LfDn9W+PxeKrLPB6P+Wrm9+amYXcbY4z5wyVPm68+XuFT55s5q8wffvlMddnV17xq1q3f47OO9Rt+MlOuPdjJa8I9r5s1W3/yqbd2614z8Z5pPuW/fvZtM2fVBuPxVHqtv9J8tnqjmfDc2z71L5r+H/PRxvXGU+lVv7LSfLJpg7l4+lsBP4uL5041H+1c67vcrnXmkrm+ndr8uWbJ42Z9/g6f8g35O821Sx6vVxv+PL/pDrMq93/GU+n1HVR6zOrc78wLm+9ocLu1yT/wo1m651dmYdq5Zvneq8zyvVeZhWljzdI9403BgTVNss6mUlm+y3hybzaejKHGs+8c55U+xHhybzaV5btdt08znM6cDWw3xmQaY8qB94FTcFJnVh3heOeRqc49Y+fHA9kEzknTYIW5xZz561MI8uqEFRQUxBm/Gk6hPRUpKijhjPNO9Klz2tgTKMo/eMhfVlZB717tfdbRu1d7ysoqDtYr99Cnazufen2Oa0dpeYVPeUHJAcae0KvGeCVBQcKY/j3J3+/7iH5+aSnn9+xdYwyTIBHG9ehFXmngR/rzyw5w/jF9fJfrfDx5ZfUbCqCsspxesb5Hfr3ijqGs0nff6qvEU8QJCafW6FgWJEH0TxjB/orAp4xurM+6m57J/8ewjp9yUrupnNRuKsM6zqFn8r2sy7qr7gZaEJN/KxI+GmnzPUFt5hLUZi7SdiEScS4m/9a6G2gkbjqb7QKGiUgUzunMKGAZ8CUwHpiOb06aycBCO3+BMcaIyGzgbRH5B9ABJ6n3EhfbRfeTuvCv61/mnMln0qaT04kpMy2bedO+oseg4wDo1qcDz/zlfc6+eFB1Z7Os9Hy++OAHuvc9eDY1dMhx3H3Pu36fnRky5OADeCNO6MqtT33AuBF9qjub7cspZM736zjFzwN4x3doy4MfzufCgX1Ijbedx/KL+O/ydfTxM7ZI3zZt+L8vv+CXvfvQLsZ2GisqZNaGdfRrG/iwv09iKvcu+5RLupxAuyi73P5C3t+xmn5JvkHPnyFJx3PPj69wTurgms/OZCxjSFLD+1Wqho5YAAARTUlEQVR0iOzKf9Ne4qTEM4m3nc3yy7NZkfsVHSObZni/SlNKbFh/n/K48BOoNGVNss4mU1mARJ5Xo0gkCCLGYgr95/RpCq76iYjIA8CvgQpgBXANzvWM6UCSLbvcGFMqIhHAm8BJOAnCJxhjttl27gGutu3cYoypM/dMbf1EKsor+PS1BSycvbRGZ7MRFw5h7JSRBIcEU1Hu4fN3l7Bofs2neIef3Zcxvzq5uns8wJKl2/ju+801OpudOqIHJ598XI31Lly9na9XbD3kKd7uDO/fxWcbyz0e3l+2hgXrtlV3f28bF8OoPt24ZHA/QoKDfOrPXLeGL7ZtqdHZ7JzjunNp3/4Bu76XV3p4d9tK5u3ZREaJ7dQWGcM5HXtx6XEn1tpl3tvS7PV8n72GrFL7WYXHMyKlPycn9a7X8v54TAXLcuazvmBJjc5mfeKGMCjpbIJdDHgUyKbsv1JSsZN20RcRHuIE69KKDNKLPyQypAs9k2sfh6Ylqcz7A0gCEnkxBNt/CJ50TMkHYAoISnAXSLSzmVIBZO//hqyS+TU6m6VEnU1y5GlHeMsOjzHlUPIe5sAXUGmTgAW1dYZMjBxfM6FVA9Q3iByVz86UlZbz2WsLWPjRIUcivziZc6eMJDQslLLSCua+5xyJVD2Il5wax7BRzpFIaJjz0RQUlDB9xiK+X7iFHHuEkZgYzSnDezBxwnBi7Qho+UUlvDFnKf9bsY3sfOfWblJcFKcP7Mbk84YQFx1RcxsrKpi1dK0dT+TgkchIeyQSFlLzv3BpRQUz161h/ratpBfbI5HoGM4+rhuX9u1PWLD//9qlngre3baSL37aXKPb+zkde/KrricGXM5bQXkxM3YtYGH2GnLLnDYSQ2MZntKPCceMIja0YYMwVVSWsyznCzYULK1+EC8uJInj44cwKHEUIUG1p91oqOSo00mOcp9K40gTCYWoiUhUoPxyzbQdR+ORyEO//gfxybHONRE7nkjmbueaSGFuEfe8cyuP3Pwf4hKjOfviQbRp71wTydzrXBMpzN/P3f90hjj80x3vMGhgV8aM6V/9vExObjFz5/7ID8t38LfHJwBOqoghfY/h/FP7Vj8vk51fzCffrWPp2l08c3vNsTH++PbHJERFctHAPrRLsNcq8gr5cPk6CkoO8OTEmue6N835iMTISH55fF/a22sie4sKmbV+LfkHDvD02PP9fha///59EsMiuaTrCTXGE3l/+2ryyg7UK2XEHSufZ2BSL8a0G1L9vExuWSFz05ewPHcTjw+4vs42/Hln55NEhcQyMPGsGuOJLM/9kpKKIiYce1uD2q1NuSePnfmvkFWygDKPkyUuLCiZlKhRHBt/HaHB/m+Vt0TGlEHJu5gDC8COJ0JQKhIxyh6JNPwJa/i5jydixwypbZ73mCGH8p53xZUvBaznPe+Xdwa+Xepvnr8xQ2qbV+uYIbXMq23MkPqOJ3LloocbNK8utY0Z0lTjiSzfO9nsyHvJlFZkVpeVVmSaHXkvm+V7r2qSdTYVT+7NxpN3n6ksXWkqK9KdV+lK48m7z3hy/+C6fX7O44nEJcfy5fTv8Hg81WUej4evZnxHvH2cPzYhmq8+XonHU+lVp5KvP1lJnNcTuqlt43hn+qIavVRzcop4Z/oi2qUe/K/VPjmOaZ8sqdFLNSuvmGmfLKF9SpzPNiZERTBn1cYa+WA8lZV8unojiVG+j6QnRkTy0aYNPvU/3rSBpFrGXE0Mi+SjXWt9l9u1jqTw+p2GtI1IZPqu+TV6qeaUFjB913xSI5Lq1YY/USExrM77lkpz8HuqNB5W531LVEjTjHVaNdq79/MyVaO9l1a46lnQ/MrXERR/v/O8THCq8wobQFD8/VDe8ETrh+uoPJ1J37GP1+5+i5UL1hAVFwkiFOcVc9Ko/lz7+OW0PaYNGWk5TH3yU1Yt3EJUTAQiUFRwgBNHdGfKHeNo28Hp+VdUdIC3py/i++83k5dXjAESEqI4ZURPfjNxGDH2Wkdh8QGmfbKUb1ZsJadgP2BIjI3ijIHdufL8IcREhdfYxj25To/VxVt3Ex0ehohQeKCUYd0688exp9EhoWbgSSvI54nv/sfCtF3EhDltFZaWMqLzMdx56ul0jPUNVABpxXk8sepLFu7bQUxoOAIUlpcyvG0X7howkg7RdR++F5WXMH3XF3yfvYa8MieYJoTGMCKlHxOPPZvokIaNw5Fbto/P977JtqIfCQ+OQhDbY7U/57afREJY4/ciXZF+NUmRw2gXfRSM9p79KyRqMkSMReydLGM8cOAzzP5pBCXPdNW+3p2xCmwGurikwP/ZCvOczmWxCU03Sntt8mznsoSoiDpq2voHnOdmEiIO7483r9QuF97yBt+pem6mqY5Aqvgd7T0oiZSos1vdaO+mIg1T9CSULQKJAQQqCyB8OBJ7BxLcwVX7GkSU+hkxlc5gURJU97NQ9fWzvsWr1M9NYwaPw3VUXlhVSjUfDSJKKVc0iCilXDkqg8htZ93HbWfdV2udOy57kTsue7HOtm697S1uve2tOuv97tGZ/O7R+t1Su/Lld7ny5XfrVbfKxFkzmDhrxmEtA/CbBW/ymwVvHvZy3m5b8Sy3rXjWVRv+vLr1Xl7dem+jt1ub5XuvYPneK5p1nU2lMvtyKrMvP9KbcXQGEaVU89EgopRyRYOIUsoVt3lnEkTkPRHZICLrRWS4iCSJyDwR2Wx/Jtq6IiJP2/wyq0VkoFc7k239zSIy2e1OKaWaj9sjkX8BnxljegMDgPXAXcB8Y0wPYL59DzAWZ+jDHjgjtr8AICJJwH3AUJxUEfdVBR6lVMvX4CAiIvHA6Tijs2OMKTPG5FEzv8yheWfesE8ZL8IZ0Lk9MAaYZ4zJMcbkAvOA+mWaVkodcW6ORLoCmcC/RWSFiLwqItFAqjFmr62TDlSNOlydd8aqyi8TqNyHiFwnIstEZFlmZqaLTVdKNRY3QSQEGAi8YIw5CSjm4KkLAHZgk0Z7ws8Y87IxZrAxZnCbNq0v2ZBSRyM3QSQNSDPGLLbv38MJKhn2NAX7044gGzC/TKPnnVFKNZ8GBxFjTDqwW0R62aJRwDoO5pcB37wzk+xdmmFAvj3t+RwYLSKJ9oLqaFumlGoF3A4F8HvgLXFGhN0GXIUTmGaKyBRgJ3CprTsHGIeTsHu/rYsxJkdEHgKW2noPGmNyXG6XUqqZuAoixpiVgL9BS0b5qWuAGwO0MxWY6mZblFJHhvZYVUq5okFEKeWKBhGllCsaRJRSrmgQUUq5okFEKeWKBhGllCsaRJRSrmgQUUq5okFEKeWKBhGllCsaRJRSrmgQUUq5okFEKeWKBhGllCuug4iIBNuBmj+277uKyGKbX2aGHbAIEQm377fY+V282rjblm8UkTFut0kp1Xwa40jkDzj5Zqo8DjxljOkO5AJTbPkUINeWP2XrISJ9gAlAX5xUEc+LSHAjbJdSqhm4zYDXCTgPeNW+F2AkzqDN4Jt3piofzXvAKFv/QmC6MabUGLMdZ/jEIW62SynVfNweifwTuAOotO+TgTxjTIV9751Dpjq/jJ2fb+tr3hmlWjE3GfDOB/YZY35oxO2pleadUarlcTNQ8ynAL0RkHBABxOHk5k0QkRB7tOGdQ6Yqv0yaiIQA8UA2mndGqVbNTd6Zu40xnYwxXXAujC4wxlwGfAmMt9UOzTtTlY9mvK1vbPkEe/emK07C7yUN3S6lVPNym3fGnzuB6SLyV2AFNuG3/fmmiGwBcnACD8aYtSIyEyfxVQVwozHG0wTbpZRqAo0SRIwxXwFf2elt+Lm7Yow5APwqwPIPAw83xrYopZqX9lhVSrmiQUQp5YoGEaWUKxpElFKuaBBRSrmiQUQp5YoGEaWUKxpElFKuaBBRSrmiQUQp5YoGEaWUKxpElFKuaBBRSrmiQUQp5YoGEaWUK27GWO0sIl+KyDoRWSsif7DlSSIyT0Q225+JtlxE5GmbX2a1iAz0amuyrb9ZRCYHWqdSquVxcyRSAdxmjOkDDANutDlk7gLmG2N6APPte4CxOEMf9gCuA14AJ+gA9wFDcQYzuq8q8CilWj43Y6zuNcYst9OFOAmsOlIzv8yheWfeMI5FOAM6twfGAPOMMTnGmFxgHk4SK6VUK9Ao10RsSsyTgMVAqjFmr52VDqTa6UD5Zeqdd0Yp1fI0Ri7eGGAWcIsxpsB7nh3N3bhdh9e6NHmVUi2M2zSaoTgB5C1jzPu2OMOepmB/7rPlgfLL1DvvjCavUqrlcXN3RnDSQKw3xvzDa5Z3fplD885MsndphgH59rTnc2C0iCTaC6qjbZlSqhVwmwHvCuBHEVlpy/4MPAbMFJEpwE7gUjtvDjAOJ2H3fuAqAGNMjog8BCy19R40xuS42C6lVDNqcBAxxnwLSIDZo/zUN8CNAdqaCkxt6LYopY4c7bGqlHJFg4hSyhUNIkopVzSIKKVc0SCilHJFg4hSyhUNIkopVzSIKKVc0SCilHJFg4hSyhUNIkopVzSIKKVc0SCilHJFg4hSyhUNIkopVzSIKKVcaTFBRETOFZGNNrnVXXUvoZRqCVpEEBGRYOA5nARXfYCJNhGWUqqFaxFBBCfz3RZjzDZjTBkwHSfZlVKqhRNn6NMjvBEi44FzjTHX2PdXAEONMTcdUu86nBScAP2ANc26oUdWCpB1pDeimf3c9rml7e+xxpg6c7O4Ge292RljXgZeBhCRZcaYwUd4k5rNz21/4ee3z611f1vK6Uy9E1gppVqWlhJElgI9RKSriIQBE3CSXSmlWrgWcTpjjKkQkZtwMt8FA1ONMWvrWOzlpt+yFuXntr/w89vnVrm/LeLCqlKq9WoppzNKqVZKg4hSypVWF0SOpu7xIrJDRH4UkZUissyWJYnIPBHZbH8m2nIRkaftfq8WkYFe7Uy29TeLyOQjtT/+iMhUEdknImu8yhptH0VkkP0Mt9hlA+WHbhYB9vd+Edljv+eVIjLOa97ddts3isgYr3K/v+f25sNiWz7D3og4sowxreaFc9F1K3AcEAasAvoc6e1ysT87gJRDyp4A7rLTdwGP2+lxwKc4SdSHAYtteRKwzf5MtNOJR3rfvPbndGAgsKYp9hFYYuuKXXZsC9zf+4E/+anbx/4OhwNd7e92cG2/58BMYIKdfhG4/kh/x63tSOTn0D3+QmCanZ4GXORV/oZxLAISRKQ9MAaYZ4zJMcbkAvOAc5t7owMxxnwD5BxS3Cj7aOfFGWMWGeev6g2vto6IAPsbyIXAdGNMqTFmO7AF53fc7++5PcoaCbxnl/f+7I6Y1hZEOgK7vd6n2bLWygBzReQH26UfINUYs9dOpwOpdjrQvrfGz6Sx9rGjnT60vCW6yZ6iTa06fePw9zcZyDPGVBxSfkS1tiBytDnVGDMQ5+nlG0XkdO+Z9r/rUX0P/uewj8ALQDfgRGAv8PcjuzmNq7UFkaOqe7wxZo/9uQ/4AOcwNsMepmN/7rPVA+17a/xMGmsf99jpQ8tbFGNMhjHGY4ypBF7B+Z7h8Pc3G+cUL+SQ8iOqtQWRo6Z7vIhEi0hs1TQwGuep5NlA1d2HycB/7fRsYJK9gzEMyLenBJ8Do0Uk0R4mj7ZlLVmj7KOdVyAiw+z1gklebbUYVQHTupiDT5/PBiaISLiIdAV64Fwo9vt7bo/avgTG2+W9P7sj50hf2W3A1e9xwCacq9f3HOntcbEfx+FcdV8FrK3aF5zz3vnAZuALIMmWC87ATVuBH4HBXm1djXNRbgtw1ZHet0P28x2cQ/hynHP4KY25j8BgnD/KrcCz2F7YLWx/37T7sxoncLT3qn+P3faNeN1ZCvR7bn9vltjP4V0g/Eh/x9rtXSnlSms7nVFKtTAaRJRSrmgQUUq5okFEKeWKBhGllCsaRJRSrmgQUUq58v/Beq3h84QpaAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "h = hds.get_data(kstpkper=(0, 0))\n",
    "x = y = np.linspace(0, delr * nrow, nrow)\n",
    "\n",
    "c = plt.contour(x, y, h[0], np.arange(100, 300, 10))\n",
    "plt.clabel(c, fmt='%1.2f')\n",
    "plt.axis('scaled')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}