Basic Example of Groundwater Modeling in MODFLOW 6 and Visualization with Paraview and Flopy

Basic tutorial to learn the procedure to build, simulate and represent a MODFLOW 6 model. The tutorial shows a introduction to the model file system on steady state conditions. The model for this tutorial is implemented with the following boundary conditions: Drains, Recharge, Wells, and Constant Head. The grid is regular with a width of 50 meters and it has 30 rows and 24 columns; the model has 4 layers and a total thickness of 130 meters. The model is called "hatari01" and is inspired in the "twri" model from the MODFLOW 2005 documentation adapted to MODFLOW 6. The model defines a constant horizontal hydraulic conductivity as well as vertical conductivity. After the simulation a Python code is run on a Jupyter Notebook to create the Unstructured VTK files for the heads, water table and boundary conditions representation as 3D objects in Paraview.

Code

This is a part of the Python code to create the VTK files, the rest you will find on the input data part of this tutorial.

```# # Import packages, read files and create empty dicts

# In[1]:

import os, re
import numpy as np
from workFunctions import workFunctions
from vtkFunctions import vtkFunctions
from transFunctions import transFunctions
from listFunctions import listFunctions

# In[2]:

#open the DIS, BAS and FHD and DRN files

# In[3]:

#create a empty dictionay to store the model features
modChd  = {}
modDis  = {}
modDrn  = {}
modIc   = {}
modNpf  = {}
modRch  = {}
modWel  = {}
modHds  = {}

# <br/>
# <br/>
# <br/>
# ___
#
# # Working with the DIS (Discretization Data) data

# ### General model features as modDis dict

# In[4]:

########################
### General model features as modDis dict
#get the number of layers, rows, columns, cell and vertex numbers
for line in disLines:
if 'NLAY' in line:
modDis['cellLays'] = int(line.split()[1])
elif 'NROW' in line:
modDis['cellRows'] = int(line.split()[1])
elif 'NCOL' in line:
modDis['cellCols'] = int(line.split()[1])

modDis["vertexLays"] = modDis["cellLays"] + 1
modDis["vertexRows"] = modDis["cellRows"] + 1
modDis["vertexCols"] = modDis["cellCols"] + 1
modDis["vertexPerLay"] = modDis["vertexRows"] * modDis["vertexCols"]
modDis["cellsPerLay"] = modDis["cellRows"] * modDis["cellCols"]

########################
### Get the DIS Breakers
modDis['DELRArray1D'] = workFunctions.getListFromDel('DELR',modDis,disLines)
modDis['DELCArray1D'] = workFunctions.getListFromDel('DELC',modDis,disLines)

modDis['cellZVertexGrid']={}
modDis['cellZVertexGrid']['lay0']=workFunctions.getUniLayerListFromTerm(modDis,disLines,'TOP').reshape(modDis['cellRows'],modDis['cellCols'])

listFromBottom = workFunctions.getListinDictxLayFromGriddataLayered(modDis,disLines,'BOTM',modDis)
#i = 0
for lay in range(1,modDis['vertexLays']):
modDis['cellZVertexGrid']['lay'+str(lay)]=np.asarray(listFromBottom['lay'+str(lay-1)]).reshape(modDis['cellRows'],modDis['cellCols'])
#    i+=1

########################
### Geolocation model data
modDis["vertexXmin"]=0
modDis["vertexYmin"]=0
modDis["vertexXmax"]=sum(modDis['DELRArray1D'])
modDis["vertexYmax"]=sum(modDis['DELCArray1D'])

########################
### List of arrays of cells and vertex coord
modDis['vertexEastingArray1D']  = np.array([modDis['vertexXmin']+np.sum(modDis['DELRArray1D'][:col]) for col in range(modDis['vertexCols'])])
modDis['vertexNorthingArray1D'] = np.array([modDis['vertexYmax']-np.sum(modDis['DELCArray1D'][:row]) for row in range(modDis['vertexRows'])])

modDis['cellEastingArray1D']    = np.array([modDis['vertexXmin']+np.sum(modDis['DELRArray1D'][:col])+modDis['DELRArray1D'][col]/2 for col in range(modDis['cellCols'])])
modDis['cellNorthingArray1D']   = np.array([modDis['vertexYmax']-np.sum(modDis['DELCArray1D'][:row])-modDis['DELCArray1D'][row]/2 for row in range(modDis['cellRows'])])

########################
### Grid of XYZ Vertex Coordinates
modDis['vertexXGrid'] = np.repeat(modDis['vertexEastingArray1D'].reshape(modDis['vertexCols'],1),modDis['vertexRows'],axis=1).T
modDis['vertexYGrid'] = np.repeat(modDis['vertexNorthingArray1D'],modDis['vertexCols']).reshape(modDis['vertexRows'],modDis['vertexCols'])
modDis['vertexZGrid'] = transFunctions.interpolateCelltoVertex(modDis,'cellZVertexGrid')

# <br/>
# <br/>
# <br/>
# ___
#
# # Get the Info for Boundary Conditions and Cell Heads

# In[5]:

# Get the NPF Info
modNpf['iCellTypeList'] = workFunctions.getListinDictxLayFromGriddataLayered(modNpf,npfLines,'ICELLTYPE',modDis)
modNpf['kList'] = workFunctions.getListinDictxLayFromGriddataLayered(modNpf,npfLines,'k LAYERED',modDis)
modNpf['K33List'] = workFunctions.getListinDictxLayFromGriddataLayered(modNpf,npfLines,'K33 LAYERED',modDis)

# Get the IC Info
modIc['strtList'] = workFunctions.getListinDictxLayFromGriddataLayered(modIc,icLines,'STRT',modDis)

# Get the DRN Info
modDrn['maxBound'] = workFunctions.getTermFromKeyword(drnLines,'MAXBOUND','DIMENSIONS')
modDrn['drnCells'] = workFunctions.getCellsforBoundary(drnLines,'drn',modDrn['maxBound'],1)

# Get the CHD Info
modChd['maxBound'] = workFunctions.getTermFromKeyword(chdLines,'MAXBOUND','DIMENSIONS')
modChd['chdCells'] = workFunctions.getCellsforBoundary(chdLines,'chd',modChd['maxBound'],1)

# Get the RCH Info
modRch['rchCellList'] = workFunctions.getUniLayerListFromTerm(modDis,rchLines,'RECHARGE')

# Get the WEL Info
modWel['maxBound'] = workFunctions.getTermFromKeyword(welLines,'MAXBOUND','DIMENSIONS')
modWel['welCells'] = workFunctions.getCellsforBoundary(welLines,'wel',modWel['maxBound'],1)

# Get the HDS info
import flopy.utils.binaryfile as bf

for lay in range(modDis['cellLays']):

# <br/>
# <br/>
# <br/>
# ___
#
# # VTK file of Model Geometry, Model Results and Boundary Conditions
#

# ## Point Data

# In[6]:

### Water Tables Vextex
listWaterTableVertex = listFunctions.listWaterTableVertexFunction(modDis,modHds)

# ## Point Definition

# In[7]:

### Definition of XYZ points for All Vertex
vertexXYZPoints = listFunctions.vertexXYZPointsFunction(modDis)

### Definition of XYZ points for Water Table
vertexWaterTableXYZPoints = listFunctions.vertexWaterTableXYZPointsFunction(listWaterTableVertex,modDis)

# ## Quad and Hexa Sequences

# In[8]:

### List of Layer Quad Sequences (Works only for a single layer)

### List of Hexa Sequences for All Cells
listHexaSequence = listFunctions.listHexaSequenceFunction(modDis)

### List of Hexa Sequences for DRN Cells
listDrnCellsHexaSecuence = listFunctions.bcCellsListFunction(modDrn,'drnCells',listHexaSequence,modDis)[1]

### List of Hexa Sequences for CHD Cells
listChdCellsHexaSecuence = listFunctions.bcCellsListFunction(modChd,'chdCells',listHexaSequence,modDis)[1]

### List of Hexa Sequences for wEL Cells
listWelCellsHexaSecuence = listFunctions.bcCellsListFunction(modWel,'welCells',listHexaSequence,modDis)[1]

# ## Cell Data

# In[9]:

### Definition of DRN cells values '1' as List
listDrnCellsIO = listFunctions.bcCellsListFunction(modDrn,'drnCells',listHexaSequence,modDis)[0]

### Definition of CHD cells values '1' as List
listChdCellsIO = listFunctions.bcCellsListFunction(modChd,'chdCells',listHexaSequence,modDis)[0]

### Definition of WEL cells values '1' as List
listWelCellsIO = listFunctions.bcCellsListFunction(modWel,'welCells',listHexaSequence,modDis)[0]

### Water Tables on Cell
listWaterTableCell = listFunctions.listWaterTableCellFunction(modDis,modHds)

# <br/>
# <br/>
# <br/>
# ___
#
# # VTK Creation

# ## Heads on Vertex and Cells VTK

# In[10]:

vtkFunctions.printPointDefinition(vtkText,vertexXYZPoints)

vtkFunctions.printCellHexaConnectivityOffsetType(vtkText,listHexaSequence)

vtkFunctions.printFooter(vtkText)

vtkText.close()

# ## Water Table VTK

# In[11]:

vtkText = open('../vtuFiles/hatari01_WaterTable.vtu','w')

vtkFunctions.printCellData(vtkText,'WaterTableElev',listWaterTableCell)

vtkFunctions.printPointDefinition(vtkText,vertexWaterTableXYZPoints)

vtkFunctions.printFooter(vtkText)

vtkText.close()

# ## DRN Package VTK

# In[12]:

vtkText = open('../vtuFiles/hatari01_DRNCells.vtu','w')

vtkFunctions.printCellData(vtkText,'DRNCells',listDrnCellsIO)

vtkFunctions.printPointDefinition(vtkText,vertexXYZPoints)

vtkFunctions.printCellHexaConnectivityOffsetType(vtkText,listDrnCellsHexaSecuence)

vtkFunctions.printFooter(vtkText)

vtkText.close()

# ## CHD Package VTK

# In[13]:

vtkText = open('../vtuFiles/hatari01_CHDCells.vtu','w')

vtkFunctions.printCellData(vtkText,'CHDCells',listChdCellsIO)

vtkFunctions.printPointDefinition(vtkText,vertexXYZPoints)

vtkFunctions.printCellHexaConnectivityOffsetType(vtkText,listChdCellsHexaSecuence)

vtkFunctions.printFooter(vtkText)

vtkText.close()

# ## WEL Package VTK

# In[14]:

vtkText = open('../vtuFiles/hatari01_WELCells.vtu','w')

vtkFunctions.printCellData(vtkText,'WELCells',listWelCellsIO)

vtkFunctions.printPointDefinition(vtkText,vertexXYZPoints)

vtkFunctions.printCellHexaConnectivityOffsetType(vtkText,listWelCellsHexaSecuence)

vtkFunctions.printFooter(vtkText)

vtkText.close()```

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.

Subscribe to our free e-newsletter for tutorials, articles, webminars, courses and more.