Modelamiento de pluma contaminante de presa de relaves con MODFLOW6 y mf6Voronoi - Tutorial

MODFLOW6 con mf6Voronoi puede manejar la simulación de tareas específicas relacionadas con la industria minera. Ya hemos cubierto la simulación de entradas de agua en el tajo mediante mallas de Voronoi y ahora vamos a modelar la filtración desde una instalación de almacenamiento de relaves. El caso aplicado abarca todos los pasos, desde la construcción de la malla, el modelo de flujo y el modelo de transporte. Los resultados se representan en 2D y se exportan a 3D en formato Vtk.

Tutorial

Codigo

Parte 1

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='tailingsStorage',maxRef = 100, multiplier=1.4, overlapping=False) ## Org

#Open limit layers and refinement definition layers
vorMesh.addLimit('basin','../shp/tsf/modelLimit.shp') ## Org
vorMesh.addLayer('tailingsseepage','../shp/tsf/tailingsSeepage.shp',25) ## Org
vorMesh.addLayer('tailings','../shp/tsf/tailingsEnvelope_v2.shp',40)
vorMesh.addLayer('river','../hatariUtils/river_network.shp',25)
#Generate point pair array
vorMesh.generateOrgDistVertices() ## Org

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

build faster, analyze more

Follow us:

Hatarilabs Hatarilabs Hatarilabs Hatarilabs Hatarilabs Hatarilabs
/--------Layer tailingsseepage discretization-------/
Progressive cell size list: [25, 60.0] m.

/--------Layer tailings discretization-------/
Progressive cell size list: [40, 96.0] m.

/--------Layer river discretization-------/
Progressive cell size list: [25, 60.0] m.

/----Sumary of points for voronoi meshing----/
Distributed points from layers: 3
Points from layer buffers: 1686
Points from max refinement areas: 2768
Points from min refinement areas: 1967
Total points inside the limit: 7361
/--------------------------------------------/

Time required for point generation: 1.60 seconds 


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

Time required for voronoi generation: 0.74 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: 7360 with length = 0.00791
Short side on polygon: 7360 with length = 0.00791
vorMesh.fixVoronoiShortSides()
vorMesh.generateVoronoi()
vorMesh.checkVoronoiQuality(threshold=0.01)
/----Generation of the voronoi mesh----/

Time required for voronoi generation: 0.67 seconds 


/----Performing quality verification of voronoi mesh----/
Your mesh has no edges shorter than your threshold
#Export generated voronoi mesh
initiateOutputFolder('../output') ## Org
getVoronoiAsShp(vorMesh.modelDis, shapePath='../output/'+vorMesh.modelDis['meshName']+'.shp') ## Org
The output folder ../output exists and has been cleared

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

Time required for voronoi shapefile: 1.86 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() ## Org
Creating a unique list of vertices [[x1,y1],[x2,y2],...]


100%|███████████████████████████████████████████████████████████████████████████| 7362/7362 [00:00<00:00, 17041.77it/s]



Extracting cell2d data and grid index


100%|████████████████████████████████████████████████████████████████████████████| 7362/7362 [00:02<00:00, 2571.43it/s]
#create folder
initiateOutputFolder('../json') ## Org

#export disv
mesh.save_properties('../json/disvDict.json') ## Org
The output folder ../json exists and has been cleared

Parte 2

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 ## <==== updated
from mf6Voronoi.tools.cellWork import getLayCellElevTupleFromRaster, getLayCellElevTupleFromElev  ## <==== updated
C:\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) ## 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
#first 4 layers
finalSurf = rasterio.open('../rst/demWithDamAndTailings.tif')  ## <==== updated
finalSurfElev = [x for x in finalSurf.sample(centroids)] ## <==== updated

#second 4 layers
botTailings = rasterio.open('../rst/demWithDam.tif') ## <==== updated
botTailingsElev = [x for x in botTailings.sample(centroids)] ## <==== updated

orgSurf = rasterio.open('../rst/tsfDem10m.tif') ## <==== updated
orgSurfElev = [x for x in orgSurf.sample(centroids)] ## <==== updated
nlay = 19   ## Org

fsurf=np.array([elev[0] for i,elev in enumerate(finalSurfElev)]) ## <==== updated
#4 layers of 2m for tailings
btail=np.array([elev[0] - 8 for i,elev in enumerate(botTailingsElev)]) ## <==== updated
#bottom of 1m geomembrane
geob=np.array([elev[0] - 9 for i,elev in enumerate(botTailingsElev)])
#4 layers of 2m for dam 
osurf=np.array([elev[0] - 17 for i,elev in enumerate(orgSurfElev)]) ## <==== updated

zbot=np.zeros((nlay,ncpl)) ## Org
mtop = fsurf 
AcuifInf_Bottom = 170 ## <==== updated

zbot[0,] = btail + 0.75*(fsurf - btail) ## <==== updated
zbot[1,] = btail + 0.5*(fsurf - btail) ## <==== updated
zbot[2,] = btail + 0.25*(fsurf - btail) ## <==== updated
zbot[3,] = btail ## <==== updated
zbot[4,] = geob ## <==== updated
zbot[5,] = osurf + 0.75*(geob - osurf) ## <==== updated
zbot[6,] = osurf + 0.5*(geob - osurf) ## <==== updated
zbot[7,] = osurf + 0.25*(geob - osurf) ## <==== updated
zbot[8,] = osurf ## <==== updated
zbot[9,] = AcuifInf_Bottom + (0.95 * (osurf - AcuifInf_Bottom)) ## <==== updated
zbot[10,] = AcuifInf_Bottom + (0.90 * (osurf - AcuifInf_Bottom)) ## <==== updated
zbot[11,] = AcuifInf_Bottom + (0.85 * (osurf - AcuifInf_Bottom)) ## <==== updated 85%
zbot[12,] = AcuifInf_Bottom + (0.78 * (osurf - AcuifInf_Bottom)) ## <==== updated 
zbot[13,] = AcuifInf_Bottom + (0.71 * (osurf - AcuifInf_Bottom)) ## <==== updated 
zbot[14,] = AcuifInf_Bottom + (0.64 * (osurf - AcuifInf_Bottom)) ## <==== updated 
zbot[15,] = AcuifInf_Bottom + (0.57 * (osurf - AcuifInf_Bottom)) ## <==== updated 
zbot[16,] = AcuifInf_Bottom + (0.50 * (osurf - AcuifInf_Bottom)) ## <==== updated 50%
zbot[17,] = AcuifInf_Bottom + (0.25 * (osurf - AcuifInf_Bottom)) ## <==== updated
zbot[18,] = AcuifInf_Bottom ## <==== updated

Create 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 = [(80*365*86400, 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
waterTableRst = rasterio.open('../rst/regWaterTable.tif')  ## <==== updated
waterTable = [x for x in waterTableRst.sample(centroids)] ## <==== updated
#ic = flopy.mf6.ModflowGwfic(gwf, strt=np.stack([mtop for i in range(nlay)])) ## Org
ic = flopy.mf6.ModflowGwfic(gwf, strt=np.stack([waterTable for i in range(nlay)])) ## <==== updated
ic.plot(mflay=0)
ncplArray = np.ones(ncpl)
Kx =[4E-4*ncplArray for i in range(9)] + [5e-6*ncplArray for i in range(3)] + \
    [8e-7*ncplArray for i in range(5)] + [5e-8*ncplArray for i in range(2)] ## <==== updated
kxArray = np.stack(Kx)
icelltype = [1 for i in range(14)] + [0 for i in range(5)] ## <=== updated

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

damGeom = gpd.read_file('../shp/tsf/damEnvelope_v2.shp').iloc[0].geometry
damCells = interIx.intersect(damGeom).cellids
kxArray[0:9,damCells.tolist()] = 2e-5


tailingsGeom = gpd.read_file('../shp/tsf/tailingsEnvelope_v2.shp').iloc[0].geometry
tailingsCells = interIx.intersect(tailingsGeom).cellids
kxArray[:4,tailingsCells.tolist()] = 5e-6

#falta la geomembrana


# 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=kxArray, ## <==== inserted
                              k33=kxArray) ## <==== inserted
## <==== inserted
crossSection = gpd.read_file('../shp/tsf/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
                              ) ## 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_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] ## Org
array([279, 289, 316, 322, 349, 350, 358, 381, 384, 428], 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/regWaterTable.tif',
                                                               '../shp/tsf/regionalFlow.shp') ## # <===== Inserted
regionalSpd = {} ## # <===== Inserted
regionalSpd[0] = [] ## # <===== Inserted
for index, layCellTuple in enumerate(layCellTupleList): ## Org
    regionalSpd[0].append([layCellTuple,cellElevList[index],0.01]) # <===== Inserted
The cell 322 has a elevation of 519.06 outside the model vertical domain
The cell 323 has a elevation of 518.79 outside the model vertical domain
The cell 326 has a elevation of 519.16 outside the model vertical domain
The cell 5531 has a elevation of 715.13 outside the model vertical domain
The cell 5538 has a elevation of 717.05 outside the model vertical domain
The cell 5539 has a elevation of 715.39 outside the model vertical domain
The cell 5541 has a elevation of 718.16 outside the model vertical domain
ghb = flopy.mf6.ModflowGwfghb(gwf, stress_period_data=regionalSpd) ## <==== modified
#regional flow plot
ghb.plot(mflay=12, 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() ## 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 533 based on size of stress_period_data
    writing package ghb_0...
INFORMATION: maxbound in ('gwf6', 'ghb', 'dimensions') changed to 289 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/01 15:40:34
 
 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/01 15:43:47
 Elapsed run time:  3 Minutes, 13.222 Seconds
 

WARNING REPORT:

  1.  Simulation convergence failure occurred 6 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] ## Org
array([544.32303871, 548.36449819, 541.34403696, 539.43147163,
       546.02745945])
# 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='../output/tailingsStorage.shp', 
                        rasterRes=10, 
                        epsg=32614, 
                        outputPath='../output/waterTable.tif')

Parte 3

#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 ## Org
C:\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) ## Org
loading 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 

layCellTupleList, cellElevList = getLayCellElevTupleFromRaster(gwf,
                                                               interIx,
                                                               '../output/waterTable.tif',
                                                               '../shp/tsf/tailingsEnvelope_v2.shp') 

for lay, cell in layCellTupleList:
    strtConc[lay,cell] = 1200

gwtIc = flopy.mf6.ModflowGwtic(gwt, strt=strtConc) ## Org
The cell 6283 has a elevation of 670.71 outside the model vertical domain
The cell 6287 has a elevation of 670.88 outside the model vertical domain
The cell 6288 has a elevation of 670.24 outside the model vertical domain
fig = 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
## <==== inserted
crossSection = gpd.read_file('../shp/tsf/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(strtConc,masked_values=[-1e+30], cmap=plt.cm.summer) ## <====== inserted
modelxsect.plot_grid()
ax.grid() ## <==== inserted
# 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) ## Org
cncSpd = {} ## Org
cncSpd[0] = [] ## Org

layCellTupleList, cellElevList = getLayCellElevTupleFromRaster(gwf,
                                                               interIx,
                                                               '../output/waterTable.tif',
                                                               '../shp/tsf/tailingsEnvelope_v2.shp') 

for layCell in layCellTupleList:
    cncSpd[0].append([layCell[0],layCell[1],1200])
    
cnc = flopy.mf6.ModflowGwtcnc(gwt,stress_period_data=cncSpd) ## Org
cnc.plot(mflay=2) ## Org
The cell 6283 has a elevation of 670.71 outside the model vertical domain
The cell 6287 has a elevation of 670.88 outside the model vertical domain
The cell 6288 has a elevation of 670.24 outside the model vertical domain
The cell 6295 has a elevation of 670.47 outside the model vertical domain
The cell 6336 has a elevation of 670.22 outside the model vertical domain
obsList = []
nameList, obsLayCellList = getLayCellElevTupleFromObs(gwf, ## Org
                  interIx, ## Org
                  '../shp/tsf/piezometers.shp', ## Org
                  'Name', ## Org
                  'Elev') ## 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 1359
Well screen elev of 529.00 found at layer 9
Working for cell 1973
Well screen elev of 533.00 found at layer 9
Working for cell 2441
Well screen elev of 537.00 found at layer 9
#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"), ## Org
    ("GWFBUDGET", "../modelFiles/mf6Model.cbc"), ## Org
]

fmi = flopy.mf6.ModflowGwtfmi(gwt, ## Org
                              flow_imbalance_correction=True, ## Org
                              packagedata=pd) ## Org
sim.write_simulation() ## Org
success, buff = sim.run_simulation() ## Org
writing 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 1658 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/01 15:50:53
 
 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/01 15:54:04
 Elapsed run time:  3 Minutes, 10.532 Seconds
 

WARNING REPORT:

  1.  Simulation convergence failure occurred 3 time(s).
 Normal termination of simulation.
concObj = gwt.output.concentration() ## Org
concObj.get_times() ## Org
[2522880000.0]
conc = concObj.get_data() ## Org
conc.shape ## Org
(19, 1, 7362)
transAoi = gpd.read_file('../shp/tsf/tailingsSeepage.shp') ## Org
xMin, yMin, xMax, yMax = transAoi.bounds.iloc[0].values ## Org
mflay = 7 ## 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,xMax)
ax.set_ylim(yMin,yMax)

plt.colorbar(cellConc, shrink=0.25) 

plt.show()

Part 4

#Vtk generation
import flopy ## Org
from mf6Voronoi.tools.vtkGen import Mf6VtkGenerator ## Org
from mf6Voronoi.utils import initiateOutputFolder ## Org
C:\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) ## Org
loading 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) ## Org
The output folder ../vtk exists and has been cleared
Hatarilabs

build faster, analyze more

Follow us:

Hatarilabs Hatarilabs Hatarilabs Hatarilabs Hatarilabs Hatarilabs
/---------------------------------------/

The Vtk generator engine has been started

/---------------------------------------/
#list models on the simulation
mf6Vtk.listModels() ## Org
Models in simulation: ['mf6model', 'gwtmodel']
mf6Vtk.loadModel(modelName) ## Org
Package 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() ## Org
Parameter 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.0802 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.0780 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.9370 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 1.6860 seconds to be generated.
/--------GHB_0 vtk generated-------/
mf6Vtk.generateHeadVtk(nper=0, crop=True) ## Org
mf6Vtk.generateWaterTableVtk(nper=0) ## Org
gwt = sim.get_model('gwtmodel')
concObj = gwt.output.concentration()
concObj.get_times()
[2522880000.0]
conc = concObj.get_data(totim=2522880000.0)
conc[4][0]
array([0., 0., 0., ..., 0., 0., 0.])
mf6Vtk.generateArrayVtk(conc, 'conc80y', nper=0,nstp=9, crop=True)

Datos de entrada

Puede descargar los datos de entrada desde este link:

owncloud.hatarilabs.com/s/xWc7GNXhSRhEFRx

Password para descargar: Hatarilabs

 

Suscríbete a nuestro boletín electrónico

Suscríbase a nuestro boletín gratuito para recibir noticias, datos interesantes y fechas de nuestros cursos en recursos hídricos.

 

Posted on August 4, 2025 and filed under TutorialPython, TutorialModflow, Modelamiento, Minería.