Introduction

The landscape of an orogenic belt is constantly shaped by competition between tectonic uplift and bedrock river incision1. It is generally accepted that tectonic uplift increases stream gradient, which conveys more water and sediment via orographic precipitation that facilitates faster incision2,3,4,5, driving the headward propagation of river knickpoint6,7,8. However, why rivers have not incised further or more deeply into active orogenic plateaus to destroy uplifted low-relief terrain over geologic time remains uncertain. The traditional explanation is that pronounced aridity in the interior plateau defeats river incision, and concomitant formation of internally drained basins9,10. Sustained internal drainage and sedimentary infilling of basins11,12 reduce elevation contrasts, raising local base levels13. In addition, the bidirectional coupling between bedrock river incision and landslides and the resulted drop in landslides frequency through time could lead to a corresponding decrease in the rate of fluvial incision14. Further, highly local rock uplift15,16 or derformation17, glacier and/or landslide damming18,19,20, or dynamic feedbacks between tectonics uplift and erosion (i.e. ‘tectonic aneurysm’ model)21,22 on the plateau margins have been also viewed as constraining features for the dissection of the plateau morphology by rivers. Thus, the controls on the erosional decay of orogenic plateaus remain controversial.

This debate exists especially in the southern Tibetan Plateau, where is characterized by a relatively low-relief and high-elevation landscape (Fig. 1a, b). The Yarlung River, the largest river in the southern Tibetan Plateau, flows eastward along the Indus-Yarlung Suture zone and then turns southward around the Eastern Himalayan Syntaxis (EHS), creating one of the deepest gorges (Fig. 1b) on Earth16,21,22,23. Evidence from paleo-altimetry data24 and the sediments in the Himalaya foreland25 and Bengal basins26 show that similar to present-day external drainage systems and high elevation in the southern Tibetan plateau persisted since at least 15 Ma. Meanwhile, the Tibetan plateau since the middle Miocene also underwent active east-west extension23, resulting in a series of N-S trending rifts developed across the plateau, especially in southern Tibetan Plateau. However, it is uncertain how high elevations in the southern plateau have been sustained in the face of river erosion and fault activity.

Fig. 1: Regional background.
figure 1

a Topography of the Tibetan Plateau and its major rivers. b Major faults, rivers and topography of the southern Tibetan Plateau. The red dots represent large knickpoints in the Yarlung River. Dashed rectangle shows the extent of Fig. 1c. The whites square shows the location of the Eastern Himalayan Syntaxis. BNS, Bangong-Nujiang Suture zone; IYS, Indus-Yarlung Suture zone. c Simplified geological map76 and sample locations along the eastern Yarlung River. d New and published thermochronological ages33,35,36 (Supplementary Table 4) along the eastern Yarlung River. AHe, apatite (U-Th-Sm)/He; AFT, apatite fission track; ZHe, zircon (U-Th)/He. Vertical error bars indicate 1σ uncertainty for thermochronological ages. The black and blue lines represent river profile and valley floor width, respectively. The river profile was extracted from ~30 m Shuttle Radar Topography Mission (SRTM) digital elevation model data, smoothing profiles used a LOESS filter to decrease noise and highlight the general shape of the river profile. The valley floor width was measured from Google Earth Image (Supplementary Note 1 and Supplementary Fig. 9).

Full size image

Thermochronological ages show a younging trend from the central part of the Tibetan Plateau to the margins (i.e., northeastern and southeastern margins), suggesting that the plateau may have progressively grown and propagated northeastwards and southeastwards27,28. This expanding uplift has been thought to accelerate river incision from the southeastern margin to the interior plateau29,30,31,32, driving the propagation of river knickpoints into the plateau33. However, some large knickpoints along river valleys coincide in space with local tectonic deformation. For example, knickpoints along the Yarlung River locally coincide with the late Cenozoic N-S trending tectonic rifts (Fig. 1b) that accommodate E-W extension in southern Tibetan Plateau34. This suggests that formation of these knickpoints may have been controlled by active fault-flank uplift. Whether the diverse activity of regional tectonics could have also inhibited knickpoint migration and thus reduced river erosion to maintain topography in southern Tibetan Plateau remains unclear.

To further address this issue, the incision history of the eastern Yarlung River valley in the southern Tibetan Plateau was constrained by low-temperature thermochronology. We collected 23 granite samples along the Gyaca gorge (knickpoint) and the broad valley of the Yarlung River (Supplementary Note 1). A total of 22 apatite fission track (AFT), 6 apatite and 8 zircon (U-Th)/He (AHe and ZHe) ages were obtained (Supplementary Tables 1–3), and previously published thermochronological ages33,35,36 in the eastern Yarlung River valley were also compiled (Fig. 1c, Supplementary Table 4 and Supplementary Fig. 1a). Three-dimensional (3D) thermokinematic modeling was then undertaken to explore the effect of the N-S trending rift on the formation and evolution of the knickpoint and the incision of the Yarlung River. Combined with previous studies, we proposed that fault systems control knickpoints migration and impede river incision in the southern Tibetan Plateau since the late Miocene.

Results

Thermochronological data

New and published thermochronological ages (Supplementary Table 4) along the eastern Yarlung River (from west to east spanning ~330 km) are presented in Fig. 1d (see Supplementary Fig. 1a for further details). Zircon (U-Th)/He (ZHe), apatite fission track (AFT), and apatite (U-Th)/He (AHe) ages range from 23.1 ± 1.4 to 7 ± 1.1 Ma, 17.5 ± 2.6 to 4.2 ± 0.4 Ma, and 9.25 ± 1.49 to 1.75 ± 0.22 Ma, respectively (Supplementary Table 4). The ZHe and AFT ages, except for the relatively young ages in the Gyaca gorge, generally show older ages in the upstream of this gorge compared to those downstream (Fig. 1d). The AHe ages show remarkable youngest ages in the Gyaca gorge and older ages downstream and upstream of this gorge (Fig. 1d).

Thermal history

Cooling histories of the new and previous samples in the study area (divided into 6 Zones, see Fig. 2) were reconstructed using inverse thermal history modelling by the QTQt software37 (Figs. 2a-2f). The spatial and temporal variations of the rock cooling rate since ~12 Ma were calculated from the modeling results (Fig. 2g). The results indicate rapid cooling of the basement rocks in the eastern Yarlung valley between ~12 Ma and ~7 Ma, but the cooling rate has increased obviously in the Gyaca gorge and decreased dramatically in its upstream and downstream since ~7 Ma (Fig. 2g).

Fig. 2: Inverse thermal history models from QTQt based on available thermochronological data (Supplementary Table 4) along the eastern Yarlung River.
figure 2

a–f Cooling histories of bedrock samples from west to east along the Yarlung River. Zone 1, upstream of the Gyaca gorge (hanging wall of the Woca normal fault); Zone 2, the Gyaca gorge (footwall of the Woca normal fault); Zone 3-6, downstream of the Gyaca gorge. Different colored lines indicate best-fit modelled thermal paths for different samples. g Spatial and temporal pattern of cooling rates calculated from thermal history modeling results. Black error bars denote the standard deviation using three different cooling rates (see Supplementary Note 4 for further details).

Full size image

3D thermo-kinematic modeling

To further explore possible mechanisms governing the regional erosion history, 3D thermo-kinematic modeling strategy using Pecube was applied (detailed setting of the models described in Supplementary Table 5). We constructed two scenarios, i.e. one without (Scenario A), and then with the influence of the Woca fault that bounds the rift structure (Scenario B) (Supplementary Fig. 2). Based on our thermal history modeling results above, the first model (Scenario A) was set-up with a two-stage (12-7 Ma and 7-0 Ma) evolution history without Woca fault activity. More than 20,000 iterations were run to find the best fitting model that constrains the incision rate of the Yarlung River. The best-fitting results exhibit a first stage of rapid incision rates of 1.98 km Ma−1, followed by a second stage with a decreased rate of 0.45 km Ma−1 (Supplementary Table 6 and Supplementary Fig. 3). The model predicted thermochronological ages, which are generally consistent with those actually observed for samples located > 25 km away from the Gyaca gorge. However, for those samples located nearby the Gyaca gorge (in the footwall of the Woca fault), the predicted AFT and AHe ages are much older than those measured (Fig. 3a, e). These large differences suggest that Scenario A, which did not integrate normal fault movements, was unable to fit the spatial distribution of observed AHe and AFT ages (Fig. 3a, e).

Fig. 3: Results of thermo-kinematic inversion and forward modeling.
figure 3

a–d The vertical velocity of thermo-kinematic inversion during the recent exhumation phase (7-0 Ma) and the predicted AHe ages for different scenarios. The arrows reflect the direction of movement of rock particles with (a) and without a normal fault (b, c, d). The yellow dots represent inflection points defining fault geometry. e–h Comparison between observed and predicted cooling ages along the eastern Yarlung River from the corresponding forward modeling of different scenarios since 12 Ma. Vertical error bars indicate 1σ uncertainty for observed age.

Full size image

For Scenario B, we added a normal fault to simulate a more complex uplift velocity field. Based on our thermal history modeling results, the time of initial normal fault activity was set at 7 Ma (Fig. 2g), and the fault trace at the surface follows that of the Woca fault (see Fig. 1c). Due to the lack of geophysical data constraining the deeper sections of the fault, we ran three distinct inversions in which we fixed the deepest segment to be horizontal at depths of 20, 30, and 50 km, respectively. The results show that the three models for Scenario B successfully reproduced the AHe and AFT ages in the Gyaca gorge (Fig. 3b–d, f–h), and that the spatial distribution of predicted ages is generally consistent with measured ages (Figs. 3b–d, f–h). This confirms that incision of the eastern Yarlung River was indeed strongly influenced by a normal fault since ~7 Ma.

Comparing the misfit values of the inversion for the three models in scenario B (Fig. 3f–h, Supplementary Table 6), we found that values for the model coupled with a 30 (3.92) or 50 km-depth fault (4.07) seem a little lower than that for a 20 km-depth fault (4.17). This would indicate that the depth of the Woca fault could reach up to ~30 km or more. Furthermore, these models predict a steep (68°∼75°) fault in the shallow crust, near the surface (Fig. 3c, d), and a fault-slip velocity rate of 1.19 ± 0.03 km Ma−1 since 7 Ma (Supplementary Table 6 and Supplementary Figs. 4–6). This finding is in very good agreement with the observed high-angle fault geometry of 60° ~ 75°38 and the Quarternary slip rate of 1.2 ± 0.6 km Ma−139. The model also shows that the incision rate along the eastern Yarlung River was up to 1.73 ± 0.07 km Ma−1 between 12 Ma and 7 Ma, followed by a rather slow incision rate of 0.18 ± 0.03 km Ma−1 upstream and downstream of the Gyaca gorge since 7 Ma (Supplementary Table 6 and Supplementary Figs. 4–6). This is also broadly consistent with the thermochronology results.

Discussion

Fault activity controls the formation of knickpoint

In a previous study33, Schmidt et al. proposed that the Gyaca knickpoint resulted from upstream migration of erosional waves starting from the Yarlung Tsangpo Grand Canyon (Fig. 1b) in response to uplift of the southeastern Tibetan Plateau prior to ∼10 Ma33. However, no geomorphic evidence (e.g., terraces) for downstream knickpoint retreat has been reported, although the decrease of base level at the Yarlung Tsangpo Grand Canyon might result in fluvial incision33. Our modeling results suggest that both headward incision of the Yarlung River since ~12 Ma due to uplift of the southeastern Tibetan Plateau33 or intensified monsoon precipitation40,41 are not compatible with the remarkable young AHe and AFT ages in the Gyaca gorge (Fig. 3e). In addition, rock outcrops in the Gyaca gorge and its downstream valley are mainly granitoids (Fig. 1c), thus lithological variation cannot be attributed as a major controlling factor for the formation of the Gyaca gorge. In contrast, as indicated by the well constrained Pecube model, a high-angle normal fault is required (Scenario B), to fit with an increased rock cooling rate at the Gyaca gorge since ~7 Ma (Fig. 3g). Therefore, we argue that movement of the Woca fault has controlled the formation of the Gyaca knickpoint since ~7 Ma, and that a similar mechanism might also explain the formation of other large knickpoints in the middle and upper reaches of the Yarlung River (Fig. 1b).

Based on QTQt and Pecube modeling results, a relatively high incision rate ( ~ 1.73 ± 0.07 km Ma−1) of the eastern Yarlung River initiated at ∼12 Ma, with a slightly decreasing trend in thermochronological ages from west to east along the Yarlung River. Published thermochronometric data from the externally drained portion of the eastern15,42,43 and central44,45 Lhasa terrane, as well as from large rivers in the southeastern Tibetan Plateau41,46, demonstrate that rapid exhumation rates ( > 1 km Ma−1) were pervasive across the southern and southeastern plateau between ∼17 and 10 Ma. These spatially large-scale synchronous rapid incision events most likely reflect enhanced Asian summer monsoon precipitation in the mid-Miocene40,41 that promoted the headward erosion of the Yarlung River channels (Fig. 4a). However, the incision rate upstream and downstream of the Gyaca gorge has decreased dramatically with the rapid activity of the Woca fault since 7 Ma (Fig. 2g and Supplementary Table 6).

Fig. 4: Schematic diagram illustrating the mechanisms for tectonic activity impeding the river incision.
figure 4

a Headwater incision of the river related to the uplift of the southeastern Tibetan Plateau along with intensified monsoon precipitation. b Regional coupled tectonic systems with accelerated extension of the plateau driving the formation of river knickpoints, which subsequently adjusts and reduces the regional river gradient, keeping river knickpoints relatively stable and slowing down upstream river incision.

Full size image

Onset of EHS exhumation and its link to rifting

Thermochronometric data from rifts in the Tibetan Plateau47,48 suggest that southern Tibetan Plateau has experienced rapid late Miocene to Pliocene rift acceleration34 (see Supplementary Fig. 7). This rapid rift activity with accelerated extension of southern Tibetan Plateau controls knickpoints, such as the Gyaca knickpoint, in the plateau interior. Further, the high rate of rifts extension in southern Tibetan Plateau also facilitate thinning of the upper crust while its lower crust is thickened by ongoing compression34,49,50. This could contribute to accelerating eastward crustal flow51 at the Eastern Himalayan Syntaxis (EHS) driving localized deformation and uplift52,53,54, leading to active coupling between crustal rock advection and river erosion55 since ~7 Ma (Fig. 4b). Detrital thermochronological data from foreland basin sediments downstream of the EHS also suggest coupling between tectonic uplift and erosion starting at 8 Ma55, 7-5 Ma56 and/or 6-4 Ma57.

Currently, the mechanism of rapid exhumation in the EHS remains controversial. A current paradigm is the Tectonic Aneurysm model, in which spatially focused surface erosion driven by the Yarlung River locally might accelerate rock uplift and exhumation of hot and weak crust at the syntaxes21,22,58. Evidence from sediments in the Himalaya foreland25, Bengal basins26 and our thermochronometric data (Fig. 2) show that the Yarlung River was definitely set in its course before ~12 Ma. However, this does not specifically address the role of the Yarlung River in driving initial rapid exhumation of the EHS, which is thought to have begun 8-6 Ma55,56,57. Although Wang et al. argued for a rapid uplift of the EHS at ∼2.5 Ma based on sediment fill immediately upstream of the Yarlung Zangbo gorge16,59,60, we consider this as only one of several stages of EHS uplift, rather than the initial uplift episode. The work by King et al. shows that the optically stimulated luminescence thermochronology to the northeast of the Namche Barwa supports the idea of northward migration of the high exhumation locus (as proposed by Seward and Burg61) instead of control by localized very fast fluvial incision62. Here, we find evidence for a synchronous exhumation pulse in the rifts (Supplementary Fig. 7) and EHS55 since the late Miocene. This rapid exhumation implies that the tectonic system resulting from the accelerated late Miocene extension of southern Tibetan Plateau drove regional fault activity to control rapid exhumation and the formation of stable knickpoints in the rifts and EHS.

Contribution of stable knickpoints to the stabilization of southern Tibetan Plateau

With the stabilization of knickpoints, the landscape equilibrates towards topographic steady-state equilibrium whereby for a given climatic condition, the slope of rivers adjust toward a geometry that allows an incision rate equal to rock uplift rate63. Therefore, the surface uplift rate might be zero, which would sustain a high-elevation base level for the Yarlung River in the rifts and EHS. Meanwhile, coeval tectonic evolution with the accelerated late Miocene extension of southern Tibetan Plateau contribute to the stabilization of southern Tibetan Plateau. If only local fault activity of the rift controlled the Gyaca knickpoint (Fig. 1b, c), and the localized deformation and uplift of EHS downstream were not particularly fast, so that EHS knickpoint was not fixed (Fig. 1b), it would migrate headwards, gradually towards the Gyaca gorge. This would cause the base level to fall in the region and the river would incise more deeply into the plateau interior.

With the stabilization of the knickpoints, river gradients above the knickpoints are commonly so low that fluvial incision may generally have difficulty in keeping up with increase in tectonic uplift rates. As a result, channel slope and stream power would decrease transiently above these Yarlung River knickpoints (Fig. 4b). Simultaneously, the reduced river gradient and sediment flux could promote upstream aggradation, burial of bedrock valley floors19,64,65, and valley widening (Fig. 1d, b), which could cause widespread backwater aggradation, forming broad valley trains occupied today by braided river systems (Fig. 1d and Supplementary Fig. 1b–f). Field survey also found large paleo-lakes66 upstream and valley-fill sediments16 downstream of Gyaca gorge (Supplementary Fig. 8). This braided river only incised when the stripping of deposited alluvial material exposed bedrock to processes of abrasion, weathering and plucking67. Therefore, the erosion rate of the river above the knickpoints will be reduced obviously due to the reduction of stream power since the late Miocene. This is supported by the predicted low erosion rates (0.18 ± 0.03 km Ma−1) (Supplementary Table 6 and Supplementary Figs. 4–6), which is in agreement with the rates (0.04–0.20 km Ma−1) determined from cosmogenic nuclide (10Be) data in the middle reaches of the Yarlung River68.

Hence, coeval tectonic systems (rifting upstream and uplift and exhumation of the EHS downstream) with the accelerated extension of southern Tibetan Plateau resulted not only in upstream reduction in river gradient and stream power, but also stabilized downstream base levels, i.e. from a stable knickpoint at the EHS. This would subsequently drive the onset of stable knickpoints formation in the Gyaca and EHS, thereby impeding a wave of upstream erosion and facilitating the stabilization of southern Tibetan Plateau topography since the late Miocene.

Processes controlling the stability/instability of knickpoints in active orogens

When uplift rates are spatially variable for a given climatic condition, two fundamentally distinctive types of knickpoints (i.e. transient and stationary knickpoints) might develop under different conditions3,4,69. Transient knickpoints, which may be produced by a change in the background regional uplift rate or a discrete base-level lowering event69. Whereas stationary knickpoints may reflect abrupt spatial changes in uplift/incision rate69, and/or a hard basement3. Based on our results (Fig. 4) and published literature69, we suggest that the transient migration and stability of river knickpoints might coexist during the evolution of orogenic belts. For example, many other active orogenic belts also first grow to a certain height, and then experience laterally outward propagation70. This expanding uplift generally increases the steepness of rivers, and promotes upstream erosion30,31,32. But this expansion is not a simple one-dimensional process, it might also be accompanied by diverse strong regional tectonic deformation of the upper crust17,71,72 (e.g., rifting, subduction, and strike-slip faulting) and lower crustal flow51. Lower crustal flow in turn will promote and further drive regional tectonic deformation of the upper crust, leading to the stability of knickpoints. An excellent example of channel adjustment to tectonic forcing is the Eastern-Western Himalaya Syntaxis52,53. However, clearly not all knickpoints associated with active faulting are stalled, such as the locations of the Hatay Graben in Turkey73, a normal fault system in the eastern California3 and the headwaters of the Yellow River in the Tibetan Plateau74, which might be mainly related to the dip and slip rate of faults73,74. The Woca normal fault in this study is dipping in the upstream direction (Fig. 1c), which might reduce the amount of base-level fall experienced by that part of the Tibetan Plateau which is drained by the Yarlung River. Thus, we propose that the activity of the local fault systems can limit the migration of these river knickpoints under specific conditions, such as those on the eastern Yarlung River, which are still located in the area of regional tectonic deformation (Fig. 5). This process should subsequently reduce fluvial incision in the upstream and prevent base-level fall downstream of the knickpoints. This might stabilize the topography of the Plateau and protect high-elevation, low-relief surfaces over geologic time in active orogenic belts.

Fig. 5: Spatial distribution of regional tectonic features and large knickpoints along major rivers in the Tibetan Plateau.
figure 5

a Longitudinal profiles of Indus River. b Mekong River. c Salween River. d Yangtze River. e Yellow River. f Syr Darya River. g Amu Darya River.

Full size image

Our findings suggest that tectonic activity does not always necessarily enhance fluvial incision, and that the diverse activity can also decrease these rates by the adjustment of fault systems in active orogenic belts. This finding also provides a new mechanism for explaining the universal stability of topography and high plateaus in some other orogenic belts.

Methods

Thermochronological data and thermal history

Apatite and zircon separates were concentrated using standard heavy liquid and magnetic separation techniques. Analyses were performed at the School of Geography, Earth and Atmospheric Sciences, University of Melbourne (see Supplementary Note 2, 3 for further details). The thermal histories based on these thermochronological data were determined using the QTQt software37 (details in Supplementary Note 4), and spatial and temporal cooling pattern and rates along the Yarlung River were calculated based on the model results.

3D thermo-kinematic modeling

To identify detailed effects of rift activity on Yarlung River incision, we used the thermokinematic finite-element Pecube75 code, which has been developed to interpret thermochronological data by solving a 3D crustal heat-transport equation. We constructed two scenarios in Pecube, i.e. without (Scenario A) and with the influences of the Woca fault (Scenario B) (detailed model parameters are described in Supplementary Note 5). In Scenario A, the Yarlung River was incised from a flat plateau surface to form the topography, ignoring the influence of the Woca normal fault in the plateau interior. This scenario mainly reflects the upstream migration of erosional waves of the Yarlung River (Supplementary Fig. 2). In Scenario B, we considered not only the incision of the Yarlung River, but also incorporated the effect of the Woca normal fault (see Supplementary Fig. 2).