Under-thrusting of continental lithospheric mantle controls the architecture of crustal deformation: numerical modelling of the northern margin of the Tibetan Plateau Open Access

The northern margin of the Tibetan Plateau is defined by the active Qilian Shan fold–thrust belt (e.g. Yin et al. 2008; Yuan et al. 2013). This Cenozoic tectonic deformation belt is c. 350 km wide and is bound by the Hexi Corridor foreland basin and the Alxa block to the north and the Qaidam Basin to the south. The Alxa block was previously considered to be the western extension of the North China Craton, but has recently been shown to be a complex of Paleozoic continental and oceanic magmatic arcs comprising the southern part of the Central Asian Orogenic Belt (e.g. Hui et al. 2021). The Qilian Shan Cenozoic contractional structures overprint a complex pre-Cenozoic geological history in the Qilian Shan that dominantly involved an early Paleozoic orogeny (e.g. Yin et al. 2008; Li et al. 2021; Wang et al. 2022). It is generally believed that, during the early Paleozoic, northward-dipping subduction systems of the Prototethyan oceanic branches were developed in the Qilian regime, specifically along the southern margins of the Alxa, North Qilian and South Qilian blocks, which resulted in these blocks being characterized by intense arc magmatism during the early Paleozoic (Zhang 2002; Song et al. 2013).

During the Silurian, these Prototethyan oceanic branches in the Qilian regime closed and the Qaidam–Qilian continental blocks were amalgamated with the southern Asian margin (Zhang 2002; Song et al. 2013). Notably, following the final collision between the Qaidam and Qilian continental blocks during the Silurian, a large-scale ultrahigh-pressure orogenic belt formed along the northern margin of the Qaidam block due to the deep-seated subduction of this continental block underlain by oceanic plateau basement (Zhang et al. 2024) (e.g. Zhang 2002; Song et al. 2013; Li et al. 2021; Zhang et al. 2024). During the late Paleozoic to Mid-Triassic, the Qilian regime was subject to magmatism and rift sedimentation related to the collapse of the early Paleozoic orogens and/or back-arc extension due to the subduction of the Palaeotethyan Ocean along the southern Kunlun margin (e.g. Zeng et al. 2018; Tian et al. 2022).

The Qilian regime is therefore characterized by the impacts of long-lived, multi-stage and intense arc- and/or rifting-related magmatism and sedimentation during the Paleozoic to Mid-Triassic and is thus underlain by a strongly reworked lithosphere characterized by a weak lithospheric strength and several discrete pre-existing lithospheric-scale weak zones (Song et al. 2013; Zuza et al. 2018). Consequently, this regime was readily reactivated during the Cenozoic convergence caused by the India–Asia collision to the south (e.g. Yin et al. 2008; Allen et al. 2017; Zuza et al. 2018, 2019). This rejuvenation led to SW under-thrusting of the Asian lithospheric mantle (ALM) beneath the Qilian Shan lithosphere, as evidenced in the Cenozoic crustal shortening record (Zuza et al. 2016, 2019) and geophysical observations (Ye et al. 2015, 2021; Gao et al. 2022).

The Qilian Shan fold–thrust belt has a broad bivergent kinematic structure, with NE-directed thrusts against the Hexi Corridor and SW-directed faults against the Qaidam Basin. The thermochronology and basin records suggest that the Qilian Shan started deforming and uplifting shortly after India collided with Asia (Yin et al. 2008; Bovet et al. 2009; Li et al. 2019, 2020; Wang et al. 2021). This rapid onset of deformation supports the view that the pre-existing early Paleozoic Qilian suture was readily reactivated in the Cenozoic. Explicit kinematic details of early Cenozoic deformation in the Qilian Shan are poorly resolved due to renewed and accelerated shortening in the Miocene. However, the low-temperature thermochronology reveals an outward-younging pattern of fault activity from the central Qilian Shan since the Miocene (Yuan et al. 2013; W.J. Zheng et al. 2013; D. Zheng et al. 2017; Zhang et al. 2022). The spatial and temporal pattern of the fault system is consistent with its recent kinematic development as a bivergent wedge (Sun et al. 2022). The culmination of Cenozoic crustal shortening has thickened the Qilian Shan crust to 60–70 km (e.g. Yue et al. 2012; Allen et al. 2017), with surface elevations >3.5 km.

The active deformation in the Quaternary indicated by the geomorphic record and seismicity suggests that the external flanks of the fold–thrust belt are the most active (W.J. Zheng et al. 2013; D. Zheng et al. 2017). The central-southern parts of the tectonic belt may be experiencing crustal-scale buckling based on a comparison of the topography and thermochronology datasets (Jolivet et al. 2022), potentially due to the development of a mid- to lower crust duplex at depth (Zuza et al. 2019).

In summary, the unique structural characteristics of the Qilian Shan fold–thrust belt include: (1) a bivergent orogenic wedge with a relatively aseismic and tectonically quiescent central segment; (2) a relatively large crustal thickness (c. 60+ km) and elevations (c. 3.5 km); and (3) focused shortening along the marginal thrust systems. We hypothesize that these characteristics are uniquely controlled by the under-thrusting of the ALM and we set out to test this via numerical simulations, which provide unique constraints on the upper crustal strain patterns linked with deep crust and mantle processes.

We structured our model to be broadly similar to observations from the Qilian Shan fold–thrust belt. The finite element model had a width of c. 360 km based on the north–south dimensions of the Qilian Shan, from the Hexi Corridor in the north to the northern margin of the Qaidam Basin in the south (Fig. 1). The model is 120 km deep and is divided into the upper crust, the lower crust and the Tibetan lithospheric mantle of an ALM. The upper crust starts with a thickness of 15 km, 240 km away from the southern end of the model. The upper surface is inclined to the north at 1.5°. The bottom surface of the upper crust is inclined to the south at an angle of 3° based on existing geophysical constraints (Fig. 1b; Gao et al. 2013; Chen et al. 2022). The Moho depth is 65 km in the south, which gradually decreases to 50 km in the centre, from where it remains constant to the northern end (Fig. 1a, b). The material parameters of the under-thrusting model are listed in Table 1.

The model is first loaded by gravity and then the two lateral boundaries converge towards the middle of the model. The first control model involved pure shear shortening of a simple rectangular model geometry without any under-thrusting geometry. The second control model was subjected to pure shear tectonic shortening with an under-thrusting architecture, but there was no décollement between the ALM and the crust. The third model was our under-thrusting model, in which ALM has a relative displacement movement with the crust along the décollement (the pink line in our model) (Fig. 2), which represents the southbound subduction channel following the reactivation of the ALM.

The primary models involved c. 15 km of convergence for each lateral boundary for c. 30 km of total shortening across the model. The simulation time started from the far-field stress transfer from the collision of the India–Asia plate to the Qilian Shan at the NE edge of the Tibetan Plateau. The ALM under-thrust to the left through the décollement layer (Fig. 2, pink) and is relatively obstructed by the Tibetan lithospheric mantle on the left-hand side. The under-thrust displacement is reduced to 12 km. The total modelled shortening is c. 30 km, which is substantially lower than Cenozoic shortening estimates across the Qilian Shan of >100 km (e.g. Allen et al. 2017; Zuza et al. 2019). We argue that our models of the deformational state are diagnostic of the predicted strain patterns and the kinematic architecture. The deformation of the final model can be compared with recent to present-day observations across the Qilian Shan of the seismicity, topography, Quaternary landforms and strain distribution. Furthermore, seismic reflection profiles across the Qilian Shan reveal that it mostly involves trans-crustal, basement-involved thrust faults rather than ramp-flat, thin-skinned fault propagation (Gao et al. 2022), which broadly matches models of plastic deformation.

In addition to the reference model in Figure 2, we developed two models to study the effects of the initial geometry and boundary conditions. Figure 3a shows a pure shear model without the southward-subducting ALM, featuring a horizontal, frictionless interface between the crust and lithospheric mantle. Figure 3b shares the same geometry as the reference model, but differs by having a coupled, non-frictional contact between the crust and lithospheric mantle. To observe the plastic strain that occurs during the loading of the boundary conditions, we quantified the EPS, which describes the permanent deformation when a material is under a stress state beyond the elastic limit. The results of the control model involving pure shear tectonic shortening with or without an under-thrusting architecture are shown in Figure 3.

Before the tectonic loading was applied to the model, the initial stress state of the model was in equilibrium with the gravity-induced stress generated by the applied boundary conditions. In the earliest stage of tectonic shortening, the model body did not experience plastic strain. After continuous tectonic shortening, plastic strain concentrated as a series of high-strain zones that look like convergent thrust ramps. The local plastic strain areas represent faulting and fracture (e.g. Strayer et al. 2001; Panian and Wiltschko 2004, 2007). We therefore express the spatial location and time series of the development of the EPS concentration zone in the model results as the evolution of faults.

The two control models (Fig. 3) show similar EPS distributions with discrete conjugate thrust faults in the middle of the models and more distributed strain between these structures. A discrete north-dipping fault is observed at the southern edge of the model and more distributed strain is observed at the north end of the model. The plastic strain bands in the models dip c. 39–42°, as expected for models of perfectly plastic deformation. This control model illustrates that the fault dip angle is controlled by the given friction angle, which is slightly higher than the theoretical value (c. 37–39°) (see Fig. S2 for details). The two pure shear models are subtly different, with the EPS in the simple pure shear model showing more uniform, distributed shortening, whereas the under-thrusting model shows reduced strain south of the subduction inflection point (Fig. 3). This strain reduction hints at the influence of the subducting plate.

Figure 4 shows our primary lithospheric mantle subduction model. Our modelling of elastic–plastic materials does not adequately capture large strains or complex fold–thrust belt evolution over time, but resolves the deformation and change in state from the initial state to the final state. After 12 km of tectonic shortening (3.33% shortening strain) (6 km at each end), minor deformation appears in the middle of the model above the subduction inflection point (Fig. 4a). With continued shortening to 16 km (4.44% stain), a robust EPS concentration zone appears in the middle of the model with discrete, connected fault structures (Fig. 4b). Minor strain appears in upper crust at the southern end of the model.

After 20 km of total shortening (5.56% strain), three discrete conjugate shear zones with a high EPS develop in the middle of the model (Fig. 4c). These faults merge near the subduction inflection point (Fig. 4c). Continued shortening (24 km shortening or 6.67% strain) amplifies the previously established EPS and a major discrete north-dipping EPS shear zone links from the southern edge of the model to the lower point of the under-thrust slab (Fig. 4d). There is a pronounced low-strain area between these north-dipping structures and the main central region of conjugate faulting. The plastic strain value of this area is almost the lowest in the model (Fig. 4d). In the final stage of our modelling with 29 km shortening (8.06% strain), the main EPS geometries remain unchanged and the low-strain regions persist (Fig. 4e). Several distributed EPS regions develop in the shallow upper crust (Fig. 4e).

Bivergent orogenic wedges have been found to operate at several spatial scales, ranging from accretionary prisms (Byrne and Hibbard 1987), such as the Aleutian Arc (Fliedner and Klemperer 2000) and the Sunda Arc (Silver and Reed 1988), to the scale of fold–thrust belts, such as the Sierras Pampeanas (Cristallini et al. 2004), and to the scale of entire orogens, such as Pyrenees (Choukroune 1989). A bivergent wedge is traditionally divided into two asymmetrical opposing thrust wedges that comply with critical taper theory (e.g. Davis et al. 1983; Dahlen 1984). These wedges can be further divided into five components: the pro-foreland basin, the pro-wedge, the retro-wedge, the retro-foreland basin and the central uplift (Willett et al. 1993; Hoth et al. 2006; Fig. S8).

In an idealized model, a velocity discontinuity (singularity) separates the down-going plate from the overriding plate, and the subduction–accretion process forms a bivergent wedge (Willett et al. 1993). The asymmetry of the convergent geometry then leads to the polarity of the crustal mass transfer. Through time, a fully developed bivergent orogen grows and propagates towards its respective forelands with a symmetrical structure (Beaumont et al. 1996).

Our simulated lithospheric mantle subduction system resembles a classical bivergent wedge, which develops to both sides from a subduction inflection point. However, a key difference is that the doubly vergent orogen formed in our models is extremely asymmetrical. Furthermore, the crustal material that comprises the bivergent wedge emanates from the upper crust and is not accreted from the down-going plate. Our model does not uniquely imply a single S-point singularity for deformation to emanate in both directions. Instead, the model shows a more complex zone in the lower crust, with a foreland-propagating thrust fault emanating from the first concave inflection of the south-directed Asian lithosphere and a second hindward/retro-wedge thrust emanating from the second convex inflection. The dimensions of this wedge are therefore initially set by these two parameters and the crustal wedge grows upwards and laterally outwards with time. Our results suggest that shallower slabs will have wider orogenic wedges and steeper subducting slabs will have narrower orogenic wedges because the distance between these inflections points is less (Fig. 5). This simple prediction is probably enhanced by the fact that shallower subduction imparts greater basal traction to the upper plate/crust to grow a wider orogenic wedge, which supports our assertion.

We recognize that deformation models are influenced by a variety of factors (Beaumont et al. 1996; Beaudoin et al. 2023). We systematically investigated the effects of the initial geometry, elastic–plastic parameters and bottom-slip layer friction coefficients on fold–thrust fault zones (Chen et al. 2022). In addition, we examined the effects of different boundary conditions (extrusion, shear, blocking and rotation) on the spatial and temporal evolution of arcuate tectonic zones (Chen et al. 2023). In this paper, due to space constraints, we have chosen to focus on the effect of changes in the dip angle of the subduction interface on deformation modelling.

The oldest structures are found in the centre of the wedge in the Qilian Shan and the structures young outwards to the north and south. This is consistent with the spatial and temporal pattern of the fault system since the Miocene as a broad bivergent wedge (Sun et al. 2022), as revealed by low-temperature thermochronology (Yuan et al. 2013; H. Zhang et al. 2017; B. Zhang et al. 2022). The orogen consists of crustal-scale, basement-cored thrust faults (Gao et al. 2022) and does not exhibit thin-skinned, ramp-flat, fold–thrust belt geometries. The culmination of Cenozoic crustal shortening has thickened the Qilian Shan crust to 60–70 km (e.g. Yue et al. 2012; Zuza et al. 2016; Allen et al. 2017) with surface elevations >3.5 km. The Quaternary active deformation observed in the geomorphic record and the seismicity both suggest that the external flanks of the fold–thrust belt are the most active (W.J. Zheng et al. 2013; D. Zheng et al. 2017). The central-southern parts of this tectonic belt may be experiencing crustal-scale buckling based on a comparison of the topography and thermochronology datasets (Jolivet et al. 2022). To better explore the deformation and strain features in the Qilian Shan, we compared these observations with our modelling results. In this section, we discuss the seismicity, fault age and evolution trends, and geomorphic proxies.

Major thrust faults across the Qilian Shan cut through the entire crust (Gao et al. 2013, 2022, Fig. S5) to develop a wedge-like geometry (Ye et al. 2021; Fig. S10). This overall geometry matches our EPS strain patterns (Fig. 4e) and therefore confirms the appropriateness of our modelling approach. Existing records of the Qilian Shan thrust propagation show a broad outward-younging pattern of Miocene deformation from c. 20 to <5 Ma along the frontal thrust margins (e.g. W.J. Zheng et al. 2013; D. Zheng et al. 2017; Zhang et al. 2022; Fig. 6a). Our models predict a similar outward growth from the lithospheric mantle subduction interface, as observed through differentiation of the surface EPS through time (i.e. the loading step) (Fig. 6b). The middle of the model above the subduction interface first experiences an increase in EPS at loading step 20 and deformation progressively migrates outwards toward the orogen margins. We believe that the model accurately produces deformation of the bivergent orogenic wedge pattern, which is consistent with observations from the Qilian Shan.

Geomorphic records and Quaternary deposits serve as proxies for upper crustal strain on shorter timescales. We compared the Quaternary geology synthesis of Zhang et al. (2017) with our model results (Fig. 7). The southern Qilian Shan has relatively flat erosional surfaces that suggest erosional planation and a relative lack of shortening strain since at least the Miocene (Zhang et al. 2017). Conversely, the central-northern parts of the Qilian Shan display the most topographic relief and have widely distributed Quaternary basin deposits. These observations suggest that the central-northern Qilian Shan has experienced relatively active shortening, range exhumation and thrust-induced loading to allow intermontane basin sedimentation. Comparison of these geomorphic records with our numerical simulations reveals similarities, with limited deformation in the central-southern Qilian Shan and focused high-frequency thrust faulting in the central-northern Qilian Shan (Fig. 7c).

Jolivet et al. (2022) demonstrated how the landforms and thermochronology records of the central Qilian Shan are consistent with orogen-scale buckling between the marginal fold–thrust belts. Our modelling, the observed geomorphology and the present-day seismicity are consistent with this broad kinematic model. However, our results highlight how the kinematics of mid- to deep crustal deformation influences the patterns of upper crustal strain. In particular, the geometry and location of the southward-directed under-thrusting ramp appear to control the locations of focused strain and regions of local tectonic quiescence.

A compilation of historical seismicity shows regions of enhanced deformation and relative tectonic quiescence. Our earthquake compilation includes historical earthquakes of M > 5 (from ad 870 to 1970, although most of them lack depth information) and of M > 3 (from 1970 to 2021) (Fig. 7a). This compilation reveals earthquakes distributed along the northern and southern margins of the Qilian Shan and a prominent seismic gap in the south-central Qilian Shan (Fig. 7). This pattern is similar to our modelling results, where the EPS is concentrated in the northern and southern flanks of the growing orogen, with a prominent lack of strain in the central-southern parts of the model (Fig. 7c).

The pure shear simulations with or without under-thrusting do not accurately predict the structural framework of the Qilian Shan (Fig. 3). The models without under-thrusting have more uniformly distributed shortening strain patterns. This pattern is at odds with the overall bivergent structure of the Qilian Shan. These models are unable to produce deep crustal-scale ramps and the strain patterns, including areas of focused and negligible deformation, do not match the observations from the Qilian Shan.

In summary, we have demonstrated that pure shear deformation yields deformation patterns that are inconsistent with the observations from the Qilian Shan, including crustal seismicity, fault kinematics and landscape evolution. Instead, a model of southward under-thrusting of the mantle lithosphere better reproduces the strain characteristics of the Qilian Shan, which supports previous geophysical and geological interpretations of the northern Tibetan Plateau.

Our models focus on the elastic–plastic behaviour of the Earth's lithosphere and the frictional contacts within the décollement layers, effectively simulating tectonic deformation at the upper crustal level and capturing the spatial and temporal evolution of the fault system in the Qilian orogenic belt. However, the models do not account for the viscosity of the lower crust and upper mantle, which has a role in stress relaxation at depth. Future studies should incorporate the effects of viscous flow within the Earth's lithosphere. Although the southward thrusting of the Asian lithosphere is well documented by geophysical data, the distribution of north–south shortening across the boundaries of the Qilian orogenic belt requires further quantitative analysis. The boundary conditions for north–south shortening presented in this study must be refined with additional observational data.

To assess how the observed and modelled strain patterns might result in surface uplift of the Qilian Shan, we compared our results with the topography of the range. The final elevation of the model is the free uplift displacement of the upper surface plus the initial elevation of the model. We did not consider the additional effect of erosion or sedimentation in the calculation, but only the thickening and uplift of the surface under the effect of bidirectional shortening of the north–south boundary. We assumed an initial elevation of c. 3 km across the Qilian Shan, with lower elevations in the Hexi Corridor to the north (e.g. Chen et al. 2022) (Fig. 8, initial terrain 1).

We first examined how pure shear shortening could impact the topography. Pure shear shortening results in broad uplift that is most concentrated in the southernmost part of the model and the topography mostly decreases toward the north (Fig. 8, model result 1). Next, we considered the under-thrusting model using the same initial topography as the pure shear shortening model (Fig. 8, initial terrain 1). These results show a final topography that is highest in the middle part of the mode (i.e. above the subduction interface) and lower to the north and south (Fig. 8, model result 2).

Both results led to surface elevations that were higher than the modern-day topography. However, these models do not consider erosion and we would expect the highest elevation and relief surfaces to be eroded, thus reducing the final elevations of these thickened locations. The thermochronology and geological records are consistent with <5 km of upper crustal erosion (e.g. Zuza et al. 2016; Jolivet et al. 2022), which could reasonably reduce our modelled elevations to the present-day observed elevations of the Qilian Shan. The pure shear model resulted in surface uplift patterns that were more dissimilar from the modern topography than the under-thrusting models (Fig. 8) – namely, the pure shear result generated the highest topography in the south, which gradually tapered to the north, whereas the under-thrusting model resulted in a bivergent topographic profile with a higher topography in the central-north regions. The under-thrusting model therefore resulted in a modelled surface uplift that may better match the modern topography, especially from its spatial derivative, which shows a greater topographic slope in the northern part of the Qilian Shan and a smoother topography to the south, close to the derivative of the real elevation (Fig. 8, derivative of model result 2, derivative of real elevation).

There are two possibilities that may explain the greater crustal thickness in the Qilian Shan. First, this region may have had a pre-existing topography due to poorly resolved Cretaceous shortening events (e.g. Li et al. 2023; Wang et al. 2022). This configuration with a pre-existing thickness is speculative, but consistent with recent geological investigations suggesting that collision of the Lhasa block with Eurasia led to deformation in the Qilian Shan (Vincent and Allen 1999; Jolivet et al. 2001).

Alternatively, an earlier, pre-subduction phase of distributed crustal shortening in the North Qaidam and southern Qilian Shan thrust belts may have thickened these regions (Sun et al. 2022). As discussed earlier, the phase of Paleogene deformation in the Qilian Shan is supported by field studies and thermochronology, but the extent of this deformation is not well constrained. It is possible that earlier deformation thickened the central and southern Qilian Shan, similar to the pre-existing topography condition described here. The model result 2 in Figure 8 recreates the modern topography, thus demonstrating the viability of this scenario. More detailed information about the terrain can be found in Figures S3 and S4.

To investigate the influence of different modes of intra-plate deformation, we compared our results and analysis of the Qilian Shan with the nearby Cenozoic Tian Shan (Fig. 9). The Tian Shan is an intra-continental orogen that formed in response to the Cenozoic India–Asia collision, located c. 1500 km north of the Indus–Tsangpo suture (Tapponnier and Molnar 1979; Burchfiel and Royden 1991; Yin et al. 1998). The Tian Shan experienced a protracted subduction and collisional history in the early Paleozoic (e.g. Windley et al. 1990; Şengör et al. 1993). After continuous plantation with the formation of a series of Jurassic intermontane basins, the Tian Shan orogen was later reactivated in the Cenozoic as a fold–thrust belt (Yin et al. 1998). Thermochronological and stratigraphic studies indicate that crustal shortening and range growth within the Tian Shan started in the Oligocene, with accelerated deformation in the Miocene (Yin et al. 1998; Sun et al. 2009; Abdulhameed et al. 2020). Although the style of deformation across the Tian Shan is debated (e.g. Zhao et al. 2003; Zhang et al. 2020; Li et al. 2022), updated constraints suggest that the Tian Shan deformed as a pure shear orogen (Li et al. 2022; Fig. S12) based on geophysical and geodetic observations, seismicity, the GPS velocity, and the widely distributed shortening and Quaternary faulting distributed across the Tian Shan with no significant increase near boundaries (Abdrakhmatov et al. 1996; Thompson et al. 2002; Zubovich et al. 2010; Charreau et al. 2017).

The Phanerozoic geological history of the Tian Shan is broadly similar to that of the Qilian Shan, including Paleozoic oceanic subduction and arc magmatism, and arc–terrane collisions (Windley et al. 1990; Şengör et al. 1993), which allows us to test whether structural differences between the two orogens reflect fundamentally different mechanisms of formation. As discussed here, the Cenozoic Qilian Shan bivergent wedge formed due to the southward subduction of the ALM. This resulted in a wide bivergent thrust wedge, focused strain along the wedge margins and a zone of tectonic quiescence within the wedge just behind the inferred subduction interface. The relatively uniform strain patterns across the Tian Shan (e.g. Zubovich et al. 2010; Charreau et al. 2017) are inconsistent with our lithospheric mantle subduction models, but are more similar to our pure shear simulations (Fig. 3). Our results therefore further support pure shear models for the Cenozoic Tian Shan (Li et al. 2022) and are evidence against under-thrusting models (e.g. Lü et al. 2019; Zhang et al. 2020).

The distinct styles of crustal deformation between the Tian Shan and Qilian Shan can be explained by major differences in how the orogens were constructed, such as pure shear orogeny v. lithospheric mantle subduction. The results from our simulations provide testable predictions for intra-plate orogens to determine the relative contributions of either of these orogenic modes. These insights can therefore be used to better interpret other bivergent orogenic wedges, both past and present.

We generated a lithospheric-scale finite element model of the Qilian Shan bivergent fold–thrust belt along the northern margin of the Tibetan Plateau. Based on geological and geophysical interpretations of southward under-thrusting of the ALM beneath the Tibetan Plateau, we tested the influence of this lithospheric mantle subduction geometry on records of crustal strain. Initial simulations of homogenous pure shear shortening did not recreate the observed faulting records of the Qilian Shan. Conversely, models that incorporated southward under-thrusting accurately predicted regions of high and low strain that matched records for the seismicity, surface geology, geomorphology and thermochronology.

Our results emphasize the roles that under-thrusting and subduction of the lithospheric mantle have on the generation of a bivergent fold–thrust belt, the topography and landscape evolution. The modelling results support the end-member mode of southward under-thrusting of the ALM. Specifically, we suggest that crustal strains are directly controlled by the geometry of the mantle subduction interface. These insights can be used to better interpret other bivergent wedges, both past and present.

We would like to thank the reviewers for their constructive comments.

QC: conceptualization (equal), methodology (equal), visualization (lead), writing – original draft (lead), writing – review and editing (equal); CH: funding acquisition (lead), methodology (lead), project administration (lead), supervision (lead), writing – original draft (equal), writing – review and editing (equal); AVZ: conceptualization (lead), formal analysis (lead), writing – original draft (lead), writing – review and editing (lead); K-JZ: visualization (equal), writing – review and editing (equal); HZ: funding acquisition (equal), writing – review and editing (equal); YS: project administration (equal), supervision (equal), writing – review and editing (equal).

This research was financially supported by the National Science Foundation of China (Grant No. 42074117) and the Fundamental Research Funds for the Central Universities. A. Zuza acknowledges the support from the National Science Foundation (Grant No. EAR 1914501).This work was funded by the Key Technologies Research and Development Program (2023YFF0804300, 2023YFF0804301).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

All data generated or analysed during this study are included in this published article.

Some maps were generated using the Generic Mapping Tools (Wessel et al. 2019) available at https://www.generic-mapping-tools.org/. The numerical modelling was performed through the General-Purpose Finite Element analysis software ABAQUS Version 2016 (ABAQUS 2016). The historical earthquake data include M_S > 5 (Department of Earthquake Disaster Prevention, China Earthquake Administration, 1995, 1999, 879 BC to1970 AD) and M > 3 (1970–2022, Centre of the China Earthquake Networks, http://data.earthquake.cn/) earthquakes. The ages of the main Qilian Shan faults are taken from W.J. Zheng et al. (2013), D. Zheng et al. (2017) and Zhang et al. (2022). The Quaternary basin fill and the erosional surfaces of the Qilian Shan are from Zhang et al. (2017).