Dane i algorytmy dla tomografii
We used the regional tomography algorithm34 with global travel time data from the International Seismological Centre35 corresponding to the period 1964-2014. Dla wybranego regionu rozpatrujemy wszelkie dane odpowiadające ścieżkom promieni przechodzących przez badany wolumen. Dotyczy to promieni pochodzących z trzęsień ziemi zlokalizowanych w badanym obszarze, zarejestrowanych przez stacje na całym świecie, oraz promieni pochodzących ze zdarzeń telesejsmicznych zarejestrowanych przez stacje zlokalizowane w badanym obszarze (Rys. S1A i B). Przed wykorzystaniem w tomografii dane zostały poddane reprocessingowi, który obejmował relokację źródeł i usunięcie wartości odstających36. Do lokalizacji zdarzeń użyto jednowymiarowego modelu prędkości AK13537.
Ten obszar był wcześniej częścią innego modelu obliczonego przy użyciu tego samego algorytmu25. Jednak poprzednie badanie obejmowało tylko północną część Indii. Ponadto nasze badanie obejmuje dane z lat 2005-2014, które nie były dostępne w poprzednim badaniu. Dziesięć lat dodatkowych zapisów dostarczyło znacznej ilości danych, zwłaszcza tych odpowiadających nowym stacjom w Indiach, które drastycznie poprawiły pokrycie promieni.
Inwersja została wykonana oddzielnie w serii nakładających się obszarów obejmujących cały badany region. Użyliśmy trzech regionów, każdy o promieniu 8° (Supplementary Materials, Fig. S1). Głębokość badanego obszaru określiliśmy na 1000 km, jednak w większości rozpatrywaliśmy wyniki do głębokości 700 km, ponieważ głębsze struktury mogły być pod wpływem anomalii zlokalizowanych poza badanym obszarem. Parametryzacja rozkładu prędkości została przeprowadzona przy użyciu zestawu węzłów rozmieszczonych na poziomach horyzontalnych na głębokościach 25, 50, 75, 100, 150, 220, 290, 360, 430, 500, 570, 640, 710, 800 i 900 km. Na każdym poziomie głębokości węzły rozmieszczono zgodnie z gęstością promieni; gęstszemu pokryciu promieniami odpowiadały mniejsze odstępy między węzłami. Minimalny odstęp został ustalony na 30 km. Aby uniknąć artefaktów związanych z geometrią siatki, obliczenia wykonano dla dwóch różnych siatek o orientacji podstawowej 0 i 45°, a następnie uśredniono wyniki.
Inwersję wykonano jednocześnie dla anomalii prędkości P i S oraz dla poprawek źródłowych. W przypadku wykorzystania danych ze zdarzeń zlokalizowanych na badanym obszarze uwzględniono cztery nieznane parametry odpowiadające przesunięciom źródeł w czasie i przestrzeni. W przypadku danych telesejsmicznych, inwertowaliśmy dla jednego parametru na zdarzenie, aby reprezentować niepewność związaną z określeniem czasu poza obszarem badań. Macierz została inwersyjna przy użyciu metody LSQR38,39. Stabilność inwersji była kontrolowana za pomocą dodatkowych równań określających amplitudę i spłaszczenie powstałych anomalii prędkości. Wartości współczynników tłumienia ustalono na podstawie kilku prób z modelami syntetycznymi.
Wyniki inwersji tomograficznej przedstawiono w pracy głównej w jednym przekroju poziomym na głębokości 100 km (Fig. 2) oraz w Materiałach uzupełniających w dwóch przekrojach poziomych na głębokości 300 i 500 km (Fig. S2) i dwóch przekrojach pionowych w tym miejscu (Fig. S3). Tutaj przedstawiamy wyniki tylko dla anomalii prędkości fali P, ponieważ dane S stanowią prawie jedną dziesiątą danych P, a wynikowy model S nie wydaje się wystarczająco stabilny.
W Materiałach Uzupełniających na Rys. S4 przedstawiamy wyniki testu szachownicy, dające informację o rozdzielczości przestrzennej dla odzyskanych modeli. Model syntetyczny składał się z naprzemiennie występujących anomalii dodatnich i ujemnych o magnitudzie 3% i rozmiarze poprzecznym 5° × 5° km. Wraz ze wzrostem głębokości, anomalie zmieniały znaki na głębokościach 200, 400 i 600 km. Dane syntetyczne zostały obliczone wzdłuż tych samych ścieżek promieni, które posłużyły do uzyskania modelu danych eksperymentalnych, i zostały zaburzone przez losowy szum o średnim odchyleniu 0.5 s. Okresowe anomalie szachownicowe zostały zdefiniowane na całej Ziemi, podczas gdy inwersja została przeprowadzona w wybranych regionach kołowych. Pozwoliło to na zbadanie wpływu anomalii zlokalizowanych poza badanym obszarem, które zostały uwzględnione przy obliczaniu danych syntetycznych. Wyniki odzyskiwania szachownicy przedstawiono w Materiałach Uzupełniających (Fig. S4). Ogólne położenie wszystkich anomalii zostało zrekonstruowane poprawnie, jednak zaobserwowano pewne ukośne rozmazania związane z dominującymi orientacjami ścieżek promieni. Stwierdziliśmy dość dobrą rozdzielczość pionową, pozwalającą na wyraźne odtworzenie zmian znaku wraz z głębokością.
Dodatkowo wykonaliśmy test syntetyczny z realistycznymi kształtami anomalii, które zostały przedstawione w przekrojach poziomych i pionowych (Rys. S5 i S6). Anomalie zdefiniowane są w ramach serii wielokątnych bloków w niektórych interwałach głębokościowych. Wyniki odzyskiwania danych pokazują, że konfiguracje poprzeczne wszystkich anomalii zostały odtworzone poprawnie. W przekrojach pionowych widzimy, że anomalie reprezentujące litosferę o zmiennej grubości są rozwiązane na prawidłowych głębokościach. Oba te testy potwierdzają wiarygodność uzyskanych wyników.
Wyniki inwersji tomograficznej przedstawiono na trzech przekrojach poziomych (Fig. 2 pracy głównej i Fig. S2) oraz dwóch przekrojach pionowych (Fig. S3). Poza wynikami związanymi z Półwyspem Indyjskim, model obejmuje niektóre obszary otaczające. Przynajmniej dwie struktury były konsekwentnie odzyskiwane w kilku poprzednich badaniach i dlatego mogą być użyte jako naturalny punkt odniesienia dla obecnego modelu. Jednym z najjaśniejszych wzorów w większości regionalnych badań tomograficznych Azji jest dobrze zbadana struktura o wysokim Vp pod Pamir-Hindukusz, która jest związana z rozkładem sejsmiczności na średnich głębokościach (do 200 km). Obrazy tej anomalii o wysokim Vp były konsekwentnie uzyskiwane przez różnych autorów przy użyciu różnych zestawów danych i algorytmów40,41,42. Drugą strukturą wzorcową jest wydłużona, skierowana na północ i południe anomalia o wysokiej Vp pod łukiem birmańskim, zaznaczona sejsmicznością na pośredniej głębokości. W naszym modelu ujawniliśmy tę anomalię w sposób opisany we wcześniejszych badaniach43,44,45,46. Te dwa przykłady silnie pokazują, że nasz obecny model jest równie stabilny dla obszarów, które nie były objęte poprzednimi badaniami.
Modelowanie kolizji kontynentów
Podejście do modelowania
Numeryczny termomechaniczny lepko-sprężysto-plastyczny 2D kod C I2ELVIS użyty do modelowania kolizji kontynentów oparty jest na metodzie różnic skończonych, zastosowany na rozłożonej siatce eulerowskiej i wykorzystuje technikę marker-w-komórce47,48. Równania zachowania pędu, masy i energii są rozwiązywane na siatce eulerowskiej, a właściwości fizyczne są transportowane przez znaczniki Lagrangiana, które poruszają się zgodnie z polem prędkości interpolowanym z siatki. W modelu zastosowano nienewtonowską reologię lepko-elasto-plastyczną opartą na doświadczalnie skalibrowanych prawach przepływu (Supplementary Materials, Table S1). Pełne szczegóły tej metody, które pozwalają na jej odtworzenie, są podane w innym miejscu47,48.
Numeryczna konstrukcja modelu. Początkowa konfiguracja modelu (Supplementary Materials, Fig. S5) jest szeroka na 6 000 km, głęboka na 300 km i rozwiązana za pomocą regularnej prostokątnej siatki od 601 × 151 do 1 201 × 151 węzłów (różne w różnych eksperymentach, Tabela S2) zawierającej 1,8 miliona losowo rozmieszczonych markerów Lagrangiana. Górna i prawa granica modelu mają mechaniczne warunki brzegowe ze swobodnym poślizgiem. Na lewej granicy zadana jest stała prędkość konwergencji 4.7 cm/rok. Dolna prędkość graniczna w dół została określona przez warunek zachowania objętości domeny obliczeniowej, a zatem była skracana i zagęszczana w każdym kroku czasowym. Warunek graniczny swobodnej powierzchni nad skorupą ziemską jest realizowany za pomocą „lepkiej” warstwy powietrza o grubości 20 km49,50 o niskiej gęstości (1 kg/m3) i lepkości (1018 Pa-s). Początkowa struktura termiczna i litologiczna modelu (Fig. S5) została zdefiniowana poprzez przypisanie kilku kontynentalnych jednostek litosferycznych odpowiadających głównym regionom zidentyfikowanym w obrębie płyty indyjskiej i euroazjatyckiej i różniących się początkowym gradientem termicznym litosfery kontynentalnej (Fig. S5). W uproszczeniu, początkowo jednolita skorupa kontynentalna o grubości 40 km składa się z 15 km felsicznej skorupy górnej, 10 km pośredniej skorupy środkowej i 15 km mafickiej skorupy dolnej (grubość warstw zmieniała się w różnych eksperymentach, Tabela S2). Przyjęta początkowo jednolita struktura skorupy jest uproszczona i nie uwzględnia np. bocznej heterogeniczności grubości skorupy kontynentu indyjskiego51,52. Uproszczenie to wynika głównie z dużej niepewności naszej wiedzy na temat początkowej grubości skorupy dla różnych kontynentalnych jednostek litosferycznych. Grubość ta będzie prawdopodobnie miała przeciwny efekt w porównaniu z początkową grubością litosfery ze względu na reologiczną słabość skorupy w porównaniu z płaszczem litosfery (Tabela S1). Początkowa grubość skorupy silnie ewoluuje podczas eksperymentów numerycznych, w których skorupa przeważnie pogrubia się w regionach początkowo cienkiej, a więc ciepłej i słabej litosfery (Tabela S2). Pochylona w prawo słaba strefa litosfery (zakończony szew subdukcyjny Oceanu Tetydy) wyznacza miejsce inicjacji kolizji na wzór indyjsko-eurazjatycki. Uproszczone liniowe gradienty geotermalne zastosowano w różnych przekrojach litosfery (grubości różniły się w różnych eksperymentach, Tabela S2) w temperaturze 273 K (na powierzchni) i 1,573 K (temperatura potencjalna płaszcza). Adiabatyczny gradient termiczny 0.5 K/km był początkowo zadany w płaszczu astenosfery. Zależne od temperatury przewodnictwo cieplne zostało użyte dla płaszcza i skorupy (Tabela S1). Termiczne warunki brzegowe to 273 K (górna), 1713 K (dolna) oraz zerowy strumień ciepła na lewej i prawej granicy. Aby zapewnić efektywny transfer ciepła z powierzchni skorupy, temperatura „lepkiego” powietrza/wody jest utrzymywana na stałym poziomie 273 K. W modelu zastosowano przyspieszenie grawitacyjne 9.81 m/s2. Należy zauważyć, że model 2D użyty w naszych badaniach nie uwzględnia bocznej zmienności trójwymiarowego modelu deformacji strefy kolizji Indie-Azja. Sądzimy jednak, że ten uproszczony model jest wystarczający dla celów naszych badań, które koncentrują się na transferze deformacji z czasem z początkowo słabszych do początkowo silniejszych regionów litosfery.
Visco-elasto-plastyczny model reologiczny
Właściwości lepkie, sprężyste i kruche (plastyczne) (Tabela S1) zostały zaimplementowane poprzez ocenę efektywnej lepkości materiału. Dla materiałów ciągliwych, uwzględniono udział różnych praw przepływu, takich jak pełzanie dyslokacyjne i dyfuzyjne, poprzez obliczenie odwrotnej średniej lepkości plastycznej ηductile
gdzie ηnewt i ηpowl są efektywnymi lepkościami dla dyfuzji i pełzania dyslokacyjnego, odpowiednio, obliczane jako
gdzie P to ciśnienie, T to temperatura (w K), \({}dot{varepsilon }}_{II}= \sqrt{1/2{({}dot{varepsilon }}_{ij})}^{2}}) jest drugim niezmiennikiem tensora szybkości odkształcania, σcr jest naprężeniem przejściowym dyfuzja-dyslokacja36, a AD, E, V, i n są eksperymentalnie wyznaczonymi parametrami prawa przepływu (Tabela S1), które oznaczają odpowiednio stałą materiałową, energię aktywacji, objętość aktywacji i wykładnik naprężenia.
Reologia plastyczna jest łączona z reologią kruchą (plastyczną) w celu uzyskania efektywnej reologii lepko-plastycznej poprzez zastosowanie następującej górnej granicy dla lepkości plastycznej:
gdzie P jest ciśnieniem, ϕ jest współczynnikiem tarcia wewnętrznego (Tabela S1), a C jest wytrzymałością skały na rozciąganie przy P = 0 (Tabela S1). Sprężystość wprowadzono w oparciu o nieściśliwy lepkosprężysto-plastyczny model Maxwella47,48. Moduły ścinania (μ) różnych materiałów są podane w Tabeli S1.
Wyniki numeryczne
Wykonaliśmy 12 eksperymentów numerycznych, zmieniając warstwowanie skorupy oraz początkowe długości i grubości litosfery w różnych przekrojach modelu (Tabela S2). Wyniki numeryczne pokazują systematyczną migrację deformacji z początkowo słabszych (tj. cieńszych) do początkowo silniejszych (tj. grubszych) sekcji litosfery, związaną ze stopniowym wzrostem wewnątrzpłytowych naprężeń ściskających zarówno w indyjskich, jak i euroazjatyckich domenach modelowych. Na tę ogólną tendencję nie mają wpływu zmiany początkowej geometrii modelu, które zostały zbadane; wpływa ona jedynie na dynamikę deformacji w najsłabszych (Tybet, Tien-Shan) sekcjach litosfery (Tabela S2).
.