Getting started with Flopy

Using MODFLOW 6 to simulate a Darcy Tube.

MODFLOW 6 and Flopy

  • MODFLOW
    • Software developed by the USGS, first USGS documentation in 1983
    • Written in FORTRAN, has evolved adding new processes using new FORTRAN features
    • Current version (in 2021) is MODFLOW-6, a significant change from previous versions
    • Command line program
      • Reads text files defining model, Calculates hydraulic head & water fluxes
      • Modular structure defining different features
      • Multiple commercial interfaces, USGS developed ModelMuse and Flopy
  • Flopy
    • Python package, leverages python capabilities for MODFLOW input/output
    • Allows use of plotting package, data manipulation
    • can use to run multiple models (via python packages)

Available at url{https://www.usgs.gov/mission-areas/water-resources/science/modflow-and-related-programs}

A First Flopy Model

To provide a Flopy overview, the simplest groundwater flow model I can conceive of, a Darcy tube, will be created.

  • one-dimension
  • homogeneous materials with a K = 1⋅10 − 5(m)/(sec)
  • steady state
  • constant head boundaries on ends

In the text below, I present a brief descriptions of the packages used followed by the Python code needed to implement them.

Flopy packages

The packages imported from Flopy to create this model include:

  • import flopy.mf6.modflow : to access flopy MODFLOW-6 packages, will set this to mf
  • mf.MFSimulation : creates simulations instance
    • mf. mf.ModflowTdis : handles simulation time
    • mf.ModflowIms : solves matrix equations
    • mf.ModflowGwf : creates instance of groundwater processes
      • mf.ModflowGwfdis : spatial parameters and grid (are other options)
      • mf.ModflowGwfic : set initial conditions in model (e.g. hydr. head)
      • mf.ModflowGwfnpf : set cell hydraulic properties
      • mf.ModflowGwfchd : I.D. const. head cells
      • mf.mfgwfoc.ModflowGwfoc : control what is output and its formatting

Import Modules, Create Base Model Object

The python text used to begin this modeling example follow:

  • Import appropriate modules to use flopy with MODFLOW6 and post process output:
1
2
3
import flopy.mf6.modflow as mf
import flopy.utils as utils
import numpy as np
  • Instantiate the base model, a python object, to be defined with other methods:
1
2
3
4
5
6
7
8
9
work_dir = "./"
mod_name = "DarcyTube"

model = mf.MFSimulation(
    sim_name=mod_name,  # stem of file names that will be created
    sim_ws=work_dir,  # directory were modflow input and output files stored
    verbosity_level=1,
    exe_name="mf6", # name of modflow executable, may need to include path
)

Model Time, Solver, and Basic Groundwater Processes

The number of stress periods and time steps in MODFLOW are set with the 'Tdis' package; the equations that characterize the groundwater flow problem are solved with the 'Ims' package, and the 'Gwf' object that defines the different flow processes (incudes grid, hydraulic parameters, and source/sink in model).

  • MFmodflowTdis : Time Discretization, stress periods and time steps
  • ModflowIms : Iterative Model Solver, define matrix solver
  • ModflowGwf : Groundwater flow, object to store flow processes
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
time_dis = mf.ModflowTdis(
    model, pname="TimeDiscret",
    time_units="seconds", nper=1,
    perioddata=[(1.0, 1, 1.0)] # period length, num steps, mutliplier
)
# steady state model in one stress period, these times will not be used

solver = mf.ModflowIms(model, pname="solver", complexity="SIMPLE")

gwtr_flow = mf.ModflowGwf(model, modelname=mod_name, model_nam_file=f"{mod_name}.nam")

Spatial Grid, Initial Conditions, Node Parameters

Several options are available in MODFLOW6 to define the finite volume grid. In this model the traditional, and simplest, grid structure (a rectilinear grid) is defined with the 'dis' package. Initial conditions or starting hydraulic head values are defined with the 'ic' MODFLOW package. Hydraulic conductance (K), dewatering attributes, and cell types (unconfined, confined, perched etc.) are defined in the 'nps' package. Below, ' icelltype=1' indicates the saturated thickness is calculated from the heads when cell water level is below the top of a cell.

  • ModflowGwfdis : Discretise grid, grid size and cell sizes for rectilinear (tradiational) grid, other grids are available in MODFLOW 6
  • ModflowGwfic : Initial Conditions, set the initial hydr. heads
  • ModflowGwfnpf : Node-Property Flow, set cell properties (hydr params, rewetting, dewatering)
1
2
3
4
5
6
7
8
9
grid_dis = mf.ModflowGwfdis(
    gwtr_flow, nlay=1, nrow=1, ncol=20, delr=0.1, delc=0.1, top=0.1, botm=0.0
)
# 3-D numpy grid used to define 1-D initial hydr. heads
init_cond = mf.ModflowGwfic(
    gwtr_flow, pname="init_cond", strt=10.0 * np.ones((1, 1, 20))
)
# icelltype=1, 'non convertable' cell, am assuming won't de-water
node_prop = mf.ModflowGwfnpf(gwtr_flow, icelltype=1, k=1.0e-5, save_flows=True)

Set Boundary Conditions

In this model, only constant hydraulic head boundaries are set. Below, a loop is used to do this. For this simple model, it would have been easier to explicitly write out this list of parameters.

  • Only use constant head boundaries, assigned in first and last columns in row
  • No-flow assumed along edges (default behavior of MODFLOW).
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
const_hd_lst = [] # somewhat convoluted way of setting two const head cells
for col in range(grid_dis.ncol.get_data()):  # loop over ncol
    if (col == 0) or (col == 19):
                    const_hd_lst.append(((0, 0, col), 10.0 - col / 38.0))

const_hd = mf.ModflowGwfchd(
            gwtr_flow,
            maxbound=2,  # max number of const hd cells in model
            stress_period_data=const_hd_lst,
            save_flows=True,
)

Defining Model Output

MODFLOW output is controlled by parameters in the 'oc' package. This includes the frequency that hydraulic head values and water budget data are saved to the listing file (or separate head of budget files).

  • Set output file names to record head values and discharge rates
  • Set frequently (stress periods & time steps) output is saved
  • Set frequency of output displayed in listing file
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
head_file = f"{mod_name}.hds"
head_filerecord = [head_file]
budgetfile = f"{mod_name}.cbb"
budget_filerecord = [budgetfile]
saverecord = [("HEAD", "ALL"), ("BUDGET", "ALL")]
printrecord = [("HEAD", "LAST")]
oc = mf.mfgwfoc.ModflowGwfoc(
    gwtr_flow, pname="oc",
    saverecord=saverecord,
    head_filerecord=head_filerecord,
    budget_filerecord=budget_filerecord,
    printrecord=printrecord,
)

Create MODFLOW Input Files and Run Model

Once the desired packages from MODFLOW6 are defined in the Python script, the input files are written to the working directory (defined in the 'Simulation' package). The model object is then run, and the modflow executable is called to process the input files written to the computer storage. Both MODFLOW and flopy need to be installed and the proper MODFLOW executable file name must be defined in 'Simulation' package for this to successfully run.

1
2
model.write_simulation()
success, buff = model.run_simulation()
  • Should see a set of MODFLOW input files after calling write_simulation
  • run_simulation runs model and produces output files including listing (\*.lst) file
1
2
3
4
5
6
# open MODFLOW binary files
heads=utils.binaryfile.HeadFile('./DarcyTube.hds')
flows=utils.binaryfile.CellBudgetFile('DarcyTube.cbb')
heads.list_records() # return records in 'heads' file
heads.get_data(kstpkper=(0,0)) # get head data for timestep=0, stress period=0
flows.get_data(kstpkper=(0,0)) # flopy will change indices to match FORTRAN conventions in output

Looking at the Data

Below is a screenshot after running this script. The working directory should have about 25 files in it with names like 'DarcyTube2.ic' and 'DarcyTube2.lst'. The file extension indicates the package associated with the file (e.g., '*.ic' for input control files). The text should indicate 'Normal termination of simulation' (i.e. no obvious errors are in the script). This does not mean the model is fine. There could still be large water balance errors, incorrect data entry, or a poorly conceptualized groundwater flow model. Users need to check the model to verify it is producing reasonable results and has properly converged on a solution.

Image of ipython session

Groundwater flow model grid for Darcy tube simulation
  • Image of ipython session running the flopy script simulating a Darcy tube
    • shows what files are written to computer storage
    • info about MODFLOW version
    • info about stress periods/time steps completed
    • if model terminated normally
  • in line 4, output head data using Flopy utilities:
    • 0.5 m total head drop, K= 1⋅10 − 5m⋅sec − 1, A=0.01 m2, over 2 m length
    • applying Darcy's Law, Q = 2.5⋅10 − 8m⋅sec − 1
    • applying to cells in col 1 and 19 (1.7 m apart), Q = 2.63⋅10 − 8m⋅sec − 1
  • Why the difference?
    • The constant head boundary conditions are averaged over cell
    • Should be treated at occurring at cell center, not at outer edge of model
    • If length is set to 1.9 m over full model length, get 'correct' answer. This distance is from the center of the first to the center of the last cell, instead of from edge to edge of the model.

Looking at the Data

Here the flux or discharge data from the compact budget are shown. These data can be accessed by opening the appropriate file (with a 'cbb' extension in this example) using the Flopy utility functions.

Image of ipython session

Screenshot of MODFLOW compact budget file data.
  • Another screenshot of ipython session, showing water budget file
    • Fluxes shown for all faces (only y-z or col-lay faces in this model)
    • Fluxes shown for:
      • source-sinks (none in this model)
      • for water from/to storage (none in this model)
      • for fluxes across faces (only faces perpendicular to x (or row) axis in this model
      • fluxes to/from constant head boundary conditions
    • Fluxes here match what was calculated from heads (with const. heads assigned at cell centers)

Looking at the Data

If errors are encountered when running MODFLOW, good starting points to correct errors are to 1) look at the listing file ('*.lst' fine in the work directory) and look at the input files (are they all present? do they end prematurely?).

MODFLOW listing file

Screenshot of Modflow listing file showing water budget.
  • MODFLOW listing files
    • file with 'lst' extension (the_model_name.lst), plain text file
    • Track input as MODFLOW reads data file
    • Will provide clues where errors are if MODFLOW run is unsuccessful, can sort out what files are causing problems
    • Water Budget (shown in image) should have water-balance errors very close to zero (<1%)
    • These files can get very large, may be slow to open in text editors. Can use more and less shell commands to scroll through and search larger files.

I hope this is helpful for people to get started with MODFLOW and Flopy. This is is a great tool for groundwater professionals and Flopy provides a great interface for making hard things easier to do. I believe this also forces the users to think more about what is being done when compared to a GUI interface.