Triangular Mesh for Groundwater Models with MODFLOW 6 and Flopy - Tutorial

TriangularMeshModelingMODFLOW6.png

MODFLOW 6 is the latest version of the groundwater modeling code MODFLOW developed by the USGS. This version includes some innovate tools and a complete rearrange of the model construction, but to date this code doesn’t have any commercial / open source Graphical User Interface (GUI) that improve the code usage for beginner and intermediate groundwater modelers.

One of the most exceptional new features from MODFLOW 6 is the different discretization options for the model mesh generation. Options range from regular grid (same as MODFLOW 2005), triangular mesh and unstructured grid. Flopy that is a Python library for the build and simulation of MODFLOW 6 and other models has tools for triangular mesh generation. The workflow on groundwater modeling with MODFLOW 6 and Flopy for triangular mesh models is pretty fluid and we see a lot of potential for local and regional groundwater flow modeling.

This tutorial shows the complete process to create a triangular mesh with the utilities from Flopy and incorporate it to a MODFLOW 6 model. The model is simulated and results are represented as colored mesh and contour lines.


Tutorial


Python code

This is the script for the tutorial:

Import required libraries

This tutorial requires some Python core libraries, Numpy, Matplotlib and Flopy.

import os
import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict

import flopy
from flopy.utils.triangle import Triangle as Triangle

flopy is installed in E:\Software\Anaconda3\lib\site-packages\flopy

Define model file path and executable

The folder for the model files and simulation output is defined as well as executables for the triangular mesh generator and Modflow 6.

workspace = '../Model/'
triExeName = '../Exe/triangle.exe'
mf6ExeName = '../Exe/mf6.exe'

Define domain and drain location

The active model domain is a trapezoid with the smaller parallel side located at the model east side. The drain was conceptualized as a rectangle on the center left part of the model.

active_domain = [(0, 0), (100, 10), (100, 55), (0, 65)]
drain_domain= [(80, 31), (95, 31), (95, 34), (80, 34)]

Tringle mesh creation

We use the Triangle tool from the flopy utils that actually implement the Triangle program to generate meshes. This tutorial comes with the Triangle executable for Windows, if you use another operating system it is needed to compile the source code that you can find from https://www.cs.cmu.edu/~quake/triangle.html.

tri = Triangle(angle=30, model_ws=workspace, exe_name=triExeName)
tri.add_polygon(active_domain)
tri.add_polygon(drain_domain)
tri.add_region((10,10),0,maximum_area=30) #coarse discretization
tri.add_region((88,33),1,maximum_area=5) #fine discretization
tri.build()

Generated triangular grid and boundary representation

The next step is the mesh representation and a identification of the boundaries create from the mesh generation

fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(edgecolor='gray')

<matplotlib.collections.PatchCollection at 0x9005c88>
output_9_1.png
### Indentification and representation of grid boundaries

fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(edgecolor='gray')
numberBoundaries = tri.get_boundary_marker_array().max()+1

cmap = plt.cm.get_cmap("hsv", numberBoundaries)

labelList = list(range(1,numberBoundaries))
i = 0
for ibm in range(1,numberBoundaries):
    tri.plot_boundary(ibm=ibm, ax=ax,marker='o', color=cmap(ibm), label= ibm)

handles, labels = plt.gca().get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys(), loc='lower right')

<matplotlib.legend.Legend at 0x91c6ac8>
output_10_1.png

Representation of triangle features

By using the available methods for the triangle object "tri" we can represent the vertex and the centroid with their position and label.

fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(ax=ax, edgecolor='gray')
tri.plot_vertices(ax=ax, marker='o', color='blue')
tri.label_vertices(ax=ax, fontsize=10, color='blue')
tri.plot_centroids(ax=ax, marker='o', color='red')
tri.label_cells(ax=ax, fontsize=10, color='red')
output_12_0.png

Configuration of a MODFLOW 6 Model

Once we have the triangular mesh we can create a MODFLOW 6 model. As you might know, MODFLOW 6 has a difference in between a model and a simulation since it implements the concept of exchange in between models. A simulation can have many models, for this case we have a simulation of just one model.

name = 'mf'
sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
                             exe_name=mf6ExeName,
                             sim_ws=workspace)
tdis = flopy.mf6.ModflowTdis(sim, time_units='SECONDS',
                             perioddata=[[1.0, 1, 1.]])
gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
ims = flopy.mf6.ModflowIms(sim, print_option='SUMMARY', complexity='complex', 
                           outer_hclose=1.e-8, inner_hclose=1.e-8)

Configuration of the triangular discretization DISV package

MODFLOW 6 differenciates the spatial discretization file and the temporal discretization file. On the spatial discretization options we have three grid types: regular (.dis), triangular (.disv) and unstructured (*.disu). This study case uses a triangular discretization.

cell2d = tri.get_cell2d()
vertices = tri.get_vertices()
xcyc = tri.get_xcyc()
nlay = 1
ncpl = tri.ncpl
nvert = tri.nvert
top = 1.
botm = [0.]

dis = flopy.mf6.ModflowGwfdisv(gwf, length_units='METERS',
                               nlay=nlay, ncpl=ncpl, nvert=nvert,
                               top=top, botm=botm, 
                               vertices=vertices, cell2d=cell2d)

Configuration of the flow package and the initial conditions

npf = flopy.mf6.ModflowGwfnpf(gwf, k=0.0001)
ic = flopy.mf6.ModflowGwfic(gwf, strt=10.0)

Definition of the Constant Head CHD boundary conditions

The left and right sides of the model were defined as CHD boundary condition to represent a regional flow right to left from an hydraulic head of 30m to 10m. See that the edge vertices are obtained from the model grid edges.

chdList = []

leftcells = tri.get_edge_cells(4)
rightcells = tri.get_edge_cells(2)

for icpl in leftcells: chdList.append([(0, icpl), 30])
for icpl in rightcells: chdList.append([(0, icpl), 10])

chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdList)

Definition of the Drain DRN boundary conditions

At the center left part of the model there is a drain conceptualized as two lines on the border from a previously declared region. The drain has a elevation of 2.5 meters lower that the eastern constant head boundary.

drnList = []

drnCells = tri.get_edge_cells(5) + tri.get_edge_cells(7)

for icpl in drnCells: drnList.append([(0, icpl), 7.5, 0.10])

chd = flopy.mf6.ModflowGwfdrn(gwf, stress_period_data=drnList)

Define output options, write input files and run the model

The tutorial defines the data type that will be stored and the records on the listing file. Then Modflow 6 model files are written on the model folder and simulation is performed.

oc = flopy.mf6.ModflowGwfoc(gwf,
                            budget_filerecord='{}.cbc'.format(name),
                            head_filerecord='{}.hds'.format(name),
                            saverecord=[('HEAD', 'LAST'),
                                        ('BUDGET', 'LAST')],
                            printrecord=[('HEAD', 'LAST'),
                                         ('BUDGET', 'LAST')])

sim.write_simulation()
success, buff = sim.run_simulation()

writing simulation...
  writing simulation name file...
  writing simulation tdis package...
  writing ims package ims_-1...
  writing model mf...
    writing model name file...
    writing package disv...
    writing package npf...
    writing package ic...
    writing package chd_0...
INFORMATION: maxbound in ('gwf6', 'chd', 'dimensions') changed to 19 based on size of stress_period_data
    writing package drn_0...
INFORMATION: maxbound in ('gwf6', 'drn', 'dimensions') changed to 32 based on size of stress_period_data
    writing package oc...
FloPy is using the following  executable to run the model: ../Exe/mf6.exe
                                   MODFLOW 6
                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
                            VERSION 6.0.3 08/09/2018

   MODFLOW 6 compiled Aug 09 2018 13:40:32 with IFORT compiler (ver. 18.0.3)

This software has been approved for release by the U.S. Geological 
Survey (USGS). Although the software has been subjected to rigorous 
review, the USGS reserves the right to update the software as needed 
pursuant to further analysis and review. No warranty, expressed or 
implied, is made by the USGS or the U.S. Government as to the 
functionality of the software and related material nor shall the 
fact of release constitute any such warranty. Furthermore, the 
software is released on condition that neither the USGS nor the U.S. 
Government shall be held liable for any damages resulting from its 
authorized or unauthorized use. Also refer to the USGS Water 
Resources Software User Rights Notice for complete use, copyright, 
and distribution information.

 Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/03/27 10:24:15

 Writing simulation list file: mfsim.lst
 Using Simulation name file: mfsim.nam
 Solving:  Stress period:     1    Time step:     1
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/03/27 10:24:15
 Elapsed run time:  0.041 Seconds

 Normal termination of simulation.

Output data representation

We import the computed heads and represent them on the model extension. A representation of the head equipotential is done by the interpolation of heads at each cell centroid.

fname = os.path.join(workspace, name + '.hds')
hdobj = flopy.utils.HeadFile(fname, precision='double')
head = hdobj.get_data()

fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
img=tri.plot(ax=ax, a=head[0, 0, :], cmap='Spectral')
fig.colorbar(img, fraction=0.02)

<matplotlib.colorbar.Colorbar at 0xa6eacf8>
output_28_1.png
from scipy.interpolate import griddata
x, y = tri.get_xcyc()[:,0],tri.get_xcyc()[:,1]
xG = np.linspace(x.min(),x.max(),100)
yG = np.linspace(y.min(),y.max(),100)
X, Y = np.meshgrid(xG, yG)
z = head[0,0,:]

fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
img=tri.plot(ax=ax, edgecolor='lightsteelblue', a=head[0, 0, :], cmap='Spectral', alpha=0.8)
fig.colorbar(img, fraction=0.02)
Ti = griddata((x, y), z, (X, Y), method='cubic')
CS = ax.contour(X,Y,Ti, colors='Cyan', linewidths=2.0)
ax.clabel(CS, inline=1, fontsize=15, fmt='%1.0f')
plt.show()
output_30_0.png


Input files

You can download the input files for this tutorial on this link.

2 Comments

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.