Dati e algoritmi per la tomografia
Abbiamo usato l’algoritmo di tomografia regionale34 con i dati sul tempo di percorrenza globale dell’International Seismological Centre35 corrispondenti al periodo 1964-2014. Per la regione selezionata, consideriamo tutti i dati corrispondenti ai percorsi dei raggi che attraversano il volume di studio. Questo include i raggi provenienti da terremoti situati nell’area di studio registrati da stazioni di tutto il mondo, e quelli da eventi telesismici registrati da stazioni situate nell’area di studio (Fig. S1A e B). Prima del loro utilizzo nella tomografia, i dati sono stati rielaborati; questo includeva la ricollocazione delle fonti e lo scarto degli outlier36. Per localizzare gli eventi, abbiamo usato il modello di velocità unidimensionale AK13537.
Questa zona è stata precedentemente parte di un altro modello calcolato utilizzando lo stesso algoritmo25. Tuttavia, lo studio precedente copriva solo la parte settentrionale dell’India. Inoltre, il nostro studio include i dati del 2005-2014, che non erano disponibili nello studio precedente. Dieci anni di registrazioni aggiuntive hanno fornito una quantità significativa di dati, soprattutto quelli corrispondenti a nuove stazioni in India che hanno migliorato drasticamente la copertura dei raggi.
L’inversione è stata eseguita separatamente in una serie di aree sovrapposte che coprono l’intera regione di studio. Abbiamo usato tre regioni, ciascuna con un raggio di 8° (Materiali supplementari, Fig. S1). Abbiamo definito la profondità del volume di studio a 1.000 km; tuttavia, abbiamo considerato principalmente i risultati fino a una profondità di 700 km perché le strutture più profonde potrebbero essere state influenzate da anomalie situate al di fuori dell’area di studio. La parametrizzazione della distribuzione della velocità è stata eseguita utilizzando un insieme di nodi distribuiti su livelli orizzontali a profondità di 25, 50, 75, 100, 150, 220, 290, 360, 430, 500, 570, 640, 710, 800 e 900 km. Ad ogni livello di profondità, i nodi sono stati distribuiti in base alla densità dei raggi; una copertura più densa dei raggi corrispondeva ad una minore spaziatura dei nodi. La spaziatura minima è stata fissata a 30 km. Per evitare artefatti legati alla geometria della griglia, abbiamo eseguito i calcoli per due diverse griglie con orientamenti di base di 0 e 45° e poi abbiamo fatto la media dei risultati.
L’inversione è stata eseguita simultaneamente per le anomalie di velocità P e S e per le correzioni della sorgente. Quando sono stati utilizzati i dati degli eventi situati nell’area di studio, abbiamo considerato quattro parametri incogniti corrispondenti agli spostamenti delle sorgenti nello spazio e nel tempo. Per i dati telesismici, abbiamo invertito per un parametro per evento per rappresentare l’incertezza della determinazione del tempo al di fuori del volume di studio. La matrice è stata invertita utilizzando il metodo LSQR38,39. La stabilità dell’inversione è stata controllata utilizzando equazioni aggiuntive che determinano l’ampiezza e l’appiattimento delle anomalie di velocità risultanti. I valori dei coefficienti di smorzamento sono stati impostati secondo diverse prove con modelli sintetici.
I risultati dell’inversione tomografica sono mostrati nell’articolo principale in una sezione orizzontale a 100 km di profondità (Fig. 2) e nei Materiali Supplementari in due sezioni orizzontali a 300 e 500 km di profondità (Fig. S2) e due sezioni verticali qui (Fig. S3). Qui, presentiamo i risultati per le anomalie di velocità delle onde P solo, perché i dati S sono quasi un decimo dei dati P, e il risultante S-modello non sembra sufficientemente stabile.
Nei Materiali Supplementari in Fig. S4, presentiamo i risultati del test checkerboard, dando le informazioni sulla risoluzione spaziale per i modelli recuperati. Il modello sintetico consisteva di anomalie positive e negative alternate con una magnitudine del 3% e dimensioni laterali di 5° × 5° km. Con l’aumentare della profondità, le anomalie hanno cambiato segno a 200, 400 e 600 km. I dati sintetici sono stati calcolati lungo gli stessi percorsi dei raggi per derivare il modello dei dati sperimentali, e questi sono stati perturbati da rumore casuale con una deviazione media di 0,5 s. Le anomalie periodiche a scacchiera sono state definite in tutta la Terra, mentre l’inversione è stata eseguita in regioni circolari selezionate. Questo ci ha permesso di esplorare l’effetto delle anomalie situate al di fuori dell’area di studio che sono state prese in considerazione durante il calcolo dei dati sintetici. I risultati del recupero a scacchiera sono mostrati nei materiali supplementari (Fig. S4). Le posizioni generali di tutte le anomalie sono state ricostruite correttamente; tuttavia, abbiamo osservato alcune sbavature diagonali legate agli orientamenti predominanti dei percorsi dei raggi. Abbiamo trovato una risoluzione verticale abbastanza buona, che ci ha permesso di recuperare chiaramente i cambiamenti di segno con la profondità.
Inoltre, abbiamo eseguito un test sintetico con forme realistiche di anomalie, che sono presentate in sezioni orizzontali e verticali (Figg. S5 e S6). Le anomalie sono definite all’interno di una serie di blocchi poligonali in alcuni intervalli di profondità. I risultati del recupero mostrano che le configurazioni laterali di tutte le anomalie sono ripristinate correttamente. Nelle sezioni verticali, vediamo che le anomalie che rappresentano la litosfera di spessore variabile sono risolte a profondità corrette. Entrambi questi test supportano l’affidabilità dei risultati derivati.
I risultati dell’inversione tomografica sono mostrati in tre sezioni orizzontali (Fig. 2 del documento principale e Fig. S2) e due sezioni verticali (Fig. S3). Oltre ai risultati relativi alla penisola indiana, il modello include alcune aree circostanti. Almeno due strutture sono state recuperate in modo coerente in diversi studi precedenti e potrebbero quindi essere utilizzate come un punto di riferimento naturale per il presente modello. Uno dei modelli più brillanti nella maggior parte degli studi di tomografia regionale asiatica è la ben studiata struttura ad alta Vp sotto il Pamir-Hindu Kush, che è associata alla distribuzione della sismicità a profondità intermedie (fino a 200 km). Le immagini di questa anomalia ad alta Vp sono state ottenute da diversi autori utilizzando diversi set di dati e algoritmi40,41,42. La seconda struttura di riferimento è un’alta anomalia Vp allungata e diretta da nord a sud sotto l’arco birmano contrassegnata dalla sismicità a profondità intermedia. Nel nostro modello, riveliamo questa anomalia come è stato riportato in studi precedenti43,44,45,46. Questi due esempi dimostrano fortemente che il nostro modello attuale è ugualmente stabile per aree che non sono state coperte da studi precedenti.
Modellazione della collisione continentale
Approccio alla modellazione
Il codice numerico termomeccanico visco-elasto-plastico 2D C-code I2ELVIS usato per modellare la collisione continentale è basato su un metodo delle differenze finite, applicato su una griglia euleriana sfalsata, e usa una tecnica marker-in-cell47,48. Le equazioni di conservazione di quantità di moto, massa ed energia sono risolte sulla griglia euleriana, e le proprietà fisiche sono trasportate da marcatori lagrangiani che si muovono secondo il campo di velocità interpolato dalla griglia. Nel modello viene utilizzata una reologia visco-elasto-plastica non newtoniana basata su leggi di flusso calibrate sperimentalmente (materiali supplementari, tabella S1). I dettagli completi di questo metodo, che permettono la sua riproduzione, sono forniti altrove47,48.
Progettazione modello numerico. La configurazione iniziale del modello (Materiali supplementari, Fig. S5) è 6.000 km di larghezza, 300 km di profondità, e risolto con una griglia rettangolare regolare di 601 × 151 a 1.201 × 151 nodi (varia in diversi esperimenti, Tabella S2) contenente 1,8 milioni di marcatori lagrangiani distribuiti in modo casuale. I confini superiore e destro del modello hanno condizioni meccaniche di libero scorrimento. Una velocità di convergenza costante di 4,7 cm/anno è prescritta al confine sinistro. La velocità del confine inferiore verso il basso è stata definita dalla condizione di conservazione del volume del dominio di calcolo, ed è stata quindi accorciata e ispessita ad ogni passo temporale. La condizione di superficie libera al di sopra della crosta è implementata usando uno strato d’aria “appiccicoso” di 20 km di spessore49,50 con bassa densità (1 kg/m3) e viscosità (1018 Pa-s). La struttura termica e litologica iniziale del modello (Fig. S5) è stata definita prescrivendo diverse unità litosferiche continentali corrispondenti alle principali regioni identificate all’interno delle placche indiane ed eurasiatiche e che differiscono nel gradiente termico iniziale della litosfera continentale (Fig. S5). In termini semplici, la crosta continentale inizialmente uniforme di 40 km di spessore è composta da 15 km di crosta felsica superiore, 10 km di crosta intermedia e 15 km di crosta mafica inferiore (lo spessore degli strati varia nei diversi esperimenti, Tabella S2). La struttura crostale inizialmente uniforme utilizzata è semplificata e trascura, per esempio, l’eterogeneità laterale dello spessore della crosta del continente indiano51,52. Questa semplificazione è dovuta principalmente alle grandi incertezze della nostra conoscenza dello spessore iniziale della crosta per le diverse unità litosferiche continentali. Questo spessore avrà probabilmente un effetto opposto rispetto allo spessore iniziale della litosfera a causa della debolezza reologica della crosta rispetto al mantello litosferico (Tabella S1). Lo spessore iniziale della crosta evolve fortemente durante gli esperimenti numerici, in cui la crosta si ispessisce prevalentemente nelle regioni di litosfera inizialmente sottile e quindi calda e debole (Tabella S2). Una zona debole su scala litosferica inclinata verso destra (una sutura di subduzione terminata dell’Oceano Tetide) segna il sito di inizio di una collisione India-Eurasia-like. I gradienti geotermici lineari semplificati sono stati utilizzati in diverse sezioni litosferiche (lo spessore varia nei diversi esperimenti, Tabella S2) a 273 K (in superficie) e 1.573 K (temperatura potenziale del mantello). Un gradiente termico adiabatico di 0.5 K/km è stato inizialmente prescritto nel mantello astenosferico. La conduttività termica dipendente dalla temperatura è stata usata per il mantello e la crosta (Tabella S1). Le condizioni termiche al contorno sono 273 K (in alto), 1.713 K (in basso), e zero flussi di calore sui confini sinistro e destro. Per assicurare un efficiente trasferimento di calore dalla superficie della crosta, la temperatura dell’aria/acqua “appiccicosa” è mantenuta costante a 273 K. Un’accelerazione gravitazionale di 9,81 m/s2 è stata usata per il modello. Va notato che il modello 2D usato nel nostro studio trascura la variabilità laterale del modello di deformazione 3D della zona di collisione India-Asia. Tuttavia pensiamo che questo modello semplificato sia sufficiente per lo scopo del nostro studio incentrato sul trasferimento della deformazione con il tempo da regioni litosferiche inizialmente più deboli a regioni inizialmente più forti.
Modello reologico visco-elasto-plastico
Le proprietà viscose, elastiche e fragili (plastiche) (Tabella S1) sono state implementate tramite la valutazione della viscosità effettiva del materiale. Per i materiali duttili, i contributi di diverse leggi di flusso come la dislocazione e lo scorrimento per diffusione sono stati considerati tramite il calcolo della viscosità duttile media inversa ηductile
dove ηnewt e ηpowl sono viscosità efficaci per la diffusione e lo scorrimento delle dislocazioni, rispettivamente, calcolate come
dove P è la pressione, T è la temperatura (in K), \({{dot{{varepsilon}_{II}=sqrt{1/2{({dot{varepsilon}{ij})}^{2}}}) è il secondo invariante del tensore della velocità di deformazione, σcr è lo stress di transizione diffusione-dislocazione36, e AD, E, V, e n sono parametri della legge di flusso determinati sperimentalmente (Tabella S1), che indicano rispettivamente la costante del materiale, l’energia di attivazione, il volume di attivazione e l’esponente di stress.
La reologia duttile è combinata con una reologia fragile (plastica) per produrre una reologia viscoso-plastica effettiva utilizzando il seguente limite superiore per la viscosità duttile:
dove P è la pressione, ϕ è il coefficiente di attrito interno (Tabella S1), e C è la resistenza alla trazione della roccia a P = 0 (Tabella S1). L’elasticità è implementata sulla base del modello visco-elasto-plastico incomprimibile di Maxwell47,48. I moduli di taglio (μ) dei diversi materiali sono specificati nella Tabella S1.
Risultati numerici
Abbiamo eseguito 12 esperimenti numerici, variando la stratificazione crostale e le lunghezze iniziali e gli spessori litosferici delle diverse sezioni del modello (Tabella S2). I risultati numerici mostrano una migrazione sistematica delle deformazioni da sezioni litosferiche inizialmente più deboli (cioè più sottili) a sezioni inizialmente più forti (cioè più spesse) associate alla crescita graduale delle tensioni compressive intra-placca sia nel dominio del modello indiano che eurasiatico. Questa tendenza generale non è influenzata dalle variazioni della geometria iniziale del modello, che sono state esplorate; essa influenza solo le dinamiche di deformazione nelle sezioni litosferiche più deboli (Tibet, Tien-Shan) (Tabella S2).