Coupling geology from a shapefile to a Voronoi mesh MODFLOW6 model - Tutorial

/

A practical approach to set up hydraulic conductivity to a MODFLOW 6 model based on geological polygons. The tutorial show all steps from mesh and model creation, the setup of hydraulic parameters and the 3d representation of the final distribution in Paraview.

Tutorial

Código

from mf6Voronoi.utils import initiateOutputFolder
import flopy
from mf6Voronoi.utils import listTemplates, copyTemplate
#get the latest modflow binary
binDir = '../bin'
initiateOutputFolder(binDir)
The output folder ../bin exists and has been cleared
flopy.utils.get_modflow(binDir, subset=['mf6'])
fetched release '23.0' info from MODFLOW-ORG/executables
using previous download 'C:\Users\saulm\Downloads\modflow_executables-23.0-win64.zip' (use 'force=True' to re-download)
extracting 1 file to 'C:\Users\saulm\Documents\howToHydraulicCondFromGeology\bin'
mf6.exe (6.6.3)
updated flopy metadata file: 'C:\Users\saulm\AppData\Local\flopy\get_modflow.json'
#listTemplates()
copyTemplate('generateVoronoi','geo')
copyTemplate('modelCreation','geo')
copyTemplate('vtkGeneration','geo')

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='coupleGeology',maxRef = 250, multiplier=1, overlapping=False) ## Org

#Open limit layers and refinement definition layers
vorMesh.addLimit('basin','../shp/hatariUtils/catchment.shp') ## Org
vorMesh.addLayer('river','../shp/hatariUtils/river_basin.shp',40) ## Org
vorMesh.addLayer('gunsight','../shp/gunsightFormation.shp',80) ## Org
#Generate point pair array
vorMesh.generateOrgDistVertices() ## Org

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

mf6Voronoi will have a web version in 2028

Follow us:

 Hatarilabs Hatarilabs Hatarilabs Hatarilabs Hatarilabs Hatarilabs 
/--------Layer river discretization-------/
Progressive cell size list: [40, 80, 120, 160, 200, 240] m.

/--------Layer gunsight discretization-------/
Progressive cell size list: [80, 160, 240] m.

/----Sumary of points for voronoi meshing----/
Distributed points from layers: 2
Points from layer buffers: 5478
Points from max refinement areas: 551
Points from min refinement areas: 1242
Total points inside the limit: 8355
/--------------------------------------------/

Time required for point generation: 1.97 seconds 


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

Time required for voronoi generation: 1.44 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.005)
/----Performing quality verification of voronoi mesh----/
Short side on polygon: 8358 with length = 0.00036
Short side on polygon: 8358 with length = 0.00036
Short side on polygon: 8358 with length = 0.00029
Short side on polygon: 8358 with length = 0.00029
Short side on polygon: 8358 with length = 0.00194
Short side on polygon: 8358 with length = 0.00194
Short side on polygon: 8358 with length = 0.00395
Short side on polygon: 8358 with length = 0.00395
vorMesh.fixVoronoiShortSides()
vorMesh.generateVoronoi()
vorMesh.checkVoronoiQuality(threshold=0.005)
/----Generation of the voronoi mesh----/

Time required for voronoi generation: 1.59 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 has been generated.

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

Time required for voronoi shapefile: 1.94 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%|███████████████████████████████████████████████████████████████████████████| 8363/8363 [00:00<00:00, 18441.82it/s]



Extracting cell2d data and grid index


100%|████████████████████████████████████████████████████████████████████████████| 8363/8363 [00:03<00:00, 2451.03it/s]
#create folder
initiateOutputFolder('../json') ## Org

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

Model creation on steady state

Part 2a: generate disv properties

import sys, json, os ## Org
import rasterio, flopy ## Org
import numpy as np ## Org
import matplotlib.pyplot as plt ## Org
import geopandas as gpd ## Org
from mf6Voronoi.meshProperties import meshShape ## Org
from shapely.geometry import MultiLineString ## Org
# 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
src = rasterio.open('../rst/modelDem.tif')  ## Org
elevation=[x for x in src.sample(centroids)] ## Org
nlay = 15   ## Org

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


AcuifInf_Bottom = 1300 ## Org
zbot[0,] = AcuifInf_Bottom + (0.95 * (mtop - AcuifInf_Bottom)) ## Org
zbot[1,] = AcuifInf_Bottom + (0.9 * (mtop - AcuifInf_Bottom)) ## Org
zbot[2,] = AcuifInf_Bottom + (0.85 * (mtop - AcuifInf_Bottom)) ## Org
zbot[3,] = AcuifInf_Bottom + (0.8 * (mtop - AcuifInf_Bottom)) ## Org
zbot[4,] = AcuifInf_Bottom + (0.75 * (mtop - AcuifInf_Bottom)) ## Org
zbot[5,] = AcuifInf_Bottom + (0.7 * (mtop - AcuifInf_Bottom)) ## Org
zbot[6,] = AcuifInf_Bottom + (0.65 * (mtop - AcuifInf_Bottom)) ## Org
zbot[7,] = AcuifInf_Bottom + (0.6 * (mtop - AcuifInf_Bottom)) ## Org
zbot[8,] = AcuifInf_Bottom + (0.55 * (mtop - AcuifInf_Bottom)) ## Org
zbot[9,] = AcuifInf_Bottom + (0.5 * (mtop - AcuifInf_Bottom)) ## Org
zbot[10,] = AcuifInf_Bottom + (0.4 * (mtop - AcuifInf_Bottom)) ## Org
zbot[11,] = AcuifInf_Bottom + (0.3 * (mtop - AcuifInf_Bottom)) ## Org
zbot[12,] = AcuifInf_Bottom + (0.2 * (mtop - AcuifInf_Bottom)) ## Org
zbot[13,] = AcuifInf_Bottom + (0.1 * (mtop - AcuifInf_Bottom)) ## Org
zbot[14,] = AcuifInf_Bottom

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
                             sim_ws=modelWs) ## Org
# create tdis package
tdis_rc = [(1000.0, 1, 1.0)] ## Org
tdis = flopy.mf6.ModflowTdis(sim, pname='tdis', time_units='SECONDS', ## Org
                             perioddata=tdis_rc) ## Org
# create gwf model
gwf = flopy.mf6.ModflowGwf(sim, ## Org
                           modelname=modelName, ## Org
                           save_flows=True, ## Org
                           newtonoptions="NEWTON UNDER_RELAXATION") ## Org
# create iterative model solution and register the gwf model with it
ims = flopy.mf6.ModflowIms(sim, ## Org
                           complexity='COMPLEX', ## Org
                           outer_maximum=50, ## Org
                           inner_maximum=30, ## Org
                           linear_acceleration='BICGSTAB') ## Org
sim.register_ims_package(ims,[modelName]) ## Org
# disv
disv = flopy.mf6.ModflowGwfdisv(gwf, nlay=nlay, ncpl=ncpl, ## Org
                                top=mtop, botm=zbot, ## Org
                                nvert=nvert, vertices=vertices, ## Org
                                cell2d=cell2d) ## Org
# initial conditions
ic = flopy.mf6.ModflowGwfic(gwf, strt=np.stack([mtop for i in range(nlay)])) ## Org
KxArray = np.ones((nlay, ncpl)) * 4e-4
KxArray[1:6] = 5e-6
KxArray[6:15] = 3e-6
KxArray[15:] = 7e-7
icelltype = [1 for x in range(7)] + [0 for x in range(nlay - 7)] 

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

geoDf = gpd.read_file('../shp/geology.shp')

#looping over the geology polygons
for index, row in geoDf.iterrows():
    geoCells=interIx.intersect(row.geometry).cellids

    # for aluvial layers, that has 'alluvium' or 'deposit' on the notes
    if 'alluvium' in row.Notes or 'deposits' in row.Notes:
        for cell in geoCells:
            KxArray[0,cell] = row.Cond #only apply to the first layer
    #for the rest of layers
    else:
        for cell in geoCells:
            KxArray[1:,cell] = row.Cond #apply to the remaining layers below layer 1
        
# node property flow
npf = flopy.mf6.ModflowGwfnpf(gwf, ## Org
                              save_specific_discharge=True, ## Org
                              icelltype=icelltype, ## Org
                              k=KxArray,
                             k33=KxArray/10) ## Org
#cross section over the main river
crossSection = gpd.read_file('../shp/crossSection.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
modelxsect.plot_array(np.log(npf.k.array), cmap='viridis', alpha=0.5)
ax.grid() ## Org
sim.write_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...

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 = '../modelFiles' ## 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 solution package mf6model...
vtkDir = '../vtk' ## Org
initiateOutputFolder(vtkDir) ## Org

mf6Vtk = Mf6VtkGenerator(sim, vtkDir) ## Org
The output folder ../vtk has been generated.
Hatarilabs

mf6Voronoi will have a web version in 2028

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']
mf6Vtk.loadModel(modelName) ## Org
Package list: ['DISV', 'IC', 'NPF']
#generate model geometry as vtk and parameter array
mf6Vtk.generateGeometryArrays() ## Org
#generate parameter vtk
mf6Vtk.generateParamVtk() ## Org
Parameter Vtk Generated

Input data

You can download the input data from this link:

owncloud.hatarilabs.com/s/PBswAAxpBCU7VyJ

/ / Comment

Saul Montoya

Saul Montoya es Ingeniero Civil graduado de la Pontificia Universidad Católica del Perú en Lima con estudios de postgrado en Manejo e Ingeniería de Recursos Hídricos (Programa WAREM) de la Universidad de Stuttgart con mención en Ingeniería de Aguas Subterráneas y Hidroinformática.

 

Suscribe to our online newsletter

Subscribe for free newsletter, receive news, interesting facts and dates of our courses in water resources.