This example develops a groundwater model in MODFLOW6 DISV that implements a Voronoi mesh generated from the basin boundary and river network with the mf6Voronoi packages that allow high performance meshing fully coupled with geospatial data. On the Voronoi mesh, the refinement levels are defined by a minimum cell size, maximum cell size and a multiplier. The applied case covers all steps on model discretization, construction, simulation and 2D visualization.

Tutorial

Code

Part 1 : Voronoi mesh generation

import warnings
warnings.filterwarnings('ignore')

import os
import geopandas as gpd
from mf6Voronoi.geoVoronoi import createVoronoi
from mf6Voronoi.geoUtils import plotOrgDistPoints, plotCirclesPoints, plotKeyList

#Create mesh object specifying the coarse mesh and the multiplier
vorMesh = createVoronoi(meshName='regionalModel',maxRef = 250, multiplier=1)

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

#Generate point pair array
vorMesh.generateOrgDistVertices()

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

/--------Layer river discretization-------/
Progressive cell size list: [100, 200] m.

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

Time required for point generation: 52.49 seconds 


/----Generation of the voronoi mesh----/

Time required for voronoi generation: 17.16 seconds

#Export generated voronoi mesh
vorMesh.getVoronoiAsShp(outputPath='output')

/----Generation of the voronoi shapefile----/
The output folder output exists

Time required for voronoi shapefile: 1.57 seconds

# Show the resulting voronoi mesh

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

Part 2 generate disv properties

import sys
import os
from mf6Voronoi.meshProperties import meshShape

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

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

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


100%|███████████████████████████████████████████████████████████████████████████| 7121/7121 [00:00<00:00, 17982.60it/s]



Extracting cell2d data and grid index


100%|████████████████████████████████████████████████████████████████████████████| 7121/7121 [00:02<00:00, 2879.62it/s]

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 3: Model construction and simulation

import rasterio
import numpy as np
import flopy
import matplotlib.pyplot as plt

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

nlay = 5

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


AcuifInf_Bottom = 2800
zbot[0,] = mtop - 30
zbot[1,] = AcuifInf_Bottom + (0.85 * (mtop - AcuifInf_Bottom))
zbot[2,] = AcuifInf_Bottom + (0.70 * (mtop - AcuifInf_Bottom))
zbot[3,] = AcuifInf_Bottom + (0.50 * (mtop - AcuifInf_Bottom))
zbot[4,] = AcuifInf_Bottom

Create simulation and model

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

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

# create gwf model
gwf = flopy.mf6.ModflowGwf(sim, modelname=modelName, save_flows=True)

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

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

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

Kx =[4E-4,5E-6,1E-6,9E-7,5E-7]
icelltype = [1,1,0,0,0]

# node property flow
npf = flopy.mf6.ModflowGwfnpf(gwf, xt3doptions=[('xt3d')],
                              save_specific_discharge=True,
                              icelltype=icelltype, 
                              k=Kx)

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

Working with rechage, evapotranspiration

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

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)

from shapely.geometry import MultiLineString

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

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

array([462, 466, 518, 530, 541, 567, 569, 570, 578, 579], dtype=object)

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

#river plot
riv.plot(mflay=0)

[<Axes: title={'center': ' drn_0 location stress period 1 layer 1'}>]

Set the Output Control and run simulation

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

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

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 818 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/19/2024

   MODFLOW 6 compiled Dec 19 2024 21:59:25 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/02/05 11:56:22
 
 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/02/05 11:58:02
 Elapsed run time:  1 Minutes, 40.269 Seconds
 
 Normal termination of simulation.

Model output visualization

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

[(0, 0)]

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

array([3740.65615377, 3741.51000605, 3741.07284441, 3739.15569422,
       3784.00173451])

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

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


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

# Boundary conditions
shpRiv = flopy.plot.plot_shapefile('Shp/river_basin', ax=ax)

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

plt.colorbar(cellhead, shrink=0.75)

plt.show()