MODFLOW 6 es la última versión del código de modelamiento de agua subterránea MODFLOW desarrollado por la USGS. Esta versión incluye herramientas innovadoras y un reestructuramiento completo de la construcción del modelo, pero hasta la fecha no tienen ninguna Interfaz Gráfica de Usuario (GUI de sus siglas en inglés) que sea comercial o de código abierto que permita mejorar la manipulación del código para modeladores de agua subterránea principiantes o intermedios.
Una de las novedades más excepcionales de MODFLOW 6 son las diferentes opciones de discretización para la generación de la grilla del modelo. Las opciones van desde grillas regulares (igual que MODFLOW 2005), grillas triangulares y grillas no estructuradas. Flopy que es una librería de Python para la construcción y simulación de MODFLOW 6 y otros modelos tiene herramientas para la generación de grillas triangulares. El flujo de trabajo en modelamiento de aguas subterráneas con MODFLOW 6 y Flopy para modelos con grillas triangulares es muy fluido y observamos mucho potencial para el modelamiento de flujo de agua subterránea en escala local y regional.
Este tutorial muestra el procedimiento completo para crear una grilla triangular con las utilidades de Flopy e incorporarlo al modelo de MODFLOW 6. El modelo es simulado y los resultados son representados como una grilla colorida e isolíneas.
Tutorial
Código de Python
Este es el script en Python para el tutorial:
Importar librerías requeridas
Este tutorial requiere algunas librerías de Python como Numpy, Matplotlib y Flopy.
import os
import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
import flopy
from flopy.utils.triangle import Triangle as Triangleflopy is installed in E:\Software\Anaconda3\lib\site-packages\flopyDefinir la el folder del modelo y el ejecutable
Se define el folder para los archivos del modelo y los resultados de la simulación, así como también los ejecutables para el generador de la grilla triangular y Modflow 6.
workspace = '../Model/'
triExeName = '../Exe/triangle.exe'
mf6ExeName = '../Exe/mf6.exe'Definir el dominio y la ubicación del dren
El dominio activo del modelo es un trapezoide con el lado paralelo más pequeño localizado al lado este. El dren fue conceptualizado como un rectángulo en la parte central izquierda del modelo.
active_domain = [(0, 0), (100, 10), (100, 55), (0, 65)]
drain_domain= [(80, 31), (95, 31), (95, 34), (80, 34)]Creación de la grilla triangular
Usamos la herramienta Triángulo de las utilidades de Flopy que implementan el programa Triángulo para generar grillas. Este tutorial viene con el ejecutable de Triángulo para Windows, si tu usas otro sistema operativo es necesario compilar el código fuente que puedes encontrar de https://www.cs.cmu.edu/~quake/triangle.html.
tri = Triangle(angle=30, model_ws=workspace, exe_name=triExeName)
tri.add_polygon(active_domain)
tri.add_polygon(drain_domain)
tri.add_region((10,10),0,maximum_area=30) #coarse discretization
tri.add_region((88,33),1,maximum_area=5) #fine discretization
tri.build()Grilla triangular generada y representación de límites
El siguiente paso es la representación de la grilla e identificación de los límites creados de la generación de la grilla.
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(edgecolor='gray')<matplotlib.collections.PatchCollection at 0x9727518>
### Indentification and representation of grid boundaries
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(edgecolor='gray')
numberBoundaries = tri.get_boundary_marker_array().max()+1
cmap = plt.cm.get_cmap("hsv", numberBoundaries)
labelList = list(range(1,numberBoundaries))
i = 0
for ibm in range(1,numberBoundaries):
tri.plot_boundary(ibm=ibm, ax=ax,marker='o', color=cmap(ibm), label= ibm)
handles, labels = plt.gca().get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys(), loc='lower right')
<matplotlib.legend.Legend at 0x9c45978>
Representación de las características triangulares
Usando los métodos disponibles para el objeto triángulo “tri” podemos representar el vértice y el centroide con su posición y etiqueta.
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(ax=ax, edgecolor='gray')
tri.plot_vertices(ax=ax, marker='o', color='blue')
tri.label_vertices(ax=ax, fontsize=10, color='blue')
tri.plot_centroids(ax=ax, marker='o', color='red')
tri.label_cells(ax=ax, fontsize=10, color='red')
Configuración de un modelo MODFLOW 6
Una vez tengamos la grilla triangular podemos crear un modelo de MODFLOW 6. Cómo pueden saber, MODFLOW 6 tiene una diferencia entre un modelo y una simulación desde que este implementa el concepto de intercambio entre modelos, para este caso tenemos una simulación de un solo modelo.
name = 'mf'
sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
exe_name=mf6ExeName,
sim_ws=workspace)
tdis = flopy.mf6.ModflowTdis(sim, time_units='SECONDS',
perioddata=[[1.0, 1, 1.]])
gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
ims = flopy.mf6.ModflowIms(sim, print_option='SUMMARY', complexity='complex',
outer_hclose=1.e-8, inner_hclose=1.e-8)Configuración del paquete de discretización triangular DISV
MODFLOW 6 diferencia el archivo de discretización espacial y el archivo de discretización temporal. En las opciones de discretización espacial tenemos tres tipos de grillas: regular (.dis), triangular (.disv) y no estructurada (.disu). Este caso de estudio usa una discretización triangular.
cell2d = tri.get_cell2d()
vertices = tri.get_vertices()
xcyc = tri.get_xcyc()
nlay = 1
ncpl = tri.ncpl
nvert = tri.nvert
top = 1.
botm = [0.]
dis = flopy.mf6.ModflowGwfdisv(gwf, length_units='METERS',
nlay=nlay, ncpl=ncpl, nvert=nvert,
top=top, botm=botm,
vertices=vertices, cell2d=cell2d)Configuración del paquete de flujo y condiciones iniciales
npf = flopy.mf6.ModflowGwfnpf(gwf, k=0.0001)
ic = flopy.mf6.ModflowGwfic(gwf, strt=10.0)Definición de Cargas Constantes (CHD) y condiciones de borde
El lado izquierdo y derecho del modelo fueron definidos como una condición de borde CHD para representar un flujo regional de derecha a izquierda desde una carga hidráulica de 30m a 10m. Observar que los bordes del vértices son obtenidos de los bordes de la grilla del modelo.
chdList = []
leftcells = tri.get_edge_cells(4)
rightcells = tri.get_edge_cells(2)
for icpl in leftcells: chdList.append([(0, icpl), 30])
for icpl in rightcells: chdList.append([(0, icpl), 10])
chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdList)Definición de la condición de borde de Dren DRN
En la parte izquierda central del modelo hay un dren conceptualizado como 2 líneas en el borde de una región declarada previamente. El dren tiene una elevación de 2.5 menor que la condición de carga constante oriental.
drnList = []
drnCells = tri.get_edge_cells(5) + tri.get_edge_cells(7)
for icpl in drnCells: drnList.append([(0, icpl), 7.5, 0.10])
chd = flopy.mf6.ModflowGwfdrn(gwf, stress_period_data=drnList)Definir las opciones de salida, escribir archivos de entrada y correr el modelo.
El tutorial define el tipo de dato que será almacenado y los registros en el archivo list. Luego los archivos del modelo MODFLOW 6 son escritos en el folder del modelo y la simulación es desarrollada.
oc = flopy.mf6.ModflowGwfoc(gwf,
budget_filerecord='{}.cbc'.format(name),
head_filerecord='{}.hds'.format(name),
saverecord=[('HEAD', 'LAST'),
('BUDGET', 'LAST')],
printrecord=[('HEAD', 'LAST'),
('BUDGET', 'LAST')])
sim.write_simulation()
success, buff = sim.run_simulation()
writing simulation...
writing simulation name file...
writing simulation tdis package...
writing ims package ims_-1...
writing model mf...
writing model name file...
writing package disv...
writing package npf...
writing package ic...
writing package chd_0...
INFORMATION: maxbound in ('gwf6', 'chd', 'dimensions') changed to 19 based on size of stress_period_data
writing package drn_0...
INFORMATION: maxbound in ('gwf6', 'drn', 'dimensions') changed to 32 based on size of stress_period_data
writing package oc...
FloPy is using the following executable to run the model: ../Exe/mf6.exe
MODFLOW 6
U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
VERSION 6.0.3 08/09/2018
MODFLOW 6 compiled Aug 09 2018 13:40:32 with IFORT compiler (ver. 18.0.3)
This software has been approved for release by the U.S. Geological
Survey (USGS). Although the software has been subjected to rigorous
review, the USGS reserves the right to update the software as needed
pursuant to further analysis and review. No warranty, expressed or
implied, is made by the USGS or the U.S. Government as to the
functionality of the software and related material nor shall the
fact of release constitute any such warranty. Furthermore, the
software is released on condition that neither the USGS nor the U.S.
Government shall be held liable for any damages resulting from its
authorized or unauthorized use. Also refer to the USGS Water
Resources Software User Rights Notice for complete use, copyright,
and distribution information.
Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/04/08 16:10:56
Writing simulation list file: mfsim.lst
Using Simulation name file: mfsim.nam
Solving: Stress period: 1 Time step: 1
Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/04/08 16:10:56
Elapsed run time: 0.071 Seconds
Normal termination of simulation.Output data representation
Importamos las cargas calculadas y las representamos en la extensión del modelo. Una representación de las cargas equipotenciales es realizada por la interpolación de las cargas en cada centroide de cada celda.
fname = os.path.join(workspace, name + '.hds')
hdobj = flopy.utils.HeadFile(fname, precision='double')
head = hdobj.get_data()
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
img=tri.plot(ax=ax, a=head[0, 0, :], cmap='Spectral')
fig.colorbar(img, fraction=0.02)
<matplotlib.colorbar.Colorbar at 0xbe44898>
from scipy.interpolate import griddata
x, y = tri.get_xcyc()[:,0],tri.get_xcyc()[:,1]
xG = np.linspace(x.min(),x.max(),100)
yG = np.linspace(y.min(),y.max(),100)
X, Y = np.meshgrid(xG, yG)
z = head[0,0,:]
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
img=tri.plot(ax=ax, edgecolor='lightsteelblue', a=head[0, 0, :], cmap='Spectral', alpha=0.8)
fig.colorbar(img, fraction=0.02)
Ti = griddata((x, y), z, (X, Y), method='cubic')
CS = ax.contour(X,Y,Ti, colors='Cyan', linewidths=2.0)
ax.clabel(CS, inline=1, fontsize=15, fmt='%1.0f')
plt.show()
Archivos de entrada
Puedes descargar los archivos de entrada para este tutorial de este link.
