Data a algoritmy pro tomografii
Použili jsme algoritmus regionální tomografie34 s globálními časovými údaji z Mezinárodního seismologického centra35 odpovídajícími období 1964-2014. Pro vybraný region jsme uvažovali všechna data odpovídající drahám paprsků procházejících studovaným objemem. To zahrnuje paprsky ze zemětřesení nacházejících se ve studované oblasti zaznamenané stanicemi po celém světě a paprsky z teleseismických událostí zaznamenaných stanicemi nacházejícími se ve studované oblasti (obr. S1A a B). Před jejich použitím v tomografii byla data znovu zpracována; to zahrnovalo přemístění zdrojů a odmítnutí odlehlých hodnot36. K lokalizaci událostí jsme použili jednorozměrný rychlostní model AK13537.
Tato oblast byla již dříve součástí jiného modelu vypočteného pomocí stejného algoritmu25. Předchozí studie však pokrývala pouze severní část Indie. Dále naše studie zahrnuje údaje z let 2005-2014, které v předchozí studii nebyly k dispozici. Deset let dodatečných záznamů poskytlo značné množství dat, zejména těch, která odpovídají novým stanicím v Indii, což výrazně zlepšilo pokrytí paprsků.
Inverze byla provedena odděleně v řadě překrývajících se oblastí pokrývajících celou studovanou oblast. Použili jsme tři oblasti, každou o poloměru 8° (Doplňkové materiály, obr. S1). Hloubku studovaného objemu jsme definovali na 1 000 km; většinou jsme však uvažovali výsledky do hloubky 700 km, protože hlubší struktury mohly být ovlivněny anomáliemi nacházejícími se mimo studovanou oblast. Parametrizace rozložení rychlostí byla provedena pomocí sady uzlů rozmístěných na horizontálních úrovních v hloubkách 25, 50, 75, 100, 150, 220, 290, 360, 430, 500, 570, 640, 710, 800 a 900 km. Na každé hloubkové úrovni byly uzly rozmístěny podle hustoty paprsků; hustšímu pokrytí paprsky odpovídaly menší rozestupy uzlů. Minimální vzdálenost byla stanovena na 30 km. Abychom se vyhnuli artefaktům souvisejícím s geometrií mřížky, provedli jsme výpočty pro dvě různé mřížky se základní orientací 0° a 45° a poté jsme výsledky zprůměrovali.
Inverze byla provedena současně pro anomálie rychlostí P a S a pro korekce zdrojů. Při použití dat z událostí nacházejících se ve studované oblasti jsme uvažovali čtyři neznámé parametry odpovídající posunům zdrojů v prostoru a čase. U teleseismických dat jsme provedli inverzi pro jeden parametr na událost, abychom reprezentovali nejistotu určení času mimo studovaný objem. Matice byla invertována pomocí metody LSQR38,39. Stabilita inverze byla řízena pomocí dalších rovnic určujících amplitudu a zploštění výsledných rychlostních anomálií. Hodnoty koeficientů tlumení byly nastaveny podle několika pokusů se syntetickými modely.
Výsledky tomografické inverze jsou uvedeny v hlavním článku na jednom horizontálním řezu v hloubce 100 km (obr. 2) a v doplňkových materiálech na dvou horizontálních řezech v hloubce 300 a 500 km (obr. S2) a dvou vertikálních řezech zde (obr. S3). Zde uvádíme výsledky pouze pro anomálie rychlosti vln P, protože S-data jsou téměř desetinou P-dat a výsledný S-model se nezdá být dostatečně stabilní.
V Doplňkových materiálech na obr. S4 uvádíme výsledky šachovnicového testu, které poskytují informace o prostorovém rozlišení pro získané modely. Syntetický model se skládal ze střídavě kladných a záporných anomálií o velikosti 3 % a laterální velikosti 5° × 5° km. S rostoucí hloubkou měnily anomálie znaménka ve vzdálenostech 200, 400 a 600 km. Syntetická data byla vypočtena podél stejných drah paprsků, z nichž byl odvozen model experimentálních dat, a ta byla narušena náhodným šumem s průměrnou odchylkou 0,5 s. Periodické šachovnicové anomálie byly definovány po celé Zemi, zatímco inverze byla provedena ve vybraných kruhových oblastech. To nám umožnilo prozkoumat vliv anomálií nacházejících se mimo studovanou oblast, které byly zohledněny při výpočtu syntetických dat. Výsledky obnovy šachovnice jsou uvedeny v doplňkových materiálech (obr. S4). Obecné polohy všech anomálií byly rekonstruovány správně; pozorovali jsme však určité diagonální rozmazání související s převládajícími orientacemi drah paprsků. Zjistili jsme poměrně dobré vertikální rozlišení, které nám umožnilo jasně obnovit změny znamének s hloubkou.
Dále jsme provedli syntetický test s realistickými tvary anomálií, které jsou prezentovány v horizontálních a vertikálních řezech (obr. S5 a S6). Anomálie jsou definovány v rámci řady polygonálních bloků v některých hloubkových intervalech. Výsledky obnovy ukazují, že příčné konfigurace všech anomálií jsou obnoveny správně. Ve vertikálních řezech vidíme, že anomálie představující litosféru proměnlivé tloušťky jsou vyřešeny ve správných hloubkách. Oba tyto testy podporují spolehlivost odvozených výsledků.
Výsledky tomografické inverze jsou zobrazeny ve třech horizontálních řezech (obr. 2 hlavního článku a obr. S2) a dvou vertikálních řezech (obr. S3). Kromě výsledků týkajících se Indického poloostrova zahrnuje model i některé okolní oblasti. Nejméně dvě struktury byly konzistentně získány v několika předchozích studiích, a proto mohly být použity jako přirozené měřítko pro tento model. Jedním z nejjasnějších vzorců ve většině asijských regionálních tomografických studií je dobře prozkoumaná struktura s vysokým Vp pod Pamírem-Hindukúšem, která souvisí s rozložením seismicity ve středních hloubkách (do 200 km). Snímky této anomálie s vysokým Vp byly konzistentně získávány různými autory za použití různých datových souborů a algoritmů40,41,42. Druhou referenční strukturou je protáhlá severojižně orientovaná anomálie s vysokým Vp pod barmským obloukem vyznačená seismicitou ve středních hloubkách. V našem modelu odhalujeme tuto anomálii, jak bylo uvedeno v předchozích studiích43,44,45,46. Tyto dva příklady jasně ukazují, že náš současný model je stejně stabilní i pro oblasti, které nebyly pokryty předchozími studiemi.
Modelování kontinentální kolize
Přístup k modelování
Numerický termomechanický viskoelasto-plastický 2D C-kód I2ELVIS použitý k modelování kontinentální kolize je založen na metodě konečných diferencí, je aplikován na odstupňované eulerovské síti a používá techniku marker-in-cell47,48 . Rovnice zachování hybnosti, hmoty a energie se řeší na eulerovské síti a fyzikální vlastnosti se přenášejí pomocí Lagrangeových značek, které se pohybují podle rychlostního pole interpolovaného ze sítě. V modelu je použita nenewtonovská viskoelastická reologie založená na experimentálně kalibrovaných zákonech proudění (doplňkové materiály, tabulka S1). Úplné podrobnosti této metody, které umožňují její reprodukci, jsou uvedeny jinde47,48.
Numerický návrh modelu. Výchozí nastavení modelu (Doplňkové materiály, obr. S5) je 6 000 km široké, 300 km hluboké a řešené pravidelnou pravoúhlou sítí 601 × 151 až 1 201 × 151 uzlů (liší se v různých experimentech, tab. S2) obsahující 1,8 milionu náhodně rozmístěných Lagrangeových značek. Horní a pravá hranice modelu mají mechanické okrajové podmínky volného skluzu. Na levé hranici je předepsána konstantní rychlost konvergence 4,7 cm/rok. Rychlost dolní hranice byla definována podmínkou zachování objemu výpočetní oblasti, a proto byla v každém časovém kroku zkrácena a zhuštěna. Hraniční podmínka volného povrchu nad kůrou je implementována pomocí 20 km silné „lepivé“ vrstvy vzduchu49,50 s nízkou hustotou (1 kg/m3) a viskozitou (1018 Pa-s). Počáteční tepelná a litologická struktura modelu (obr. S5) byla definována předepsáním několika kontinentálních litosférických jednotek, které odpovídají hlavním oblastem identifikovaným v rámci Indické a Euroasijské desky a liší se počátečním tepelným gradientem kontinentální litosféry (obr. S5). Zjednodušeně řečeno, původně jednotná kontinentální kůra o tloušťce 40 km se skládá z 15 km felsické svrchní kůry, 10 km střední kůry a 15 km mafické spodní kůry (tloušťka vrstev se v různých experimentech lišila, tabulka S2). Použitá původně jednotná struktura kůry je zjednodušená a zanedbává např. laterální heterogenitu tloušťky kůry indického kontinentu51,52. Toto zjednodušení je způsobeno především velkou nejistotou našich znalostí počáteční tloušťky kůry pro různé jednotky kontinentální litosféry. Tato tloušťka bude mít pravděpodobně opačný vliv než počáteční tloušťka litosféry v důsledku reologické slabosti kůry ve srovnání s litosférickým pláštěm (tab. S1). Počáteční tloušťka kůry se během numerických experimentů silně vyvíjí, přičemž kůra se převážně zesiluje v oblastech původně tenké, a tedy teplé a slabé litosféry (tab. S2). Vpravo skloněná slabá zóna litosférického měřítka (ukončený subdukční šev oceánu Tethys) označuje místo iniciace kolize podobné kolizi Indie a Eurasie. Zjednodušené lineární geotermální gradienty byly použity v různých řezech litosféry (tloušťky se v různých experimentech lišily, tab. S2) při teplotách 273 K (na povrchu) a 1 573 K (potenciální teplota pláště). V astenosférickém plášti byl původně předepsán adiabatický teplotní gradient 0,5 K/km. Pro plášť a kůru byla použita teplotně závislá tepelná vodivost (tab. S1). Tepelné okrajové podmínky jsou 273 K (horní), 1 713 K (dolní) a nulové tepelné toky na levé a pravé hranici. Pro zajištění účinného přenosu tepla z povrchu kůry je teplota „lepivého“ vzduchu/vody udržována na konstantní hodnotě 273 K. Pro model bylo použito gravitační zrychlení 9,81 m/s2 . Je třeba poznamenat, že 2D model použitý v naší studii zanedbává laterální variabilitu 3D deformačního modelu indicko-asijské kolizní zóny. Domníváme se však, že tento zjednodušený model je pro účely naší studie zaměřené na přenos deformace v čase z původně slabších do původně silnějších litosférických oblastí dostačující.
Visko-elasto-plastický reologický model
Viskozní, elastické a křehké (plastické) vlastnosti (tabulka S1) byly implementovány prostřednictvím vyhodnocení efektivní viskozity materiálu. Pro tvárné materiály, byly uvažovány příspěvky různých zákonů toku, jako je dislokační a difuzní creep, prostřednictvím výpočtu inverzní průměrné tvárné viskozity ηductile
kde ηnewt a ηpowl jsou efektivní viskozity pro difuzní a dislokační tečení, vypočtené jako
kde P je tlak, T je teplota (v K), \({\dot{\varepsilon }}_{II}=\sqrt{1/2{({\dot{\varepsilon }}}_{ij})}^{2}}\) je druhý invariant tenzoru rychlosti deformace, σcr je difuzně-dislokační přechodové napětí36 a AD, E, V a n jsou experimentálně stanovené parametry zákona toku (tabulka S1), které označují materiálovou konstantu, aktivační energii, aktivační objem a exponent napětí.
Duktilní reologie se kombinuje s křehkou (plastickou) reologií, aby se získala efektivní viskózně-plastická reologie pomocí následující horní meze pro duktální viskozitu:
kde P je tlak, ϕ je koeficient vnitřního tření (tabulka S1) a C je pevnost horniny v tahu při P = 0 (tabulka S1). Pružnost je implementována na základě nestlačitelného viskoelasto-plastického Maxwellova modelu47,48. Smykové moduly (μ) různých materiálů jsou uvedeny v tabulce S1.
Numerické výsledky
Provedli jsme 12 numerických experimentů, přičemž jsme měnili vrstvení kůry a počáteční délky a tloušťky litosféry různých modelových úseků (tabulka S2). Numerické výsledky ukazují systematickou migraci deformací od původně slabších (tj. tenčích) k původně silnějším (tj. tlustším) litosférickým úsekům spojenou s postupným růstem vnitroplošných tlakových napětí v indické i euroasijské modelové oblasti. Tato obecná tendence není ovlivněna variacemi počáteční geometrie modelu, které byly zkoumány; ovlivňuje pouze dynamiku deformace v nejslabších (tibetských, tien-šanských) litosférických úsecích (tab. S2).
.