Applications de la géophysique avec Python : débit des eaux souterraines et modélisation sismique

Les phénomènes géophysiques – du débit des eaux souterraines à la propagation des ondes sismiques – sont régis par des équations aux dérivées partielles qui décrivent comment les grandeurs physiques changent dans l’espace et dans le temps. Le poste 242 a introduit le cadre PDE plus large pour la modélisation des systèmes terrestres, couvrant les méthodes de volume fini et les outils disponibles pour s’attaquer à ces équations à l’échelle. Cet article se concentre sur deux domaines pratiques où Python a transformé les flux de travail géophysiques : modélisation des flux d’eau souterraine et modélisation des ondes sismiques.

Dans les deux domaines, Python a déplacé le paradigme d’un logiciel propriétaire en boîte noire à une recherche transparente et reproductible. Que vous configuriez une simulation d’eau souterraine Modflow avec Flopy ou que vous implémentiez un solveur d’équations à différence finie avec DeVito, Python offre la même combinaison de clarté mathématique, de flexibilité des écosystèmes et de soutien communautaire qui a remodelé le calcul scientifique dans toutes les disciplines.

Nous allons parcourir des exemples de code Python concrets, comparer les écosystèmes d’outils (Modflow/Flopy vs Devito/Sweep) et montrer comment sélectionner le bon outil pour votre problème géophysique spécifique.

Modélisation de l’écoulement des eaux souterraines avec Modflow et Flopy

L’écoulement des eaux souterraines à travers les milieux poreux est l’une des applications les plus matures de la géophysique informatique. Les équations qui gouvernent associent la loi de Darcy à la conservation de masse, produisant une PDE de type diffusion qui, malgré sa relative simplicité par rapport à la propagation des ondes sismiques, présente des défis pratiques dans la paramétrisation, les conditions aux limites et la stabilité numérique.

L’équation régissant

L’écoulement des eaux souterraines suit la Darcy’s Law combiné avec la conservation de masse, produisant la PDE de type diffusion :

∂∂x(kxx∂H∂x)+∂∂y(Kyy∂H∂y) +∂∂Z(KZZ∂H∂Z)+W=SS∂H∂t
frac{partial {partial x}left(k_{xx}frac{partial h}{partial x}right) + frac{partial {partial y}left(k_{yy}frac{partial h}{partial y}right) + frac{partial {partial z}left(k_{zz}frac{partial z}right) + w = s_sfrac{partielle h}{partielle t}

où

HH

est la tête hydraulique,

KK

est la conductivité hydraulique,

WW

représente les sources et les puits (puits, recharge), et

ss_s

est un stockage spécifique [1].

La discrétisation des différences finies sous-jacente à MODFLOW est bien adaptée aux réseaux structurés, ce qui en fait la norme de l’industrie pour les études régionales sur les eaux souterraines. Mais travailler directement avec les fichiers d’entrée de Modflow est fastidieux – c’est là que Flopy entre.

Premiers pas avec Flopy

Flopy est un package Python qui automatise les flux de travail ModFlow 6 : il génère des fichiers d’entrée, exécute des simulations et analyse la sortie pour la visualisation. L’écosystème s’est étendu au-delà de Flopy lui-même pour inclure modflow-setup (USGS, 2022) pour l’automatisation du flux de travail, SFRmaker pour le routage de l’eau de surface et Linesink-Maker pour le routage des flux. Ensemble, ils forment un écosystème complet pour la construction et l’automatisation des modèles d’eau souterraine [2].

Voici un démarrage rapide pratique de Flopy pour un modèle simple d’eau souterraine :

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

Le principal avantage de Flopy est qu’il remplace des dizaines de fichiers d’entrée Modflow par une seule API Python. Vous définissez la géométrie de la grille, les propriétés hydrauliques et les conditions aux limites en tant qu’objets Python, et Flopy traite la génération, la validation et l’exécution des fichiers.

Visualisation et résultats

Les capacités de visualisation de Flopy facilitent l’interprétation des résultats de la simulation. Après avoir exécuté une simulation :

# 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()

Cela produit des tracés de contour de la tête et une visualisation des flux vectoriels [3] – le type de sortie que les hydrogéologues utilisent pour interpréter les régimes d’écoulement, identifier les chemins d’écoulement et valider le comportement du modèle.

La sortie de visualisation de Flopy montre les contours de la tête hydrauliques recouverts de vecteurs de flux, ce qui facilite la compréhension des schémas de mouvement des eaux souterraines dans les systèmes aquifères complexes. La visualisation de démarrage rapide de Flopy démontre clairement cette capacité.

Attribution : Flopy QuickStart.png — ModFlowPy

Flux de travail pratique avec modflow-setup

Pour les projets plus importants, modflow-setup fournit une approche de configuration basée sur YAML. Construit automatiquement des grilles à partir de données source. Cela représente le passage à l’automatisation de workflow moderne. Vous définissez la structure du modèle dans les fichiers de configuration déclarative et modflow-setup gère la géométrie fastidieuse et l’affectation des propriétés.

Le flux de travail pratique ressemble à ceci :

1. Définir les propriétés de la grille et de l’aquifère (conductivité, coefficients de stockage, épaisseurs de couche)
2. Définir les conditions aux limites (têtes constantes, puits, taux de recharge)
3. Exécutez la simulation via l’exécution automatique de ModFlow 6 de Flopy
4. Visualisez les résultats avec les utilitaires de traçage de Flopy ou exportez vers Paraview pour le rendu 3D

Pour les hydrogéologues travaillant avec des modèles réglementaires d’eau souterraine, ce flux de travail est essentiel car MODFLOW reste la norme de l’industrie pour la régulation et la gestion des eaux souterraines appliquées [4]. Flopy comble le fossé entre le prototypage de la recherche et les flux de production. Vous obtenez la reproductibilité et la flexibilité de Python tout en maintenant l’acceptation réglementaire des résultats de ModFlow.

Modélisation des ondes sismiques avec Python

La modélisation des ondes sismiques est fondamentalement différente de l’écoulement des eaux souterraines. Là où les eaux souterraines suivent des équations de type diffusion qui lissent les caractéristiques nettes au fil du temps, la propagation des ondes sismiques suit la équation des ondes élastiques – une PDE hyperbolique du second ordre qui préserve les fronts d’onde nets et nécessite un traitement numérique minutieux pour éviter la dispersion artificielle.

L’équation de l’onde

L’équation des ondes élastiques régit la propagation des ondes sismiques :

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

où

𝐮mathbf{u}

est un déplacement,

𝐂mathbf{c}

est le tenseur de rigidité, et

𝐟mathbf{f}

représente les forces corporelles.

Contrairement au débit des eaux souterraines (où la pression diffuse et les gradients nets lisses), les ondes sismiques se propagent avec une dispersion minimale sur des milliers de kilomètres, ce qui signifie que les schémas numériques doivent préserver la netteté du front d’onde pour obtenir des résultats précis. Cela rend les méthodes de domaine temporel (FDTD) à différence finie particulièrement intéressantes car elles sont à la fois efficaces et précises pour les problèmes de propagation des ondes.

Devito : DSL symbolique pour les solveurs d’équations d’onde

devito (Georgia Tech Slim Group) génère un code C optimisé pour les simulations de différences finies par équation des ondes à l’aide de mathématiques symboliques symphoniques. Le framework est conçu pour les chercheurs qui souhaitent une propagation d’ondes différentiables, essentiel pour les flux de travail d’inversion de formes d’ondes complètes (FWI) et de migration inversée (RTM).

Voici l’exemple de modélisation avant de Devito TLE (ensemble de décalage temporel) du didacticiel mince :

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()

Ce code illustre le flux de travail clé : définissez une grille, spécifiez un modèle de vitesse, configurez l’opérateur d’équation d’onde et exécutez-le. Le moteur symbolique basé sur Sympy de DeVito traduit l’équation de haut niveau en un code C optimisé qui fonctionne à des performances quasi natives – une caractéristique essentielle lors de la simulation de grands ensembles de données sismiques.

L’animation de propagation d’ondes avant/arrière du groupe mince dans le modèle de Marmousi montre le type de résultats que DeVito produit : propagation claire du front d’onde grâce à un modèle de vitesse réaliste [5]. Cette animation illustre à la fois l’évolution vers l’avant et vers l’arrière des champs d’ondes, les mêmes champs d’onde qui sous-tendent les flux de travail FWI.

Attribution : Slim Forward.gif — Slim Group, Georgia Tech

Sweep 2026 : différenciation automatique pour l’inversion sismique

Le développement récent le plus important de la modélisation sismique python est sweep, un article ARXIV publié en avril 2026 qui introduit la différenciation automatique comme une alternative principale aux dérivations manuelles de l’état adjoint pour l’inversion sismique. Cela représente un changement de paradigme.

L’inversion sismique traditionnelle utilise la méthode à l’état adjoint, où les chercheurs dérivent manuellement des équations adjointes et les implémentent comme des codes de propagation à temps inverse distincts. Le balayage automatise ce processus : en utilisant la différenciation automatique (AD) avec des schémas de pas de temps natifs, les praticiens peuvent calculer les gradients directement à partir de leur modèle avant sans dériver les équations adjointes à la main.

Sweep prend en charge :

• Équations d’ondes acoustiques, élastiques et anisotropes (VTI, TTI)
• Différenciation automatique avec les backends Jax et PyTorch
• Modélisation par lots multi-GPU pour une inversion à grande échelle
• Cadre d’équations d’ondes unifiées qui remplace les codes avant/adjoint séparés

Le contraste avec DeVito est instructif : DeVito utilise des expressions DSL symboliques où vous écrivez l’équation d’onde comme formule symbolique et DeVito en génère un code C optimisé. Sweep adopte une approche différente – vous écrivez du code Python natif pour le modèle avancé et des frameworks AD (Jax, PyTorch) calculent automatiquement les gradients grâce à la différenciation en mode inverse.

Pour les praticiens, cela signifie :

• DeVito (DSL symbolique) : code C explicite et optimisé ; Nécessite un apprentissage de la syntaxe DSL ; Gradients calculés via la génération symbolique adjointe
• Sweep (différenciation automatique) : code Python natif ; Dégradés calculés automatiquement par AD ; backends flexibles (Jax, PyTorch); Potentiellement plus intuitif pour les développeurs Python

DeepWave : gradients adjoint-état accélérés GPU

deepwave (Kit Geophysical Institute) Implémente la propagation d’ondes différentiables dans PyTorch, permettant l’inversion de formes d’onde pleine accélérée par GPU et le calcul du gradient d’état adjoint. Contrairement à l’approche symbolique de DeVito, DeepWave fonctionne directement avec les tenseurs de PyTorch et calcule les gradients grâce à la différenciation automatique, mais nécessite une implémentation adjointe explicite plutôt qu’une publicité de bout en bout.

DeepWave est particulièrement utile lorsque vous souhaitez une accélération du GPU avec l’écosystème de PyTorch, et que vous êtes à l’aise pour mettre en œuvre des méthodes de l’état adjoint explicitement.

SimWave : propagation des ondes acoustiques hautes performances

Simwave, publié en 2025, comble une lacune dans l’écosystème sismique python : propagation des ondes acoustiques hautes performances avec backends C et GPU vérifiés. Il comble le fossé entre les prototypes de recherche et les outils géophysiques de qualité de production, utiles lorsque vous avez besoin de performances qui correspondent aux codes propriétaires.

Conseils pratiques : Choisir votre outil

La sélection du bon outil dépend de votre problème spécifique : domaine, échelle et rôle de votre modèle dans le flux de travail de recherche. Voici une comparaison pratique :

Guide de sélection d’outils

Meilleures pratiques pour la modélisation géophysique

Que vous modélisiez les eaux souterraines ou les ondes sismiques, ces directives pratiques s’appliquent dans les deux domaines :

1. Start with analytical benchmarks: Before running complex models, verify your implementation against known analytical solutions. Les problèmes de référence publiés par Modflow et les carnets de notes de Devito fournissent cette base.
2. Faites attention à la condition CFL : en cas de pas de temps sismique explicite, le pas de temps doit satisfaire à

   Δt≤Δx2Vmaxdelta t leq frac{delta X}{2V_{max}}où

   VMAXV_{max}est la vitesse maximale. La violation de la LCF entraîne une instabilité numérique, une erreur courante dans la modélisation sismique vers l’avant.

3. Utilisez les conditions aux limites appropriées : pour les eaux souterraines, ajoutez des zones tampons avec des limites à tête fixe loin de votre région d’intérêt. Pour les ondes sismiques, utilisez des conditions aux limites absorbantes (couches parfaitement adaptées) pour empêcher les réflexions artificielles.
4. Valider à plusieurs échelles : effectuez des études de raffinement de la grille pour garantir la convergence de vos résultats. Ceci est essentiel pour la modélisation géophysique crédible et doit être documenté pour la reproductibilité.
5. Versions et paramètres du modèle de document : enregistrez les versions logicielles exactes, les fichiers de paramètres et les graines aléatoires utilisées. Ceci est essentiel pour la reproductibilité et l’acceptation réglementaire dans la modélisation des eaux souterraines.

Erreurs courantes à éviter

Les erreurs les plus fréquentes que je vois dans les projets de modélisation géophysique sont :

• Commencer avec le mauvais solveur pour votre échelle de problème : ModFlow est optimisé pour les études régionales sur les eaux souterraines. Pour la multiphysique couplée (par exemple, flux de fluide + transport de chaleur), FIPY ou DeVito peut être plus approprié.
• Ignorer la condition de LFC dans un pas de temps explicite : cela conduit à une instabilité numérique dans les simulations sismiques. Vérifiez toujours votre pas de temps par rapport au critère de LFC avant d’exécuter de grandes simulations.
• Des suppositions initiales médiocres pour les problèmes d’inversion : les problèmes inverses sont mal posés ; Les modèles de départ médiocres conduisent à des minima locaux ou des résultats non physiques. Utilisez la régularisation (lisse, bornes), l’inversion hiérarchique (commencez avec des modèles plus simples) et les contraintes géologiques à partir de données indépendantes.
• Test de convergence inadéquat : accepter les tolérances du solveur donne des résultats inexacts de manière trop lâche ; Trop gaspille trop de calculs. Effectuez des études de convergence de maillage — Affinez votre grille jusqu’à ce que les principales sorties se stabilisent.
• Oublier de vérifier les conditions aux limites : les limites manquantes ou mal choisies introduisent des artefacts qui contaminent les résultats. Effectuez une analyse du domaine d’influence pour vous assurer que vos limites sont suffisamment éloignées de votre région d’intérêt.

Résumé et étapes suivantes

Cet article a couvert deux domaines géophysiques où Python a transformé les flux de travail de modélisation :

Débit d’eau souterraine avec modflow/flopy est la norme de l’industrie pour la modélisation des eaux souterraines appliquées. Flopy automatise les flux de travail ModFlow 6, de la construction de modèles à la visualisation, tandis que l’écosystème étendu (ModFlow-Setup, SFRMaker) fournit une automatisation modulaire pour les projets à l’échelle de la production. Le workflow pratique est le suivant : définir les propriétés de grille et d’aquifère, définir les conditions aux limites, exécuter la simulation et visualiser les résultats.

La modélisation des ondes sismiques est passée des codes propriétaires C++/Fortran aux écosystèmes Python. Le DSL symbolique de DeVito génère un code C optimisé pour une propagation d’ondes différenciées, la différenciation automatique de Sweep élimine les dérivations manuelles et DeepWave permet des flux de travail à l’état accéléré GPU dans PyTorch. Le choix entre ces outils dépend si vous préférez le DSL symbolique (DeVito) ou le Python natif avec différenciation automatique (Sweep).

Recommandation : Pour les eaux souterraines, utilisez Flopy si vous automatisez les flux de travail ModFlow. Pour une modélisation avancée sismique avec différenciation, utilisez DeVito. Pour les workflows d’inversion où vous souhaitez un calcul de gradient automatique, considérez l’approche basée sur la publicité comme une alternative moderne à la dérivation manuelle.

Clé de compromis : le DSL symbolique (devito) vous donne un code C explicite et optimisé mais nécessite l’apprentissage des expressions DSL. La différenciation automatique (Sweep, DeepWave) fonctionne avec du code natif et calcule les gradients automatiquement, avec des caractéristiques de performances différentes en fonction de votre backend (Jax vs PyTorch).

Commencer

Si vous débutez dans la modélisation géophysique avec Python, voici le chemin recommandé :

1. Définissez votre physique : notez les PDE gouvernantes avec tous les paramètres et les conditions aux limites.
2. Sélectionnez le logiciel : faites correspondre les capacités de l’outil à la complexité du problème (voir le tableau de sélection d’outils ci-dessus).
3. Créez un cas de test simple : vérifiez la mise en œuvre par rapport aux solutions analytiques ou aux repères publiés.
4. Effectuez une analyse de sensibilité : identifiez les paramètres et les hypothèses qui affectent le plus les résultats.
5. Document à fond : enregistrez les versions du modèle, les paramètres et les graines aléatoires pour la reproductibilité.

Pour une vue plus large de la façon dont les PDE sous-tendent la modélisation géophysique sur tous les systèmes terrestres, consultez le Guide d’aperçu géophysique. Pour obtenir des conseils sur le choix entre les solveurs Python PDE plus généralement, lisez l’analyse PDE Solver et l’analyse FEM vs FVM vs FDM.

Si vous avez besoin d’aide pour la conception de flux de travail géophysiques ou la sélection des bons outils pour vos recherches, nous sommes heureux de discuter de vos besoins spécifiques de projet – Demandez une consultation.

Références

[1] USGS. (2017). Documentation pour le modèle d’écoulement des eaux souterraines MODFLOW 6. Techniques et méthodes 6-A55.

[2] Bakker, M., et al. (2016). Développement de modèles de scripts ModFlow à l’aide de Python et de Flopy. eau souterraine, 54(5), 656-663.

[3] Lisez-moi de flopy github. https://github.com/modflowpy/flopy

[4] NIST. (2024). FIPY : solveur PDE de volume fini. https://pages.nist.gov/fipy/

[5] Slim Group, Georgia Tech. Modélisation sismique et solveurs d’équations d’ondes. https://slim.gatech.edu/research/modélisation

[6] Sweep 2026. Différenciation automatique pour la modélisation des ondes sismiques. arxiv : 2604.14189.

[7] Devito-example le carnet de commandes vers l’avant. https://github.com/slimgroup/devito-examples/blob/master/seismic/tutorials/tle_forward.ipynb

[8] Introduction technique DeepWave. Kit Institut géophysique. https://www.gpi.kit.edu/anglais/518_951.php