El desarrollo de hidrogeológicos puede llevar mucho tiempo debido a todos los pasos involucrados como la construcción del modelo, la calibración y la visualización de salida. Es importante utilizar herramientas que puedan optimizar estas tareas y permitir que el tiempo ahorrado se utilice en el análisis del sistema de flujo y el estado de la calidad del agua subterránea.

Flopy es un conjunto versátil de scripts de Python que se pueden usar para ejecutar MODFLOW y MT3D, entre otros programas de agua subterránea relacionados con MODFLOW de una manera simple y eficiente. Este tutorial desarrollará un caso completo aplicado de modelamiento de flujo y transporte con MODFLOW, MT3D-USDS y Flopy. El caso de estudio describe un flujo de agua subterránea regional que tiene una fuente puntual con un esquema de remediación simulado en condiciones de estado estacionario para condiciones de flujo y transitorio para el modelamiento de transporte.

Contenido

El tutorial tiene el siguiente contenido:

Modelo de flujo:

Configuración de un modelo MODFLOW con Flopy
Establecimiento de discretización espacial y temporal.
Definición de parámetros hidráulicos.
Modelo de ejecución y representación de resultados

Modelo de transporte:

Configuración de un modelo MT3D-USGS con Flopy
Implementación de advección, dispersión y fuentes de masa.
Ejecución de transporte y simulación de salida

Tutorial

Código

Este es el código en Flopy para la simulación de transporte de contaminantes:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import flopy
import flopy.modflow as mf
import flopy.mt3d as mt
import flopy.utils as fu

flopy is installed in E:\Software\Anaconda3\lib\site-packages\flopy

#MODFLOW 2005
modelname = 'flowModel'
modPath= '../Model/'
mfModel = mf.Modflow(modelname = modelname, model_ws=modPath, exe_name="../Exe/MODFLOW-NWT_64.exe")

#DIS file
Lx = 610
Ly = 310
nrow = 31
ncol = 61
nlay = 1
delr = Lx / ncol
delc = Ly / nrow
top = np.ones((nrow, ncol))*20
botm = np.ones((nrow, ncol))*-40
nper = 1
perlen = 2700
nstp = 5

dis = mf.ModflowDis(mfModel, nlay, nrow, ncol, delr = delr, delc = delc, 
                    top = top, botm = botm, laycbd = 0, itmuni=4, 
                    nper = nper, perlen = perlen, nstp = nstp, tsmult=1.4)

# Output Control: Create a flopy output control object
oc = mf.ModflowOc(mfModel,stress_period_data={(nper-1, nstp-1): ['save head']},)

#BCF file
laycon=1 #confined
tran=1600.0 #transmissivity
bcf = flopy.modflow.mfbcf.ModflowBcf(mfModel,laycon=laycon, tran=tran)

#BAS file
strt=15 #starting head
ibound=np.ones((nrow, ncol))
bas = mf.ModflowBas(mfModel, ibound = ibound, strt = strt)

#PCG file
pcg = flopy.modflow.mfpcg.ModflowPcg(mfModel, mxiter=20, iter1=30, hclose=1e-03, rclose=1e-03, relax=1.0)

#CHD
chd_data = []
for c in range(mfModel.dis.nrow):
    lChd = np.array([0, c, 0, 20, 20])
    rChd = np.array([0, c, mfModel.dis.ncol-1, 12, 15])
    chd_data.append(lChd)
    chd_data.append(rChd)
stress_period_data = 

stress_period_data[0][:5]

[array([ 0,  0,  0, 20, 20]),
 array([ 0,  0, 60, 12, 15]),
 array([ 0,  1,  0, 20, 20]),
 array([ 0,  1, 60, 12, 15]),
 array([ 0,  2,  0, 20, 20])]

chd = mf.mfchd.ModflowChd(mfModel, stress_period_data=stress_period_data)

#WELL
inyectingWell = 200 #m3/d
pumpingWell = -400
wel_sp1 = []
wel_sp1.append([0, 15, 15, inyectingWell])
wel_sp1.append([0, 5, 45, pumpingWell])
stress_period_data = {0: wel_sp1}

wel = flopy.modflow.ModflowWel(mfModel, stress_period_data=stress_period_data)

#LMT Linkage with MT3DMS for multi-species mass transport modeling
lmt = flopy.modflow.ModflowLmt(mfModel, output_file_name='mt3d_link.ftl')

#Write input files
mfModel.write_input()

# run the model
mfModel.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: flowModel.nam 
 Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/11/14 16:17:46

 Solving:  Stress period:     1    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     1    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     1    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     1    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     1    Time step:     5    Groundwater-Flow Eqn.
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/11/14 16:17:46
 Elapsed run time:  0.048 Seconds

  Normal termination of simulation





(True, [])

#Plot model results
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf

# Create the headfile object
headobj = bf.HeadFile(modPath + modelname+'.hds')
times = headobj.get_times()
print(times)
head = headobj.get_data(totim=times[-1])


# Setup contour parameters
levels = np.arange(10, 20, 1)
extent = (delr/2., Lx - delr/2., delc/2., Ly - delc/2.)

# Make the plots
fig = plt.figure(figsize=(18,8))
plt.subplot(1, 1, 1, aspect='equal')
plt.title('Head distribution (m)')
plt.imshow(head[0, :, :], extent=extent, cmap='YlGnBu', vmin=15., vmax=20.)
plt.colorbar()

contours = plt.contour(np.flipud(head[0, :, :]), levels=levels, extent=extent, zorder=10)
plt.clabel(contours, inline=1, fontsize=10, fmt='%d')# zorder=11)

plt.show()

[2699.9998]

#MT3D-USGS
namemt3d='transModel'
mt_model = mt.Mt3dms(modelname=namemt3d,  model_ws=modPath,version='mt3d-usgs', 
                     exe_name='../Exe/mt3d-usgs_1.1.0_64.exe', modflowmodel=mfModel)

#BTN file
btn = flopy.mt3d.Mt3dBtn(mt_model, sconc=0.0, prsity=0.3, thkmin=0.01, munit='g')#, icbund=icbund)

#ADV file
mixelm = -1 #Third-order TVD scheme (ULTIMATE)
percel = 1 #Courant number PERCEL is also a stability constraint
adv = flopy.mt3d.Mt3dAdv(mt_model, mixelm=mixelm, percel=percel)

#GCG file
mxiter = 1 #Maximum number of outer iterations
iter1 = 200 #Maximum number of inner iterations
isolve = 3 #Preconditioner = Modified Incomplete Cholesky
gcg = flopy.mt3d.Mt3dGcg(mt_model, mxiter=mxiter, iter1=iter1, isolve=isolve)

#DSP file
al = 10 #longitudinal dispersivity
dmcoef = 0 #effective molecular diffusion coefficient
trpt = 0.1 #ratio of the horizontal transverse dispersivity to the longitudinal dispersivity
trpv = 0.01 #ratio of the vertical transverse dispersivity to the longitudinal dispersivity

dsp = mt.Mt3dDsp(mt_model, al=al, dmcoef=dmcoef, trpt=trpt, trpv=trpv)

#SSM file
itype = flopy.mt3d.Mt3dSsm.itype_dict()
itype

{'CHD': 1,
 'BAS6': 1,
 'PBC': 1,
 'WEL': 2,
 'DRN': 3,
 'RIV': 4,
 'GHB': 5,
 'MAS': 15,
 'CC': -1}

#[K,I,J,CSS,iSSType] = layer, row, column, source concentration, type of sink/source: well-constant concentration cell 
ssm_data = {}
ssm_data[0] = [(0, 15, 15, 65.0, 2)]

ssm = flopy.mt3d.Mt3dSsm(mt_model, stress_period_data=ssm_data)

#Write model input
mt_model.write_input()

#Run the model
mt_model.run_model(silent=True)

(False, [])

#Plot concentration results
conc = fu.UcnFile(modPath+'MT3D001.UCN')
conc.get_times()

[2699.9998]

fig = plt.figure(figsize=(18,8))
conc.plot(totim=times[-1], colorbar='Concentration (mg/l)', cmap='Blues')
plt.title('Concentration distribution (mg/l)')
plt.show()