Example for FloPy$_3$ methods note

Import the {\tt modflow} and {\tt utilities} subpackages of FloPy and give them the aliases {\tt fpm} and {\tt fpu}, respectively


In [5]:
import numpy as np
import flopy.modflow as fpm
import flopy.utils as fpu

Create a MODFLOW model object. Here, the MODFLOW model object is stored in a Python variable called {\tt model}, but this can be an arbitrary name. This object name is important as it will be used as a reference to the model in the remainder of the FloPy script. In addition, a {\tt modelname} is specified when the MODFLOW model object is created. This {\tt modelname} is used for all the files that are created by FloPy for this model.


In [15]:
exe = "mf2005"
model = fpm.Modflow(modelname='gwexample',exe_name=exe)

The discretization of the model is specified with the discretization file (DIS) of MODFLOW. The aquifer is divided into 201 cells of length 10 m and width 1 m. The first input of the discretization package is the name of the model object. All other input arguments are self explanatory.


In [16]:
fpm.ModflowDis(model, nlay=1, nrow=1, ncol=201,
               delr=10, delc=1, top=50, botm=0)


Out[16]:
    MODFLOW Discretization Package Class.

    Parameters
    ----------
    model : model object
        The model object (of type :class:`flopy.modflow.Modflow`) to which
        this package will be added.
    nlay : int
        Number of model layers (the default is 1).
    nrow : int
        Number of model rows (the default is 2).
    ncol : int
        Number of model columns (the default is 2).
    nper : int
        Number of model stress periods (the default is 1).
    delr : float or array of floats (ncol), optional
        An array of spacings along a row (the default is 1.0).
    delc : float or array of floats (nrow), optional
        An array of spacings along a column (the default is 0.0).
    laycbd : int or array of ints (nlay), optional
        An array of flags indicating whether or not a layer has a Quasi-3D
        confining bed below it. 0 indicates no confining bed, and not zero
        indicates a confining bed. LAYCBD for the bottom layer must be 0. (the
        default is 0)
    top : float or array of floats (nrow, ncol), optional
        An array of the top elevation of layer 1. For the common situation in
        which the top layer represents a water-table aquifer, it may be
        reasonable to set Top equal to land-surface elevation (the default is
        1.0)
    botm : float or array of floats (nlay, nrow, ncol), optional
        An array of the bottom elevation for each model cell (the default is
        0.)
    perlen : float or array of floats (nper)
        An array of the stress period lengths.
    nstp : int or array of ints (nper)
        Number of time steps in each stress period (default is 1).
    tsmult : float or array of floats (nper)
        Time step multiplier (default is 1.0).
    steady : boolean or array of boolean (nper)
        true or False indicating whether or not stress period is steady state
        (default is True).
    itmuni : int
        Time units, default is days (4)
    lenuni : int
        Length units, default is meters (2)
    extension : string
        Filename extension (default is 'dis')
    unitnumber : int
        File unit number (default is 11).
    xul : float
        x coordinate of upper left corner of the grid, default is None
    yul : float
        y coordinate of upper left corner of the grid, default is None
    rotation : float
        clockwise rotation (in degrees) of the grid about the upper left
        corner. default is 0.0
    proj4_str : str
        PROJ4 string that defines the xul-yul coordinate system
        (.e.g. '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ').
        Can be an EPSG code (e.g. 'EPSG:4326'). Default is 'EPSG:4326'
    start_dateteim : str
        starting datetime of the simulation. default is '1/1/1970'

    Attributes
    ----------
    heading : str
        Text string written to top of package input file.

    Methods
    -------

    See Also
    --------

    Notes
    -----

    Examples
    --------

    >>> import flopy
    >>> m = flopy.modflow.Modflow()
    >>> dis = flopy.modflow.ModflowDis(m)

     _ModflowDis__thickness = <flopy.utils.util_array.Util3d object at 0x108747dd0> ('flopy.utils.util_array.Util3d)
 acceptable_dtypes (list, items = 3
 allowDuplicates = False (bool)
 botm = <flopy.utils.util_array.Util3d object at 0x108747a90> ('flopy.utils.util_array.Util3d)
 delc = <flopy.utils.util_array.Util2d object at 0x108747990> ('flopy.utils.util_array.Util2d)
 delr = <flopy.utils.util_array.Util2d object at 0x108747910> ('flopy.utils.util_array.Util2d)
 extra = 
 file_name = gwexample.dis
 fn_path = ./gwexample.dis (str)
 itmuni = 4 (int)
 itmuni_dict = {0: 'undefined', 1: 'seconds', 2: 'minutes', 3: 'hours', 4: 'days', 5: 'years'} (dict)
 laycbd = <flopy.utils.util_array.Util2d object at 0x1087478d0> ('flopy.utils.util_array.Util2d)
 lenuni = 2 (int)
 ncol = 201 (int)
 nlay = 1 (int)
 nper = 1 (int)
 nrow = 1 (int)
 nstp = <flopy.utils.util_array.Util2d object at 0x108747bd0> ('flopy.utils.util_array.Util2d)
 perlen = <flopy.utils.util_array.Util2d object at 0x108747ad0> ('flopy.utils.util_array.Util2d)
 sr = xul:0, yul:1, rotation:0, proj4_str:EPSG:4326 ('flopy.utils.reference.SpatialReference)
 start_datetime = 1/1/1970 (str)
 steady = <flopy.utils.util_array.Util2d object at 0x108747cd0> ('flopy.utils.util_array.Util2d)
 top = <flopy.utils.util_array.Util2d object at 0x108747a10> ('flopy.utils.util_array.Util2d)
 tsmult = <flopy.utils.util_array.Util2d object at 0x108747c50> ('flopy.utils.util_array.Util2d)
 unit_number = 11

Active cells and the like are defined with the Basic package (BAS), which is required for every MODFLOW model. It contains the {\tt ibound} array, which is used to specify which cells are active (value is positive), inactive (value is 0), or fixed head (value is negative). The {\tt numpy} package (aliased as {\tt np}) can be used to quickly initialize the {\tt ibound} array with values of 1, and then set the {\tt ibound} value for the first and last columns to -1. The {\tt numpy} package (and Python, in general) uses zero-based indexing and supports negative indexing so that row 1 and column 1, and row 1 and column 201, can be referenced as [0, 0], and [0, -1], respectively. Although this simulation is for steady flow, starting heads still need to be specified. They are used as the head for fixed-head cells (where {\tt ibound} is negative), and as a starting point to compute the saturated thickness for cases of unconfined flow.


In [17]:
ibound = np.ones((1, 201))
ibound[0, 0] = ibound[0, -1] = -1
fpm.ModflowBas(model, ibound=ibound, strt=20)


Out[17]:
    MODFLOW Basic Package Class.

    Parameters
    ----------
    model : model object
        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which
        this package will be added.
    ibound : array of ints, optional
        The ibound array (the default is 1).
    strt : array of floats, optional
        An array of starting heads (the default is 1.0).
    ixsec : bool, optional
        Indication of whether model is cross sectional or not (the default is
        False).
    ichflg : bool, optional
        Flag indicating that flows between constant head cells should be
        calculated (the default is False).
    stoper : float
        percent discrepancy that is compared to the budget percent discrepancy
        continue when the solver convergence criteria are not met.  Execution
        will unless the budget percent discrepancy is greater than stoper
        (default is None). MODFLOW-2005 only
    hnoflo : float
        Head value assigned to inactive cells (default is -999.99).
    extension : str, optional
        File extension (default is 'bas'.
    unitnumber : int, optional
        FORTRAN unit number for this package (default is 13).

    Attributes
    ----------
    heading : str
        Text string written to top of package input file.
    options : list of str
        Can be either or a combination of XSECTION, CHTOCH or FREE.
    ifrefm : bool
        Indicates whether or not packages will be written as free format.

    Methods
    -------

    See Also
    --------

    Notes
    -----

    Examples
    --------

    >>> import flopy
    >>> m = flopy.modflow.Modflow()
    >>> bas = flopy.modflow.ModflowBas(m)

     acceptable_dtypes (list, items = 3
 allowDuplicates = False (bool)
 extra = 
 file_name = gwexample.bas
 fn_path = ./gwexample.bas (str)
 hnoflo = -999.99 (float)
 ibound = <flopy.utils.util_array.Util3d object at 0x10874c810> ('flopy.utils.util_array.Util3d)
 ichflg = False (bool)
 ifrefm = True (bool)
 ixsec = False (bool)
 options =  (str)
 stoper = None (NoneType)
 strt = <flopy.utils.util_array.Util3d object at 0x10874c950> ('flopy.utils.util_array.Util3d)
 unit_number = 13

The hydraulic properties of the aquifer are specified with the Layer Properties Flow (LPF) package (alternatively, the Block Centered Flow (BCF) package may be used). Only the hydraulic conductivity of the aquifer and the layer type ({\tt laytyp}) need to be specified. The latter is set to 1, which means that MODFLOW will calculate the saturated thickness differently depending on whether or not the head is above the top of the aquifer.


In [18]:
fpm.ModflowLpf(model, hk=10, laytyp=1)


Out[18]:
    MODFLOW Layer Property Flow Package Class.

    Parameters
    ----------
    model : model object
        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which
        this package will be added.
    ipakcb : int
        A flag that is used to determine if cell-by-cell budget data should be
        saved. If ipakcb is non-zero cell-by-cell budget data will be saved.
        (default is 53)
    hdry : float
        Is the head that is assigned to cells that are converted to dry during
        a simulation. Although this value plays no role in the model
        calculations, it is useful as an indicator when looking at the
        resulting heads that are output from the model. HDRY is thus similar
        to HNOFLO in the Basic Package, which is the value assigned to cells
        that are no-flow cells at the start of a model simulation. (default
        is -1.e30).
    laytyp : int or array of ints (nlay)
        Layer type (default is 0).
    layavg : int or array of ints (nlay)
        Layer average (default is 0).
        0 is harmonic mean
        1 is logarithmic mean
        2 is arithmetic mean of saturated thickness and logarithmic mean of
        of hydraulic conductivity
    chani : float or array of floats (nlay)
        contains a value for each layer that is a flag or the horizontal
        anisotropy. If CHANI is less than or equal to 0, then variable HANI
        defines horizontal anisotropy. If CHANI is greater than 0, then CHANI
        is the horizontal anisotropy for the entire layer, and HANI is not
        read. If any HANI parameters are used, CHANI for all layers must be
        less than or equal to 0. Use as many records as needed to enter a
        value of CHANI for each layer. The horizontal anisotropy is the ratio
        of the hydraulic conductivity along columns (the Y direction) to the
        hydraulic conductivity along rows (the X direction).
    layvka : float or array of floats (nlay)
        a flag for each layer that indicates whether variable VKA is vertical
        hydraulic conductivity or the ratio of horizontal to vertical
        hydraulic conductivity.
    laywet : float or array of floats (nlay)
        contains a flag for each layer that indicates if wetting is active.
    wetfct : float
        is a factor that is included in the calculation of the head that is
        initially established at a cell when it is converted from dry to wet.
        (default is 0.1).
    iwetit : int
        is the iteration interval for attempting to wet cells. Wetting is
        attempted every IWETIT iteration. If using the PCG solver
        (Hill, 1990), this applies to outer iterations, not inner iterations.
        If IWETIT  less than or equal to 0, it is changed to 1.
        (default is 1).
    ihdwet : int
        is a flag that determines which equation is used to define the
        initial head at cells that become wet. (default is 0)
    hk : float or array of floats (nlay, nrow, ncol)
        is the hydraulic conductivity along rows. HK is multiplied by
        horizontal anisotropy (see CHANI and HANI) to obtain hydraulic
        conductivity along columns. (default is 1.0).
    hani : float or array of floats (nlay, nrow, ncol)
        is the ratio of hydraulic conductivity along columns to hydraulic
        conductivity along rows, where HK of item 10 specifies the hydraulic
        conductivity along rows. Thus, the hydraulic conductivity along
        columns is the product of the values in HK and HANI.
        (default is 1.0).
    vka : float or array of floats (nlay, nrow, ncol)
        is either vertical hydraulic conductivity or the ratio of horizontal
        to vertical hydraulic conductivity depending on the value of LAYVKA.
        (default is 1.0).
    ss : float or array of floats (nlay, nrow, ncol)
        is specific storage unless the STORAGECOEFFICIENT option is used.
        When STORAGECOEFFICIENT is used, Ss is confined storage coefficient.
        (default is 1.e-5).
    sy : float or array of floats (nlay, nrow, ncol)
        is specific yield. (default is 0.15).
    vkcb : float or array of floats (nlay, nrow, ncol)
        is the vertical hydraulic conductivity of a Quasi-three-dimensional
        confining bed below a layer. (default is 0.0).
    wetdry : float or array of floats (nlay, nrow, ncol)
        is a combination of the wetting threshold and a flag to indicate
        which neighboring cells can cause a cell to become wet.
        (default is -0.01).
    storagecoefficient : boolean
        indicates that variable Ss and SS parameters are read as storage
        coefficient rather than specific storage. (default is False).
    constantcv : boolean
         indicates that vertical conductance for an unconfined cell is
         computed from the cell thickness rather than the saturated thickness.
         The CONSTANTCV option automatically invokes the NOCVCORRECTION
         option. (default is False).
    thickstrt : boolean
        indicates that layers having a negative LAYTYP are confined, and their
        cell thickness for conductance calculations will be computed as
        STRT-BOT rather than TOP-BOT. (default is False).
    nocvcorrection : boolean
        indicates that vertical conductance is not corrected when the vertical
        flow correction is applied. (default is False).
    novfc : boolean
         turns off the vertical flow correction under dewatered conditions.
         This option turns off the vertical flow calculation described on p.
         5-8 of USGS Techniques and Methods Report 6-A16 and the vertical
         conductance correction described on p. 5-18 of that report.
         (default is False).
    extension : string
        Filename extension (default is 'lpf')
    unitnumber : int
        File unit number (default is 15).


    Attributes
    ----------

    Methods
    -------

    See Also
    --------

    Notes
    -----

    Examples
    --------

    >>> import flopy
    >>> m = flopy.modflow.Modflow()
    >>> lpf = flopy.modflow.ModflowLpf(m)

     acceptable_dtypes (list, items = 3
 allowDuplicates = False (bool)
 chani = <flopy.utils.util_array.Util2d object at 0x108747310> ('flopy.utils.util_array.Util2d)
 extra = 
 file_name = gwexample.lpf
 fn_path = ./gwexample.lpf (str)
 hani = <flopy.utils.util_array.Util3d object at 0x10874c6d0> ('flopy.utils.util_array.Util3d)
 hdry = -1e+30 (float)
 hk = <flopy.utils.util_array.Util3d object at 0x10874c650> ('flopy.utils.util_array.Util3d)
 ihdwet = 0 (int)
 ipakcb = 53 (int)
 iwetit = 1 (int)
 layavg = <flopy.utils.util_array.Util2d object at 0x108747510> ('flopy.utils.util_array.Util2d)
 laytyp = <flopy.utils.util_array.Util2d object at 0x108747750> ('flopy.utils.util_array.Util2d)
 layvka = <flopy.utils.util_array.Util2d object at 0x108747fd0> ('flopy.utils.util_array.Util2d)
 laywet = <flopy.utils.util_array.Util2d object at 0x10874c7d0> ('flopy.utils.util_array.Util2d)
 nplpf = 0 (int)
 options =   (str)
 ss = <flopy.utils.util_array.Util3d object at 0x10874c510> ('flopy.utils.util_array.Util3d)
 sy = <flopy.utils.util_array.Util3d object at 0x1086e20d0> ('flopy.utils.util_array.Util3d)
 unit_number = 15
 vka = <flopy.utils.util_array.Util3d object at 0x108747810> ('flopy.utils.util_array.Util3d)
 vkcb = <flopy.utils.util_array.Util3d object at 0x1086e2190> ('flopy.utils.util_array.Util3d)
 wetdry = <flopy.utils.util_array.Util3d object at 0x1086e2210> ('flopy.utils.util_array.Util3d)
 wetfct = 0.1 (float)

Aquifer recharge is simulated with the Recharge package (RCH) and the extraction of water at the two ditches is simulated with the Well package (WEL). The latter requires specification of the layer, row, column, and injection rate of the well for each stress period. The layers, rows, columns, and the stress period are numbered (consistent with Python's zero-based numbering convention) starting at 0. The required data are stored in a Python dictionary ({\tt lrcQ} in the code below), which is used in FloPy to store data that can vary by stress period. The {\tt lrcQ} dictionary specifies that two wells (one in cell 1, 1, 51 and one in cell 1, 1, 151), each with a rate of -1 m$^3$/m/d, will be active for the first stress period. Because this is a steady-state model, there is only one stress period and therefore only one entry in the dictionary.


In [19]:
fpm.ModflowRch(model, rech=0.001)
lrcQ = {0: [[0, 0, 50, -1], [0, 0, 150, -1]]}
fpm.ModflowWel(model, stress_period_data=lrcQ)


Out[19]:
    MODFLOW Well Package Class.

    Parameters
    ----------
    model : model object
        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which
        this package will be added.
    ipakcb : int
        A flag that is used to determine if cell-by-cell budget data should be
        saved. If ipakcb is non-zero cell-by-cell budget data will be saved.
        (default is 0).
    stress_period_data : list of boundaries, or recarray of boundaries, or
        dictionary of boundaries
        Each well is defined through definition of
        layer (int), row (int), column (int), flux (float).
        The simplest form is a dictionary with a lists of boundaries for each
        stress period, where each list of boundaries itself is a list of
        boundaries. Indices of the dictionary are the numbers of the stress
        period. This gives the form of::

            stress_period_data =
            {0: [
                [lay, row, col, flux],
                [lay, row, col, flux],
                [lay, row, col, flux]
                ],
            1:  [
                [lay, row, col, flux],
                [lay, row, col, flux],
                [lay, row, col, flux]
                ], ...
            kper:
                [
                [lay, row, col, flux],
                [lay, row, col, flux],
                [lay, row, col, flux]
                ]
            }

        Note that if the number of lists is smaller than the number of stress
        periods, then the last list of wells will apply until the end of the
        simulation. Full details of all options to specify stress_period_data
        can be found in the flopy3 boundaries Notebook in the basic
        subdirectory of the examples directory
    options : list of strings
        Package options. (default is None).
    extension : string
        Filename extension (default is 'wel')
    unitnumber : int
        File unit number (default is 11).

    Attributes
    ----------
    mxactw : int
        Maximum number of wells for a stress period.  This is calculated
        automatically by FloPy based on the information in
        stress_period_data.

    Methods
    -------

    See Also
    --------

    Notes
    -----
    Parameters are not supported in FloPy.

    Examples
    --------

    >>> import flopy
    >>> m = flopy.modflow.Modflow()
    >>> lrcq = {0:[[2, 3, 4, -100.]], 1:[[2, 3, 4, -100.]]}
    >>> wel = flopy.modflow.ModflowWel(m, stress_period_data=lrcq)

     acceptable_dtypes (list, items = 3
 allowDuplicates = False (bool)
 dtype = [('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('flux', '<f4')] (numpy.dtype)
 extra = 
 file_name = gwexample.wel
 fn_path = ./gwexample.wel (str)
 ipakcb = 0 (int)
 np = 0 (int)
 options (list, items = 0
 specify = False (bool)
 stress_period_data = <flopy.utils.util_list.MfList object at 0x1086e2a10> ('flopy.utils.util_list.MfList)
 unit_number = 20

The Preconditioned Conjugate-Gradient (PCG) solver, using the default settings, is specified to solve the model.


In [20]:
fpm.ModflowPcg(model)


Out[20]:
    MODFLOW Pcg Package Class.

    Parameters
    ----------
    model : model object
        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which
        this package will be added.
    mxiter : int
        maximum number of outer iterations. (default is 50)
    iter1 : int
        maximum number of inner iterations. (default is 30)
    npcond : int
        flag used to select the matrix conditioning method. (default is 1).
        specify npcond = 1 for Modified Incomplete Cholesky.
        specify npcond = 2 for Polynomial.
    hclose : float
        is the head change criterion for convergence. (default is 1e-5).
    rclose : float
        is the residual criterion for convergence. (default is 1e-5)
    relax : float
        is the relaxation parameter used with npcond = 1. (default is 1.0)
    nbpol : int
        is only used when npcond = 2 to indicate whether the estimate of the
        upper bound on the maximum eigenvalue is 2.0, or whether the estimate
        will be calculated. nbpol = 2 is used to specify the value is 2.0;
        for any other value of nbpol, the estimate is calculated. Convergence
        is generally insensitive to this parameter. (default is 2).
    iprpcg : int
        solver print out interval. (default is 0).
    mutpcg : int
        If mutpcg = 0, tables of maximum head change and residual will be
        printed each iteration.
        If mutpcg = 1, only the total number of iterations will be printed.
        If mutpcg = 2, no information will be printed.
        If mutpcg = 3, information will only be printed if convergence fails.
        (default is 3).
    damp : float
        is the steady-state damping factor. (default is 1.)
    dampt : float
        is the transient damping factor. (default is 1.)
    ihcofadd : int
        is a flag that determines what happens to an active cell that is
        surrounded by dry cells.  (default is 0). If ihcofadd=0, cell
        converts to dry regardless of HCOF value. This is the default, which
        is the way PCG2 worked prior to the addition of this option. If
        ihcofadd<>0, cell converts to dry only if HCOF has no head-dependent
        stresses or storage terms.
    extension : list string
        Filename extension (default is 'pcg')
    unitnumber : int
        File unit number (default is 27).

    Attributes
    ----------

    Methods
    -------

    See Also
    --------

    Notes
    -----

    Examples
    --------

    >>> import flopy
    >>> m = flopy.modflow.Modflow()
    >>> pcg = flopy.modflow.ModflowPcg(m)

     acceptable_dtypes (list, items = 3
 allowDuplicates = False (bool)
 damp = 1.0 (float)
 dampt = 1.0 (float)
 extra = 
 file_name = gwexample.pcg
 fn_path = ./gwexample.pcg (str)
 hclose = 1e-05 (float)
 ihcofadd = 0 (int)
 iprpcg = 0 (int)
 iter1 = 30 (int)
 mutpcg = 3 (int)
 mxiter = 50 (int)
 nbpol = 0 (int)
 npcond = 1 (int)
 rclose = 1e-05 (float)
 relax = 1.0 (float)
 unit_number = 27

The frequency and type of output that MODFLOW writes to an output file is specified with the Output Control (OC) package. In this case, the budget is printed and heads are saved (the default), so no arguments are needed.


In [21]:
fpm.ModflowOc(model)


Out[21]:
    MODFLOW Output Control Package Class.

    Parameters
    ----------
    model : model object
        The model object (of type :class:`flopy.modflow.mf.Modflow`) to which
        this package will be added.
    ihedfm : int
        is a code for the format in which heads will be printed.
        (default is 0).
    iddnfm : int
        is a code for the format in which heads will be printed.
        (default is 0).
    chedfm : string
        is a character value that specifies the format for saving heads.
        The format must contain 20 characters or less and must be a valid
        Fortran format that is enclosed in parentheses. The format must be
        enclosed in apostrophes if it contains one or more blanks or commas.
        The optional word LABEL after the format is used to indicate that
        each layer of output should be preceded with a line that defines the
        output (simulation time, the layer being output, and so forth). If
        there is no record specifying CHEDFM, then heads are written to a
        binary (unformatted) file. Binary files are usually more compact than
        text files, but they are not generally transportable among different
        computer operating systems or different Fortran compilers.
        (default is None)
    cddnfm : string
        is a character value that specifies the format for saving drawdown.
        The format must contain 20 characters or less and must be a valid
        Fortran format that is enclosed in parentheses. The format must be
        enclosed in apostrophes if it contains one or more blanks or commas.
        The optional word LABEL after the format is used to indicate that
        each layer of output should be preceded with a line that defines the
        output (simulation time, the layer being output, and so forth). If
        there is no record specifying CDDNFM, then drawdowns are written to a
        binary (unformatted) file. Binary files are usually more compact than
        text files, but they are not generally transportable among different
        computer operating systems or different Fortran compilers.
        (default is None)
    cboufm : string
        is a character value that specifies the format for saving ibound.
        The format must contain 20 characters or less and must be a valid
        Fortran format that is enclosed in parentheses. The format must be
        enclosed in apostrophes if it contains one or more blanks or commas.
        The optional word LABEL after the format is used to indicate that
        each layer of output should be preceded with a line that defines the
        output (simulation time, the layer being output, and so forth). If
        there is no record specifying CBOUFM, then ibounds are written to a
        binary (unformatted) file. Binary files are usually more compact than
        text files, but they are not generally transportable among different
        computer operating systems or different Fortran compilers.
        (default is None)
    stress_period_data : dictionary of of lists
        Dictionary key is a tuple with the zero-based period and step 
        (IPEROC, ITSOC) for each print/save option list. 
        (default is {(0,0):['save head']})
        
        The list can have any valid MODFLOW OC print/save option:
            PRINT HEAD
            PRINT DRAWDOWN
            PRINT BUDGET
            SAVE HEAD
            SAVE DRAWDOWN
            SAVE BUDGET
            SAVE IBOUND
            
            The lists can also include (1) DDREFERENCE in the list to reset 
            drawdown reference to the period and step and (2) a list of layers 
            for PRINT HEAD, SAVE HEAD, PRINT DRAWDOWN, SAVE DRAWDOWN, and
            SAVE IBOUND.
        
        The list is used for every stress period and time step after the 
        (IPEROC, ITSOC) tuple until a (IPEROC, ITSOC) tuple is entered with
        and empty list.
    compact : boolean
        Save results in compact budget form. (default is True).
    extension : list of strings
        (default is ['oc','hds','ddn','cbc']).
    unitnumber : list of ints
        (default is [14, 51, 52, 53]).

    Attributes
    ----------

    Methods
    -------

    See Also
    --------

    Notes
    -----
    The "words" method for specifying output control is the only option 
    available.  Also, the "compact" budget should normally be used as it 
    produces files that are typically much smaller.  The compact budget form is 
    also a requirement for using the MODPATH particle tracking program.

    Examples
    --------

    >>> import flopy
    >>> m = flopy.modflow.Modflow()
    >>> spd = {(0, 0): ['print head'],
    ...   (0, 1): [],
    ...   (0, 249): ['print head'],
    ...   (0, 250): [],
    ...   (0, 499): ['print head', 'save ibound'],
    ...   (0, 500): [],
    ...   (0, 749): ['print head', 'ddreference'],
    ...   (0, 750): [],
    ...   (0, 999): ['print head']}
    >>> oc = flopy.modflow.ModflowOc3(m, stress_period_data=spd, cboufm='(20i5)')

     acceptable_dtypes (list, items = 3
 allowDuplicates = False (bool)
 cboufm = None (NoneType)
 cddnfm = None (NoneType)
 chedfm = None (NoneType)
 compact = True (bool)
 extra (list, items = 4
 file_name (list, items = 4
 fn_path = ./gwexample.oc (str)
 ibouun = 0 (int)
 iddnfm = 0 (int)
 ihedfm = 0 (int)
 stress_period_data = {(0, 0): ['save head']} (dict)
 unit_number (list, items = 4

Finally the MODFLOW input files are written (eight files for this model) and the model is run. This requires, of course, that MODFLOW is installed on your computer and FloPy can find the executable in your path.


In [22]:
model.write_input()
model.run_model()


FloPy is using the following executable to run the model: /Users/jwhite/bin/mf2005

                                  MODFLOW-2005     
    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL
                             Version 1.11.00 8/8/2013                        

 Using NAME file: gwexample.nam 
 Run start date and time (yyyy/mm/dd hh:mm:ss): 2016/01/19  9:13:09

 Solving:  Stress period:     1    Time step:     1    Ground-Water Flow Eqn.
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2016/01/19  9:13:09
 Elapsed run time:  0.007 Seconds

  Normal termination of simulation
Out[22]:
(True, [])

After MODFLOW has responded with the positive {\tt Normal termination of simulation}, the calculated heads can be read from the binary output file. First a file object is created. As the modelname used for all MODFLOW files was specified as {\tt gwexample} in step 1, the file with the heads is called {\tt gwexample.hds}. FloPy includes functions to read data from the file object, including heads for specified layers or time steps, or head time series at individual cells. For this simple mode, all computed heads are read.


In [24]:
hfile = fpu.HeadFile('gwexample.hds')
h = hfile.get_data(totim=1.0)

The heads are now stored in the Python variable {\tt h}. FloPy includes powerful plotting functions to plot the grid, boundary conditions, head, etc. This functionality is demonstrated later. For this simple one-dimensional example, a plot is created with the matplotlib package


In [ ]:
import matplotlib.pyplot as plt
ax = plt.subplot(111)
x = model.dis.sr.xcentergrid[0]
ax.plot(x,h[0,0,:])
ax.set_xlim(0,x.max())
ax.set_xlabel("x(m)")
ax.set_ylabel("head(m)")
plt.show()

In [ ]: