In Problem P4.2 from page 172-173 in Anderson, Woessner and Hunt (2015), we are asked to construct a profile model to compute groundwater flow under a dam and solve the model using either an FD or FE code. An impermeable concrete dam 60 m long and 40 m wide is constructed over an isotropic and homogeneous silty sand (K = 10 m/d) that is 26.5 m thick and fills a river valley underlain by impermeable bedrock (Fig. P4.2). The depth of the reservoir pool is 30 m and the tail water stream elevation below the dam is 10 m. The crosssectional area of interest is 200 m long centered on the dam. As you formulate your numerical model think carefully about how to represent boundary conditions.
In this notebook, we will work through the problem using MODFLOW and the Python tool set Flopy. Notice how much code is reused from P4.1 because the variable names (all CAPS) remained the same.
<br > [Acknowledgements: This tutorial was created by Randy Hunt and all failings are mine. The exercise here has benefited greatly from the online Flopy tutorial and example notebooks developed by Chris Langevin and Joe Hughes for the USGS Spring 2015 Python Training course GW1774]
In this example, we will create a simple groundwater flow model by following the tutorial included on the Flopy website. We will make a few small changes so that the tutorial works with our file structure.
Visit the tutorial website here.
In [1]:
%matplotlib inline
import sys
import os
import shutil
import numpy as np
from subprocess import check_output
# Import flopy
import flopy
For this tutorial, we will work in the 21_FlopyIntro directory, which is located up one folder and over to the Data folder. We can use some fancy Python tools to help us manage the directory creation. Note that if you encounter path problems with this workbook, you can stop and then restart the kernel and the paths will be reset.
In [2]:
# Set the name of the path to the model working directory
dirname = "P4-2_DamProfile"
datapath = os.getcwd()
modelpath = os.path.join(datapath, dirname)
print 'Name of model path: ', modelpath
# Now let's check if this directory exists. If not, then we will create it.
if os.path.exists(modelpath):
print 'Model working directory already exists.'
else:
print 'Creating model working directory.'
os.mkdir(modelpath)
In [3]:
# model domain and grid definition
# for clarity, user entered variables are all caps; python syntax are lower case or mixed case
# we will use a layer orientation profile for easy plotting (see Box 4.2 on page 126)
LX = 200.
LY = 26.5
ZTOP = 1. # the "thickness" of the profile will be 1 m (= ZTOP - ZBOT)
ZBOT = 0.
NLAY = 1
NROW = 53
NCOL = 200
DELR = LX / NCOL # recall that MODFLOW convention is DELR is along a row, thus has items = NCOL; see page XXX in AW&H (2015)
DELC = LY / NROW # recall that MODFLOW convention is DELC is along a column, thus has items = NROW; see page XXX in AW&H (2015)
DELV = (ZTOP - ZBOT) / NLAY
BOTM = np.linspace(ZTOP, ZBOT, NLAY + 1)
HK = 10.
VKA = 1.
print "DELR =", DELR, " DELC =", DELC, ' DELV =', DELV
print "BOTM =", BOTM
Create a flopy MODFLOW object: flopy.modflow.Modflow.
In [4]:
# Assign name and create modflow model object
modelname = 'P4-2'
exe_name = os.path.join(datapath, 'mf2005')
print 'Model executable: ', exe_name
MF = flopy.modflow.Modflow(modelname, exe_name=exe_name, model_ws=modelpath)
Create a flopy discretization package object: flopy.modflow.ModflowDis.
In [5]:
# Create the discretization object
TOP = np.ones((NROW, NCOL),dtype=np.float)
In [6]:
DIS_PACKAGE = flopy.modflow.ModflowDis(MF, NLAY, NROW, NCOL, delr=DELR, delc=DELC,
top=TOP, botm=BOTM[1:], laycbd=0)
# print DIS_PACKAGE #uncomment this on far left to see information about the flopy object
Create a flopy basic package object: flopy.modflow.ModflowBas.
In [7]:
# Variables for the BAS package
IBOUND = np.ones((NLAY, NROW, NCOL), dtype=np.int32) # all nodes are active (IBOUND = 1)
# make the top of the profile specified head by setting the IBOUND = -1
IBOUND[:, 0, 0:81] = -1 #don't forget arrays are zero-based!
IBOUND[:, 0, 121:200] = -1 #don't forget arrays are zero-based!
print IBOUND
In [8]:
STRT = 40 * np.ones((NLAY, NROW, NCOL), dtype=np.float32) # set starting head to 40 through out model domain
STRT[:, 0, 0:81] = 76.5 # pool elevation upstream of the dam
STRT[:, 0, 121:200] = 36.5 # pool elevation downstream of the dam
print STRT
In [9]:
BAS_PACKAGE = flopy.modflow.ModflowBas(MF, ibound=IBOUND, strt=STRT)
# print BAS_PACKAGE # uncomment this at far left to see the information about the flopy BAS object
Create a flopy layer property flow package object: flopy.modflow.ModflowLpf.
In [10]:
LPF_PACKAGE = flopy.modflow.ModflowLpf(MF, hk=HK, vka=VKA) # we defined the K and anisotropy at top of file
# print LPF_PACKAGE # uncomment this at far left to see the information about the flopy LPF object
Create a flopy output control object: flopy.modflow.ModflowOc.
In [11]:
OC_PACKAGE = flopy.modflow.ModflowOc(MF) # we'll use the defaults for the model output
# print OC_PACKAGE # uncomment this at far left to see the information about the flopy OC object
Create a flopy pcg package object: flopy.modflow.ModflowPcg.
In [12]:
PCG_PACKAGE = flopy.modflow.ModflowPcg(MF) # we'll use the defaults for the PCG solver
# print PCG_PACKAGE # uncomment this at far left to see the information about the flopy PCG object
In [13]:
#Before writing input, destroy all files in folder
#This will prevent us from reading old results
modelfiles = os.listdir(modelpath)
for filename in modelfiles:
f = os.path.join(modelpath, filename)
if modelname in f:
try:
os.remove(f)
print 'Deleted: ', filename
except:
print 'Unable to delete: ', filename
In [14]:
#Now write the model input files
MF.write_input()
Yup. It's that simple, the model datasets are written using a single command (mf.write_input).
Check in the model working directory and verify that the input files have been created. Or if you might just add another cell, right after this one, that prints a list of all the files in our model directory. The path we are working in is returned from this next block.
In [15]:
# return current working directory
modelpath
Out[15]:
In [16]:
silent = False #Print model output to screen?
pause = False #Require user to hit enter? Doesn't mean much in Ipython notebook
report = True #Store the output from the model in buff
success, buff = MF.run_model(silent=silent, pause=pause, report=report)
To read heads from the MODFLOW binary output file, we can use the flopy.utils.binaryfile module. Specifically, we can use the HeadFile object from that module to extract head data arrays.
In [17]:
#imports for plotting and reading the MODFLOW binary output file
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
#Create the headfile object and grab the results for last time.
headfile = os.path.join(modelpath, modelname + '.hds')
headfileobj = bf.HeadFile(headfile)
#Get a list of times that are contained in the model
times = headfileobj.get_times()
print 'Headfile (' + modelname + '.hds' + ') contains the following list of times: ', times
In [18]:
#Get a numpy array of heads for totim = 1.0
#The get_data method will extract head data from the binary file.
HEAD = headfileobj.get_data(totim=1.0)
#Print statistics on the head
print 'Head statistics'
print ' min: ', HEAD.min()
print ' max: ', HEAD.max()
print ' std: ', HEAD.std()
In [19]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(15,15))
#setup contour levels and plot extent
LEVELS = np.arange(36., 77., 1.)
EXTENT = (DELR/2., LX - DELR/2., DELC/2., LY - DELC/2.)
print 'Contour Levels: ', LEVELS
print 'Extent of domain: ', EXTENT
#Make a contour plot on the first axis
AX1 = FIG.add_subplot(1, 2, 1, aspect='equal')
AX1.contour(np.flipud(HEAD[0, :, :]), levels=LEVELS, extent=EXTENT)
#Make a color flood on the second axis
AX2 = FIG.add_subplot(1, 2, 2, aspect='equal')
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.25)
Look at the bottom of the MODFLOW output file (ending with a *.list) and note the water balance reported. It's that simple, the model automates the calculation of the water budget so the modeler does not have to.
In [ ]: