A total of 15 unconsolidated modern fluvial samples were collected over the study area, among which 12 tributary samples were taken at the outlet of their respective upstream basins (Figure 1). These tributary basins cover an area of 12 to 22,788 km² (Figure 1 and Table 1). Three samples were taken from the mainstream of the Yarlung Zangbo River. All samples were collected from fresh fluvial sandbars away from fluvial terraces, landslides, and/or debris flows, aiming to obtain representative erosional product for upstream drainage basins.

Samples were processed for in situ ¹⁰Be and ²⁶Al separation at the cosmogenic nuclide sample preparation laboratory of the Institute of Geochemistry, Chinese Academy of Sciences at Guiyang. Samples were first sieved to isolate the 250-550 μm fraction for analysis. After sieving, quartz was pre-concentrated and purified through magnetic separation and sequential selective dissolutions in 5%-10% HF/HNO₃. Dissolution of quartz and extraction of the Be and Al fractions were carried out following procedures proposed by Kohl and Nishiizumi (1992). Prior to completing HF dissolution, ~250 μg of ⁹Be carrier were added to each sample. After ion exchange column chemistry and hydroxides precipitation, Be(OH)₂ and Al(OH)₃ were oxidized in furnace at 800 ℃ and then mixed with Nb (BeO:Nb=1:6 by weight) and Ag (Al₂O₃:Ag=1:2 by weight), respectively. The mixed samples were pressed into a sample holder and loaded into the ion source of accelerator mass spectrometer (AMS).

Both ¹⁰Be/⁹Be ratios and ²⁶Al/²⁷Al ratios were then determined using the 0.5 MV AMS facility at School of Earth System Science, Tianjin University (Dong et al., 2018). This AMS system, designed by National Electrostatic Corporation (NEC), is an extension of the compact ¹⁴C AMS system for ¹⁰Be and ²⁶Al. Repeated measurements of secondary standards indicted ~3% precision on samples with ¹⁰Be/⁹Be and ²⁶Al/²⁷Al ratios of 10^-13-10^-12. The detction limit for ¹⁰Be/⁹Be and ²⁶Al/²⁷Al was 3.8 × 10^-15 and 6.9 × 10^-15, respectively (Dong et al., 2019). Our AMS performance is similar to the similar system installed at GNS Science, New Zealand (Zondervan et al., 2015).

¹⁰Be/⁹Be ratios were calibrated against the NIST SRM4325 standard and ²⁶Al/²⁷Al ratios were calibrated against the Al01-4-1 and Al01-4-2 standards (Nishiizumi, 2004). Independent chemistry procedural blanks were prepared from the same carrier solution in order to run for the correction of the measured ¹⁰Be or ²⁶Al concentrations. In this study, we took 1.387 ± 0.012 Myr for¹⁰Be half-life (Chmeleff et al., 2010) and 0.708 ± 0.017 Myr for ²⁶Al half-life (Nishiizumi, 2004).

We used the online erosion rate calculator (formerly known as the CRONUS-Earth online erosion rate calculator Version 3.0; http://hess.ess.washington.edu/math/v3/v3_erosion_in.html; Balco et al., 2008) to calculate ¹⁰Be basin-scale erosion rates based on the time-independent scaling model of Lal (1991)/Stone (2000) (St scaling model), which is calibrated to a ¹⁰Be surface production rate of 4.01 ± 0.33 atoms g^-1 yr^-1 at sea level high latitude (SLHL) (Borchers et al., 2016). The calculator consists of several MATLAB functions regarding the calculation of neutronic and muonic nuclide production rates. For the effective attenuation length of neutron in rock, the calculator internally took the value of 160 g cm^-2. The muonic production rate of ¹⁰Be and ²⁶Al were computed based on the depth below the surface and on-site atmospheric pressure, using the MATLAB function ‘P_mu_total' from the calculator (Balco et al., 2008). The topographic shielding factor for each basin was computed after Li (2013). As recommended by Balco et al. (2008), we substituted the latitude and elevation at sample site with mean basin latitude and mean basin elevation. Mean basin latitude and mean basin elevation were calculated with ArcGIS by performing the 90 m resolution Shuttle Radar Topography Mission (SRTM) DEM. The rock density was assigned with 2.7 g cm^-3 for each sample.

Repeated measurements of standards indicate 1σ uncertainties of < 3% for samples with ¹⁰Be/⁹Be and ²⁶Al/²⁷Al ratios of 10^-13 in this study. Measured ¹⁰Be/⁹Be ratios vary from (1.67 ± 0.06) × 10^-13 to (2.25 ± 0.06) × 10^-12, whereas ²⁶Al/²⁷Al ratios vary from (3.89 ± 0.11) × 10^-13 to (4.01 ± 0.10) × 10^-12. ¹⁰Be and ²⁶Al concentrations show a 10-fold spread from (2.28 ± 0.07) × 10⁵ to (2.26 ± 0.07) × 10⁶ atoms g^-1 and from (1.44 ± 0.05) × 10⁶ to (2.03 ± 0.05) × 10⁷ atoms g^-1, respectively (Table 2). Correspondingly, ²⁶Al/¹⁰Be ratios range from 6.33 ± 0.29 to 7.67 ± 0.68 (exceptionally high ratio for sample T1) (Table 2). In Figure 2, it can be seen that all samples, except sample T1, fall within the constant exposure-steady-state erosion island, satisfying the precondition of steady-state erosion. Note that sample T1 yields the highest ²⁶Al/¹⁰Be ratio (8.96 ± 0.37), which is extremely higher than the canonical spallation production ratio of 6.75 at surface. The high ²⁶Al/¹⁰Be ratio may have resulted from incorporation of deeply sourced materials or suggest overestimation of ²⁷Al measurement. Since sample T1 yields the highest native ²⁷Al concentration (287.60 ppm), it seems that overestimation of ²⁷Al measurement should be accountable for the high ²⁶Al/¹⁰Be ratio of T1. In any case, ¹⁰Be concentration and thus resulted erosion rate of sample T1 are still reliable given the high quality of determination of ¹⁰Be/⁹Be ratios (discussed in Section 5.1).

The ¹⁰Be concentrations give the basin-scale erosion rates which vary greatly from 20.60 ± 1.79 to 154.00 ± 13.60 m Myr^-1 (Table 2 and Figure 2). Spatially, samples with high erosion rates are mainly from the upper Nyang River basin. The highest value (154.00 ± 13.60 m Myr^-1) was derived from basin T8, which is situated at the source of the Nyang River (Figure 1). Erosion rates in the Nyang River basin are apparently higher than those in the Lhasa River basin. Relatively low erosion rates (20.60 ± 1.79 to 83.90 ± 7.14 m Myr^-1) were observed in drainage basins sampled from the upstream areas of the Gyaca knickpoint (samples T1-T7, M1-M3).

The χ° values decrease from upstream basins to downstream basins, inversely consistent with the increase in erosion rates (Figure 3a). In the downstream direction, the hypsometric curves shift from S- and concave-shaped to convex-shaped (Figure 4). Accordingly, the basin mean HI values generally increase downstream (Table 1). It can be deduced that stages of landscape in the downstream areas of the Gyaca knickpoint are younger than that in the upstream areas of the Gyaca knickpoint. Compared to the areas in the southern Himalayas, our 15 basins and the basins in the central TP generally yield higher HI values, indicating younger stages of landscape (Figures 5d-9d and 11). In addition, basin-averaged hillslope angles and local reliefs for the 15 basins are also generally lower than those in the southern Himalayan areas. As shown in Figures 5-6 and 8-9, erosion rates roughly show positive correlations with variance in hillslope angles, local reliefs and precipitation. In Figure 7, however, these correlations are unable to be tested due to the high risk of perturbation by lanslides (Abrahami et al., 2016). The correlation between relief and erosion rate is not so significant among the 15 basins (Figure 3b). As indicated in Figure 3c, precipitation exerts a positive influence on erosion for most of the sampled basins. The mean EI value varies from 1.0-1.4 for the 15 basins (Table 1) and shows no systematic correlations with ¹⁰Be-derived basin-scale erosion rates (Figure 3d). Vegetation coverage (measured by NDVI) seems to exert different controls on erosion rates in the upper and middle drainage basins (Figure 3e). In the upper and middle drainage basins, the increased vegetation density accelarates surface erosion processes. However, landscapes with sparse plants in the middle drainage basin seem to be prone to erosion.

Before the establishment of the cosmogenic nuclide method, traditional methods of estimating basin-scale erosion rates depended on the measurement of modern sediment flux. However, it is difficult to measure sediment flux accurately, because this method depends on centennial river gauging records which cannot normally document rare flood events transporting large volumes of sediment. Moreover, significant sediment storage or remobilization in upstream valleys can strongly alter the precision of measuring modern sediment flux. In these cases, recorded sediment flux cannot be directly translated into secular erosion rates. Basically, ¹⁰Be-derived erosion rates can integrate over a period of 10³-10⁵ yr (Granger et al., 1996). By comparing ¹⁰Be derived basin-scale erosion rate with present-day hydrologic erosion rate, we can assess the fluctuations in erosion rate and provide information about anthropogenic impacts on basin-scale erosion processes.

Previous studies reported similar results between ¹⁰Be-derived basin-scale erosion rates and present-day hydrologic erosion rates even though they represent different time-scales (Von Blanckenburg, 2005). There are four major hydrologic stations along the Yarlung Zangbo River (Figure 1). The Lhaze hydrologic station and the Nugesha hydrologic station are close to the mainstream samples M1 and M3. The Yangcun station records sediment flux from upstream basins above the Gyaca knickpoint. The Nuxia station is located a little far from downstream at the confluence between the Nyang River and the Yarlung Zangbo River. It records sediment fluxes from the upstream and midstream of the Yarlung Zangbo River. Basic information about the four stations is shown in Table 3.

Based on the annual sedement load recorded by these four hydrologic stations (Shi et al., 2018), we calculated the present-day hydrologic erosion rates for the above drainage areas. It can be seen that the present-day hydrologic erosion rates increase downstream for areas above the Gyaca knickpoint (Table 3). Similarly, ¹⁰Be derived basin-scale erosion rates also show a downstream increase in this region (Figure 3 and Table 2). Note that sample M3 yielded a lower erosion rate than the upstream Nugesha station (61.3 m Myr^-1 and 80.4 m Myr^-1, respectively). The relatively high erosion rate for the Nugesha hydrologic station may have resulted from anthropogenic influences such as agriculture in upstream areas. Several studies have reported the anthropogenic impacts on basin-scale erosion processes. For instance, Hewawasam et al. (2003) evaluated the influences of human agricultural activities on soil production and bedrock weathering processes. They found that present-day soil erosion rates have been strongly accelerated in agriculturally utilized regions. Anthropogenic disturbances have become a severe problem and could even alter the cosmogenic nuclide concentrations from large river basins. Human tillage activities can sometimes be deep enough into the regolith. Schmidt et al. (2016) demonstrated that these activities could excavate materials with low ¹⁰Be concentrations from great depth to the surface. Addition of these materials to streams will consequently raise the apparent ¹⁰Be basin-scale erosion rates. However, this may not be the case for the drainage area above the Nugesha hydrologic station since sample M3 yielded a ²⁶Al/¹⁰Be ratio of 6.74 ± 0.60 and indicates steady-state erosion.

The Yangcun station and the Nuxia station are located within the narrow gorges where knickpoints developed (Wang et al., 2015). Before these knickpoints, Wang et al. (2015) reported a huge amount of sediment storage. The decrease in present-day hydrologic erosion rate (from 122.1 m Myr^-1 to 57.6 m Myr^-1) may result from sediment storage within the broad valley sections between these two stations. As discussed in many previous studies, sediments that have been stored for prolonged periods in upstream valleys could influence the downstream erosion rate derived using cosmogenic nuclide when these buried deposits are tapped episodically by rare flood events or by river network segmentation triggered by tectonic activity. Given the high frequency and high magnitude of flood events triggered by breaks of glacial or landslide dammed lakes (Lang et al., 2013), the formerly buried sediments between these two stations are very likely to be remobilized and channelled downstream. In the downstream of the Yarlung Zangbo River, the risk of overestimating ¹⁰Be basin-scale erosion rate would increase due to incorporation of low-¹⁰Be -concentration sediments.

In the use of ¹⁰Be in deriving basin-scale erosion rate, a steady-state erosion assumption for upstream basins is strictly required (Granger et al., 1996). However, during extreme erosion events such as landsliding, debris flow and rock avalanches, large numbers of materials with lower ¹⁰Be concentration might mix into river sediment budget and cause significant biases in derived erosion rates (West et al., 2014). These catastrophic erosion processes are common across our study area due to poor soil and water conservation capacity. To evaluate the influence of these non-steady-state erosion processes, we conducted a paired analysis of ¹⁰Be and ²⁶Al. As shown in Figure 2, except for sample T8 which failed in ²⁶Al measurement, 13 out of the remaining 14 samples were plotted within the constant exposure—steady state erosion island, satisfying the assumption of steady-state erosion. However, sample T1 has a greatly high ²⁶Al/¹⁰Be ratio of 8.96 ± 0.37 and was plotted in the “forbidden area” (area above the steady state erosion island) (Lal, 1991). By analysing the relationship between ²⁶Al/¹⁰Be, depth, and rock density, Akçar et al. (2017) demonstrated that the ²⁶Al/¹⁰Be ratio would increase with depth in bedrock profile (from 6.75-7.25 in the upper 2 m to ~8.4 between 5-10 m at depth) due to discrepancy in muonic production for ¹⁰Be and ²⁶Al. Materials originating from a greater depth are very likely to be exhumed to the surface in landscapes vulnerable to deep erosion, thus increasing the ²⁶Al/¹⁰Be ratios in sampled river sediment. Addition of deeply sourced materials can not only increase the ²⁶Al/¹⁰Be ratio in river sediment but will also lower their ¹⁰Be concentration due to low nuclide production rate by muons. However, sample T1 yields the highest ¹⁰Be concentration over the other 14 samples, implying no dilution of ¹⁰Be by deep-sourced materials. Alternatively, overestimation of ²⁷Al can also lead to a high ²⁶Al/¹⁰Be ratio. Noting that, sample T1 has the highest native ²⁷Al concentration. In this case, overestimation of ²⁷Al seems to be accountable for the high ²⁶Al/¹⁰Be ratio of sample T1. Given that the adjacent basins (T2 and T3) that are undergoing steady-state erosion have similar topographic characteristics to basin T1 (Figure 2 and Table 1), it is reasonable to deduce that basin T1 satisfies the assumption of steady-state erosion. Nevertheless, ¹⁰Be concentration and thus, derived erosion rate of basin T1 are still reliable.

The addition of ²⁶Al can suggest whether deep-sourced sediments have mixed in. However, we didn't conduct extensive sampling to cover all the basins located in the two flanks of the Yarlung Zangbo River. Basins that are not undergoing steady-state erosion may actively supply sediments to the mainstream. Thus, it may be erroneous to constrain their erosion rate with just ¹⁰Be and the derived erosion rate is very likely to be overestimated. At least, it can still be deduced that contributions of these basins are less important to dominate the mainstream sediment fluxes since the ²⁶Al/¹⁰Be ratios of samples M1, M2 and M3 have recovered close to 6.75 within uncertainty (Figure 2). These mainstream samples might have received significant amounts of sediment from other areas undergoing steady-state erosion.

Topography and climate are two interactional factors controlling the erosion rate. They may alternately exert leading control on erosion in different regions. Abrahami et al. (2016) reported that in the Sikkim Himalayas, erosion processes are strongly controlled by glacial inheritance while distribution of precipitation plays only a second-order role. In central Nepal, variations in precipitation also exert only a second-order control in modulating the erosion signal (Godard et al., 2014). Instead, the primary force that shows a significant co-variation with erosion is the long-term rock uplift rate. Further east in the Bhutanese Himalayas, rapid erosion in steep river basins along the Himalayan front is thought to be driven by high rainfall associated with the Indian summer monsoon. In these regions, climate may play the decisive role on erosion processes (Hodges et al., 2004).

Based on our analysis, annual precipitation and vegetation coverage (measured by NDVI) vary largely across the study area (Figures 3c and 3e). Erosion rates are spatially well coincident with an increase in precipitation and NDVI for basins T1-T7 and M1-M3 in the upstream areas of the Gyaca knickpoint, emphasizing the control of climate in this region. With respect to vegetation coverage, distribution of vegetation across different lithologies may lead to different geomorphologic processes, and as a result, different landscapes form. Previous cosmogenic nuclide studies have shown that dense vegetation coverage will accelarate surface erosion processes for most landscapes (Hahm et al., 2014). Increase in vegetation coverage may result in intense variations in soil conservation and rock weathering through different rooting and organic ligand release. As a result, it will result in differences in erosion by chemical, physical, and biological processes. This may be the case for the upstream areas of the Gyaca knickpoint (basins T1-T7 and M1-M3). In this region, increase of vegetation density accelarates surface erosion processes (Figure 3e). However, vegetation exerts a negligible role in the downstream areas of the Gyaca knickpoint since vegetation shows no positive influence on erosion (basins T8-T12; Figure 3e).

Our results revealed that basin-scale erosion rates generally fit well with indices that deal with the terrain gradient of sampled basins, i.e., mean basin χ° values (Figure 3a) and mean hypsometric integral (HI) (Figures 5d-9d), but almost no significant correlations with local relief (Figure 3b). Since the HI values are positively correlated with tectonic uplift rate, basins with relative uplift and basins with relative subsidence can be clearly discriminated with the help of the scatter plot of HI values in relation to average slope (Lifton and Chase, 1992). As shown in Figure 10, drainage basins (T1-T7, M1-M3) were associated with low average slope values (average=17.70°) and low HI values (average=0.41), indicating relatively low values of uplift rate and mature or final (old) erosional stages (group HI-1). These drainage basins are positioned in the upstream areas of the Gyaca knickpoint. Drainage areas (T8-T12) located in the Nyang River basin have high average slope values (average=26.20°) and intermediate HI values (average=0.51), representing relatively young transient landscapes with high uplift rates (group HI-2). We demonstrated that the different erosion rates for these two regions can be interpreted by their different tectonic uplift background.

The χ° values of the 15 basins fit well with derived basin-scale erosion rates, with lower χ° value indicating higher erosion rates (Figure 3a). The 10 basins (basins T1-T7, M1-M3) located in the upstream areas of the Gyaca knickpoint have relatively higher mean χ° values than those located in the downstream areas (Figures 3a and 11a). The apparent disequilibrium of the river channel geometry might have resulted from tectonics or changing climate. As discussed above, drainage basins (T1-T7, M1-M3) with low average slope values and low HI values imply low uplift rates, whereas basins (T8-T12) located in the downstream areas of the Gyaca knickpoint have high average slope values and intermediate HI values, representing a relatively younger more transient landscape with high uplift rates. We demonstrated that tectonics should be responsible for the difference in erosion rates between these two regions.

When compiling all the reported erosion data shown in Figure 1 using cosmogenic nuclide in the southern Tibet, we found that both topography and climate play significant roles in erosion in southern Tibet. In Figure 5, it can be seen that erosion rates show roughly positive correlations with variance in hillslope angles, local reliefs and precipitation, indicating that both topography and climate play significant roles in erosion in this region. In other regions shown in Figures 6 and 8-9, similar conclusions can also be drawn. In this situation, it is hard to tell whether topography or climate is the main driving force on erosion across southern Tibet. However, at least, we can make a conclusion that topography plays a leading role in erosion over climate in our study area. Erosion rates are modulated by precipitation and/or vegetation (second-order control).

Moreover, erosion rates are also proved to be strongly related to rock erodibility in mountainous areas. Landscapes with high rock strength could evolve steep topography. In our study area, the averaged rock erodibility lies between 1.0-1.4 and indicates low to medium erodibility (Table 1). The negligible correlation between rock erodibility and erosion rate revealed that lithological complexity has no significant influence on erosion rate for the 15 basins in this study (Figure 3d).

By far, hundreds of local or basin-scale erosion rates derived using cosmogenic ¹⁰Be have been reported on the Tibetan Plateau (Godard et al., 2014; Abrahami et al., 2016; Lupker et al., 2017). All of this data together can depict a general erosion picture of the plateau. It is worth mentioning that biases may exist when sampling rock from surfaces that are more resistant to erosion than the surrounding landscapes. However, erosion rates from different lithologies across the plateau yielded similar and low results, indicating low rates in secular erosion-planation on the central plateau surfaces (Li et al., 2014). Our results are similar to those from adjacent basins repotred by Lupker et al. (2017) but much lower than those from the marginal Himalayas (Figures 5-9), implying that the plateau is being dissected at its boundaries by river headward erosion.

The χ° value highlights that the disequilibrium of the river channel geometry resulted from tectonics or changing climate. Watersheds with lower χ° values (higher erosion rates ) tend to migrate towards areas with higher χ° values (lower erosion rates) through river headward capture. As shown, drainage river systems in the southern Himalayas exhibit much lower χ° values thus indicating higher erosion rates (Figure 11a). Towards the north and into the central plateau, the χ° values generally increase with increasing elevations and decreasing local reliefs (Figures 11b-11l). The northward increase in river χ° values revealed that drainage divides will migrate northward in the long-run under these circumstances. This trend is consistent with the spacial pattern of erosion that landscapes in the southern Himalayas are northwardly eroding faster than areas in this study and areas in the central plateau. It can be predicted that generally northward erosion in the southern Himalayas and the YZRB will ultimately consume the stable central TP surfaces.