Tutorial básico para aprender el procedimiento de construcción, simulación y representación de un modelo hidrogeológico en MODFLOW 6. El tutorial muestra una introducción al sistema de ficheros de un modelo en condiciones de régimen uniforme. El modelo de este tutorial tiene implementado las siguientes condiciones de borde: Drenes, Recarga, y Carga Constante. La grilla del modelo es regular con un ancho de 50 metros y tiene 30 filas con 24 columnas; el modelo consta de 4 capas y un espesor total de 130 metros.
El modelo es llamado "hatari01" y esta inspirado en el modelo "twri" de la documentación de MODFLOW 2005 adaptado a MODFLOW 6. Dentro de la conceptualización hidrogeológica se especifica la conductividad hidráulica en el plano horizontal y vertical. Luego de la simulación se ejecuta un script de Python dentro de un Jupyter Notebook para crear archivos VTK de grilla no estructurada para las cargas hidráulicas, napa freática y condiciones de borde como objetos 3D en Paraview.
Galería de Fotos
Tutorial
Datos de Ingreso
Ustedes pueden descargar los datos de ingreso del siguiente enlace.
Código
Esta es una parte del código en Python para crear los archivos VTK, el resto de código lo puede ubica en la parte de datos de ingreso de este 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()
