In Problem P5.3 from pages 246-247 in Anderson, Woessner and Hunt (2015), we are asked to construct an areal 2D model to represent an industrial facility in an arid setting. It is disposing of fluids in a 900 m by 900 m pond that is leaking at a rate of 0.2 m/d (Fig. P5.2). Recharge from precipitation in this area is negligible. The pond is located in the center of the horizontal problem domain and is underlain by a sequence of sediment layers consisting of sand, clay, and sand and gravel. Wet areas around the pond perimeter at the land surface are causing some local water-logging of the soils and impacting vegetation. The owners of the pond believe that the water-logging is caused by seepage out of the pond through the berms around the sides of the pond. The state regulatory agency, however, suspects that leakage through the bottom of the pond has created a water table mound that intersects the land surface. The objective of modeling is to determine whether the groundwater mound beneath the pond reaches the land surface and is water-logging the soil.
The consulting firm hired by the industrial facility recommends a 2D areal steady-state unconfined model as a quick and easy way to address the modeling objective. As a newly hired hydrogeologist of the consulting firm, you are instructed to construct the model. The width of the problem domain is 11,700 m. Use no flow boundary conditions at the north and south ends of the problem domain and specified heads along the sides (Fig. P5.2). Use a uniform nodal spacing of 900 m. Use Eqns (B5.3.2) and (B5.3.3) in Box 5.3 to compute the average horizontal and vertical hydraulic conductivity for the layer. Although the vertical hydraulic conductivity is not used in a one-layer 2D areal model, the vertical anisotropy ratio of the layer is of interest. Produce a contour map of the water table (use 1-m contour intervals) using the computed heads. Under this representation does the water table intersect the land surface?
Below is an iPython Notebook that builds a Python MODFLOW model for this problem and plots results. See the Github wiki associated with this Chapter for information on one suggested installation and setup configuration for Python and iPython Notebook.
[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 modelusing the Flopy website approach. 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 a new subdirectory underneath the directory where the notebook is located. 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 = "P5-3_Industrial_pond"
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 typically all caps except for layer thickness (b) and K direction (Kh or Kv); python syntax are lower case or mixed case
# This is an unconfined areal 2D model.
LX = 13500. # aquifer width of 11700 m + 2 900m nodes for constant head boundary condition on each end
LY = 9000. # height of aquifer of 9000 m, no flow boundaries do not require explicit cells in MODFLOW
ZTOP = 130. # the system is unconfined
ZBOT = 0.
NLAY = 1
NROW = 10
NCOL = 15
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)
# RCH = 0.0 #not needed for Problem P5.3
# WELLQ = 0. #not needed for Problem P5.3
POND_SEEP=0.2
print "DELR =", DELR, " DELC =", DELC, ' DELV =', DELV
print "BOTM =", BOTM
#print "Recharge =", RCH
print "Pond Seepage =", POND_SEEP, "m/d"
#print "Pumping well rate =", WELLQ
In this problem we are required to calculate the equivalent 1-layer hydraulic conductivity from the 3-layer system showed in Figure P5.2 using the correct equation from Box 5.3 (equation B5.3.2)
In [4]:
LAY1b=((120.-80)+(90-80))/2 #this calcultes the average saturated thickness of Layer 1
print "The average Layer 1 thickness =", LAY1b, "m"
This agrees with the value given in the caption in figure P5.2
In [5]:
LAY2b = 80.-40
LAY3b = 40.-0
print "Layer 2 thickness =", LAY2b, "m"
print "Layer 3 thickness =", LAY3b, "m"
In [6]:
#In equation B5.3.2 we also need the total thickess of the aquifer B
TOT_THICK = LAY1b+LAY2b+LAY3b
print "The total aquifer thickness B =", TOT_THICK, "m"
In [7]:
#Now we can assign the Kh for each layer from Figure P5.2
LAY1Kh=44.
LAY2Kh=0.13
LAY3Kh=113.
#now we can calculate the equivalent Kh from Box 5.3's equation B5.3.2
EQUIV_Kh=(LAY1Kh*LAY1b/TOT_THICK)+(LAY2Kh*LAY2b/TOT_THICK)+(LAY3Kh*LAY3b/TOT_THICK)
print "The 1-layer equivlant Kh =", EQUIV_Kh, "m/d"
Don't forget that we are asked to do a similar calculation for equivalent Kv. I will leave that for you to calculate because it is not needed for a 1-layer MODFLOW model.
Create a flopy MODFLOW object: flopy.modflow.Modflow.
In [8]:
# Assign name and create modflow model object
modelname = 'P5-3'
#exe_name = os.path.join(datapath, 'mf2005.exe') # for Windows OS
exe_name = os.path.join(datapath, 'mf2005') # for Mac OS
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 [9]:
# Create the discretization object
TOP = ZTOP * np.ones((NROW, NCOL),dtype=np.float)
In [10]:
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 [11]:
# 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] = -1 #don't forget arrays are zero-based! Sets first column
IBOUND[:, :, -1] = -1 # Sets last column
print IBOUND
In [12]:
STRT = 130 * np.ones((NLAY, NROW, NCOL), dtype=np.float32) # set starting head to landsurface (130 m) throughout model domain
STRT[:, :, 0] = 90. # leftmost constant head
STRT[:, :, -1] = 120. # rightmost constant head
print STRT
In [13]:
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 [14]:
LPF_PACKAGE = flopy.modflow.ModflowLpf(MF, laytyp=1, hk=EQUIV_Kh) # 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
In [15]:
#WEL_PACKAGE = flopy.modflow.ModflowWel(MF, stress_period_data=[0,0,0,WELLQ]) # remember python 0 index, layer 0 = layer 1 in MF
#print WEL_PACKAGE # uncomment this at far left to see the information about the flopy WEL object
Create a flopy output control object: flopy.modflow.ModflowOc.
In [16]:
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 [17]:
PCG_PACKAGE = flopy.modflow.ModflowPcg(MF, mxiter=500, iter1=100, hclose=1e-04, rclose=1e-03, relax=0.98, damp=0.5)
# print PCG_PACKAGE # uncomment this at far left to see the information about the flopy PCG object
Because the pond seepage is given as a flux over the cell (L/T) rather than a volumetric flow rate (L3/T), the Recharge Package is the easier to use than MODFLOW's Well Package. However, we only want to add pond seepage to the cells that represent the pond. Therefore, we'll need to make an array like we do for starting heads.
In [18]:
SEEP_ARRAY = 0 * np.ones((NROW, NCOL), dtype=np.float32) # set seepage to 0 over model grid
SEEP_ARRAY[4, 7] = POND_SEEP # add Pond seepage at the pond location
print SEEP_ARRAY
Create a flopy pcg package object: flopy.modflow.ModflowRch.
In [19]:
RCH_PACKAGE = flopy.modflow.ModflowRch(MF, rech=SEEP_ARRAY)
# print RCH_PACKAGE # uncomment this at far left to see the information about the flopy RCH object
In [20]:
#Before writing input, destroy all files in folder to prevent reusing old files
#Here's the working directory
print modelpath
#Here's what's currently in the working directory
modelfiles = os.listdir(modelpath)
print modelfiles
In [21]:
#delete these files to 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 [22]:
#Now write the model input files
MF.write_input()
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 [23]:
# return current working directory
print "You can check the newly created files in", modelpath
In [24]:
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 [25]:
#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 [26]:
#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()
Land surface elevation = 130 m, which is significantly higher than the maximum head in the model.
In [27]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(12,10))
#setup contour levels and plot extent
LEVELS = np.arange(90., 126., 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.set_xlabel("x")
AX1.set_ylabel("y")
YTICKS = np.arange(0, 9000, 1000)
AX1.set_yticks(YTICKS)
AX1.set_title("P5.3 1-layer Industrial Pond Problem")
AX1.text(6000, 5500, r"pond", fontsize=10, color="blue")
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')
AX2.set_xlabel("x")
AX2.set_ylabel("y")
AX2.set_yticks(YTICKS)
AX2.set_title("P5.3 color flood")
AX2.text(6000, 5500, r"pond", fontsize=10, color="white")
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.45)
In [28]:
#Using code from P4.4, let's plot a cross section of head in row = 4 from headobj, recall that Python is zero based
#so that MODFLOW row 5 with the pond is equal to Python row 4
#define Y as HEAD along the row and then print; ROW is a variable that allows us to change the row plotted easily
ROW = 4
Y = HEAD[0,ROW,:]
print Y
In [29]:
#in order to plot the cross section we'll need to create X-coordinates to match with heads at the node centers
XCOORD = np.arange(0, 13500, 900) + 450
print XCOORD
In [30]:
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(1, 1, 1)
TITLE = '1-layer model 900m x 900m grid: cross section of head along Row = ({0})'.format(ROW) #this allows the title to be updated as ROW changes
ax.set_title(TITLE)
ax.set_xlabel('X')
ax.set_ylabel('head')
ax.set_xlim(0, 13500.)
ax.set_ylim(0.,150.)
ax.text(6400,125, r"pond", fontsize=12, color="black")
ax.text(2030,132, r"land surface", fontsize=15, color="green")
ax.text(13200, 80, r"constant head = 120m", fontsize=10, color="blue",rotation='vertical')
ax.text(150, 80, r"constant head = 90m", fontsize=10, color="blue",rotation='vertical')
ax.plot(XCOORD, Y)
ax.plot(XCOORD, XCOORD*0+130) #land surface
Out[30]:
The state regulatory agency insists that a 3D model be developed to examine how vertical flow and anisotropy influence the height of the groundwater mound. They point out that the low hydraulic conductivity of layer 2 andthe vertical anisotropy present in the layered sequence of units might cause the mound to rise to the surface. Construct a three- layer steady-state model based on the information in Fig. P5.2. Specified head and no flow boundaries extend to all layers. Generate an equipotential surface (using 1-m contour intervals) for each layer. Also show the head distribution in a cross section that passes through the specified head boundaries and the pond.
In [31]:
#we have to redefine the layering
ZTOP = 130. # the system is unconfined so set the top above land surface so that the water table never > layer top
ZBOT = 0.
NLAY = 3
BOTM = np.zeros((NLAY, NROW, NCOL), dtype=np.float)
BOTM[0,:,:] = 80.
BOTM[1,:,:] = 40.
BOTM[2,:,:] = 0.
print BOTM
In [32]:
# Create the discretization object
TOP = ZTOP * np.ones((NROW, NCOL),dtype=np.float)
DIS_PACKAGE = flopy.modflow.ModflowDis(MF, NLAY, NROW, NCOL, delr=DELR, delc=DELC,
top=TOP, botm=BOTM, laycbd=0)
Now we must do the same for MODFLOW's Basic Package
In [33]:
IBOUND = np.ones((NLAY, NROW, NCOL), dtype=np.int)
IBOUND[:,:,0] = -1
IBOUND[:,:,-1] = -1
print IBOUND
In [34]:
STRT = 130. * np.ones((NLAY, NROW, NCOL), dtype=np.float)
STRT[:,:,0] = 90.
STRT[:,:,-1] = 120.
print STRT
In [35]:
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
Now we must do the same for hydraulic conductivity.
In [36]:
#Now we can assign the Kh for each layer from Figure P5.2
LAY1Kh=44.
LAY2Kh=0.13
LAY3Kh=113.
KH_ARRAY = np.zeros((NLAY, NROW, NCOL), dtype=np.float)
KH_ARRAY[0,:,:] = LAY1Kh
KH_ARRAY[1,:,:] = LAY2Kh
KH_ARRAY[2,:,:] = LAY3Kh
print KH_ARRAY
In [37]:
#Now we can assign the Kv for each layer from Figure P5.2
LAY1Kv=4.4
LAY2Kv=0.013
LAY3Kv=11.3
KV_ARRAY = np.zeros((NLAY, NROW, NCOL), dtype=np.float)
KV_ARRAY[0,:,:] = LAY1Kv
KV_ARRAY[1,:,:] = LAY2Kv
KV_ARRAY[2,:,:] = LAY3Kv
print KV_ARRAY
In [38]:
LPF_PACKAGE = flopy.modflow.ModflowLpf(MF, laytyp=1, hk=KH_ARRAY, vka=KV_ARRAY) # 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
In [39]:
#Before writing input, destroy all files in folder to prevent reusing old files
#Here's the working directory
print modelpath
#Here's what's currently in the working directory
modelfiles = os.listdir(modelpath)
print modelfiles
#delete these files to 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 [40]:
#Now write the model input files and run MODFLOW
MF.write_input()
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)
In [41]:
#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
#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()
Note the new maximum head value (recall land surface elevation = 130).
In [42]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(12,10))
#setup contour levels and plot extent
LEVELS = np.arange(90., 146., 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.set_xlabel("x")
AX1.set_ylabel("y")
YTICKS = np.arange(0, 9000, 1000)
AX1.set_yticks(YTICKS)
AX1.set_title("3-layer P5.3 Industrial Pond Problem 900x900m grid")
AX1.text(6000, 5500, r"pond", fontsize=10, color="blue")
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')
AX2.set_xlabel("x")
AX2.set_ylabel("y")
AX2.set_yticks(YTICKS)
AX2.set_title("P5.3 color flood")
AX2.text(6000, 5500, r"pond", fontsize=10, color="white")
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.45)
In [43]:
#Using code from P4.4, let's again plot a cross section of head in row = 4 from headobj, recall that Python is zero based
#so that MODFLOW row 5 with the pond is equal to Python row 4
#define Y as HEAD along the row and then print; ROW is a variable that allows us to change the row plotted easily
ROW = 4
Y = HEAD[0,ROW,:]
print Y
In [44]:
#in order to plot the cross section we'll need to create X-coordinates to match with heads at the node centers
#(we could have just used the XCOORD calculated above as the X-Y spacing did not change, but we'll do it again
#to make sure we have the correct X coordinates)
XCOORD = np.arange(0, 13500, 900) + 450
print XCOORD
In [45]:
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(1, 1, 1)
TITLE = '3-layer model 900m x 900m grid: cross section of head along Row = ({0})'.format(ROW) #this allows the title to be updated as ROW changes
ax.set_title(TITLE)
ax.set_xlabel('X')
ax.set_ylabel('head')
ax.set_xlim(0, 13500.)
ax.set_ylim(0.,150.)
ax.text(6150,125, r"pond", fontsize=12, color="black")
ax.text(2030,132, r"land surface", fontsize=15, color="green")
ax.text(2030,75, r"Layer 1 bottom", fontsize=10, color="red")
ax.text(2030,35, r"Layer 2 bottom", fontsize=10, color="teal")
ax.text(13200, 80, r"constant head = 120m", fontsize=10, color="blue",rotation='vertical')
ax.text(150, 80, r"constant head = 90m", fontsize=10, color="blue",rotation='vertical')
ax.plot(XCOORD, Y)
ax.plot(XCOORD, XCOORD*0+130) #land surface
ax.plot(XCOORD, XCOORD*0+80) #Bottom of Layer 1
ax.plot(XCOORD, XCOORD*0+40) #Bottom of Layer 2
Out[45]:
Examine your results and answer the following:
i. Explain why the results of the 2D model are different from the 3D model. What are the main factors that control the height of the water table mound under the pond? Discuss whether the 2D model is appropriate for this problem.
ii. Is it likely that the water table intersects the land surface away from the pond? If so, use shading on a map of the land surface in the vicinity of the pond to show the area affected by leakage.
When the modeling report is sent out for review, reviewers question whether the large nodal spacing of 900 m sufficiently captures the head gradient that defines the mound. They say that the surface area affected by the mound is underestimated. Use the three-layer model developed in (b) to assess the effect of nodal spacing on the solution. Reduce the nodal spacing uniformly over the grid/mesh to 300 m, or construct an irregular FD grid, unstructured FD grid, or FE mesh with fine nodal spacing in the vicinity of the pond. Run the model and generate equipotential maps for each layer (use a 1-m contour interval) and a cross section that passes through the pond and constant head boundaries. If the mound intersects the land surface, show the area impacted by shading on a map of the land surface. Compare and contrast results with those of parts (a) and (b).
In [46]:
# let's redefine the grid
LX = 12300. # aquifer width of 11700 m + 2 300m nodes for constant head boundary condition on each end
LY = 9000. # height of aquifer of 9000 m, no flow boundaries do not require explicit cells in MODFLOW
ZTOP = 130. # the system is unconfined
NROW = 30 #these two lines are the only things we need to change to refine
NCOL = 41 # the grid
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)
POND_SEEP=0.2
print "DELR =", DELR, " DELC =", DELC
print "Pond Seepage =", POND_SEEP, "m/d"
In [47]:
#repeat the modeling building steps from P5.3b to make sure all model definition is current
NLAY = 3
ZBOT = 0.
BOTM = np.zeros((NLAY, NROW, NCOL), dtype=np.float)
BOTM[0,:,:] = 80.
BOTM[1,:,:] = 40.
BOTM[2,:,:] = 0.
#print BOTM
# Create the discretization object
TOP = ZTOP * np.ones((NROW, NCOL),dtype=np.float)
DIS_PACKAGE = flopy.modflow.ModflowDis(MF, NLAY, NROW, NCOL, delr=DELR, delc=DELC,
top=TOP, botm=BOTM, laycbd=0)
IBOUND = np.ones((NLAY, NROW, NCOL), dtype=np.int)
IBOUND[:,:,0] = -1
IBOUND[:,:,-1] = -1
#print IBOUND
STRT = 130. * np.ones((NLAY, NROW, NCOL), dtype=np.float)
STRT[:,:,0] = 90.
STRT[:,:,-1] = 120.
#print STRT
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
#Now we can assign the Kh for each layer from Figure P5.2
LAY1Kh=44.
LAY2Kh=0.13
LAY3Kh=113.
KH_ARRAY = np.zeros((NLAY, NROW, NCOL), dtype=np.float)
KH_ARRAY[0,:,:] = LAY1Kh
KH_ARRAY[1,:,:] = LAY2Kh
KH_ARRAY[2,:,:] = LAY3Kh
#print KH_ARRAY
#Now we can assign the Kv for each layer from Figure P5.2
LAY1Kv=4.4
LAY2Kv=0.013
LAY3Kv=11.3
KV_ARRAY = np.zeros((NLAY, NROW, NCOL), dtype=np.float)
KV_ARRAY[0,:,:] = LAY1Kv
KV_ARRAY[1,:,:] = LAY2Kv
KV_ARRAY[2,:,:] = LAY3Kv
#print KV_ARRAY
LPF_PACKAGE = flopy.modflow.ModflowLpf(MF, laytyp=1, hk=KH_ARRAY, vka=KV_ARRAY) # 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
#We have to place the pond in the new grid too via MODFLOW's Recharge Package (like in P5.3a)
SEEP_ARRAY = 0 * np.ones((NROW, NCOL), dtype=np.float32) # set seepage to 0 over model grid
SEEP_ARRAY[12, 20] = POND_SEEP # add Pond seepage at the pond location (have 9 nodes now with finer grid)
SEEP_ARRAY[13, 20] = POND_SEEP
SEEP_ARRAY[14, 20] = POND_SEEP
SEEP_ARRAY[12, 21] = POND_SEEP
SEEP_ARRAY[13, 21] = POND_SEEP
SEEP_ARRAY[14, 21] = POND_SEEP
SEEP_ARRAY[12, 22] = POND_SEEP
SEEP_ARRAY[13, 22] = POND_SEEP
SEEP_ARRAY[14, 22] = POND_SEEP
RCH_PACKAGE = flopy.modflow.ModflowRch(MF, rech=SEEP_ARRAY)
In [48]:
#Before writing input, destroy all files in folder to prevent reusing old files
#Here's the working directory
print modelpath
#Here's what's currently in the working directory
modelfiles = os.listdir(modelpath)
print modelfiles
#delete these files to 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 [49]:
#Now we are ready to write input and re-run MODFLOW with the 300 x 300 m grid
MF.write_input()
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)
In [50]:
#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
#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 [51]:
#Create a contour plot of heads
FIG = plt.figure(figsize=(12,10))
#setup contour levels and plot extent
LEVELS = np.arange(90., 146., 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.set_xlabel("x")
AX1.set_ylabel("y")
YTICKS = np.arange(0, 9000, 1000)
AX1.set_yticks(YTICKS)
AX1.set_title("3-layer P5.3 Industrial Pond Problem 300x300m grid")
AX1.text(6000, 5500, r"pond", fontsize=10, color="blue")
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')
AX2.set_xlabel("x")
AX2.set_ylabel("y")
AX2.set_yticks(YTICKS)
AX2.set_title("P5.3 color flood")
AX2.text(6000, 5500, r"pond", fontsize=10, color="white")
cax = AX2.imshow(HEAD[0, :, :], extent=EXTENT, interpolation='nearest')
cbar = FIG.colorbar(cax, orientation='vertical', shrink=0.45)
In [52]:
#Using code from P4.4, let's look at the head in column = 13 from headobj, and then plot it
#print HEAD along a column; COL is a variable that allows us to change this easily
ROW = 13
Y = HEAD[0,ROW,:]
print Y
In [53]:
#for our cross section create X-coordinates to match with heads
XCOORD = np.arange(0, 12300, 300) + 150
print XCOORD
In [54]:
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(1, 1, 1)
TITLE = '3-layer model 300m x 300m grid: cross section of head along Row = ({0})'.format(ROW) #this allows the title to be updated as ROW changes
ax.set_title(TITLE)
ax.set_xlabel('X')
ax.set_ylabel('head')
ax.set_xlim(0, 12300.)
ax.set_ylim(0.,150.)
ax.text(6150,125, r"pond", fontsize=12, color="black")
ax.text(2030,132, r"land surface", fontsize=15, color="green")
ax.text(2030,75, r"Layer 1 bottom", fontsize=10, color="red")
ax.text(2030,35, r"Layer 2 bottom", fontsize=10, color="teal")
ax.text(12000, 80, r"constant head = 120m", fontsize=10, color="blue",rotation='vertical')
ax.text(150, 80, r"constant head = 90m", fontsize=10, color="blue",rotation='vertical')
ax.plot(XCOORD, Y)
ax.plot(XCOORD, XCOORD*0+130) #land surface
ax.plot(XCOORD, XCOORD*0+80) #Bottom of Layer 1
ax.plot(XCOORD, XCOORD*0+40) #Bottom of Layer 2
Out[54]:
Compare and contrast results with those of parts (a) and (b).
What is your conclusion about the cause of water-logging? Support your answer with modeling results.
In this problem we assumed that all parameters were known, and we examined how layering (and anisotropy) and nodal spacing affected the height and configuration of the water table. Make a list of factors (in addition to dimensionality and nodal spacing) that would add uncertainty to a forecast of the height of the water table mound in this problem. Speculate how the uncertainty might affect modeling results and conclusions.
In [ ]: