Données et algorithmes pour la tomographie
Nous avons utilisé l’algorithme de tomographie régionale34 avec les données de temps de parcours global du Centre sismologique international35 correspondant à la période 1964-2014. Pour la région sélectionnée, nous considérons toutes les données correspondant aux trajets des rayons passant par le volume d’étude. Cela inclut les rayons provenant de séismes situés dans la zone d’étude et enregistrés par des stations du monde entier, ainsi que ceux provenant d’événements télésismiques enregistrés par des stations situées dans la zone d’étude (Fig. S1A et B). Avant d’être utilisées pour la tomographie, les données ont été retraitées, ce qui comprenait la relocalisation des sources et le rejet des valeurs aberrantes36. Pour localiser les événements, nous avons utilisé le modèle de vitesse unidimensionnel AK13537.
Cette zone a déjà fait partie d’un autre modèle calculé à l’aide du même algorithme25. Cependant, l’étude précédente ne couvrait que la partie nord de l’Inde. De plus, notre étude inclut des données de 2005 à 2014, qui n’étaient pas disponibles dans l’étude précédente. Dix années d’enregistrements supplémentaires ont fourni une quantité importante de données, en particulier celles correspondant aux nouvelles stations en Inde qui ont considérablement amélioré la couverture des rayons.
L’inversion a été réalisée séparément dans une série de zones se chevauchant et couvrant toute la région d’étude. Nous avons utilisé trois régions, chacune avec un rayon de 8° (matériaux supplémentaires, figure S1). Nous avons défini la profondeur du volume d’étude à 1 000 km ; cependant, nous avons surtout considéré les résultats jusqu’à une profondeur de 700 km car les structures plus profondes ont pu être affectées par des anomalies situées en dehors de la zone d’étude. La paramétrisation de la distribution des vitesses a été réalisée en utilisant un ensemble de nœuds répartis sur des niveaux horizontaux à des profondeurs de 25, 50, 75, 100, 150, 220, 290, 360, 430, 500, 570, 640, 710, 800 et 900 km. À chaque niveau de profondeur, les nœuds ont été répartis en fonction de la densité des rayons ; une couverture de rayons plus dense correspondait à un espacement plus petit des nœuds. L’espacement minimal a été fixé à 30 km. Pour éviter les artefacts liés à la géométrie de la grille, nous avons effectué les calculs pour deux grilles différentes avec des orientations de base de 0 et 45°, puis nous avons calculé la moyenne des résultats.
L’inversion a été réalisée simultanément pour les anomalies de vitesse P et S et pour les corrections de source. Lorsque les données des événements situés dans la zone d’étude ont été utilisées, nous avons considéré quatre paramètres inconnus correspondant aux décalages des sources dans l’espace et dans le temps. Pour les données télésismiques, nous avons inversé pour un paramètre par événement pour représenter l’incertitude de la détermination du temps en dehors du volume d’étude. La matrice a été inversée en utilisant la méthode LSQR38,39. La stabilité de l’inversion a été contrôlée à l’aide d’équations supplémentaires déterminant l’amplitude et l’aplatissement des anomalies de vitesse résultantes. Les valeurs des coefficients d’amortissement ont été fixées selon plusieurs essais avec des modèles synthétiques.
Les résultats de l’inversion tomographique sont présentés dans l’article principal dans une section horizontale à 100 km de profondeur (Fig. 2) et dans les matériaux supplémentaires dans deux sections horizontales à 300 et 500 km de profondeur (Fig. S2) et deux sections verticales ici (Fig. S3). Ici, nous présentons les résultats pour les anomalies de vitesse des ondes P uniquement, car les données S représentent presque un dixième des données P, et le modèle S résultant ne semble pas suffisamment stable.
Dans les matériaux supplémentaires, à la Fig. S4, nous présentons les résultats du test du damier, donnant les informations sur la résolution spatiale pour les modèles récupérés. Le modèle synthétique consistait en une alternance d’anomalies positives et négatives avec une magnitude de 3% et une taille latérale de 5° × 5° km. Avec l’augmentation des profondeurs, les anomalies ont changé de signe à 200, 400 et 600 km. Les données synthétiques ont été calculées le long des mêmes trajets de rayons que ceux utilisés pour obtenir le modèle de données expérimentales, et elles ont été perturbées par un bruit aléatoire avec une déviation moyenne de 0,5 s. Les anomalies périodiques en damier ont été définies sur l’ensemble de la Terre, tandis que l’inversion a été réalisée dans des régions circulaires sélectionnées. Cela nous a permis d’explorer l’effet des anomalies situées en dehors de la zone d’étude qui ont été prises en compte lors du calcul des données synthétiques. Les résultats de la récupération en damier sont présentés dans les Matériaux supplémentaires (Fig. S4). Les emplacements généraux de toutes les anomalies ont été reconstruits correctement ; cependant, nous avons observé un certain étalement diagonal lié aux orientations prédominantes des trajectoires des rayons. Nous avons constaté une assez bonne résolution verticale, nous permettant de récupérer clairement les changements de signe avec la profondeur.
En outre, nous avons effectué un test synthétique avec des formes réalistes d’anomalies, qui sont présentées en sections horizontales et verticales (Figs S5 et S6). Les anomalies sont définies dans une série de blocs polygonaux dans certains intervalles de profondeur. Les résultats de la récupération montrent que les configurations latérales de toutes les anomalies sont restaurées correctement. Dans les sections verticales, on voit que les anomalies représentant la lithosphère d’épaisseur variable sont résolues à des profondeurs correctes. Ces deux tests soutiennent la fiabilité des résultats dérivés.
Les résultats de l’inversion tomographique sont présentés dans trois sections horizontales (Fig. 2 de l’article principal et Fig. S2) et deux sections verticales (Fig. S3). Outre les résultats liés à la péninsule indienne, le modèle inclut certaines zones environnantes. Au moins deux structures ont été retrouvées de manière cohérente dans plusieurs études précédentes et pourraient donc être utilisées comme référence naturelle pour le présent modèle. L’une des structures les plus brillantes dans la plupart des études tomographiques régionales asiatiques est la structure à haute Vp bien étudiée sous le Pamir-Hindu Kush, qui est associée à la distribution de la sismicité à des profondeurs intermédiaires (jusqu’à 200 km). Les images de cette anomalie à haute Vp ont été obtenues de manière cohérente par différents auteurs utilisant différents ensembles de données et algorithmes40,41,42. La deuxième structure de référence est une anomalie haute Vp allongée, orientée nord-sud, située sous l’arc birman et marquée par la sismicité de profondeur intermédiaire. Dans notre modèle, nous révélons cette anomalie comme cela a été rapporté dans des études précédentes43,44,45,46. Ces deux exemples montrent fortement que notre modèle actuel est également stable pour les zones qui n’étaient pas couvertes par les études précédentes.
Modélisation de la collision continentale
Approche de modélisation
Le code numérique thermomécanique visco-élasto-plastique 2D C-code I2ELVIS utilisé pour modéliser la collision continentale est basé sur une méthode de différences finies, appliquée sur une grille eulérienne décalée, et utilise une technique de marqueur dans la cellule47,48. Les équations de conservation de la quantité de mouvement, de la masse et de l’énergie sont résolues sur la grille eulérienne, et les propriétés physiques sont transportées par des marqueurs lagrangiens qui se déplacent selon le champ de vitesse interpolé à partir de la grille. Une rhéologie visco-élasto-plastique non-newtonienne basée sur des lois d’écoulement calibrées expérimentalement est utilisée dans le modèle (Matériaux supplémentaires, Tableau S1). Les détails complets de cette méthode, qui permettent sa reproduction, sont fournis ailleurs47,48.
Conception numérique du modèle. La configuration initiale du modèle (matériaux supplémentaires, figure S5) a une largeur de 6 000 km, une profondeur de 300 km et est résolue avec une grille rectangulaire régulière de 601 × 151 à 1 201 × 151 nœuds (variable dans différentes expériences, tableau S2) contenant 1,8 million de marqueurs lagrangiens distribués de façon aléatoire. Les limites supérieure et droite du modèle ont des conditions de limites mécaniques de glissement libre. Une vitesse de convergence constante de 4,7 cm/an est prescrite à la frontière gauche. La vitesse de la frontière inférieure descendante est définie par la condition de conservation du volume du domaine de calcul, et est donc raccourcie et épaissie à chaque pas de temps. La condition de surface libre au-dessus de la croûte est implémentée en utilisant une couche d’air « collante » de 20 km d’épaisseur49,50 avec une faible densité (1 kg/m3) et une faible viscosité (1018 Pa-s). La structure thermique et lithologique initiale du modèle (Fig. S5) a été définie en prescrivant plusieurs unités lithosphériques continentales correspondant aux principales régions identifiées au sein des plaques indienne et eurasienne et différant par le gradient thermique initial de la lithosphère continentale (Fig. S5). En termes simples, la croûte continentale initialement uniforme de 40 km d’épaisseur est composée de 15 km de croûte supérieure felsique, de 10 km de croûte moyenne intermédiaire et de 15 km de croûte inférieure mafique (l’épaisseur des couches varie dans les différentes expériences, tableau S2). La structure crustale initialement uniforme utilisée est simplifiée et néglige, par exemple, l’hétérogénéité latérale de l’épaisseur de la croûte du continent indien51,52. Cette simplification est principalement due aux grandes incertitudes de notre connaissance de l’épaisseur crustale initiale pour les différentes unités lithosphériques continentales. Cette épaisseur aura probablement un effet opposé à celui de l’épaisseur initiale de la lithosphère en raison de la faiblesse rhéologique de la croûte par rapport au manteau lithosphérique (Tableau S1). L’épaisseur initiale de la croûte évolue fortement au cours des expériences numériques, dans lesquelles la croûte s’épaissit de manière prédominante dans les régions de lithosphère initialement mince et donc chaude et faible (tableau S2). Une zone faible d’échelle lithosphérique inclinée vers la droite (une suture de subduction terminée de l’océan Téthys) marque le site d’initiation d’une collision de type Inde-Eurasie. Des gradients géothermiques linéaires simplifiés ont été utilisés dans différentes sections lithosphériques (les épaisseurs varient selon les expériences, tableau S2) à 273 K (à la surface) et 1 573 K (température potentielle du manteau). Un gradient thermique adiabatique de 0,5 K/km a été initialement prescrit dans le manteau asthénosphérique. La conductivité thermique dépendant de la température a été utilisée pour le manteau et la croûte (Tableau S1). Les conditions aux limites thermiques sont 273 K (supérieur), 1 713 K (inférieur), et des flux de chaleur nuls sur les limites gauche et droite. Pour assurer un transfert de chaleur efficace depuis la surface de la croûte, la température de l’air/eau » collante » est maintenue constante à 273 K. Une accélération gravitationnelle de 9,81 m/s2 a été utilisée pour le modèle. Il convient de noter que le modèle 2D utilisé dans notre étude néglige la variabilité latérale du modèle de déformation 3D de la zone de collision Inde-Asie. Cependant, nous pensons que ce modèle simplifié est suffisant pour l’objectif de notre étude axée sur le transfert de la déformation avec le temps des régions lithosphériques initialement plus faibles aux régions lithosphériques initialement plus fortes.
Modèle rhéologique Visco-élasto-plastique
Les propriétés visqueuses, élastiques et fragiles (plastiques) (tableau S1) ont été implémentées via l’évaluation de la viscosité effective du matériau. Pour les matériaux ductiles les contributions des différentes lois d’écoulement telles que le fluage par dislocation et diffusion ont été prises en compte via le calcul de la viscosité ductile moyenne inverse ηductile
où ηnewt et ηpowl sont des viscosités effectives pour la diffusion et le fluage des dislocations, respectivement, calculées comme suit
où P est la pression, T est la température (en K), \({\dot{\varepsilon }}_{II}=\sqrt{1/2{({\dot{\varepsilon }}_{ij})}^{2}}\) est le deuxième invariant du tenseur de vitesse de déformation, σcr est la contrainte de transition diffusion-dislocation36, et AD, E, V et n sont des paramètres de loi d’écoulement déterminés expérimentalement (tableau S1), qui désignent respectivement la constante de matière, l’énergie d’activation, le volume d’activation et l’exposant de contrainte.
La rhéologie ductile est combinée avec une rhéologie fragile (plastique) pour donner une rhéologie visqueuse-plastique effective en utilisant la limite supérieure suivante pour la viscosité ductile :
où P est la pression, ϕ est le coefficient de friction interne (tableau S1), et C est la résistance à la traction de la roche à P = 0 (tableau S1). L’élasticité est mise en œuvre sur la base du modèle visco-élasto-plastique incompressible de Maxwell47,48. Les modules de cisaillement (μ) des différents matériaux sont précisés dans le tableau S1.
Résultats numériques
Nous avons réalisé 12 expériences numériques, en faisant varier la stratification crustale et les longueurs initiales et épaisseurs lithosphériques des différentes sections du modèle (tableau S2). Les résultats numériques montrent une migration systématique des déformations des sections lithosphériques initialement plus faibles (c’est-à-dire plus minces) vers les sections lithosphériques initialement plus fortes (c’est-à-dire plus épaisses), associée à la croissance progressive des contraintes de compression intra-plaque dans les domaines modèles indien et eurasien. Cette tendance générale n’est pas affectée par les variations de la géométrie initiale du modèle, qui ont été explorées ; elle n’influence la dynamique de déformation que dans les sections lithosphériques les plus faibles (Tibet, Tien-Shan) (tableau S2).