Daten und Algorithmen für die Tomographie
Wir haben den regionalen Tomographie-Algorithmus34 mit globalen Laufzeitdaten des Internationalen Seismologischen Zentrums35 aus dem Zeitraum 1964-2014 verwendet. Für die ausgewählte Region berücksichtigen wir alle Daten, die den Strahlenpfaden entsprechen, die durch das Untersuchungsvolumen verlaufen. Dazu gehören Strahlen von Erdbeben im Untersuchungsgebiet, die von Stationen auf der ganzen Welt aufgezeichnet wurden, und Strahlen von teleseismischen Ereignissen, die von Stationen im Untersuchungsgebiet aufgezeichnet wurden (Abb. S1A und B). Vor der Verwendung in der Tomographie wurden die Daten neu verarbeitet; dabei wurden die Quellen verlagert und Ausreißer entfernt36. Zur Lokalisierung von Ereignissen verwendeten wir das eindimensionale Geschwindigkeitsmodell AK13537.
Dieses Gebiet war bereits Teil eines anderen Modells, das mit demselben Algorithmus berechnet wurde25. Die frühere Studie bezog sich jedoch nur auf den nördlichen Teil Indiens. Außerdem enthält unsere Studie Daten aus den Jahren 2005-2014, die in der vorherigen Studie nicht verfügbar waren. Zehn Jahre zusätzlicher Aufzeichnungen haben eine beträchtliche Menge an Daten geliefert, insbesondere solche, die neuen Stationen in Indien entsprechen, die die Strahlenabdeckung drastisch verbessert haben.
Die Inversion wurde separat in einer Reihe von sich überlappenden Gebieten durchgeführt, die die gesamte Studienregion abdecken. Wir verwendeten drei Regionen mit einem Radius von jeweils 8° (Ergänzende Materialien, Abb. S1). Die Tiefe des Untersuchungsvolumens wurde auf 1.000 km festgelegt; wir haben jedoch hauptsächlich Ergebnisse bis zu einer Tiefe von 700 km berücksichtigt, da tiefere Strukturen durch Anomalien außerhalb des Untersuchungsgebiets beeinflusst worden sein könnten. Die Parametrisierung der Geschwindigkeitsverteilung erfolgte anhand einer Reihe von Knoten, die auf horizontalen Ebenen in den Tiefen 25, 50, 75, 100, 150, 220, 290, 360, 430, 500, 570, 640, 710, 800 und 900 km verteilt waren. Auf jeder Tiefenebene wurden die Knoten entsprechend der Dichte der Strahlen verteilt; eine dichtere Strahlenabdeckung entsprach einem geringeren Knotenabstand. Der Mindestabstand wurde auf 30 km festgelegt. Um Artefakte im Zusammenhang mit der Gittergeometrie zu vermeiden, wurden die Berechnungen für zwei verschiedene Gitter mit Grundorientierungen von 0 und 45° durchgeführt und die Ergebnisse anschließend gemittelt.
Die Inversion wurde gleichzeitig für die P- und S-Geschwindigkeitsanomalien und für Quellenkorrekturen durchgeführt. Bei der Verwendung von Daten aus den Ereignissen im Untersuchungsgebiet wurden vier unbekannte Parameter berücksichtigt, die den Verschiebungen der Quellen in Raum und Zeit entsprechen. Bei den teleseismischen Daten wurde für jedes Ereignis ein Parameter invertiert, um die Unsicherheit der Zeitbestimmung außerhalb des Untersuchungsgebiets darzustellen. Die Matrix wurde mit der LSQR-Methode38,39 invertiert. Die Stabilität der Inversion wurde durch zusätzliche Gleichungen kontrolliert, die die Amplitude und Abflachung der resultierenden Geschwindigkeitsanomalien bestimmen. Die Werte der Dämpfungskoeffizienten wurden nach mehreren Versuchen mit synthetischen Modellen festgelegt.
Die Ergebnisse der tomographischen Inversion sind im Hauptartikel in einem horizontalen Schnitt in 100 km Tiefe (Abb. 2) und in den ergänzenden Materialien in zwei horizontalen Schnitten in 300 und 500 km Tiefe (Abb. S2) und zwei vertikalen Schnitten hier (Abb. S3) dargestellt. Hier präsentieren wir nur die Ergebnisse für die P-Wellen-Geschwindigkeitsanomalien, da die S-Daten fast ein Zehntel der P-Daten ausmachen und das resultierende S-Modell nicht ausreichend stabil erscheint.
In den ergänzenden Materialien in Abb. S4 präsentieren wir die Ergebnisse des Schachbrett-Tests, der Aufschluss über die räumliche Auflösung der gewonnenen Modelle gibt. Das synthetische Modell bestand aus abwechselnd positiven und negativen Anomalien mit einer Größenordnung von 3 % und einer lateralen Größe von 5° × 5° km. Mit zunehmender Tiefe änderten die Anomalien ihre Vorzeichen bei 200, 400 und 600 km. Die synthetischen Daten wurden entlang derselben Strahlenpfade berechnet, aus denen das Modell der experimentellen Daten abgeleitet wurde, und diese wurden durch zufälliges Rauschen mit einer durchschnittlichen Abweichung von 0,5 s gestört. Die periodischen Schachbrettanomalien wurden auf der gesamten Erde definiert, während die Inversion in ausgewählten kreisförmigen Regionen durchgeführt wurde. Dies ermöglichte es uns, die Auswirkungen von Anomalien außerhalb des Untersuchungsgebiets zu untersuchen, die bei der Berechnung der synthetischen Daten berücksichtigt wurden. Die Ergebnisse der Checkerboard-Wiederherstellung sind in den ergänzenden Materialien (Abb. S4) dargestellt. Die allgemeine Lage aller Anomalien wurde korrekt rekonstruiert; wir beobachteten jedoch einige diagonale Verschmierungen, die mit den vorherrschenden Orientierungen der Strahlengänge zusammenhängen. Die vertikale Auflösung war recht gut, so dass wir die Vorzeichenänderungen in der Tiefe deutlich erkennen konnten.
Zusätzlich haben wir einen synthetischen Test mit realistischen Formen von Anomalien durchgeführt, die in horizontalen und vertikalen Schnitten dargestellt sind (Abb. S5 und S6). Die Anomalien sind in einigen Tiefenintervallen durch eine Reihe von polygonalen Blöcken definiert. Die Ergebnisse der Wiederherstellung zeigen, dass die seitlichen Konfigurationen aller Anomalien korrekt wiederhergestellt werden. In vertikalen Schnitten sehen wir, dass die Anomalien, die die Lithosphäre mit variabler Mächtigkeit darstellen, in korrekten Tiefen aufgelöst werden. Diese beiden Tests bestätigen die Zuverlässigkeit der abgeleiteten Ergebnisse.
Die Ergebnisse der tomographischen Inversion sind in drei horizontalen Schnitten (Abb. 2 des Hauptbeitrags und Abb. S2) und zwei vertikalen Schnitten (Abb. S3) dargestellt. Abgesehen von den Ergebnissen, die sich auf die indische Halbinsel beziehen, umfasst das Modell auch einige umliegende Gebiete. Mindestens zwei Strukturen wurden in mehreren früheren Studien übereinstimmend gefunden und konnten daher als natürliche Referenz für das vorliegende Modell verwendet werden. Eines der auffälligsten Muster in den meisten regionalen Tomographiestudien Asiens ist die gut untersuchte Struktur mit hohem Vp-Wert unter dem Pamir-Hindukusch, die mit der Verteilung der Seismizität in mittleren Tiefen (bis zu 200 km) verbunden ist. Bilder dieser Anomalie mit hohem Vp-Wert wurden von verschiedenen Autoren mit unterschiedlichen Datensätzen und Algorithmen erstellt40,41,42. Die zweite Referenzstruktur ist eine längliche, in Nord-Süd-Richtung verlaufende Anomalie mit hohem Vp-Wert unter dem burmesischen Bogen, die durch die Seismizität in mittlerer Tiefe gekennzeichnet ist. In unserem Modell zeigen wir diese Anomalie, wie sie in früheren Studien berichtet wurde43,44,45,46. Diese beiden Beispiele zeigen deutlich, dass unser gegenwärtiges Modell auch für Gebiete stabil ist, die in früheren Studien nicht erfasst wurden.
Modellierung der Kontinentalkollision
Modellierungsansatz
Der zur Modellierung der Kontinentalkollision verwendete numerische thermomechanische visko-elasto-plastische 2D-C-Code I2ELVIS basiert auf einer Finite-Differenzen-Methode, die auf einem gestaffelten Euler-Gitter angewendet wird, und verwendet eine Marker-in-Cell-Technik47,48. Die Impuls-, Massen- und Energieerhaltungsgleichungen werden auf dem Eulerschen Gitter gelöst, und die physikalischen Eigenschaften werden durch Lagrangesche Marker transportiert, die sich entsprechend dem vom Gitter interpolierten Geschwindigkeitsfeld bewegen. Im Modell wird eine nicht-newtonsche visko-elasto-plastische Rheologie verwendet, die auf experimentell kalibrierten Fließgesetzen beruht (Ergänzende Materialien, Tabelle S1). Die vollständigen Einzelheiten dieser Methode, die ihre Reproduktion ermöglichen, werden an anderer Stelle47,48 beschrieben.
Numerischer Modellaufbau. Der ursprüngliche Modellaufbau (Ergänzende Materialien, Abb. S5) ist 6.000 km breit, 300 km tief und wird durch ein regelmäßiges rechteckiges Gitter mit 601 × 151 bis 1.201 × 151 Knoten (variiert in verschiedenen Experimenten, Tabelle S2) aufgelöst, das 1,8 Millionen zufällig verteilte Lagrangesche Marker enthält. Für die oberen und rechten Grenzen des Modells gelten mechanische Randbedingungen mit freiem Schlupf. An der linken Grenze ist eine konstante Konvergenzgeschwindigkeit von 4,7 cm/Jahr vorgegeben. Die nach unten gerichtete untere Grenzgeschwindigkeit wurde durch die Volumenerhaltungsbedingung des Rechengebiets definiert und wurde daher bei jedem Zeitschritt verkürzt und verdickt. Die Randbedingung der freien Oberfläche über der Kruste wird durch eine 20 km dicke „klebrige“ Luftschicht49,50 mit geringer Dichte (1 kg/m3) und Viskosität (1018 Pa-s) realisiert. Die anfängliche thermische und lithologische Struktur des Modells (Abb. S5) wurde durch die Vorgabe mehrerer kontinentaler Lithosphäreneinheiten definiert, die den wichtigsten Regionen innerhalb der indischen und eurasischen Platte entsprechen und sich durch den anfänglichen thermischen Gradienten der kontinentalen Lithosphäre unterscheiden (Abb. S5). Vereinfacht ausgedrückt besteht die anfänglich einheitliche 40 km dicke kontinentale Kruste aus 15 km felsischer Oberkruste, 10 km mittlerer Kruste und 15 km mafischer Unterkruste (die Dicke der Schichten variiert in den verschiedenen Experimenten, Tabelle S2). Die ursprünglich einheitliche Krustenstruktur ist vereinfacht und vernachlässigt z. B. die seitliche Heterogenität der Krustendicke des indischen Kontinents51,52. Diese Vereinfachung ist in erster Linie auf die großen Unsicherheiten in unserem Wissen über die anfängliche Krustendicke für verschiedene kontinentale lithosphärische Einheiten zurückzuführen. Aufgrund der rheologischen Schwäche der Kruste im Vergleich zum lithosphärischen Mantel (Tabelle S1) hat diese Dicke wahrscheinlich einen entgegengesetzten Effekt im Vergleich zur anfänglichen Dicke der Lithosphäre. Die anfängliche Krustendicke entwickelt sich während der numerischen Experimente stark, wobei sich die Kruste vorwiegend in Regionen mit ursprünglich dünner und damit warmer und schwacher Lithosphäre verdickt (Tabelle S2). Eine nach rechts geneigte lithosphärische Schwächezone (eine abgeschlossene Subduktionsnaht des Tethys-Ozeans) markiert den Initiationsort einer Indien-Eurasien-ähnlichen Kollision. Vereinfachte lineare geothermische Gradienten wurden in verschiedenen Lithosphärenabschnitten (die Dicken variierten in verschiedenen Experimenten, Tabelle S2) bei 273 K (an der Oberfläche) und 1.573 K (potenzielle Manteltemperatur) verwendet. Im asthenosphärischen Mantel wurde zunächst ein adiabatischer Temperaturgradient von 0,5 K/km vorgegeben. Für den Mantel und die Kruste wurde eine temperaturabhängige Wärmeleitfähigkeit verwendet (Tabelle S1). Die thermischen Randbedingungen sind 273 K (oben), 1.713 K (unten) und null Wärmeströme an der linken und rechten Grenze. Um eine effiziente Wärmeübertragung von der Oberfläche der Kruste zu gewährleisten, wird die Temperatur der „klebrigen“ Luft/Wasser konstant bei 273 K gehalten. Für das Modell wurde eine Schwerkraftbeschleunigung von 9,81 m/s2 verwendet. Es sei darauf hingewiesen, dass das in unserer Studie verwendete 2D-Modell die laterale Variabilität des 3D-Deformationsmusters der indisch-asiatischen Kollisionszone vernachlässigt. Wir sind jedoch der Meinung, dass dieses vereinfachte Modell für den Zweck unserer Studie, die sich auf den Transfer der Deformation mit der Zeit von anfänglich schwächeren zu anfänglich stärkeren lithosphärischen Regionen konzentriert, ausreichend ist.
Visco-elasto-plastisches rheologisches Modell
Die viskosen, elastischen und spröden (plastischen) Eigenschaften (Tabelle S1) wurden durch Auswertung der effektiven Viskosität des Materials implementiert. Für die duktilen Materialien, wurden die Beiträge verschiedener Fließgesetze wie Versetzungs- und Diffusionskriechen über die Berechnung der inversen durchschnittlichen duktilen Viskosität ηductile
wobei ηnewt und ηpowl effektive Viskositäten für Diffusion und Versetzungskriechen sind, die wie folgt berechnet werden:
wobei P der Druck, T die Temperatur (in K) ist, \({\dot{\varepsilon }}_{II}=\sqrt{1/2{({\dot{\varepsilon }}_{ij})}^{2}}\) die zweite Invariante des Dehnungsgeschwindigkeitstensors ist, σcr die Diffusionsversetzungsübergangsspannung36 ist und AD, E, V und n experimentell bestimmte Fließgesetzparameter (Tabelle S1) sind, die die Materialkonstante, die Aktivierungsenergie, das Aktivierungsvolumen bzw. den Spannungsexponenten bezeichnen.
Die duktile Rheologie wird mit einer spröden (plastischen) Rheologie kombiniert, um eine effektive viskos-plastische Rheologie zu erhalten, indem die folgende Obergrenze für die duktile Viskosität verwendet wird:
wobei P der Druck, ϕ der innere Reibungskoeffizient (Tabelle S1) und C die Gesteinszugfestigkeit bei P = 0 ist (Tabelle S1). Die Elastizität wird auf der Grundlage des inkompressiblen visko-elasto-plastischen Maxwell-Modells47,48 implementiert. Die Schermodule (μ) der verschiedenen Materialien sind in Tabelle S1 angegeben.
Numerische Ergebnisse
Wir haben 12 numerische Experimente durchgeführt und dabei die Krustenschichtung sowie die Anfangslängen und Lithosphärendicken der verschiedenen Modellabschnitte variiert (Tabelle S2). Die numerischen Ergebnisse zeigen eine systematische Migration der Deformationen von anfänglich schwächeren (d.h. dünneren) zu anfänglich stärkeren (d.h. dickeren) lithosphärischen Abschnitten, die mit dem allmählichen Anstieg der Kompressionsspannungen innerhalb der Platte sowohl in der indischen als auch in der eurasischen Modelldomäne einhergehen. Diese allgemeine Tendenz wird von den untersuchten Variationen der anfänglichen Modellgeometrie nicht beeinflusst; sie wirkt sich nur auf die Deformationsdynamik in den schwächsten (Tibet, Tien-Shan) Lithosphärenabschnitten aus (Tabelle S2).