Geophysikanwendungen mit Python: Grundwasserströmung und seismische Modellierung

Geophysikalische Phänomene – von der Grundwasserströmung bis zur Ausbreitung seismischer Wellen – werden durch partielle Differentialgleichungen bestimmt, die beschreiben, wie sich physikalische Größen über Raum und Zeit ändern. Post 242 führte das breitere PDE-Framework für die Modellierung von Erdsystemen ein, das die Methoden endlicher Volumen und die verfügbaren Werkzeuge für die Bewältigung dieser Gleichungen im Maßstab abdeckt. Dieser Artikel geht in zwei praktische Bereiche ein, in denen Python geophysikalische Workflows transformiert hat: Grundwasserflussmodellierung und seismische Wellenmodellierung .

In beiden Bereichen hat Python das Paradigma von proprietärer Black-Box-Software auf transparente, reproduzierbare Forschung verlagert. Unabhängig davon, ob Sie eine Modflow-Grundwassersimulation mit Flopy einrichten oder einen Finite-Differenz-Wellen-Gleichungslöser mit DeVito implementieren, bietet Python die gleiche Kombination aus mathematischer Klarheit, Ökosystemflexibilität und Community-Unterstützung, die das wissenschaftliche Rechnen über alle Disziplinen hinweg umgestaltet hat.

Wir gehen durch konkrete Python-Codebeispiele, vergleichen Werkzeug-Ökosysteme (Modflow / Flopy vs Devito / Sweep) und zeigen, wie Sie das richtige Werkzeug für Ihr spezifisches geophysikalisches Problem auswählen.

Grundwasser-Strömungsmodellierung mit Modflow und Flopy

Der Grundwasserfluss durch poröse Medien ist eine der ausgereiftesten Anwendungen der Computergeophysik. Die maßgeblichen Gleichungen kombinieren das Darcysche Gesetz mit der Massenerhaltung und ergeben eine PDE vom Diffusionstyp, die trotz ihrer relativen Einfachheit im Vergleich zur Ausbreitung der seismischen Wellen praktische Herausforderungen in Bezug auf Parametrisierung, Randbedingungen und numerische Stabilität mit sich bringt.

die maßgebliche Gleichung

Der Grundwasserfluss folgt dem Darcy’s Law in Kombination mit der Massenerhaltung, wobei die PDE vom Diffusionstyp erzeugt wird:

∂∂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}rechts) + frac{partial}{partial}{partial} y}left(k_{yy}frac{partial h}{partial y}rechts) + frac{partial{partial Z}left(k_{zz}frac{partial h}{partial z}rechts) + w = s_sfrac{partial h}{partial t}

wo

HH

ist hydraulischer Kopf,

kk

ist hydraulische Leitfähigkeit,

ww

repräsentiert Quellen und Senken (Brunnen, Aufladen) und

sss_s

ist ein spezifischer Speicher [1].

Die finite-Differenz-Diskretisierung, die dem Modflow zugrunde liegt, eignet sich gut für strukturierte Netze und ist damit der Industriestandard für regionale Grundwasserstudien. Die direkte Arbeit mit den Eingabedateien von Modflow ist jedoch mühsam – hier tritt Flopy ein.

Erste Schritte mit Flopy

Flopy ist ein Python-Paket, das Workflows von Modflow 6 automatisiert: Es generiert Eingabedateien, führt Simulationen aus und analysiert die Ausgabe für die Visualisierung. Das Ökosystem hat sich über Flopy selbst hinaus erweitert und umfasst modflow-setup (USGS, 2022) für die Workflow-Automatisierung, SFRmaker für das Routing von Oberflächenwasser und Linesink-Maker für das Streamflow-Routing. Zusammen bilden sie ein umfassendes Ökosystem für den Bau und die Automatisierung des Grundwassermodells [2].

Hier ist ein praktischer Flopy-Schnellstart für ein einfaches Grundwassermodell:

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

Der Hauptvorteil von Flopy ist, dass es Dutzende von Modflow-Eingabedateien durch eine einzelne Python-API ersetzt. Sie definieren Gittergeometrie, hydraulische Eigenschaften und Randbedingungen als Python-Objekte und Flopy-Handles für Dateigenerierung, Validierung und Ausführung.

Visualisierung und Ergebnisse

Die Visualisierungsfunktionen von Flopy machen es einfach, Simulationsergebnisse zu interpretieren. Nach dem Ausführen einer 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()

Dies erzeugt Kopfkonturdiagramme und Vektorflussvisualisierung [3] – die Art von Ausgabe, die Hydrogeologen verwenden, um Flussregime zu interpretieren, Flusspfade zu identifizieren und das Modellverhalten zu validieren.

Die Visualisierungsausgabe von Flopy zeigt hydraulische Kopfkonturen, die mit Strömungsvektoren überlagert sind, was es leichter macht, Grundwasserbewegungsmuster in komplexen Aquifersystemen zu verstehen. Die Flopy-Quickstart-Visualisierung demonstriert diese Fähigkeit deutlich.

Attribution: Flopy Quickstart.png — Modflowpy

Praktischer Workflow mit Modflow-Setup

Bei größeren Projekten bietet Modflow-Setup einen auf YAML basierenden Konfigurationsansatz, der Erstellt automatisch Gitter aus Quelldaten. Dies stellt die moderne Verschiebung in Richtung Workflow-Automatisierung dar. Sie definieren die Modellstruktur in deklarativen Konfigurationsdateien, und Modflow-Setup übernimmt die mühsame Geometrie und Eigenschaftszuweisung.

Der praktische Arbeitsablauf sieht so aus:

1. Eigenschaften des Gitters und der Grundwasserleiter (Leitfähigkeit, Speicherkoeffizienten, Schichtdicken)
2. Randbedingungen festlegen (konstante Köpfe, Brunnen, Wiederaufladeraten)
3. Simulation ausführen durch die automatisierte Modflow 6-Ausführung von Flopy
4. Ergebnisse visualisieren Mit Flopys Plot-Utilities oder Export nach ParaView für 3D-Rendering

Für Hydrogeologen, die mit regulatorischen Grundwassermodellen arbeiten, ist dieser Workflow unerlässlich, da Modflow der Industriestandard für angewandte Grundwasserregulierung und -management [4] bleibt. Flopy überbrückt die Lücke zwischen Prototyping- und Produktionsabläufen – Sie erhalten die Reproduzierbarkeit und Flexibilität von Python, während die regulatorische Akzeptanz der Modflow-Ergebnisse beibehalten wird.

Seismische Wellenmodellierung mit Python

Die Modellierung seismischer Wellen unterscheidet sich grundlegend von der Grundwasserströmung. Wenn das Grundwasser Diffusionsgleichungen folgt, die im Laufe der Zeit scharfe Merkmale glätten, folgt die Ausbreitung seismischer Wellen der Elastikwellengleichung – einer hyperbolischen PDE zweiter Ordnung, die scharfe Wellenfronten beibehält und eine sorgfältige numerische Behandlung erfordert, um eine künstliche Dispersion zu vermeiden.

Die Wellengleichung

Die elastische Wellengleichung regelt die Ausbreitung der seismischen Wellen:

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

wo

𝐮Mathbf{u}

ist Verschiebung,

𝐂Mathbf{c}

ist der Steifigkeitstensor und

𝐟Mathbf{f}

stellt Körperkräfte dar.

Im Gegensatz zu Grundwasserströmungen (bei denen der Druck diffundiert und scharfe Gradienten glätten) breiten sich seismische Wellen mit minimaler Streuung über Tausende von Kilometern aus – was bedeutet, dass numerische Schemata die Schärfe der Wellenfront beibehalten müssen, um genaue Ergebnisse zu erzielen. Dies macht die FDTD-Methode (Finite Difference Time-Domain) besonders attraktiv, da sie sowohl effizient als auch für Wellenausbreitungsprobleme sind.

Devito: Symbolische DSL für Wellengleichungslöser

Devito (Georgia Tech Slim Group) generiert optimierten C-Code für Wellengleichung Finite-Differenz-Simulationen mit symptischer symbolischer Mathematik. Das Framework richtet sich an Forscher, die eine differenzierbare Wellenausbreitung wünschen – unerlässlich für Workflows (Full-Waveform Inversion, FWI) und Reverse Time Migration (RTM).

Hier ist das Devito TLE (Time-Lag Ensemble) Vorwärtsmodellierungsbeispiel aus dem schlanken Tutorial:

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

Dieser Code demonstriert den wichtigsten Workflow: Definieren Sie ein Gitter, geben Sie ein Geschwindigkeitsmodell an, richten Sie den Operator der Wellengleichung ein und führen Sie ihn aus. Die sympy-basierte symbolische Engine von DeVito übersetzt die Gleichung auf hoher Ebene in optimierten C-Code, der mit nahezu nativer Leistung ausgeführt wird – ein kritisches Merkmal bei der Simulation großer seismischer Datensätze.

Die Vorwärts-/Rückwärts-Wellenanimation der schlanken Gruppe im Marmousi-Modell zeigt die Art der Ergebnisse, die DeVito liefert: klare Wellenfrontausbreitung durch ein realistisches Geschwindigkeitsmodell [5]. Diese Animation demonstriert sowohl vorwärts als auch rückwärts Wavefield Evolution – dieselben Wellenfelder, die FWI-Workflows untermauern.

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

Sweep 2026: Automatische Differenzierung für seismische Inversion

Die bedeutendste jüngste Entwicklung in der seismischen Python-Modellierung ist sweep , ein im April 2026 veröffentlichtes ARXIV-Papier, das die automatische Differenzierung als primäre Alternative zu manuellen Ableitungen des Adjoint-State für die seismische Inversion einführt. Dies stellt einen Paradigmenwechsel dar.

Die traditionelle seismische Inversion verwendet die Adjoint-State-Methode, bei der Forscher adjungierte Gleichungen manuell ableiten und sie als separate Ausbreitungscodes für die umgekehrte Zeit implementieren. Sweep automatisiert diesen Prozess: Durch die Verwendung der automatischen Differenzierung (AD) mit nativen Zeitschrittschemata können die Praktiker Gradienten direkt aus ihrem Vorwärtsmodell berechnen, ohne von Hand adjungierte Gleichungen abzuleiten.

Sweep unterstützt:

• Akustische, elastische und anisotrope Wellengleichungen (VTI, TTI)
• Automatische Differenzierung mit Jax- und PyTorch-Backends
• Multi-GPU-Batchmodellierung für großflächige Inversion
• Unified Wave Equation Framework, das separate Vorwärts-/Adjoint-Codes ersetzt

Der Kontrast zu DeVito ist aufschlussreich: DeVito verwendet symbolische DSL-Ausdrücke, bei denen Sie die Wellengleichung als symbolische Formel schreiben und DeVito daraus optimierten C-Code generiert. Sweep verfolgt einen anderen Ansatz – Sie schreiben nativen Python-Code für das Vorwärtsmodell, und AD-Frameworks (JAX, PyTorch) berechnen automatisch Gradienten durch die Differenzierung des Reverse-Modus.

Für Praktiker bedeutet das:

• Devito (symbolisches DSL): expliziter, optimierter C-Code; Erfordert das Lernen von DSL-Syntax; Gradienten, die über symbolische Adjoint-Generierung berechnet werden
• Sweep (automatische Differenzierung): nativer Python-Code; Gradienten, die automatisch von AD berechnet werden; Flexible Backends (JAX, PyTorch); Potenziell intuitiver für Python-Entwickler

DeepWave: GPU-beschleunigte Steigungen im Adjoint-State

DeepWave (KIT Geophysical Institute) Implementiert eine differenzierbare Wellenausbreitung in PyTorch und ermöglicht eine GPU-beschleunigte Vollwellenform-Inversion und eine Gradientenberechnung mit Adjoint-State-Gradienten. Im Gegensatz zu DeVitos symbolischem Ansatz arbeitet DeepWave direkt mit PyTorch-Tensoren und berechnet Gradienten durch automatische Differenzierung – erfordert jedoch eine explizite Adjoint-Implementierung anstelle einer End-to-End-Anzeige.

DeepWave ist besonders nützlich, wenn Sie die GPU-Beschleunigung mit dem PyTorch-Ökosystem wünschen und die explizite Implementierung von Adjoint-State-Methoden.

SimWave: Hochleistungs-Akustikwellenausbreitung

SimWave, veröffentlicht im Jahr 2025, schließt eine Lücke im seismischen Python-Ökosystem: Hochleistungs-Akustikwellenausbreitung mit verifizierten C- und GPU-Backends. Es überbrückt die Lücke zwischen Forschungsprototypen und geophysikalischen Werkzeugen in Produktionsqualität – nützlich, wenn Sie eine Leistung benötigen, die proprietäre Codes übereinstimmt.

Praktische Anleitung: Auswahl Ihres Werkzeugs

Die Auswahl des richtigen Tools hängt von Ihrem spezifischen Problem ab: Domäne, Skalierung und die Rolle Ihres Modells im Forschungsworkflow. Hier ein praktischer Vergleich:

Handbuch zur Werkzeugauswahl

Best Practices für geophysikalische Modellierung

Ob Sie Grundwasser oder seismische Wellen modellieren, diese praktischen Richtlinien gelten für beide Bereiche:

1. Starten Sie mit analytischen Benchmarks: Verifizieren Sie Ihre Implementierung anhand bekannter analytischer Lösungen, bevor Sie komplexe Modelle ausführen. Die von Modflow veröffentlichten Benchmark-Probleme und Devitos Tutorial-Notizbücher bieten diese Grundlage.
2. Auf die CFL-Bedingung achten: Beim expliziten seismischen Zeitschritt muss der Zeitschritt

   Δt≤Δx2vmaxdelta t leq frac{delta erfüllen x}{2v_{max}}wobei

   vmaxv_{max} die maximale Geschwindigkeit ist. Verstöße gegen die CFL führt zu numerischer Instabilität – ein häufiger Fehler in der seismischen Vorwärtsmodellierung.

3. Geeignete Randbedingungen verwenden: Fügen Sie für Grundwasser Pufferzonen mit festen Kopfgrenzen hinzu, die weit entfernt von Ihrer Region von Interesse sind. Verwenden Sie bei seismischen Wellen absorbierende Randbedingungen (perfekt abgestimmte Schichten), um künstliche Reflexionen zu verhindern.
4. Über mehrere Skalen validieren: Führen Sie Studien zur Rasterverfeinerung durch, um sicherzustellen, dass Ihre Ergebnisse konvergieren. Dies ist für eine glaubwürdige geophysikalische Modellierung unerlässlich und sollte für die Reproduzierbarkeit dokumentiert werden.
5. Dokumentenmodellversionen und -parameter: Zeichnen Sie die verwendeten Softwareversionen, Parameterdateien und zufälligen Seeds auf. Dies ist für die Reproduzierbarkeit und regulatorische Akzeptanz bei der Grundwassermodellierung von wesentlicher Bedeutung.

Häufige Fehler zu vermeiden

Die häufigsten Fehler, die ich in geophysikalischen Modellierungsprojekten sehe, sind:

• Angefangen mit dem falschen Solver für Ihre Problemskala: Modflow ist für regionale Grundwasserstudien optimiert. Für gekoppelte Multiphysik (z. B. Fluidfluss + Wärmetransport) kann FIPY oder DeVito geeigneter sein.
• Ignorieren der CFL-Bedingung in expliziter Zeitschritt : Dies führt zu numerischen Instabilitäten bei seismischen Simulationen. Überprüfen Sie Ihren Zeitschritt immer anhand des CFL-Kriteriums, bevor Sie große Simulationen ausführen.
• Schlechte anfängliche Vermutungen für Inversionsprobleme : Inverse Probleme sind schlecht gestellt; Schlechte Startmodelle führen zu lokalen Minima oder nicht-physischen Ergebnissen. Verwenden Sie Regularisierung (Glätte, Grenzen), hierarchische Inversion (beginnen mit einfacheren Modellen) und geologische Einschränkungen aus unabhängigen Daten.
• Unadäquate Konvergenztests: Akzeptieren von Solvertoleranzen zu lose ergibt ungenaue Ergebnisse; Zu eng wird die Berechnung verschwendet. Führen Sie Studienkonvergenzstudien durch – Verfeinern Sie Ihr Raster, bis sich die Schlüsselausgaben stabilisieren.
• Vergessen von Randbedingungen: Fehlende oder schlecht gewählte Grenzen führen Artefakte ein, die die Ergebnisse kontaminieren. Führen Sie eine Domain-of-Influence-Analyse durch, um sicherzustellen, dass Ihre Grenzen weit genug von Ihrer Region entfernt sind.

Zusammenfassung und nächste Schritte

Dieser Artikel behandelte zwei geophysikalische Domänen, in denen Python die Modellierungsworkflows transformiert hat:

Grundwasserströmung mit Modflow/Flopy ist der Industriestandard für die Modellierung der angewandten Grundwassermodellierung. Flopy automatisiert Modflow 6-Workflows – von der Modellkonstruktion über die Ausführung bis zur Visualisierung – während das erweiterte Ökosystem (ModFlow-Setup, SFRMaker) modulare Automatisierung für Projekte im Produktionsmaßstab bereitstellt. Der praktische Workflow ist: Gitter- und Aquifereigenschaften definieren, Randbedingungen festlegen, Simulation ausführen und Ergebnisse visualisieren.

Seismische Wellenmodellierung hat sich von proprietären C ++ / Fortran-Codes zu Python-Ökosystemen verlagert. Die symbolische DSL von DeVito generiert optimierten C-Code für die differenzierbare Wellenausbreitung, die automatische Differenzierung von Sweep eliminiert manuelle Adjoint-Ableitungen und DeepWave ermöglicht GPU-beschleunigte Workflows mit Adjoint-State in PyTorch. Die Wahl zwischen diesen Werkzeugen hängt davon ab, ob Sie symbolisches DSL (DeVito) oder Native Python mit automatischer Differenzierung (Sweep) bevorzugen.

Empfehlung: Verwenden Sie für Grundwasser Flopy, wenn Sie Modflow-Workflows automatisieren. Verwenden Sie für die seismische Vorwärtsmodellierung mit Differenzierbarkeit DeVito. Betrachten Sie bei Inversions-Workflows, bei denen Sie eine automatische Farbverlaufsberechnung wünschen, den AD-basierten Ansatz von Sweep als moderne Alternative zur manuellen adjunkten Ableitung.

Schlüsselkompromiss: Die symbolische DSL (DeVito) bietet Ihnen expliziten, optimierten C-Code, erfordert jedoch das Erlernen von DSL-Ausdrücken. Die automatische Differenzierung (Sweep, DeepWave) arbeitet mit nativem Code und berechnet Farbverläufe automatisch mit unterschiedlichen Leistungsmerkmalen, abhängig von Ihrem Backend (JAX vs Pytorch).

Loslegen

Wenn Sie mit Python noch nicht in der geophysikalischen Modellierung sind, finden Sie hier den empfohlenen Pfad:

1. Definieren Sie Ihre Physik: Schreiben Sie die maßgeblichen PDEs mit allen Parametern und Randbedingungen auf.
2. Auswählen von Software: Tool-Funktionen an Problemkomplexität anpassen (siehe Tabelle der Werkzeugauswahl oben).
3. Einen einfachen Testfall erstellen: Überprüfen Sie die Implementierung anhand von analytischen Lösungen oder veröffentlichten Benchmarks.
4. Sensitivitätsanalyse: Identifizieren Sie, welche Parameter und Annahmen die Ergebnisse am stärksten beeinflussen.
5. Dokument gründlich: Aufzeichnung von Modellversionen, Parametern und zufälligen Startwerten für die Reproduzierbarkeit.

Eine umfassendere Ansicht darüber, wie PDEs die geophysikalische Modellierung über alle Erdsysteme hinweg untermauern, finden Sie im Geophysik-Übersichtshandbuch . Weitere Informationen zur Auswahl zwischen Python-PDE-Lösern allgemeiner finden Sie in der Analyse PDE-Solver-Vergleich und der FEM vs FVM vs FDM -Analyse.

Wenn Sie Unterstützung beim Entwerfen geophysikalischer Workflows oder bei der Auswahl der richtigen Tools für Ihre Forschung benötigen, besprechen wir gerne Ihre spezifischen Projektanforderungen – Fordern Sie eine Beratung an .

Referenzen

[1] USGS. (2017). Dokumentation für das Modflow 6-Grundwasserflussmodell. Techniken und Methoden 6-A55.

[2] Bakker, M. et al. (2016). Scripting Modflow-Modellentwicklung mit Python und Flopy. Grundwasser, 54(5), 656-663.

[3] Flopy GitHub Readme. https://github.com/modflowpy/flopy

[4] NIST. (2024). FIPY: Endlicher PDE-Solver-Volume. https://pages.nist.gov/fipy/

[5] Slim Group, Georgia Tech. Seismische Modellierung und Wellengleichungslöser. https://slim.gatech.edu/research/modelling

[6] Sweep 2026. Automatische Differenzierung für die Modellierung seismischer Wellen. Arxiv: 2604.14189.

[7] Devito-Beispiele für das Forward-Notebook. https://github.com/slimgroup/devito-examples/blob/master/seismic/tutorials/tle_forward.ipynb

[8] Technische Einführung in die DeepWave. KIT Geophysikalisches Institut. https://www.gpi.kit.edu/english/518_951.php