Import the required libraries

This tutorial requires some Python core libraries, Scipy libraries (numpy and matplotlib) as well as the Flopy library. For the management and data extraction from raster files we use the GDAL library.

import flopy
import os
import flopy.utils.binaryfile as bf
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt

Create a basic MODFLOW model object

As well as any other model generated by Flopy, first we have to setup the model name and the working directory. We strongly recommend to follow the folder/file estructure provided on the "Input Files" of this tutorial. Based on the complex geometry and the differences among layers elevations and layer thickness, the NWT solver has been applied to the model.

modelname = "model1"
modelpath = "../Model/"
#Modflow Nwt
mf1 = flopy.modflow.Modflow(modelname, exe_name="../Exe/MODFLOW-NWT_64.exe", version="mfnwt",model_ws=modelpath)
nwt = flopy.modflow.ModflowNwt(mf1 ,maxiterout=150,linmeth=2)

Open and read raster files

In this part we will use the GDAL library to open the rasters of elevation and boundary conditions and get their geotransformation parameters. The tutorial requires that the raster resolution is the same as the grid resolution. Raster bands are read as Numpy arrays and inserted to the model as model top.

#Raster paths
demPath ="../Rst/DEM_200.tif"
crPath = "../Rst/CR.tif"

#Open files
demDs =gdal.Open(demPath)
crDs = gdal.Open(crPath)
geot = crDs.GetGeoTransform() #Xmin, deltax, ?, ymax, ?, delta y

geot

(613882.5068, 200.0, 0.0, 8375404.9524, 0.0, -200.0)

# Get data as arrays
demData = demDs.GetRasterBand(1).ReadAsArray()
crData = crDs.GetRasterBand(1).ReadAsArray()

# Get statistics
stats = demDs.GetRasterBand(1).GetStatistics(0,1) 
stats

[3396.0, 4873.0, 4243.0958193041, 287.90182492684]

Spatial discretization

Since this tutorial deals with a regional scale, the differences in elevation on the area of study can be more than a hundred meters or more than a thousand meters if we deal with andean basins. The hydrogeological conceptualization for this tutorial is the definition of two upper thin layers (Layer1 and Layer1_2) and thicker layers below.

#Elevation of layers
Layer1 = demData - 7.5
Layer1_2 = demData - 15
Layer2 = (demData - 3000)*0.8 + 3000
Layer3 = (demData - 3000)*0.5 + 3000
Layer4 = (demData>0)*3000

Column and row number and dimension come from the raster attributes

#Boundaries for Dis = Create discretization object, spatial/temporal discretization
ztop = demData
zbot = [Layer1, Layer1_2, Layer2, Layer3, Layer4]
nlay = 5
nrow = demDs.RasterYSize
ncol = demDs.RasterXSize
delr = geot[1]
delc = abs(geot[5])

Definition of the flow packages for the MODFLOW model

In this part we define the packages related to groundwater flow on the MODFLOW model. First we define the DIS package that has the geometry as well as the simulation type (steady / transient). The model run on steady conditions.

dis = flopy.modflow.ModflowDis(mf1, nlay,nrow,ncol,delr=delr,delc=delc,top=ztop,botm=zbot,itmuni=1)

Then we define another MODEL packages. Some boundary conditions are setup from raster valures, while other packages use constant value. These are the required packages for the simulation:

• the BAS package for setting the constant head cells and inactive cells,
• the UPW that defines the vertical / horizontal hydraulic conductivity,
• the RCH package that applies rechage from precipitation
• the EVT package for the simulation of groundwater consumption from roots
• the DRN package for the groundwater discharge to the channel network
• the OC packages for the output record

# Variables for the BAS package
iboundData = np.zeros(demData.shape, dtype=np.int32)
iboundData[crData > 0 ] = 1
bas = flopy.modflow.ModflowBas(mf1,ibound=iboundData,strt=ztop, hnoflo=-2.0E+020)

# Array of hydraulic heads per layer
hk = [1E-4, 1E-5, 1E-7, 1E-8, 1E-9]

# Add UPW package to the MODFLOW model
upw = flopy.modflow.ModflowUpw(mf1, laytyp = [1,1,1,1,0], hk = hk)

#Add the recharge package (RCH) to the MODFLOW model
rch = np.ones((nrow, ncol), dtype=np.float32) * 0.120/365/86400
rch_data = {0: rch}
rch = flopy.modflow.ModflowRch(mf1, nrchop=3, rech =rch_data)

#Add the evapotranspiration package (EVT) to the MODFLOW model
evtr = np.ones((nrow, ncol), dtype=np.float32) * 1/365/86400
evtr_data = {0: evtr}
evt = flopy.modflow.ModflowEvt(mf1,nevtop=1,surf=ztop,evtr=evtr_data, exdp=0.5)

#Add the drain package (DRN) to the MODFLOW model
river = np.zeros(demData.shape, dtype=np.int32)
river[crData == 3 ] = 1
list = []
for i in range(river.shape[0]):
    for q in range(river.shape[1]):
        if river[i,q] == 1:
            list.append([0,i,q,ztop[i,q],0.001]) #layer,row,column,elevation(float),conductanceee
rivDrn = {0:list}

drn = flopy.modflow.ModflowDrn(mf1, stress_period_data=rivDrn)

# Add OC package to the MODFLOW model
oc = flopy.modflow.ModflowOc(mf1) #ihedfm= 1, iddnfm=1

Write files of the MODFLOW model and run simulation

Write the .nam file and the files declared on the .nam file. Then runs the simulation and shows the simulation information on the screen

#Write input files -> write file with extensions
mf1.write_input()

#run model -> gives the solution
mf1.run_model()

FloPy is using the following  executable to run the model: ../Exe/MODFLOW-NWT_64.exe

                                  MODFLOW-NWT-SWR1 
    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUNDWATER-FLOW MODEL
                             WITH NEWTON FORMULATION
                             Version 1.1.4 4/01/2018                         
                    BASED ON MODFLOW-2005 Version 1.12.0 02/03/2017                       

                    SWR1 Version 1.04.0 09/15/2016                       

 Using NAME file: model1.nam 
 Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/02/22 10:19:24

 Solving:  Stress period:     1    Time step:     1    Groundwater-Flow Eqn.
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/02/22 10:20:04
 Elapsed run time: 39.475 Seconds

  Normal termination of simulation





(True, [])

Model results post-processing

This tutorial has some model features and output representation done with Flopy and Matplotlib. Another 3D representation on the VTK format is performed on another notebook from this tutorial

#Import model
ml = flopy.modflow.Modflow.load('../Model/model1.nam')

Model grid representation

# First step is to set up the plot
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')

# Next we create an instance of the ModelMap class
modelmap = flopy.plot.ModelMap(sr=ml.dis.sr)

# Then we can use the plot_grid() method to draw the grid
# The return value for this function is a matplotlib LineCollection object,
# which could be manipulated (or used) later if necessary.
linecollection = modelmap.plot_grid(linewidth=0.4)

Cross section of model grid representation

fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1)
# Next we create an instance of the ModelCrossSection class
#modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Column': 5})
modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Row': 99})

# Then we can use the plot_grid() method to draw the grid
# The return value for this function is a matplotlib LineCollection object,
# which could be manipulated (or used) later if necessary.
linecollection = modelxsect.plot_grid(linewidth=0.4)
t = ax.set_title('Column 6 Cross-Section - Model Grid')

Active/inactive cells on model extension

fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=0)
quadmesh = modelmap.plot_ibound(color_noflow='cyan')
linecollection = modelmap.plot_grid(linewidth=0.4)

Cross sections of active/inactive cells

fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1)
modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Column': 5})
patches = modelxsect.plot_ibound(color_noflow='cyan')
linecollection = modelxsect.plot_grid(linewidth=0.4)
t = ax.set_title('Column 6 Cross-Section with IBOUND Boundary Conditions')

Channel network as drain (DRN) package

fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=-14)
quadmesh = modelmap.plot_ibound(color_noflow='cyan')
quadmesh = modelmap.plot_bc('DRN', color='blue')
linecollection = modelmap.plot_grid(linewidth=0.4)

Model grid and heads representation

fname = os.path.join(modelpath, 'model1.hds')
hdobj = flopy.utils.HeadFile(fname)
head = hdobj.get_data()

fig = plt.figure(figsize=(30, 30))
ax = fig.add_subplot(1, 2, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=-14)
#quadmesh = modelmap.plot_ibound()
quadmesh = modelmap.plot_array(head, masked_values=[-2.e+20], alpha=0.8)
linecollection = modelmap.plot_grid(linewidth=0.2)