Flopy is the Python library to create, run and represent results from MODFLOW models, including support for MODFLOW 6. With Flopy, the model definition and set up of boundary conditions can be done with lines of Python code, and some Python tools of optimization and machine learning can interact with the MODFLOW models. This tutorial show a basic example of model creation, configuration and simulation in Python 3 on a Jupyter Notebook. The tutorial also shows the procedure to import the MODFLOW results done with Flopy in Model Muse.

If you want to install Flopy on the Anaconda 3 framework please have a look on this tutorial.

Video

Code

This is the code used to create, setup and run the MODFLOW model:

import flopy, os
import numpy as np

modelname = "flopymodel"
os.chdir('ModelFiles')
mf = flopy.modflow.Modflow(modelname, exe_name='../mf2005')

# Model domain and grid definition
Lx = 1000.
Ly = 1000.
ztop = 100.
zbot = -100.
nlay = 4
nrow = 10
ncol = 10
delr = Lx/ncol
delc = Ly/nrow
delv = (ztop - zbot) / nlay
botm = np.linspace(ztop, zbot, nlay + 1)

# Create the discretization object
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc,
                               top=ztop, botm=botm[1:])

# Variables for the BAS package
ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
ibound[0, :, 0] = -1
ibound[0, :, -1] = -1
strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
strt[0, :, 0] = 100.
strt[0, :, -1] = 80.
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)

# Add LPF package to the MODFLOW model
lpf = flopy.modflow.ModflowLpf(mf, hk=1e-7, vka=1e-7)

# Add OC package to the MODFLOW model
oc = flopy.modflow.ModflowOc(mf)

# Add PCG package to the MODFLOW model
pcg = flopy.modflow.ModflowPcg(mf)

# Write the MODFLOW model input files
mf.write_input()

# Run the MODFLOW model
success, buff = mf.run_model()

#Plot results
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
plt.subplot(1,2,1,aspect='equal')
hds = bf.HeadFile(modelname+'.hds')
head = hds.get_data(totim=1.0)
levels = np.arange(80,100,1)
extent = (delr/2., Lx - delr/2., Ly - delc/2., delc/2.)
plt.contour(head[0, :, :], levels=levels, extent=extent)

plt.subplot(1,2,2,aspect='equal')
hds = bf.HeadFile(modelname+'.hds')
head = hds.get_data(totim=1.0)
levels = np.arange(80,100,1)
extent = (delr/2., Lx - delr/2., Ly - delc/2., delc/2.)
plt.contour(head[3, :, :], levels=levels, extent=extent)
plt.show()

#Export to Model Muse
import shutil
shutil.copyfile(modelname+'.hds',modelname+'.bhd')