Tutorial de Modelamiento Regional de Aguas Subterráneas con MODFLOW y Flopy

RegionalModelingMODFLOWandFlopyTutorial.jpg

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 plt

Crear 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)*3000

El 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=1




Escribir 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)
output_25_0.png

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')
output_27_0.png

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)
output_29_0.png

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')
output_31_0.png

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)
output_33_0.png

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)
output_35_0.png

Datos de entrada

Usted puedes descargar los datos de entrada de este enlace.



Smiley face

Suscríbase a nuestro boletín gratuito para recibir noticias, datos interesantes y fechas de nuestros cursos en recursos hídricos.

Posted on February 28, 2019 and filed under TutorialModflow.