Tutorial Básico de Modelamiento en MODFLOW con Python usando Flopy

Flopy es la librería creada por el Servicio Geológico de los Estados Unidos (USGS) para la creación, configuración y representación de resultados de modelos en MODFLOW. Flopy es una herramienta avanzada que tiene soporte incluso para la creación de grillas no estructuradas en MODFLOW 6. Con el uso conjunto de Python y MODFLOW mediante Flopy se extienden las posibilidades de modelamiento y de gestión de aguas subterráneas al permitir la configuración de nuevos esquemas de optimización y la aplicación de algoritmos de inteligencia artificial.

Este tutorial muestra el proceso completo de construcción de un modelo numérico en condiciones transientes con Flopy. El modelo tiene un pozo con caudal variable y un flujo regional determinado por la gradiente en las cargas hidráulicas. El código en Python ha sido realizado de manera interactiva en Jupyter Notebook.

Usted puede acceder a un tutorial para instalar Flopy en Anaconda 3 en este enlace.

Video

Código

Este es el código completo en Python:

import flopy, os
import numpy as np

modelname = "EjemploFlopy"
workspace = "ModelFiles"
mf = flopy.modflow.Modflow(modelname, exe_name='mf2005', model_ws="ModelFiles")

# Model domain and grid definition
Lx = 1000.
Ly = 1000.
ztop = 100.
zbot = -100.
#SpatialDiscretization
nlay = 4
nrow = 10
ncol = 10
delr = Lx/ncol
delc = Ly/nrow
delv = (ztop - zbot) / nlay
botm = np.linspace(ztop, zbot, nlay + 1)
#TemporalDiscretization
nper = 4
perlen = [2.0, 3.0, 4.0, 3.0]
perlen = [x * 86400 for x in perlen]
nstp = [3,1,2,1]
steady=False
#Unitsystem
itmuni=1
lenuni=2

# Create the discretization object
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc,
                               top=ztop, botm=botm[1:], 
                               nper=nper, perlen=perlen, nstp=nstp, steady=False,
                               itmuni=1, lenuni=2)

# Variables for the BAS package
ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
ibound[0, :, 0] = -1
ibound[0, :, -1] = -1
strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
strt[0, :, 0] = 100.
strt[0, :, -1] = 80.
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)

#Configuration of well package
sp_data = {0:[2, 6, 6, -0.01], 1:[2, 6, 6, -0.025], 2:[2, 6, 6, -0.036], 3:[2, 6, 6, 0]}
wel1 = flopy.modflow.ModflowWel(mf, stress_period_data=sp_data)

# Add LPF package to the MODFLOW model
lpf = flopy.modflow.ModflowLpf(mf, hk=2e-5, vka=2e-5)

# Add OC package to the MODFLOW model
stressperioddata={(0,2):['save head'],(1,0):['save head'],(2,1):['save head'],(3,0):['save head']}
oc = flopy.modflow.ModflowOc(mf, stress_period_data=stressperioddata)

# Add PCG package to the MODFLOW model
pcg = flopy.modflow.ModflowPcg(mf)

# Write the MODFLOW model input files
mf.write_input()

# Run the MODFLOW model
success, buff = mf.run_model()

import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf

totaltimes = np.cumsum(perlen)

plt.figure(figsize=(12,8))

plt.subplot(1,2,1,aspect='equal')
hds = bf.HeadFile(os.path.join(workspace,modelname+'.hds'))
head = hds.get_data(totim=totaltimes[2])
levels = np.arange(0,100,1)
extent = (delr/2., Lx - delr/2., Ly - delc/2., delc/2.)
plt.contour(head[0, :, :], levels=levels, extent=extent)

plt.subplot(1,2,2,aspect='equal')
hds = bf.HeadFile(os.path.join(workspace,modelname+'.hds'))
head = hds.get_data(totim=totaltimes[2])
levels = np.arange(0,100,1)
extent = (delr/2., Lx - delr/2., Ly - delc/2., delc/2.)
plt.contour(head[1, :, :], levels=levels, extent=extent)

plt.show()

import shutil
shutil.copyfile(os.path.join(workspace,modelname+'.hds'),os.path.join(workspace,modelname+'.bhd'))