Groundwater seepage from a water reservoir modeling with Modflow6 and mf6Voronoi - Tutorial
/
Applied case of local groundwater modeling based on a regional groundwater flow to evaluate the interaction of the water reservoir with the groundwater flow system and the seepage rate from the reservoir. The model runs on steady state and shows the procedure to modify original/final surfaces without error on layer thickness.
Tutorial
Code
#!pip install -U mf6Voronoifrom mf6Voronoi.utils import listTemplates, copyTemplate#listTemplates()copyTemplate('generateVoronoi','reservoir')copyTemplate('modelCreation','reservoir')copyTemplate('vtkGeneration','reservoir')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='reservoirModel',maxRef = 250, multiplier=1, overlapping=False) ## Org
#Open limit layers and refinement definition layers
vorMesh.addLimit('basin','../../shp/reservoirModelLimit.shp') ## Org
vorMesh.addLayer('dam','../../shp/damEnvelopeDownstream.shp',20) ## Org
vorMesh.addLayer('reservoir','../../shp/reservoirEnvelope.shp',80) ## 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
Follow us: |
 |
 |
 |
 |
 |
 |
/--------Layer dam discretization-------/
Progressive cell size list: [20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240] m.
/--------Layer reservoir discretization-------/
Progressive cell size list: [80, 160, 240] m.
/--------Layer river discretization-------/
Progressive cell size list: [80, 160, 240] m.
/----Sumary of points for voronoi meshing----/
Distributed points from layers: 3
Points from layer buffers: 4855
Points from max refinement areas: 2111
Points from min refinement areas: 3713
Total points inside the limit: 12484
/--------------------------------------------/
Time required for point generation: 2.87 seconds
/----Generation of the voronoi mesh----/
Time required for voronoi generation: 1.15 seconds#Uncomment the next two cells if you have strong differences on discretization or you have encounter an FORTRAN error while running MODFLOW6vorMesh.checkVoronoiQuality(threshold=0.01)/----Performing quality verification of voronoi mesh----/
Short side on polygon: 12483 with length = 0.00380
Short side on polygon: 12483 with length = 0.00380
Short side on polygon: 12483 with length = 0.00221
Short side on polygon: 12483 with length = 0.00221
Short side on polygon: 12483 with length = 0.00673
Short side on polygon: 12483 with length = 0.00673
Short side on polygon: 12483 with length = 0.00007
Short side on polygon: 12483 with length = 0.00007
Short side on polygon: 12483 with length = 0.00264
Short side on polygon: 12483 with length = 0.00264
Short side on polygon: 12483 with length = 0.00154
Short side on polygon: 12483 with length = 0.00154
Short side on polygon: 12483 with length = 0.00003
Short side on polygon: 12483 with length = 0.00003
Short side on polygon: 12483 with length = 0.00003
Short side on polygon: 12483 with length = 0.00008
Short side on polygon: 12483 with length = 0.00008
Short side on polygon: 12483 with length = 0.00003
Short side on polygon: 12483 with length = 0.00008
Short side on polygon: 12483 with length = 0.00008vorMesh.fixVoronoiShortSides()
vorMesh.generateVoronoi()
vorMesh.checkVoronoiQuality(threshold=0.01)/----Generation of the voronoi mesh----/
Time required for voronoi generation: 1.21 seconds
/----Performing quality verification of voronoi mesh----/
Your mesh has no edges shorter than your threshold#Export generated voronoi mesh
initiateOutputFolder('../../reservoirFiles/output') ## Org
getVoronoiAsShp(vorMesh.modelDis, shapePath='../../reservoirFiles/output/'+vorMesh.modelDis['meshName']+'.shp') ## OrgThe output folder ../../reservoirFiles/output has been generated.
/----Generation of the voronoi shapefile----/
Time required for voronoi shapefile: 2.54 seconds# Show the resulting voronoi mesh
#open the mesh file
mesh=gpd.read_file('../../reservoirFiles/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('../../reservoirFiles/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%|█████████████████████████████████████████████████████████████████████████| 12494/12494 [00:00<00:00, 17930.04it/s]
Extracting cell2d data and grid index
100%|██████████████████████████████████████████████████████████████████████████| 12494/12494 [00:03<00:00, 3265.71it/s]#create folder
initiateOutputFolder('../../reservoirFiles/json') ## Org
#export disv
mesh.save_properties('../../reservoirFiles/json/disvDict.json') ## OrgThe output folder ../../reservoirFiles/json exists and has been clearedModel 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
from mf6Voronoi.tools.cellWork import getLayCellElevTupleFromRaster, getLayCellElevTupleFromElev# open the json file
with open('../../reservoirFiles/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/n33w111_wgs84_int32_20m.tif') ## Org
elevation=[x for x in src.sample(centroids)] ## Org
finalSrc = rasterio.open('../../rst/dam/demwithReservoir.tif')
finalElevation=[x for x in finalSrc.sample(centroids)]#fix elevation to stay allways below final elevation
for index, elev in enumerate(elevation):
if elev > finalElevation[index]:
elevation[index] = finalElevation[index]
else:
passnlay = 14 #10 ## Org
#original and final surface
orgMtop = np.array([elev[0] for i,elev in enumerate(elevation)]) ## Bottom of layer 1
mtop = np.array([elev[0] for i,elev in enumerate(finalElevation)])
zbot=np.zeros((nlay,ncpl)) ## Org
AcuifInf_Bottom = 700 ## Org
#bottom original layer 1
botLay1 = np.copy(AcuifInf_Bottom + (0.95 * (orgMtop - AcuifInf_Bottom)))
#added layers
zbot[0,] = botLay1 + 0.8 * (mtop - botLay1)
zbot[1,] = botLay1 + 0.6 * (mtop - botLay1)
zbot[2,] = botLay1 + 0.4 * (mtop - botLay1)
zbot[3,] = botLay1 + 0.2 * (mtop - botLay1)
#original layers but different index
zbot[4,] = AcuifInf_Bottom + (0.95 * (orgMtop - AcuifInf_Bottom)) ## <==== updated
zbot[5,] = AcuifInf_Bottom + (0.90 * (orgMtop - AcuifInf_Bottom)) ## <==== updated
zbot[6,] = AcuifInf_Bottom + (0.85 * (orgMtop - AcuifInf_Bottom)) ## <==== updated 85%
zbot[7,] = AcuifInf_Bottom + (0.78 * (orgMtop - AcuifInf_Bottom)) ## <==== updated
zbot[8,] = AcuifInf_Bottom + (0.71 * (orgMtop - AcuifInf_Bottom)) ## <==== updated
zbot[9,] = AcuifInf_Bottom + (0.64 * (orgMtop - AcuifInf_Bottom)) ## <==== updated
zbot[10,] = AcuifInf_Bottom + (0.57 * (orgMtop - AcuifInf_Bottom)) ## <==== updated
zbot[11,] = AcuifInf_Bottom + (0.50 * (orgMtop - AcuifInf_Bottom)) ## <==== updated 50%
zbot[12,] = AcuifInf_Bottom + (0.25 * (orgMtop - AcuifInf_Bottom)) ## <==== updated
zbot[13,] = AcuifInf_Bottom ## <==== updatedCreate simulation and model
# create simulation
simName = 'mf6Sim' ## Org
modelName = 'mf6Model' ## Org
modelWs = '../../reservoirFiles/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) ## Orgdisv.top.plot(figsize=(12,8), alpha=0.8) ## Org
crossSection = gpd.read_file('../../shp/crossSectionReservoir.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.set_ylim(900,1300)
ax.set_xlim(1500,2500)
# initial conditions
waterTableSrc = rasterio.open('../../regionalFiles/output/waterTable.tif')
regWaterTable = [x for x in waterTableSrc.sample(centroids)]
ic = flopy.mf6.ModflowGwfic(gwf, strt=np.stack([regWaterTable for i in range(nlay)])) ## Orgic.plot(mflay=0)
kxOnes = np.ones([nlay,ncpl])
kxOnes.shape
#original matrix of Kx values
Kx =[4E-4, 4E-4, 4E-4, 4E-4, 4E-4, 5E-5, 3E-6, 3E-6, 2.5E-6, 2.5E-6, 2.5E-6, 1E-6, 9E-7, 5E-7] ## Org
layerValues = np.array(Kx).reshape(nlay,1)
kxArray = kxOnes*layerValues
kxArrayarray([[4.e-04, 4.e-04, 4.e-04, ..., 4.e-04, 4.e-04, 4.e-04],
[4.e-04, 4.e-04, 4.e-04, ..., 4.e-04, 4.e-04, 4.e-04],
[4.e-04, 4.e-04, 4.e-04, ..., 4.e-04, 4.e-04, 4.e-04],
...,
[1.e-06, 1.e-06, 1.e-06, ..., 1.e-06, 1.e-06, 1.e-06],
[9.e-07, 9.e-07, 9.e-07, ..., 9.e-07, 9.e-07, 9.e-07],
[5.e-07, 5.e-07, 5.e-07, ..., 5.e-07, 5.e-07, 5.e-07]])icelltype = [1,1,1,1,1,1,1,1,1,1,0,0,0,0] ## Org
# Define intersection object
interIx = flopy.utils.gridintersect.GridIntersect(gwf.modelgrid) ## Org
#reservoir apply kx to cells in layers from 1 to 5
resDf = gpd.read_file('../../shp/reservoirEnvelope.shp')
resCells = interIx.intersect(resDf.iloc[0].geometry).cellids
for cell in resCells:
for layer in range(5):
kxArray[layer,cell] = 1
#dam apply kx to cells in layers from 1 to 5 from cells downstream
damDf = gpd.read_file('../../shp/damEnvelopeDownstream.shp')
damCells = interIx.intersect(damDf.iloc[0].geometry).cellids
for cell in damCells:
for layer in range(6):
kxArray[layer,cell] = 1e-8
#apply dam kx to layers upstream
layCellTupleList, cellElevList = getLayCellElevTupleFromRaster(gwf,
interIx,
'../../rst/dam/damElev.tif',
'../../shp/damEnvelopeUpstream.shp') ## # <=====
#looping over celltuples to get cells from intersected layer to layer 5
for layCell in layCellTupleList:
#filter layers below layer 5
if layCell[0] <=5:
for extendedLayer in range(layCell[0],6):
kxArray[extendedLayer,layCell[1]] = 1e-8The cell 1869 has a elevation of 1221.23 outside the model vertical domain
The cell 1877 has a elevation of 1228.92 outside the model vertical domain
The cell 1886 has a elevation of 1229.89 outside the model vertical domain
The cell 1897 has a elevation of 1222.93 outside the model vertical domain
The cell 1903 has a elevation of 1228.70 outside the model vertical domain
The cell 1927 has a elevation of 1228.41 outside the model vertical domain# node property flow
npf = flopy.mf6.ModflowGwfnpf(gwf, ## Org
save_specific_discharge=True, ## Org
icelltype=icelltype, ## Org
k=kxArray) ## Org## <==== inserted
crossSection = gpd.read_file('../../shp/crossSectionReservoir.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
ax.set_ylim(900,1300)
ax.set_xlim(1600,2600)
# 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
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) ## OrgDefinition 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_clip.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([48, 55, 59, 60, 72, 97, 101, 110, 116, 121], 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
#general boundary condition
#regional flow
layCellTupleList, cellElevList = getLayCellElevTupleFromRaster(gwf,
interIx,
'../../regionalFiles/output/waterTable.tif',
'../../shp/regionalFlow.shp') ## # <===== Inserted
regionalSpd = {} ## # <===== Inserted
regionalSpd[0] = [] ## # <===== Inserted
for index, layCellTuple in enumerate(layCellTupleList): ## Org
regionalSpd[0].append([layCellTuple,cellElevList[index],0.01,'regionalFlow']) # <===== InsertedThe cell 55 has a elevation of 1208.29 outside the model vertical domain
The cell 383 has a elevation of 1176.29 outside the model vertical domain
The cell 408 has a elevation of 1175.19 outside the model vertical domain
The cell 1042 has a elevation of 1098.56 outside the model vertical domain
The cell 1056 has a elevation of 1098.94 outside the model vertical domain
The cell 1072 has a elevation of 1099.35 outside the model vertical domain
The cell 7764 has a elevation of 1324.99 outside the model vertical domain
The cell 10935 has a elevation of 1242.91 outside the model vertical domain
The cell 11478 has a elevation of 1241.94 outside the model vertical domainlayCellTupleList = getLayCellElevTupleFromElev(gwf,
interIx,
1220,
'../../shp/reservoirEnvelope.shp')
for layCellTuple in layCellTupleList:
regionalSpd[0].append([layCellTuple, 1220, 10, 'reservoir'])
ghb = flopy.mf6.ModflowGwfghb(gwf, stress_period_data=regionalSpd, boundnames=True) ## <==== modified
# Observation package for Drain
obsDict = { # <===== Inserted
"{}.ghb.obs.csv".format(modelName): [ # <===== Inserted
("regional_flow", "ghb", "regionalFlow"), # <===== Inserted
("reservoir", "ghb", "reservoir") # <===== Inserted
] # <===== Inserted
} # <===== Inserted
# Attach observation package to DRN package
ghb.obs.initialize( # <===== Inserted
filename=gwf.name+".ghb.obs", # <===== Inserted
digits=10, # <===== Inserted
print_input=True, # <===== Inserted
continuous=obsDict # <===== Inserted
) # <===== Inserted
#ghb plot
ghb.plot(mflay=0, kper=0) # <===== InsertedYou have inserted a fixed elevation
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() ## 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 500 based on size of stress_period_data
writing package ghb_0...
INFORMATION: maxbound in ('gwf6', 'ghb', 'dimensions') changed to 4811 based on size of stress_period_data
writing package obs_0...
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/23 12:22:26
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/23 12:22:50
Elapsed run time: 24.049 Seconds
ERROR REPORT:
1. Simulation convergence failure occurred 1 time(s).
Premature termination of simulation.Model output visualization
headObj = gwf.output.head() ## Org
headObj.get_kstpkper() ## Org[(np.int32(0), np.int32(0))]heads = headObj.get_data() ## Org
heads[2,0,:5] ## Orgarray([1210.357455 , 1203.58311951, 1211.3184162 , 1204.7309277 ,
1207.04141144])# 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 the heads for a defined layer and boundary conditions
boundsDf = gpd.read_file('../../shp/damEnvelope.shp')
bounds = boundsDf.total_bounds
reservoirDf = gpd.read_file('../../shp/reservoirEnvelope.shp')
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.5) ## Org
linecollection = modelmap.plot_grid(linewidth=0.3, alpha=0.5, color='cyan', ax=ax) ## Org
ax.set_ylim(bounds[1]-500,bounds[3]+500)
ax.set_xlim(bounds[0]-1500,bounds[2]+1500)
reservoirDf.plot(ax=ax, alpha=0.5)
plt.colorbar(cellhead, shrink=0.75) ## Org
plt.show() ## Org
3d geometry generation on Vtk format
#Vtk generation
import flopy ## Org
from mf6Voronoi.tools.vtkGen import Mf6VtkGenerator ## Org
from mf6Voronoi.utils import initiateOutputFolder ## Org# load simulation
simName = 'mf6Sim' ## Org
modelName = 'mf6Model' ## Org
modelWs = '../../reservoirfiles/modelFiles' ## 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 solution package mf6model...vtkDir = '../../reservoirFiles/vtk' ## Org
initiateOutputFolder(vtkDir) ## Org
mf6Vtk = Mf6VtkGenerator(sim, vtkDir) ## OrgThe output folder ../../reservoirFiles/vtk has been generated.mf6Voronoi will have a web version in 2028
Follow us: |
|
|
|
|
|
|
/---------------------------------------/
The Vtk generator engine has been started
/---------------------------------------/#list models on the simulation
mf6Vtk.listModels() ## OrgModels in simulation: ['mf6model']mf6Vtk.loadModel(modelName) ## OrgPackage list: ['DISV', 'IC', 'NPF', 'STO', 'RCHA_0', 'EVTA_0', 'DRN_0', 'GHB_OBS', 'GHB_0', 'OC']#show output data
headObj = mf6Vtk.gwf.output.head() ## Org
headObj.get_kstpkper() ## Org[(np.int32(0), np.int32(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.2023 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.2495 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 3.4729 seconds to be generated.
/--------DRN_0 vtk generated-------/
/--------GHB_0 vtk generation-------/
Working for GHB_0 package, creating the datasets: ('bhead', 'cond', 'boundname')
Vtk file took 37.4578 seconds to be generated.
/--------GHB_0 vtk generated-------/mf6Vtk.generateHeadVtk(nper=0, crop=True) ## Orgmf6Vtk.generateWaterTableVtk(nper=0) ## OrgInput files
You can download the input files from this link:
owncloud.hatarilabs.com/s/PrAOTk6ppYkAyON
Password: Hatarilabs
