Aplicaciones de geofísica con Python: flujo de agua subterránea y modelado sísmico

Los fenómenos geofísicos, desde el flujo de agua subterránea hasta la propagación de las ondas sísmicas, se rigen por ecuaciones diferenciales parciales que describen cómo las cantidades físicas cambian a través del espacio y el tiempo. La publicación 242 introdujo el marco PDE más amplio para el modelado de sistemas de tierra, que cubre los métodos de volumen finito y las herramientas disponibles para abordar estas ecuaciones a escala. Este artículo profundiza en dos dominios prácticos donde Python ha transformado los flujos de trabajo geofísicos: Modelado de flujo de agua subterránea y Modelado de ondas sísmicas.

En ambos dominios, Python ha cambiado el paradigma de un software patentado de caja negra a una investigación transparente y reproducible. Ya sea que esté configurando una simulación de agua subterránea Modflow con Flopy o implementando un solucionador de ecuaciones de ondas de diferencia finita con DeVito, Python proporciona la misma combinación de claridad matemática, flexibilidad del ecosistema y apoyo comunitario que ha remodelado la computación científica en todas las disciplinas.

Recorreremos ejemplos de códigos de Python concretos, compararemos ecosistemas de herramientas (modflow/flopy vs Devito/Sweep) y mostraremos cómo seleccionar la herramienta adecuada para su problema geofísico específico.

Modelado de flujo de agua subterránea con ModFlow y Flopy

El flujo de agua subterránea a través de medios porosos es una de las aplicaciones más maduras de la geofísica computacional. Las ecuaciones gobernantes combinan la ley de Darcy con la conservación de masas, lo que produce una PDE de tipo difusión que, a pesar de su relativa simplicidad en comparación con la propagación de ondas sísmicas, presenta desafíos prácticos en la parametrización, las condiciones de contorno y la estabilidad numérica.

La ecuación de gobierno

El flujo de agua subterránea sigue La ley de Darcy combinada con la conservación de la masa, produciendo el PDE de tipo difusión:

∂∂X(Kxx∂H∂x)+∂∂Y(Kyy∂H∂Y) +∂∂Z(KZZ∂H∂Z)+W=SS∂H∂T
frac{parcial}{parcial x}left(k_{xx}frac{parcial h}{parcial x}derecho) + frac{parcial}{parcial y}left(k_{yy}frac{parcial h}{parcial y}derecho) + frac{parcial}{parcial z}left(k_{zz}frac{parcial h}{parcial z}derecho) + w = s_sfrac{parcial h}{parcial t}

dónde

charl

es la cabeza hidráulica,

kar

es conductividad hidráulica,

uu

representa fuentes y sumideros (pozos, recarga), y

SS_S

es un almacenamiento específico [1].

La discretización de la discretización de diferencias finitas subyacente a Modflow es adecuada para redes estructuradas, lo que lo convierte en el estándar de la industria para estudios regionales de aguas subterráneas. Pero trabajar directamente con los archivos de entrada de ModFlow es tedioso: ahí es donde entra Flopy.

Primeros pasos con Flopy

Flopy es un paquete de Python que automatiza los flujos de trabajo de Modflow 6: genera archivos de entrada, ejecuta simulaciones y analiza la salida para la visualización. El ecosistema se ha expandido más allá de Flopy mismo para incluir modflow-setup (USGS, 2022) para la automatización del flujo de trabajo, SFRmaker para el enrutamiento del agua superficial y Linesink-Maker para el enrutamiento de flujo de flujo. Juntos, forman un ecosistema integral para la construcción y automatización de modelos de aguas subterráneas [2].

Aquí hay un práctico flopy de inicio rápido para un modelo de agua subterránea simple:

import flopy
import numpy as np

# Define model parameters
model_name = 'simple_groundwater'
top = 100.0 # model top elevation
botm = 0.0 # bottom elevation
nlay = 10 # number of layers
nrow = 50 # number of rows
ncol = 50 # number of columns

# Create model simulation
sim = flopy.mf6.MFSimulation(sim_name=model_name, exe_name='mf6')
ibs = flopy.mf6.ModIbs(sim)
sts = flopy.mf6.ModSts(sim)

# Create model domain (discharge package)
dis = flopy.mf6.MFDischarge(sim, print_flows=True, save_first=True)

# Set up numerical grid
grid = flopy.mf6.MFModelGrid((nlay, nrow, ncol), top=top)

# Define hydraulic properties
k = np.full((nlay, nrow, ncol), 1e-3) # hydraulic conductivity (m/s)
hsy = np.full((nrow, ncol), 95.0) # head at surface (head boundary)
ss = np.full((nlay, nrow, ncol), 1e-4) # specific storage

# Create simulation and run
sim.create_simulation()
# ... set up stress packages, run simulation

La ventaja clave de Flopy es que reemplaza docenas de archivos de entrada de ModFlow con una sola API de Python. Define la geometría de la cuadrícula, las propiedades hidráulicas y las condiciones de los límites como objetos de Python, y Flopy maneja la generación, validación y ejecución de archivos.

Visualización y resultados

Las capacidades de visualización de Flopy hacen que sea sencillo interpretar los resultados de la simulación. Después de ejecutar una simulación:

# Load model results
mf = flopy.modflow.Modflow(sim)
h = mf.hds.get_data()

# Plot head contours and flow vectors
model.plot_head()
model.plot_vector_flow()

Esto produce gráficos de contorno de cabeza y visualización de flujo vectorial [3] — el tipo de salida que utilizan los hidrogeólogos para interpretar regímenes de flujo, identificar trayectorias de flujo y validar el comportamiento del modelo.

La salida de visualización de Flopy muestra los contornos de la cabeza hidráulico superpuestos con vectores de flujo, lo que facilita la comprensión de los patrones de movimiento del agua subterránea en sistemas acuíferos complejos. La visualización de inicio rápido de flopy demuestra claramente esta capacidad.

Atribución: Flopy QuickStart.png — Modflowpy

Flujo de trabajo práctico con ModFlow-Setup

Para proyectos más grandes, modflow-setup Proporciona un enfoque de configuración basado en YAML que Construye automáticamente cuadrículas a partir de datos de origen. Esto representa el cambio moderno hacia la automatización del flujo de trabajo: se define la estructura del modelo en los archivos de configuración declarativos y la configuración de Modflow se encarga de la geometría tediosa y la asignación de propiedades.

El flujo de trabajo práctico se ve así:

1. Definir las propiedades de la red y del acuífero (conductividad, coeficientes de almacenamiento, espesores de capa)
2. Establecer condiciones de contorno (cabezas constantes, pozos, tasas de recarga)
3. Ejecutar simulación a través de la ejecución de Modflow 6 automatizada de Flopy
4. Visualizar resultados con las utilidades de trazado de Flopy o exportar a ParaView para renderizado en 3D

Para los hidrogeólogos que trabajan con modelos regulatorios de aguas subterráneas, este flujo de trabajo es esencial porque ModFlow sigue siendo el estándar de la industria para la regulación y gestión de aguas subterráneas aplicadas [4]. Flopy cierra la brecha entre la creación de prototipos de investigación y los flujos de trabajo de producción: obtienes la reproducibilidad y flexibilidad de Python mientras mantienes la aceptación regulatoria de los resultados de ModFlow.

Modelado de ondas sísmicas con Python

El modelado de ondas sísmicas es fundamentalmente diferente del flujo de agua subterránea. Donde el agua subterránea sigue las ecuaciones de tipo difusión que suavizan las características nítidas a lo largo del tiempo, la propagación de ondas sísmicas sigue la ecuación de ondas elásticas, una PDE hiperbólica de segundo orden que preserva los frentes de onda nítidos y requiere un cuidado numérico cuidadoso para evitar la dispersión artificial.

La ecuación de onda

La ecuación de onda elástica gobierna la propagación de ondas sísmicas:

∂2𝐮∂t2=∇⋅(𝐂:∇𝐮)+𝐟
frac{parcial^2 mathbf{u}}{parcial t^2} = nabla cdot (mathbf{c} : nabla mathbf{u}) + mathbf{f}

dónde

𝐮mathbf{u}

es el desplazamiento,

𝐂Matemáticas{c}

es el tensor de rigidez, y

𝐟mathbf{f}

representa las fuerzas del cuerpo.

A diferencia del flujo de agua subterránea (donde la presión se difunde y los gradientes agudos), las ondas sísmicas se propagan con una dispersión mínima en miles de kilómetros, lo que significa que los esquemas numéricos deben preservar la nitidez del frente de onda para producir resultados precisos. Esto hace que los métodos de dominio del tiempo de diferencia finita (FDTD) sean particularmente atractivos porque son eficientes y precisos para los problemas de propagación de ondas.

DeVito: DSL simbólico para solucionadores de ecuaciones de onda

devito (Grupo Slim de Georgia Tech) genera un código C optimizado para Simulaciones de diferencia finita de onda-ecuación utilizando matemáticas simbólicas simpy. El marco está diseñado para investigadores que desean propagación de ondas diferenciables, esencial para los flujos de trabajo de inversión de forma de onda completa (FWI) y de migración de tiempo inverso (RTM).

Aquí está el ejemplo de modelado directo de Devito TLE (conjunto de tiempo de retraso) del tutorial delgado:

from devito import Grid, TimeAxis, Function, TimeFunction
from devito import Eq, Operator, where
from scipy import sparse

# Define grid dimensions
nx, nz = 200, 200
space_order = 4
time_order = 2
time_steps = 100

# Create computational grid
grid = Grid(shape=(nx, nz))

# Define velocity model (m/s)
v = Function(name='velocity', grid=grid, data=1500.0 * np.ones((nx, nz)))

# Initialize wavefield
u = TimeFunction(name='u', grid=grid, time_order=time_order)

# Set up initial conditions
src = TimeFunction(name='source', grid=grid)

# Define wave equation operator
eq = Eq(u.forward, u.back + src + v**2 * u.laplacian)

# Create and run operator
op = Operator(eq)
op()

# Analyze results
import matplotlib.pyplot as plt
plt.imshow(u.data[-1].T, cmap='viridis')
plt.title('Wavefield snapshot')
plt.show()

Este código demuestra el flujo de trabajo clave: Defina una cuadrícula, especifique un modelo de velocidad, configure el operador de ecuación de onda y ejecútelo. El motor simbólico basado en Sympy de DeVito traduce la ecuación de alto nivel en un código C optimizado que se ejecuta a un rendimiento casi nativo, una característica crítica cuando se simulan grandes conjuntos de datos sísmicos.

La animación de propagación de ondas hacia adelante y hacia atrás del grupo delgado en el modelo Marmousi muestra el tipo de resultados que produce DeVito: clara propagación del frente de onda a través de un modelo de velocidad realista [5]. Esta animación demuestra la evolución de campo de onda hacia adelante y hacia atrás, los mismos campos de onda que sustentan los flujos de trabajo de FWI.

Atribución: Slim Forward.gif — Slim Group, Georgia Tech

Barrido 2026: Diferenciación automática para inversión sísmica

El desarrollo reciente más significativo en el modelado sísmico de Python es Sweep, un artículo de ARXIV publicado en abril de 2026 que introduce la diferenciación automática como alternativa principal a las derivaciones manuales adjuntas para la inversión sísmica. Esto representa un cambio de paradigma.

La inversión sísmica tradicional utiliza el método del estado adjunto, donde los investigadores derivan manualmente ecuaciones adjuntas y las implementan como códigos de propagación de tiempo inverso separados. Sweep automatiza este proceso: al usar la diferenciación automática (AD) con esquemas de pasos de tiempo nativos, los profesionales pueden calcular gradientes directamente desde su modelo directo sin derivar ecuaciones adjuntas a mano.

Soportes de barrido:

• Ecuaciones de onda acústica, elástica y anisotrópica (VTI, TTI)
• Diferenciación automática con backends JAX y PyTorch
• Modelado por lotes multi-GPU para inversión a gran escala
• Marco de ecuaciones de onda unificado que reemplaza los códigos directos/adjuntos separados

El contraste con DeVito es instructivo: DeVito usa expresiones DSL simbólicas donde escribes la ecuación de onda como fórmula simbólica y DeVito genera un código C optimizado a partir de ella. Sweep adopta un enfoque diferente: escribe código de Python nativo para el modelo de avance y los marcos de anuncios (JAX, PyTorch) calculan gradientes automáticamente a través de la diferenciación en modo inverso.

Para los practicantes, esto significa:

• Devito (DSL simbólico): código C explícito y optimizado; Requiere sintaxis de aprendizaje de DSL; Gradientes calculados a través de la generación adjunta simbólica
• Sweep (diferenciación automática): código de Python nativo; gradientes calculados automáticamente por AD; Backends flexibles (JAX, PyTorch); Potencialmente más intuitivo para los desarrolladores de Python

Deepwave: gradientes de estado adjunto acelerados por GPU

onda profunda (Kit Geophysical Institute) Implementa la propagación de ondas diferenciables en PyTorch, lo que permite la inversión de forma de onda completa acelerada por GPU y el cálculo del gradiente de estado adjunto. A diferencia del enfoque simbólico de DeVito, DeepWave trabaja directamente con los tensores de PyTorch y calcula gradientes a través de la diferenciación automática, pero requiere una implementación adjunta explícita en lugar de un anuncio de extremo a extremo.

DeepWave es particularmente útil cuando desea una aceleración de GPU con el ecosistema de PyTorch, y se siente cómodo implementando explícitamente métodos de estado adjunto.

SIMWAVE: propagación de ondas acústicas de alto rendimiento

SimWave, publicado en 2025, llena un vacío en el ecosistema sísmico de Python: propagación de ondas acústicas de alto rendimiento con backends verificados de C y GPU. Cierra la brecha entre los prototipos de investigación y las herramientas geofísicas de producción, útiles cuando se necesita un rendimiento que coincida con los códigos propietarios.

Orientación práctica: elegir su herramienta

Seleccionar la herramienta adecuada depende de su problema específico: dominio, escala y el papel de su modelo en el flujo de trabajo de investigación. Aquí hay una comparación práctica:

Guía de selección de herramientas

Mejores prácticas para el modelado geofísico

Ya sea que esté modelando aguas subterráneas u ondas sísmicas, estas pautas prácticas se aplican en ambos dominios:

1. Comienza con benchmarks analíticos: antes de ejecutar modelos complejos, verifica tu implementación con soluciones analíticas conocidas. Los problemas de referencia publicados de ModFlow y los cuadernos tutoriales de DeVito brindan esta base.
2. Preste atención a la condición CFL: en los pasos de tiempo sísmicos explícitos, el paso de tiempo debe satisfacer

   Δt≤Δx2vmaxdelta t leq frac{delta x}{2v_{max}}donde

   vmaxv_{max}es la velocidad máxima. Violar la CFL conduce a la inestabilidad numérica, un error común en el modelado de avance sísmico.

3. Utilice las condiciones de contorno apropiadas: para las aguas subterráneas, agregue zonas de amortiguamiento con límites de cabeza fija lejos de su región de interés. Para ondas sísmicas, utilice condiciones de contorno absorbentes (capas perfectamente adaptadas) para evitar reflejos artificiales.
4. Validar a múltiples escalas: Realice estudios de refinamiento de la red para garantizar que sus resultados converjan. Esto es esencial para el modelado geofísico creíble y debe documentarse para la reproducibilidad.
5. Versiones y parámetros del modelo de documento: registre las versiones exactas del software, los archivos de parámetros y las semillas aleatorias utilizadas. Esto es esencial para la reproducibilidad y la aceptación regulatoria en el modelado de aguas subterráneas.

Errores comunes a evitar

Los errores más frecuentes que veo en los proyectos de modelado geofísico son:

• Comenzando con el solucionador incorrecto para su escala de problemas: ModFlow está optimizado para estudios regionales de aguas subterráneas. Para la multifísica acoplada (por ejemplo, flujo de fluido + transporte de calor), puede ser más apropiado Fipy o DeVito.
• Ignorar la condición CFL en pasos de tiempo explícitos: esto conduce a la inestabilidad numérica en las simulaciones sísmicas. Siempre verifique su paso de tiempo con el criterio CFL antes de ejecutar grandes simulaciones.
• Pobres conjeturas iniciales para problemas de inversión: los problemas inversos no están planteados; Los modelos de arranque deficientes conducen a mínimos locales o resultados no físicos. Use regularización (agilidad, límites), inversión jerárquica (comienza con modelos más simples) y restricciones geológicas a partir de datos independientes.
• Pruebas de convergencia inadecuadas: Aceptar las tolerancias del solucionador produce resultados inexactos de manera demasiado vaga; Demasiado desperdicios de cálculo. Realice estudios de convergencia de malla: refine su cuadrícula hasta que las salidas clave se estabilicen.
• Olvidar verificar las condiciones de los límites: los límites faltantes o mal elegidos introducen artefactos que contaminan los resultados. Realice un análisis de dominio de influencia para asegurarse de que sus límites estén lo suficientemente lejos de su región de interés.

Resumen y próximos pasos

Este artículo cubrió dos dominios geofísicos donde Python ha transformado los flujos de trabajo de modelado:

Flujo de agua subterránea con Modflow/Flopy es el estándar de la industria para el modelado de aguas subterráneas aplicadas. Flopy automatiza los flujos de trabajo de ModFlow 6, desde la construcción del modelo hasta la ejecución hasta la visualización, mientras que el ecosistema expandido (Modflow-Setup, SFRMaker) proporciona automatización modular para proyectos a escala de producción. El flujo de trabajo práctico es: Defina las propiedades de la cuadrícula y del acuífero, establezca las condiciones de los límites, ejecute la simulación y visualice los resultados.

Modelado de ondas sísmicas ha pasado de los códigos C++/Fortran patentados a los ecosistemas de Python. El DSL simbólico de DeVito genera un código C optimizado para la propagación de ondas diferenciables, la diferenciación automática de Sweep elimina las derivaciones adjuntas manuales y DeepWave habilita flujos de trabajo adjuntos acelerados por GPU en PyTorch. La elección entre estas herramientas depende de si prefiere DSL simbólico (DeVito) o Python nativo con diferenciación automática (Sweep).

Recomendación: para las aguas subterráneas, use Flopy si está automatizando los flujos de trabajo de Modflow. Para el modelado sísmico de avance con diferenciabilidad, use DeVito. Para los flujos de trabajo de inversión en los que desea calcular el gradiente automático, considere el enfoque basado en anuncios de Sweep como una alternativa moderna a la derivación adjunta manual.

Equipo de claves: DSL simbólico (DeVito) le brinda código C explícito y optimizado, pero requiere el aprendizaje de expresiones DSL. La diferenciación automática (Sweep, DeepWave) funciona con código nativo y calcula los gradientes automáticamente, con diferentes características de rendimiento dependiendo de su backend (Jax vs PyTorch).

Primeros pasos

Si es nuevo en el modelado geofísico con Python, aquí está la ruta recomendada:

1. Definir su física: anote las PDES que gobiernan con todos los parámetros y condiciones de contorno.
2. Seleccione software: haga coincidir las capacidades de la herramienta con la complejidad del problema (consulte la tabla de selección de herramientas anterior).
3. Construya un caso de prueba simple: verifique la implementación con soluciones analíticas o puntos de referencia publicados.
4. Realizar análisis de sensibilidad: Identificar qué parámetros y supuestos afectan más a los resultados.
5. Documento a fondo: Registre versiones, parámetros y semillas aleatorias para la reproducibilidad.

Para obtener una vista más amplia de cómo las PDEs respaldan el modelado geofísico en todos los sistemas terrestres, consulte la Guía de descripción general de la geofísica. Para obtener orientación sobre cómo elegir entre los solucionadores de PDE de Python en general, lea el comparación del solucionador de PDE y el fem vs fvm vs fdm.

Si necesita soporte para el diseño de flujos de trabajo geofísicos o la selección de las herramientas adecuadas para su investigación, nos complacerá discutir las necesidades de su proyecto específico: buscar una consulta.

referencias

[1] USGS. (2017). Documentación para el modelo de flujo de agua subterránea Modflow 6. Técnicas y Métodos 6-A55.

[2] Bakker, M., et al. (2016). Desarrollo de modelos de scripting Modflow usando Python y flopy. Agua subterránea, 54(5), 656-663.

[3] Flopy GitHub Léame. https://github.com/modflowpy/flopy

[4] NIST. (2024). fipy: solucionador de PDE de volumen finito. https://pages.nist.gov/fipy/

[5] Slim Group, Georgia Tech. Modelado sísmico y solucionadores de ecuación de ondas. https://slim.gatech.edu/research/modeling

[6] Sweep 2026. Diferenciación automática para modelado de ondas sísmicas. Arxiv:2604.14189.

[7] Devito-Ejemplos TLE Notebook Forward. https://github.com/slimgroup/devito-examples/blob/master/sismic/tutorials/tle_forward.ipynb

[8] Introducción técnica de onda profunda. Kit Instituto Geofísico. https://www.gpi.kit.edu/english/518_951.php