A continuación ponemos los códigos y simulaciones con datos iniciales.
A través de las páginas de Jupyter, utilizando la biblioteca FloPy, podemos realizar modelos hidrogeológicos de manera más eficiente y visual. FloPy, una interfaz de Python para MODFLOW y otros modelos de flujo de agua subterránea, permite a los usuarios construir, ejecutar y analizar modelos de manera interactiva dentro del entorno de Jupyter Notebook. Esto facilita la manipulación de datos, la visualización de resultados y la modificación de parámetros del modelo en tiempo real, lo cual es crucial para el desarrollo iterativo de modelos precisos y confiables. La capacidad de ejecutar bloques de código de manera secuencial y de integrar gráficos y tablas directamente en las celdas del notebook proporciona una plataforma robusta para la enseñanza y la investigación en hidrogeología.
Además, Jupyter Notebooks ofrece un entorno colaborativo donde los resultados del modelado pueden ser fácilmente compartidos y reproducidos. La visualización dinámica, mediante gráficos interactivos y mapas, mejora significativamente la interpretación de los resultados del modelo, permitiendo a los usuarios identificar patrones y comportamientos del agua subterránea con mayor claridad. Esta integración de FloPy con Jupyter no solo aumenta la productividad de los modeladores, sino que también abre nuevas posibilidades para la comunicación de hallazgos y la toma de decisiones informadas basadas en datos. En resumen, el uso de Jupyter y FloPy juntos proporciona una herramienta poderosa y flexible para el modelado hidrogeológico, combinando la capacidad de procesamiento de Python con la interactividad y visualización de Jupyter.
document.addEventListener(‘DOMContentLoaded’, () => { try { if (window.MathJax) { MathJax.Hub.Config({ TeX: { equationNumbers: { autoNumber: ‘AMS’, useLabelIds: true } }, tex2jax: { inlineMath: [[‘$’, ‘$’], [‘\\(‘, ‘\\)’]], displayMath: [[‘$$’, ‘$$’], [‘\\[‘, ‘\\]’]], processEscapes: true, processEnvironments: true }, displayAlign: ‘center’, CommonHTML: { linebreaks: { automatic: true } } }); MathJax.Hub.Queue([‘Typeset’, MathJax.Hub]); } } catch (error) { console.error(‘Error initializing MathJax:’, error); // Puedes registrar el error en la consola sin mostrarlo en la página // También puedes enviar el error a un servidor de registro de errores si es necesario } });` is generated by Mermaid in place of `\n`, * but _any_ “malformed” tag will break the SVG rendering entirely. */ const RE_VOID_ELEMENT = /]*?)\s*>/gi; /** * Ensure a void element is closed with a slash, preserving any attributes. */ function replaceVoidElement(match, tag, rest) { rest = rest.trim(); if (!rest.endsWith(‘/’)) { rest = `${rest} /`; } return “; } void Promise.all([…diagrams].map(renderOneMarmaid)); });
In [1]:
from pathlib import Path
from tempfile import TemporaryDirectory
In [2]:
import matplotlib.pyplot as plt
import numpy as np
In [ ]:
import flopy
In [4]:
h1 = 100
Nlay = 10
N = 101
L = 400.0
H = 50.0
k = 1.0
q = -1000.0
In [5]:
temp_dir = TemporaryDirectory()
workspace = temp_dir.name
name = "tutorial01_mf6"
In [6]:
temp_dir
Out[6]:
<TemporaryDirectory '/var/folders/k8/cxrxpj813s3bxcy0ry98xtxh0000gn/T/tmpg40d3a1_'>
In [7]:
sim = flopy.mf6.MFSimulation(
sim_name=name, exe_name="/Users/roberto/MODFLOW/mac/mf6", version="mf6", sim_ws=workspace
)
In [8]:
tdis = flopy.mf6.ModflowTdis(
sim, pname="tdis", time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)]
)
In [9]:
ims = flopy.mf6.ModflowIms(
sim,
pname="ims",
complexity="SIMPLE",
linear_acceleration="BICGSTAB",
)
In [10]:
model_nam_file = f"{name}.nam"
gwf = flopy.mf6.ModflowGwf(
sim,
modelname=name,
model_nam_file=model_nam_file,
save_flows=True,
newtonoptions="NEWTON UNDER_RELAXATION",
)
In [11]:
bot = np.linspace(-H / Nlay, -H, Nlay)
delrow = delcol = L / (N - 1)
dis = flopy.mf6.ModflowGwfdis(
gwf,
nlay=Nlay,
nrow=N,
ncol=N,
delr=delrow,
delc=delcol,
top=0.0,
botm=bot,
)
In [12]:
dis
Out[12]:
package_name = dis
filename = tutorial01_mf6.dis
package_type = dis
model_or_simulation_package = model
model_name = tutorial01_mf6
Block dimensions
--------------------
nlay
{internal}
(10)
nrow
{internal}
(101)
ncol
{internal}
(101)
Block griddata
--------------------
delr
{constant 4.0}
delc
{constant 4.0}
top
{constant 0.0}
botm
Layer_1{constant -5.0}
Layer_2{constant -10.0}
Layer_3{constant -15.0}
Layer_4{constant -20.0}
Layer_5{constant -25.0}
Layer_6{constant -30.0}
Layer_7{constant -35.0}
Layer_8{constant -40.0}
Layer_9{constant -45.0}
Layer_10{constant -50.0}
In [13]:
start = h1 * np.ones((Nlay, N, N))
ic = flopy.mf6.ModflowGwfic(gwf, pname="ic", strt=start)
In [14]:
npf = flopy.mf6.ModflowGwfnpf(
gwf,
icelltype=1,
k=k,
)
In [15]:
chd_rec = []
layer = 0
for row_col in range(0, N):
chd_rec.append(((layer, row_col, 0), h1))
chd_rec.append(((layer, row_col, N - 1), h1))
if row_col != 0 and row_col != N - 1:
chd_rec.append(((layer, 0, row_col), h1))
chd_rec.append(((layer, N - 1, row_col), h1))
chd = flopy.mf6.ModflowGwfchd(
gwf,
stress_period_data=chd_rec,
)
In [16]:
iper = 0
ra = chd.stress_period_data.get_data(key=iper)
ra
Out[16]:
rec.array([((0, 0, 0), 100.), ((0, 0, 100), 100.), ((0, 1, 0), 100.),
((0, 1, 100), 100.), ((0, 0, 1), 100.), ((0, 100, 1), 100.),
((0, 2, 0), 100.), ((0, 2, 100), 100.), ((0, 0, 2), 100.),
((0, 100, 2), 100.), ((0, 3, 0), 100.), ((0, 3, 100), 100.),
((0, 0, 3), 100.), ((0, 100, 3), 100.), ((0, 4, 0), 100.),
((0, 4, 100), 100.), ((0, 0, 4), 100.), ((0, 100, 4), 100.),
((0, 5, 0), 100.), ((0, 5, 100), 100.), ((0, 0, 5), 100.),
((0, 100, 5), 100.), ((0, 6, 0), 100.), ((0, 6, 100), 100.),
((0, 0, 6), 100.), ((0, 100, 6), 100.), ((0, 7, 0), 100.),
((0, 7, 100), 100.), ((0, 0, 7), 100.), ((0, 100, 7), 100.),
((0, 8, 0), 100.), ((0, 8, 100), 100.), ((0, 0, 8), 100.),
((0, 100, 8), 100.), ((0, 9, 0), 100.), ((0, 9, 100), 100.),
((0, 0, 9), 100.), ((0, 100, 9), 100.), ((0, 10, 0), 100.),
((0, 10, 100), 100.), ((0, 0, 10), 100.), ((0, 100, 10), 100.),
((0, 11, 0), 100.), ((0, 11, 100), 100.), ((0, 0, 11), 100.),
((0, 100, 11), 100.), ((0, 12, 0), 100.), ((0, 12, 100), 100.),
((0, 0, 12), 100.), ((0, 100, 12), 100.), ((0, 13, 0), 100.),
((0, 13, 100), 100.), ((0, 0, 13), 100.), ((0, 100, 13), 100.),
((0, 14, 0), 100.), ((0, 14, 100), 100.), ((0, 0, 14), 100.),
((0, 100, 14), 100.), ((0, 15, 0), 100.), ((0, 15, 100), 100.),
((0, 0, 15), 100.), ((0, 100, 15), 100.), ((0, 16, 0), 100.),
((0, 16, 100), 100.), ((0, 0, 16), 100.), ((0, 100, 16), 100.),
((0, 17, 0), 100.), ((0, 17, 100), 100.), ((0, 0, 17), 100.),
((0, 100, 17), 100.), ((0, 18, 0), 100.), ((0, 18, 100), 100.),
((0, 0, 18), 100.), ((0, 100, 18), 100.), ((0, 19, 0), 100.),
((0, 19, 100), 100.), ((0, 0, 19), 100.), ((0, 100, 19), 100.),
((0, 20, 0), 100.), ((0, 20, 100), 100.), ((0, 0, 20), 100.),
((0, 100, 20), 100.), ((0, 21, 0), 100.), ((0, 21, 100), 100.),
((0, 0, 21), 100.), ((0, 100, 21), 100.), ((0, 22, 0), 100.),
((0, 22, 100), 100.), ((0, 0, 22), 100.), ((0, 100, 22), 100.),
((0, 23, 0), 100.), ((0, 23, 100), 100.), ((0, 0, 23), 100.),
((0, 100, 23), 100.), ((0, 24, 0), 100.), ((0, 24, 100), 100.),
((0, 0, 24), 100.), ((0, 100, 24), 100.), ((0, 25, 0), 100.),
((0, 25, 100), 100.), ((0, 0, 25), 100.), ((0, 100, 25), 100.),
((0, 26, 0), 100.), ((0, 26, 100), 100.), ((0, 0, 26), 100.),
((0, 100, 26), 100.), ((0, 27, 0), 100.), ((0, 27, 100), 100.),
((0, 0, 27), 100.), ((0, 100, 27), 100.), ((0, 28, 0), 100.),
((0, 28, 100), 100.), ((0, 0, 28), 100.), ((0, 100, 28), 100.),
((0, 29, 0), 100.), ((0, 29, 100), 100.), ((0, 0, 29), 100.),
((0, 100, 29), 100.), ((0, 30, 0), 100.), ((0, 30, 100), 100.),
((0, 0, 30), 100.), ((0, 100, 30), 100.), ((0, 31, 0), 100.),
((0, 31, 100), 100.), ((0, 0, 31), 100.), ((0, 100, 31), 100.),
((0, 32, 0), 100.), ((0, 32, 100), 100.), ((0, 0, 32), 100.),
((0, 100, 32), 100.), ((0, 33, 0), 100.), ((0, 33, 100), 100.),
((0, 0, 33), 100.), ((0, 100, 33), 100.), ((0, 34, 0), 100.),
((0, 34, 100), 100.), ((0, 0, 34), 100.), ((0, 100, 34), 100.),
((0, 35, 0), 100.), ((0, 35, 100), 100.), ((0, 0, 35), 100.),
((0, 100, 35), 100.), ((0, 36, 0), 100.), ((0, 36, 100), 100.),
((0, 0, 36), 100.), ((0, 100, 36), 100.), ((0, 37, 0), 100.),
((0, 37, 100), 100.), ((0, 0, 37), 100.), ((0, 100, 37), 100.),
((0, 38, 0), 100.), ((0, 38, 100), 100.), ((0, 0, 38), 100.),
((0, 100, 38), 100.), ((0, 39, 0), 100.), ((0, 39, 100), 100.),
((0, 0, 39), 100.), ((0, 100, 39), 100.), ((0, 40, 0), 100.),
((0, 40, 100), 100.), ((0, 0, 40), 100.), ((0, 100, 40), 100.),
((0, 41, 0), 100.), ((0, 41, 100), 100.), ((0, 0, 41), 100.),
((0, 100, 41), 100.), ((0, 42, 0), 100.), ((0, 42, 100), 100.),
((0, 0, 42), 100.), ((0, 100, 42), 100.), ((0, 43, 0), 100.),
((0, 43, 100), 100.), ((0, 0, 43), 100.), ((0, 100, 43), 100.),
((0, 44, 0), 100.), ((0, 44, 100), 100.), ((0, 0, 44), 100.),
((0, 100, 44), 100.), ((0, 45, 0), 100.), ((0, 45, 100), 100.),
((0, 0, 45), 100.), ((0, 100, 45), 100.), ((0, 46, 0), 100.),
((0, 46, 100), 100.), ((0, 0, 46), 100.), ((0, 100, 46), 100.),
((0, 47, 0), 100.), ((0, 47, 100), 100.), ((0, 0, 47), 100.),
((0, 100, 47), 100.), ((0, 48, 0), 100.), ((0, 48, 100), 100.),
((0, 0, 48), 100.), ((0, 100, 48), 100.), ((0, 49, 0), 100.),
((0, 49, 100), 100.), ((0, 0, 49), 100.), ((0, 100, 49), 100.),
((0, 50, 0), 100.), ((0, 50, 100), 100.), ((0, 0, 50), 100.),
((0, 100, 50), 100.), ((0, 51, 0), 100.), ((0, 51, 100), 100.),
((0, 0, 51), 100.), ((0, 100, 51), 100.), ((0, 52, 0), 100.),
((0, 52, 100), 100.), ((0, 0, 52), 100.), ((0, 100, 52), 100.),
((0, 53, 0), 100.), ((0, 53, 100), 100.), ((0, 0, 53), 100.),
((0, 100, 53), 100.), ((0, 54, 0), 100.), ((0, 54, 100), 100.),
((0, 0, 54), 100.), ((0, 100, 54), 100.), ((0, 55, 0), 100.),
((0, 55, 100), 100.), ((0, 0, 55), 100.), ((0, 100, 55), 100.),
((0, 56, 0), 100.), ((0, 56, 100), 100.), ((0, 0, 56), 100.),
((0, 100, 56), 100.), ((0, 57, 0), 100.), ((0, 57, 100), 100.),
((0, 0, 57), 100.), ((0, 100, 57), 100.), ((0, 58, 0), 100.),
((0, 58, 100), 100.), ((0, 0, 58), 100.), ((0, 100, 58), 100.),
((0, 59, 0), 100.), ((0, 59, 100), 100.), ((0, 0, 59), 100.),
((0, 100, 59), 100.), ((0, 60, 0), 100.), ((0, 60, 100), 100.),
((0, 0, 60), 100.), ((0, 100, 60), 100.), ((0, 61, 0), 100.),
((0, 61, 100), 100.), ((0, 0, 61), 100.), ((0, 100, 61), 100.),
((0, 62, 0), 100.), ((0, 62, 100), 100.), ((0, 0, 62), 100.),
((0, 100, 62), 100.), ((0, 63, 0), 100.), ((0, 63, 100), 100.),
((0, 0, 63), 100.), ((0, 100, 63), 100.), ((0, 64, 0), 100.),
((0, 64, 100), 100.), ((0, 0, 64), 100.), ((0, 100, 64), 100.),
((0, 65, 0), 100.), ((0, 65, 100), 100.), ((0, 0, 65), 100.),
((0, 100, 65), 100.), ((0, 66, 0), 100.), ((0, 66, 100), 100.),
((0, 0, 66), 100.), ((0, 100, 66), 100.), ((0, 67, 0), 100.),
((0, 67, 100), 100.), ((0, 0, 67), 100.), ((0, 100, 67), 100.),
((0, 68, 0), 100.), ((0, 68, 100), 100.), ((0, 0, 68), 100.),
((0, 100, 68), 100.), ((0, 69, 0), 100.), ((0, 69, 100), 100.),
((0, 0, 69), 100.), ((0, 100, 69), 100.), ((0, 70, 0), 100.),
((0, 70, 100), 100.), ((0, 0, 70), 100.), ((0, 100, 70), 100.),
((0, 71, 0), 100.), ((0, 71, 100), 100.), ((0, 0, 71), 100.),
((0, 100, 71), 100.), ((0, 72, 0), 100.), ((0, 72, 100), 100.),
((0, 0, 72), 100.), ((0, 100, 72), 100.), ((0, 73, 0), 100.),
((0, 73, 100), 100.), ((0, 0, 73), 100.), ((0, 100, 73), 100.),
((0, 74, 0), 100.), ((0, 74, 100), 100.), ((0, 0, 74), 100.),
((0, 100, 74), 100.), ((0, 75, 0), 100.), ((0, 75, 100), 100.),
((0, 0, 75), 100.), ((0, 100, 75), 100.), ((0, 76, 0), 100.),
((0, 76, 100), 100.), ((0, 0, 76), 100.), ((0, 100, 76), 100.),
((0, 77, 0), 100.), ((0, 77, 100), 100.), ((0, 0, 77), 100.),
((0, 100, 77), 100.), ((0, 78, 0), 100.), ((0, 78, 100), 100.),
((0, 0, 78), 100.), ((0, 100, 78), 100.), ((0, 79, 0), 100.),
((0, 79, 100), 100.), ((0, 0, 79), 100.), ((0, 100, 79), 100.),
((0, 80, 0), 100.), ((0, 80, 100), 100.), ((0, 0, 80), 100.),
((0, 100, 80), 100.), ((0, 81, 0), 100.), ((0, 81, 100), 100.),
((0, 0, 81), 100.), ((0, 100, 81), 100.), ((0, 82, 0), 100.),
((0, 82, 100), 100.), ((0, 0, 82), 100.), ((0, 100, 82), 100.),
((0, 83, 0), 100.), ((0, 83, 100), 100.), ((0, 0, 83), 100.),
((0, 100, 83), 100.), ((0, 84, 0), 100.), ((0, 84, 100), 100.),
((0, 0, 84), 100.), ((0, 100, 84), 100.), ((0, 85, 0), 100.),
((0, 85, 100), 100.), ((0, 0, 85), 100.), ((0, 100, 85), 100.),
((0, 86, 0), 100.), ((0, 86, 100), 100.), ((0, 0, 86), 100.),
((0, 100, 86), 100.), ((0, 87, 0), 100.), ((0, 87, 100), 100.),
((0, 0, 87), 100.), ((0, 100, 87), 100.), ((0, 88, 0), 100.),
((0, 88, 100), 100.), ((0, 0, 88), 100.), ((0, 100, 88), 100.),
((0, 89, 0), 100.), ((0, 89, 100), 100.), ((0, 0, 89), 100.),
((0, 100, 89), 100.), ((0, 90, 0), 100.), ((0, 90, 100), 100.),
((0, 0, 90), 100.), ((0, 100, 90), 100.), ((0, 91, 0), 100.),
((0, 91, 100), 100.), ((0, 0, 91), 100.), ((0, 100, 91), 100.),
((0, 92, 0), 100.), ((0, 92, 100), 100.), ((0, 0, 92), 100.),
((0, 100, 92), 100.), ((0, 93, 0), 100.), ((0, 93, 100), 100.),
((0, 0, 93), 100.), ((0, 100, 93), 100.), ((0, 94, 0), 100.),
((0, 94, 100), 100.), ((0, 0, 94), 100.), ((0, 100, 94), 100.),
((0, 95, 0), 100.), ((0, 95, 100), 100.), ((0, 0, 95), 100.),
((0, 100, 95), 100.), ((0, 96, 0), 100.), ((0, 96, 100), 100.),
((0, 0, 96), 100.), ((0, 100, 96), 100.), ((0, 97, 0), 100.),
((0, 97, 100), 100.), ((0, 0, 97), 100.), ((0, 100, 97), 100.),
((0, 98, 0), 100.), ((0, 98, 100), 100.), ((0, 0, 98), 100.),
((0, 100, 98), 100.), ((0, 99, 0), 100.), ((0, 99, 100), 100.),
((0, 0, 99), 100.), ((0, 100, 99), 100.), ((0, 100, 0), 100.),
((0, 100, 100), 100.)],
dtype=[('cellid', 'O'), ('head', '<f8')])
In [17]:
wel_rec = [(Nlay - 1, int(N / 4), int(N / 4), q)]
wel = flopy.mf6.ModflowGwfwel(
gwf,
stress_period_data=wel_rec,
)
In [18]:
headfile = f"{name}.hds"
head_filerecord = [headfile]
budgetfile = f"{name}.cbb"
budget_filerecord = [budgetfile]
saverecord = [("HEAD", "ALL"), ("BUDGET", "ALL")]
printrecord = [("HEAD", "LAST")]
oc = flopy.mf6.ModflowGwfoc(
gwf,
saverecord=saverecord,
head_filerecord=head_filerecord,
budget_filerecord=budget_filerecord,
printrecord=printrecord,
)
In [19]:
sim.write_simulation()
writing simulation...
writing simulation name file...
writing simulation tdis package...
writing ims package ims...
writing model tutorial01_mf6...
writing model name file...
writing package dis...
writing package ic...
writing package npf...
writing package chd_0...
INFORMATION: maxbound in ('gwf6', 'chd', 'dimensions') changed to 400 based on size of stress_period_data
writing package wel_0...
INFORMATION: maxbound in ('gwf6', 'wel', 'dimensions') changed to 1 based on size of stress_period_data
writing package oc...
In [20]:
success, buff = sim.run_simulation()
assert success, "MODFLOW 6 did not terminate normally."
FloPy is using the following executable to run the model: /Users/roberto/MODFLOW/mac/mf6
MODFLOW 6
U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
VERSION 6.4.1 Release 12/09/2022
MODFLOW 6 compiled Dec 10 2022 04:47:08 with Intel(R) Fortran Intel(R) 64
Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
Build 20220726_000000
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): 2023/03/29 9:35:03
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): 2023/03/29 9:35:05
Elapsed run time: 1.559 Seconds
Normal termination of simulation.
In [20]:
h = gwf.output.head().get_data(kstpkper=(0, 0))
x = y = np.linspace(0, L, N)
y = y[::-1]
vmin, vmax = 90.0, 100.0
contour_intervals = np.arange(90, 100.1, 1.0)
In [21]:
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(1, 1, 1, aspect="equal")
c = ax.contour(x, y, h[0], contour_intervals, colors="black")
plt.clabel(c, fmt="%2.1f")
Out[21]:
<a list of 6 text.Text objects>
In [22]:
fig, axes = plt.subplots(2, 1, figsize=(6, 12), constrained_layout=True)
# first subplot
ax = axes[0]
ax.set_title("Model Layer 1")
modelmap = flopy.plot.PlotMapView(model=gwf, ax=ax)
pa = modelmap.plot_array(h, vmin=vmin, vmax=vmax)
quadmesh = modelmap.plot_bc("CHD")
linecollection = modelmap.plot_grid(lw=0.5, color="0.5")
contours = modelmap.contour_array(
h,
levels=contour_intervals,
colors="black",
)
ax.clabel(contours, fmt="%2.1f")
cb = plt.colorbar(pa, shrink=0.5, ax=ax)
# second subplot
ax = axes[1]
ax.set_title(f"Model Layer {Nlay}")
modelmap = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=Nlay - 1)
linecollection = modelmap.plot_grid(lw=0.5, color="0.5")
pa = modelmap.plot_array(h, vmin=vmin, vmax=vmax)
quadmesh = modelmap.plot_bc("CHD")
contours = modelmap.contour_array(
h,
levels=contour_intervals,
colors="black",
)
ax.clabel(contours, fmt="%2.1f")
cb = plt.colorbar(pa, shrink=0.5, ax=ax)
In [24]:
bot
Out[24]:
array([ -5., -10., -15., -20., -25., -30., -35., -40., -45., -50.])
In [ ]:
