Contaminant plume modeling of a mining waster rock dump with MODFLOW 6 and mf6Voronoi - Tutorial
/
This tutorial shows an applied case that covers the discretization of a model using Voronoi grids, taking into consideration the water network, the predictive flow area, and the model boundaries. Based on this mesh, a flow and transport model under steady-state conditions is constructed to simulate the contaminant plume over a period of 50 years, evaluated every 5 years. The model includes piezometers and the development of concentrations over time is evaluated using the Pandas library with Matplotlib. Finally, the entire model geometry, boundary conditions, hydraulic loads, and concentrations are visualized in 3D with Paraview.
Tutorial
Code
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='wasteDump',maxRef = 200, multiplier=1.5, overlapping=False) ## Org
#Open limit layers and refinement definition layers
vorMesh.addLimit('limit','../shp/wasteDump/modelLimit.shp') ## <==== updated
vorMesh.addLayer('refinement','../shp/wasteDump/modelRefinement.shp',25) ## <==== updated
vorMesh.addLayer('river','../hatariUtils/river_network.shp',50) ## <==== updated
vorMesh.addLayer('regFlow','../shp/wasteDump/regionalFlow.shp',100) ## <==== updated#Generate point pair array
vorMesh.generateOrgDistVertices() ## Org
#Generate the point cloud and voronoi
vorMesh.createPointCloud() ## Org
vorMesh.generateVoronoi() ## Org
build faster, analyze more
Follow us: |
 |
 |
 |
 |
 |
 |
/--------Layer refinement discretization-------/
Progressive cell size list: [25, 62.5, 118.75] m.
/--------Layer river discretization-------/
Progressive cell size list: [50, 125.0] m.
/--------Layer regFlow discretization-------/
Progressive cell size list: [100] m.
/----Sumary of points for voronoi meshing----/
Distributed points from layers: 3
Points from layer buffers: 2356
Points from max refinement areas: 785
Points from min refinement areas: 17979
Total points inside the limit: 22158
/--------------------------------------------/
Time required for point generation: 2.58 seconds
/----Generation of the voronoi mesh----/
Time required for voronoi generation: 2.49 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)#vorMesh.fixVoronoiShortSides()
#vorMesh.generateVoronoi()
#vorMesh.checkVoronoiQuality(threshold=0.01)#Export generated voronoi mesh
initiateOutputFolder('../output') ## Org
getVoronoiAsShp(vorMesh.modelDis, shapePath='../output/'+vorMesh.modelDis['meshName']+'.shp') ## OrgThe output folder ../output exists and has been cleared
/----Generation of the voronoi shapefile----/
Time required for voronoi shapefile: 4.95 seconds# Show the resulting voronoi mesh
#open the mesh file
mesh=gpd.read_file('../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('../output/'+vorMesh.modelDis['meshName']+'.shp') ## Org# get the list of vertices and cell2d data
gridprops=mesh.get_gridprops_disv() ## OrgCreating a unique list of vertices [[x1,y1],[x2,y2],...]
100%|█████████████████████████████████████████████████████████████████████████| 22158/22158 [00:01<00:00, 17508.92it/s]
Extracting cell2d data and grid index
100%|██████████████████████████████████████████████████████████████████████████| 22158/22158 [00:09<00:00, 2270.88it/s]#create folder
initiateOutputFolder('../json') ## Org
#export disv
mesh.save_properties('../json/disvDict.json') ## OrgThe output folder ../json exists and has been clearedPart 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.cellWork import getCellFromGeom, getLayCellElevTupleFromRaster ## <==== updated
from mf6Voronoi.tools.graphs2d import generateRasterFromArray ## <==== updatedC:\Users\saulm\anaconda3\Lib\site-packages\geopandas\_compat.py:7: DeprecationWarning: The 'shapely.geos' module is deprecated, and will be removed in a future version. All attributes of 'shapely.geos' are available directly from the top-level 'shapely' namespace (since shapely 2.0.0).
import shapely.geos# open the json file
with open('../json/disvDict.json') as file: ## Org
gridProps = json.load(file) ## Orgcell2d = 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 xyPart 2b: Model construction and simulation
#Extract dem values for each centroid of the voronois
src = rasterio.open('../rst/clipDem100m.tif') ## Org
elevation=[x for x in src.sample(centroids)] ## Orgnlay = 11 ## <==== updated
mtop=np.array([elev[0] for i,elev in enumerate(elevation)]) ## Org
zbot=np.zeros((nlay,ncpl)) ## Org
AcuifInf_Bottom = 170 ## <==== updated
zbot[0,] = AcuifInf_Bottom + (0.96 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[1,] = AcuifInf_Bottom + (0.92 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[2,] = AcuifInf_Bottom + (0.88 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[3,] = AcuifInf_Bottom + (0.84 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[4,] = AcuifInf_Bottom + (0.80 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[5,] = AcuifInf_Bottom + (0.76 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[6,] = AcuifInf_Bottom + (0.72 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[7,] = AcuifInf_Bottom + (0.68 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[8,] = AcuifInf_Bottom + (0.64 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[9,] = AcuifInf_Bottom + (0.50 * (mtop - AcuifInf_Bottom)) ## <==== updated
zbot[10,] = AcuifInf_Bottom ## OrgCreate simulation and model
# create simulation
simName = 'mf6Sim' ## Org
modelName = 'mf6Model' ## Org
modelWs = '../modelFiles' ## Org
sim = flopy.mf6.MFSimulation(sim_name=modelName, version='mf6', ## Org
exe_name='../bin/mf6.exe', ## Org
continue_=True,
sim_ws=modelWs) ## Org# create tdis package
tdis_rc = [(50*86400*365.0, 10, 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# initial conditions
ic = flopy.mf6.ModflowGwfic(gwf, strt=np.stack([mtop for i in range(nlay)])) ## OrgKx =[4E-4,5E-6,5E-6,5E-6,5E-6,1E-6,1E-6,1E-6,1E-6,9E-7,5E-7] ## <===== updated
icelltype = [1,1,1,1,0,0,0,0,0,0,0] ## <===== updated
# node property flow
npf = flopy.mf6.ModflowGwfnpf(gwf, ## Org
save_specific_discharge=True, ## Org
save_flows=True, ## <==== inserted
save_saturation=True, ## <==== inserted
icelltype=icelltype, ## Org
k=Kx) ## Org## <==== inserted
crossSection = gpd.read_file('../shp/wasteDump/crossSection.shp') ## <==== inserted
sectionLine =list(crossSection.iloc[0].geometry.coords) ## <==== inserted
fig, ax = plt.subplots(figsize=(12,8)) ## <==== inserted
modelxsect = flopy.plot.PlotCrossSection(model=gwf, line={'Line': sectionLine}) ## <==== inserted
linecollection = modelxsect.plot_grid(lw=0.5) ## <==== inserted
modelxsect.plot_array(np.log(npf.k.array), alpha=0.5) ## <====== inserted
ax.grid() ## <==== inserted
# define storage and transient stress periods
sto = flopy.mf6.ModflowGwfsto(gwf, ## Org
iconvert=1, ## Org
steady_state={ ## Org
0:True, ## Org
} ## Org
) ## OrgWorking with rechage, evapotranspiration
interIx = flopy.utils.gridintersect.GridIntersect(gwf.modelgrid)
rchr = (0.1486*mtop + 45.72)*0.5/1000/365/86400 ## <===== update
#tailings initial concentration
cellList = getCellFromGeom(gwf,
interIx,
'../shp/wasteDump/wasteDump.shp')
for cell in cellList: ## <===== inserted
rchr[cell] = 5e-9 ## <===== inserted
rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=rchr) ## Org
rch.plot()
evtr =(-0.0633*mtop + 1336.3)/1000/365/24/86400 ## <===== update
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) ## Org#open the river shapefile
rivers =gpd.read_file('../hatariUtils/river_network.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] ## Orgarray([10, 11, 24, 25, 133, 164, 170, 177, 309, 312], 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
#river package # <===== Inserted
layCellTupleList, cellElevList = getLayCellElevTupleFromRaster(gwf,
interIx,
'../rst/waterTable.tif',
'../shp/wasteDump/regionalFlow.shp') ## # <===== Inserted
regionalSpd = {} ## # <===== Inserted
regionalSpd[0] = [] ## # <===== Inserted
for index, layCellTuple in enumerate(layCellTupleList): ## Org
regionalSpd[0].append([layCellTuple,cellElevList[index],0.01]) # <===== InsertedThe cell 132 has a elevation of 478.76 outside the model vertical domain
The cell 4711 has a elevation of 546.37 outside the model vertical domain
The cell 4714 has a elevation of 547.37 outside the model vertical domain
The cell 4999 has a elevation of 546.99 outside the model vertical domainghb = flopy.mf6.ModflowGwfghb(gwf, stress_period_data=regionalSpd)#regional flow plot
ghb.plot(mflay=2, kper=0) ## <===== modified
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", "ALL"),("BUDGET","ALL")]) ## Org# Run the simulation
sim.write_simulation() ## Org
success, buff = sim.run_simulation() ## Orgwriting 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 326 based on size of stress_period_data
writing package ghb_0...
INFORMATION: maxbound in ('gwf6', 'ghb', 'dimensions') changed to 304 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/08/13 18:47:03
Writing simulation list file: mfsim.lst
Using Simulation name file: mfsim.nam
Solving: Stress period: 1 Time step: 1
Solving: Stress period: 1 Time step: 2
Solving: Stress period: 1 Time step: 3
Solving: Stress period: 1 Time step: 4
Solving: Stress period: 1 Time step: 5
Solving: Stress period: 1 Time step: 6
Solving: Stress period: 1 Time step: 7
Solving: Stress period: 1 Time step: 8
Solving: Stress period: 1 Time step: 9
Solving: Stress period: 1 Time step: 10
Run end date and time (yyyy/mm/dd hh:mm:ss): 2025/08/13 18:49:14
Elapsed run time: 2 Minutes, 10.872 Seconds
WARNING REPORT:
1. Simulation convergence failure occurred 1 time(s).
Normal termination of simulation.Model output visualization
headObj = gwf.output.head() ## Org
headObj.get_kstpkper() ## Org[(0, 0),
(1, 0),
(2, 0),
(3, 0),
(4, 0),
(5, 0),
(6, 0),
(7, 0),
(8, 0),
(9, 0)]heads = headObj.get_data() ## Org
heads[2,0,:5] ## Orgarray([483.96540049, 483.88786761, 483.93685517, 484.14018628,
479.36809215])# 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
#plot cross section
crossSection = gpd.read_file('../shp/wasteDump/crossSection.shp') ## <==== inserted
sectionLine =list(crossSection.iloc[0].geometry.coords) ## <==== inserted
waterTable = flopy.utils.postprocessing.get_water_table(heads)
fig = plt.figure(figsize=(18, 5))
ax = fig.add_subplot(1, 1, 1)
ax.set_title("plot_array()")
xsect = flopy.plot.PlotCrossSection(model=gwf, line={"line": sectionLine})
patch_collection = xsect.plot_array(heads, alpha=0.5)
xsect.plot_surface(waterTable)
line_collection = xsect.plot_grid()
cb = plt.colorbar(patch_collection, shrink=0.75)
waterTable = flopy.utils.postprocessing.get_water_table(heads)
generateRasterFromArray(gwf,
waterTable,
rasterRes=10,
epsg=32614,
outputPath='../output/waterTable.tif',
limitLayer='../shp/wasteDump/modelLimit.shp')Raster X Dim: 6600.00, Raster Y Dim: 7900.00
Number of cols: 661, Number of rows: 791
[483.9779653 483.90023243 483.94941999 ... 579.09532336 579.4460603
580.50993036]#Basic lines for transport modeling
import flopy ## Org
import json, os ## Org
import numpy as np ## Org
import geopandas as gpd ## Org
import matplotlib.pyplot as plt ## Org
from mf6Voronoi.tools.cellWork import getLayCellElevTupleFromRaster, getLayCellElevTupleFromObs ## OrgC:\Users\saulm\anaconda3\Lib\site-packages\geopandas\_compat.py:7: DeprecationWarning: The 'shapely.geos' module is deprecated, and will be removed in a future version. All attributes of 'shapely.geos' are available directly from the top-level 'shapely' namespace (since shapely 2.0.0).
import shapely.geos# load simulation
simName = 'mf6Sim' ## Org
modelName = 'mf6Model' ## Org
modelWs = os.path.abspath('../modelFiles') ## Org
sim = flopy.mf6.MFSimulation.load(sim_name=modelName, version='mf6', ## Org
exe_name='../bin/mf6.exe', ## Org
sim_ws=modelWs) ## Org
transWs = os.path.abspath('../transFiles') ## Org
#change working directory
sim.set_sim_path(transWs) ## Org
sim.write_simulation(silent=True) ## Orgloading simulation...
loading simulation name file...
loading tdis package...
loading model gwf6...
loading package disv...
loading package ic...
loading package npf...
loading package sto...
loading package rch...
loading package evt...
loading package drn...
loading package ghb...
loading package oc...
loading solution package mf6model...#list model names
sim.model_names ## Org['mf6model']#select the flow model
gwf = sim.get_model('mf6model') ## Org# open the json file
with open('../json/disvDict.json') as file: ## Org
gridProps = json.load(file) ## Org
cell2d = gridProps['cell2d'] #cellid, cell centroid xy, vertex number and vertex id list ## Org
vertices = gridProps['vertices'] #vertex id and xy coordinates ## Org
ncpl = gridProps['ncpl'] #number of cells per layer ## Org
nvert = gridProps['nvert'] #number of verts ## Org
centroids=gridProps['centroids'] ## Org#define the transport model ## Org
gwt = flopy.mf6.ModflowGwt(sim, ## Org
modelname='gwtModel', ## Org
save_flows=True) ## Org#register solver for transport model
imsGwt = flopy.mf6.ModflowIms(sim, print_option='SUMMARY', ## Org
outer_dvclose=1e-4, ## Org
inner_dvclose=1e-4, ## Org
linear_acceleration='BICGSTAB') ## Org
sim.register_ims_package(imsGwt,[gwt.name]) ## Org# apply discretization to transport model
disv = flopy.mf6.ModflowGwtdisv(gwt, ## Org
nlay=gwf.modelgrid.nlay, ## Org
ncpl=ncpl, ## Org
top=gwf.modelgrid.top, ## Org
botm=gwf.modelgrid.botm, ## Org
nvert=nvert, ## Org
vertices=vertices, ## Org
cell2d=cell2d) ## Org#define starting concentrations
strtConc = np.zeros((gwf.modelgrid.nlay, ncpl), dtype=np.float32) ## Org
interIx = flopy.utils.gridintersect.GridIntersect(gwf.modelgrid) ## Org
#waste dump recharge
layCellTupleList, cellElevList = getLayCellElevTupleFromRaster(gwf, ## <===== inserted
interIx, ## <===== inserted
'../output/waterTable.tif', ## <===== inserted
'../shp/wasteDump/wasteDump.shp') ## <===== inserted
for layCell in layCellTupleList: ## <===== inserted
strtConc[layCell[0],layCell[1]] = 1200 ## <===== inserted
gwtIc = flopy.mf6.ModflowGwtic(gwt, strt=strtConc) ## OrgThe cell 11204 has a elevation of 558.07 outside the model vertical domain
The cell 11206 has a elevation of 558.26 outside the model vertical domain
The cell 11207 has a elevation of 558.34 outside the model vertical domain
The cell 11209 has a elevation of 558.42 outside the model vertical domain
The cell 11210 has a elevation of 558.45 outside the model vertical domainfig = plt.figure(figsize=(12, 12)) ## Org
ax = fig.add_subplot(1, 1, 1, aspect = 'equal') ## Org
mapview = flopy.plot.PlotMapView(model=gwf,layer = 2) ## Org
plot_array = mapview.plot_array(strtConc,masked_values=[-1e+30], cmap=plt.cm.summer) ## Org
plt.colorbar(plot_array, shrink=0.75,orientation='horizontal', pad=0.08, aspect=50) ## Org
# set advection, dispersion
adv = flopy.mf6.ModflowGwtadv(gwt, scheme='UPSTREAM') ## Org
dsp = flopy.mf6.ModflowGwtdsp(gwt, alh=1.0, ath1=0.1) ## Org
#define mobile storage and transfer
porosity = 0.05 ## Org
sto = flopy.mf6.ModflowGwtmst(gwt, porosity=porosity) ## Org#define sink and source package
sourcerecarray = [()] ## Org
ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray) ## OrgcncSpd = {} ## Org
cncSpd[0] = [] ## Org
for layCell in layCellTupleList: ## <===== inserted
cncSpd[0].append([layCell[0],layCell[1],1200]) ## <===== inserted
cnc = flopy.mf6.ModflowGwtcnc(gwt,stress_period_data=cncSpd) ## Org
cnc.plot(mflay=2) ## Org
obsList = []
nameList, obsLayCellList = getLayCellElevTupleFromObs(gwf, ## Org
interIx, ## Org
'../shp/wasteDump/piezometers.shp', ## Org
'Name', ## Org
'screenElev') ## Org
for obsName, obsLayCell in zip(nameList, obsLayCellList): ## Org
obsList.append((obsName,'concentration',obsLayCell[0]+1,obsLayCell[1]+1)) ## Org
obs = flopy.mf6.ModflowUtlobs( ## Org
gwt,
filename=gwt.name+'.obs', ## Org
digits=10, ## Org
print_input=True, ## Org
continuous={gwt.name+'.obs.csv': obsList} ## Org
)Working for cell 14012
Well screen elev of 452.00 found at layer 6
Working for cell 7673
Well screen elev of 528.00 found at layer 5
Working for cell 5061
Well screen elev of 435.00 found at layer 6
Working for cell 2682
Well screen elev of 628.00 found at layer 4
Working for cell 7804
Well screen elev of 478.00 found at layer 6#define output control
oc_gwt = flopy.mf6.ModflowGwtoc(gwt, ## Org
budget_filerecord='%s.cbc'%gwt.name, ## Org
concentration_filerecord='%s.ucn'%gwt.name, ## Org
saverecord=[('CONCENTRATION', 'LAST'), ('BUDGET', 'LAST')], ## Org
printrecord=[('CONCENTRATION', 'LAST'), ('BUDGET', 'LAST')], ## Org
)pd = [ ## Org
("GWFHEAD", "../modelFiles/mf6Model.hds"), ## <===== updated
("GWFBUDGET", "../modelFiles/mf6Model.cbc"), ## <===== updated
]
fmi = flopy.mf6.ModflowGwtfmi(gwt, ## Org
flow_imbalance_correction=True, ## Org
packagedata=pd) ## Orgsim.write_simulation() ## Org
success, buff = sim.run_simulation() ## Orgwriting simulation...
writing simulation name file...
writing simulation tdis package...
writing solution package ims_0...
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...
writing package ghb_0...
writing package oc...
writing model gwtModel...
writing model name file...
writing package disv...
writing package ic...
writing package adv...
writing package dsp...
writing package mst...
writing package ssm...
writing package cnc_0...
INFORMATION: maxbound in ('gwt6', 'cnc', 'dimensions') changed to 6140 based on size of stress_period_data
writing package obs_0...
writing package oc...
writing package fmi...
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/08/13 18:59:18
Writing simulation list file: mfsim.lst
Using Simulation name file: mfsim.nam
Solving: Stress period: 1 Time step: 1
Solving: Stress period: 1 Time step: 2
Solving: Stress period: 1 Time step: 3
Solving: Stress period: 1 Time step: 4
Solving: Stress period: 1 Time step: 5
Solving: Stress period: 1 Time step: 6
Solving: Stress period: 1 Time step: 7
Solving: Stress period: 1 Time step: 8
Solving: Stress period: 1 Time step: 9
Solving: Stress period: 1 Time step: 10
Run end date and time (yyyy/mm/dd hh:mm:ss): 2025/08/13 19:03:32
Elapsed run time: 4 Minutes, 14.194 Seconds
WARNING REPORT:
1. Simulation convergence failure occurred 1 time(s).
Normal termination of simulation.concObj = gwt.output.concentration() ## Org
concObj.get_times() ## Org[1576800000.0]conc = concObj.get_data() ## Org
conc.shape ## Org(11, 1, 22158)transAoi = gpd.read_file('../shp/wasteDump/modelRefinement.shp') ## Org
xMin, yMin, xMax, yMax = transAoi.bounds.iloc[0].values ## Orgmflay = 2 ## Org
# 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(0,conc.max()/2,num=10) ## Org
quadmesh = modelmap.plot_bc('DRN', color='crimson') ## Org
contour = modelmap.contour_array(conc[mflay],ax=ax,levels=levels,cmap='summer') ## Org
ax.clabel(contour) ## Org
linecollection = modelmap.plot_grid(linewidth=0.1, alpha=0.8, color='cyan', ax=ax) ## Org
cellConc = modelmap.plot_array(conc[mflay],ax=ax,cmap='Blues') ## Org
quadmesh = modelmap.plot_bc('DRN', color='slategrey') ## Org
#dump1 = modelmap.plot_shapefile('../shp/wasteDump1.shp')
#piezo = modelmap.plot_shapefile('../shp/piezometers2.shp', radius=10)
ax.set_xlim(xMin-500,xMax+500)
ax.set_ylim(yMin-500,yMax+500)
plt.colorbar(cellConc, shrink=0.25)
plt.show()
#Vtk generation
import flopy ## Org
from mf6Voronoi.tools.vtkGen import Mf6VtkGenerator ## Org
from mf6Voronoi.utils import initiateOutputFolder ## Org
from mf6Voronoi.tools.graphs2d import generateRasterFromArray, generateContoursFromRaster ## <=== insertedC:\Users\saulm\anaconda3\Lib\site-packages\geopandas\_compat.py:7: DeprecationWarning: The 'shapely.geos' module is deprecated, and will be removed in a future version. All attributes of 'shapely.geos' are available directly from the top-level 'shapely' namespace (since shapely 2.0.0).
import shapely.geos# load simulation
simName = 'mf6Sim' ## Org
modelName = 'mf6Model' ## Org
modelWs = '../transFiles' ## Org
sim = flopy.mf6.MFSimulation.load(sim_name=modelName, version='mf6', ## Org
exe_name='../bin/mf6.exe', ## Org
sim_ws=modelWs) ## Orgloading simulation...
loading simulation name file...
loading tdis package...
loading model gwf6...
loading package disv...
loading package ic...
loading package npf...
loading package sto...
loading package rch...
loading package evt...
loading package drn...
loading package ghb...
loading package oc...
loading model gwt6...
loading package disv...
loading package ic...
loading package adv...
loading package dsp...
loading package mst...
loading package ssm...
loading package cnc...
loading package obs...
loading package oc...
loading package fmi...
loading solution package mf6model...vtkDir = '../vtk' ## Org
initiateOutputFolder(vtkDir) ## Org
mf6Vtk = Mf6VtkGenerator(sim, vtkDir) ## OrgThe output folder ../vtk exists and has been clearedbuild faster, analyze more
Follow us: |
|
|
|
|
|
|
/---------------------------------------/
The Vtk generator engine has been started
/---------------------------------------/#list models on the simulation
mf6Vtk.listModels() ## OrgModels in simulation: ['mf6model', 'gwtmodel']mf6Vtk.loadModel(modelName) ## OrgPackage list: ['DISV', 'IC', 'NPF', 'STO', 'RCHA_0', 'EVTA_0', 'DRN_0', 'GHB_0', 'OC']#show output data
headObj = mf6Vtk.gwf.output.head() ## Org
headObj.get_kstpkper() ## Org[(0, 0),
(1, 0),
(2, 0),
(3, 0),
(4, 0),
(5, 0),
(6, 0),
(7, 0),
(8, 0),
(9, 0)]#generate model geometry as vtk and parameter array
mf6Vtk.generateGeometryArrays() ## Org#generate parameter vtk
mf6Vtk.generateParamVtk() ## OrgParameter Vtk Generated#generate bc and obs vtk
mf6Vtk.generateBcObsVtk(nper=0) ## Org/--------RCHA_0 vtk generation-------/
Working for RCHA_0 package, creating the datasets: dict_keys(['irch', 'recharge', 'aux'])
Vtk file took 0.1880 seconds to be generated.
/--------RCHA_0 vtk generated-------/
/--------EVTA_0 vtk generation-------/
Working for EVTA_0 package, creating the datasets: dict_keys(['ievt', 'surface', 'rate', 'depth', 'aux'])
Vtk file took 0.2043 seconds to be generated.
/--------EVTA_0 vtk generated-------/
/--------DRN_0 vtk generation-------/
Working for DRN_0 package, creating the datasets: ('elev', 'cond')
Vtk file took 2.1800 seconds to be generated.
/--------DRN_0 vtk generated-------/
/--------GHB_0 vtk generation-------/
Working for GHB_0 package, creating the datasets: ('bhead', 'cond')
Vtk file took 2.1950 seconds to be generated.
/--------GHB_0 vtk generated-------/mf6Vtk.generateHeadVtk(nper=0, crop=True) ## Orgmf6Vtk.generateWaterTableVtk(nper=0) ## Orggwt = sim.get_model('gwtmodel')
concObj = gwt.output.concentration()
concObj.get_times()[1576800000.0]conc = concObj.get_data(totim=1576800000.0)
conc[4][0]array([4.55627840e-03, 6.32667844e-03, 5.88916975e-03, ...,
4.33387581e-23, 2.19599450e-21, 2.17184460e-23])mf6Vtk.generateArrayVtk(conc, 'conc50y', nper=0,nstp=9, crop=True)#generate contours
concPath = '../output/conc50y_lay4.tif'
generateRasterFromArray(mf6Vtk.gwf, conc[3][0],
rasterRes=10, epsg=32614,
outputPath=concPath,
limitLayer='../shp/wasteDump/modelLimit.shp')Raster X Dim: 6600.00, Raster Y Dim: 7900.00
Number of cols: 661, Number of rows: 791
[0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 5.49645745e-26
0.00000000e+00 6.82984931e-26]generateContoursFromRaster('../output/conc50y_lay4.tif',5,
'../output/conc50y_lay4.shp')import pandas as pd
import matplotlib.pyplot as pltgwtDf = pd.read_csv('../transFiles/gwtModel.obs.csv', na_values=3e30)
gwtDf.head()| time | OBS1 | OBS2 | OBS3 | OBS4 | OBS5 | |
|---|---|---|---|---|---|---|
| 0 | 157680000.0 | 107.286321 | 60.666416 | 1.349865 | 0.664076 | 0.227713 |
| 1 | 315360000.0 | 339.175837 | 246.724857 | 9.548925 | 21.047711 | 45.011523 |
| 2 | 473040000.0 | 558.399591 | 452.015263 | 31.501198 | 71.914114 | 143.841551 |
| 3 | 630720000.0 | 772.058932 | 626.672568 | 71.062098 | 138.455888 | 265.740493 |
| 4 | 788400000.0 | 902.509562 | 757.446171 | 126.441138 | 207.338829 | 391.997757 |
obsDf = pd.DataFrame()
gwtDf['time'] = gwtDf['time']/86400/365
gwtDf = gwtDf.set_index('time')
gwtDf| OBS1 | OBS2 | OBS3 | OBS4 | OBS5 | |
|---|---|---|---|---|---|
| time | |||||
| 5.0 | 107.286321 | 60.666416 | 1.349865 | 0.664076 | 0.227713 |
| 10.0 | 339.175837 | 246.724857 | 9.548925 | 21.047711 | 45.011523 |
| 15.0 | 558.399591 | 452.015263 | 31.501198 | 71.914114 | 143.841551 |
| 20.0 | 772.058932 | 626.672568 | 71.062098 | 138.455888 | 265.740493 |
| 25.0 | 902.509562 | 757.446171 | 126.441138 | 207.338829 | 391.997757 |
| 30.0 | 972.281811 | 849.016516 | 192.773815 | 264.863908 | 509.393987 |
| 35.0 | 1007.302296 | 910.365695 | 263.757513 | 306.632472 | 613.392793 |
| 40.0 | 1024.571750 | 950.030272 | 333.965065 | 334.390513 | 702.609107 |
| 45.0 | 1034.216752 | 975.147890 | 400.164664 | 352.186816 | 777.010479 |
| 50.0 | 1039.716124 | 990.717337 | 460.434205 | 363.492754 | 837.433407 |
fig, ax = plt.subplots(figsize=(12,8))
gwtDf.plot(ax=ax)
ax.grid()
Input data
You can download the input data from this link:
owncloud.hatarilabs.com/s/PvfQ8SGGk68BUS2
Password to download data: Hatarilabs
