The role of continental lithospheric thermal structure in the evolution of orogenic systems: application to the Himalayan–Tibetan collision zone

Continental collision is a crucial process in plate tectonics. However, in terms of the evolution and the controlling parameters of its lateral heterogeneity, our understanding of the tectonic complexities at such a convergent plate boundary remains largely unclear. In this study, we conducted a series of two-dimensional numerical experiments to investigate how continental lithospheric thermal structure influences the development of lateral heterogeneity along the continental collision zone. The following two end-members were achieved. First, continuous subduction mode, which prevails when the model has a cold procontinental Moho temperature (≤450 ^∘C). In this case, a narrow collision orogen develops, and the subducting angle steepens with the increasing retrocontinental Moho temperature. Second, continental subduction with a slab break-off, which generates a relative wide collision orogen and dominates when the model has a relatively hot procontinental Moho temperature (≥500 ^∘C), especially when the Moho temperature ≥ 550 ^∘C. Radioactive heat production is the second-order controlling parameter in varying the continental collision mode, while it prefers to enhance strain localization in the upper part of the continental lithosphere and promote the growth of shear zones there. By comparing the model results with geological observations, we suggest that the discrepant evolutionary paths from the continuous subduction underlying the Hindu Kush to the continental subduction after slab break-off beneath eastern Tibet may originate from the inherited lateral inhomogeneity of the Indian lithospheric thermal structure. Besides, the high content of crustal radioactive elements may be one of the most important factors that controls the formation of large thrust fault zones in the Himalayas.

1 Introduction

The collision between the Indian and Asian continents, following the closure of the Neo-Tethys Ocean, is regarded as one of the key geodynamical processes of plate tectonics (Toussaint et al., 2004), which creates the world's largest orogen (the Alpine–Himalaya orogen) and highest plateau (the Tibetan Plateau). It would be really hard to maintain a homogeneous plate morphology and the crustal–mantle deformation features during such a large-scale collision stretching thousands of kilometers. Many of the previous observations have already provided insights into the significant lateral heterogeneity of Himalayan–Tibetan orogen (Zhou and Murphy, 2005; Li et al., 2008; Chen et al., 2015). Specifically, the horizontal underthrusting distance of the Indian lithosphere decreases, and the subducting angle increases laterally from west to east. Despite various mechanisms that have been invoked, including the geometry of the subducting plate (Nettesheim et al., 2018; Koptev et al., 2019), inhomogeneous inherit lithospheric structure, variations in rock composition, characters of terrane assembly, and the complex tectonic settings (Beaumont et al., 1994; Chen et al., 2016, 2019; Vogt et al., 2018; Liu and Yang, 2022), the development of the collision zone's lateral heterogeneity remains ambiguous and requires further investigation.

Many of the laboratory (Willingshofer and Sokoutis, 2009; Luth et al., 2010) and numerical modeling experiments (Beaumont et al., 1994; Pysklywec, 2001; Toussaint et al., 2004; Huangfu et al., 2017; Liao and Gerya, 2017; Vogt et al., 2018; S. Liu et al., 2022) have suggested that the lithospheric thermal structure is one of the important controlling parameters contributing to the regulation the evolution of continental collision (Pysklywec, 2001; Faccenda et al., 2008; Ueda et al., 2012; Tang et al., 2020; M. Liu et al., 2022). Pysklywec (2001) once used a series of thermomechanical models that focused on the evolutionary modes of the continental lithospheric mantle during collision. The experiments suggested that an increase in the lithospheric geotherm resulted in a transition of the lithospheric mantle from an asymmetric, subduction-like mode to ablative subduction. Similarly, on the basis of two-dimensional mechanical models, Toussaint et al. (2004) investigated the potential scenarios for the evolution of continental collision. They concluded that models with a relatively low (<550 ^∘C) Moho temperature and a fast convergence rate (>5 cm yr⁻¹) are more likely to maintain stable subduction. On the contrary, among warmer and slower convergent models, lithospheric folding, pure-shear thickening, and Rayleigh–Taylor (RT) instabilities prevailed. Ghazian and Buiter (2013) used 2D numerical models to discuss the sensitivity of various collision modes to velocity, crust–mantle temperature structure, lithospheric rheology, and the density difference between the lithospheric mantle and asthenosphere. They recognized that velocity, lithospheric rheology, and temperature have important controls. More precisely, stable continental subduction evolves over a large range of values for convergent velocity and lithospheric temperature. Fast and cold models tend to fold, while slow and warm ones can generate RT-type dripping. Heron and Pysklywec (2016) subsequently presented dynamic thermomechanical experiments to explore the role of lithospheric inherited weakness in the orogenic process. According to their model results, lithospheric temperature exerts a strong control on lithospheric strength; it may influence the brittle–ductile transition depth to help determine which layer (upper crust, lower crust, or lithospheric mantle) dominantly controls the lithospheric deformation style in a collision system. Recently, Huangfu et al. (2019) employed systematic numerical simulations to study the dynamics of the lithospheric-scale subduction and crustal-scale underthrusting of the continental collision zone. They found that models with low convergence velocity and a cold upper plate are inclined to lithospheric-scale subduction, whereas models with high velocity and a hot overriding plate contribute to crustal-scale underthrusting. Though numerous numerical simulations are employed to qualify the effects of lithospheric thermal structure, it is noteworthy, however, that earlier studies consistently and mainly focused on the influences of the overriding lithospheric thermal structure. A fairly large number of them neglected oceanic subduction prior to the continental collision for simplicity. Thus, our understanding of tectonic complexities that may emerge during progressive crust–mantle deformation at the convergent plate boundary remains largely unclear.

Here we present a series of high-resolution numerical experiments to gain new insights to how continental lithospheric thermal structure influences the evolutionary path of the continental collision system. Temperature heterogeneities are incorporated by altering the continental Moho temperature (T_moho) and crustal radioactive heat production (H_r), respectively. In the end, the results of our models are used to draw comparisons with some first-order characters of the Indian–Asian collision zone.

2 Materials and methods

2.1 Governing equation

We use the parallel code Advanced Solver for Problems in Earth's ConvecTion (ASPECT; Kronbichler et al., 2012), an extensible code of the C$++$ program library deal.ii (Differential Equations Analysis Library; https://www.dealii.org/, last access: 24 June 2023), which is targeted at the computational solution of partial differential equations, using adaptive finite elements (Arndt et al., 2021), to solve three conservation equations. The methodology is similar to our recent work (Liu and Yang, 2022; M. Liu et al., 2022). The models are calculated by solving the following equations of conservation of mass, momentum, and internal energy for an incompressible medium and adopting the extended Boussinesq approximation (van Zelst et al., 2022).

$$\begin{array}{}\text{(1)}& {\displaystyle }-\mathrm{\nabla }\cdot \mathrm{2}\mathit{\eta }\dot{\mathit{\epsilon }}\left(u\right)+\mathrm{\nabla }p=\mathit{\rho }g\mathrm{in}\mathrm{\Omega },\text{(2)}& {\displaystyle }\mathrm{\nabla }\cdot u=\mathrm{0}\mathrm{in}\mathrm{\Omega },\text{(3)}& {\displaystyle }\begin{array}{rl}{\mathit{\rho }}_{\mathrm{0}}{C}_p\left({\displaystyle \frac{\partial T}{\partial t}}+u\cdot \mathrm{\nabla }T\right)\mathrm{\nabla }& \cdot k\mathrm{\nabla }T={\mathit{\rho }}_{\mathrm{0}}H+\mathrm{2}\mathit{\eta }\dot{\mathit{\epsilon }}\left(u\right):\dot{\mathit{\epsilon }}\left(u\right)\\ & +\mathit{\alpha }T(u\cdot \mathrm{\nabla }P)\mathrm{in}\mathrm{\Omega }\end{array}\text{(4)}& {\displaystyle }{\displaystyle \frac{\partial c_i}{\partial t}}+u\cdot \mathrm{\nabla }{c}_i=\mathrm{0}\end{array}$$

On the right-hand side of Eq. (3), the three terms represent the internal heat production that includes radioactive decay, frictional heating, and adiabatic compression of material. η is the viscosity, $\dot{\mathit{\epsilon }}\left(u\right)=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla }u+\mathrm{\nabla }{u}^T)$ is the deviator of the strain rate tensor, u is the velocity, p is the pressure, g is the gravitational acceleration, Ω is the interesting domain, C_p is the heat capacity, k is the heat conductivity, H is the intrinsic specific heat production, α is the thermal expansion coefficient, and ρ₀ is the adiabatic reference density. c_i is named compositional fields here (e.g., upper crust and lower crust; Kronbichler et al., 2012). As the medium was considered incompressible, we assumed that the density (ρ) in Eq. (1) satisfied the following equation:

$$\begin{array}{}\text{(5)}& \mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{1}-\mathit{\alpha }\left(T-T_{\mathrm{0}}\right)\right),\end{array}$$

where ρ₀ is the reference density at reference temperature T₀ (293 K).

Figure 1 (a) Reference model configuration and boundary conditions. Different colors reflect different lithologies, where 1 and 8 are the continental upper crust, 2 and 9 are the continental lower crust, 3 and 10 are the continental lithospheric mantle, 4 is the weak zone, 5 is the sediment, 6 is the oceanic crust, 7 is the oceanic lithospheric mantle, and 11 is the sub-lithospheric mantle (Table S1). The white lines are isotherms of 200 ^∘C increments in the range of 200 to 1200 ^∘C. T_surf is the surface temperature, and T_lab is the bottom temperature of continental lithosphere. V_in and V_out denote where material flows in and out. (b) Initial distributions of procontinental (T_{cp_pro}), retrocontinental (T_{cp_retro}), and oceanic (T_op) lithospheric temperatures.

2.2 Rheology

Our models are based on the common assumption that solid Earth materials are treated as highly viscous fluids. Thus, we apply a viscoplastic rheology (Glerum et al., 2018). The viscous regime uses a composite of diffusion and dislocation creep that can be conveniently formulated, as follows (Karato and Wu, 1993; Karato, 2008):

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_{\mathrm{diff}/\mathrm{disl}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{diff}/\mathrm{disl}}^{-\frac{\mathrm{1}}{n}}{d}^{\frac{m}{n}}{\dot{{\mathit{\epsilon }}_e}}^{\frac{\mathrm{1}-n}{n}}\exp \left({\displaystyle \frac{Q_{\mathrm{diff}/\mathrm{disl}}+{\mathrm{PV}}_{\mathrm{diss}/\mathrm{disl}}}{n\mathrm{RT}}}\right)\text{(7)}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_{\mathrm{comp}}={\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{diff}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{disl}}}}\right)}^{-\mathrm{1}},\end{array}$$

where ${\dot{\mathit{\epsilon }}}_e$ is the square root of second invariant of the deviatoric strain rate. In the case of diffusion creep, n=1 and m>0, while for dislocation creep n>1 and m=0. Definitions and values of other symbols are shown in Table S1 in the Supplement. Surface erosion and sedimentation are neglected.

The plastic yielding is defined by the Drucker–Prager criterion as Eq. (8) (Davis and Selvadurai, 2002):

$$\begin{array}{}\text{(8)}& {\mathit{\eta }}_{\mathrm{yield}}={\displaystyle \frac{C\cos \left(\mathit{\varphi }\right)+P\sin \left(\mathit{\varphi }\right)}{\mathrm{2}\dot{{\mathit{\epsilon }}_e}}},\end{array}$$

where C is the cohesion, and ϕ is the initial friction angle. We also include strain weakening in the plastic regime, through which C and ϕ are linearly weakened by a factor of 2 or 4 between plastic strains of 0.5 and 1.5.

The effective viscosity is eventually given by

$$\begin{array}{}\text{(9)}& {\mathit{\eta }}_{\mathrm{eff}}=\min \left(\max \left({\mathit{\eta }}_i,{\mathit{\eta }}_{\min }\right),{\mathit{\eta }}_{\max }\right),\end{array}$$

where i is one of the subscripts among diff, disl, comp, and yield. η_max and η_min are user-defined viscosity cutoffs (Table S1).

2.3 Initial model configuration and boundary conditions

The numerical model domain is 2000 km × 660 km in the horizontal and vertical directions. It contains a retrocontinental lithosphere and a procontinental lithosphere, with an oceanic lithosphere in between (Fig. 1). A weak zone with a width of 10 km (gabbro; Wilks and Carter, 1990) abuts the oceanic plate to facilitate subduction. The continental lithosphere is 120 km thick, consisting of a 25 km thick upper crust (wet quartzite; Gleason and Tullis, 1995), a 10 km lower crust (wet anorthite; Rybacki et al., 2006), and an 85 km lithospheric mantle. The oceanic lithosphere (∼90 Ma) contains a 4 km thick sediment layer (gabbro; Wilks and Carter, 1990), an 8 km thick oceanic crust layer (gabbro; Wilks and Carter, 1990), and an 88 km thick lithospheric mantle. Dry olivine (Hirth and Kohlstedt, 2003) is used for the mantle. All material properties are listed in Table S1. The numerical resolution incrementally increases along the depth direction from 2 km × 2 km to 8 km × 8 km.

Temperature is fixed at the model surface (0 ^∘C) and the base of the lithosphere (1300 ^∘C), while various continental T_moho is defined ranging from 450–600 ^∘C. The initial continental lithospheric temperature profile follows Chapman (1986), taking into account the thickness of each lithologic layer (Eqs. 11–13). According to these equations, T_moho and H_r are not independent parameters; we thus fixed T_moho when investigating the influence of H_r. Specially, we first choose a surface (T_T), Moho temperature (T_B), and radioactive heat production (H) and then calculate the surface heat flow (q_T), based on Eq. (12). After that, we substitute the above T_T, q_T, and H into Eq. (11) to obtain the continental crustal temperature structure. The oceanic lithospheric temperature distribution follows the half-space cooling model (Turcotte and Schubert, 2002), and in all models, we use the same initial oceanic lithospheric thermal structure. Besides, a temperature gradient of 0.5 ^∘C km⁻¹ is assumed for the sublithospheric mantle (Huangfu et al., 2019).

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{T}_Z=T_{\mathrm{T}}-{\displaystyle \frac{q_{\mathrm{T}}}{kZ}}-{\displaystyle \frac{H{Z}^{\mathrm{2}}}{\mathrm{2}k}}\text{(11)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{B}}=T_{\mathrm{T}}+{\displaystyle \frac{q_{\mathrm{T}}}{k}}\mathrm{\Delta }Z-{\displaystyle \frac{H\mathrm{\Delta }{Z}^{\mathrm{2}}}{\mathrm{2}k}}\text{(12)}& {\displaystyle }& {\displaystyle }{q}_{\mathrm{B}}=q_{\mathrm{T}}-H\mathrm{\Delta }Z\end{array}$$

T_Z is the temperature at depth Z. T is the temperature, q is the heat flux, and the subscripts “T” and “B” denote the top and bottom boundaries of each lithologic layer, respectively. ΔZ is the layer thickness, H is the volumetric heat production, and k is the thermal conductivity.

We apply a free-surface and free-slip lower boundary and impose a 5 cm yr⁻¹ convergence rate to the trailing end of the procontinental lithosphere. Equivalent material flows out through both side walls of the mantle domain to keep mass balance within the whole computational domain (Fig. 1).

3 Results

We conducted 48 numerical experiments by varying the Moho temperature, upper crustal H_{r_uc}, and lower crustal H_{r_lc} of the retrocontinent and procontinent, respectively, to capture different continental thermal structures. Two distinct continental collision evolutionary paths were recognized, namely (i) continuous subduction and (ii) continental subduction with a slab break-off. All the simulations are summarized in Table S2.

3.1 Continuous subduction without slab break-off (Mode I)

Figure 2 shows the development of the continuous subduction mode, and it takes the Model m2 (reference model; Table S2) as an example. The model started with a relatively shallow angle subduction characterized by the oceanic plate dips downward along the low-viscosity weak zone (Fig. 2a1; 3 Myr). During this stage, a portion of oceanic sediment was scraped off and stacked at the surface, accompanied by the slight uplift of the retrocontinental foreland (Fig. 2a2 and c). As the subduction goes on, the dipping angle of the slab in depth gradually steepens due to the increasing slab pull (Fig. 2a2; 9 Myr). At ∼11.5 Myr (Fig. 2a3), continents collide with each other after the full consumption of the oceanic plate. Under the continued oceanic slab pulling, the procontinental lithosphere inherited the subduction, which is characterized by further subducting the angle steepening that resulted in the superposition of retrocontinental upper crust and oceanic sediment on the procontinental forepart. With the proceeding of procontinental subduction, the collision wedge uplifted prominently (Fig. 2c). After that, a large part of the oceanic sediment previously scraped off was then entrained and subducted with the procontinental crust, followed by the rapid uplift at the retro side of the collision zone and the retro-ward advance of the suture (Fig. 2c). From ∼17.5 Myr on, the buoyant procontinental upper crust began to detach from the underlying lithosphere, mixed with the exhumed oceanic sediment, and accumulated together in the subduction channel (Figs. 2a4 and a5). During this time, uplift within the collision zone gradually expands to the retrocontinental interior (Fig. 2c).

Figure 2 Time evolution of Model m2. Snapshots of compositional fields (a1–a5), viscosity (b1–b5), and topography (c) for selected model times are shown. The lithologies and isotherms are the same as in Fig. 1.

3.2 Continental subduction with a slab break-off (Mode II)

In Model m7, both the pro- and retrocontinent have higher T_moho than Model m2 (Table S2). The lithospheric deformation behaviors and the evolution of the model during the oceanic subduction phase are quite similar to those observed in Model m2 (Fig. 3a1). The only difference is that the slab steepened more rapidly at a later stage, which gave rise to the decoupling between the retrocontinental lithosphere and the oceanic plate. Such a process created an ideal space for the asthenospheric upwelling and resulted in the significant weakening of the retrocontinental lithosphere (Fig. 3a2; 8.5 Myr). After ∼10.2 Myr, the combined effect of the increasing slab pull and asthenosphere upwelling promoted necking around the oceanic–procontinental lithospheric transition zone (Fig. 3a3; 11 Myr). The procontinental lithosphere then arrived at the trench and initiated subducting along the inclined weak zone, thus accommodating most of the ongoing convergence. During this stage, the continuous compression and rebound after slab break-off significantly uplifted the collision zone. Later on, the upper part of the retrocontinental lithosphere indented into the procrust, scraping off most of the procontinental upper crust and causing significant thickening and localized shear zones there, while leaving the small residual portion to subduct with the lithosphere. Subsequently, the accumulated orogenic wedge gradually grew into a relatively wide surface uplift.

Figure 3 Results of Model m7. Snapshots of compositional fields (a1–a5), viscosity (b1–b5), and topography (c) for selected model times are shown. The lithologies and isotherms are the same as in Fig. 1.

3.3 Sensitivity tests of continental lithospheric thermal structure on the evolution of continental convergence

A regime diagram (Fig. 4) of all the models summarizes the template simulations investigating the effects of continental T_moho and H_r, and two contrasting end members of continental collision modes are identified. Models with a cold procontinental T_moho (≤450 ^∘C) generally exhibit a continuous subduction mode without slab break-off. Meanwhile, the hotter the retrocontinental T_moho, the steeper the subducting angle is (Figs. 4a and 5a). In comparison, Mode II dominates among the models with a relative hot procontinental T_moho (≥500 ^∘C), especially when the retrocontinental T_moho is greater than 550 ^∘C (Fig. 4a). In addition, H_r is a second-order controlling parameter compared to the T_moho, as it has a greater propensity to alter the upper part of lithospheric deformation styles than the continental collision mode (Fig. 6c–i).

Figure 4 Schematic figures of models with different (a) continental T_moho, (b) upper, and (c) lower crustal radioactive heat production (Table S2). The diagram in each panel shows the dependence of collision mode on the corresponding controlling parameters of both colliding plates.

4 Discussions

4.1 How does continental Moho temperature influence the trajectory of the collision system?

As mentioned above, models incorporating a cold prolithosphere generally evolve into continuous subduction without slab break-off. It may come from the fact that a cold prolithosphere is strong enough to maintain its strength and that it also keeps coupling with the retrolithosphere during plate convergence. Thus, the prolithosphere suffers moderate roll-back, under the condition that the stress is not sufficient to yield the lithosphere and generate break-off (Fig. 2 and S1 in the Supplement). Based on further analysis of this mode of models, we notice that the subducting angle increases as the retrocontinental T_moho increases, accompanied by much more procrust accumulation at shallow and intense retrolithospheric mantle destruction (Fig. 5a–d). The mechanism behind this is that increasing the retrocontinental T_moho can weaken the retrolithosphere and increase the thermal structural difference between the two continents, which may lead to plate decoupling. As a consequence, it is more likely for the prolithosphere to roll back as the magnitude of decoupling increases. Besides, weakening of the retrolithosphere may also offer a favorable condition for the intrusion of accumulated procrustal material into it.

Figure 5 (a) Geometries of subducting plate in models with different retrocontinental T_moho outlined by 800 ^∘C isotherms. Compositional fields of models with decreasing retrocontinental T_moho, (b) m4, retro T_moho=600 ^∘C, (c) m12, retro T_moho=550 ^∘C, (d) m8, and retro T_moho=500 ^∘C. White lines denote the isotherms (200 ^∘C increments). Panel (e) shows the relationship between the times of slab break-off and continental T_moho. Panels (f)–(h) are strain rate of models with various procontinental T_moho.

Figure 6 Temperature profiles (a), topography (b), and strain rate (c–i; at 20 Myr) of models with different crustal H_r.

Regarding Mode II, models with hot procontinental T_moho (≥500 ^∘C) always evolve into continental subduction, following a slab break-off. In these models, when the retrocontinental T_moho is determined (1) as the procontinental T_moho decreases, the procontinental lithosphere becomes cold, which results in a rheologically strong lithosphere that can be subjected to greater deformation; the slab thus breaks off much later (Fig. 5e). This is consistent with many of the previous numerical studies (van de Zedde and Wortel, 2001; van Hunen and Allen, 2011; Duretz et al., 2011). However, it appears that Duretz and Gerya (2013) proposed an apparent opposite tendency. It may come from the fact that these models do not use a layered crust. As a consequence, a strong crust is more likely to result in strong crust–mantle coupling and a deeper slab break-off under slower convergence, while a weak crust leads to crust–mantle decoupling that may evolve into delamination. (2) The depth of the slab break-off in our models is between ∼30–60 km, which is in the range indicated by Davies and von Blanckenburg (1995), Li et al. (2002), and Duretz et al. (2011, 2012). (3) After slab break-off, the buoyant continental lithosphere experiences a rebound accompanied by strong surface uplift within a short duration (Fig. 3c), which coincides with Duretz and Gerya (2013) and Magni et al. (2017). Slab break-off is often considered to be an early process of continental collision, and numerous numerical simulations have been conducted to investigate it. Comparing their results with our models, we suggest that slab break-off is closely related to the rheological strength of the lithosphere. Therefore, the parameters that have prominent impacts on it, such as oceanic plate age, convergence velocity, continental crustal structure, layered crustal rheological strength, and the strength of the interface between the subducting and overriding plates, etc. (Gerya et al., 2004; Burov, 2011; Duretz et al., 2011; Magni et al., 2017; Koptev et al., 2022), may significantly influence the evolutionary path of slab break-off.

In conclusion, we propose that continental T_moho is an important parameter in shaping the continental collision system and especially influences the fate of subduction. The hotter the procontinent, the weaker and easier it is to evolve into a slab break-off. The hotter either the procontinent or the retrocontinent, the less coupled it is with the adjacent plate. In addition, as a result of ongoing subduction, the slab rolls back rapidly, which may generate a gap between the adjacent plates. The asthenosphere then upwells into it, which in turn aggravates the roll-back. This contributes to the development of intense strain localization in the prolithosphere and finally leads to slab break-off (Figs. 3, S2 and 5f–h).

4.2 How crustal radioactive heat production influences the trajectory of the collision system

Radioactive elements are thought to mainly exist in the continental crust, and the radioactive heat production that originates from their decay is one of the important internal heat sources in the Earth, which influences the continental lithospheric thermal structure significantly (Turcotte and Schubert, 2002; Faccenda et al., 2008). According to our model results, increasing the upper crustal H_r can distinctly increase the thermal gradient of the lithospheric upper part, under the effect of which the steep surface topography built at the early stage of continental collision would tend to be relatively gentle. In comparison, varying lower crustal H_r has much less influence on the variation in the continental lithospheric thermal gradient and the evolution of surface topography. Moreover, as the H_r in the crust increases, the lithospheric rheological strength decreases, which facilitates the growth of crustal shear zones. That is, the higher the crustal H_r, the more shear zones it may generate in the crust (Fig. 6c–i). In conclusion, crustal H_r is a crucial parameter in enhancing strain localization in the lithospheric upper part that is closely related to the growth of shear zones, while it exerts a second-order influence on altering the continental collision mode.

4.3 Implications for the development of lateral heterogeneity of Indian–Asian collision

The ∼2000 km long east–west Himalayan–Tibetan orogenic system was created by continental collision from ∼60 to 50 Ma, with the Nanga Parbat syntaxis as its west margin and the Namcha Barwa syntaxis in the east (Chen et al., 2015). It has been proven by contrasting along-strike lithospheric structures (Li et al., 2008; Chen et al., 2015; Kufner et al., 2021). By comparison, our models are capable of providing several first-order fits of the crust–mantle deformation behaviors to this natural collision system, especially at the Hindu Kush profile and the eastern Tibet profile.

Hindu Kush is the westernmost end of the Himalayan–Tibetan orogen. Previous seismic tomography studies have confirmed that there is a north-dipping, downward-steepening, and thinning high-velocity anomaly (HVA), accompanied by increasingly intense seismicity, underneath the region (Kufner et al., 2016, 2021). The inclined HVA is believed to be a subducting plate of Indian provenance and is part of the Indian–Asian collision system. The thinning of the HVA is in the extent between 180–220 km and was suggested to be at the transition zone between Indian continental lithosphere and Neo-Tethys oceanic lithosphere. In agreement with the observations, Model m2 has a similar tectonic background and shows an analogous continuous subduction characterized by a nearly vertical dipping angle and a slab necking at a depth of ∼220 km. In addition, the relative narrow collision orogen in the model is also consistent with the natural topography in this region (Figs. 7a and 2c). It is worth noticing that the observed inclined distribution of earthquakes seems to coincide with the subducting continental crust. This may offer a possible explanation for the active intermediate earthquakes beneath this region.

Figure 7 (a) Geographic setting of Himalayan–Tibetan orogen. Schematic geologic cross sections across the (b) western (Hindu Kush; Kufner et al., 2021) and eastern Himalayan–Tibetan orogen (Li et al., 2008; Wang et al., 2019). Panel (c) shows the compositional field of Model m2 (mode I). Panels (e) and (g) are the strain rate and compositional field of Model m7 (Mode II), respectively.

The main part of the Himalayan–Tibetan orogen is the product of a collision between the cold Indian craton and the relative warm Asian continent after the closure of Neo-Tethys, the geological settings of which are quite analogous with our Model m7. During the continental collision phase, the Indus–Yarlung suture zone bordering the Indian and Asian lithospheres moved significantly, creating an ∼8 km mountain (the Himalayas) and a ∼4 km plateau (the Tibetan Plateau; Yin and Harrison, 2000). Unlike continuous subduction underlying the Hindu Kush, geological and seismological evidence shows that the eastern part of the profile is underlain by a shallower-angle subducting Indian lithosphere, following a slab break-off; these are also similar to the model results of Model m7 (Figs. 7b and 2c). Furthermore, previous field observations have shown that the northern part of the Indian continental crust is abundant with radioactive elements (Vidal et al., 1982; Scaillet et al., 1990; Macfarlane, 1992; Faccenda et al., 2008). Our sensitivity tests have recognized that continental crust with high radioactive heat production is more prone to result in brittle fractures and intense deformation of the lithospheric upper part, which may lead to the development of new shear zones (Fig. 6). This resembles the tectonic characteristics of large thrust fault zones in the Himalayas (Fig. 7b).

As a consequence, we speculate that the continental lithospheric thermal structural difference may be one of the key parameters that control the evolution of lateral heterogeneity along the Himalayan–Tibetan orogen and that the high content of crustal radioactive elements is one of the significant factors that dominates the growth of large thrust fault zones in the Himalayas.

5 Conclusions

In this work, we systematically discuss the influences of the continental lithospheric thermal structure on the evolution of the continental collision system based on high-resolution thermomechanical numerical experiments. The model results demonstrate that two end-members of continental collision are obtained because the continuous subduction mainly occurs with a relative cold overriding lithosphere (T_moho≤450^∘), and as the retrocontinental T_moho increases, the subducting angle steepens. Slab break-off dominates when the model has a relatively hot procontinental T_moho (≥500 ^∘C), especially when the retrocontinental T_moho is greater than 550 ^∘C. H_r is a second-order controlling parameter compared with T_moho in shaping the continental collision mode, while it is more prone to facilitating the growth of crustal shear zones.

The continental lithospheric thermal structure may have played a significant role in the development of lateral heterogeneity along the Himalayan–Tibetan orogenic belt. We suggest that the different evolutionary paths between the continuous subduction underlying the Hindu Kush and the continental subduction beneath eastern Tibet may come from the inherited lateral variation in the Indian lithospheric geothermal gradient. In addition, the high content of radioactive elements in the continental crust may be one of the important reasons for the development of deep and large thrust fault zones in this collision orogen.