Special Issue: Climate Change and Its Regional Response

Erosions on the southern Tibetan Plateau: Evidence from in-situ cosmogenic nuclides 10Be and 26Al in fluvial sediments

  • ZHANG Xiaolong , 1 ,
  • XU Sheng , 1, * ,
  • CUI Lifeng 1 ,
  • ZHANG Maoliang 1 ,
  • ZHAO Zhiqi 2 ,
  • LIU Congqiang 1
  • 1. Institute of Surface-Earth System Science, School of Earth System Science, Tianjin University, Tianjin 300072, China
  • 2. School of Earth Science and Resources, Chang'An University, Xi'an 710054, China
*Xu Sheng (1963-), Professor, specialized in isotopic geochemistry and geochronology. E-mail:

Zhang Xiaolong (1989-), PhD Candidate, specialized in cosmogenic nuclide applications. E-mail:

Received date: 2021-02-02

  Accepted date: 2021-09-14

  Online published: 2022-04-25

Supported by

Second Tibetan Plateau Scientific Expedition and Research (STEP) Program(2019QZKK0707)

National Key Research and Development Program of China(2020YFA0607700)

National Natural Science Foundation of China(41930863)

China Seismic Experimental Site(2019CSES0104)


Investigating topographic and climatic controls on erosion at variable spatial and temporal scales is essential to our understanding of the topographic evolution of the orogen. In this work, we quantified millennial-scale erosion rates deduced from cosmogenic 10Be and 26Al concentrations in 15 fluvial sediments from the mainstream and major tributaries of the Yarlung Zangbo River draining the southern Tibetan Plateau (TP). The measured ratios of 26Al/10Be range from 6.33 ± 0.29 to 8.96 ± 0.37, suggesting steady-state erosion processes. The resulted erosion rates vary from 20.60 ± 1.79 to 154.00 ± 13.60 m Myr-1, being spatially low in the upstream areas of the Gyaca knickpoint and high in the downstream areas. By examining the relationships between the erosion rate and topographic or climatic indices, we found that both topography and climate play significant roles in the erosion process for basins in the upstream areas of the Gyaca knickpoint. However, topography dominantly controls the erosion processes in the downstream areas of the Gyaca knickpoint, whereas variations in precipitation have only a second-order control. The marginal Himalayas and the Yarlung Zangbo River Basin (YZRB) yielded significantly higher erosion rates than the central plateau, which indicated that the landscape of the central plateau surface is remarkably stable and is being intensively consumed at its boundaries through river headward erosion. In addition, our 10Be erosion rates are comparable to present-day hydrologic erosion rates in most cases, suggesting either weak human activities or long-term steady-state erosion in this area.

Cite this article

ZHANG Xiaolong , XU Sheng , CUI Lifeng , ZHANG Maoliang , ZHAO Zhiqi , LIU Congqiang . Erosions on the southern Tibetan Plateau: Evidence from in-situ cosmogenic nuclides 10Be and 26Al in fluvial sediments[J]. Journal of Geographical Sciences, 2022 , 32(2) : 333 -357 . DOI: 10.1007/s11442-022-1950-4

1 Introduction

In general, hillslope gradients are thought to be the predominant regressor in accommodating surface erosion processes (Portenga and Bierman, 2011). Up to a threshold slope, increased frequency of landsliding can facilitate erosion processes without further steepening hillslopes. Furthermore, an increased hillslope gradient will increase stream power which in turn facilitates the transportation of river sediments and raises the potential of fluvial incision. Several studies proposed that erosion rates are spatially dominated by regional gradients in precipitation more than hillslope steepness (Montgomery and Brandon, 2002), whereas others stressed the predominant control on spatial erosion patterns by intense tectonic deformation (Finnegan et al., 2008). This long-lasting debate has attracted more attention in Earth sciences over the last several decades.
The Tibetan Plateau (TP) is characterized by active continental collision as well as highly variable precipitation (Abrahami et al., 2016). As a result, it is ideally suited as a natural laboratory for the study of feedbacks between erosion rates, tectonics and climate. Previous studies suggested that during the last few million years, climate change may have acted as a decisive force in controlling erosion processes in the Himalayan Mountains belt (Vance et al., 2003). On millennial time scales (e.g., Holocene), a number of studies documented an increase in erosion rates during intensified monsoon phases (Bookhagen et al., 2005). Decadal erosion patterns also exhibit a strong correlation with precipitation gradients (Burbank et al., 2012). However, in the central Nepalese Himalayas, the spatial patterns of long-term exhumation revealed by low-temperature thermochronology seem to be decoupled with precipitation, but well correlated with tectonics (Burbank et al., 2003). In the central Himalayas, recent studies quantified millennial-scale erosion rates with cosmogenic 10Be (Godard et al., 2014). Their results revealed a similar spatial pattern for millennial-scale erosion and long-term exhumation (derived from thermochronology), as well as indicating a predominant control by tectonics.
In this paper we quantify the basin-scale erosion rates across the middle Yarlung Zangbo River Basin (YZRB), deduced from cosmogenic 10Be and 26Al in river-borne sediments. By combining topographic analysis and basin-scale erosion rates, we aim to investigate the relative contributions of topography (e.g., relief) and climate (e.g., precipitation) to erosion processes of the southern TP. Furthermore, the 10Be-derived erosion rates are compared to the present-day hydrological erosion rates in order to assess whether or not there is any temporal variations of surface erosion.

2 Study area

The Yarlung Zangbo River is one of the longest and most powerful rivers in the world (Figure 1). It originates at an elevation of ~5200 m from Gyima Yangzoin Glacier on the south-central TP and drains a large area of ~2.4 × 105 km2 with average elevation above 4600 m in the central and eastern Himalayas and the southern TP. The river flows eastward along the Indus-Yalu suture (IYS) zone before sharply bending to the south at the Namcha Barwa-Gyala Peri massif (NBGPm). It carries a large volume of sediment to the ocean, serving as the second largest river in China (Milliman and Meade, 1983). The middle reaches of the Yarlung Zangbo River are bounded to the south by the Himalayas and to the north by the Nyainqentanglha Mountains. This part of the river is characterized by alternating broad valleys and narrow gorges (Shi et al., 2018). Compared to the narrow gorges, the broad valley sections are usually favoured for sediment storage. It has been reported that the maximum thickness of these deposits could reach up to ~800 m (Wang et al., 2015).
Figure 1 Shaded relief map showing sampling sites (square) in this work and sample locations (circles) in previous studies. Green squares refer to tributary samples and white squares to mainstream samples in this study; colored circles refer to basin or bedrock samples in previous studies; only tributary basins in this study are shown in light yellow. Sample IDs are consistent with those listed in Tables 1 and 2.
Owing to the influence of the Indian Ocean monsoon system as well as the plateau topography, the climate of the YZRB varies significantly in the longitudinal direction. Precipitation distributes unevenly across the YZRB, being concentrated in the downstream area but gradually decreases westwards. The west portion of the drainage basin is shielded from significant monsoonal rainfall due to the rain shadow effect from the Himalayas (annual precipitation < 300 mm). River discharge and sediment flux in the downstream area show more intensive peaks, which are almost consistent with fluctuations of more frequent rainfall events (Shi et al., 2018). There are a total of 10,816 glaciers developed within the YZRB, accounting for ~1.6% of the total drainage area. Most of the glaciers (~80% in the glacial zone) are distributed on the northern flank of the Yarlung Zangbo River (Guo et al., 2015).

3 Samples and methods

3.1 Sampling strategy

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 km2 (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.
Table 1 Sample information and parameters characterized for each catchment used in Figures 3 and4 -8

'Derived with the 90 m resolution SRTM DEM in ArcGIS.

"Calculated from the WorldClim V. 2 monthly climate dataset (1 km2 resolution), covering precipitation from 1970 to 2000 (Fick and Hijmans, 2017).

* Calculated following Perez- Pea et al. (2016)..

**** Calculated based on the Global Lithological Map database (GLiM) V. 1.0 (Hartmann and Moosdorf, 2012) after Moosdorf et al. (2018).

**Derived from the MODIS dataset with a spatial resolution of 500 m (https://pdaac.ugs.gov).

Samples were processed for in situ 10Be and 26Al 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/HNO3. 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 9Be carrier were added to each sample. After ion exchange column chemistry and hydroxides precipitation, Be(OH)2 and Al(OH)3 were oxidized in furnace at 800 ℃ and then mixed with Nb (BeO:Nb=1:6 by weight) and Ag (Al2O3: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 10Be/9Be ratios and 26Al/27Al 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 14C AMS system for 10Be and 26Al. Repeated measurements of secondary standards indicted ~3% precision on samples with 10Be/9Be and 26Al/27Al ratios of 10-13-10-12. The detction limit for 10Be/9Be and 26Al/27Al 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).
10Be/9Be ratios were calibrated against the NIST SRM4325 standard and 26Al/27Al 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 10Be or 26Al concentrations. In this study, we took 1.387 ± 0.012 Myr for10Be half-life (Chmeleff et al., 2010) and 0.708 ± 0.017 Myr for 26Al half-life (Nishiizumi, 2004).

3.2 Calculation of basin-scale erosion rate

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 10Be basin-scale erosion rates based on the time-independent scaling model of Lal (1991)/Stone (2000) (St scaling model), which is calibrated to a 10Be 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 10Be and 26Al 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.

3.3 Hydromorphologic analysis

3.3.1 Hillslope angles and local relief

The topography of each basin was analysed with the 90 m resolution DEM in MATLAB™, ArcGIS and the TopoToolbox Version 2.0 (Schwanghart and Scherler, 2014a). After the delineation of each drainage basin, we computed their corresponding area (km2). Other topographic parameters such as mean basin slope (degrees) and mean basin local relief (m) were subsequently determined following Jarvis et al. (2008) (Table 1).

3.3.2 Annual precipitation

The annual precipitation was derived from the WorldClim Version 2 monthly climate dataset (approximately 1 km2 resolution (Fick and Hijmans, 2017) (Table 1), which covers precipitation data from 1970 to 2000. The dataset was produced by interpolating monthly meteoric data from worldwide weathering stations with covariates obtained with the MODIS satellite platform.

3.3.3 Modified Chi value (χ°)

Drainage systems evolve with tectonic deformation and climate change. Hence, their spatial patterns can be used to interpret the potential interactions between tectonics and climate (Giachetta and Willett, 2018). The dynamics of these interactions usually control the migration of watersheds through river capture. More often, the chi (χ) value, derived by integrating the contributing drainage area with upstream channel length, is used to interpret the geometric characteristics of fluvial channels. With the help of chi (χ) plot, disequilibrium along the river channels resulted from tectonic activities or climate changing can be directely highlighted (Winterberg and Willett, 2019). For basins with similar elevation difference, lower χ value means larger steepness index or higher erosion rates. Field work and numerical models have both shown that watersheds migrate towards areas with lower erosion rates (Scherler et al., 2014). Therefore, this feature can be considered as a good indicator to basin-scale erosion. However, chi (χ) value might be influenced by varations in precipitation distribution on river discharge. To eliminate this influence in this study, we applied the modified Chi (χ) value (χ') proposed by Giachetta and Willett (2018) which also included additional correction of base level in closed basins.

3.3.4 Hypsometric integral (HI)

Hypsometry is another useful indicator that describes the topographic surface area distribution with respect to elevation. It is proved to be strongly corrolated with tectonic uplift rates and lithology differences (Pérez-Peña et al., 2009). It is often assessed with the the hypsometric integral (HI) or the hypsometric curve. The shape of hypsometric curves gives information about the volume of original basin remains, and indicates the erosional stage of drainage basins as an outcome of interaction between tectonic activity, lithological heterogeneity and climatic conditions (Huang and Niemann, 2006). According to Pérez-Peña et al. (2009), convex hypsometric curves typically represent young and slightly eroded areas, while S- and concave-shaped curves characterize moderately to highly eroded areas with mature or final (old) stage erosion. The value of hypsometric integral (HI) is equal to the area below the curve. Therefore, HI values near to 1 (convex curve) means that the mean elevations are close to the maximum elevations thus indicating a young transient landscape. In contrast, values of HI close to 0 (S- and concave-shaped curve) specifically describe mature or final (old) erosional stage areas (Pérez-Pea et al., 2016). High values of HI are usually addressed to erosional processes associated with dissected fluvial areas, while low values of HI are common in areas with prevailing deflation and fluvial aggradation processes (Ruszkiczay-Rüdiger et al., 2009). Moreover, Lifton and Chase (1992) reported that the HI values are positively correlated with tectonic uplift rate. To evaluate the stage of geomorphic development of the landscapes in our study area, we extracted the hypsometric curves and computed corresponding HI values for each basin using ArcGIS (Pérez-Peña et al., 2009).

3.3.5 Rock erodibility index (EI)

Rock erodibility index (EI) was often used to describe basin-scale erodibility potential of different lithologies (Moosdorf et al., 2018). Based on the Global Erodibility Index (Moosdorf et al., 2018) and the Global Lithological Map database (GLiM) Version 1.0 (Hartmann and Moosdorf, 2012), we calculated the mean EI values for each sampled basin (Table 1). Basically, rocks with different EI values can be divided into three groups. Rocks with EI values ranging from 1.0 to 1.2 exhibit low erodibility (e.g., acid igneous rocks). EI values of 1.4-1.5 refer to rocks of medium erodibility (e.g., basic igneous rocks). EI values near to 3.2 indicate rocks of high erodibility (e.g., unconsolidated sediments).

3.3.6 Normalized difference vegetation index (NDVI)

The normalized difference vegetation index (NDVI) is a remotely-sensed measure of vegetation density. Studies on connections between lithology and vegetation reveals that erosion rates might vary significantly across lanscapes with different vegetation densities. Hahm et al. (2014) analysed global cosmogenic nuclide data (n=691, covering all the 7 continents) and found that 10Be derived erosion rates in vegitation-covered lanscapes are generally faster than that in bare rocks. To explore the possible impact of vegetation on derived erosion rates, we extracted the NDVI values for sampled basins from the MODIS dataset (500 m spatial resolution; https://lpdaac.usgs.gov). Valid values of NDVI fall within -1 to 1 with negative values indicating cloud, water or snow areas, 0 for bare rock, and 1 for significantly dense vegetation cover (Table 1).

4 Results

4.1 Basin-scale erosion rate

Repeated measurements of standards indicate 1σ uncertainties of < 3% for samples with 10Be/9Be and 26Al/27Al ratios of 10-13 in this study. Measured 10Be/9Be ratios vary from (1.67 ± 0.06) × 10-13 to (2.25 ± 0.06) × 10-12, whereas 26Al/27Al ratios vary from (3.89 ± 0.11) × 10-13 to (4.01 ± 0.10) × 10-12. 10Be and 26Al concentrations show a 10-fold spread from (2.28 ± 0.07) × 105 to (2.26 ± 0.07) × 106 atoms g-1 and from (1.44 ± 0.05) × 106 to (2.03 ± 0.05) × 107 atoms g-1, respectively (Table 2). Correspondingly, 26Al/10Be 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 26Al/10Be ratio (8.96 ± 0.37), which is extremely higher than the canonical spallation production ratio of 6.75 at surface. The high 26Al/10Be ratio may have resulted from incorporation of deeply sourced materials or suggest overestimation of 27Al measurement. Since sample T1 yields the highest native 27Al concentration (287.60 ppm), it seems that overestimation of 27Al measurement should be accountable for the high 26Al/10Be ratio of T1. In any case, 10Be concentration and thus resulted erosion rate of sample T1 are still reliable given the high quality of determination of 10Be/9Be ratios (discussed in Section 5.1).
Figure 2 Plot of 26Al/10Be ratios versus 10Be concentration. To allow data from samples with different production rates to be compared on the same diagram, 10Be concentrations are normalized to sea level high latitude (SLHL).
Table 2 Analytical results of cosmogenic radionuclide 1'Be and 26A1 for river sediments from southern Tibet

*ID = Sample number shown in Figure 1; T refers to tributary and M to mainstream.

**"Topographic shielding factor calculated following Li (2013) with 5° in elevation angles and 10° interval in azimuth.

***"Timescale of integration for erosion (T=z*/erosion rate, with z* the mean attenuation path length; a value of 60 cm was assigned for all the samples according to Lal, 1991 and Von Blanckenburg, 2005_ ENREF 10_ENREF_14)

The 10Be 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).

4.2 Hydromorphologic indices

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 10Be-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.
Figure 3 Characteristics of sampled basins. (a) χ-transformed river profiles, derived using TopoToolbox v2 (Schwanghart and Scherler, 2014b). The black line in each river network denotes mainstream. Corresponding tributaries are shown in grey. The red line represents the linear fitting for the drainage network. Each drainage network is positioned along the abscissa based on the N90-projected distance from their outlet to sample T1. (b) 3 km radius local relief; (c) annual precipitation; (d) rock erodibility index (EI), calculated from the Global Lithological Map (GLiM) v1.0 after Moosdorf et al. (2018); (e) normalized difference vegetation index (NDVI); and (f) 10Be-derived basin-scale erosion rate. Errors in Figures 3b-3f are reported at 1σ level.
Figure 4 Hypsometric curves for the 15 drainage basins in this study. Total elevation (H) denotes the relief within the drainage basin; total area (A) denotes the total area of the drainage basin; area (a) equals to the area within the drainage basin above a certain altitude (h). Convex hypsometric curves and high HI values typically represent young and slightly eroded areas, while S- and concave-shaped curves together with low HI values characterize moderately to highly eroded areas with mature or final (old) erosional stage. The value of hypsometric integral (HI) equals to the area below the curve (Pérez-Pea et al., 2016).
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, 10Be-derived erosion rates can integrate over a period of 103-105 yr (Granger et al., 1996). By comparing 10Be 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 10Be-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.
Table 3 Information about the mainstream hydrologic stations and calculated hydrologic erosion rate*
Station Latitude
Drainage area above
the hydrologic
station (km2)
(108 m3)
sediment load
(t km-2 yr-1)
erosion rate
(m Myr-1)
Lhaze 29.1167 87.5833 4000 47832 273 52 39.4 14.6
Nugesha 29.3000 89.7500 3720 106060 378 158 217.0 80.4
Yangcun 29.2667 91.8167 3500 151507 361 306 329.8 122.1
Nuxia 29.4667 94.6500 2780 191235 480 576 155.5 57.6

*Data in columns 2-8 are from Shi et al. (2018). Hydrologic erosion rates in column 9 are calculated by dividing annual sediment load with a density of 2.7 g cm-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, 10Be 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 10Be concentrations from great depth to the surface. Addition of these materials to streams will consequently raise the apparent 10Be 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 26Al/10Be 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 10Be basin-scale erosion rate would increase due to incorporation of low-10Be -concentration sediments.

5 Discussion

5.1 Influences of deep-sourced materials on 10Be-derived basin-scale erosion rates

In the use of 10Be 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 10Be 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 10Be and 26Al. As shown in Figure 2, except for sample T8 which failed in 26Al 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 26Al/10Be 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 26Al/10Be, depth, and rock density, Akçar et al. (2017) demonstrated that the 26Al/10Be 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 10Be and 26Al. Materials originating from a greater depth are very likely to be exhumed to the surface in landscapes vulnerable to deep erosion, thus increasing the 26Al/10Be ratios in sampled river sediment. Addition of deeply sourced materials can not only increase the 26Al/10Be ratio in river sediment but will also lower their 10Be concentration due to low nuclide production rate by muons. However, sample T1 yields the highest 10Be concentration over the other 14 samples, implying no dilution of 10Be by deep-sourced materials. Alternatively, overestimation of 27Al can also lead to a high 26Al/10Be ratio. Noting that, sample T1 has the highest native 27Al concentration. In this case, overestimation of 27Al seems to be accountable for the high 26Al/10Be 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, 10Be concentration and thus, derived erosion rate of basin T1 are still reliable.
The addition of 26Al 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 10Be 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 26Al/10Be 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.

5.2 Topographic, climatic and lithologic controls on the spatial distribution of erosion rate

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).
Figure 5 Swath profile from SW to NE in Figure 1, showing the variations of (a) elevation; (b) slope; (c) annual precipitation; (d) hypsometric integral (HI) for sampled basins in this paper (squares) and in previous studies (cycles). Colors of different gradient refer to erosion rate of different magnitude.
Figure 6 Swath profile from SW to NE in Figure 1: (a) elevation; (b) slope; (c) annual precipitation; (d) hypsometric integral (HI) for sampled basins in this paper (squares) and in previous studies (cycles). Colors of different gradient refer to erosion rate of different magnitudes.
Figure 7 Swath profile from SW to NE in Figure 1: (a) elevation; (b) slope; (c) annual precipitation; (d) hypsometric integral (HI) for sampled basins in this paper (squares) and in previous studies (cycles). Colors of different gradient refer to erosion rate of different magnitude.
Figure 8 Swath profile from SW to NE in Figure 1: (a) elevation; (b) slope; (c) annual precipitation; (d) hypsometric integral (HI) for sampled basins in this paper (squares) and in previous studies (cycles). Colors of different gradient refer to erosion rate of different magnitude.
Figure 9 Swath profile from SE to NW in Figure 1: (a) elevation; (b) slope; (c) annual precipitation; (d) hypsometric integral (HI) for sampled basins in this paper (squares) and in previous studies (cycles). Colors of different gradient refer to erosion rate of different magnitude
Figure 10 Scatter plot of average slope and HI values. Drainage basins are clustered in two groups (HI-1 and HI-2); higher values of HI and average slopes indicate more dissected surfaces and younger erosional stage (Ruszkiczay-Rüdiger et al., 2009; see Table 1 for HI and average slope values).
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).

5.3 Implications for landscape evolution

By far, hundreds of local or basin-scale erosion rates derived using cosmogenic 10Be 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.
Figure 11 Map of χ' for southern Tibetan Plateau, showing geometric disequilibrium in drainage networks. Difference in χ' values between two adjacent basins can indicate the migration of basin divide and exchange of basin area. In Figures b-k, the white thick line starts from A to B and crosses the divide between basin A and basin B along the valley bottoms. The maps below Figures b-k show the variations of topography following these white thick lines. Values of χ' are inversely consistent with topographic asymmetry between basin A and basin B.

6 Conclusions

This study presents 15 new basin-scale erosion rates derived from cosmogenic 10Be and 26Al in fluvial sediments across the southern Tibetan Plateau. The erosion rates range from 20.60 ±1.79 to 154.00 ± 13.60 m Myr-1, being low in the upstream areas of the Gyaca knickpoint and high in the downstream areas. Paired analysis of 10Be and 26Al reveals that all basins satisfy the assumption of steady-state erosion and thus are free from the perturbation of non-steady-state erosion processes. Both topographic and climatic indices play significant roles in erosion for basins in the upstream areas of the Gyaca knickpoint and for basins in the southern Himalayas. Variations in precipitation have only a second-order control in the downstream areas of the Gyaca knickpoint. Instead, topography takes the leading role on erosion processes for this region. The marginal Himalayas and the YZRB (this study) exhibit significantly higher erosion rates than the central plateau. In other words, the landscape of the central plateau is remarkably stable and is being intensively consumed at its boundaries through river headwards erosion. In addition, comparison of erosion rates between 10Be and hydrological records suggests that present-day hydrologic erosion rates have been accelerated in areas under severe anthropogenic perturbations.
Abrahami R, van der Beek P, Huyghe P et al., 2016. Decoupling of long-term exhumation and short-term erosion rates in the Sikkim Himalaya. Earth and Planetary Science Letters, 433: 76-88.


Akçar N, Ivy-Ochs S, Alfimov V et al., 2017. Isochron-burial dating of glaciofluvial deposits: First results from the Swiss Alps. Earth Surface Processes and Landforms, 42(14): 2414-2425.


Balco G, Stone J O, Lifton N A et al., 2008. A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements. Quaternary Geochronology, 3(3): 174-195.


Bookhagen B, Thiede R C, Strecker M R, 2005. Late Quaternary intensified monsoon phases control landscape evolution in the northwest Himalaya. Geology, 33(2): 149-152.


Borchers B, Marrero S, Balco G et al., 2016. Geological calibration of spallation production rates in the CRONUS-Earth project. Quaternary Geochronology, 31: 188-198.


Burbank D W, Blythe A E, Putkonen J et al., 2003. Decoupling of erosion and precipitation in the Himalayas. Nature, 426(6967): 652-655.


Burbank D W, Bookhagen B, Gabet E J et al., 2012. Modern climate and erosion in the Himalaya. Comptes Rendus Geoscience, 344(11/12): 610-626.


Chmeleff J, von Blanckenburg F, Kossert K et al., 2010. Determination of the 10Be half-life by multicollector ICP-MS and liquid scintillation counting. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 268(2): 192-199.

Dong K, Lang Y, Xu S, 2019. Progress in XCAMS at Tianjin University. Radiation Detection Technology and Methods, 3(4): 1-8.


Dong K J, Lang Y C, Hu N et al., 2018. The new AMS facility at Tianjin University. Radiation Detection Technology and Methods, 2(1): 30.


Fick S E, Hijmans R J, 2017. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12): 4302-4315.


Finnegan N J, Hallet B, Montgomery D R et al., 2008. Coupling of rock uplift and river incision in the Namche Barwa-Gyala Peri massif, Tibet. Geological Society of America Bulletin, 120(1/2): 142-155.


Giachetta E, Willett S D, 2018. A global dataset of river network geometry. Scientific Data, 5: 180127.


Godard V, Bourles D L, Spinabella F et al., 2014. Dominance of tectonics over climate in Himalayan denudation. Geology, 42 (3): 243-246.


Granger D E, Kirchner J W, Finkel R, 1996. Spatially averaged long-term erosion rates measured from in situ-produced cosmogenic nuclides in alluvial sediment. The Journal of Geology, 104(3): 249-257.


Guo W Q, Liu S Y, Xu J L et al., 2015. The Second Chinese Glacier Inventory: Data, methods and results. Journal of Glaciology, 61 (226): 357-372.


Hahm W J, Riebe C S, Lukens C E et al., 2014. Bedrock composition regulates mountain ecosystems and landscape evolution. Proceedings of the National Academy of Sciences of the United States of America, 111(9): 3338-3343.

Hartmann J, Moosdorf N, 2012. The new global lithological map database GLiM: A representation of rock properties at the Earth surface. Geochemistry, Geophysics, Geosystems, 13(12).

Hewawasam T, von Blanckenburg F, Schaller M et al., 2003. Increase of human over natural erosion rates in tropical highlands constrained by cosmogenic nuclides. Geology, 31(7): 597-600.


Hodges K V, Wobus C, Ruhl K et al., 2004. Quaternary deformation, river steepening, and heavy precipitation at the front of the Higher Himalayan ranges. Earth and Planetary Science Letters, 220(3/4): 379-389.


Huang X J, Niemann J D, 2006. An evaluation of the geomorphically effective event for fluvial processes over long periods. Journal of Geophysical Research: Earth Surface, 111(F3).

Jarvis A, Reuter H I, Nelson A, et al., 2008. Hole-filled SRTM for the globe Version 4. available from the CGIAR-CSI SRTM 90m Database(15: 25-54.

Kohl C P, Nishiizumi K, 1992. Chemical isolation of quartz for measurement of in-situ-produced cosmogenic nuclides. Geochimica et Cosmochimica Acta, 56(9): 3583-3587.


Lal D, 1991. Cosmic ray labeling of erosion surfaces: In situ nuclide production rates and erosion models. Earth and Planetary Science Letters, 104(91): 424-439.


Lang K A, Huntington K W, Montgomery D R, 2013. Erosion of the Tsangpo Gorge by megafloods, eastern Himalaya. Geology, 41(9): 1003-1006.


Li Y K, 2013. Determining topographic shielding from digital elevation models for cosmogenic nuclide analysis: A GIS approach and field validation. Journal of Mountain Science, 10(3): 355-362.


Li Y K, Li D W, Liu G N et al., 2014. Patterns of landscape evolution on the central and northern Tibetan Plateau investigated using in-situ produced 10Be concentrations from river sediments. Earth and Planetary Science Letters, 398: 77-89.


Lifton N A, Chase C G, 1992. Tectonic, climatic and lithologic influences on landscape fractal dimension and hypsometry: Implications for landscape evolution in the San Gabriel Mountains, California. Geomorphology, 5(1/2): 77-114.


Lupker M, Lavé J, France-Lanord C et al., 2017. 10Be systematics in the Tsangpo-Brahmaputra catchment: The cosmogenic nuclide legacy of the eastern Himalayan syntaxis. Earth Surface Dynamics, 5(3): 429-449.


Milliman J D, Meade R H J J o G, 1983. World-wide delivery of river sediment to the oceans. Journal of Geology, 91(1): 1-21.


Montgomery D R, Brandon M T, 2002. Topographic controls on erosion rates in tectonically active mountain ranges. Earth and Planetary Science Letters, 201(3/4): 481-489.


Moosdorf N, Cohen S, von Hagke C, 2018. A global erodibility index to represent sediment production potential of different rock types. Applied Geography, 101: 36-44.


Nishiizumi K, 2004. Preparation of 26Al AMS standards. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 223/224: 388-392.

Pérez-Pea J V, Al-Awabdeh M, Azaón J M et al., 2016. SwathProfiler and NProfiler: Two new ArcGIS Add-ins for the automatic extraction of swath and normalized river profiles. Computers & Geosciences, 104: 135-150.


Pérez-Peña J V, Azañón J M, Booth-Rea G et al., 2009. Differentiating geology and tectonics using a spatial autocorrelation technique for the hypsometric integral. Journal of Geophysical Research, 114(F2).

Portenga E W, Bierman P R, 2011. Understanding Earth's eroding surface with 10Be. GSA Today, 21 (8): 4-10.

Ruszkiczay-Rüdiger Z, Fodor L, Horváth E et al., 2009. Discrimination of fluvial, eolian and neotectonic features in a low hilly landscape: A DEM-based morphotectonic analysis in the Central Pannonian Basin, Hungary. Geomorphology, 104(3/4): 203-217.


Scherler D, Bookhagen B, Strecker M R, 2014. Tectonic control on 10Be-derived erosion rates in the Garhwal Himalaya, India. Journal of Geophysical Research: Earth Surface, 119(2): 83-105.


Schmidt A H, Neilson T B, Bierman P R et al., 2016. Influence of topography and human activity on apparent in situ 10Be-derived erosion rates in Yunnan, SW China. Earth Surface Dynamics, 4(4): 819-830.


Schwanghart W, Scherler D, 2014a. Short Communication: TopoToolbox 2 - MATLAB-based software for topographic analysis and modeling in Earth surface sciences. Earth Surface Dynamics, 2(1): 1-7.


Schwanghart W, Scherler D, 2014b. TopoToolbox 2-MATLAB-based software for topographic analysis and modeling in Earth surface sciences. Earth Surface Dynamics, 2(1): 1-7.


Shi X N, Zhang F, Lu X X et al., 2018. Spatiotemporal variations of suspended sediment transport in the upstream and midstream of the Yarlung Tsangpo River (the upper Brahmaputra), China. Earth Surface Processes and Landforms, 43(2): 432-443.


Stone J O, 2000. Air pressure and cosmogenic isotope production. Journal of Geophysical Research-Solid Earth, 105(B10): 23753-23759.


Vance D, Bickle M, Ivy-Ochs S et al., 2003. Erosion and exhumation in the Himalaya from cosmogenic isotope inventories of river sediments. Earth and Planetary Science Letters, 206(3/4): 273-288.


Von Blanckenburg F, 2005. The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment. Earth and Planetary Science Letters, 237(3): 462-479.


Wang Z, Yu G, Wang X et al., 2015. Sediment storage and morphology of the Yalu Tsangpo valley due to uneven uplift of the Himalaya. Science China Earth Sciences, 58(8): 1440-1445.


West A J, Hetzel R, Li G et al., 2014. Dilution of 10Be in detrital quartz by earthquake-induced landslides: Implications for determining denudation rates and potential to provide insights into landslide sediment dynamics. Earth and Planetary Science Letters, 396: 143-153.


Winterberg S, Willett S D, 2019. Greater Alpine river network evolution, interpretations based on novel drainage analysis. Swiss Journal of Geosciences, 112(1): 3-22.


Zondervan A, Hauser T, Kaiser J et al., 2015. XCAMS: The compact 14C accelerator mass spectrometer extended for 10Be and 26Al at GNS Science, New Zealand. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 361: 25-33.