El modelamiento regional de aguas subterráneas es una tarea importante en la gestión estratégica del agua que involucra a todos los usuarios, actividades y ecosistemas involucrados y proporciona un uso sostenible para las condiciones actuales y futuras. Existen algunas consideraciones específicas sobre el modelamiento regional con respecto a la línea base y la discretización espacial; un modelo regional no pretende proporcionar la respuesta del acuífero para un área determinada, sino busca la evaluación del flujo regional del agua subterránea y la cuantificación de la recarga, descarga y otros procesos del balance hídrico.
Este tutorial es la versión de Flopy de un ejemplo numérico de la cuenca de Angascancha realizada con MODFLOW y Model Muse que se puede encontrar en este enlace. El ejemplo está en régimen uniforme y se resuelve con el solucionador NWT. Las representaciones de salida del modelo se han realizado con las herramientas Flopy / Matplotlib, así como también con algunos códigos Python para crear archivos VTU y representarlas en Paraview.
Video
Código de Python
Este es el código de Python descrito para este tutorial:
Importar las bibliotecas requeridas
Este tutorial requiere algunas bibliotecas centrales de Python, bibliotecas Scipy (Numpy y Matplotlib), así como la biblioteca Flopy. Para la administración y extracción de datos de archivos raster usamos la biblioteca GDAL.
import flopy
import os
import flopy.utils.binaryfile as bf
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as pltCrear un objeto de modelo MODFLOW
Así como para cualquier modelo generado por Flopy, primero tenemos que configurar el nombre del modelo y el directorio de trabajo. Recomendamos seguir la estructura de las carpetas / archivos que se proporciona en la parte de "Datos de entrada" de este tutorial. Debido a la geometría compleja y a las diferencias de elevaciones y espesores de la capas, se aplicó el solucionador NWT.
modelname = "model1"
modelpath = "../Model/"
#Modflow Nwt
mf1 = flopy.modflow.Modflow(modelname, exe_name="../Exe/MODFLOW-NWT_64.exe", version="mfnwt",model_ws=modelpath)
nwt = flopy.modflow.ModflowNwt(mf1 ,maxiterout=150,linmeth=2)Abrir y leer archivos raster
En esta parte utilizaremos la biblioteca GDAL para abrir los rásters de elevación y las condiciones de contorno y obtener sus parámetros de geotransformación. El tutorial requiere que la resolución ráster sea la misma que la resolución de la grilla. Las bandas del raster de elevación se lee como matriz Numpy y se inserta en la parte superior del modelo.
#Raster paths
demPath ="../Rst/DEM_200.tif"
crPath = "../Rst/CR.tif"
#Open files
demDs =gdal.Open(demPath)
crDs = gdal.Open(crPath)
geot = crDs.GetGeoTransform() #Xmin, deltax, ?, ymax, ?, delta y
geot(613882.5068, 200.0, 0.0, 8375404.9524, 0.0, -200.0)# Get data as arrays
demData = demDs.GetRasterBand(1).ReadAsArray()
crData = crDs.GetRasterBand(1).ReadAsArray()
# Get statistics
stats = demDs.GetRasterBand(1).GetStatistics(0,1)
stats[3396.0, 4873.0, 4243.0958193041, 287.90182492684]Discretización espacial
Dado que este tutorial se ocupa de una escala regional, las diferencias de elevación en el área de estudio pueden ser de más de cien metros o más de mil metros si tratamos con cuencas andinas. La conceptualización hidrogeológica para este tutorial define dos capas delgadas superiores (Layer1 y Layer1_2) y capas más gruesas debajo.
#Elevation of layers
Layer1 = demData - 7.5
Layer1_2 = demData - 15
Layer2 = (demData - 3000)*0.8 + 3000
Layer3 = (demData - 3000)*0.5 + 3000
Layer4 = (demData>0)*3000El número y dimensión de columnas y filas provienen de los atributos ráster.
#Boundaries for Dis = Create discretization object, spatial/temporal discretization
ztop = demData
zbot = [Layer1, Layer1_2, Layer2, Layer3, Layer4]
nlay = 5
nrow = demDs.RasterYSize
ncol = demDs.RasterXSize
delr = geot[1]
delc = abs(geot[5])Definición de los paquetes de flujo para el modelo MODFLOW
En esta parte definimos los paquetes relacionados con el flujo de agua subterránea en el modelo MODFLOW. Primero definimos el paquete DIS que tiene la geometría y el tipo de simulación (uniforme / transitorio). El modelo funciona en condiciones de flujo uniforme.
dis = flopy.modflow.ModflowDis(mf1, nlay,nrow,ncol,delr=delr,delc=delc,top=ztop,botm=zbot,itmuni=1)Luego definimos otros paquetes del modelo. Algunas condiciones de contorno se configuran a partir de los valores ráster mientras que otros paquetes utilizan un valor constante. Estos son los paquetes requeridos para la simulación:
- BAS para configurar las celdas constantes y las celdas inactivas,
- UPW que define la conductividad hidráulica vertical / horizontal,
- RCH que aplica recarga de precipitación.
- EVT para la simulación del consumo de agua subterránea desde las raíces.
- DRN para la descarga de agua subterránea a la red de canales.
- OC para el registro de salida
# Variables for the BAS package
iboundData = np.zeros(demData.shape, dtype=np.int32)
iboundData[crData > 0 ] = 1
bas = flopy.modflow.ModflowBas(mf1,ibound=iboundData,strt=ztop, hnoflo=-2.0E+020)# Array of hydraulic heads per layer
hk = [1E-4, 1E-5, 1E-7, 1E-8, 1E-9]
# Add UPW package to the MODFLOW model
upw = flopy.modflow.ModflowUpw(mf1, laytyp = [1,1,1,1,0], hk = hk)#Add the recharge package (RCH) to the MODFLOW model
rch = np.ones((nrow, ncol), dtype=np.float32) * 0.120/365/86400
rch_data = {0: rch}
rch = flopy.modflow.ModflowRch(mf1, nrchop=3, rech =rch_data)#Add the evapotranspiration package (EVT) to the MODFLOW model
evtr = np.ones((nrow, ncol), dtype=np.float32) * 1/365/86400
evtr_data = {0: evtr}
evt = flopy.modflow.ModflowEvt(mf1,nevtop=1,surf=ztop,evtr=evtr_data, exdp=0.5)#Add the drain package (DRN) to the MODFLOW model
river = np.zeros(demData.shape, dtype=np.int32)
river[crData == 3 ] = 1
list = []
for i in range(river.shape[0]):
for q in range(river.shape[1]):
if river[i,q] == 1:
list.append([0,i,q,ztop[i,q],0.001]) #layer,row,column,elevation(float),conductanceee
rivDrn = {0:list}
drn = flopy.modflow.ModflowDrn(mf1, stress_period_data=rivDrn)# Add OC package to the MODFLOW model
oc = flopy.modflow.ModflowOc(mf1) #ihedfm= 1, iddnfm=1Escribir los archivos del modelo MODFLOW y ejecutar la simulación.
Esta parte del script escribe el archivo .nam y los archivos declarados en el archivo .nam. Luego se ejecuta la simulación y se muestra la información de la simulación en la pantalla.
#Write input files -> write file with extensions
mf1.write_input()
#run model -> gives the solution
mf1.run_model()FloPy is using the following executable to run the model: ../Exe/MODFLOW-NWT_64.exe
MODFLOW-NWT-SWR1
U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUNDWATER-FLOW MODEL
WITH NEWTON FORMULATION
Version 1.1.4 4/01/2018
BASED ON MODFLOW-2005 Version 1.12.0 02/03/2017
SWR1 Version 1.04.0 09/15/2016
Using NAME file: model1.nam
Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/02/28 9:28:21
Solving: Stress period: 1 Time step: 1 Groundwater-Flow Eqn.
Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/02/28 9:29:01
Elapsed run time: 39.992 Seconds
Normal termination of simulation
(True, [])Model results post-processing
Este tutorial tiene algunas representaciones de las características del modelo y los datos de salida realizadas con Flopy y Matplotlib. La representación 3D es en el formato VTK que se realiza en otro tutorial.
#Import model
ml = flopy.modflow.Modflow.load('../Model/model1.nam')Representación de la grilla del modelo
# First step is to set up the plot
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
# Next we create an instance of the ModelMap class
modelmap = flopy.plot.ModelMap(sr=ml.dis.sr)
# Then we can use the plot_grid() method to draw the grid
# The return value for this function is a matplotlib LineCollection object,
# which could be manipulated (or used) later if necessary.
linecollection = modelmap.plot_grid(linewidth=0.4)
Sección transversal de la grilla del modelo
fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1)
# Next we create an instance of the ModelCrossSection class
#modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Column': 5})
modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Row': 99})
# Then we can use the plot_grid() method to draw the grid
# The return value for this function is a matplotlib LineCollection object,
# which could be manipulated (or used) later if necessary.
linecollection = modelxsect.plot_grid(linewidth=0.4)
t = ax.set_title('Column 6 Cross-Section - Model Grid')
Celdas activas / inactivas en la extensión del modelo
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=0)
quadmesh = modelmap.plot_ibound(color_noflow='cyan')
linecollection = modelmap.plot_grid(linewidth=0.4)
Secciones transversales de celdas activas / inactivas
fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1)
modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Column': 5})
patches = modelxsect.plot_ibound(color_noflow='cyan')
linecollection = modelxsect.plot_grid(linewidth=0.4)
t = ax.set_title('Column 6 Cross-Section with IBOUND Boundary Conditions')
Red de canales como paquete de drenaje (DRN)
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=-14)
quadmesh = modelmap.plot_ibound(color_noflow='cyan')
quadmesh = modelmap.plot_bc('DRN', color='blue')
linecollection = modelmap.plot_grid(linewidth=0.4)
Representación de grilla y cargas hidráulicas
fname = os.path.join(modelpath, 'model1.hds')
hdobj = flopy.utils.HeadFile(fname)
head = hdobj.get_data()
fig = plt.figure(figsize=(30, 30))
ax = fig.add_subplot(1, 2, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=-14)
#quadmesh = modelmap.plot_ibound()
quadmesh = modelmap.plot_array(head, masked_values=[-2.e+20], alpha=0.8)
linecollection = modelmap.plot_grid(linewidth=0.2)
