examples/Notebooks/flopy3_PEST.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FloPy\n",
    "## Parameter Estimation with FloPy\n",
    "This notebook demonstrates the current parameter estimation functionality that is available with FloPy.  The capability to write a simple template file for PEST is the only capability implemented so far.  The plan is to develop functionality for creating PEST instruction files as well as the PEST control file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import flopy\n",
    "import flopy.pest.templatewriter as tplwriter\n",
    "import flopy.pest.params as params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook will work with a simple model using the dimensions below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define the model dimensions\n",
    "nlay = 3\n",
    "nrow = 20\n",
    "ncol = 20\n",
    "\n",
    "# Create the flopy model object and add the dis and lpf packages\n",
    "m = flopy.modflow.Modflow(modelname='mymodel', model_ws='./data')\n",
    "dis = flopy.modflow.ModflowDis(m, nlay, nrow, ncol)\n",
    "lpf = flopy.modflow.ModflowLpf(m, hk=10.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple One Parameter Example\n",
    "In order to create a PEST template file, we first need to define a parameter.  For example, let's say we want to parameterize hydraulic conductivity, which is a static variable in flopy and MODFLOW.  As a first step, let's define a parameter called HK_LAYER_1 and assign it to all of layer 1.  We will not parameterize hydraulic conductivity for layers 2 and 3 and instead leave HK at its value of 10. (as assigned in the block above this one). We can do this as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mfpackage = 'lpf'\n",
    "partype = 'hk'\n",
    "parname = 'HK_LAYER_1'\n",
    "idx = np.empty((nlay, nrow, ncol), dtype=np.bool)\n",
    "idx[0] = True\n",
    "idx[1:] = False\n",
    "\n",
    "# The span variable defines how the parameter spans the package\n",
    "span = {'idx': idx}\n",
    "\n",
    "# These parameters have not affect yet, but may in the future\n",
    "startvalue = 10.\n",
    "lbound = 0.001\n",
    "ubound = 1000.\n",
    "transform='log'\n",
    "\n",
    "p = params.Params(mfpackage, partype, parname, startvalue, \n",
    "                  lbound, ubound, span)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, we have enough information to the write a PEST template file for the LPF package.  We can do this using the following statement:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The following packages will be parameterized: ['LPF']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "tw = tplwriter.TemplateWriter(m, [p])\n",
    "tw.write_template()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, the lpf template file will have been created.  The following block will print the template file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ptf ~\n",
      "# LPF for MODFLOW, generated by Flopy.\n",
      "53    -1E+30         0\n",
      "0         0         0\n",
      "0         0         0\n",
      "1.000000E+00   1.000000E+00   1.000000E+00\n",
      "0         0         0\n",
      "0         0         0\n",
      "CONSTANT ~ HK_LAYER_1  ~    #hk Layer 1\n",
      "CONSTANT    1.000000E+00                           #vka layer 1\n",
      "CONSTANT 10.0    #hk Layer 2\n",
      "CONSTANT    1.000000E+00                           #vka layer 2\n",
      "CONSTANT 10.0    #hk Layer 3\n",
      "CONSTANT    1.000000E+00                           #vka layer 3\n"
     ]
    }
   ],
   "source": [
    "lines = open('./data/mymodel.lpf.tpl', 'r').readlines()\n",
    "for l in lines:\n",
    "    print(l.strip())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The span variable will also accept 'layers', in which the parameter applies to the list of layers, as shown next.  When 'layers' is specifed in the span dictionary, then the original hk value of 10. remains in the array, and the multiplier is specified on the array control line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The following packages will be parameterized: ['LPF']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "mfpackage = 'lpf'\n",
    "partype = 'hk'\n",
    "parname = 'HK_LAYER_1-3'\n",
    "\n",
    "# Span indicates that the hk parameter applies as a multiplier to layers 0 and 2 (MODFLOW layers 1 and 3)\n",
    "span = {'layers': [0, 2]}\n",
    "\n",
    "# These parameters have not affect yet, but may in the future\n",
    "startvalue = 10.\n",
    "lbound = 0.001\n",
    "ubound = 1000.\n",
    "transform='log'\n",
    "\n",
    "p = params.Params(mfpackage, partype, parname, startvalue, \n",
    "                  lbound, ubound, span)\n",
    "tw = tplwriter.TemplateWriter(m, [p])\n",
    "tw.write_template()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ptf ~\n",
      "# LPF for MODFLOW, generated by Flopy.\n",
      "53    -1E+30         0\n",
      "0         0         0\n",
      "0         0         0\n",
      "1.000000E+00   1.000000E+00   1.000000E+00\n",
      "0         0         0\n",
      "0         0         0\n",
      "INTERNAL ~ HK_LAYER_1-3  ~ (FREE) -1      #hk Layer 1\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "CONSTANT    1.000000E+00                           #vka layer 1\n",
      "CONSTANT 10.0    #hk Layer 2\n",
      "CONSTANT    1.000000E+00                           #vka layer 2\n",
      "INTERNAL ~ HK_LAYER_1-3  ~ (FREE) -1      #hk Layer 3\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "CONSTANT    1.000000E+00                           #vka layer 3\n"
     ]
    }
   ],
   "source": [
    "lines = open('./data/mymodel.lpf.tpl', 'r').readlines()\n",
    "for l in lines:\n",
    "    print(l.strip())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple Parameter Zoned Approach\n",
    "\n",
    "The params module has a helper function called zonearray2params that will take a zone array and some other information and create a list of parameters, which can then be passed to the template writer.  This next example shows how to create a slightly more complicated LPF template file in which both HK and VKA are parameterized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a zone array\n",
    "zonearray = np.ones((nlay, nrow, ncol), dtype=int)\n",
    "zonearray[0, 10:, 7:] = 2\n",
    "zonearray[0, 15:, 9:] = 3\n",
    "zonearray[1] = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create a list of parameters for HK\n",
    "mfpackage = 'lpf'\n",
    "parzones = [2, 3, 4]\n",
    "parvals = [56.777, 78.999, 99.]\n",
    "lbound = 5\n",
    "ubound = 500\n",
    "transform = 'log'\n",
    "plisthk = params.zonearray2params(mfpackage, 'hk', parzones, lbound, \n",
    "                                  ubound, parvals, transform, zonearray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, Flopy will create three parameters: hk_2, hk_3, and hk_4, which will apply to the horizontal hydraulic conductivity for cells in zones 2, 3, and 4, respectively.  Only those zone numbers listed in parzones will be parameterized.  For example, many cells in zonearray have a value of 1.  Those cells will not be parameterized.  Instead, their hydraulic conductivity values will remain fixed at the value that was specified when the Flopy LPF package was created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a list of parameters for VKA\n",
    "parzones = [1, 2]\n",
    "parvals = [0.001, 0.0005]\n",
    "zonearray = np.ones((nlay, nrow, ncol), dtype=int)\n",
    "zonearray[1] = 2\n",
    "plistvk = params.zonearray2params(mfpackage, 'vka', parzones, lbound, \n",
    "                                  ubound, parvals, transform, zonearray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hk_2 lpf 56.777\n",
      "hk_3 lpf 78.999\n",
      "hk_4 lpf 99.0\n",
      "vka_1 lpf 0.001\n",
      "vka_2 lpf 0.0005\n"
     ]
    }
   ],
   "source": [
    "# Combine the HK and VKA parameters together\n",
    "plist = plisthk + plistvk\n",
    "for p in plist:\n",
    "    print(p.name, p.mfpackage, p.startvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The following packages will be parameterized: ['LPF']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Write the template file\n",
    "tw = tplwriter.TemplateWriter(m, plist)\n",
    "tw.write_template()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ptf ~\n",
      "# LPF for MODFLOW, generated by Flopy.\n",
      "53    -1E+30         0\n",
      "0         0         0\n",
      "0         0         0\n",
      "1.000000E+00   1.000000E+00   1.000000E+00\n",
      "0         0         0\n",
      "0         0         0\n",
      "INTERNAL 1.0 (FREE) -1      #hk Layer 1\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0            10.0\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~ ~    hk_2     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_3     ~\n",
      "~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_3     ~\n",
      "~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_3     ~\n",
      "~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_3     ~\n",
      "~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~\n",
      "10.0            10.0            10.0            10.0            10.0            10.0            10.0 ~    hk_2     ~ ~    hk_2     ~ ~    hk_3     ~\n",
      "~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~ ~    hk_3     ~\n",
      "CONSTANT ~    vka_1    ~    #vka Layer 1\n",
      "CONSTANT ~    hk_4     ~    #hk Layer 2\n",
      "CONSTANT ~    vka_2    ~    #vka Layer 2\n",
      "CONSTANT 10.0    #hk Layer 3\n",
      "CONSTANT ~    vka_1    ~    #vka Layer 3\n"
     ]
    }
   ],
   "source": [
    "# Print contents of template file\n",
    "lines = open('./data/mymodel.lpf.tpl', 'r').readlines()\n",
    "for l in lines:\n",
    "    print(l.strip())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-Dimensional Transient Arrays\n",
    "\n",
    "Flopy supports parameterization of transient two dimensional arrays, like recharge.  This is similar to the approach for three dimensional static arrays, but there are some important differences in how span is specified.  The parameter span here is also a dictionary, and it must contain a 'kper' key, which corresponds to a list of stress periods (zero based, of course) for which the parameter applies.  The span dictionary must also contain an 'idx' key.  If span['idx'] is None, then the parameter is a multiplier for those stress periods.  If span['idx'] is a tuple (iarray, jarray), where iarray and jarray are a list of array indices, or a boolean array of shape (nrow, ncol), then the parameter applies only to the cells specified in idx."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define the model dimensions (made smaller for easier viewing)\n",
    "nlay = 3\n",
    "nrow = 5\n",
    "ncol = 5\n",
    "nper = 3\n",
    "\n",
    "# Create the flopy model object and add the dis and lpf packages\n",
    "m = flopy.modflow.Modflow(modelname='mymodel', model_ws='./data')\n",
    "dis = flopy.modflow.ModflowDis(m, nlay, nrow, ncol, nper=nper)\n",
    "lpf = flopy.modflow.ModflowLpf(m, hk=10.)\n",
    "rch = flopy.modflow.ModflowRch(m, rech={0: 0.001, 2: 0.003})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plist = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a multiplier parameter for recharge\n",
    "mfpackage = 'rch'\n",
    "partype = 'rech'\n",
    "parname = 'RECH_MULT'\n",
    "startvalue = None\n",
    "lbound = None\n",
    "ubound = None\n",
    "transform = None\n",
    "\n",
    "# For a recharge multiplier, span['idx'] must be None\n",
    "idx = None\n",
    "span = {'kpers': [0, 1, 2], 'idx': idx}\n",
    "p = params.Params(mfpackage, partype, parname, startvalue,\n",
    "                  lbound, ubound, span)\n",
    "plist.append(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The following packages will be parameterized: ['RCH']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Write the template file\n",
    "tw = tplwriter.TemplateWriter(m, plist)\n",
    "tw.write_template()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ptf ~\n",
      "# RCH for MODFLOW, generated by Flopy.\n",
      "3         0\n",
      "1        -1 # Stress period 1\n",
      "INTERNAL ~   RECH_MULT   ~ (FREE) -1      #rech_1\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "1        -1 # Stress period 2\n",
      "INTERNAL ~   RECH_MULT   ~ (FREE) -1      #rech_1\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "1        -1 # Stress period 3\n",
      "INTERNAL ~   RECH_MULT   ~ (FREE) -1      #rech_3\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n"
     ]
    }
   ],
   "source": [
    "# Print the results\n",
    "lines = open('./data/mymodel.rch.tpl', 'r').readlines()\n",
    "for l in lines:\n",
    "    print(l.strip())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplier parameters can also be combined with index parameters as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plist = []\n",
    "\n",
    "# Create a multiplier parameter for recharge\n",
    "mfpackage = 'rch'\n",
    "partype = 'rech'\n",
    "parname = 'RECH_MULT'\n",
    "startvalue = None\n",
    "lbound = None\n",
    "ubound = None\n",
    "transform = None\n",
    "\n",
    "# For a recharge multiplier, span['idx'] must be None\n",
    "span = {'kpers': [1, 2], 'idx': None}\n",
    "p = params.Params(mfpackage, partype, parname, startvalue,\n",
    "                  lbound, ubound, span)\n",
    "plist.append(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now create an index parameter\n",
    "mfpackage = 'rch'\n",
    "partype = 'rech'\n",
    "parname = 'RECH_ZONE'\n",
    "startvalue = None\n",
    "lbound = None\n",
    "ubound = None\n",
    "transform = None\n",
    "\n",
    "# For a recharge index parameter, span['idx'] must be a boolean array or tuple of array indices\n",
    "idx = np.empty((nrow, ncol), dtype=np.bool)\n",
    "idx[0:3, 0:3] = True\n",
    "span = {'kpers': [1], 'idx': idx}\n",
    "p = params.Params(mfpackage, partype, parname, startvalue,\n",
    "                  lbound, ubound, span)\n",
    "plist.append(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The following packages will be parameterized: ['RCH']\n",
      "\n",
      "ptf ~\n",
      "# RCH for MODFLOW, generated by Flopy.\n",
      "3         0\n",
      "1        -1 # Stress period 1\n",
      "CONSTANT    1.000000E-03                           #rech_1\n",
      "1        -1 # Stress period 2\n",
      "INTERNAL ~   RECH_MULT   ~ (FREE) -1      #rech_1\n",
      "~  RECH_ZONE  ~ ~  RECH_ZONE  ~ ~  RECH_ZONE  ~ 0.0010000000474974513 0.0010000000474974513\n",
      "~  RECH_ZONE  ~ ~  RECH_ZONE  ~ ~  RECH_ZONE  ~ 0.0010000000474974513 0.0010000000474974513\n",
      "~  RECH_ZONE  ~ ~  RECH_ZONE  ~ ~  RECH_ZONE  ~ 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513 0.0010000000474974513\n",
      "1        -1 # Stress period 3\n",
      "INTERNAL ~   RECH_MULT   ~ (FREE) -1      #rech_3\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n",
      "0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032 0.003000000026077032\n"
     ]
    }
   ],
   "source": [
    "# Write the template file\n",
    "tw = tplwriter.TemplateWriter(m, plist)\n",
    "tw.write_template()\n",
    "\n",
    "# Print the results\n",
    "lines = open('./data/mymodel.rch.tpl', 'r').readlines()\n",
    "for l in lines:\n",
    "    print(l.strip())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}