Contrasting exhumation histories and relief development within the Three Rivers Region (south-east Tibet)

The Three Rivers Region in south-east Tibet represents a transition between the strongly deformed zone around the Eastern Himalayan Syntaxis (EHS) and the less deformed south-east Tibetan Plateau margin in Yunnan and Sichuan. In this study, we compile and model published thermochronometric ages for two massifs facing each other across the Mekong River in the core of the Three Rivers Region (TRR), using the thermo-kinematic code Pecube to constrain their exhumation and relief history. Modelling results for the low-relief (< 600 m), moderate-elevation (∼ 4500 m) Baima Xueshan massif, east of the Mekong River, suggest regional rock uplift at a rate of 0.25 km/Myr since ∼ 10 Ma, following slow exhumation at a rate of 0.01 km/Myr since at least 22 Ma. Estimated Mekong River incision accounts for 30 % of the total exhumation since 10 Ma. We interpret exhumation of the massif as a response to regional uplift around the EHS and conclude that the low relief of the massif was acquired at high elevation (> 4500 m), probably in part due to glacial “buzzsaw-like” processes active at such high elevation and particularly efficient during Quaternary glaciations. Exhumation of the Baima Xueshan is significantly higher (2.5 km since ∼ 10 Ma) than that estimated for the most emblematic low-relief “relict” surfaces of eastern Tibet, where apatite (U–Th) $/$ He (AHe) ages > 50 Ma imply only a few hundreds of metres of exhumation since the onset of the India–Asia collision. The low-relief Baima Xueshan massif, with its younger AHe ages (< 50 Ma) that record significant rock uplift and exhumation, thus cannot be classified as a relict surface. Modelling results for the high-relief, high-elevation Kawagebo massif, to the west of the Mekong, imply a similar contribution of Mekong River incision (25 %) to exhumation but much stronger local rock uplift at a rate of 0.45 km/Myr since at least 10 Ma, accelerating to 1.86 km/Myr since 1.6 Ma. We show that the thermochronometric ages are best reproduced by a model of rock uplift on a kinked westward-dipping thrust striking roughly parallel to the Mekong River, with a steep shallow segment flattening out at depth. Thus, the strong differences in elevation and relief of two massifs are linked to variable exhumation histories due to strongly differing tectonic imprint.

1 Introduction

Despite its high mean elevation of ∼ 5000 m, Tibet is mainly characterized by low relief (< 1 km), with an average slope of ∼ 5^∘ (Fielding et al., 1994). However, topographic characteristics vary throughout the Tibetan Plateau. Central and north Tibet are characterized by extensive low-relief surfaces, with narrow, linear, and parallel mountain ranges that rise slightly above the mean elevation of the plateau, defining “positive topography” where elevation is positively correlated with relief and mean slope (Liu-Zeng et al., 2008). In those regions, structural relief that was tectonically generated during the collision between India and Asia has been smoothed out locally by internal drainages, which transport sediments from the nearby ranges to fill local intermontane basins but cannot evacuate sediments out of Tibet (Meyer et al., 1998). In contrast, south-east Tibet is characterized by “negative topography”, where deeply incised valleys are separated by patches of low-relief surfaces at a mean elevation of ∼ 4500 m, slightly lower than the average plateau height (Fig. 1a). Although only scattered remnants of these surfaces occur in the Three Rivers Region of south-east Tibet, they gradually become more continuous northward and eventually integrate into the extensive low-relief surfaces of the plateau interior (Fig. 1b, c). The most emblematic surface with remarkable low relief occurs on the Triassic Daocheng granite (Fig. 1c, d), located between the Yangtze and the Yalong rivers draining the south-eastern plateau margin in Yunnan and Sichuan (Fig. 1c). These low-relief surfaces, with < 600 m relief in a radius of 5 km, have previously been interpreted as “relict surfaces”, i.e. remnants of a paleo-landscape that was formed at low elevation and subsequently uplifted (Clark et al., 2005, 2006). No such low-relief relict surfaces are observed further west, between the Salween River and the Eastern Himalayan Syntaxis (Fig. 1b, c).

Samples collected from these low-relief surfaces in south-east Tibet show much older thermochronometric ages than samples from the deep river valleys that dissect them (Fig. 1b). The Daocheng plateau and other low-relief surfaces of the Yunnan and Sichuan margin are characterized by apatite (U–Th) $/$ He (AHe) ages > 50 Ma and apatite fission track (AFT) ages > 100 Ma, whereas samples collected from the gorges of the Dadu, Yalong, and Yangtze rivers show AHe ages < 15 Ma and AFT ages < 50 Ma (e.g. Clark et al., 2005; Ouimet et al., 2010). These data have been interpreted as recording a regional rapid incision phase starting between ∼ 13 and 9 Ma in eastern Tibet, which was considered a proxy for widespread plateau uplift shortly preceding this incision (Clark et al., 2005; Ouimet et al., 2010). However, the direct link between river incision and regional tectonic uplift has been challenged by others, who argued for a delay between crustal shortening, occurring mostly in the Eocene, and significant incision mostly during the Miocene in south-east Tibet (Liu-Zeng et al., 2008). A purely climatic forcing on river incision by intensified monsoonal precipitation at ∼ 17 Ma was proposed by Nie et al. (2018). Other thermochronological studies have also provided evidence for earlier phases of rapid exhumation, the timing of which varies regionally between 30–20 Ma in the Longmenshan (Wang et al., 2012; Tan et al., 2014) to 40–30 Ma in the Yalong thrust belt (Zhang et al., 2016) and ∼ 60 to ∼ 40 Ma in the Baima Xueshan massif (Liu-Zeng et al., 2018). The latter has been linked to uplift of the south-east Tibetan Plateau, consistent with paleo-elevation data implying that the plateau has been close to its present-day elevation since the late Eocene–Oligocene (Hoke et al., 2014; Li et al., 2015; Wu et al., 2018).

An inherent problem with thermochronology data is that they do not provide direct information on (either rock or surface) uplift (England and Molnar, 1990; Reiners, 2007). However, such data do potentially provide constraints on paleo-relief (Braun, 2002; Reiners, 2007) and careful thermo-kinematic modelling of thermochronological datasets may allow differentiating between regional exhumation, presumably related to rock uplift, and relief change through valley incision (e.g. Valla et al., 2010, 2011). Spatial patterns of thermochronological ages also potentially allow constraining the kinematics of rock exhumation and the underlying fault geometry at depth (e.g. Robert et al., 2011; Braun et al., 2012).

In this paper, we use the thermo-kinematic modelling code Pecube (Braun et al., 2012) to quantify both river incision and structurally controlled rock exhumation, based on a dense thermochronological dataset from the Three Rivers Region. The Three Rivers Region forms the transition between the Eastern Himalayan Syntaxis and the south-east margin of the Tibetan Plateau in Yunnan and Sichuan, where the Salween (Nu), Mekong (Lancang), and Yangtze (Jinsha) rivers run closely parallel to each other over hundreds of kilometres and deeply dissect the plateau margin (Fig. 1b). The narrow spacing between the rivers has been interpreted as resulting from strong lateral shortening in response to indentation of the Indian-plate corner (Hallet and Molnar, 2001; Yang et al., 2015). This transition region, separating areas without low-relief surfaces to the west and south-west from areas with extensive low-relief surfaces to the east and north, is key to better understanding plateau-growth mechanisms and the geodynamic processes operating both within the high strain zone around the syntaxis and within the lower strain zone of the Yunnan and Sichuan margin (Fig. 1b).

In the core of this region, the Mekong River separates the Kawagebo massif (max. elevation 6740 m) to the west from the Baima Xueshan massif (max. elevation ∼ 5400 m) to the east. The Baima Xueshan massif corresponds to the southern prolongation of the elongated low-relief surface of Markam (mean elevation ∼ 4500 m), located upstream between the Yangtze and the Mekong rivers; both were mapped as relict surfaces by Clark et al. (2006) (Fig. 1c, e). The Kawagebo massif is thought to have been rapidly exhumed recently (< 10 Ma) due to motion on a local thrust fault in a restraining bend between two regional-scale strike-slip faults (Fig. 1b, c; Replumaz et al., 2020). Here, we quantify the exhumation history of the low-relief Baima Xueshan massif during the collision period to identify the different roles played by regional rock uplift and river incision, and we compare it with the exhumation history of the high-relief Kawagebo massif, while constraining the geometry of the crustal fault responsible for uplifting the latter.

Figure 1 Tectonic and geomorphic setting of south-east Tibet. (a) Topography and tectonic framework of the central and eastern Tibetan Plateau, based on the GTOPO30 digital elevation model (DEM). White box shows extent of Fig. 1b; red box shows extent of Fig. 1c. (b) Main faults, rivers, relict surfaces, and apatite (U–Th) $/$ He (AHe) data, overlain on Landsat image mosaic of south-east Tibet (USGS). (c) Relief map (modified from Zhang et al., 2015), derived from 30 m resolution DEM; 1 km radius relief is shown for pixels at elevations > 3500 m above sea level. Black outlines delimit low-relief relict surfaces with relief < 600 m, as mapped by Clark et al. (2006). Note that some averaging of the relief across the relict surface and active landscape occurs at the edges of the relict surface. In these areas local relief can be as high as 1 km. Two red rectangles represent topographic profile and relief shown in (d) and (e). (d, e) Topographic profiles across Daocheng granite plateau and Kawagebo and Baima Xueshan massifs, respectively, showing local topography and relief with available thermochronometric ages. Note that these two panels share the same legend. Apatite (U–Th) $/$ He (AHe) and apatite fission track (AFT) ages are from [1] Yang et al. (2016), [2] Liu-Zeng et al. (2018), [3] Replumaz et al. (2020), [4] Ouimet et al. (2010), [5] Clark et al. (2005), [6] Zhang et al. (2016), [7] Tian et al. (2014), [8] Nie et al. (2018), [9] Godard et al. (2009), [10] Wang et al. (2012), [11] Xiao et al. (2015), [12] Zeitler et al. (2014), [13] Gourbet et al. (2019), [14] Cao et al. (2020), [15] Shen et al. (2016), [16] Wang et al. (2018), [17] Zhang et al. (2015), [18] Reid et al. (2005), [19] Wilson and Fowler (2011), [20] Lai et al. (2007), [21] Lei et al. (2006), [22] Zhang et al. (2017), [23] Dai et al. (2013), and [24] Xu and Kamp (2000). Abbreviations: ASRRSZ – Ailaoshan-Red River Shear Zone; B – Baima Xueshan; EHS – Eastern Himalaya Syntaxis; JB – Jianchuan Basin; RRF – Red River Fault; TRR – Three Rivers Region; XSHF – Xianshuihe Fault.

2 Geologic setting

2.1 Cenozoic tectonic evolution of south-east Tibet

At the regional scale, the mean elevation between the Eastern Himalayan Syntaxis and the Mekong River (> 5000 m) is higher than the mean elevation of the low-relief surfaces east of the river (∼ 4500 m), suggesting stronger shortening and crustal thickening around the syntaxis (Fig. 1a, b). Within this zone, the Kawagebo massif, near the city of Deqing, peaks at 6740 m and towers over the plateau by more than 1700 m. Across the Mekong River in this area, a remarkable topographic and geomorphologic contrast is observed between the Kawagebo massif to the west, showing a high-topography, high-relief, and rugged landscape (Figs. 1e, 2c), and the Baima Xueshan massif to the east, characterized by an elevation close to the mean plateau elevation (∼ 4500 m) and low relief (Figs. 1e, 2b). Furthermore, deeper and larger U-shaped glacial valleys are observed in the Kawagebo massif than in the Baima Xueshan.

From the late Eocene to the early Miocene, the dominant driver of deformation in south-east Tibet was the extrusion of the Indochina block along the left-lateral Ailao Shan–Red River shear zone following the Red River (Fig. 1a), subsequently inverted along the right-lateral Red River Fault since ∼ 5–10 Ma (Leloup et al., 1995, 2001; Replumaz et al., 2001; Fyhn and Phach, 2015) (Fig. 1b). High-temperature thermochronology data from south-east Tibet shear zones show ages mainly between ∼ 34 and ∼ 17 Ma, along the Red, the Salween and the Mekong rivers, which have been linked to this extrusion process (Leloup et al., 2001; Wang et al., 2016, 2018). Compared to the timing of extrusion, no clear estimate for the timing of plateau uplift in the Three Rivers Region has been obtained. Rapid sediment filling of the Jianchuan basin (Fig. 1b), located downstream of the Three Rivers Region, around 37–35 Ma demonstrates significant erosion and suggests uplift in the source region just predating extrusion (Gourbet et al., 2017). The Jianchuan basin subsequently experienced significant exhumation along thrust faults between ∼ 28 and 20 Ma (Cao et al., 2019), suggesting regional uplift at that time.

In the Three Rivers Region, the Ailao Shan–Red River shear zone is prolongated along a distinctive and intensively sheared Jurassic red-wine-coloured clastic formation, following the Mekong River (Figs. 1b, 2a). Intense shear and deformation is attested by mostly steep to subvertical stratification and fault planes along the Mekong River, particularly in the Jurassic unit (Replumaz et al., 2020). Mostly right-lateral sub-horizontal striations have been recognized, incompatible with either motion along the left-lateral Ailao Shan–Red River ductile shear zone or thrusting. Parallel to this shear zone along the Mekong River, several thrusts affecting Eocene basins have been observed, but the timing and amount of exhumation related to activity on these thrusts have not been fully quantified (Liu-Zeng et al., 2018; Replumaz et al., 2020). The roughly north–south-trending thrusts and shear zones have been left-laterally offset by regional strike-slip faults, with the most prominent offset occurring north of the Kawagebo massif (Fig. 2a). Some of these left-lateral strike-slip faults have been reactivated as right-lateral faults, such as the Zhongdian Fault crossing the Baima Xueshan massif near Deqing.

Present-day deformation of south-east Tibet is dominated by active right-lateral strike-slip faults, such as the Jiali Fault, splitting eastward into the Parlung and Po Qu faults (Fig. 1b), which accommodate overall clockwise rotation around the Eastern Himalayan Syntaxis (Gan et al., 2007; Bai et al., 2018). In the Three Rivers Region, the Parlung Fault jumps north-eastward to the Zhongdian Fault, uplifting the Kawagebo massif in a restraining overstep since at least ∼ 10 Ma (Fig. 2a), and then jumps south-eastward to the Red River Fault, opening the Lijiang pull-apart basin (Replumaz et al., 2020).

Figure 2 (a) Geological map of the Kawagebo area in the Three Rivers Region. Active faults are in red, other major faults in black. Thermochronology data are shown with colour according to age and symbol according to method. ASRRSZ: Ailao Shan–Red River Shear Zone; K: Kawagebo massif; B: Baima Xueshan massif. (b) View of the low-relief Baima Xueshan massif at ∼ 4500 m mean elevation, with summits peaking at ∼ 5400 m visible at the background. (c) View of high-relief Kawagebo massif, with summit peaking at 6740 m, while the Mekong River valley lies at ∼ 2000 m elevation, creating almost 5000 m of relief.

2.2 Timing and quantification of exhumation and incision along the Mekong River

In the upper and lower reaches of the Mekong River, AHe ages from the valley bottom cluster around 20–15 Ma and have been interpreted as recording > 700 m of river incision in the interval ∼ 17–14 Ma, linked to the intensification of monsoonal precipitation (Nie et al., 2018). However, Oligocene or earlier (> 34 Ma) entrenchment of the lower reach of the Mekong River has previously been proposed, based on river offsets by structures related to the extrusion of the Indochina block in Burma (Lacassin et al., 1998). Regionally, a northward increase in erosion rate along the Mekong River and a westward increase toward the syntaxis since ∼ 10 Ma have been proposed (Yang et al., 2016).

In the Baima Xueshan massif, located east of the Mekong River in its middle reach, AHe ages between 5 and 16 Ma, AFT ages between 20 and 60 Ma, and zircon (U–Th) $/$ He (ZHe) ages between ∼ 80 and 120 Ma have been reported by Yang et al. (2016) and Liu-Zeng et al. (2018). Age–elevation profiles have been interpreted as recording two rapid phases of exhumation: one between 60 and 40 Ma, interpreted as being linked to the main regional uplift and crustal thickening phase; the other since ∼ 20 Ma, speculated to be linked to regional Miocene river incision (Liu-Zeng et al., 2018). However, quantitative time–temperature modelling of these age–elevation profiles using the 1D QTQt model (Gallagher, 2012) suggests a different scenario, with rapid exhumation only starting at 10 Ma and no indication of a previous phase of rapid exhumation (Replumaz et al., 2020).

In the Kawagebo massif, also along the middle reach but to the west of the Mekong River, all thermochronometric ages are much younger: AHe ages between 1 and 4 Ma, AFT ages between 3 and 7 Ma, and a single ZHe age at ∼ 8 Ma have been reported from the Mekong Valley and its tributary Yanzhi valley (Yang et al., 2016; Replumaz et al., 2020). Quantitative time–temperature inversion of the Kawagebo data suggests rapid exhumation since at least 8 Ma, with no clear estimate of the onset timing and an acceleration after ∼ 1.5 Ma (Replumaz et al., 2020). Rapid exhumation of the Kawagebo massif has been interpreted to be linked to tectonic uplift along a local thrust fault located along the Mekong River, in a restraining step over between the Parlung and Zhongdian strike-slip faults (Replumaz et al., 2020; Fig. 2). However, this local thrust fault, inferred to be related to reactivation of regional north–south-trending thrusts, has not been documented in the field. Therefore, additional work is needed in this region to resolve the exhumation history of the low-relief mean-elevation Baima Xueshan and the high-relief high-elevation Kawagebo massifs during the India–Asia collision period, in order to distinguish the effects of regional plateau uplift, incision of the Mekong River, and uplift along local tectonic structures.

Table 1 Thermal and flexural parameter values for Pecube modelling.

3 Methods: thermo-kinematic modelling of multiple thermochronology data

We synthesized published AHe, AFT, and ZHe ages in the Three Rivers Region around Deqing (Table S1, Supplement) to constrain our thermo-kinematic modelling using Pecube (Braun et al., 2012). This model is designed to solve the three-dimensional heat-transport equation in a crustal block undergoing lateral and vertical rock-particle transport, exhumation, and surface change. The flexural response to surface change is also taken into account, leading to additional rock uplift and exhumation. Default thermal and flexural parameters are shown in Table 1; some of these are optimized during the inversion. In forward mode, this set of model parameters together with a tectonic/geomorphic scenario (exhumation rates, transition times, fault geometry, surface evolution) results in rock-particle paths and an evolving 3D crustal temperature field through time, which are used to predict cooling ages for different thermochronometers. These predicted cooling ages are then compared to measured ages, defining a misfit value ϕ for this set of parameters as expressed by

$$\begin{array}{}\text{(1)}& \mathit{\varphi }=\sqrt{{\sum }_{i=\mathrm{1}}^N{\displaystyle \frac{{\left({\mathit{\alpha }}_{i,\mathrm{obs}}-{\mathit{\alpha }}_{i,\mathrm{pre}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_i^{\mathrm{2}}}}/N},\end{array}$$

where N is the number of data points, α_i,obs and α_i,pre are the observed and predicted ages for data point i, respectively, and σ_i is the uncertainty on the age for data point i. Inversion using the neighbourhood algorithm is employed to search for a best-fitting geological scenario, i.e. the set of model parameters leading to the minimum misfit value, as well as to define the resolution with which these parameter values can be constrained (Braun et al., 2012). Modelling results are represented in scatterplots for different parameter couples and as 1D or 2D marginal probability density functions (pdf's) of each parameter (Braun et al., 2012).

Geological scenarios in Pecube take into account temporal variations in exhumation rates and in paleo-topography, separately or combined. Different transition times mark the timing of exhumation or topographic changes between any two different phases. Here, we define a “steady-state topography scenario” as a scenario considering only spatially and temporally varying rock exhumation rates but sustaining the modern topography without any evolution through time. In contrast, we define an “incision scenario” as a scenario considering a simple topographic evolution from a plateau at a prescribed elevation to the modern topography, with possibly different phases separated by transition times. In order to explore the influence of incision of the Mekong River, we set up an initial plateau at 4500 m elevation, which corresponds to the average elevation of low-relief surfaces in south-east Tibet (e.g. Markam and Daocheng surfaces), evolving toward the present-day topography characterized by the deeply incised Mekong River valley (Fig. 2c). A combined scenario starting with a plateau but allowing additional rock exhumation is defined as a “plateau scenario”, where transition times mark the onset timing of river incision from the initial plateau as well as of different exhumation phases. For the Kawagebo massif, in order to model spatial variations of uplift above a crustal-scale fault rooting at a mid-crustal depth, a 3D fault geometry is implemented following Robert et al. (2011), defining a “tectonic scenario”. In this scenario, the topography is in steady state, and rock exhumation is due to motion along the fault. The fault trace at the surface and the number of segments at depth are fixed parameters in this inversion. The geometry of the fault is represented by different segments, separated by several inflection points, which are defined by their coordinates (r, s), corresponding to the depth and distance relative to the surface trace of fault, respectively. Those coordinates and velocities are inverted to constrain the fault geometry and exhumation history best fitting the thermochronology dataset.

Figure 3 Model results for the Baima Xueshan massif. (a) Age–elevation profiles of measured (solid blue symbols; triangles, squares, and diamonds corresponding to AHe, AFT, and ZHe ages, respectively) and predicted (open symbols coloured according to scenario) ages for the steady-state topography and incision scenarios. (b) Age–elevation profile of measured AHe and AFT ages (< 22 Ma) and predicted ages for the plateau scenario and predictions of the scenarios shown in (a) for comparison. (c) Three-dimensional view of the present-day topography (at resolution of 900 m) of the Baima Xueshan massif above the Mekong River, coloured according to predicted AHe ages for the plateau scenario, the best-fit model for Baima Xueshan transect (orange symbols in b). Note that the plateau is shown slightly higher than its real elevation for clarity. The squares indicate the sample locations.

Figure 4 Scatterplots of Pecube inversion misfits for the plateau scenario model of the Baima Xueshan massif. Coloured dots represent individual forward-model runs plotted in two-dimensional projections of the parameter space, with colours corresponding to misfit values, shown on a log scale. Posterior probability density functions for parameter values are plotted along the axes. The best-fit solution is represented by a yellow star with 2σ (dark red) and 1σ (light red) confidence contours. (a) Transition time (between first and second exhumation phase) versus exhumation rate during the first phase. (b) Transition time versus exhumation rate during the second phase.

4 Result

The strong contrasts in elevation, relief, and thermochronological ages between the Baima Xueshan and Kawagebo massifs on the two flanks of Mekong River led us to model their exhumation histories separately. We tested steady-state, incision, and plateau scenarios for both massifs, to resolve the influence of Mekong River incision on their exhumation history (Table 2). For the Kawagebo massif, we additionally tested a tectonic scenario implementing a thrust fault with a surface trace along the Mekong River, as recently suggested by Replumaz et al. (2020), aiming to determine its geometry and the contribution of its activity to exhumation in the massif.

Table 2 Inversion parameters for the different scenarios tested for the Baima Xueshan and Kawagebo massifs. The best-fit scenario, with the lowest misfit for each massif, is in bold. A range of exhumation rates is explored for different phases, separated by transition times. Coordinates of the fault inflection points define the geometry of the fault plane, with r_i and s_i (in km) corresponding to the depth and distance relative to the surface trace of fault, respectively. Note that the tectonic scenario includes a slip rate on the fault rather than a vertical exhumation velocity and that upward slip on a thrust is indicated by negative values (see Braun et al., 2012, for more details).

4.1 Baima Xueshan massif

The exhumation history of the Baima Xueshan massif has been modelled using the available data for three thermochronometers (AHe, AFT, ZHe; Fig. 3a, Table S1) from Yang et al. (2016) and Liu-Zeng et al. (2018). Models were run for 110 Ma to encompass the full history recorded by these data. We initially ran steady-state topography scenarios (Table 2). A model with two exhumation phases, with protracted slow exhumation at a rate of 0.04 km/Myr from 110 to 7 Ma followed by rapid exhumation at a rate 0.42 km/Myr since 7 Ma, fits the entire dataset with a relatively low misfit (3.29). However, the break in slope in AFT ages around 40 Ma, marked by three ages at elevations between ∼ 4000 and ∼ 4700 m and interpreted by Liu-Zeng et al. (2018) as recording a rapid phase of exhumation, is not reproduced (Fig. 3a). Moreover, AHe ages > 10 Ma, marking a subtle break in slope around 4500 m elevation, and AFT ages < 25 Ma are not well modelled either in this steady-state topography scenario (Fig. 3a, b). To test the existence of a rapid exhumation phase around 40 Ma, as proposed by Liu-Zeng et al. (2018), we modelled more complex scenarios with three or four exhumation phases and specified transition times. The best-fitting of these models comprises a four-phase scenario (Table 2, Fig. S1), resulting in a slightly lower misfit (3.13) than the previous two-phase scenario. However, none of the more complex models resolve a rapid exhumation phase around 40 Ma; they all result in best-fit models characterized by slow exhumation that accelerated in the last few million years (< 10 Ma), very similar to the two-phase exhumation scenario (Table 2). We conclude that the data do not require rapid exhumation prior to the last few million years and we do not further consider these more complex scenarios.

Second, we ran an incision scenario, comprising a plateau at an elevation of 4500 m at the beginning of the model (110 Ma), which linearly evolves toward the present-day topography after a transition time. The best-fitting model shows a transition time at ∼ 10 Ma (Table 2). For this endmember incision scenario, only the lowest-elevation AHe and ZHe ages are relatively well predicted, while the predicted AHe ages at > ∼ 3300 m elevation and the AFT ages do not fit the observed ages (Fig. 3a), as expressed by a high overall misfit of 12.6. The high misfit and overprediction of all but the lowest-elevation ages imply that the total amount of exhumation is underestimated in this scenario.

We finally tested a plateau scenario, including both regional rock uplift and incision, fitting only the AHe and AFT ages younger than 22 Ma to concentrate on the Neogene history, considering that the steady-state models, including the old ZHe ages, adequately predict the early history. The best-fit model includes initial protracted slow exhumation at a rate of 0.01 km/Myr, a transition time at ∼ 10 Ma, marking the onset time of river incision, and more rapid exhumation at a rate of 0.25 km/Myr since ∼ 10 Ma. This model has the lowest misfit (1.16), with well-resolved parameter values (Fig. 4). Note that this lowest misfit (1.16) cannot be directly compared to the misfit of the best steady-state scenario (3.29), which models a longer timescale (110 Ma) with more data (AFT ages > 22 Ma and ZHe ages). However, this plateau scenario better reproduces the AHe and AFT age–elevation trends, including the slope break around ∼ 10 Ma in the AHe age–elevation plot (Fig. 3b). Maximum river incision in this scenario from a plateau of 4500 m is 2.5 km for the Mekong Valley bottom (currently at 2000 m elevation); the average contribution of river incision to exhumation of the modelled samples in the valley (average elevation 3500 m) is ∼ 1 km, accounting for ∼ 30 % of the total exhumation of ∼ 3.5 km (∼ 1 km topographic change and 2.5 km of rock uplift) since 10 Ma (see Fig. 7).

4.2 Kawagebo massif

Published thermochronology data from the Kawagebo massif are dispersed over more than 40 km along the Mekong River valley (Fig. 2). To the north, several samples yield ages with large error bars, with AHe ages clustering between 2 and 5 Ma, AFT between 5 and 30 Ma, and ZHe ages between 7 and 50 Ma (Wilson and Fowler, 2011; Liu-Zeng et al., 2018; Replumaz et al., 2020). To the south, near the village of Yanzhi, along a west-bank tributary incising the Kawagebo massif, samples show more coherent ages, with AHe ages between 1 and 4 Ma, AFT ages between 3 and 7 Ma, and a single ZHe age at ∼ 8 Ma (Yang et al., 2016; Replumaz et al., 2020). We use this dataset from the Yanzhi valley to constrain thermo-kinematic modelling of the massif evolution since 10 Ma. The steady-state topography scenario predicts a transition at 1.2 Ma between slower exhumation (0.54 km/Myr) and rapid exhumation (1.75 km/Myr), fitting the data reasonably well (misfit = 1.05). However, this scenario does not reproduce the increasing AFT ages with elevation nor the observed slope break in AHe ages in the age–elevation plot (Fig. 5a). As the high-elevation samples were collected toward the core of the massif, farther from the Mekong River, these characteristics of the age–elevation plot could potentially be explained by spatially varying exhumation rates due to a curved fault at depth and/or incision of the Mekong River. To explore these possibilities, we include tectonic particle transport along a thrust fault striking along the Mekong River and incision in the following inversions.

Figure 5 Modelling results for the Kawagebo massif. Age–elevation profiles of measured and predicted ages (same legend as Fig. 4) for (a) the steady-state topography scenario and the tectonic scenario with a planar (one-segment) thrust; (b) the tectonic and the tectonic + plateau scenarios, both with a kinked (two-segment) thrust. (c) Three-dimensional view of the present-day topography of the Kawagebo massif and the Mekong River valley, coloured according to predicted AHe ages of the tectonic + plateau scenario with a two-segment thrust. The geometry of the thrust and rock-particle trajectories are shown at depth, with squares indicating the sample locations. The thrust geometry is defined by P₁, P₂, and P₃ and their coordinates in distance relative to the surface trace of thrust. This scenario corresponds to the best-fit model for the Kawagebo massif (orange symbols in b). Note that the plateau is shown slightly higher than its real elevation for clarity. The arrows show the present-day velocity field generated by movement along the fault, during the recent rapid exhumation phase since 1.5 Ma; the two-segment geometry leads to spatially varying exhumation rates (shown by colours, scale on the right).

Figure 6 Scatterplots of PECUBE inversion results for two-segment thrust scenario model of the Kawagebo massif. Coloured dots represent single forward-model runs plotted in two-dimensional projections of the parameter space, with colours corresponding to misfit values, shown on a log scale. Posterior probability density functions for parameter values are plotted along the axes. The best-fit solution is represented by a star with 2σ (dark red) and 1σ (light red) confidence contours. (a) Depth of the base of the flat (s₁) versus the horizontal location of the base of the ramp, relative to the surface trace of the thrust (r₂). (b) The horizontal distance of the base of the ramp (r₂) versus the depth of the base of the ramp (s₂). For explanation of these parameters, see Fig. 5c. (c) Transition time versus slip rate along the fault during the first phase. (d) Transition time versus slip rate along the fault during the second phase.

A tectonic scenario with a simple planar (one-segment) fault converges to a similar misfit (1.08) as the steady-state topography scenario (1.05). The best-fit model shows a fault dipping 85^∘ to the west and a transition time at 1.3 Ma, with a slow slip rate of 0.58 km/Myr before and a faster rate of 1.87 km/Myr after. Due to the steep dip of this reverse fault, this tectonic scenario is very similar to the steady-state topography scenario considering purely vertical motion (Fig. 5a).

A more complex tectonic scenario with a kinked (two-segment) thrust fits the data slightly better (misfit = 0.96) but much better reproduces the observed age–elevation relationships (Fig. 5b). In this scenario, the best-fit fault geometry is a thrust with a steep shallow segment (dipping at 65^∘) and a near-horizontal (1.5^∘ dip) deeper segment, below 15 km (Fig. 5c). The transition time is 1.5 Ma, with a slip rate of 0.69 km/Myr before, increasing to 1.84 km/Myr after. The velocity field predicted by this model implies decreasing exhumation rates away from the Mekong Valley, as high exhumation rates above the steep shallow segment of the thrust give way to lower rates above the deep flat segment (Fig. 5c). The low slope in the AFT age–elevation plot and the slope break in AHe ages are reproduced due to this spatial variation in exhumation rates. Input parameters of this scenario are relatively well constrained, as shown by the misfit scatterplots and the posterior pdf's of parameter values (Fig. 6).

For consistency with the modelling results of the Baima Xueshan massif across the Mekong River, we also ran a “tectonic + plateau” scenario, adding the incision of a 4500 m high plateau to the tectonic scenario (two-segment thrust). This tectonic + plateau scenario converges to the lowest misfit (0.66), with a similar transition time (1.6 Ma) but somewhat smaller fault-slip rates (0.45 km/Myr until 1.6 Ma, 1.86 km/Myr thereafter) than the tectonic scenario with two thrust segments (Fig. 5b). The average contribution of river incision from a plateau of 4500 m to exhumation of the modelled samples (average elevation 2500 m) is ∼ 2 km or ∼ 25 % of the total exhumation of ∼ 8 km (∼ 2 km topographic change and 6.1 km of rock uplift) since 10 Ma (Fig. 7), a proportion similar to that of Baima Xueshan. The similarity between scenarios with and without topographic change suggests that river incision plays a minor role in the exhumation history of the Kawagebo massif, characterized mainly by tectonic forcing.

Figure 7 Exhumation since ∼ 10 Ma for the Kawagebo (red arrows) and Baima Xueshan (green arrows) massifs; with inferred Mekong river incision (downward red/green arrows) from scenarios with a schematic plateau at 4.5 km (horizontal grey thick line). Purple dotted thick line marks the glacial equilibrium line altitude (ELA) during the Last Glacial Maximum (LGM; ∼ 20 ka) (Fu et al., 2013), at which glacial erosion is most effective. Glaciation allows the formation of low-relief surfaces at high elevation through buzzsaw-like processes in Baima Xueshan, where rock exhumation is relatively low compared to the Kawagebo massif. Note that the geometry of the west-dipping thrust along Mekong river is not scaled at depth.

5 Discussion: tectonic and climatic forcing on exhumation

5.1 Tectonic and climatic forcing on the evolution of the low-relief Baima Xueshan massif

For the Baima Xueshan massif, east of the Mekong River, our best-fit model, exploring the exhumation history since 110 Ma, shows rapid regional exhumation at a rate of 0.42 km/Myr since ∼ 7 Ma, succeeding a phase of slow regional exhumation at a rate of 0.04 km/Myr before that. The onset of rapid exhumation occurred slightly later compared to published thermal history models using subsets of the data, which showed well-defined rapid exhumation at a rate of 0.26 km/Myr since 10 Ma (Replumaz et al., 2020). Our results do not support the history inferred by Liu-Zeng et al. (2018), obtained by using pseudo-elevation profiles (Reiners and Brandon, 2006), with rapid exhumation between 100 and 80 Ma at rate of 0.15 km/Myr, between 60 and 40 Ma at a rate of 0.6 km/Myr, and since ∼ 20 Ma at a rate between 0.22 and 0.65 km/Myr. Surprisingly, both the previous 1D thermal-history modelling of Replumaz et al. (2020) and the 3D thermo-kinematic modelling reported here imply that none of the tectonic events occurring in south-east Tibet before 7 Ma have been recorded as exhumation phases in the Baima Xueshan massif. Those events include extrusion of the the Indochina block between ∼ 34 and ∼ 17 Ma (e.g. Leloup et al., 2001), significant erosion in the source region recorded by rapid sediment filling in Eocene basins downstream of the Three Rivers Region between 37–35 Ma (Gourbet et al., 2017), and shortening of these basins between 28 and 20 Ma (Cao et al., 2019). These results suggest that tectonic forcing in this massif has been negligible in driving its exhumation history before 7 Ma.

Our best-fit model focusing on the Neogene history of the massif since 22 Ma implies the incision of a plateau at an initial elevation of 4500 m since ∼ 10 Ma, associated with regional exhumation at a rate of 0.25 km/Myr, succeeding a phase of slow regional exhumation at a rate of 0.01 km/Myr until 10 Ma. In this model, about 1 km of river incision (at the elevation of the samples) is added to 2.5 km of regional uplift-driven exhumation, resulting in a total exhumation of ∼ 3.5 km since 10 Ma (Fig. 7). Using such a simplified topographic evolution scenario mimicking the incision of a plateau, we show that differential erosion between the low-relief, mean-elevation Baima Xueshan massif and the deeply incised Mekong River valley is necessary to reproduce the break in slope in the AHe age–elevation profile, between samples on the plateau and samples from the Mekong River valley (Fig. 3c). This scenario of incision of the Mekong River since ∼ 10 Ma is compatible with the onset of incision estimated between 13 and 9 Ma for the Dadu and Yalong rivers gorges to the north-east (Clark et al., 2006; Ouimet et al., 2010).

Our models show that river incision alone is not sufficient to reproduce the AHe ages at higher elevations, as well as the AFT ages from the Baima Xueshan massif (Fig. 3a), and accounts for ∼ 30 % of the total exhumation since ∼ 10 Ma. For the highest samples on the plateau, our results imply that ∼ 2.5 km of overburden was removed from the plateau surface since 10 Ma, which is an order of magnitude higher than the exhumation (∼ 0.23 km) estimated for the low-relief Daocheng plateau to the north-east during the same time span (Clark et al., 2005). This higher exhumation is required by the much younger AHe and AFT ages in the Baima Xueshan massif (< 18 and < 56 Ma, respectively; Liu-Zeng et al., 2018) compared to the Daocheng massif (> 50 and > 100 Ma, respectively; Clark et al., 2005). Therefore, the Baima Xueshan massif cannot be considered to represent the same relict surface as the Daocheng granite, despite the fact that both have low relief (Fig. 1). We conclude that a relict surface sensu stricto, i.e. a remnant of a paleo-landscape that was barely affected by exhumation during the India–Asia collision (Clark et al., 2005), should show AHe ages older than 50 Ma, and that not all low-relief surfaces in south-east Tibet can be classified as a relict surface.

Figure 8 Three-dimensional view of the Kawagebo (K) and Baima Xueshan (B) massifs (Landsat image draped on DEM). Active faults are in red, other major faults in black. Thermochronological ages are shown with colour representing age and symbol representing method. MK samples are from Yang et al. (2016), DQ from Liu-Zeng et al. (2018) (Baima Xueshan transect), and KW from Replumaz et al. (2020) (Kawagebo transect). Also shown are the crustal geometries of faults (in red) active in the last 10 Myr. Arrows annotated B or C: positions of field photos shown in Fig. 2.

Figure 9 (a) Simplified geological structures of south-east Tibet, with GPS velocity field (modified from Wang et al., 2017), showing a N–S component of motion north of India. Blue patch is the area of average elevation above 5000 m on the Tibetan Plateau. The reconstructed position of Indian continent at 10 Ma (green contour from Replumaz and Tapponnier, 2003) shows continuing convergence since that time (green north-directed arrows). (c) Correlation between thermochronological (AHe and AFT) ages in south-east Tibet and topography above 5000 m (blue contour). White square represents our study area. More additional data are from Burg et al. (1998), Ding and Wang (1995), Seward and Burg (2008), Tu et al. (2015), Yu et al. (2011), and Zhang et al. (2015) (see Fig. 1 for references of other thermochronological ages). Abbreviations: GF, Ganzi Fault; ITS, Indus–Tsangpo Suture; JF, Jiali Fault; LMS, Longmenshan thrust system; LSS, Longmucuo–Shuanghu Suture; MBT, Main Boundary Thrust; MCT, Main Central Thrust; PF, Parlung Fault; RRF, Red River Fault; SF, Sagaing Fault; STDS, South Tibet Detachment System; XJF, Xiaojiang Fault; XSHF, Xianshuihe Fault; ZF, Zhongdian Fault.

The significant removal of overburden (∼ 2.5 km) from the Baima Xueshan massif surface (currently at a mean elevation of ∼ 4500 m and peaking at ∼ 5400 m) since ∼ 10 Ma implies the efficient erosion processes at high elevation have generated the currently observed low-relief landscape (Fig. 1). Glacial erosion and vigorous periglacial processes have been shown to be efficient in smoothing high-elevation regions above the glacial equilibrium line altitude (ELA) (Brozovic et al., 1997; Egholm et al., 2009, 2017; Hales and Roering, 2009). Zhang et al. (2016) have proposed that these processes could be active in south-eastern Tibet and could lead to seemingly continuous low-relief high-elevation surfaces despite spatially differential and diachronous tectonic exhumation of 2–4 km between the surfaces. The mean ELA across south-east Tibet during the Last Glacial Maximum (LGM, ∼ 20 ka) was estimated to be ∼ 4.6 km, much lower than the present-day ELA at ∼ 5.4 km (Fig. 7; Fu et al., 2013). The ELA during previous glacial phases of the Quaternary would have been comparable to that of the LGM, with the average Quaternary ELA lying somewhere between these two values. Numerous cirques, moraines, and U-shaped valleys are observed across the Baima Xueshan massif, providing evidence for significant glacial erosion (Fig. 8). Based on this observation as well as the maximum elevation of the massif (5400 m), which is close to the present-day ELA (5400 m), we suggest that buzzsaw-like (peri-) glacial processes could be active in the Baima Xueshan massif, such that any tectonic uplift bringing the elevation of the plateau above the ELA could trigger glacial processes that would erode and smooth these highlands, as previously proposed by Zhang et al. (2016). These glacial processes since at least ∼ 2 Ma could thus contribute, in part, to removing the overburden of 2.5 km of materials, coupled with other erosional processes since 10 Ma.

Increasing regional rock uplift toward the eastern Himalayan Syntaxis has been inferred from relatively young thermochronological ages (AHe/AFT ages < 15 Ma e.g. Burg et al., 1998; Replumaz et al., 2020; Ouimet et al., 2010) in a region between the South Tibet Detachment System–Parlung Fault and the Longmucuo–Shuanghu Suture (Fig. 9), which roughly follows the Mekong River. This zone is also characterized by relatively high regional topography, with the plateau surface generally exceeding 5000 m (Fig. 9). In contrast, the area to the east of the Longmucuo–Shuanghu Suture, including the low-relief surfaces of Markam and Daocheng, show lower average elevations and older thermochronological ages (AHe/AFT > 15 Ma). This zone between the South Tibet Detachment System–Parlung Fault and the Longmucuo–Shuanghu Suture is thus potentially associated with regionally enhanced late Miocene uplift north and east of the Eastern Himalayan Syntaxis (EHS), as argued for instance by Zeitler et al. (2014) and Schmidt et al. (2015), possibly due to the continuous northward advance of the eastern indenter corner of the Indian plate. In this zone, a north–south component of motion is also observed from the current GPS velocity field, which shows ongoing north–south shortening north of India (Fig. 9) associated with surface uplift (Liang et al., 2013). The Baima Xueshan massif is located in the south-eastern-most extremity of this regional uplift zone. Within the massif, N–S-oriented thrusts like the Zigaishi thrust, east of and parallel to the Mekong River, are associated with a clear topographic crest bounding the massif to the west and could have been reactivated at a moderate rate (0.25 km/Myr) to accommodate shortening and thickening of the Three Rivers Region since ∼ 10 Ma (Figs. 2a and 8). Moderate uplift along those thrusts could be due to the distance from the EHS, with a west-to-east gradient of decreasing exhumation and erosion rates across the Three Rivers Region (Yang et al., 2016).

5.2 Dominant tectonic forcing on exhumation in the Kawagebo massif

Our best-fit model for the Kawagebo massif combines incision of a 4500 m high plateau and tectonically controlled rock uplift above a segmented thrust (Figs. 5 and 7). Two phases have been modelled in this tectonic + plateau scenario, with an initial slip rate on the fault of 0.45 km/Myr from at least 10 Ma until 1.6 Ma, followed by rapid slip at a rate of 1.86 km/Myr, leading to 6.1 km of tectonically driven exhumation since 10 Ma (Fig. 7). The available dataset only covers a period since 10 Ma and cannot constrain the onset timing of the initial exhumation phase. Our model outcomes corroborate earlier thermal-history modelling of the samples from the Yanzhi valley in the Kawagebo massif using QTQt, which suggested rapid cooling from 8 Ma and a significant acceleration from 1.5 Ma to the present (Replumaz et al., 2020). The average topographic lowering due to river incision at the elevation of the samples is ∼2 km, accounting for ∼ 25 % of the total exhumation of the massif (Fig. 7). This result implies that tectonic forcing has been dominant in exhuming the Kawagebo massif since at least 10 Ma. Potentially, the inferred acceleration of exhumation since 1.6 Ma could be due in part to more efficient erosion related to the onset of high-altitude glaciations and/or monsoon intensification. Efficient glacial erosion in the Kawagebo massif is attested to by deep U-shaped valleys, like the Yanzhi valley, suggesting > 1 km of glacial erosion (Fig. 8). In contrast to the Baima Xueshan massif, glacial erosion has not led to relief reduction in the Kawagebo massif; this difference could be due to the much higher rock uplift rate of up to ∼1.7 km/Myr in Kawagebo versus 0.25 km/Myr in Baima Xueshan, consistent with worldwide comparisons of glacial relief in mountain belts (Pedersen et al., 2010).

Combining the 3D distribution of sample ages relative to the fault trace, fault-controlled exhumation is predicted to be linked to a steep (∼ 65^∘) thrust in the upper crust, which flattens at depth and strikes roughly parallel to the Mekong River (Fig. 5c). This fault geometry reproduces the low slope in the AFT age–elevation profile as well as the observed slope break for AHe ages (Fig. 5b, c). Replumaz et al. (2020) suggested that late Miocene uplift and exhumation of the Kawagebo massif, which peaks at ∼ 2300 m above the mean plateau elevation of the plateau, was accommodated by such a fault. The fault was suggested to have activated at ∼ 10 Ma in a large-scale restraining left-stepping overstep between the right-lateral Parlung and Zhongdian strike-slip faults (Fig. 8), resulting in more pronounced local uplift compared to the regional uplift recorded in the Baima Xueshan massif. Steep to subvertical stratification and fault planes are observed along the Mekong River in the sampling area (Replumaz et al., 2020), in agreement with the modelled steeply westward-dipping thrust fault. However, observed striations are mostly sub-horizontal rather than dip slip, suggesting that the active Mekong thrust is probably not outcropping in the field.

6 Conclusions

Using 3D thermo-kinematic models constrained by a relatively dense low-temperature thermochronology dataset, we have compared different exhumation scenarios for massifs along the middle reach of the Mekong River in the Three Rivers Region. We are able to discriminate between tectonically controlled rock uplift and river incision for both the Baima Xueshan and Kawagebo massifs. This method allows refining and validating the first-order conclusions on exhumation history inferred from age–elevation relationships and assessing the relative contributions of tectonic or climatic forcing to inferred phases of rapid exhumation.

We show that different tectonic processes dominate the exhumation histories of the two flanks of the Mekong River. To the east, the low-relief Baima Xueshan massif is not a relict surface as defined by Clark et al. (2006), as it has experienced an amount of exhumation (∼ 2.5 km since ∼ 10 Ma) that is an order of magnitude higher than that determined for the emblematic Daocheng relict surface (∼ 0.23 km). None of the tectonic events affecting south-east Tibet before 10 Ma has been recorded as exhumation phases in the Baima Xueshan massif, suggesting that tectonic forcing on the massif's exhumation has been negligible during that time. Moderate rock exhumation at a rate of 0.25 km/Myr since ∼ 10 Ma was mainly triggered by the northward advance of the indenting Indian plate, leading to higher regional elevation around the EHS including the Baima Xueshan, with the incision of the Mekong River contributing only ∼ 30 % of the exhumation recorded by the samples. The contrasting exhumation history of the Kawagebo massif to the west is dominated by strong local uplift driven by a segmented thrust dipping westward in a regional compressional overstep. Total tectonically driven exhumation since ∼ 10 Ma at a rate of 0.45 km/My, with an acceleration to 1.86 km/My since 1.6 Ma, is estimated to be 6.1 km, with the incision of the Mekong River contributing only ∼ 25 % of the total exhumation.