Data och algoritmer för tomografi
Vi använde oss av den regionala tomografialgoritmen34 med globala restidsdata från det internationella seismologiska centret35 som motsvarar perioden 1964-2014. För den valda regionen tar vi hänsyn till alla data som motsvarar de strålbanor som passerar genom undersökningsvolymen. Detta inkluderar strålar från jordbävningar som ligger i undersökningsområdet och som registrerats av stationer över hela världen och strålar från teleseismiska händelser som registrerats av stationer som ligger i undersökningsområdet (fig. S1A och B). Innan uppgifterna användes för tomografi bearbetades de på nytt, vilket innebar att man flyttade källor och avvisade outliers36. För att lokalisera händelser använde vi den endimensionella hastighetsmodellen AK13537.
Detta område har tidigare ingått i en annan modell som beräknats med samma algoritm25. Den tidigare studien omfattade dock endast den norra delen av Indien. Vidare omfattar vår studie data från 2005-2014, som inte fanns tillgängliga i den tidigare studien. Tio års ytterligare registreringar har gett betydande mängder data, särskilt de som motsvarar nya stationer i Indien som drastiskt förbättrade stråltäckningen.
Inversionen utfördes separat i en rad överlappande områden som täcker hela undersökningsområdet. Vi använde tre områden, vart och ett med en radie på 8° (Supplementary Materials, Fig. S1). Vi definierade studievolymens djup till 1 000 km; vi tog dock främst hänsyn till resultat ner till ett djup på 700 km eftersom djupare strukturer kan ha påverkats av anomalier som ligger utanför studieområdet. Parametriseringen av hastighetsfördelningen utfördes med hjälp av en uppsättning noder fördelade på horisontella nivåer på 25, 50, 75, 100, 150, 220, 290, 360, 430, 500, 570, 640, 710, 800 och 900 km djup. På varje djupnivå fördelades noderna i enlighet med stråltätheten; tätare stråltäckning motsvarade mindre avstånd mellan noderna. Det minsta avståndet sattes till 30 km. För att undvika artefakter relaterade till rutnätsgeometri utförde vi beräkningarna för två olika rutnät med grundorienteringarna 0 och 45° och medelvärdesberäknade sedan resultaten.
Inversionen utfördes samtidigt för P- och S-hastighetsanomalier och för källkorrigeringar. När data från de händelser som ligger i undersökningsområdet användes, tog vi hänsyn till fyra okända parametrar som motsvarar källornas förskjutningar i rum och tid. För de teleseismiska uppgifterna inverterade vi för en parameter per händelse för att representera osäkerheten i tidsbestämningen utanför undersökningsområdet. Matrisen inverterades med hjälp av LSQR-metoden38,39. Stabiliteten hos inverteringen kontrollerades med hjälp av ytterligare ekvationer som bestämde amplituden och utplaningen av de resulterande hastighetsanomalierna. Värdena på dämpningskoefficienterna fastställdes enligt flera försök med syntetiska modeller.
Resultaten av den tomografiska inversionen visas i huvudartikeln i ett horisontellt snitt på 100 km djup (fig. 2) och i kompletterande material i två horisontella snitt på 300 och 500 km djup (fig. S2) och två vertikala snitt här (fig. S3). Här presenterar vi endast resultaten för P-vågshastighetsanomalier, eftersom S-data är nästan en tiondel av P-data, och den resulterande S-modellen verkar inte tillräckligt stabil.
I Supplementary Materials i Fig. S4 presenterar vi resultaten av checkerboard-testet, vilket ger information om den rumsliga upplösningen för de hämtade modellerna. Den syntetiska modellen bestod av omväxlande positiva och negativa anomalier med en magnitud på 3 % och en lateral storlek på 5° × 5° km. Med ökande djup bytte anomalierna tecken vid 200, 400 och 600 km. De syntetiska uppgifterna beräknades längs samma strålbanor för att härleda den experimentella datamodellen, och dessa stördes av slumpmässigt brus med en genomsnittlig avvikelse på 0,5 s. De periodiska schackbrädavanomalierna definierades över hela jorden, medan inversionen utfördes i utvalda cirkulära områden. Detta gjorde det möjligt för oss att utforska effekten av anomalier som ligger utanför undersökningsområdet och som togs i beaktande vid beräkningen av de syntetiska uppgifterna. Resultaten av återvinningen av rutnätet visas i Supplementary Materials (Fig. S4). De allmänna platserna för alla anomalier rekonstruerades korrekt, men vi observerade en viss diagonal utsmetning som hänger samman med de dominerande orienteringarna för strålbanorna. Vi fann en ganska god vertikal upplösning, vilket gjorde att vi tydligt kunde återskapa teckenförändringarna med djupet.
Därutöver har vi utfört ett syntetiskt test med realistiska former av anomalier, som presenteras i horisontella och vertikala sektioner (fig. S5 och S6). Anomalierna definieras inom en serie polygonala block i vissa djupintervall. Återvinningsresultaten visar att alla anomaliernas sidokonfigurationer återställs korrekt. I vertikala sektioner ser vi att de anomalier som representerar litosfären med varierande tjocklek är upplösta på korrekta djup. Båda dessa tester stöder tillförlitligheten hos de härledda resultaten.
Resultaten av den tomografiska inversionen visas i tre horisontella sektioner (fig. 2 i huvudartikeln och fig. S2) och två vertikala sektioner (fig. S3). Förutom de resultat som rör den indiska halvön omfattar modellen vissa omgivande områden. Minst två strukturer återfanns konsekvent i flera tidigare studier och kunde därför användas som ett naturligt riktmärke för den nuvarande modellen. Ett av de ljusaste mönstren i de flesta asiatiska regionala tomografistudier är den välstuderade strukturen med hög Vp under Pamir-Hindu Kush, som är förknippad med fördelningen av seismicitet på mellandjup (ner till 200 km). Bilder av denna höga Vp-anomali har konsekvent erhållits av olika författare med hjälp av olika datamängder och algoritmer40,41,42. Den andra riktmärkesstrukturen är en långsträckt nord-sydlig riktad hög Vp-anomali under den burmesiska bågen som är markerad med seismicitet på mellandjup. I vår modell avslöjar vi denna anomali som rapporterats i tidigare studier43,44,45,46. Dessa två exempel visar starkt att vår nuvarande modell är lika stabil för områden som inte täcktes av tidigare studier.
Modellering av kontinentala kollisioner
Modelleringsmetod
Den numeriska termomekaniska visco-elasto-plastiska 2D C-koden I2ELVIS som används för att modellera kontinentala kollisioner bygger på en finita differensmetod, som tillämpas på ett förskjutet eulerskt rutnät och som använder sig av en markör-i-cell teknik47,48. Ekvationerna för bevarande av rörelsemängd, massa och energi löses i det eulerska rutnätet, och fysiska egenskaper transporteras av lagrangska markörer som rör sig i enlighet med det hastighetsfält som interpoleras från rutnätet. I modellen används icke-newtonsk visko-elasto-plastisk reologi baserad på experimentellt kalibrerade flödeslagar (kompletterande material, tabell S1). Alla detaljer om denna metod, som gör det möjligt att reproducera den, finns på annat håll47,48.
Numerisk modellutformning. Den ursprungliga modelluppställningen (Supplementary Materials, Fig. S5) är 6 000 km bred, 300 km djup och upplöses med ett regelbundet rektangulärt rutnät med 601 × 151 till 1 201 × 151 noder (varieras i olika experiment, Tabell S2) som innehåller 1,8 miljoner slumpmässigt fördelade lagrangska markörer. Modellens övre och högra gränser har mekaniska randvillkor för fri glidning. En konstant konvergenshastighet på 4,7 cm/år föreskrivs vid den vänstra gränsen. Den nedåtgående nedre gränshastigheten bestämdes av beräkningsdomänens volymbevaringsvillkor och förkortades och förtjockades således vid varje tidssteg. Gränsvillkoret för den fria ytan ovanför jordskorpan genomförs med hjälp av ett 20 km tjockt ”klibbigt” luftskikt49,50 med låg densitet (1 kg/m3) och viskositet (1018 Pa-s). Modellens initiala termiska och litologiska struktur (fig. S5) definierades genom att föreskriva flera kontinentala litosfäriska enheter som motsvarar större regioner som identifierats inom de indiska och eurasiska plattorna och som skiljer sig åt i den initiala termiska gradienten hos den kontinentala litosfären (fig. S5). Enkelt uttryckt består den initialt enhetliga 40 km tjocka kontinentala skorpan av 15 km felsisk övre skorpa, 10 km intermediär mellanskorpa och 15 km mafisk nedre skorpa (skikttjockleken varierade i de olika experimenten, tabell S2). Den använda initialt enhetliga skorpestrukturen är förenklad och försummar t.ex. den indiska kontinentens laterala heterogenitet i skorpans tjocklek51,52. Denna förenkling beror främst på de stora osäkerheterna i vår kunskap om ursprunglig skorpetjocklek för olika kontinentala litosfäriska enheter. Denna tjocklek kommer sannolikt att ha en motsatt effekt jämfört med litosfärens initiala tjocklek på grund av skorpans reologiska svaghet jämfört med den litosfäriska manteln (tabell S1). Den initiala skorpans tjocklek utvecklas kraftigt under de numeriska experimenten, där skorpan övervägande tjocknar i områden med initialt tunn och därmed varm och svag litosfär (tabell S2). En högerinriktad svag zon i lithosfärisk skala (en avslutad subduktionssutur i Tethyshavet) markerar startplatsen för en Indien-Eurasien-liknande kollision. Förenklade linjära geotermiska gradienter användes i olika litosfäriska sektioner (tjockleken varierade i olika experiment, tabell S2) vid 273 K (vid ytan) och 1 573 K (mantelns potentiella temperatur). En adiabatisk termisk gradient på 0,5 K/km föreskrevs inledningsvis i den asthenosfäriska manteln. Temperaturberoende värmeledningsförmåga användes för manteln och skorpan (tabell S1). De termiska randvillkoren är 273 K (övre), 1 713 K (nedre) och nollvärmeflöden på vänster och höger gräns. För att säkerställa en effektiv värmeöverföring från skorpans yta hålls temperaturen för den ”klibbiga” luften/vattnet konstant på 273 K. En gravitationsacceleration på 9,81 m/s2 användes för modellen. Det bör noteras att den 2D-modell som används i vår studie försummar den laterala variabiliteten i 3D-deformationsmönstret i kollisionszonen mellan Indien och Asien. Vi anser dock att denna förenklade modell är tillräcklig för syftet med vår studie som fokuserar på överföringen av deformation med tiden från initialt svagare till initialt starkare litosfäriska regioner.
Visco-elasto-plastisk reologisk modell
De viskösa, elastiska och spröda (plastiska) egenskaperna (tabell S1) implementerades via utvärdering av materialets effektiva viskositet. För de duktila materialen, beaktades bidragen från olika flödeslagar, t.ex. dislokations- och diffusionskrypning, genom beräkning av den omvända genomsnittliga duktila viskositeten ηductile
där ηnewt och ηpowl är effektiva viskositeter för diffusion och dislokationskryp, som beräknas som
där P är trycket, T är temperaturen (i K), \({\dot{\varepsilon }}_{II}=\sqrt{1/2{({\dot{\varepsilon }}_{ij})}^{2}}}}\) är den andra invarianten i tensorn för töjningshastighet, σcr är diffusion-dislokationsövergångsspänningen36 och AD, E, V och n är experimentellt bestämda flödeslagsparametrar (tabell S1), som betecknar materialkonstanten, aktiveringsenergin, aktiveringsvolymen respektive spänningsexponenten.
Den duktila reologin kombineras med en spröd (plastisk) reologi för att ge en effektiv viskös-plastisk reologi genom att använda följande övre gräns för den duktila viskositeten:
där P är trycket, ϕ är den inre friktionskoefficienten (tabell S1) och C är bergets draghållfasthet vid P = 0 (tabell S1). Elasticiteten implementeras baserat på den inkompressibla visko-elasto-plastiska Maxwellmodellen47,48. Skermodulerna (μ) för olika material anges i tabell S1.
Numeriska resultat
Vi utförde 12 numeriska experiment där vi varierade skiktningen av jordskorpan och de initiala längderna och litosfäriska tjocklekarna för olika modellavsnitt (tabell S2). De numeriska resultaten visar systematisk migration av deformationer från initialt svagare (dvs. tunnare) till initialt starkare (dvs. tjockare) litosfäriska sektioner i samband med den gradvisa tillväxten av kompressionsspänningar inom plattan i både de indiska och eurasiska modellområdena. Denna allmänna tendens påverkas inte av variationerna i den ursprungliga modellgeometrin, som har undersökts; den påverkar endast deformationsdynamiken i de svagaste (Tibet, Tien-Shan) litosfäriska sektionerna (tabell S2).