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
Uploaded by None on 2018-01-31.
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()
