To constrain the exhumation of the Three Rivers Region, especially the section of the knickzone, sampling from a vertical transect was performed from the western margin of the central Three Rivers Region (Figures 1B, 3). Six rock samples were collected from Mesozoic granitic intrusions from the near peak of the Heipushan to the deeply incised valley bottom of the Dulong River (Kongdang Village) (Figure 3). Sample’s elevations range from 3,326 to 1,562 m, forming a vertical profile spanning ∼1,760 m relief over a lateral extent of ∼18 km (Figure 3B). The intrusions, where the samples were collected, are undeformed Mesozoic plutons with intrusive contact, in which no faulting has been observed during field investigations. Previous AFT (closure temperature, ∼110 ± 20°C; Reiners and Brandon, 2006) ages reported by Lei et al. (2006) for the same transect (Figure 3) are between 4 and 6.8 Ma (Figure 4). To gain more detailed information for the cooling history since the Pliocene and new insights into surface processes, we report new AHe (closure temperature, ∼60 ± 20°C; Farley et al., 1996) data for the Dulong batholith.

Apatite (U-Th)/He analyses for the Dulong transect were conducted at the National Institute of Natural Hazards, Ministry of Emergency Management of China (NINH-MEMC). Apatite concentrates were extracted using standard crushing, sieving, electromagnetic, and heavy liquid mineral separation techniques. Apatite grains with euhedral morphology and no visible inclusions were selected under a microscope and only grains >70 μm in both length and width were considered suitable for (U-Th)/He dating. Grain dimensions were measured from digital photographs for the calculation of the equivalent spherical grain radius and the α-ejection correction factor. Each grain was then wrapped in a 1 mm × 1 mm platinum capsule and loaded into the laser chamber. Each grain was thermally outgassed under vacuum at ∼900°C for 5 min, using a diode laser (970 nm wavelength) with 8 A current. Then, spiked with ³He, gas volumes were determined using a PrismaPLus QME 220 quadrupole mass analyzer at NINH-MEMC. We checked that gas released during replicate heating yielded approximately the same as hot blanks to ensure total extraction for each grain. After degassing, molar abundances of U and Th were determined by isotope dilution using a mixed ²³⁵U-²³⁰Th spike. U-Th analyses were carried out on an inductively coupled plasma quadrupole massspectrometer at NINH-MEMC. The age calculation was processed by applying the α-ejection correction factor (F_T) (Farley et al., 1996) to each crystal to derive a corrected (U-Th)/He age (Table 1). The age error was derived from the analytical uncertainties in U and Th measurements, and the variance of the single grain ages. Six fragments of Durango apatite were run as reference standards together with and identically to our samples to verify analytical accuracy. A weighted mean average age of 31.7 ± 0.5 Ma (Table 1) was obtained for these fragments, which is in consistent with the nominal age of the Durango apatite (McDowell et al., 2005).

To investigate the thermal evolution of the Dulong vertical transect, we modeled the thermal history using the program QTQt, which has been developed to invert thermochronological ages for multiple samples with a known altitudinal relationship implementing a Markov chain Monte Carlo method (Gallagher, 2012) Figure 5. The modeling approach employs an alpha-damage-dependent kinetic model of helium diffusion in apatite (Flowers et al., 2009) and a multikinetic AFT annealing model (Ketcham et al., 2007). The AHe data in this study and AFT data from Lei et al. (2006) are modeled jointly. The input parameters used to model the thermal history for individual samples are as follows: (1) present-day mean surface temperature of 10 ± 10°C; (2) the prior for the paleotemperature offsets, or temperature difference between the uppermost and lowermost samples in a vertical profile, were defined as ∼53 ± 53°C equivalent to temperature gradient prior of 30 ± 30°C/km (Clark et al., 2005) and the temperature offsets were also allowed to vary over time; (3) an initial time-temperature constraint is set at 100–200°C at a time span slightly older than the oldest AFT age. These prior settings were always included with a large uncertainty so as to give the modeling enough freedom to search for a wide range of data-constrained thermal histories. The final thermal history models were sampled 400,000 iterations: 200,000 used to stabilize or burn-in the inversion, and the second 200,000 used to form the posterior ensemble (Gallagher, 2012). Exploratory runs using larger numbers did not appreciably change model outcomes.

Bedrock river profiles are often described using the stream power incision model (Whipple and Tucker, 1999), which expresses the erosion rate in terms of channel slope and drainage area ∂ z ( x , t ) ∂ t = U ( x , t ) − K ( x ) A ( x , t ) m S n (1) where ( ∂ z ( x , t ) ∂ t ) is the change in elevation of the channel bed with respect to time, U is rock uplift rate relative to the base level, K is rock erodibility, A is drainage area, S is channel slope, m and n are constants. Under the assumption of a topographic steady state ( ∂ z ( x , t ) ∂ t = 0 ) and U and K are spatially and temporally uniform, the equilibrium slope is then a function of S = ( U K ) 1 n A − m n (2) where m n is the concavity of the equilibrium profile and ( U K ) 1 n is the channel steepness which can be determined by scaling the slope and area relationship.

Deriving the channel slope data directly from the digital elevation model (DEM) can be problematic due to the noise of the DEM data. To avoid the scatter of noise during the estimation of slope, we used an alternative method (Perron and Royden, 2013) for the equilibrium river profiles by substituting the channel slope with elevation, which leads to z ( x ) = z ( x b ) + ( U KA 0 m ) 1 n χ (3) and χ = ∫ x b x ( A 0 A ( x′ ) ) m n dx′ (4) where x b is the reference of local base level, and A 0 is an arbitrary scaling factor. Then channel steepness K sn is the slope of the χ -elevation plot K sn = ( U KA 0 m ) 1 n (5) which is proportional to the rock uplift rate.

We used the SRTM DEM, which has a resolution of ∼90 m, to extract the longitudinal profiles and steepness index of the Dulong and Salween rivers (Figure 6). A threshold drainage area of 5 km² was used to exclude regions that are potentially dominated by debris flows or hillslope processes. We selected a concavity, m / n , of 0.45 and a scaling area, A 0 , of 1 m² (Wobus et al., 2006). The channel steepness was then estimated from the slope of the χ -plot with the linear regression method by using a χ interval of 1.

Four to eight single-grain AHe age analyses were performed for each of the six Dulong samples, as summarized in Table 1. The samples yield mostly consistent AHe ages except the uppermost sample (G18-1) has two abnormally old ages. The two AHe outliers of sample G18-1 do not show clear relationships with eU and grain size (Table 1), indicating radiation damage and grain size variation do not appear to be controlling the distribution of ages (Gautheron et al., 2012). U-zoning in the core leads to overestimate of the alpha-ejection correction, but cannot explain the abnormally old ages in our study, because even the uncorrected ages (11.5 and 13.4 Ma) of the two grains are older than the AFT age (6.8 ± 0.5 Ma) at the same elevation. Additional sources of ⁴He other than the analyzed apatite, such as U-rich mineral inclusions in apatite, U-rich neighbouring minerals (Spiegel et al., 2009) may be possible explanations for the outliners. Excluding outliers, all remaining AHe data show a strong positive relationship with elevation (Figure 4). The age-elevation relationship has an inflection point at the elevation of ∼2,500 m, and the AHe ages below this point are generally less than 3 Ma, while the AHe age above are significantly older (3–7 Ma). Excluding outliers, the calculated weighted mean AHe ages range from 6.18 ± 0.9 to 2.36 ± 0.43 Ma and show a positive correlation with elevation. The regression of the age-elevation relationship suggests a significant increase in erosion rate from ∼0.18–0.3 km/Myr to ∼1.3–3.0 km/Myr after ∼2.6 Ma (see below for the timing from the thermal history modeling).

The modeling results show a thermal history with two phases of rapid cooling since the late Miocene (Figure 5A). The first episode commenced at ∼7–8 Ma; all the samples passed through the AFT partial annealing zone (PAZ) rapidly and some upper samples might have reached the AHe partial retention zone (PRZ) during this cooling event. This phase of fast cooling also revealed by the overlap of the AHe and AFT ages (∼6–7 Ma) in the uppermost elevation (Figure 4). However, the current available data cannot provide a precise constraint on the timing of onset. The duration of this episode of rapid cooling, the induced mechanism and the potential links to tectonics or climate change need further work that are well beyond the scope of this study. After the first phase of fast cooling, a period of slow cooling or isothermal holding lasted for ∼5 Myr; then the cooling rate increased at ∼2.6 Ma, and all the samples exhumed to the near Earth’s surface (Figure 5A). Such a thermal history is generally consistent with our thermochronological observations (Figure 5B). Note that the AHe observations are very well fitted by the modeled values, supporting the validity of the Quaternary enhanced cooling and exhumation (Figure 5A). In summary, the inverse modeling results from the vertical transect suggest that it experienced two episodes of rapid cooling commenced before Pliocene and at the beginning of the Quaternary, which are in accordance with the age-elevation profile (Figure 4). In the sections below, we focus on the Quaternary enhanced cooling and expand its implications to regional exhumation and tectonics.

Our new AHe data and thermal modeling suggest increased exhumation rates in the upper reach of the Dulong River, central Three Rivers Region, at the beginning of the Quaternary (∼2.6 Ma) (Figure 5). Although our data cannot provide detailed information for the exhumation processes after 2.4 Ma (Figure 4), the mean exhumation rate of ∼0.83 mm/year since ∼2.4 Ma can be estimated given the ∼2 km magnitude of erosion derived from the closure temperature of AHe (∼60 ± 20°C; Farley et al., 1996) and the recommended geothermal gradient (∼30°C/km; Clark et al., 2005) in the region. Thus, we conclude that the study area should have experienced faster exhumation during the Quaternary than before (Figures 4, 5). This is similar with previous findings from thermochronological studies at about the same latitude in the gorges of the Salween and Mekong (Figure 1B). Pre-existing thermochronological data from the valley bottoms of the Salween and Mekong have suggested enhanced exhumation (>0.75 mm/year) near 28°N in the past 2 Myr (Yang et al., 2016). A recently reported set of AHe and AFT data from Kawagebo massif have also revealed rapid Quaternary exhumation (>1 mm/year) at the valley bottom of the Mekong (Replumaz et al., 2020). Our results suggest that this increase in exhumation rate has also occurred in the upper Dulong River, the western margin of the central Three Rivers Region. Together with previous studies, we infer that an enhanced Quaternary exhumation with significant magnitude may exist in the central Three Rivers Region. This conclusion is supported by the increase in sedimentary flux to the marginal sea basins in the past 2 Myr (Métivier et al., 1999; Clift, 2006).

It is worth noting that all the young thermochronological ages younger than 2.6 Ma in the Three Rivers Region are located between 26 and 30°N (Figure 6A), implying that the central part of the Three Rivers Region may have experienced fast erosion during the Quaternary. The locus of rapid erosion was focused at the same area in different river gorges may suggest that the same mechanism may underline this phase of fast exhumation in the central Three Rivers Region. As mentioned above, the Quaternary enhanced exhumation could be induced by tectonic uplift, climate change, river reorganization or fault activity. Based on several lines of evidence, the rapid Quaternary exhumation in the central Three Rivers Region was most likely controlled by localized tectonic uplift. First, the locus of rapid erosion coincides with the conspicuous large-scale knickzone in the Three Rivers Region (Figure 6). In this region, the Three Rivers and the Dulong River are most closely spaced, have the highest steepness index in river long profiles, coinciding with the steepest reach in plateau edge as suggested by the maximum elevation envelop (Figures 6B,C). The pattern of the knickzones, with high steepness values limited to the knickzone region and similar lower values above and below the knickzones (Figures 6B,C), identifies they as “vertical-step” knickpoints (Kirby and Whipple, 2012), suggesting that they are related to spatially focused rock uplift given that there is no obvious variation in lithology associated with the knickzones (Replumaz et al., 2020). Second, in the central Three Rivers Region, the low-relief and high-elevation landscapes are absent (Clark et al., 2006) due to the intense dissection and high relief, which may be caused by the local uplift. Third, short-term (millennial) erosion pattern in the Three Rivers Region revealed by detrital cosmogenic nuclide was used to infer that tectonics is the primary control and the east-west enhanced erosion gradient mirrors a gradient in rock uplift rates (Henck et al., 2011). Fourth, structural and kinematic analyses reveal that the amalgamation area of the Gaoligong and Chongshan shear zones, from Fugong to Gongshan area, is just located at the neck of the large-scale boudin structure and experienced strongly partitioned dextral transpression and consequent uplift at the corner of the eastern Himalayan syntaxis (Huang et al., 2015). Finally, in nearby region, enhanced rock uplift since ca. 2.5–2 Ma in the eastern Himalayan syntaxis has been inferred from the existence of Quaternary thick alluvium sediments above the Yarlung Tsangpo gorge (Wang et al., 2014) and multiple thermochronometries and geomorphology analysis (King et al., 2016; Salvi et al., 2017; Yang et al., 2021). It is likely that the Quaternary enhanced exhumation in the Three Rivers Region is synchronous with the eastern Himalayan syntaxis, and a response to the continuous indentation of the northeast corner of the Indian plate.

It is significant that the ages we obtained for the initiation of rapid exhumation in the central Three Rivers Region (∼2.6 Ma) closely approximate the estimated timing of global cooling (Herman et al., 2013). U-shaped valleys are widely distributed above ∼3,000 m in the Dulong area (Figure 2C), indicating the imprint of glacial erosion. The past extent of glaciers in the southeast Tibet, reconstructed based on glacial landforms and sediments, indicates that this was one of the most extensively glaciated area of the Tibetan Plateau during the Quaternary (Li, 1996; Fu et al., 2013). However, only a few areas exhibit rapid Quaternary exhumation implied by low-temperature thermochronology in the vast region of the southeast Tibet (Figure 1B), suggesting the glacial erosion was not the main force for the fast exhumation in the Three Rivers Region.

The Three Rivers Region is strongly influenced by the Asian monsoon precipitation (An et al., 2001) (Figure 2A). The youngest ages in the Three Rivers Region are in the area where the precipitation rate decreases abruptly (Figure 6B,C and 2A). Along the Salween, the modern rainfall increases steadily from the immediately south of the edge-plateau to the lowland while erosion rates decrease (Figure 6C). Thus, if during the Quaternary the climatic gradient was similar to the modern one, the exhumation pattern is unlikely related to the precipitation.

Drainage area loss or gain will decrease or increase the erosion rate near the capture point based on the stream power law (Whipple and Tucker, 1999). Potential capture of the formerly northwest-to-southeast-flowing paleo-Yarlung Tsangpo-Dulong River by the Brahmaputra River was proposed as the drainage reorganization event in the region (Clark et al., 2004), although the timing of this process is still unclear. However, if this capture event took place in the headwater of the Dulong River, the loss of the drainage area would result in the decreased erosion rate in the downstream of the capture point. This case is not supported by the observation of our study. Our results cannot preclude the possibility of the rapid exhumation induced by river capture in the downstream of the knickzone of the Dulong River, but we argue that even though the river capture occurred in the Quaternary and resulted in the consequent rapid exhumation, it was possible triggered by the enhanced rock uplift.

In summary, although climatic factors or river capture may play somewhat roles on the enhanced Quaternary exhumation in the central Three Rivers Region, the localized tectonic uplift may have exerted first-order control on this exhumation, similar to the eastern Himalayan syntaxis, the central Longmen Shan and the Gongga Shan where the tectonics activity was regarded as active during the recent past.

The geodynamics of the formation of the southeastern Tibetan Plateau is hotly debated. Various models have been proposed to explain the plateau growth and the formation of the unique landscape in this region. They include: indentation and progressive crustal thickening (England and McKenzie, 1982), tectonic extrusion (Tapponnier et al., 2001), lower crustal channel flow (Clark and Royden, 2000) or whole crustal flow (Copley and McKenzie, 2007) driven by the topographic difference between the plateau and its surroundings. Although the timing of each model exerted is still controversy, all existing models have in common that the southeastern Tibetan Plateau must have grown outwards with respect to its interior. This outward expansion of the plateau is also thought to be responsible for the propagation of topography and thus the focus of erosion. However, our new results and the available datasets indicate that the locus of rapid erosion in the recent geological past was confined to the central part of Three Rivers Region rather than the plateau margin (Figure 6A), in contrast to previous plateau expansion models. Our results cannot provide constrains on the topographic evolution or surface uplift during the Quaternary. Nevertheless, significant regional rock uplift in the high-strain zone probably caused by the expansion of the eastern Himalayan syntaxis is required to explain the previous and our new thermochronological data in the Three Rivers Region.