Datos y algoritmos para la tomografía
Utilizamos el algoritmo de tomografía regional34 con datos de tiempo de recorrido global del Centro Sismológico Internacional35 correspondientes al periodo 1964-2014. Para la región seleccionada, consideramos cualquier dato correspondiente a las trayectorias de los rayos que pasan por el volumen de estudio. Esto incluye los rayos de los terremotos localizados en el área de estudio registrados por estaciones de todo el mundo, y los de los eventos telesísmicos registrados por estaciones localizadas en el área de estudio (Fig. S1A y B). Antes de su uso en la tomografía, los datos fueron reprocesados; esto incluyó la reubicación de las fuentes y el rechazo de los valores atípicos36. Para localizar los eventos, se utilizó el modelo de velocidad unidimensional AK13537.
Esta área ha sido previamente parte de otro modelo calculado utilizando el mismo algoritmo25. Sin embargo, el estudio anterior sólo abarcaba la parte norte de la India. Además, nuestro estudio incluye datos de 2005 a 2014, que no estaban disponibles en el estudio anterior. Diez años de registros adicionales han proporcionado una cantidad significativa de datos, especialmente los correspondientes a las nuevas estaciones de la India que mejoraron drásticamente la cobertura de rayos.
La inversión se realizó por separado en una serie de zonas superpuestas que cubrían toda la región de estudio. Utilizamos tres regiones, cada una con un radio de 8° (Materiales Suplementarios, Fig. S1). Definimos la profundidad del volumen de estudio en 1.000 km; sin embargo, consideramos principalmente los resultados hasta una profundidad de 700 km porque las estructuras más profundas podrían haber sido afectadas por anomalías situadas fuera del área de estudio. La parametrización de la distribución de la velocidad se realizó utilizando un conjunto de nodos distribuidos en niveles horizontales a profundidades de 25, 50, 75, 100, 150, 220, 290, 360, 430, 500, 570, 640, 710, 800 y 900 km. En cada nivel de profundidad, los nodos se distribuyeron en función de la densidad de rayos; una cobertura de rayos más densa correspondía a un espaciado de nodos menor. El espaciado mínimo se fijó en 30 km. Para evitar artefactos relacionados con la geometría de la cuadrícula, realizamos los cálculos para dos cuadrículas diferentes con orientaciones básicas de 0 y 45° y luego promediamos los resultados.
La inversión se realizó simultáneamente para las anomalías de velocidad P y S y para las correcciones de fuentes. Cuando se utilizaron los datos de los eventos localizados en el área de estudio, se consideraron cuatro parámetros desconocidos correspondientes a los desplazamientos de las fuentes en el espacio y el tiempo. Para los datos telesísmicos, invertimos para un parámetro por evento para representar la incertidumbre de la determinación del tiempo fuera del volumen de estudio. La matriz se invirtió utilizando el método LSQR38,39. La estabilidad de la inversión se controló mediante ecuaciones adicionales que determinan la amplitud y el aplanamiento de las anomalías de velocidad resultantes. Los valores de los coeficientes de amortiguación se fijaron según varias pruebas con modelos sintéticos.
Los resultados de la inversión tomográfica se muestran en el artículo principal en una sección horizontal a 100 km de profundidad (Fig. 2) y en los Materiales Complementarios en dos secciones horizontales a 300 y 500 km de profundidad (Fig. S2) y dos secciones verticales aquí (Fig. S3). Aquí, presentamos los resultados para las anomalías de velocidad de las ondas P solamente, porque los datos S son casi una décima parte de los datos P, y el modelo S resultante no parece suficientemente estable.
En los Materiales Suplementarios en la Fig. S4, presentamos los resultados de la prueba de tablero de ajedrez, dando la información sobre la resolución espacial para los modelos recuperados. El modelo sintético consistía en alternar anomalías positivas y negativas con una magnitud del 3% y un tamaño lateral de 5° × 5° km. Al aumentar la profundidad, las anomalías cambiaron de signo a 200, 400 y 600 km. Los datos sintéticos se calcularon a lo largo de las mismas trayectorias de rayos para derivar el modelo de datos experimentales, y éstos fueron perturbados por ruido aleatorio con una desviación media de 0,5 s. Las anomalías periódicas en damero se definieron en toda la Tierra, mientras que la inversión se realizó en regiones circulares seleccionadas. Esto nos permitió explorar el efecto de las anomalías situadas fuera del área de estudio que se tuvieron en cuenta al calcular los datos sintéticos. Los resultados de la recuperación del tablero de control se muestran en los Materiales Suplementarios (Fig. S4). Las localizaciones generales de todas las anomalías fueron reconstruidas correctamente; sin embargo, observamos algunas manchas diagonales relacionadas con las orientaciones predominantes de las trayectorias de los rayos. Encontramos una resolución vertical bastante buena, lo que nos permitió recuperar claramente los cambios de signo con la profundidad.
Además, hemos realizado una prueba sintética con formas realistas de las anomalías, que se presentan en secciones horizontales y verticales (Figs S5 y S6). Las anomalías se definen dentro de una serie de bloques poligonales en algunos intervalos de profundidad. Los resultados de la recuperación muestran que las configuraciones laterales de todas las anomalías se restauran correctamente. En las secciones verticales, vemos que las anomalías que representan la litosfera de espesor variable se resuelven a profundidades correctas. Ambas pruebas apoyan la fiabilidad de los resultados derivados.
Los resultados de la inversión tomográfica se muestran en tres secciones horizontales (Fig. 2 del artículo principal y Fig. S2) y dos secciones verticales (Fig. S3). Aparte de los resultados relacionados con la península de la India, el modelo incluye algunas áreas circundantes. Al menos dos estructuras fueron recuperadas de forma consistente en varios estudios anteriores y, por lo tanto, podrían ser utilizadas como un punto de referencia natural para el presente modelo. Uno de los patrones más brillantes en la mayoría de los estudios de tomografía regional de Asia es la bien estudiada estructura de alta Vp bajo el Pamir-Hindu Kush, que está asociada a la distribución de la sismicidad a profundidades intermedias (hasta 200 km). Las imágenes de esta anomalía de alto Vp han sido obtenidas de forma consistente por diferentes autores utilizando diferentes conjuntos de datos y algoritmos40,41,42. La segunda estructura de referencia es una anomalía alargada de alto Vp dirigida de norte a sur bajo el arco birmano marcada con la sismicidad de profundidad intermedia. En nuestro modelo, revelamos esta anomalía tal y como se informó en estudios anteriores43,44,45,46. Estos dos ejemplos demuestran fuertemente que nuestro modelo actual es igualmente estable para áreas que no fueron cubiertas por estudios anteriores.
Modelación de la colisión continental
Enfoque de la modelación
El código C numérico termomecánico visco-elasto-plástico 2D I2ELVIS utilizado para modelar la colisión continental se basa en un método de diferencias finitas, aplicado en una malla euleriana escalonada, y utiliza una técnica de marcador en celda47,48. Las ecuaciones de conservación del momento, la masa y la energía se resuelven en la malla euleriana, y las propiedades físicas son transportadas por marcadores lagrangianos que se mueven según el campo de velocidad interpolado desde la malla. En el modelo se utiliza una reología viscoelástica no newtoniana basada en leyes de flujo calibradas experimentalmente (Materiales Suplementarios, Tabla S1). Los detalles completos de este método, que permiten su reproducción, se proporcionan en otra parte47,48.
Diseño del modelo numérico. La configuración inicial del modelo (Materiales Suplementarios, Fig. S5) es de 6.000 km de ancho, 300 km de profundidad, y se resuelve con una malla rectangular regular de 601 × 151 a 1.201 × 151 nodos (variada en diferentes experimentos, Tabla S2) que contiene 1,8 millones de marcadores lagrangianos distribuidos aleatoriamente. Los límites superior y derecho del modelo tienen condiciones de contorno mecánicas de deslizamiento libre. En el límite izquierdo se prescribe una velocidad de convergencia constante de 4,7 cm/año. La velocidad del límite inferior hacia abajo fue definida por la condición de conservación del volumen del dominio computacional, por lo que se acortó y engrosó en cada paso de tiempo. La condición de frontera de superficie libre por encima de la corteza se implementa utilizando una capa de aire «pegajosa» de 20 km de espesor49,50 con baja densidad (1 kg/m3) y viscosidad (1018 Pa-s). La estructura térmica y litológica inicial del modelo (Fig. S5) se definió prescribiendo varias unidades litosféricas continentales correspondientes a las principales regiones identificadas dentro de las placas India y Euroasiática y que difieren en el gradiente térmico inicial de la litosfera continental (Fig. S5). En términos sencillos, la corteza continental inicialmente uniforme de 40 km de espesor se compone de 15 km de corteza superior félsica, 10 km de corteza media intermedia y 15 km de corteza inferior máfica (el espesor de las capas varía en diferentes experimentos, Tabla S2). La estructura de la corteza inicialmente uniforme utilizada está simplificada y descuida, por ejemplo, la heterogeneidad lateral del espesor de la corteza del continente indio51,52. Esta simplificación se debe principalmente a las grandes incertidumbres de nuestro conocimiento del espesor inicial de la corteza para las diferentes unidades litosféricas continentales. Es probable que este espesor tenga un efecto opuesto al del espesor inicial de la litosfera debido a la debilidad reológica de la corteza en comparación con el manto litosférico (Tabla S1). El espesor inicial de la corteza evoluciona fuertemente durante los experimentos numéricos, en los que la corteza se engrosa predominantemente en regiones de litosfera inicialmente delgada y, por tanto, cálida y débil (Tabla S2). Una zona débil a escala litosférica inclinada hacia la derecha (una sutura de subducción terminada del Océano Tethys) marca el lugar de inicio de una colisión tipo India-Eurasia. Se utilizaron gradientes geotérmicos lineales simplificados en diferentes secciones litosféricas (los espesores variaron en los diferentes experimentos, Tabla S2) a 273 K (en la superficie) y 1.573 K (temperatura potencial del manto). En el manto astenosférico se prescribió inicialmente un gradiente térmico adiabático de 0,5 K/km. Se utilizó una conductividad térmica dependiente de la temperatura para el manto y la corteza (Tabla S1). Las condiciones térmicas de contorno son 273 K (superior), 1.713 K (inferior), y flujos de calor nulos en los límites izquierdo y derecho. Para asegurar una transferencia de calor eficiente desde la superficie de la corteza, la temperatura del aire/agua «pegajosa» se mantiene constante en 273 K. Se utilizó una aceleración gravitacional de 9,81 m/s2 para el modelo. Cabe señalar que el modelo 2D utilizado en nuestro estudio no tiene en cuenta la variabilidad lateral del patrón de deformación 3D de la zona de colisión India-Asia. Sin embargo, pensamos que este modelo simplificado es suficiente para el propósito de nuestro estudio centrado en la transferencia de la deformación con el tiempo desde las regiones litosféricas inicialmente más débiles a las inicialmente más fuertes.
Modelo reológico disco-elasto-plástico
Las propiedades viscosas, elásticas y frágiles (plásticas) (Tabla S1) se implementaron mediante la evaluación de la viscosidad efectiva del material. Para los materiales dúctiles las contribuciones de las diferentes leyes de flujo, como la dislocación y la fluencia por difusión, se consideraron mediante el cálculo de la viscosidad dúctil media inversa ηductile
donde ηnewt y ηpowl son viscosidades efectivas para la difusión y la fluencia por dislocación, respectivamente, calculadas como
donde P es la presión, T es la temperatura (en K), \({\dot{varepsilon }}_{II}=\sqrt{1/2{{\dot{varepsilon }}_{ij})}^{2}}) es la segunda invariante del tensor de velocidad de deformación, σcr es la tensión de transición difusión-dislocación36, y AD, E, V y n son parámetros de la ley de flujo determinados experimentalmente (Tabla S1), que denotan la constante del material, la energía de activación, el volumen de activación y el exponente de tensión, respectivamente.
La reología dúctil se combina con una reología frágil (plástica) para obtener una reología viscoso-plástica efectiva utilizando el siguiente límite superior para la viscosidad dúctil:
donde P es la presión, ϕ es el coeficiente de fricción interna (Tabla S1), y C es la resistencia a la tracción de la roca en P = 0 (Tabla S1). La elasticidad se implementa en base al modelo visco-elasto-plástico incompresible de Maxwell47,48. Los módulos de corte (μ) de los diferentes materiales se especifican en la Tabla S1.
Resultados numéricos
Realizamos 12 experimentos numéricos, variando la estratificación de la corteza y las longitudes iniciales y espesores litosféricos de diferentes secciones del modelo (Tabla S2). Los resultados numéricos muestran una migración sistemática de las deformaciones desde las secciones litosféricas inicialmente más débiles (es decir, más delgadas) a las inicialmente más fuertes (es decir, más gruesas), asociada al crecimiento gradual de las tensiones de compresión intraplaca tanto en los dominios del modelo indio como euroasiático. Esta tendencia general no se ve afectada por las variaciones de la geometría inicial del modelo, que han sido exploradas; sólo influye en la dinámica de deformación en las secciones litosféricas más débiles (Tíbet, Tien-Shan) (Tabla S2).