Malla triangular para modelamiento de aguas subterráneas con MODFLOW6 y Flopy - Tutorial

TriangularMeshModelingMODFLOW6.png

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 Triangle
flopy is installed in E:\Software\Anaconda3\lib\site-packages\flopy

Definir 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>
output_9_1.png
### 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>
output_10_1.png

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

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>
output_28_1.png
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()
output_30_0.png

Archivos de entrada

Puedes descargar los archivos de entrada para este tutorial de este link.

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 April 8, 2019 and filed under TutorialModflow, TutorialPython.