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.
Photo Gallery
Tutorial
Input data
You can download the input data from the following link.
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
chdLines = open('../model/hatari01.chd').readlines()
disLines = open('../model/hatari01.dis').readlines()
drnLines = open('../model/hatari01.drn').readlines()
icLines = open('../model/hatari01.ic').readlines()
npfLines = open('../model/hatari01.npf').readlines()
rchLines = open('../model/hatari01.rch').readlines()
welLines = open('../model/hatari01.wel').readlines()
# 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
### Store heads per lay
import flopy.utils.binaryfile as bf
modHds['cellHeadGrid'] = {}
headObject = bf.HeadFile('..\\model\\hatari01.hds', precision='double')
headObjectList = headObject.get_data()
headObject.close()
for lay in range(modDis['cellLays']):
modHds['cellHeadGrid']['lay'+str(lay)] = headObjectList[lay]
vertexHeadGridCentroid = transFunctions.vertexHeadGridCentroidFunction(modDis,modHds)
modHds['vertexHeadGrid'] = transFunctions.vertexHeadGridFunction(vertexHeadGridCentroid,modDis,modHds)
# <br/>
# <br/>
# <br/>
# ___
#
# # VTK file of Model Geometry, Model Results and Boundary Conditions
#
# ## Point Data
# In[6]:
### Vertex Heads
listVertexHead = listFunctions.listCellHeadsFunction('vertexLays','vertexHeadGrid',modDis,modHds)
### 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)
listLayerQuadSequence = listFunctions.listLayerQuadSequenceFunction(modDis)
### 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 cellHead
listCellHead = listFunctions.listCellHeadsFunction('cellLays','cellHeadGrid',modDis,modHds)
### 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]:
vtkText = open('../vtuFiles/hatari01_Heads.vtu','w')
vtkFunctions.printHeader(vtkText,len(vertexXYZPoints),len(listHexaSequence))
vtkFunctions.printPointData(vtkText,'VertexHeads',listVertexHead)
vtkFunctions.printCellData(vtkText,'CellHeads',listCellHead)
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.printHeader(vtkText,len(vertexWaterTableXYZPoints),len(listWaterTableCell))
vtkFunctions.printCellData(vtkText,'WaterTableElev',listWaterTableCell)
vtkFunctions.printPointDefinition(vtkText,vertexWaterTableXYZPoints)
vtkFunctions.printCellQuadConnectivityOffsetType(vtkText,listLayerQuadSequence)
vtkFunctions.printFooter(vtkText)
vtkText.close()
# ## DRN Package VTK
# In[12]:
vtkText = open('../vtuFiles/hatari01_DRNCells.vtu','w')
vtkFunctions.printHeader(vtkText,len(vertexXYZPoints),len(listDrnCellsHexaSecuence))
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.printHeader(vtkText,len(vertexXYZPoints),len(listChdCellsHexaSecuence))
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.printHeader(vtkText,len(vertexXYZPoints),len(listWelCellsHexaSecuence))
vtkFunctions.printCellData(vtkText,'WELCells',listWelCellsIO)
vtkFunctions.printPointDefinition(vtkText,vertexXYZPoints)
vtkFunctions.printCellHexaConnectivityOffsetType(vtkText,listWelCellsHexaSecuence)
vtkFunctions.printFooter(vtkText)
vtkText.close()
