Determination of Flow Direction Vectors from a MODFLOW Model with Python and Flopy - Tutorial


Groundwater flow direction representation is useful to understand the actual and predicted conditions of the groundwater flow regime. The arrow direction and magnitude give a quick perspective of the main groundwater flow directions and the interconexion between sources and discharge points. Flow direction vectors wasn´t a very common feature in groundwater modeling with the open source stack until the development of Flopy library that has a special tool for quiver plot representation.

This tutorial show the complete workflow to determine the flow directions from a MODFLOW model done with Model Muse. The scripting insert a background image, georeference the model from parameters exported as comments, and export the resulting figure as a PNG file. The tutorial is done in Python 3 on a Jupyter Notebook.



This is the complete Python script for this tutorial:

Import the required libraries

From numpy, matplotlib and flopy as well as Python core libraries.

%matplotlib inline
import sys, os, re
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import flopy
import flopy.utils.binaryfile as bf

Set a background image for the figure

This is not mandatory, but it helps to locate and analyse groundwater flow with respect to geomorphological and geographical features. Backgroundmap can be generated from QGIS by exporting the canvas that actually comes with the world file (.pgw or similar), this world file allows the georeferenciation of the image. There are options to improve the graphics with contour lines and the representation of boundary condition, please check the flopy examples:

wrlFile = open('../Rst/BackgroundImage.pgw').read().split('\n')
cellsize = float(wrlFile[0])
bcgImg = mpimg.imread('../Rst/BackgroundImage.png')
left = float(wrlFile[4])
right = float(wrlFile[4]) + bcgImg.shape[1]*cellsize
bottom = float(wrlFile[5]) - bcgImg.shape[0]*cellsize
top = float(wrlFile[5])

Create our modflow model as an object

Basic setup of a modflow model in flopy

modPath = '../Model/'
modName = 'SingleLayerModel'
exeName = '../Exe/mf2005.exe'   
ml = flopy.modflow.Modflow.load(modName+'.nam', model_ws=modPath, 

Import some geospatial parameters from the model

Model Muse inserts some geospatial parameters as comments on the DIS file. We have to import it by opening the DIS file and apply the parameters to the model object. Note: The EPSG code have to be set by hand since it is not written on the DIS file.

#Open the DIS file and extract the spatial parametes
spatialpar = open(modPath + modName + '.dis').readlines()[:6]
ul =  re.findall("\d+\.\d+", spatialpar[1])
lr =  re.findall("\d+\.\d+", spatialpar[4])
xul, yul = float(ul[0]), float(ul[1])
xlr, ylr = float(ul[0]), float(ul[1])
rotation = -float(re.findall("\d+", spatialpar[5])[0])

#Assign parameters to the model object, yul=yul, rotation=-rotation) 'EPSG:32718' #### Defined by the user
xul:199558.1423; yul:8807344.096; rotation:0; proj4_str:+init=EPSG:32718; units:meters; lenuni:2; length_multiplier:1.0

Import the output heads

h = flopy.utils.formattedfile.FormattedHeadFile(modPath+modName+'.fhd',model=ml)
head = h.get_data(totim=1)

head[0,30:33,30:33] #Sample of the heads values
array([[24.3446, 26.4447, 28.8097],
       [23.7384, 25.8897, 28.1835],
       [23.2935, 25.5388, 27.3164]], dtype=float32)

Create a figure

The figure is generated as a matplotlib figure, where flow direction vectors overlay the background image. There are options to filter the number of arrows with the istep/jstep, this is highly recommended to enhance the representation and user experience.

fig = plt.figure(figsize=(16, 10))
ax = fig.add_subplot(1, 1, 1, aspect='equal')

modelmap = flopy.plot.ModelMap(model=ml)
cbb = bf.CellBudgetFile(modPath+modName+'.cbc')

frf = cbb.get_data(text='FLOW RIGHT FACE', totim=1)[0]
fff = cbb.get_data(text='FLOW FRONT FACE', totim=1)[0]

lc = modelmap.plot_grid(linewidth=1, color='gray',alpha=0.3)
vectors = modelmap.plot_discharge(frf, fff, head=head,istep=2,jstep=2,normalize=False,color='cyan')

image = plt.imshow(bcgImg, extent=(left, right, bottom, top), alpha=0.7)

ax.set_xlabel('Este', size=14)
ax.set_ylabel('Norte', size=14)


Input file

You can download the requiered files for this tutorial on the following link.

Saul Montoya

Saul Montoya es Ingeniero Civil graduado de la Pontificia Universidad Católica del Perú en Lima con estudios de postgrado en Manejo e Ingeniería de Recursos Hídricos (Programa WAREM) de la Universidad de Stuttgart con mención en Ingeniería de Aguas Subterráneas y Hidroinformática.

Smiley face

Subscribe to our free e-newsletter for tutorials, articles, webminars, courses and more.