Applied case of groundwater modeling in the area influenced by a water reservoir. This tutorial provides the regional flow conditions that serve as boundary conditions for a refined model of the reservoir's interaction with groundwater.

Tutorial

Code

#!pip install -U mf6Voronoi

from mf6Voronoi.utils import listTemplates, copyTemplate

listTemplates()

/-------- List of available mf6Voronoi templates --------/

Nr 1: generateVoronoi
    File: p1_generateVoronoi.ipynb
    Description: Template for Voronoi mesh generation

Nr 2: modelCreation
    File: p2_modelCreation.ipynb
    Description: Template for regional model creation on steady state

Nr 3: obsCalculated
    File: p3_obsCalculated.ipynb
    Description: Template for comparison between observed and calculated

Nr 4: vtkGeneration
    File: p4_vtkGeneration.ipynb
    Description: Template for 3d geometry generation on Vtk format

Nr 5: multilayeredTransient
    File: p5_multilayeredTransient.ipynb
    Description: Template for multilayer (15 layers) and transient (6 stress periods) regional model

Nr 6: basicTransport
    File: p6_basicTransport.ipynb
    Description: Template with basic lines for transport modeling based on a flow model

copyTemplate('generateVoronoi','regWk')

copyTemplate('modelCreation','regWk')

copyTemplate('vtkGeneration','regWk')

Mesh generation

Part 1 : Voronoi mesh generation

import warnings ## Org
warnings.filterwarnings('ignore') ## Org

import os, sys ## Org
import geopandas as gpd ## Org
from mf6Voronoi.geoVoronoi import createVoronoi ## Org
from mf6Voronoi.meshProperties import meshShape ## Org
from mf6Voronoi.utils import initiateOutputFolder, getVoronoiAsShp ## Org

#Create mesh object specifying the coarse mesh and the multiplier
vorMesh = createVoronoi(meshName='regionalModel',maxRef = 800, multiplier=1.5) ## Org

#Open limit layers and refinement definition layers
vorMesh.addLimit('basin','../../hatariUtils/catchment.shp') ## Org
vorMesh.addLayer('river','../../hatariUtils/river_basin.shp',80) ## Org

#Generate point pair array
vorMesh.generateOrgDistVertices() ## Org

#Generate the point cloud and voronoi
vorMesh.createPointCloud() ## Org
vorMesh.generateVoronoi() ## Org

mf6Voronoi will have a web version in 2028

/--------Layer river discretization-------/
Progressive cell size list: [80, 200.0, 380.0, 650.0] m.

/----Sumary of points for voronoi meshing----/
Distributed points from layers: 1
Points from layer buffers: 10469
Points from max refinement areas: 1028
Points from min refinement areas: 0
Total points inside the limit: 13890
/--------------------------------------------/

Time required for point generation: 14.28 seconds 


/----Generation of the voronoi mesh----/
Something went wrong

Time required for voronoi generation: 11.51 seconds

#Uncomment the next two cells if you have strong differences on discretization or you have encounter an FORTRAN error while running MODFLOW6

vorMesh.checkVoronoiQuality(threshold=0.01)

/----Performing quality verification of voronoi mesh----/
Short side on polygon: 13890 with length = 0.00218
Short side on polygon: 13890 with length = 0.00218

vorMesh.fixVoronoiShortSides()
vorMesh.generateVoronoi()
vorMesh.checkVoronoiQuality(threshold=0.01)

/----Generation of the voronoi mesh----/
Something went wrong

Time required for voronoi generation: 12.56 seconds 


/----Performing quality verification of voronoi mesh----/
Your mesh has no edges shorter than your threshold

#Export generated voronoi mesh
initiateOutputFolder('../../regionalFiles') #<==== updated
initiateOutputFolder('../../regionalFiles/output') ## Org
getVoronoiAsShp(vorMesh.modelDis, shapePath='../../regionalFiles/output/'+vorMesh.modelDis['meshName']+'.shp') ## Org

The output folder ../../regionalFiles exists and has been cleared
The output folder ../../regionalFiles/output has been generated.

/----Generation of the voronoi shapefile----/

Time required for voronoi shapefile: 1.47 seconds

# Show the resulting voronoi mesh

#open the mesh file
mesh=gpd.read_file('../../regionalFiles/output/'+vorMesh.modelDis['meshName']+'.shp') ## Org
#plot the mesh
mesh.plot(figsize=(35,25), fc='crimson', alpha=0.3, ec='teal') ## Org

Part 2 generate disv properties

# open the mesh file
mesh=meshShape('../../regionalFiles/output/'+vorMesh.modelDis['meshName']+'.shp') ## Org

# get the list of vertices and cell2d data
gridprops=mesh.get_gridprops_disv() ## Org

Creating a unique list of vertices [[x1,y1],[x2,y2],...]


100%|██████████| 13892/13892 [00:00<00:00, 44092.10it/s]



Extracting cell2d data and grid index


100%|██████████| 13892/13892 [00:01<00:00, 8038.46it/s]

#create folder
initiateOutputFolder('../../regionalFiles/json') ## Org

#export disv
mesh.save_properties('../../regionalFiles/json/disvDict.json') ## Org

The output folder ../../regionalFiles/json has been generated.

Model creation on steady state

Part 2a: generate disv properties

import sys, json, os ## Org
import rasterio, flopy ## Org
import numpy as np ## Org
import matplotlib.pyplot as plt ## Org
import geopandas as gpd ##Org
from mf6Voronoi.meshProperties import meshShape ## Org
from shapely.geometry import MultiLineString ## Org
from mf6Voronoi.tools.graphs2d import generateRasterFromArray, generateContoursFromRaster

# open the json file
with open('../../regionalFiles/json/disvDict.json') as file: ## Org
    gridProps = json.load(file) ## Org

cell2d = gridProps['cell2d']           #cellid, cell centroid xy, vertex number and vertex id list
vertices = gridProps['vertices']       #vertex id and xy coordinates
ncpl = gridProps['ncpl']               #number of cells per layer
nvert = gridProps['nvert']             #number of verts
centroids=gridProps['centroids']       #cell centroids xy

Part 2b: Model construction and simulation

#Extract dem values for each centroid of the voronois
src = rasterio.open('../../rst/n33w111_wgs84_int32_50m.tif')  ## Org
elevation=[x for x in src.sample(centroids)] ## Org

nlay = 10   ## Org

mtop=np.array([elev[0] for i,elev in enumerate(elevation)]) ## Org
zbot=np.zeros((nlay,ncpl)) ## Org

AcuifInf_Bottom = 700 ## Org
zbot[0,] = AcuifInf_Bottom + (0.95 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[1,] = AcuifInf_Bottom + (0.90 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[2,] = AcuifInf_Bottom + (0.85 * (mtop - AcuifInf_Bottom)) ## <==== updated 85%
zbot[3,] = AcuifInf_Bottom + (0.78 * (mtop - AcuifInf_Bottom)) ## <==== updated 
zbot[4,] = AcuifInf_Bottom + (0.71 * (mtop - AcuifInf_Bottom)) ## <==== updated 
zbot[5,] = AcuifInf_Bottom + (0.64 * (mtop - AcuifInf_Bottom)) ## <==== updated 
zbot[6,] = AcuifInf_Bottom + (0.57 * (mtop - AcuifInf_Bottom)) ## <==== updated 
zbot[7,] = AcuifInf_Bottom + (0.50 * (mtop - AcuifInf_Bottom)) ## <==== updated 50%
zbot[8,] = AcuifInf_Bottom + (0.25 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[9,] = AcuifInf_Bottom ## <==== updated

Create simulation and model

# create simulation
simName = 'mf6Sim' ## Org
modelName = 'mf6Model' ## Org
modelWs = '../../regionalFiles/modelFiles' ## Org
sim = flopy.mf6.MFSimulation(sim_name=modelName, version='mf6', ## Org
                             exe_name='../../bin/mf6.exe', ## Org
                             sim_ws=modelWs) ## Org

# create tdis package
tdis_rc = [(1000.0, 1, 1.0)] ## Org
tdis = flopy.mf6.ModflowTdis(sim, pname='tdis', time_units='SECONDS', ## Org
                             perioddata=tdis_rc) ## Org

# create gwf model
gwf = flopy.mf6.ModflowGwf(sim, ## Org
                           modelname=modelName, ## Org
                           save_flows=True, ## Org
                           newtonoptions="NEWTON UNDER_RELAXATION") ## Org

# create iterative model solution and register the gwf model with it
ims = flopy.mf6.ModflowIms(sim, ## Org
                           complexity='COMPLEX', ## Org
                           outer_maximum=50, ## Org
                           inner_maximum=30, ## Org
                           linear_acceleration='BICGSTAB') ## Org
sim.register_ims_package(ims,[modelName]) ## Org

# disv
disv = flopy.mf6.ModflowGwfdisv(gwf, nlay=nlay, ncpl=ncpl, ## Org
                                top=mtop, botm=zbot, ## Org
                                nvert=nvert, vertices=vertices, ## Org
                                cell2d=cell2d) ## Org

disv.top.plot(figsize=(12,8), alpha=0.8) ## Org

crossSection = gpd.read_file('../../shp/crossSectionRegional.shp') ## Org
sectionLine =list(crossSection.iloc[0].geometry.coords) ## Org

fig, ax = plt.subplots(figsize=(12,8)) ## Org
modelxsect = flopy.plot.PlotCrossSection(model=gwf, line={'Line': sectionLine}) ## Org
linecollection = modelxsect.plot_grid(lw=0.5) ## Org
ax.grid() ## Org

# initial conditions
ic = flopy.mf6.ModflowGwfic(gwf, strt=np.stack([mtop for i in range(nlay)])) ## Org

Kx =[4E-4, 5E-5, 3E-6, 3E-6, 2.5E-6, 2.5E-6, 2.5E-6, 1E-6, 9E-7, 5E-7] ## Org
icelltype = [1,1,1,1,1,1,0,0,0,0] ## Org

# node property flow
npf = flopy.mf6.ModflowGwfnpf(gwf, ## Org
                              save_specific_discharge=True, ## Org
                              icelltype=icelltype, ## Org
                              k=Kx) ## Org

# define storage and transient stress periods
sto = flopy.mf6.ModflowGwfsto(gwf, ## Org
                              iconvert=1, ## Org
                              steady_state={ ## Org
                                0:True, ## Org
                              } ## Org
                              ) ## Org

Working with rechage, evapotranspiration

rchr = 0.15/365/86400 ## Org
rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=rchr) ## Org
evtr = 1.2/365/86400 ## Org
evt = flopy.mf6.ModflowGwfevta(gwf,ievt=1,surface=mtop,rate=evtr,depth=1.0) ## Org

Definition of the intersect object

For the manipulation of spatial data to determine hydraulic parameters or boundary conditions

# Define intersection object
interIx = flopy.utils.gridintersect.GridIntersect(gwf.modelgrid) ## Org

#open the river shapefile
rivers =gpd.read_file('../../hatariUtils/river_basin.shp') ## Org
list_rivers=[] ## Org
for i in range(rivers.shape[0]): ## Org
    list_rivers.append(rivers['geometry'].loc[i]) ## Org
    
riverMls = MultiLineString(lines=list_rivers) ## Org

#intersec rivers with our grid
riverCells=interIx.intersect(riverMls).cellids ## Org
riverCells[:10] ## Org

array([110, 131, 139, 151, 155, 158, 166, 170, 185, 186], dtype=object)

#river package
riverSpd = {} ## Org
riverSpd[0] = [] ## Org
for cell in riverCells: ## Org
    riverSpd[0].append([(0,cell),mtop[cell],0.01]) ## Org
riv = flopy.mf6.ModflowGwfdrn(gwf, stress_period_data=riverSpd) ## Org

#river plot
riv.plot(mflay=0) ## Org

Set the Output Control and run simulation

#oc
head_filerecord = f"{gwf.name}.hds" ## Org
budget_filerecord = f"{gwf.name}.cbc" ## Org
oc = flopy.mf6.ModflowGwfoc(gwf, ## Org
                            head_filerecord=head_filerecord, ## Org
                            budget_filerecord = budget_filerecord, ## Org
                            saverecord=[("HEAD", "LAST"),("BUDGET","LAST")]) ## Org

# Run the simulation
sim.write_simulation() ## Org
success, buff = sim.run_simulation() ## Org

writing simulation...
  writing simulation name file...
  writing simulation tdis package...
  writing solution package ims_-1...
  writing model mf6Model...
    writing model name file...
    writing package disv...
    writing package ic...
    writing package npf...
    writing package sto...
    writing package rcha_0...
    writing package evta_0...
    writing package drn_0...
INFORMATION: maxbound in ('gwf6', 'drn', 'dimensions') changed to 2413 based on size of stress_period_data
    writing package oc...
FloPy is using the following executable to run the model: ..\..\bin\mf6.exe
                                   MODFLOW 6
                U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
                            VERSION 6.6.0 12/20/2024

   MODFLOW 6 compiled Dec 31 2024 17:10:16 with Intel(R) Fortran Intel(R) 64
   Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
                             Build 20220726_000000

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.

 
 MODFLOW runs in SEQUENTIAL mode
 
 Run start date and time (yyyy/mm/dd hh:mm:ss): 2025/09/18 16:07:33
 
 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): 2025/09/18 16:07:41
 Elapsed run time:  8.320 Seconds
 
 Normal termination of simulation.

Model output visualization

headObj = gwf.output.head() ## Org
headObj.get_kstpkper() ## Org

[(0, 0)]

heads = headObj.get_data() ## Org
heads[2,0,:5] ## Org

array([1327.73127261, 1313.44523133, 1325.54397454, 1134.81998389,
       1284.78111455])

# Plot the heads for a defined layer and boundary conditions
fig = plt.figure(figsize=(12,8)) ## Org
ax = fig.add_subplot(1, 1, 1, aspect='equal') ## Org
modelmap = flopy.plot.PlotMapView(model=gwf) ## Org

####
levels = np.linspace(heads[heads>-1e+30].min(),heads[heads>-1e+30].max(),num=50) ## Org
contour = modelmap.contour_array(heads[3],ax=ax,levels=levels,cmap='PuBu') ## Org
ax.clabel(contour) ## Org


quadmesh = modelmap.plot_bc('DRN') ## Org
cellhead = modelmap.plot_array(heads[3],ax=ax, cmap='Blues', alpha=0.8) ## Org

linecollection = modelmap.plot_grid(linewidth=0.3, alpha=0.5, color='cyan', ax=ax) ## Org

plt.colorbar(cellhead, shrink=0.75) ## Org

plt.show() ## Org

waterTable = flopy.utils.postprocessing.get_water_table(heads)

generateRasterFromArray(gwf, 
                        waterTable,
                        meshLayer='../../regionalFiles/output/regionalModel.shp', 
                        rasterRes=10, 
                        epsg=32612, 
                        outputPath='../../regionalFiles/output/waterTable.tif',
                        limitLayer='../../hatariUtils/catchment.shp')

WARNING: No head on vextex was found, or something went wrong with your mesh
WARNING: No head on vextex was found, or something went wrong with your mesh
Raster X Dim: 40700.00, Raster Y Dim: 37850.00
Number of cols:  4071, Number of rows: 3786

Input data

You can download the input files from this link:

owncloud.hatarilabs.com/s/KPOgDALbD3LCLib

Password: Hatarilabs

Blog

Regional groundwater modeling near a water reservoir using Modflow6 and mf6Voronoi - Tutorial

Tutorial

Code

Mesh generation

Part 1 : Voronoi mesh generation

Follow us:

Part 2 generate disv properties

Model creation on steady state

Part 2a: generate disv properties

Part 2b: Model construction and simulation

Create simulation and model

Working with rechage, evapotranspiration

Definition of the intersect object

Set the Output Control and run simulation

Model output visualization

Input data

Hatari Labs

Our courses:

Blog

Tutorial

Code

Mesh generation

Part 1 : Voronoi mesh generation

Follow us:

Part 2 generate disv properties

Model creation on steady state

Part 2a: generate disv properties

Part 2b: Model construction and simulation

Create simulation and model

Working with rechage, evapotranspiration

Definition of the intersect object

Set the Output Control and run simulation

Model output visualization

Input data

Hatari Labs

Our courses:

Follow Us

Join us on social media for the latest webinars, tutorials and courses.