Seit Mitte der achtziger Jahre wird die Methode der seismischen Tomographie auf die Struktur des Erdmantels angewandt. Die damit ermittelten 3-D Modelle der seismischen Geschwindigkeiten in der Erde erlaubten zum ersten Mal einen Einblick in die möglichen Konvektionsmuster des Erdmantels. Die Verfeinerung dieses snapshots der Dynamik des Erdinnern gehört heute zu den Hauptaufgaben der globalen Seismologie. Die Struktur des Erdmantels, die Größe und die Raumwellenlänge der beobachteten Geschwindigkeitsperturbationen sind wichtige Randbedingungen, die in die numerische Modellierung der Dynamik des Erdinnern eingehen. Die spezielle Form der Mantelkonvektion bleibt weiterhin umstritten. Albarede und van der Hilst (1999) bemerken zum Beispiel: "(…) seismic imaging (…) anisotropy, scattering (...) needs to be improved in order to confirm or refute the existence of compositional domains in the deep mantle (…)".

Um Fragen in der Geodynamik, wie das Funktionieren der Subduktion (z.B. Lay, 1994; Zhong and Gurnis, 1995), oder das Entstehen von hot spots, beantworten zu können, bedarf es einer strukturellen Auflösung, die die Kapazitäten der bis heute verwendeten - auf Strahlentheorie basierenden - Tomographie übersteigt. In diesem Zusammenhang ist zu bemerken, dass bis heute niemand in der Lage ist, das vollständige, breitbandige Wellenfeld für eine dreidimensionale globale Erde zu berechnen, ohne auf zum Teil schwerwiegende Näherungen zu verzichten. Das heißt, dass - nachdem ein 3-D Geschwindigkeitsmodell mittels Tomographie erstellt wurde - nicht zuverlässig nachgeprüft werden kann, ob das 3-D Modell wirklich die beobachteten Wellenformen erklärt.

Aus diesen Gründen ist das Problem der Berechnung von vollständigen synthetischen Seismogrammen für eine globale 3-D Erde heute eine der größten Herausforderungen in der theoretischen globalen Seismologie. Offensichtlich ist die Berechnung von vollständigen seismischen Wellenfeldern für eine 3-D Erde nur mit numerischen Methoden zu erreichen. Dies impliziert, dass diese Aufgabe mit einem enormen rechen-technischen Aufwand verbunden ist. Teleseismische Wellenfelder enthalten Frequenzen bis in den 1Hz-Bereich und stellen damit eine große Herausforderung für eine genaue 3-D Modellierung dar.

Welche Effekte werden beobachtet, die eine Berechnung des vollständigen Wellenfeldes suggerieren? Ich möchte mit einem Beispiel aus der regionalen Seismologie beginnen: Aktive Verwerfungszonen wie der San-Andreas Graben und die Nord-Anatolische Verwerfung zeichnen sich dadurch aus, dass die (oft nahezu vertikale) Verwerfungsfläche durch eine ca. 100-200 m dicke Niedrig-geschwindigkeitszone gekennzeichnet ist. Diese hat erheblichen Einfluss auf das emittierte Wellenfeld, sofern der Bruch innerhalb dieser Zone stattfindet. Dann werden sogenannte fault zone waves (oder Kanalwellen, geführte Wellen) erzeugt, die bei Beobachtungen direkt über der Verwerfung zu Wellenzügen mit starken Amplituden und dispersivem Wellencharakter führen (z.B. Li et al., 1990). Die regionale Tomographie übersieht solche Zonen erniedrigter Geschwindigkeit, da sie nicht aufgelöst werden können. Die Analyse der Dispersion der geführten Wellen erlaubt jedoch eine Abschätzung der Dicke und ggf. der Geschwindigkeit der Schicht. Außerdem können durch die Beobachtung der geführten Wellen Beben, die innerhalb oder außerhalb dieser Zone stattfinden, unterschieden werden. Diese erlaubt eine Kartierung der Verwerfungszone mit der Tiefe, die mit den Methoden der Tomographie nicht möglich wäre. Über diesen Effekt kann also die strukturelle Auflösung verbessert werden. Ähnliche Effekte in anderem Zusammenhang wurden kürzlich von Shapiro et al. (2000) berichtet, wobei Zonen erniedrigter Geschwindigkeiten in subduzierender Lithosphäre modelliert wurden.

In der globalen Seismologie sind schon lange Effekte bekannt, die nicht im Rahmen der kugelsymmetrischen Erdmodelle bzw. der Strahlentheorie erklärbar sind. So werden im kurzperiodischen Bereich Vorläufer zur teleseismischen PKP-Phase beobachtet, die in den vollständigen theoretischen Seismo-grammen, die für kugelsymmetrische Medien exakt berechnet werden können, nicht auftreten (z.B. Doornbos und Vlaar, 1971; Haddon und Cleary, 1974; Thomas et al., 1999). Diese Energie wird durch Streuer nahe der Kern-Mantel Grenze erzeugt. Auch Streuer im ganzen Mantel wurden als Erklärung in Erwägung gezogen (Hedlin et al., 1997). Um solche Effekte zuverlässig erklären zu können, ist die Berechnung vollständiger Seismo-gramme unerlässlich. Die ist weiter unterlegt durch Beobachtungen von an der Kern-Mantel-Grenze (KMG) diffraktierter Phasen von Beben im Südpazifik, die auf dem nordamerikanischen Kontinent beobachtet werden. Für kleine azimuthale Differenzen ergeben sich starke Variationen der Amplitude der diffraktierten Phasen. Dies deutet auf starke laterale Heterogenitäten oberhalb der KMG hin (Abb. 1).

Ein weiteres Beispiel, das die Notwendigkeit der 3-D Modellierung für teleseismische Wellenfelder zeigt, sind die Beobachtungen die von Schulte-Pelkum und Earle (1999) berichtet werden: Die im südkalifornischen Netz beobachteten direkten P-Phasen zeigen deutliche Polarisationsanomalien. Die Teilchenbewegung des P-welleneinsatzes ist bis zu 15 Grad von der senkrecht zur Wellenfront erwarteten Bewegung verschieden. Dieser Effekt ist zudem frequenzabhängig, wobei die maximale Abweichungen bei einer Periode von 15 Sekunden beobachtet wird. Dieser Effekt, der entweder mit Anisotropie oder/und mit starken Heterogenitäten zusammenhängt, stellt eine große Herausforderung für eine 3-D Modellierung dar.

Wie können wir Wellenausbreitung in einer allgemein heterogenen 3-D Kugel - mit Anwendung auf unseren Planeten - modellieren? Zur vollständigen Berechnung des Wellenfeldes bieten sich u.a. numerische Methoden wie finite Differenzen (FD), finite Elemente (FE), spektrale Elemente (SE), pseudo-spektrale Methoden (PS), finite Volumen (FV) an.

Numerische Methoden zur seismischen Wellenausbreitung wurden im Bereich der Explorationsseismologie v.a. seit Mitte der achtziger Jahre intensiv entwickelt. Anfangs wurde im Besonderen die FD Methode in 2-D angewandt (z.B. Virieux, 1986; Korn, 1987), die später auf Operatoren der höheren Ordnung (z.B. Dablain, 1985) und 3-D (z.B. Witte und Richards, 1987; Mora, 1989) sowie anisotrope Medien (z.B. Igel et al., 1995)

erweitert wurde. Auch die spektralen (Fourier oder Chebyshev) Methoden (z.B. Kosloff et al., 1984; Reshef et al., 1988; Tessmer und Kosloff, 1994) wurden auf das Problem der seismischen Wellenausbreitung angewandt. Die – mit Ausnahme der Chebyshev Methode – regelmäßigen Gitter wurden auf im Raum variierende Gitter erweitert (z.B. Moczo, 1980; Jastram und Behle, 1992; Jastram und Tessmer, 1994). Dabei wurde entweder durch analytische Funktionen das Gitter gestreckt oder das Problem in mehrere Untergitter mit unterschiedlichen Gitterdichten aufgeteilt. Damit bleibt die (Matrix-) Strukturierung des Gitters erhalten.

Dabei spielte die Entwicklung der Rechentechnik eine fundamentale Rolle: Die parallele Rechnerarchitektur ermöglicht eine höchst effiziente Implementierung der seismischen Wellengleichung unter Aus-nutzung der intrinsischen Parallelität der physikalischen Erhaltungssätze und der lokalen Operatoren der FD Methode. Dies erlaubt die Berechnung vollständiger synthetischer Seismogramme für realistische Modelle, wodurch diese Methodik zahlreiche Anwendungen fand. Allerdings wurden diese Entwicklungen hauptsächlich für kartesische Geometrien vorangetrieben. Leider können die Standard FD Methoden nicht ohne weiteres auf Probleme in Kugelkoordinaten angewandt werden. Der Grund dafür liegt in der raumabhängigen Gitterelementgröße, die bei einer regelmäßigen Diskretisierung der Kugelkoordinaten (r,j ,q ) auftritt und zu Stabilitätsproblemen führt (Abb. 2). Die Stabilität eines FD Algorithmus ist proportional zum Zeitinkrement und zum reziproken Gitterabstand. Das heißt, je kleiner der Abstand, desto kleiner das notwendige Zeitinkrement, um die Simulation numerisch stabil zu halten. Im Fall der Kugel würden die notwendigen Zeitinkremente unrealistisch klein (und damit die Simulationsdauer ernorm groß) werden, wenn keine Mehrgebiets-methode verwendet würde.

In den letzten Jahren wurde mit der Entwicklung von Algorithmen begonnen, die

darauf abzielen, globale vollständige Seismogramme für lateral heterogene Strukturen berechnen zu können. Dabei gibt es bis jetzt drei Ansätze:

(1) Die elastische Wellengleichung in Kugelkoordinaten wird in einer axi-symmetrischen Form gelöst (Igel und Weber, 1995, 1996; Chaljub und Tarantola, 1997; Igel und Gudmundsson, 1997; Thomas et. al, 2000). Damit wird das Problem auf zwei Dimensionen reduziert, wobei trotzdem die 3-D Charakteristik des Wellenfeldes korrekt simuliert wird. Nachteil: auch die Modelle und Quellterme sind rotationssymmetrisch. Allerdings können viele modellabhängige Streueffekte damit untersucht werden. Außerdem können jetzt schon Wellenfelder mit dominanten Perioden bis zu 10 Sekunden simuliert werden. Mit Erweiterungen (Verbesserung der Operatoren, Paralleli-sierung) soll die Simulation von kurz-periodischen Wellenfeldern (Periode 1-2 Sekunden) erreicht werden. Ein ähnlich scheinender Ansatz ist die Berechnung des Wellenfeldes in Zylinderkoordinaten (Furumura et al., 1998). Allerdings lässt sich mit diesem Ansatz die geometrische Divergenz nicht korrekt simulieren, was die Anwendung erheblich einschränkt.

(2) Das Wellenfeld wird in einer Kugelsektion simuliert, wobei - um Probleme mit Singularitäten zu vermeiden - die Sektion um den Äquator zentriert wird und bis zu einer gewissen Maximaltiefe - also nicht bis zum Zentrum - gerechnet wird (Abb. 3, Igel, 1999).

Abb. 3: Chebyshev Gitter zur Simulation von Wellenausbreitung in einer Kugelsektion. Zentrierung dieses Gebiets um dem Äquator verhindert Probleme mit der Singularität der Wellengleichungen entlang der Achse Q=0 (Igel, 1999).

Dieser Ansatz erlaubt es, Standardmethoden für regelmäßige Gitter wie die FD bzw. PS Methode einzusetzen. Die PS Methode ist bekannt dafür, dass Randbedingungen sehr genau implementiert werden können. Damit sollen im Besonderen 3-D Effekte von Oberflächenwellen sowie azimuthale Effekte für stark heterogene geodynamische Strukturen studiert werden. Diese Methode wird besonders für Quellnahe und mittlere Epizentraldistanzen für die Datenanalyse eine Rolle spielen.

(3) Für die ganze Erde kann nicht auf die numerischen Methoden für regelmäßige Gitter zurückgegriffen werden. Es müssen Methoden angewandt werden, die allgemeine Gitter zulassen. Chaljub und Vilotte (1999) setzen auf die SE Methode, die auf einer Diskretisierung mit gekrümmten kubischen Elementen arbeitet, was allerdings zu einer unnatürlichen Beschreibung der inneren Erde mit kubischen Elementen führt. Eine Alternative ist die Lösung der Wellen-gleichung auf unstrukturierten, kugelförmigen Gittern mittels expliziter Differential-operatoren (Käser, 1999; Käser et al., 2000). Dabei kann das Innere der Erde gemäß der geometrischen Form der Heterogenitäten diskretisiert werden ( Abb. 4 und 5). Dann werden für diese Gitter die Gewichte für die Berechnung der Differentiale mit Hilfe der Methode der natürlichen Nachbarn (Sambridge and Braun, 1995) oder der Methode der finiten Volumen berechnet.

Die Genauigkeit der numerischen Lösungen soll durch Vergleich mit anderen 3-D Methoden gewährleistet werden (COSY Project, http://www.geophysik.uni-muenchen.

de/COSY). Erste Vergleiche (Igel et al., 2000) haben gezeigt, dass zum Beispiel die Ankunftszeitperturbationen für globale 3-D Modelle berechnet mit verschiedenen Methoden erheblich variieren können. Diese Ungenauigkeiten finden sich in den tomographischen Modellen als Fehler wieder. Es zeigte sich auch, dass zur Berechnung von hochfrequenten globalen Seismogrammen mit numerischen Methoden technische Verbesserungen notwendig sind.

Welche hardware Resourcen sind zur Berechnung globaler Modelle notwendig? In Tabelle 1 werden für die bereits veröffentlichten Algorithmen (1) und (2) Beispiele gegeben. Die beschriebenen Ergebnisse für einen Prozessor stellen Größenordnungen dar, wie sie heute auf leistungsfähigen PCs erreicht werden können. Das heißt, dass – sofern die Algorithmen auf Höchsleistungsrechern installiert werden – wesentlich größere Modelle gerechnet werden und somit höhere Frequenzen erreicht werden können. Hätte man für den axisymmetrischen Algorithmus 50 Gbyte zur Verfügung könnten globale Modelle mit Gitterabständen im Subkilometerbereich berechnet werden. Mit der geeigneten Wahl der numerischen Operatoren könnten somit kurzperiodische Seismogramme (Bereich 1 Hz) simuliert werden.

Die Anwendung der oben genannten numerischen Methoden zur Berechnung der globalen Wellenausbreitung bietet sich besonders für Regionen mit starken lateralen Heterogenitäten an. Dazu gehören im Besonderen der obere und der untere Mantel. Igel und Weber (1996) untersuchten die Effekte von lateral variierender Dicke der D" Schicht auf die PcP und PdP Phasen. Mögliche Strukturen von Subduktionszonen, die mit Hilfe von numerischen Rechnungen des Subduktionsprozesses erstellt wurden (z.B. Zhong und Gurnis, 1995), wurden in den axisymmetrischen Algorithmus zur SH- Wellenausbreitung (Igel und Weber, 1995) eingegeben und Wellenfelder berechnet (Igel und Ita, 1997, Abb. 6).

Durch die Existenz der Subduktionszone werden z.B. sekundäre Oberflächenwellen, Vorläufer der Oberflächenmultiplen (z.B. SS, sSS, etc.) und zurückgestreute Phasen erzeugt. Der nächste Schritt ist hier die Untersuchung der Wellenfelder für globale Strukturen, die durch 3-D numerische Mantelkonvektions-rechnungen erhalten werden. Für kleine bis mittlere Epizentraldistanzen bietet sich hier der Kugelsektionsansatz an. Einfache Hochgeschwindigkeitsstrukturen wurden bereits mit der Chebyshevmethode auf Fokussierungseffekte untersucht (Igel, 1999). Wellenausbreitung für lateral variierende Topographie der 670 km Diskontinuität wurde von Chaljub und Tarantola (1997) mit der FD Methode berechnet. Igel und Gudmundsson (1997) berechneten synthetische Seismo-gramme für lateral variierende Strukturen des oberen Mantels mit bekannten statistischen Eigenschaften. Dabei wurde darauf hin-gewiesen, dass die an der Erdoberfläche beobachteten Laufzeitperturbationen nicht notwendigerweise die statistischen Eigen-schaften der Mantelheterogenitäten wieder-spiegeln, da finite Frequenzen Effekte (z.B. wavefront healing) je nach Verhältnis der Korrelationslänge zu Wellenlänge eine Rolle spielen.

Für die oben beschriebenen Algorithmen wird es besonders wenn die sich schnell entwickelnde hardware höhere Frequenz-bereiche erlaubt - zahlreiche weitere Anwendungen geben. So ist grundsätzlich die Berechnung vollständiger Seismogramme für heutige tomographische Modelle zu nennen. Dies betrifft sowohl die Raumwellen-tomographie wie die Oberflächenwellen-tomographie. Dabei wird die Frage im Vordergrund stehen, wie gut die 3-D Modelle die gemessenen Daten wirklich erklären. Die so ermittelten Differenzen können dann wiederum zur Verbesserung der tomo-graphischen Strukturen herangezogen werden. Zur Verbesserung der strukturellen Auflösung von Subduktionszonen, plumes, D" etc. sollten die Wellenfeldeffekte möglicher Strukturen (durch Modellrechnungen) untersucht werden, die dann im Inversionsprozess herangezogen werden können. Modellrechnungen sollten

auch als Hilfe für das Design von teleseismischen Experimenten herangezogen werden, um Experimentgeometrien etc. für die entsprechenden Fragestellungen zu optimieren.

Um die Dynamik des Erdinnern zu verstehen bedarf es höherer struktureller Auflösung wie sie die seismische Tomographie gegenwärtig bieten kann. Es sollten (wie in der regionalen Seismologie zum Beispiel bei fault zone waves) Wellenfeldeffekte zur Strukturanalyse herangezogen werden. Um diese zu simulieren, ist die numerische Berechnung des vollständigen Wellenfeldes in Kugelgeometrie notwendig. Die Anwendung numerischer Methoden auf das Problem der teleseismischen Wellenausbreitung ist ein rasch expandierendes Gebiet, das durch die rasante Entwicklung in der hardware Technologie an Bedeutung gewinnen wird.