p9_make_pred_model.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Chapter 10 problem using Flopy###\n",
      "\n",
      "requires flopy, available at <link>https://pypi.python.org/pypi/flopy/2.2.306</link>\n",
      "or via pip install: <link>https://pypi.python.org/pypi/pip</link>\n",
      "\n",
      "much of this code follows the flopy tutorial: <link>https://flopy.googlecode.com/svn/sphinx/_build/html/tutorial2.html</link>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Define general model characteristics####"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import pandas as pd\n",
      "import flopy\n",
      "\n",
      "%matplotlib inline\n",
      "\n",
      "#model domain and grid definition\n",
      "Lx = 1500.\n",
      "Ly = 1500.\n",
      "ztop = 600.\n",
      "zbot = 450.\n",
      "nlay = 1\n",
      "nrow = 15\n",
      "ncol = 15\n",
      "delr = Lx / ncol\n",
      "delc = Ly / nrow\n",
      "delv = (ztop - zbot) / nlay\n",
      "botm = np.linspace(ztop, zbot, nlay + 1)\n",
      "\n",
      "# properties\n",
      "Khvalues = {1: 44, 2: 5.5} # dictionary of K values by zone number (1 = sand, 2 = silt)\n",
      "Vani = 1.\n",
      "sy = 0.1\n",
      "ss = 1.e-4\n",
      "laytyp = 1\n",
      "\n",
      "# global BC settings\n",
      "m_riv = 2 # riverbed thickness\n",
      "w_riv = 100 # riverbed width\n",
      "R = 0.0001 # recharge rate\n",
      "Qleak = 45000 # flow through southern boundary (pos. = inflow)\n",
      "Rcond = 150000\n",
      "\n",
      "# pumping well for transient simulation\n",
      "QA2 = -20000\n",
      "pumping_well_info2 = [1, 7, 11, QA2] # l, r, c, Q\n",
      "QA3 = -30000\n",
      "pumping_well_info3 = [1, 9, 9, QA3] # l, r, c, Q\n",
      "\n",
      "# Stress Periods\n",
      "nper = 3\n",
      "perlen = [1, 3, 365]\n",
      "nstp = [1, 10, 20]\n",
      "tsmult = [1, 2, 2]\n",
      "steady = [True, False, False]\n",
      "\n",
      "# create flopy objects for packages\n",
      "modelname = 'P10'\n",
      "mf = flopy.modflow.Modflow(modelname, exe_name='mf2005')\n",
      "dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc,\n",
      "                               top=ztop, botm=botm[1:],\n",
      "                               nper=nper, perlen=perlen, tsmult=tsmult, nstp=nstp, steady=steady)\n",
      "\n",
      "# Variables for the BAS package\n",
      "# Note that changes from the previous tutorial!\n",
      "ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)\n",
      "strt = 515. * np.ones((nlay, nrow, ncol), dtype=np.float32) # starting heads\n",
      "bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)\n",
      "\n",
      "# Output Control\n",
      "words = ['head','drawdown','budget', 'phead', 'pbudget']\n",
      "save_head_every = 1\n",
      "oc = flopy.modflow.ModflowOc(mf, words=words, save_head_every=save_head_every, compact=True)\n",
      "\n",
      "# solver\n",
      "'''\n",
      "note: had trouble getting pcg to solve for some k values for the silt. Switched to NWT and it ran faster with no problems.\n",
      "'''\n",
      "#pcg = flopy.modflow.ModflowPcg(mf, relax=1)\n",
      "nwt = flopy.modflow.ModflowNwt(mf)\n",
      "\n",
      "# recharge package\n",
      "rch = flopy.modflow.mfrch.ModflowRch(mf, nrchop=3, irchcb=0, rech=R, irch=1, extension='rch', unitnumber=19)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Make the hydraulic conductivity array (upw package):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# cell numbers with silt\n",
      "silt = np.append(np.arange(78, 84), np.arange(93, 99))\n",
      "silt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 32,
       "text": [
        "array([78, 79, 80, 81, 82, 83, 93, 94, 95, 96, 97, 98])"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# now creaty Kzones array and assign second zone value to cells with silt\n",
      "Kzones = np.ones(225) * 1\n",
      "Kzones[silt] = 2\n",
      "Kzones = np.reshape(Kzones, (15, 15))\n",
      "\n",
      "# show the Kzones array\n",
      "plt.imshow(Kzones, interpolation='none')\n",
      "\n",
      "# save the Kzones for calibration with PEST\n",
      "np.savetxt('../pest/Kzones.dat', Kzones, delimiter=' ', fmt='%i')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD7CAYAAABOrvnfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD4BJREFUeJzt3X9Ilfffx/GXTSNIqhWpy2M7Uok/U8sQYqNanYIgaSaj\nXCj9Giw2thZRf92jGyqdi35s/40sY7GC/TFdOFmu2WIirZ02xhqzb3jYUctBTTd3KtN97j9aYt9b\nTa/Odaw+zwcc8Mfxen8Cn12e4/H6RBljjABYYdxYLwBA5BA8YBGCByxC8IBFCB6wCMEDFol268AT\nJ85RKPQftw4PYAiLFi1SQ0PDoJ+Lcuv38FFRUZLeHeYeDZIWuzH6MZvJ3Kd35uM6d7eGypof6QGL\nOA6+rq5OqampmjNnjsrLy8O5JgAucRR8X1+f3njjDdXV1eny5cv65JNP9Msvv4zyKF4nox/RWMxk\n7tM788mb6yj4CxcuaPbs2fJ6vYqJidHatWtVXV09yqN4nYx+RGMxk7lP78wnb66j4Nva2pSUlNT/\nvsfjUVtbm6MFAIgcR7+Wu/cM/Eg0DHjbq7H73xB4mgX+vT2co+ATExMVDAb73w8Gg/J4PIPcc7GT\nwwMYFa8ePJmeG/Kejn6kz8vL05UrVxQIBNTT06NTp06poKDAyaEARJCjM3x0dLQ+/PBDrVixQn19\nfdq0aZPS0tLCvTYAYTaGr7QD4A5eaQdABA9YheABixA8YBGCByxC8IBFCB6wCMEDFiF4wCIED1iE\n4AGLEDxgEYIHLELwgEUIHrAIwQMWIXjAIgQPWITgAYsQPGARx8EHg0EtWbJEGRkZyszM1OHDh8O5\nLgAucHSZakmKiYnRgQMHlJOTo+7ubs2fP18+n4/LVQOPMcdn+ISEBOXk5EiSYmNjlZaWpvb29rAt\nDED4heUxfCAQ0KVLl5Sfnx+OwwFwySMH393draKiIh06dEixsbHhWBMAlzh+DC9Jd+/e1Zo1a7R+\n/XqtXr16kHs0DHjbK3aPBdwQ0Eh3j3W81ZQxRqWlpZo2bZoOHDjw/w/MVlPAGHFhq6lvv/1WH3/8\nsb7++mvl5uYqNzdXdXV1jpcIwH2Of6R/4YUX9M8//4RzLQBcxivtAIsQPGARggcsQvCARQgesAjB\nAxYheMAiBA9YhOABixA8YBGCByxC8IBFCB6wCMEDFiF4wCIED1iE4AGLEDxgEYIHLELwgEUIHrDI\nIwXf19en3NxcrVq1KlzrAeCiRwr+0KFDSk9P/3fTCQCPO8fBt7a2qra2Vps3bx5ylwsAjxfHwW/b\ntk0VFRUaN46nAYAnhaNaT58+rbi4OOXm5nJ2B54gjraaamxsVE1NjWpra3X79m39+eefKikp0fHj\nx//rng0D3vaK3WMBNwTk+u6x9507d07vv/++Pv/88wcPzO6xwBhxYffYgXiWHngyPPIZfsgDc4YH\nxojLZ3gATwaCByxC8IBFHP1a7mn1rnaP9RKeCEbOn6T9X/1PGFeC0eIMD1iE4AGLEDxgEYIHLELw\ngEUIHrAIwQMWIXjAIgQPWITgAYsQPGARggcsQvCARQgesAh/HotRixKXJn9ScYYHLELwgEUcB9/Z\n2amioiKlpaUpPT1dTU1N4VwXABc4fgz/1ltvaeXKlfr000/V29urv//+O5zrAuACR8F3dXXp/Pnz\nqqqquneQ6GhNnjw5rAsDEH6OfqRvaWnR9OnTtWHDBs2bN09btmxRKBQK99oAhJmj4Ht7e+X3+7V1\n61b5/X5NnDhRZWVlg9yzYcAt4HSNAIYV0IOtDc3Rj/Qej0cej0cLFiyQJBUVFQ0R/GInhwcwKl49\nuDPzuSHv6egMn5CQoKSkJDU3N0uS6uvrlZGR4eRQACLI8bP0H3zwgV599VX19PRo1qxZOnr0aDjX\nBcAFjoPPzs7Wd999F861AHAZr7QDLELwgEUIHrAIfx47wG69O9ZLAFzFGR6wCMEDFiF4wCIED1iE\n4AGLEDxgEYIHLELwgEUIHrAIwQMWIXjAIgQPWITgAYsQPGARggcsQvCARQgesIjj4Pft26eMjAxl\nZWWpuLhYd+7cCee6ALjAUfCBQEAfffSR/H6/fvrpJ/X19enkyZPhXhuAMHN0TbtJkyYpJiZGoVBI\nzzzzjEKhkBITE8O9NgBh5ugMP3XqVG3fvl0zZ87UjBkzNGXKFC1btizcawMQZo6Cv3r1qg4ePKhA\nIKD29nZ1d3frxIkTg9yzQeweC7gtoJHuHuso+IsXL2rhwoWaNm2aoqOjVVhYqMbGxkHuuXjAzetk\nFICH8urB1obmKPjU1FQ1NTXp1q1bMsaovr5e6enpTg4FIIIcBZ+dna2SkhLl5eVp7ty5kqTXXnst\nrAsDEH5RxhjjyoGjoiR2cgHGwG4NlTWvtAMsQvCARQgesAjBAxYheMAiBA9YhOABixA8YBGCByxC\n8IBFCB6wCMEDFiF4wCIED1iE4AGLEDxgEYIHLELwgEUIHrAIwQMWGTb4jRs3Kj4+XllZWf0fu3nz\npnw+n1JSUrR8+XJ1dna6vkgA4TFs8Bs2bFBdXd0DHysrK5PP51Nzc7OWLl2qsrIyVxcIIHyGDf7F\nF1/Us88++8DHampqVFpaKkkqLS3VZ5995t7qAITVqB/Dd3R0KD4+XpIUHx+vjo6OsC8KgDse6Um7\nqKiofzecAPAkGPX+8PHx8bp+/boSEhJ07do1xcXFDXPvhgFve8WGkoAbAhrp7syjPsMXFBSoqqpK\nklRVVaXVq1cPc+/FA27e0Y4CMCJehWX32HXr1mnhwoX69ddflZSUpKNHj2rXrl06c+aMUlJSdPbs\nWe3atSs8awbgOjaTBJ46bCYJQAQPWIXgAYsQPGARggcsQvCARQgesAjBAxYheMAiBA9YhOABixA8\nYBGCByxC8IBFCB6wCMEDFiF4wCIED1iE4AGLEDxgEYIHLPLQ4AfbQXbHjh1KS0tTdna2CgsL1dXV\n5eoiAYTHQ4MfbAfZ5cuX6+eff9aPP/6olJQU7du3z7UFAgifhwY/2A6yPp9P48bd+9L8/Hy1tra6\nszoAYfXIj+ErKyu1cuXKcKwFgMseKfg9e/Zo/PjxKi4uDtd6ALho1LvH3nfs2DHV1tbqq6++GuZe\nDQPe9ooNJQE3BDTS3WMdBV9XV6eKigqdO3dOEyZMGOaei50cHsCoePXgyfTckPd86I/0/72DbGVl\npd588011d3fL5/MpNzdXW7dufeQlA3Afu8cCTx12jwUgggesQvCARQgesAjBAxYheMAiBA9YhOAB\nixA8YBGCByxC8IBFCB6wCMEDFiF4wCIED1iE4AGLEDxgEYIHLELwgEUIHrAIwQMWGTb4wXaOvW//\n/v0aN26cbt686driAITXsMEPtnOsJAWDQZ05c0bPP/+8awsDEH7DBj/YzrGS9M477+i9995zbVEA\n3DHqx/DV1dXyeDyaO3euG+sB4KJR7S0XCoW0d+9enTlzpv9jw29c0zDgba/YTBJwQ0CubCZ59epV\nBQIBZWdnS5JaW1s1f/58XbhwQXFxcYN8xeLRHB6AI16NdDPJUQWflZWljo6O/veTk5P1/fffa+rU\nqaNbH4AxMexj+Ps7xzY3NyspKUlHjx594PP3NowE8KQYw91jA4r8Y/qxmMncp3fm4zr3sdw9NmDJ\nTOY+vTOfvLm8tBawyKietButefOeG/Jz7e2xmjFj6M+7YSxmMvfpnfm4zvX7h/lC45JFixYZSdy4\ncYvwbdGiRUN26dqTdgAePzyGByxC8IBFIh58XV2dUlNTNWfOHJWXl0dkZjAY1JIlS5SRkaHMzEwd\nPnw4InMlqa+vT7m5uVq1alXEZnZ2dqqoqEhpaWlKT09XU1NTRObu27dPGRkZysrKUnFxse7cuePK\nnMGu03Dz5k35fD6lpKRo+fLl6uzsjMjcHTt2KC0tTdnZ2SosLFRXV5frM+9zdE0Kt560G0xvb6+Z\nNWuWaWlpMT09PSY7O9tcvnzZ9bnXrl0zly5dMsYY89dff5mUlJSIzDXGmP3795vi4mKzatWqiMwz\nxpiSkhJz5MgRY4wxd+/eNZ2dna7PbGlpMcnJyeb27dvGGGNeeeUVc+zYMVdmffPNN8bv95vMzMz+\nj+3YscOUl5cbY4wpKyszO3fujMjcL7/80vT19RljjNm5c2fY5w420xhjfvvtN7NixQrj9XrNjRs3\nRny8iJ7hL1y4oNmzZ8vr9SomJkZr165VdXW163MTEhKUk5MjSYqNjVVaWpra29tdn9va2qra2lpt\n3rz5IX9VGD5dXV06f/68Nm7cKEmKjo7W5MmTXZ87adIkxcTEKBQKqbe3V6FQSImJia7MGuw6DTU1\nNSotLZUklZaW6rPPPovIXJ/Pp3Hj7mWUn5+v1tZW12dKzq9JEdHg29ralJSU1P++x+NRW1tbJJeg\nQCCgS5cuKT8/3/VZ27ZtU0VFRf83RCS0tLRo+vTp2rBhg+bNm6ctW7YoFAq5Pnfq1Knavn27Zs6c\nqRkzZmjKlClatmyZ63Pv6+joUHx8vCQpPj7+gT/yipTKykqtXLnS9TmPck2KiAY/1n9s093draKi\nIh06dEixsbGuzjp9+rTi4uKUm5sbsbO7JPX29srv92vr1q3y+/2aOHGiysrKXJ979epVHTx4UIFA\nQO3t7eru7taJEydcnzuYqKioiH+v7dmzR+PHj1dxcbGrc+5fk2L37t39HxvN91dEg09MTFQwGOx/\nPxgMyuPxRGT23bt3tWbNGq1fv16rV692fV5jY6NqamqUnJysdevW6ezZsyopKXF9rsfjkcfj0YIF\nCyRJRUVF8g/70qvwuHjxohYuXKhp06YpOjpahYWFamxsdH3uffHx8bp+/bok6dq1a0Ncn8Edx44d\nU21tbUT+gxt4TYrk5OT+a1L8/vvvI/r6iAafl5enK1euKBAIqKenR6dOnVJBQYHrc40x2rRpk9LT\n0/X222+7Pk+S9u7dq2AwqJaWFp08eVIvvfSSjh8/7vrchIQEJSUlqbm5WZJUX1+vjIwM1+empqaq\nqalJt27dkjFG9fX1Sk9Pd33ufQUFBaqqqpIkVVVVReQ/deneb50qKipUXV2tCRMmuD7v/jUpWlpa\n1NLSIo/HI7/fP/L/4ML4hOKI1NbWmpSUFDNr1iyzd+/eiMw8f/68iYqKMtnZ2SYnJ8fk5OSYL774\nIiKzjTGmoaEhos/S//DDDyYvL8/MnTvXvPzyyxF5lt4YY8rLy016errJzMw0JSUlpqenx5U5a9eu\nNc8995yJiYkxHo/HVFZWmhs3bpilS5eaOXPmGJ/PZ/744w/X5x45csTMnj3bzJw5s//76vXXX3dl\n5vjx4/v/rQMlJyeP6ll6XloLWIRX2gEWIXjAIgQPWITgAYsQPGARggcsQvCARQgesMj/AaTfUDG5\nVNnYAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x11007add0>"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# assign K values based on zone\n",
      "hk = np.ones((15, 15), dtype=float) # initialize new array for Kvalues\n",
      "for z, v in Khvalues.iteritems(): hk[Kzones == z] = v # assign K value corresponding to zone for each value in Kzones array\n",
      "\n",
      "# make the lpf object\n",
      "#lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp)\n",
      "upw = flopy.modflow.ModflowUpw(mf, hk=hk, vka=1, sy=sy, ss=ss, laytyp=laytyp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Now make the river cells:####"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# bring in river cell information from csv file\n",
      "rivcells = pd.read_csv('../rivercells.csv')\n",
      "rivcells"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>layer</th>\n",
        "      <th>row</th>\n",
        "      <th>column</th>\n",
        "      <th>stage</th>\n",
        "      <th>Rcond</th>\n",
        "      <th>rbot</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>0 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td>  1</td>\n",
        "      <td> 510.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 506.0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td>  2</td>\n",
        "      <td> 509.5</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 505.5</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td>  3</td>\n",
        "      <td> 509.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 505.0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td>  4</td>\n",
        "      <td> 508.5</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 504.5</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td>  5</td>\n",
        "      <td> 508.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 504.0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>5 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td>  6</td>\n",
        "      <td> 507.5</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 503.5</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>6 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td>  7</td>\n",
        "      <td> 507.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 503.0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>7 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 2</td>\n",
        "      <td>  8</td>\n",
        "      <td> 506.5</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 502.5</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>8 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 3</td>\n",
        "      <td>  9</td>\n",
        "      <td> 506.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 502.0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>9 </th>\n",
        "      <td> 1</td>\n",
        "      <td> 3</td>\n",
        "      <td> 10</td>\n",
        "      <td> 505.5</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 501.5</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>10</th>\n",
        "      <td> 1</td>\n",
        "      <td> 3</td>\n",
        "      <td> 11</td>\n",
        "      <td> 505.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 501.0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>11</th>\n",
        "      <td> 1</td>\n",
        "      <td> 3</td>\n",
        "      <td> 12</td>\n",
        "      <td> 504.5</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 500.5</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>12</th>\n",
        "      <td> 1</td>\n",
        "      <td> 2</td>\n",
        "      <td> 13</td>\n",
        "      <td> 504.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 500.0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>13</th>\n",
        "      <td> 1</td>\n",
        "      <td> 2</td>\n",
        "      <td> 14</td>\n",
        "      <td> 503.5</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 499.5</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>14</th>\n",
        "      <td> 1</td>\n",
        "      <td> 1</td>\n",
        "      <td> 15</td>\n",
        "      <td> 503.0</td>\n",
        "      <td> 150000</td>\n",
        "      <td> 499.0</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 35,
       "text": [
        "    layer  row  column  stage   Rcond   rbot\n",
        "0       1    1       1  510.0  150000  506.0\n",
        "1       1    1       2  509.5  150000  505.5\n",
        "2       1    1       3  509.0  150000  505.0\n",
        "3       1    1       4  508.5  150000  504.5\n",
        "4       1    1       5  508.0  150000  504.0\n",
        "5       1    1       6  507.5  150000  503.5\n",
        "6       1    1       7  507.0  150000  503.0\n",
        "7       1    2       8  506.5  150000  502.5\n",
        "8       1    3       9  506.0  150000  502.0\n",
        "9       1    3      10  505.5  150000  501.5\n",
        "10      1    3      11  505.0  150000  501.0\n",
        "11      1    3      12  504.5  150000  500.5\n",
        "12      1    2      13  504.0  150000  500.0\n",
        "13      1    2      14  503.5  150000  499.5\n",
        "14      1    1      15  503.0  150000  499.0"
       ]
      }
     ],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# update rCond values with single paramter value from top of script\n",
      "rivcells['Rcond'] = Rcond\n",
      "\n",
      "# make dataframe into list for flopy input\n",
      "rivdata = rivcells.values.tolist()\n",
      "\n",
      "# need to copy river cell info for each stress period\n",
      "rivdata = [rivdata for p in perlen]\n",
      "rivdata"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 36,
       "text": [
        "[[[1.0, 1.0, 1.0, 510.0, 150000.0, 506.0],\n",
        "  [1.0, 1.0, 2.0, 509.5, 150000.0, 505.5],\n",
        "  [1.0, 1.0, 3.0, 509.0, 150000.0, 505.0],\n",
        "  [1.0, 1.0, 4.0, 508.5, 150000.0, 504.5],\n",
        "  [1.0, 1.0, 5.0, 508.0, 150000.0, 504.0],\n",
        "  [1.0, 1.0, 6.0, 507.5, 150000.0, 503.5],\n",
        "  [1.0, 1.0, 7.0, 507.0, 150000.0, 503.0],\n",
        "  [1.0, 2.0, 8.0, 506.5, 150000.0, 502.5],\n",
        "  [1.0, 3.0, 9.0, 506.0, 150000.0, 502.0],\n",
        "  [1.0, 3.0, 10.0, 505.5, 150000.0, 501.5],\n",
        "  [1.0, 3.0, 11.0, 505.0, 150000.0, 501.0],\n",
        "  [1.0, 3.0, 12.0, 504.5, 150000.0, 500.5],\n",
        "  [1.0, 2.0, 13.0, 504.0, 150000.0, 500.0],\n",
        "  [1.0, 2.0, 14.0, 503.5, 150000.0, 499.5],\n",
        "  [1.0, 1.0, 15.0, 503.0, 150000.0, 499.0]],\n",
        " [[1.0, 1.0, 1.0, 510.0, 150000.0, 506.0],\n",
        "  [1.0, 1.0, 2.0, 509.5, 150000.0, 505.5],\n",
        "  [1.0, 1.0, 3.0, 509.0, 150000.0, 505.0],\n",
        "  [1.0, 1.0, 4.0, 508.5, 150000.0, 504.5],\n",
        "  [1.0, 1.0, 5.0, 508.0, 150000.0, 504.0],\n",
        "  [1.0, 1.0, 6.0, 507.5, 150000.0, 503.5],\n",
        "  [1.0, 1.0, 7.0, 507.0, 150000.0, 503.0],\n",
        "  [1.0, 2.0, 8.0, 506.5, 150000.0, 502.5],\n",
        "  [1.0, 3.0, 9.0, 506.0, 150000.0, 502.0],\n",
        "  [1.0, 3.0, 10.0, 505.5, 150000.0, 501.5],\n",
        "  [1.0, 3.0, 11.0, 505.0, 150000.0, 501.0],\n",
        "  [1.0, 3.0, 12.0, 504.5, 150000.0, 500.5],\n",
        "  [1.0, 2.0, 13.0, 504.0, 150000.0, 500.0],\n",
        "  [1.0, 2.0, 14.0, 503.5, 150000.0, 499.5],\n",
        "  [1.0, 1.0, 15.0, 503.0, 150000.0, 499.0]],\n",
        " [[1.0, 1.0, 1.0, 510.0, 150000.0, 506.0],\n",
        "  [1.0, 1.0, 2.0, 509.5, 150000.0, 505.5],\n",
        "  [1.0, 1.0, 3.0, 509.0, 150000.0, 505.0],\n",
        "  [1.0, 1.0, 4.0, 508.5, 150000.0, 504.5],\n",
        "  [1.0, 1.0, 5.0, 508.0, 150000.0, 504.0],\n",
        "  [1.0, 1.0, 6.0, 507.5, 150000.0, 503.5],\n",
        "  [1.0, 1.0, 7.0, 507.0, 150000.0, 503.0],\n",
        "  [1.0, 2.0, 8.0, 506.5, 150000.0, 502.5],\n",
        "  [1.0, 3.0, 9.0, 506.0, 150000.0, 502.0],\n",
        "  [1.0, 3.0, 10.0, 505.5, 150000.0, 501.5],\n",
        "  [1.0, 3.0, 11.0, 505.0, 150000.0, 501.0],\n",
        "  [1.0, 3.0, 12.0, 504.5, 150000.0, 500.5],\n",
        "  [1.0, 2.0, 13.0, 504.0, 150000.0, 500.0],\n",
        "  [1.0, 2.0, 14.0, 503.5, 150000.0, 499.5],\n",
        "  [1.0, 1.0, 15.0, 503.0, 150000.0, 499.0]]]"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# make the riv package object\n",
      "riv = flopy.modflow.mfriv.ModflowRiv(mf, irivcb=59, layer_row_column_data=rivdata, extension='riv', unitnumber=18, options=None, naux=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 37
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Make the leaking ditch:####"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# designate flux cells\n",
      "flux_cells = zip(np.ones(15, dtype=int) * 15, np.arange(1, 16))\n",
      "flux_cells\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 38,
       "text": [
        "[(15, 1),\n",
        " (15, 2),\n",
        " (15, 3),\n",
        " (15, 4),\n",
        " (15, 5),\n",
        " (15, 6),\n",
        " (15, 7),\n",
        " (15, 8),\n",
        " (15, 9),\n",
        " (15, 10),\n",
        " (15, 11),\n",
        " (15, 12),\n",
        " (15, 13),\n",
        " (15, 14),\n",
        " (15, 15)]"
       ]
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q = Qleak / len(flux_cells) # flow rate in each constant flux cell"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# now make list of layer, row, column, q info for each pumping cell, for each stress period\n",
      "bflux = [[[1, c[0], c[1], q] for c in flux_cells] for p in perlen]\n",
      "\n",
      "# add pumping well to second (transient) stress period\n",
      "bflux[1].append(pumping_well_info2)\n",
      "\n",
      "# add second well at location M for prediction stress period\n",
      "bflux[2].append(pumping_well_info3)\n",
      "\n",
      "bflux"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 40,
       "text": [
        "[[[1, 15, 1, 3000],\n",
        "  [1, 15, 2, 3000],\n",
        "  [1, 15, 3, 3000],\n",
        "  [1, 15, 4, 3000],\n",
        "  [1, 15, 5, 3000],\n",
        "  [1, 15, 6, 3000],\n",
        "  [1, 15, 7, 3000],\n",
        "  [1, 15, 8, 3000],\n",
        "  [1, 15, 9, 3000],\n",
        "  [1, 15, 10, 3000],\n",
        "  [1, 15, 11, 3000],\n",
        "  [1, 15, 12, 3000],\n",
        "  [1, 15, 13, 3000],\n",
        "  [1, 15, 14, 3000],\n",
        "  [1, 15, 15, 3000]],\n",
        " [[1, 15, 1, 3000],\n",
        "  [1, 15, 2, 3000],\n",
        "  [1, 15, 3, 3000],\n",
        "  [1, 15, 4, 3000],\n",
        "  [1, 15, 5, 3000],\n",
        "  [1, 15, 6, 3000],\n",
        "  [1, 15, 7, 3000],\n",
        "  [1, 15, 8, 3000],\n",
        "  [1, 15, 9, 3000],\n",
        "  [1, 15, 10, 3000],\n",
        "  [1, 15, 11, 3000],\n",
        "  [1, 15, 12, 3000],\n",
        "  [1, 15, 13, 3000],\n",
        "  [1, 15, 14, 3000],\n",
        "  [1, 15, 15, 3000],\n",
        "  [1, 7, 11, -20000]],\n",
        " [[1, 15, 1, 3000],\n",
        "  [1, 15, 2, 3000],\n",
        "  [1, 15, 3, 3000],\n",
        "  [1, 15, 4, 3000],\n",
        "  [1, 15, 5, 3000],\n",
        "  [1, 15, 6, 3000],\n",
        "  [1, 15, 7, 3000],\n",
        "  [1, 15, 8, 3000],\n",
        "  [1, 15, 9, 3000],\n",
        "  [1, 15, 10, 3000],\n",
        "  [1, 15, 11, 3000],\n",
        "  [1, 15, 12, 3000],\n",
        "  [1, 15, 13, 3000],\n",
        "  [1, 15, 14, 3000],\n",
        "  [1, 15, 15, 3000],\n",
        "  [1, 9, 9, -30000]]]"
       ]
      }
     ],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# create the well package\n",
      "wel = flopy.modflow.ModflowWel(mf, layer_row_column_Q=bflux)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####run the model####"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#write the model input files\n",
      "mf.write_input()\n",
      "\n",
      "\n",
      "# manually append info for writing out riv package results to nam file (couldn't figure out how to make flopy do this)\n",
      "nam = open('{}.nam'.format(modelname),'a')\n",
      "nam.write('DATA(BINARY)  59 {}.rivout REPLACE\\n'.format(modelname))\n",
      "nam.close()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 42
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Show everything####"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# show the observation locations with the K zones, River cells, and pumping wells\n",
      "fig = plt.figure()\n",
      "ax = fig.add_subplot(111)\n",
      "extent = [0, 1500, 0, 1500]\n",
      "\n",
      "# add the river cells and wells to the Kzones array so we can verify where they are\n",
      "for i in range(len(rivcells)):\n",
      "    Kzones[rivcells.row[i] - 1, rivcells.column[i] - 1] = 0\n",
      "for w in bflux[1]:\n",
      "    r, c = w[1:3]\n",
      "    Kzones[r - 1, c - 1] = 3\n",
      "ax.imshow(Kzones, extent=extent, interpolation='None')\n",
      "\n",
      "# bring in observation info\n",
      "obs = pd.read_csv('../observations.csv')\n",
      "ax.scatter(obs.X, obs.Y, c='w')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 43,
       "text": [
        "<matplotlib.collections.PathCollection at 0x1101b9e10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQkAAAD7CAYAAAB5RWHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHb5JREFUeJzt3XtYVPW6B/DvDBcBL2gqEAwmIncQBNG0i5qiZVqGpqkH\nvGQ3dzd3oe5ummfDDJmnNPM8z362aFvblud00d1W8pZuU0MUb4GJ2aCAwFGQ+22A3/mD3RipS1jM\nzJoZvp/n4XlizeV9IebrmjVr/V6VEEKAiOg21Eo3QETWjSFBRJIYEkQkiSFBRJIYEkQkiSFBRJIc\nlW7gt8aMGYODBw8q3QZRlzN69GgcOHDglrdZ1Z7EwYMHIYS449fy5cvbdT9TfLEWa3WFWlL/OFtV\nSBCR9WFIEJEkmwyJMWPGsBZrsZaFqIQQVnPthkqlghW1Q9RlSL32bHJPgogshyFBRJIYEkQkiSFB\nRJKs6ozL9lKpVijdgu1TqeQ97uxy0/bRFck5Nh+xosMPGTDAHZcuLZZRrC3uSRCRJIYEEUliSBCR\nJIYEEUliSBCRJIYEEUliSBCRJIYEEUmSDIkFCxbA09MTERERN922evVqqNVqlJWVGbdptVoEBAQg\nODgYu3fvNm4/ceIEIiIiEBAQgFdeecWE7RORuUmGxPz585Genn7T9vz8fOzZswf33HOPcVtOTg4+\n//xz5OTkID09HYsWLTJeevrCCy9gw4YNuHDhAi5cuHDL5yQi6yQZEg888AD69Olz0/Y//vGPeO+9\n99ps2759O2bNmgUnJycMHDgQgwcPRkZGBoqKilBVVYXhw4cDABITE/H111+b8EcgInPq8DGJ7du3\nQ6PRYMiQIW22X7lyBRqNxvi9RqNBYWHhTdt9fHxQWFjYiZaJyJI6dIFXbW0tUlJSsGfPHuM2U68k\ntWLFCuN/jxkz5tbLdp1dcfM26hi5/9+GvNvxx5yxgYvC5Pw65PwuAEDOtXVy/uadbn/TgQMHbruE\n/u91KCQuXryIvLw8REZGAgAKCgoQExODjIwM+Pj4ID8/33jfgoICaDQa+Pj4oKCgoM12Hx+f29b4\nbUgQkXn8/h/gd9+9feB16O1GREQESkpKoNfrodfrodFokJWVBU9PTzz22GP47LPP0NjYCL1ejwsX\nLmD48OHw8vJCr169kJGRASEENm/ejKlTp8r+4YjIsiRDYtasWRg1ahRyc3Ph6+uLjRs3trld9Zs1\nCUJDQzFjxgyEhobikUcewfr16423r1+/HgsXLkRAQAAGDx6Mhx9+2Aw/ChGZg02ulq360QLN2Du5\n/9ujVnb8MTwm0ZacYxKnO/47HOAEXApq3325WjYRycaQICJJDAkiksSQICJJDAkiksSQICJJDAki\nksSQICJJDAkikmSbZ1xmW6AZurUWGX8ukTLO0rQ0C50FKZuM/gY4AZcC2/n0POOSiORiSBCRJIaE\nJVWWA5d/BgyNSndC1G4MCQtx2PwhnCcMQL/n4+A6JRC4wEtZyTYwJCzh7DH02vw+fs7JxtVLeqxb\nuQJuSTOU7oqoXRgSlnD+DMaPj4Ovry8AYN7cRNTrc/m2g2wCQ8ISfAfh8OHDqKysBADs3bsX3Tzu\nBpycFW6M6M46tBAuyTR8LErvexQDQ0Ix0D8AP+Vko27VNqW7ImqXDo/5S0pKQkhICCIjIxEfH4+K\nigrjbRzzdxsqFRqWfIDrH6fjZMIbqPvyR2D4GKW7ImqXDo/5mzBhArKzs3H69GkEBgZCq9UC4Ji/\ndgkIB0bFAX09lO6EqN06POYvLi4OanXrw0aMGGGcqcExf0T2qVMHLtPS0jBp0iQAHPNHZK9kH7hM\nTk6Gs7MzZs+ebcp+2jfmj5SjlnGlkS0sqS+HnIvCrITZxvz9atOmTdi5cyf27dtn3MYxf0S2w2xj\n/gAgPT0dq1atwvbt2+Hi4mLczjF/RPZJck9i1qxZOHjwIK5duwZfX1+8++670Gq1aGxsRFxcHABg\n5MiRWL9+fZsxf46OjjeN+Zs3bx7q6uowadIkjvkjsiFcdIbMz2r+wkzMyo9JcNEZIrIIhgQRSWJI\nEJEkhgQRSWJIEJEkhgQRSWJIEJEkhgQRSWJIEJEkLl9H5mflZyaSNO5JEJEkhgQRSWJIUOcZDMC1\nEqC5WelObq+qovWLOowhQZ3z3Q50G+uJ7vFhcH14IHA2U+mO2jI0wnXZHDiP84HzOB+4LpnFoUgd\nxJAg+YoL4Lb8aRza/S2qS69hy7o1cF38hFW9CB03pGK44TrKr11FRek13Cuq4PiXZKXbsikMCZLv\n5x8RFhmF2NhYAEB8fDzcHB2A4oI7PNBy3LIzsPiF5+Dq6goXFxf8cdHzcMs+pnRbNoUhQfJ5+SL3\nXDZKS0sBAOfPn0dVRblVzRVp8LoHe747aPx+z3cH0eg1QMGObA/PkyD5Boeh7omFCBwShcjoGGRm\nHIVh2VrArYfSnRk1PL8cGxeMxuETJ6BSqXC+sBj1aQfv/EAykly+bsGCBfjnP/8JDw8PnD17FgBQ\nVlaGmTNn4tKlSxg4cCC2bduG3r17A2gd85eWlgYHBwesXbsWEyZMANA65m/evHmor6/HpEmTsGbN\nmls3w+XrbFNOFlDwCxAQAfgFKd3NzWprgMwDrf8dO9qqQsycLLJ83a3G/Ol0OsTFxSE3Nxfjxo2D\nTqcDwDF/XVpoNDBhunUGBAC4dQdGP9r61UUCwpQ6POZvx44dmDt3LgBg7ty5xpF9HPNHZJ86fOCy\npKQEnp6eAABPT0+UlJQA4Jg/InvVqQOXKpXKOFvDVKxtzJ9K5nrwzX0cTNyJ7RIyr/ByuG7FZ3Da\nOLOO+fP09ERxcTG8vLxQVFQED4/Wj7s45o/Idph1zN9jjz2GTz75BADwySefGEf2ccwfkX3q0Ji/\nlStXYtmyZZgxYwY2bNhg/AgUAMf8Edkpjvm7Uy0THZOoqhZY/4kKxVcd8eAIA554pOusxMJjEsrg\nmD8bUlcnMPZJN5y++Dg0gSl44z0NUj/mgU2yDdyTuFMtE+xJbNsh8Jdtsdi7LwMqlQr5+fkIDh6E\nyvNNUKvtf4+CexLK4J6EDamtAzw9vYzHaDw8PGAwtFj1Gi1Ev2JIWMC4+4G9e/di8+bNyMnJwdML\n5uDR8S5wcrL/vQiyfXy7cadaJjpwefy0wOv/2RNF/wc8OKIJ/7W8Dj17dI2Q4NsNZZjq7QZD4k61\neMZlpzEklMFjEkRkEQwJIpLEkCAiSVy+zkzkHsug35BxuOydSHl/0ivP8PjH7XBPgogkMSSISBJD\noqUFaKhXugubJgRQV6d0F9KahUCz9Xzaf2v1dbLeYplblw4J1ZdpcBrVGw4j3dF93gPAtWKlW7I5\ne/8FeMe4wD3EAcGj3fDjT0p31FaLEPjW2RnvqdV4T61GurMzWqzthfhzNtweD4HDSHe4PHQ3cHi3\n0h210XVD4vQPcP/vt3H2eCYa6+rw/Lj70f2NBKW7silFJcCsF13x96070dBgwJ/eXIcp811hMCjd\n2Q0ZDg5wiIjA1dJSlF6/jm5RUTjqaEXH65ua4PbyY1jzpyQYGhrw7f9ug+uf5gAl1rMObBcOiaOY\nMW0agoKCoFarseKtN9GQ9b3SXdmUM+eAyCFhGDt2LFQqFebOm4/mFlcUFCnd2Q1Fbm5IeuMNuLu7\no2fPnkh64w0Uu7oq3dYNV6/A2VCPhU8vgEqlwoMPPogh0THAT6eU7syo64ZEPy9knMhC878vxczM\nzIRz/7sVbsq2ePUHzuf+gsrKSgBAXl4eyiuq0bfPHR5oQW4GA44ePmz8/ocjR+DW1KRgR7/j3he1\nlZX45ZdfAABVVVX4+adzgBX9LVrRfpeFTXgSF77ZjLDYEQgKDsbe3d+iNnmz0l3ZlMgwYPqkGsRE\nB2Pkvfdi77790C1rQq+eSnd2w6j6eqR9/DFOZmZCpVIhKyMDc+rrAROv8i6bW3c0/XEVokfdj4fG\nj0dGRgaqxzzeOvDISsi+wEur1WLLli1Qq9WIiIjAxo0bUVNT0+ERgG2asfQFXs3NwKFdQPk1IPp+\nYMDgm2vJPCmqpU/X2Uk7eBS4eAmICgOiI26+Xe4FXuqyjp/gdKuTqeqEwMV//7c/ANdbBITiJ1Od\nO9n6FsPHr3UUoQlCTNGrQPPy8vDQQw/h3Llz6NatG2bOnIlJkyYhOzsb/fr1w5IlS5Camorr169D\np9MhJycHs2fPRmZmJgoLCzF+/Hjk5uZCrW77QrKnq0C7UkjcidIh0R6Kh4QZKHoVaK9eveDk5ITa\n2lo0NTWhtrYW3t7eHRoBeOzYMTmlicjCZIXEXXfdhddeew0DBgyAt7c3evfujbi4uA6PACQi6ydr\n3+zixYv48MMPkZeXB3d3dzz55JPYsmVLm/vcaQTg7W6ztjF/sneVr7eYuJMuSMb78pWnreiTCytm\n1jF/AHD8+HGMGjUKffv2BQDEx8fj6NGj8PLyavcIwNuN+uOYPyLzM+uYPwAIDg7GDz/8gLq6Oggh\nsHfvXoSGhmLKlCkdGgFIRNZP1p5EZGQkEhMTMWzYMKjVakRHR+PZZ59FVVVVh0cAEpF140K4ZsJF\nZzpP1vGgW/39FOiB9M9ab3v4KcB30M33scN/tLgQLlF7XMyBa8JIPN9chBdECVwSRgIXflS6K5vS\ndU/Lpi7BbWMq3kl6DUuXLAEA+A8ciHfSdKjVbrnDI+lX3JMgu+ZYXQH/QTfeXvgPGgTH6goFO7I9\n3JMgu1b1wGQsXbHSuCTA628vR/UTLyjdlk1hSJBdE/FP41JlGUY++hggBOqnPYuWJ59Vui2bwk83\nzISfbnSeyT7daA9+usFPN4hIHoYEEUliSBCRJB64NBO5V49SJ9nhsQWlcU+CiCQxJIhIEkOC2rLk\nJ+JCWOVYO2qLIUGtGhvgsmIhHGO7o9v9d8Fh4yrz1WpqQjftS3Aa0RPOo3rD8ePlDAsrxpAgAIDz\nR29hZE0x/u9KIXJOHMfdOzYA6dvMUsvpr1pEFWSjME+Pn7N/hN/hHVB9scEstajzGBIEAHA5tg+6\nFe+gT58+GDRoEP706stwObbPLLXcMvch+e030b9/f/j6+uKdpNfQPdM8tajzGBIEAGjp0x9nzt5Y\nZ+HEmbMw9O5nllpNffrj1Jmzxu9PnjmLRnfz1KLO47Ub1Cr7BFwXPYKpU6bgWlkZDp/JRu3mI0Af\nM7x4L56Dy8KxmDxxIurq67H/6A+o+9sRwPPWiyOTPIpfu1FeXo7p06cjJCQEoaGhyMjIQFlZGeLi\n4hAYGIgJEyagvLzceH+tVouAgAAEBwdj9+7dcsuSuYTFoO7vx7BVMwx7hk1B7dZM8wQEAPiHoP6z\nE/hfv5H4Z/h41H2exYCwYrL3JObOnYvRo0djwYIFaGpqQk1NDZKTk+1uzB+RrVJ0T6KiogKHDh3C\nggULAACOjo5wd3fnmD8iOyQrJPR6Pfr374/58+cjOjoazzzzDGpqajjmj8gOybrAq6mpCVlZWVi3\nbh1iY2Px6quvQqfTtbmPvYz5I7JHZh/zp9FooNFoEBsbCwCYPn06tFotx/wR2Qizj/nz8vKCr68v\ncnNzAQB79+5FWFgYx/wR2SHZ60l89NFHmDNnDhobG+Hv74+NGzeiubmZY/6I7AxPprIUgwH4dhtw\nrRgYeh8Qea/SHXUdV4uAvV+2XkQ2Ph7w8Fa6o5tl7AdysgCfga09qjt/MrTiJ1NRBzQ1we2lyYj+\n5i94vjEfvV+fBtXXm5Tuqmu4fBGuM2Mw43ImZhacgMtTMcClC0p31YbjhlR4rFyIRYZCBG3WwfWt\nuVZ1VSz3JCxh/3aEbtbizNHDcHBwwLlz5xA54l4YjpRzuTUzc317PpZFD8Y7b70JAEjW6pCSkYPa\n5L8p3Nm/1VTB+SFv6HPPw9vbG/X19fALDUdx8qfAkBGdemruSdiS8lIEBwXBwcEBABAYGIiW+jqg\nyaBwY/bPqaIU4aEhxu8jwkLhWFGqYEe/U1kOlx494O3d+hbIxcUF9wzyB8qtp0eGhCVE34/0XTvx\n3XffobKyEq8tXQaX6PsAJ2elO7N71fdOwNspOly6dAn5+fl4888pqBkRp3RbN3h4o6nXXdC9twpV\nVVX46quvcOb0KSA0RunOjPh2w1IO7YKb9kU0Xi1Gt2EPouY/NwH9PJXuyv61tMDp43eg/vy/AQi0\nTH8OhpeTTXJg0GQK9OjxViLqfzyBbr5+qFn+VyBqZKef1lRvNxgSRHaKxySIyCIYEkQkiSFBRJIY\nEkQkiSFBRJIYEkQkiSFBRJIYEtQ5506ix6xYdHugL7o/Ox4ouqx0R2RiDAmSr7wUrn+YhHWvv4zL\nP53D6w+PgdsfHgWam5XujEyIIUHyZZ9ASHAI5iYmwMPDA8vfehPOVWVASYHSnZEJMSRIvl69cSX/\nMhoaGgAA165dQ21lJdC9l8KNkSnJXr6OCOGxqAiOQezosXhk7Bhs/fIriMTFgHsfpTsjE+rUnkRz\nczOGDh2KKVOmAADH/HU1KhXqdH/H2Sf+gPeqXZH/8ioY/rBS6a7IxDoVEmvWrEFoaKhxUVudToe4\nuDjk5uZi3LhxxlkcOTk5+Pzzz5GTk4P09HQsWrQILS0tne+elOfgAEyeAzz/NjBmstLdkBnIDomC\nggLs3LkTCxcuNF5iyjF/RPZHdkgsXrwYq1atajP0l2P+iOyPrJD45ptv4OHhgaFDh95+NRuZY/6I\nyLrI+nTjyJEj2LFjB3bu3In6+npUVlYiISEBnp6eJh3zx1mgRObRkVmgnV6+7uDBg3j//ffxj3/8\nA0uWLEHfvn2xdOlS6HQ6lJeXQ6fTIScnB7Nnz8axY8dQWFiI8ePH4+eff75pb4LL1xGZjqmWrzPJ\neRK/vtiXLVvGMX9EdoYL4RLZKS6ES0QWwZAgIkkMCSKSxJAgIkkMCSKSxJAgIkkMCSKSxJAgIkkM\nCSKSxJAgIkkMCSKSxJAgIkkMCSKSxJAg+yYE1J9+hJ5TQ9FzagjUW9YA1nPhs03g3A2ya6odf4PP\nF+vxxd//BpVKhWn/kYh8t54Q8QuUbs1mcE+C7FrP777C6j+vRGxsLIYNG4YPUv6Mnt99pXRbNoUh\nQXatybUHCgpurMxeUFCIZrceCnZke/h2g+xa7bwleOOZcSi4cgVqtRof/3UD6v6yR+m2bIqsPYn8\n/HyMHTsWYWFhCA8Px9q1awFwzB9ZoaAhqP/ke3xQ7YLVFU6o23QICI5SuiubImuNy+LiYhQXFyMq\nKgrV1dWIiYnB119/jY0bN6Jfv35YsmQJUlNTcf369TarZWdmZhpXy87NzW0z2AfgGpdEpqToGpde\nXl6IimpN4x49eiAkJASFhYUc80dkhzp94DIvLw8nT57EiBEjOOaPyA51KiSqq6sxbdo0rFmzBj17\n9mxzG8f8EdkH2Z9uGAwGTJs2DQkJCZg6dSoAcMwfkY0w+5g/IQTmzp2Lvn374oMPPjBu55g/Iuth\nqgOXskLi+++/x4MPPoghQ4YYX+harRbDhw/HjBkzcPnyZeOYv969ewMAUlJSkJaWBkdHR6xZswYT\nJ07sUKNt7seQILojRUPCXCweEi0twLHvgOvXgKiRwN0DTPTERMrjmL/Oam6G2+tPYuB/vYrx32+D\n68yY1sAgoja67mnZe76A3/VCnDyVBScnJ3z77beY9vyzqPnmgtKdEVmVrrsnUVKA++8dAScnJwDA\nfffdh4biAoWbIrI+XTckIkbgf774Enq9HkIIpL6/Gt0i71W6KyKr03XfbkTfh/L5yxAYFg4HtRqO\nfkGo+fBrpbsisjpd+9MNADAYgLoaoKc7wLNAyY6Y6tONrrsn8SsnJ8Cpt9JdEFmtrntMgojahSFB\nRJJs8u3G8ghmG9GduA8YAOTldfp5+GojIkkMCSKSxJAgIkkMCSKSxJAgIkkMCSKSxJAgIkkMCSKS\nZNGQSE9PR3BwMAICApCammrJ0kQkk8VCorm5GS+++CLS09ORk5ODrVu34ty5c5YqT0QyWSwkjh07\nhsGDB2PgwIFwcnLCU089he3bt1uqPBHJZLGQKCwshK+vr/F7jvojsg0Wu8DLlGP9wmbONNlzEdkr\nt759TfI8FguJ34/6y8/PbzNE+FftGfP35Nat5miRqMsw+5g/OZqamhAUFIR9+/bB29sbw4cPx9at\nWxESEnKjmXYuX0dEpmUVy9c5Ojpi3bp1mDhxIpqbm/H000+3CQgisk42uRAuEZkWx/wRkWwMCSKS\nZJMh0d6jsqzFWqzVeQwJ1mIt1pJkkyFBRJbDkCAiacKKjB49WgDgF7/4ZeGv0aNH3/Z1aVXnSRCR\n9eHbDSKSxJAgIkk2FRKmXv4uPz8fY8eORVhYGMLDw7F27VoAQFlZGeLi4hAYGIgJEyagvLzc+Bit\nVouAgAAEBwdj9+7dHa7Z3NyMoUOHYsqUKWatVV5ejunTpyMkJAShoaHIyMgwWy2tVouwsDBERERg\n9uzZaGhoMFmtBQsWwNPTExEREcZtcp77xIkTiIiIQEBAAF555ZV210pKSkJISAgiIyMRHx+PiooK\nk9S6Xb1frV69Gmq1GmVlZSarJ5ulDkp2VlNTk/D39xd6vV40NjaKyMhIkZOT06nnLCoqEidPnhRC\nCFFVVSUCAwNFTk6OSEpKEqmpqUIIIXQ6nVi6dKkQQojs7GwRGRkpGhsbhV6vF/7+/qK5ublDNVev\nXi1mz54tpkyZIoQQZquVmJgoNmzYIIQQwmAwiPLycrPU0uv1ws/PT9TX1wshhJgxY4bYtGmTyWr9\n61//EllZWSI8PNy4rSPP3dLSIoQQIjY2VmRkZAghhHjkkUfErl272lVr9+7dxv6WLl1qslq3qyeE\nEJcvXxYTJ04UAwcOFKWlpSarJ5fNhMSRI0fExIkTjd9rtVqh1WpNWuPxxx8Xe/bsEUFBQaK4uFgI\n0RokQUFBQgghUlJShE6nM95/4sSJ4ujRo+1+/vz8fDFu3Dixf/9+MXnyZCGEMEut8vJy4efnd9N2\nc9QqLS0VgYGBoqysTBgMBjF58mSxe/duk9bS6/VtXkgdfe4rV66I4OBg4/atW7eK5557rl21fuvL\nL78Uc+bMMVmt29WbPn26OH36dJuQMFU9OWzm7Ya5l7/Ly8vDyZMnMWLECJSUlMDT0xMA4OnpiZKS\nEgDAlStX2iyU09EeFi9ejFWrVkGtvvFrN0ctvV6P/v37Y/78+YiOjsYzzzyDmpoas9S666678Npr\nr2HAgAHw9vZG7969ERcXZ7bfIdDx39nvt/v4+Mj620lLS8OkSZPMWmv79u3QaDQYMmRIm+3m/tmk\n2ExImHL5u9+rrq7GtGnTsGbNGvTs2fOmulK129vXN998Aw8PDwwdOvS2l+SaqlZTUxOysrKwaNEi\nZGVloXv37tDpdGapdfHiRXz44YfIy8vDlStXUF1djS1btpil1u0ea86/jV8lJyfD2dkZs2fPNluN\n2tpapKSk4N133zVuu93fiiXZTEi0d/m7jjIYDJg2bRoSEhIwdepUAK3/OhUXFwMAioqK4OHhccse\nCgoK4OPj0646R44cwY4dO+Dn54dZs2Zh//79SEhIMEstjUYDjUaD2NhYAMD06dORlZUFLy8vk9c6\nfvw4Ro0ahb59+8LR0RHx8fE4evSoWWr9qiO/M41GAx8fHxQUFMiuuWnTJuzcuROffvqpcZs5al28\neBF5eXmIjIyEn58fCgoKEBMTg5KSErP9bO1i0jcvZmQwGMSgQYOEXq8XDQ0NJjlw2dLSIhISEsSr\nr77aZntSUpLx/Z9Wq73pYFVDQ4P45ZdfxKBBg4wHjzriwIEDxmMS5qr1wAMPiPPnzwshhFi+fLlI\nSkoyS61Tp06JsLAwUVtbK1paWkRiYqJYt26dSWv9/n27nOcePny4+OGHH0RLS4vkwb3f19q1a5cI\nDQ0VV69ebXM/U9S6Vb3futWBy87Wk8NmQkIIIXbu3CkCAwOFv7+/SElJ6fTzHTp0SKhUKhEZGSmi\noqJEVFSU2LVrlygtLRXjxo0TAQEBIi4uTly/ft34mOTkZOHv7y+CgoJEenq6rLoHDhwwfrphrlqn\nTp0Sw4YNE0OGDBFPPPGEKC8vN1ut1NRUERoaKsLDw0ViYqJobGw0Wa2nnnpK3H333cLJyUloNBqR\nlpYm67mPHz8uwsPDhb+/v3jppZfaVWvDhg1i8ODBYsCAAca/jxdeeMEktX5bz9nZ2fiz/Zafn58x\nJExRTy6elk1EkmzmmAQRKYMhQUSSGBJEJIkhQUSSGBJEJIkhQUSSGBJEJIkhQUSS/h9uzDAx5pUu\n1QAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10fe8b710>"
       ]
      }
     ],
     "prompt_number": 43
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Run the model####\n",
      "if using a PC, could just do this:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "    #run the model\n",
      "    success, mfoutput = mf.run_model2(silent=True, pause=False)\n",
      "    if not success:\n",
      "        raise Exception('MODFLOW did not terminate normally.')\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I couldn't get this to work on the mac, though, so I just ran from the terminal:  \n",
      "\n",
      "    $ wine mf2005.exe p9\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 43
    }
   ],
   "metadata": {}
  }
 ]
}