examples/Notebooks/flopy3_Zaidel_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FloPy\n",
    "\n",
    "## MODFLOW-USG $-$ Discontinuous water table configuration over a stairway impervious base\n",
    "\n",
    "One of the most challenging numerical cases for MODFLOW arises from drying-rewetting problems often associated with abrupt changes in the elevations of impervious base of a thin unconfined aquifer. This problem simulates a discontinuous water table configuration over a stairway impervious base and flow between constant-head boundaries in column 1 and 200. This problem is based on \n",
    "\n",
    "[Zaidel, J. (2013), Discontinuous Steady-State Analytical Solutions of the Boussinesq Equation and Their Numerical Representation by Modflow. Groundwater, 51: 952–959. doi: 10.1111/gwat.12019](http://onlinelibrary.wiley.com/doi/10.1111/gwat.12019/abstract)\n",
    "\n",
    "The model consistes of a grid of 200 columns, 1 row, and 1 layer; a bottom altitude of ranging from 20 to 0 m; constant heads of 23 and 5 m in column 1 and 200, respectively; and a horizontal hydraulic conductivity of $1x10^{-4}$ m/d.  The discretization is 5 m in the row direction for all cells.\n",
    "\n",
    "In this example results from MODFLOW-USG will be evaluated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flopy is installed in /Users/jdhughes/Documents/Development/flopy_git/flopy_fork/flopy\n",
      "3.7.3 | packaged by conda-forge | (default, Jul  1 2019, 14:38:56) \n",
      "[Clang 4.0.1 (tags/RELEASE_401/final)]\n",
      "numpy version: 1.17.3\n",
      "matplotlib version: 3.1.1\n",
      "flopy version: 3.3.1\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "import platform\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": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set name of MODFLOW exe\n",
    "#  assumes executable is in users path statement\n",
    "exe_name = 'mfusg'\n",
    "if platform.system() == 'Windows':\n",
    "    exe_name += '.exe'\n",
    "mfexe = exe_name\n",
    "\n",
    "modelpth = os.path.join('data')\n",
    "modelname = 'zaidel'\n",
    "\n",
    "#make sure modelpth directory exists\n",
    "if not os.path.exists(modelpth):\n",
    "    os.makedirs(modelpth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model dimensions\n",
    "nlay, nrow, ncol = 1, 1, 200\n",
    "delr = 50.\n",
    "delc = 1.\n",
    "\n",
    "# boundary heads\n",
    "h1 = 23.\n",
    "h2 = 5.\n",
    "\n",
    "# cell centroid locations\n",
    "x = np.arange(0., float(ncol)*delr, delr) + delr / 2.\n",
    "\n",
    "# ibound\n",
    "ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
    "ibound[:, :, 0] = -1\n",
    "ibound[:, :, -1] = -1\n",
    "\n",
    "# bottom of the model\n",
    "botm = 25 * np.ones((nlay + 1, nrow, ncol), dtype=np.float)\n",
    "base = 20.\n",
    "for j in range(ncol):\n",
    "    botm[1, :, j] = base\n",
    "    #if j > 0 and j % 40 == 0:\n",
    "    if j+1 in [40,80,120,160]:\n",
    "        base -= 5\n",
    "\n",
    "# starting heads\n",
    "strt = h1 * np.ones((nlay, nrow, ncol), dtype=np.float)\n",
    "strt[:, :, -1] = h2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create and run the MODFLOW-USG model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FloPy is using the following  executable to run the model: /Users/jdhughes/.local/bin/mfusg\n",
      "\n",
      "                                  MODFLOW-USG      \n",
      "    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUNDWATER FLOW MODEL\n",
      "                             Version 1.5.00 02/27/2019                       \n",
      "\n",
      " Using NAME file: zaidel.nam \n",
      " Run start date and time (yyyy/mm/dd hh:mm:ss): 2020/06/26 10:16:30\n",
      "\n",
      " Solving:  Stress period:     1    Time step:     1    Groundwater Flow Eqn.\n",
      " Run end date and time (yyyy/mm/dd hh:mm:ss): 2020/06/26 10:16:30\n",
      " Elapsed run time:  0.012 Seconds\n",
      "\n",
      "  Normal termination of simulation\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(True, [])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#make the flopy model\n",
    "mf = flopy.modflow.Modflow(modelname=modelname, exe_name=mfexe, model_ws=modelpth)\n",
    "dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol,\n",
    "                               delr=delr, delc=delc, \n",
    "                               top=botm[0, :, :], botm=botm[1:, :, :], \n",
    "                               perlen=1, nstp=1, steady=True)\n",
    "bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)\n",
    "lpf = flopy.modflow.ModflowLpf(mf, hk=0.0001, laytyp=4)\n",
    "oc = flopy.modflow.ModflowOc(mf, \n",
    "                             stress_period_data={(0,0): ['print budget', 'print head', \n",
    "                                                         'save head', 'save budget']})\n",
    "sms = flopy.modflow.ModflowSms(mf, nonlinmeth=1, linmeth=1,\n",
    "                               numtrack=50, btol=1.1, breduc=0.70, reslim = 0.0,\n",
    "                               theta=0.85, akappa=0.0001, gamma=0., amomentum=0.1,\n",
    "                               iacl=2, norder=0, level=5, north=7, iredsys=0, rrctol=0.,\n",
    "                               idroptol=1, epsrn=1.e-5,\n",
    "                               mxiter=500, hclose=1.e-3, hiclose=1.e-3, iter1=50)\n",
    "mf.write_input()\n",
    "\n",
    "# remove any existing head files\n",
    "try:\n",
    "    os.remove(os.path.join(model_ws, '{0}.hds'.format(modelname)))\n",
    "except:\n",
    "    pass\n",
    "\n",
    "# run the model\n",
    "mf.run_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read the simulated MODFLOW-USG model results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the mfusg headfile object\n",
    "headfile = os.path.join(modelpth, '{0}.hds'.format(modelname))\n",
    "headobj = flopy.utils.HeadFile(headfile)\n",
    "times = headobj.get_times()\n",
    "mfusghead = headobj.get_data(totim=times[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot MODFLOW-USG results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAF3CAYAAABaGRmGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8ddnQsISQBIIyhZABBFZoqKyyCZolbr2q3VpBcUWvz9F1KrfVq3V79e6Vdq6VutaW1u1WuuKK6KoKSpSEBCpqMgOgsi+ZDm/P84kDDEJM0lmbu7M+/l43MfM3Lkz95Mr5p1z77nnmHMOERERCZdI0AWIiIhI4hTgIiIiIaQAFxERCSEFuIiISAgpwEVEREJIAS4iIhJCSQ9wM+tiZtPNbKGZLTCzS6LrrzezFWY2J7qMTXYtIiIi6cKSfR+4mXUAOjjnZptZK+Aj4BTgh8AW59yUpBYgIiKShpokewfOuVXAqujzzWa2EOiU7P2KiIiks5ReAzezbsAhwPvRVZPM7GMze9jM8lJZi4iISJgl/RR65Y7MWgJvAzc6554xs32BdYADbsCfZp9QzecmAhMBcnNzD+vdu3dK6hUREQnaRx99tM45V1DdeykJcDPLBl4EXnXO/a6a97sBLzrn+tb2PQMHDnSzZs1KSo0iIiKNjZl95JwbWN17qeiFbsBDwMLY8I52bqtwKjA/2bWIiIiki6R3YgOGAucA88xsTnTd1cBZZlaEP4W+BLggBbWIiIikhVT0Qn8XsGremprsfYuIiKQrjcQmIiISQgpwERGREFKAi4iIhJACXEREJIQU4CIiUm9mxjnnnFP5urS0lIKCAk444YTKdc8++yz9+/end+/e9OvXj2effbbyvXPPPZfu3bszYMAAevXqxbhx41ixYkXl+926daNfv34UFRVRVFREcXExS5YsoW/f7w4f4pzj17/+NT179qRXr16MGjWKBQsWAHDHHXdw6aWXVm57wQUXMGbMmMrXd911F5MnT/7Od15//fVMmbLn1B3dunVj3bp1ANx4440cfPDB9O/fn6KiIt5///3K43D11VfTs2fPytpvvPHG+A7qXqTiNjIREUlzubm5zJ8/n+3bt9O8eXNef/11OnXaPe3F3LlzueKKK3j99dfp3r07X375Jccccwz7778//fv3B+C2227jtNNOwznH7bffzqhRo5g/fz45OTkATJ8+nXbt2lV+55IlS6qt5Z577qG4uJi5c+fSokULXnvtNU466SQWLFjAkCFD+Otf/1q57Zw5cygvL6esrIysrCyKi4s55ZRTEvrZ//Wvf/Hiiy8ye/ZsmjZtyrp169i1axcAv/zlL1m9ejXz5s2jWbNmbN68md/+9rcJfX9N1AIXEZEGcfzxx/PSSy8B8Pjjj3PWWWdVvjdlyhSuvvpqunfvDkD37t256qqruO22277zPWbGZZddxn777cfLL7+ccB233nord911Fy1atADg2GOPrQzuQw45hP/85z9s376djRs30qJFC4qKipg3bx4AxcXFDBkyJKH9rVq1inbt2tG0aVMA2rVrR8eOHdm2bRsPPPAAd911F82aNQOgVatWXH/99Qn/TNVRgIuISIM488wzeeKJJ9ixYwcff/wxRx55ZOV7CxYs4LDDDttj+4EDB1ae2q7OoYceyqefflr5etSoURQVFe3xvVVt2rSJrVu30qNHj2r31aRJE4qKivjwww+ZOXMmRx55JIMGDaK4uJiVK1finKNLly4J/dzHHnssy5Yto1evXlx44YW8/fbbACxevJjCwkJatWqV0PfFS6fQRUTSzJlnwurVDfd9++0HTzyx9+369+/PkiVLePzxxxk7duwe7znn8CNr176u6vuxqp5CT0TsvoYOHUpxcTHbt29n8ODB9OzZk5tuuomCgoIaW9811WlmtGzZko8++oh33nmH6dOnc8YZZ3DLLbdw6KGH7rHtI488wh133MH69espLi5O+A+FqhTgIiJpJp6wTZaTTjqJK664grfeeov169dXrj/44IOZNWtW5fVugNmzZ9OnT58av+vf//43o0ePTmj/rVu3Jjc3ly+++IL9999/j32NGDECgCFDhvDHP/6RHTt2cNFFF1FQUMAnn3xCQUEBQ4cOBfx19AceeACAqVOn0rZtW1atWrXHvjZv3kybNm0AyMrKYuTIkYwcOZJ+/frx6KOP8sMf/pClS5eyefNmWrVqxXnnncd5551H3759KSsrS+jnqo5OoYuISIOZMGECv/rVr+jXr98e66+44gpuvvnmyo5nS5Ys4aabbuLyyy//znc457jzzjtZtWoVxx13XMI1XHnllUyePJnt27cD8MYbb/Duu+9y9tlnAz7AZ86cyddff0379u0xMwoKCnjuuecqW+AXXXQRc+bMYc6cOXTs2JHhw4fz/PPPs3nzZgCeeeYZBgwYQFZWFosWLeKzzz6r3P+cOXPo2rUrLVq04Pzzz2fSpEns2LEDgLKyssoObvWlFriIiDSYzp07c8kll3xnfVFREbfeeisnnngiJSUlZGdn85vf/IaioqLKba688kpuuOEGtm3bxqBBg5g+fXplD/SaLFq0iM6dO1e+/v3vf8/FF1/Mhg0b6NevH1lZWey3334899xzNG/eHIC8vDwKCgo4+OCDKz83ePBg3nvvPQYMGFDtfvr378+kSZM46qijMDPat2/Pgw8+CMCWLVu4+OKL+fbbb2nSpAkHHHAA999/P+BvL7v22mvp27cvrVq1onnz5owfP56OHTvGeURrlpL5wBuK5gMXEZFMEuh84CIiItLwFOAiIiIhlLEB7hyUlwddhYiISN1kbCe2NWvg9NOhZUs44ggYNAiOPBLy84OuTEREZO8yNsD32w/eeQc2bYIPPoCZM+Gee2DDBujZEwYP9qHety9kZQVdrYiIyJ4yNsArtG4NY8b4Bfyp9cWL4V//gj/+EebPh6ZN4fDDfagPHapWuoiIBE+3kcVhyxb46CMoLob33oNvv4UBA2DYML/ETLgjIiLSYGq7jSzjW+DxaNkSRozwC0BZGXz8sT8Ff+mlsGIF9Orlw3z4cDjgAKhleF8REZF6Uwu8ATgHn33mA33GDP+8WzcYPdov3boFXaGIiIRRbS1wBXiSfPklTJvml6++gn79fJgffTTUcTIdERHJMArwgJWX+85w06bBm2/Cxo3+1rUxY/wp9+ic8yIiIntQgDcyu3b5W9emTYO33vIBftxxMHYsVJmDXkREMpgCvJFbtw5eeQWmToUlS/z952PH+k5xTZsGXZ2IiARFAR4ipaW+dT51qu8Q1749nHoqfP/7EJ03XkREMoRuIwuRJk1gyBC/gO8A9+yz8MMfQnY2nHyyX/bdN9g6RUQkWGqBh8jXX8Pzz8Nzz8GOHfBf/wWnnQZt2wZdmYiIJIPmA08TBQVw/vk+xJ980o/RfvbZvkX+t7/5EeNERCQzKMBDKi8PfvITePVV+MMfYPVqf538xz+Gt9/2g8uIiEj6UoCngU6d4Gc/88F91VXwz3/6SVemTPGn3UVEJP0owNPMwQfD7bf7AWM6dIAf/QjOOMPfc65WuYhI+lCAp6lmzXx4v/Ya/N//wYsv+vvKH33UDyQjIiLhpgDPAAceCL//vQ/xtWt9kN98M2zeHHRlIiJSVwrwDNKmDVx5Jbz7LnTsCMccA7/5DWzdGnRlIiKSKAV4BsrOhvHj/fSn+fl+hrQ77tCpdRGRMFGAZ7DsbH8r2owZEInAiBF+CFcREWn8FOBC06Zw8cX+Gvnzz/te66tXB12ViIjURgEuldq2hfvug0mT/AQq//hH0BWJiEhNFODyHcOGwRtvwMsvw//7f7BzZ9AViYhIVQpwqVZuLjz4oJ8V7cQT4Ztvgq5IRERiaTpRqdU550BhIZxwAjz9tL/9TEREgqcWuOzViBFw992+c9vGjUFXIyIioACXOB16KFx3nQ9x3S8uIhI8BbjEbcwYOP54uOeeoCsREZGkB7iZdTGz6Wa20MwWmNkl0fX5Zva6mX0WfcxLdi1SfxdeCE89pWlKRUSClooWeClwuXPuIGAQcJGZ9QF+AUxzzvUEpkVfSyOXnQ2/+pVfREQkOEkPcOfcKufc7OjzzcBCoBNwMvBodLNHgVOSXYs0jOOOgyVLYMWKoCsREclcKb2NzMy6AYcA7wP7OudWgQ95M2tfw2cmAhMBCgsLG7Se2bNns3Tp0gb9zkyRm3sQL7/s+MlPegddiohIRkpZgJtZS+AfwKXOuU1mFtfnnHP3A/cDDBw40DVkTdu3b6dZs2Y0a9asIb82IxQUlLNiRXnQZYiIZKyU9EI3s2x8eP/VOfdMdPUaM+sQfb8DsDYVtUjDyM8vYd26rKDLEBHJWKnohW7AQ8BC59zvYt56HhgffT4eeC7ZtUjDyc/fxfr12UGXISKSsVJxCn0ocA4wz8zmRNddDdwC/N3MzgeWAqenoBZpIG3blvDNNwpwEZGgJD3AnXPvAjVd8B6d7P1LcuTnl7BhQ07QZYiIZCyNxCZ10rJlGVu2NMG5Bu1XKCIicVKAS51EIv6kSnm5eqKLiARBAS71UlamABcRCYICXOosN7eEDRsU4CIiQVCAS521abOTlSsV4CIiQVCAS521abOTNWvUiU1EJAgKcKmzvLydrFoVdBUiIplJAS511qbNDtasCboKEZHMpACXOsvL26UAFxEJiAJc6sy3wOObVU5ERBqWAlzqrE2bnaxdqwAXEQmCAlzqrGnTErZsUYCLiARBAS51Zgag28hERIKgAJd603DoIiKppwCXemndupQNG4KuQkQk8yjApc7MjH33LeHGG+HVV1GQi4ikUJOgC5DwikQiXHzxSjZuzGPmTLjrLvj2W+jVCwYPhkGDoE8fyMoKulIRkfSjAJc6MzOaNi3lmGPgmGP8OufgP/+BmTPhD3+ABQugeXM4/HAf6IMGQbt2wdYtIpIOFOBSZ2ZGaWlplXVw4IF+GT/er9uyBT780If6Aw/AunVw0EEwbBgMHw7dulX0aBcRkXgpwKXOqgvw6rRsCaNG+QV8r/VPP4UZM+CXv4Qvv/QhPmyYX/r0gYh6Z4iI1EoBLnVmZpSVlSX8uUjEh3SfPvDf/+3XffWVD/Q774RPPvGn2YcPh9GjoV8/BbqISFUKcKmzSCQSVws8Hl27wjnn+AVg7Vp4+224916YNw86d/ZhPmYM7L9/g+xSRCTUFOBSZ3VtgcejfXs4/XS/ACxdCtOmwbXX+lPuffr4QD/6aNh336SUICLSqCnApc6SGeBVFRbCeef5xTl/mn3aNLjgAli/HoYOhbFj/e1r2dkpKUlEJFAKcKmzeDuxNfx+4eCD/TJ5MpSUwHvvwQsvwFVX+dPtY8fC8cf7lryISDpSgEudpbIFXpvsbBg50i8AS5bA1Klw/vmwebNff8opMGCAblcTkfShAJc6i0QijSLAq+rWDS680C/bt8P06X5QmY8/9qfaTz3Vn2rXCHEiEmYKcKmX8vJynHNYI23aNm/uT6ePHQtlZf5U+z/+AZdfDkVFcMYZMGKEwlxEwkcBLnVmZpgZ5eXlZIUgAbOy/L3lw4f7jnCzZ8MTT/jr5kOHwtlnw2GH6TS7iISDAlzqxTlHWVlZKAI8lpkP68MO8yPDvfsuPPggXHwxfO97MG6c7jcXkcZNAS71UtECD7NIZHfLfNcuPzXqlVfCtm3+trVTToGcnKCrFBHZkwaolHprjB3Z6ionB0480V8nf/hh+PxzH+xXXgmLFgVdnYjIbgpwqbewt8Br0qGDvz5eXOxPq//qV/7e8pdf9tfQRUSCpACXekmHU+h7E4n4MdiffBLuvx/efBOOOgoeegh27Ai6OhHJVApwqZeKTmyZoksXuO023wrftMnfgnbfff7auYhIKinApd7SvQVendat4bLL4K23/GAxw4b5FnlJSdCViUimUIBLvWVSC7yq5s19kL/5Jqxb54dtffvtoKsSkUygAJd6y8QWeFW5ufDzn/ve6w8+6O8jX7Mm6KpEJJ0pwKXeFOC77bcf/OUvMGGCH3P9gQfUY11EkkMBLvWSaZ3Y4jVypL8+/sUXfrz1DRuCrkhE0o0CXOpFAV6znBy4+WaYONFPpqKBYESkISnApV4ikQilpaVBl9GojRkDjz3mh2V9552gqxGRdKEAl3oxM7XA49CjB7zwAlxzDcyaFXQ1IpIOFOBSL2amFnic2raFv/8dJk2Czz4LuhoRCTsFuNSLmVGi0UviFttLXX/3iEh9JD3AzexhM1trZvNj1l1vZivMbE50GZvsOiQ5dAo9cT17wkkn+VvMRETqKhUt8D8Bx1Wz/vfOuaLoMjUFdUgSqBNb3VxyCfz5z/DNN0FXIiJhlfQAd87NAPRrKk2pBV43OTm+Q9uvfx10JSISVk0C3PckMxsHzAIud85VO9SFmU0EJgIUFhamsDyJh5mxceNG5s6dG3QpodO5M8yY0YNt2yK0aNEi6HJEJGSCCvB7gRsAF338LTChug2dc/cD9wMMHDhQg1I2Mjk5OWzbto1ly5YFXUoobdvWkU2byhTgIpKwQALcOVc5zYOZPQC8GEQdUn9mRm5ubtBlhFYkYpSVaSx5EUlcILeRmVmHmJenAvNr2lYknTVvXsrGjQpwEUlc0lvgZvY4MBJoZ2bLgeuAkWZWhD+FvgS4INl1iDRGbdrsYM2abPr0CboSEQmbpAe4c+6salY/lOz9ioRBmzY7WbMmO+gyRCSENBKbSIDy8nayZs3etxMRqUoBLhKgvLydrF2r/w1FJHH6zSESoLy8Er7+Wv8bikji9JtDJED5+bsU4CJSJ/rNIRKgvLxdrF8f5ICIIhJWCnCRAOXmlrF5s/43FJHE6TeHSIAiEcNpgGARqQMFuEiAIpEIzjmFuIgkTAEuEiA/lnwJGzcGXYmIhI0CXCRAZkabNjtYvTroSkQkbBTgIgEyM/bZZ4dGYxORhCnARQLkA3y7AlxEEqYAFwlQRQt81Sr1YhORxCjARQK0+xq4AlxEEqMAFwlYXt4uBbiIJEwBLhKwNm12sHix6VYyEUmIBmEWCVjz5mUcd1wpZ52Vw6ZN0Ls3DBrklz59IKI/s0WkGgpwkUZg8uSdXHNNDuXlsGgRzJwJd94JCxdCixZwxBEweDAceSS0bRt0tSLSGCjARRqB8vJywLe2DzrIL+ed59/btAk+/NCH+n33wTffQI8ecNRRMHw49OoFZgEWLyKBUICLBMzMKCsrq/H91q1h9Gi/ADgHixfDu+/Crbf6Fvt++/kwHzYMBgyArKwUFS8igVGAiwTMOVfZAo+HGfTs6ZeKVvrKlT7QH34YPv4YWrWCoUN9oB9xBDRtmqTiRSQwCnCRRiCRAK9Ox47wwx/6BWDDBnjvPXjhBbj2Wn8dfdQo34ovKlLHOJF0oAAXaQRqO4VeF3l5cMIJfgF/3fytt+Chh2DuXH/KffRoGDMGDjhA19BFwmivAW5mbYBxQLfY7Z1zk5NXlkhmqW8LfG/y8+EHP/ALwPLl8OabcMMN/nr6gQfCscfC977ntxWRxi+eFvhUYCYwD0jubxmRDOSca/AW+N507gzjxvnFOd8R7tVX4cc/hp07fet87FjfIU6tc5HGKZ4Ab+ac+1nSKxHJUEEEeCwzP3hM795wySWwZYtvnd93H8yZA337+jA/5hjfOU5EGod4AvwvZvZT4EVgZ8VK59w3SatKJIPs7TayVGvZEk46yS/Owfz5MHUq3H23D/BTT/Xv6VS7SLDiCfBdwG3ANUDFjAsO2D9ZRYlkkkgkQmlpadBlVMsM+vXzy89/7m9X++c/4ayzfE/2k0+GU07xneJEJLXiCfCfAQc459YluxiRTGRmjTbAq+rYES66yC9ffw3PPw8//Sls3+5b5WeeCe3bB12lSGaI527QBcC2ZBcikqnCFOCxCgrg/PP9veZPPw377APjx/sgf/JJH+oikjzxtMDLgDlmNp09r4HrNjKRBhDWAI/Vpo0P7/Hj/S1qf/ub7/R24IG+p/vw4erNLtLQ4mmBPwvcCBQDH8UsItIAGlsntvrq3Bn+53/gnXdg0iR45hk/rOvvfgfr1wddnUj62GsL3Dn3aCoKEclUjbkTW32YwSGH+GX7dnjqKT/Ua6dOcOmlcOihQVcoEm4aEVkkYOlwCn1vmjf3p9KnTYOf/QzuuMOP+vbii5DkQehE0pYCXCRg6XYKfW+KiuDRR+GRR6C42M9r/thjkEGHQKRBKMBFApYJLfDqdOwIN90Er7wCX37ppz79+9/VIheJ114D3Mx6mdkDZvaamb1ZsaSiOJFMYGZJn8ykMWvd2k95+tJLMG8ejBgB06cHXZVI4xfPbWRPAfcBD+BvKRORBpSundgSlZfnZ0dbs8b3Yn/wQZgyBTp0CLoykcYpngAvdc7dm/RKRDJUpl0D35t99/XXyGfMgNNOgzPO8CO/ZWUFXZlI4xLPNfAXzOxCM+tgZvkVS9IrE8kgzrmMPo1eneHD4a23/C1oY8fCihVBVyTSuMTTAh8ffbwyZp0mMxFpIBYdoqy8vJxIRP1KY2Vn+0lUxoyB00/3p9hHjw66KpHGIZ6BXLqnohCRTJbpHdn25rDD/JSmZ58NX30FEyYEXZFI8GoMcDM72jn3ppn9oLr3nXPPJK8skcyj6+C1a9MGnnvOz362dStcfHHQFYkEq7YW+AjgTeDEat5zgAJcpAGpBb532dnw0EO+JV5Q4KcvFclUNQa4c+666ON5qStHJDOpJ3r8srLgz3+Gk0+Gbt1g0KCgKxIJRtJ7zJjZw2a21szmx6zLN7PXzeyz6GNesusQaczUCz0xTZv6lvjPf64hWCVzpaLL65+A46qs+wUwzTnXE5gWfS2S0RTgienUCY47Dh5+OOhKRIKR9AB3zs0Avqmy+mSgYprSR4FTkl2HSGOnU+iJu+wyH+AbNgRdiUjqxXMf+B7MbCCwyjlXn2EV9nXOrQJwzq0ys/b1+C6RtLBp0yayNNxYwk44oSkvvBBh3LjmQZciklIJBzhwMdDfzP7jnDujoQuqyswmAhMBCgsLk707kUBkZWWxYMGCoMsIpU2b9uNf/2rFuHEHBV2KSEolHODOufEAZtaqHvtdY2Ydoq3vDsDaWvZ3P3A/wMCBA1099inSaLVs2TLoEkKrc+cIb79dl7aISLjVNpDLobV90Dk3ux77fR4/ROst0cfn6vFdIpLB2rUrYf36pjjnKoelFckEtf3Z+tvoYzNgIDAXMKA/8D5wVDw7MLPHgZFAOzNbDlyHD+6/m9n5wFLg9LoULyLStu0uNmxoRnl5ufoQSEapbSCXUQBm9gQw0Tk3L/q6L3BFvDtwzp1Vw1uakkBE6i0nx1FaGlGAS8aJ5zay3hXhDeCcmw8UJa8kEZHE6TY8yTTx9PxYaGYPAo/hx0D/MbAwqVWJiCQgO7ucrVvLadYs6EpEUieeFvh5wALgEuBS4JPoOhGRRiEvbwcrV+omFcks8cwHvgP4fXQREWl08vN9gPfrF3QlIqmz1wA3s57AzUAffI90AJxz+yexLhGRuPkAD7oKkdSK5xT6I8C9QCkwCvgz8JdkFiUikgh/Cj3oKkRSK54Ab+6cmwaYc+4r59z1wNHJLUtEJH75+TtYvVqDuEhmiacX+g4ziwCfmdkkYAWgyUdEpNHIz9/BrFmpmB1ZpPGI51/8pUALYDJwGP42svHJLEpEJBFt2uxgzRq1wCWzxNML/UMAM3POOd0+JiKNTnZ2Obt2BV2FSGrttQVuZoPN7BOig7eY2QAz+0PSKxMRiZNzDud0H7hklnhOod8OfA9YD+CcmwsMT2ZRIiKJMDNychzbtwddiUjqxNXrwzm3rMoqDTosIo1GJBKhoGAXq1YFXYlI6sQT4MvMbAjgzCzHzK5AY6GLSCNiZrRrt0v3gktGiec2sv8G7gA6AcuB14CLklmUiEgizIyDDtrKZZfBoYfC4YfDEUdAnz7QJJ7fciIhFE8v9HXAj1JQi4hInZgZw4atZ/Lk7ixcCB98APfeC598AllZMGCAD/TDD4cePcB0x5mkgRoD3Mzuwk8fWi3n3OSkVCQikiAzo6ysjKws6NvXLxMm+Pd27IC5c+HDD+GGG+Dzz6FlSxg40Af64YdDx47B1i9SF7W1wGfFPP9f4Lok1yIiUieRSITS0tJq32vWDI480i8VNm2Cjz7yof7YY7ByJXTuDEOG+KWoCLKzU1S8SB3VGODOuUcrnpvZpbGvRUQak4oWeLxat4ZRo/xSYdkyKC72gX7FFT7AjzgChg6FwYMhPz8JhYvUQ7zdOzRCgog0WokGeHW6dIEzzvALwNatvoVeXAz33QcbNkDv3r6FPnQo9Oqla+kSLPXPFJHQM7MaT6HXVW4ujBzpF4Dycli0yAf6Lbf454WFMGKE36Z3bwW6pFZtndg2s7vl3cLMNlW8BTjnXOtkFyciEo9kBHhVkQgcdJBfzj/fr/vqK3j7bbjtNvj0U9+Krwh9BbokW23XwFulshARkbqqrRNbMnXtCuPG+QVg6VIf6FOmwMKFvmPcmDFw7LHQrVvKy5M0p1PoIhJ6ZkZ5eXnQZVBYCOec4xfwgf7GG/CLX/jnhxziw3zUKN+RTqQ+FOAikhbKy8txzmGN6Lx1YaG/H33CBH8N/d//htdeg3vugbIyf/187Fg/elwkrpkpRHZTgItI6FWEdllZGU0a6dipkQgcdphfrrrK93J/6y145BGYNAn69YMTT4TRo30HOpG9aZz/0kVEEtRYTqPHKzcXvv99vzgH8+bBCy/A7bfv+V6XLkFXKo2VAlxE0kaYAjyWGfTv75drroE1a+Dll+HSS2HdOn+a/bTT/DjuIhUU4CKSNuo7mEtjse++cO65ftm61Yf5Ndf4IV+PP96Hec+eQVcpQVOAi0jaCGsLvDa5uT6wTzsNtm2DV1+F666DFSvgBz+As86C9u2DrlKCoH6PIpI20qUFXpMWLeDUU+Fvf4OXXoK8PH8P+sknw9NP+5nXJHMowEUkbaRjC7wmLVv68H7lFbj7bli82Pdgv/BC+PjjoKuTVFCAi0jayKQAj9Wlix8s5t13/SAyU6bA0UfDX/6iVnk6U4CLSLBo9ZgAABEuSURBVNpI91Poe2Pmpz7985/hqadg7Vo/Lvv//i98/XXQ1UlDU4CLSNrI1BZ4ddq2hcsvh/fegz594PTT/en1xYuDrkwaigJcRNKCcy7jW+DVycry4T19uu+xftllfja1r74KujKpLwW4iKQNtcBrZgbDhvnR3s49F37yE5g8GVavDroyqSsFuIikDbXA4zNsmJ9U5fvf963zm2+GnTuDrkoSpQAXkbRgZoHMCR5WZvC97/kJVfLzfWe3l14KuipJhEZiE5G0oACvm6wsuOAC3xK/9lq46y4/Y1rnzv72tC5d/PP8fB/60ngowEUkLSjA6yc/389TvmoVfP45LFsGn3ziT7UvWwbffOO3y8mBTp32DPeK5/vso5BPJQW4iKQFBXjD6NDBLzXZudOPw75sGSxf7kd9e+kl/3rjRr9Ns2Y+2GPDvSLsW7dOzc+RCRTgIpIWIpGIAjwFmjaF/ff3S022b98d8suWwaxZ8Oyz/vnmzX7+8xYtqm/Fd+niJ3CRvVOAi0haMDP1Qm8kmjeHAw7wS022bvUt+IqW/MyZfvS4Zcv8e2Z+vPeuXaGw0D9WPO/QwV+7z3QKcBFJCzqFHi65uXDggX6pyaZNsHSpH3Rm6VKYO9c/rloFZWXQpIlvvceGe9euvhXfokXqfpagKMBFJC0owNNP69bQt69fqlNS4lvvFSH/wQe+Fb90qT+ND1BQ8N0WfNeufqjZsHe4CzTAzWwJsBkoA0qdcwODrEdEwkun0DNPdjZ07+6X6jgH69btbsEvWQIzZvjX69b5bVq02B3qsY+dO/sWfmPWGMob5ZxbF3QRIhJu6sQmVZn5FnhBAQysoXm4dau/7l4R8q+95p+vWAGlpT7ECwuhWze/dO/uHzt1Cj7gG0OAi4jUm5lpLHRJWG4u9O7tl+rs2uVP0y9ZAl9+CW+84Z8vX+6vw2dn7xnww4fXfEagoQUd4A54zcwc8Efn3P1VNzCzicBEgMLCwhSXJyJhYWZs2bKFqVOnBl1KaBUVFdGxY8egy2hUcnJqv21u1y7fgl+yxC/ffJM5AT7UObfSzNoDr5vZp865GbEbREP9foCBAwe6IIoUkcYvEonQtm1bnNOvibrYtGkT2yt6fknccnKgRw+/pFqgAe6cWxl9XGtm/wSOAGbU/ikRkepFIpqfqa7UhyB8AvvXbma5Ztaq4jlwLDA/qHpERDJZJBKhpKQk6DIkAUG2wPcF/mn+RrwmwN+cc68EWI+ISMbSbXjhE1iAO+e+AAYEtX8REdlNA+GEjy4YiYiIAjyEFOAiIqJT6CGkABcREfVCDyEFuIiIqAUeQgpwERFRgIeQAlxERNSJLYQU4CIiohZ4CCnARUSESCSiAA8ZBbiIiABQXl6uyWBCRAEuIiKYmeZUDxkFuIiIAOCc02n0EFGAi4gIgFrgIaMAFxGRSmqBh4cCXEREKqkFHh4KcBERAXQKPWwU4CIiAqgTW9gowEVEpJJa4OGhABcRkUpqgYeHAlxERCqpBR4eCnAREamkAA8PBbiIiADqxBY2CnAREQEU4GGjABcREcBPKVpaWhp0GRInBbiIiAB+IBe1wMNDAS4iIoAPcLXAw0MBLiIigA/wkpKSoMuQOCnARUQE0Cn0sFGAi4gIoE5sYaMAFxERQC3wsFGAi4gIoE5sYaMAFxERQC3wsFGAi4gIoBZ42CjARUQE8J3Y1AIPDwW4iIgAOoUeNgpwEREBdge4cy7oUiQOCnAREQF8gDvnFOAhoQAXEZE9lJeXB12CxEEBLiIilXQdPDwU4CIisge1wMNBAS4iIpXMTAEeEgpwERGp5JzTKfSQUICLiMge1AIPBwW4iIjsQS3wcFCAi4jIHtQCDwcFuIiI7EEBHg6BBriZHWdmi8xssZn9IshaREREndjCJLAAN7Ms4B7geKAPcJaZ9QmqHhER8dQCD4cgW+BHAIudc18453YBTwAnB1iPiEjGUws8PJoEuO9OwLKY18uBI6tuZGYTgYkAhYWFqalMRCRDmRlr165lx44dQZcSSllZWfTo0QMzS/q+ggzw6n6670yB45y7H7gfYODAgZoiR0QkiXJzc1m5ciUrV64MupRQKi8vp1u3bjRpkvx4DTLAlwNdYl53BvQvRkQkQNnZ2eyzzz5BlxFa3377bcr2FeQ18A+BnmbW3cxygDOB5wOsR0REJDQCa4E750rNbBLwKpAFPOycWxBUPSIiImES5Cl0nHNTgalB1iAiIhJGGolNREQkhBTgIiIiIaQAFxERCSEFuIiISAgpwEVEREJIAS4iIhJCCnAREZEQUoCLiIiEkAJcREQkhBTgIiIiIaQAFxERCSEFuIiISAgpwEVEREJIAS4iIhJCCnAREZEQUoCLiIiEkAJcREQkhBTgIiIiIaQAFxERCSEFuIiISAgpwEVEREJIAS4iIhJCCnAREZEQUoCLiIiEkAJcREQkhBTgIiIiIaQAFxERCSEFuIiISAgpwEVEREJIAS4iIhJCCnAREZEQUoCLiIiEkAJcREQkhBTgIiIiIaQAFxERCSEFuIiISAgpwEVEREJIAS4iIhJCCnAREZEQUoCLiIiEkAJcREQkhBTgIiIiIaQAFxERCSEFuIiISAgpwEVEREIokAA3s+vNbIWZzYkuY4OoQ0REJKyaBLjv3zvnpgS4fxERkdDSKXQREZEQCrIFPsnMxgGzgMudcxuCKKKkpAQzC2LXIiIidWbOueR8sdkbwH7VvHUNMBNYBzjgBqCDc25CDd8zEZgYfXkgsKiBSmzXvHnzkuzs7BYN9H0ZadeuXS1ycnK2BV1HmOkYNgwdx4ah41g/5eXlZVu2bCnHZ1xD6OqcK6jujaQFeLzMrBvwonOub4r3O8s5NzCV+0xHOo71p2PYMHQcG4aOY/2l6hgG1Qu9Q8zLU4H5QdQhIiISVkFdA/+NmRXhT6EvAS4IqA4REZFQCiTAnXPnBLHfKu4PuoA0oeNYfzqGDUPHsWHoONZfSo5h4NfARUREJHG6D1xERCSEMjLAzew4M1tkZovN7BdB19OYmFkXM5tuZgvNbIGZXRJdn29mr5vZZ9HHvJjPXBU9lovM7Hsx6w8zs3nR9+60DLvh3syyzOzfZvZi9LWOYYLMrI2ZPW1mn0b/TQ7WcUyMmV0W/X95vpk9bmbNdAz3zsweNrO1ZjY/Zl2DHTcza2pmT0bXvx+9IysxzrmMWoAs4HNgfyAHmAv0CbquxrIAHYBDo89bAf8B+gC/AX4RXf8L4Nbo8z7RY9gU6B49tlnR9z4ABgMGvAwcH/TPl+Jj+TPgb/jbJNExrNMxfBT4SfR5DtBGxzGh49cJ+BJoHn39d+BcHcO4jt1w4FBgfsy6BjtuwIXAfdHnZwJPJlpjJrbAjwAWO+e+cM7tAp4ATg64pkbDObfKOTc7+nwzsBD/S+Bk/C9Too+nRJ+fDDzhnNvpnPsSWAwcEb1VsLVz7l/O/wv9c8xn0p6ZdQa+DzwYs1rHMAFm1hr/S/QhAOfcLufct+g4JqoJ0NzMmgAtgJXoGO6Vc24G8E2V1Q153GK/62lgdKJnNTIxwDsBy2JeL4+ukyqip3QOAd4H9nXOrQIf8kD76GY1Hc9O0edV12eK24H/Acpj1ukYJmZ/4GvgkeiliAfNLBcdx7g551YAU4ClwCpgo3PuNXQM66ohj1vlZ5xzpcBGoG0ixWRigFf3F4664ldhZi2BfwCXOuc21bZpNetcLevTnpmdAKx1zn0U70eqWZfRxzCqCf4U5r3OuUOArfjTljXRcawieo32ZPxp3Y5Arpn9uLaPVLMuo49hnOpy3Op9TDMxwJcDXWJed8afUpIoM8vGh/dfnXPPRFevqRhBL/q4Nrq+puO5PPq86vpMMBQ4ycyW4C/RHG1mj6FjmKjlwHLn3PvR10/jA13HMX5jgC+dc18750qAZ4Ah6BjWVUMet8rPRC9v7MN3T9nXKhMD/EOgp5l1N7McfOeB5wOuqdGIXoN5CFjonPtdzFvPA+Ojz8cDz8WsPzPao7I70BP4IHp6abOZDYp+57iYz6Q159xVzrnOzrlu+H9fbzrnfoyOYUKcc6uBZWZ2YHTVaOATdBwTsRQYZGYtoj/7aHy/Fh3DumnI4xb7Xafhf08kdlYj6J5+QSzAWHzv6s+Ba4KupzEtwFH40zgfA3Oiy1j8tZlpwGfRx/yYz1wTPZaLiOmZCgzEj3P/OXA30YGDMmkBRrK7F7qOYeLHrwg/5fDHwLNAno5jwsfwf4FPoz//X/A9pXUM937cHsf3GyjBt5bPb8jjBjQDnsJ3ePsA2D/RGjUSm4iISAhl4il0ERGR0FOAi4iIhJACXEREJIQU4CIiIiGkABcREQkhBbhIgsxsS5XX55rZ3Ql+x0nWgDPhmZ+168I4t90SxzbXm9kV0ef/Z2Zjatn2FDPrE3+1yWFmA83szqDrEEkVBbhIiplZE+fc8865Wxrwa9vgZzdqcM65Xznn3qhlk1PwszEFyjk3yzk3Oeg6RFJFAS7SgMysq5lNM7OPo4+F0fV/MrPfmdl04NbYVruZzYlZtpvZiOi8w89Gv2emmfWPbnt9dJ7it8zsCzOrCKxbgB7R77jNzFpG9z87OhfxXmfcM7NrzM9l/AZwYMz6P5nZadHnt5jZJ9G6ppjZEOAk4LbovnuY2U/N7EMzm2tm/zCzFjHfc6eZFUdrPy1mH/8TrXOumd0SXdfDzF4xs4/M7B0z672X+kfa7rnXazpOVT+zxcxuje7jDTM7IuYzJ+3tmIkEKujRbrRoCdsClLF7lLo5+OEq746+9wIwPvp8AvBs9PmfgBfZPUfwuRWfifneE4F3gGzgLuC66PqjgTnR59cDxfjRtNoB66Pbd2PPeYub4KcxJLrdYnaPALWlmp/pMGAefrrJ1tHtr4ip/TQgHz/KVMX3tIl9P+a72sY8/zVwccx2T+EbDn3w0/oCHB/9mVpEX+dHH6cBPaPPj8QPNVnbf5eR7B71rtrjVM1nHLvnZ/4n8Fr0eA6oOOZatDTWpQkikqjtzrmiihdmdi5+uESAwcAPos//Avwm5nNPOefKqvtCM+sJ3AYc7ZwrMbOjgP8CcM69aWZtzWyf6OYvOed2AjvNbC2wb3VfCdxkZsPxU5p2im63uoafaRjwT+fctmg91c0PsAnYATxoZi/h/yCpTl8z+zX+tH5L4NWY9551zpUDn5hZRd1jgEcq9u2c+8b8bHhDgKds9xTJTWvYX02qO07Lq2yzC3gl+nwesDN6/Ofh/ygSabQU4CLJFTtW8dbqNjA/x/XfgZ865ypmKqptqsGdMevKqP7/4x8BBcBh0UBagh97Od5av/umc6VmdgR+QowzgUn4swNV/Qk4xTk3N/rHzciY92Jrt5jHqvuOAN/G/qFUB/EcpxLnXMW+yys+45wrNz9DlEijpWvgIg2rGB9u4EP03Tg+8wi+BfpOzLoZ0c9jZiOBda72edk3A61iXu+Dn5O8xMxGAV33UsMM4FQza25mrfCn8/cQbRXv45ybClyKn2ikun23AlaZn5b2R3vZL/jT1hNirpXnR3/WL83s9Og6M7MB0eenmtnNcXyvSFrTX5giDWsy8LCZXQl8DZxX28Zm1hV/fbmXmU2Irv4J/hruI2b2MbCN3dMOVss5t97M3jOz+cDLwK3AC2Y2C3+d/tO9fH62mT0Z3fYr/LX4qloBz5lZM3yr+bLo+ieAB6IdxU4DrgXej37PPPYM9+r2/YqZFQGzzGwXMBW4Gh/+95rZL/HXpZ8A5gI98KfzRTKaZiMTkVAxs8eAy5xzXwddi0iQFOAiIiIhpGvgIiIiIaQAFxERCSEFuIiISAgpwEVEREJIAS4iIhJCCnAREZEQUoCLiIiE0P8HqPWtpWkjAlIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "fig.subplots_adjust(left=None, bottom=None, right=None, top=None,\n",
    "                    wspace=0.25, hspace=0.25)\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.plot(x, mfusghead[0, 0, :], linewidth=0.75, color='blue', label='MODFLOW-USG')\n",
    "ax.fill_between(x, y1=botm[1, 0, :], y2=-5, color='0.5', alpha=0.5)\n",
    "leg = ax.legend(loc='upper right')\n",
    "leg.draw_frame(False)\n",
    "ax.set_xlabel('Horizontal distance, in m')\n",
    "ax.set_ylabel('Head, in m')\n",
    "ax.set_ylim(-5,25);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}