The geomorphologic evolution of catchments is the basis for analyzing the changes in the surface environment. The incision rate and spatiotemporal variation of a river network are core and difficult points in the analysis of geomorphic characteristics given a long-timescale background, which is helpful to understand the historical evolution of a basin and the sediment-contributing area (Demoulin et al., 2017; Wu et al., 2020). The advancement of geochronological techniques has also made it possible to measure the formation ages of river terraces and long-term erosion rates (Guo et al., 2012; Hu et al., 2016; Zhang et al., 2018). Reconstructing ancient channel profiles through river terraces and using sediment dating results to estimate river erosion rates are considered effective methods of achieving this. Common dating methods for sedimentary strata include thermoluminescence (TL), photoluminescence (OSL), electron spin resonance dating (ESR), cosmogenic nuclide (CN), and others (Kreutzer et al., 2018; Roda-Boluda et al., 2019; Su et al., 2019; Tecsa et al., 2020). In recent years, there have been breakthroughs in low-temperature thermochronology, typified by apatite and zircon dating, which have been made, and they have been used in the analysis of plateau residual geomorphology and inversion of paleo-topography (Reiners et al., 2006; Wang et al., 2018). Using these dating methods to determine the age of sediments on landforms such as river terraces, alluvial-flood fans, and fault cliffs, we could reconstruct the formation and evolution history and quantify the erosion rate of landforms (Wu et al., 2020). However, basins in areas with an intense orogenic belt and humid climate have active erosion characteristics in which river valleys are mostly V-shaped and terrace sediments are rarely preserved; thus, the quality and collection process of sediment samples can be adversely affected. Besides, the drainage reorganization in the basin, which is accompanied by river erosion that includes drainage divide migration and river capture (Reiners et al., 2006; Portenga et al., 2011; Singh et al., 2017), is often neglected because of the large spatial scale of the river basin and the sparse field relics, such as beheaded rivers and pirate bays.

The Buyuan River basin, located in the upper Lancang-Mekong, is a large branch basin with such characteristics and is an important source of sediment for Lancang River trunk streams (Gu et al., 2019). For this basin, it is not always possible to fully comprehend the erosion rate of the whole channel by studying river terraces preserved in the local river reaches. The determination of the incision rate of a river and its divide migration are related to the stable evaluation of existing river engineering and the practical need to select the best site for major projects, and also indicate the characteristics of Quaternary geomorphologic evolution. In order to analyze the incision rate and river reorganization trend in this basin, this study starts from the perspective of geomorphology and reveals the geomorphologic characteristics of this basin through the calculation and interpretation of quantitative geomorphologic indicators, including the steep index (Shi, 2019), morphological dating index (Demoulin, 2011, 2012), and Chi-plots (χ-z) (Willett et al., 2014; Sean, 2018). We determined the incision rate of the river based on a reconstruction of the equilibrium longitudinal profile of the main channel and a determination of the time of the last base-level decline. Thereafter, we discuss the migration trend of sub-basins by calculating the spatial distribution characteristics of the Chi-value distribution.

The Buyuan River is an important tributary of the upper Lancang-Mekong, the river basin covers an area of about 8931 km² with an elevation range of 500-2290 m. The area was characterized by the interaction of coastal continental and transgression facies in the Early Paleozoic period. During the Late Paleozoic period, the basin area was a marine environment, mainly composed of limestone and dolomite strata, evolving into a terrestrial environment with strata and argillaceous sediments, during the early Mesozoic era (Ma, 2002). The majority of the faults in the basin are strike-slip faults, some of which are normal faults with a strike-slip rate of about 0.2-4 mm/yr, and the main fault activities commenced in the Pliocene epoch (Wu, 2016). Although the modern surface outcrop in the study area is a sedimentary layer formed in the Mesozoic period, the lithological characteristics mainly differ in age. A large area of Cretaceous strata is exposed in the north, and older strata are exposed in the south; later Paleozoic strata are exposed in the local area downstream (Figure 1). Currently, the basin is affected by a subtropical monsoon climate, with an annual average temperature of around 21.5℃ and annual average precipitation of 1488 mm, mainly concentrated between June and November. River discharge is recharged by precipitation, and the annual average value is about 178 m³/s, accounting for 1.2% of the total discharge of the Lancang-Mekong River (He, 1995).

The research is mainly based on Digital Elevation Models (DEMs) with a~30 m spatial resolution and a~5 m height error (Elkhrachy, 2018; Zhang et al., 2019) provided by the Shuttle Radar Topography Mission (SRTM) from the National Aeronautics and Space Administration (http://earthexplorer.usgs.gov).

3.2.1 Morphometric dating

The morphological dating method was proposed by Demoulin (2011, 2012). The response of erosion to a base-level decline begins in the main channel and is transmitted through the river network to the entire basin due to the traceable migration of knickpoints. The response time of the longitudinal profile of the mainstream, the river network and the whole basin to the decline of base-level corresponds to a short, medium and long-time scale, respectively. Therefore, the metric R^* designed by Demoulin (2011, 2012) can be presented as follows:

where $\mathop{\int }^{}{{H}_{\text{r}}}$, $\mathop{\int }^{}{{H}_{\text{n}}}$ and $\mathop{\int }^{}{{H}_{\text{b}}}$ represent the hypsometric integral of the mainstream, the hypsometric integral of drainage network and the hypsometric integral of the basin, respectively (Demoulin, 2011; Demoulin et al., 2017; Shi, 2019). These can be calculated using the following formula:

Where$\bar{H},$ H_min, and H_max represent the mean, minimum, and maximum elevation of longitudinal profiles, the river network and the basin, respectively (Pike et al., 1971; Shi, 2019). Catchments of different sizes can be extracted along the main channel and the R^* value of each catchment can be calculated by Eqs. (1) and (2) (Demoulin et al., 2017; Shi, 2019). As the influences of lithology and climate on a river network and the whole basin are similar (Demoulin, 2012), the ratio operation in Eq. (1) cancels out the influences to the greatest extent and R^* is mainly affected by the irregularity of the basin plane morphology (Demoulin et al., 2017). However, R^* is corrected by multiplying a factor of $1/\sqrt{E}$ to eliminate the effect of the river network on basin elongation and branching level. E is calculated as 4A/(πL_b²), in which A (km²) is the catchment area and L_b (km) is the maximal length of a basin, measured from its outlet. An empirical relation established according to known dates of uplift episodes in many regions of the world is:

where t represents the starting time (Ma), S_R is the slope of the linear fitting function of R^* = S_R (Ln A) + b in a region of homogeneous uplift (Demoulin, 2012). The morphological dating method requires that the area of a tectonic uplift block should not be less than 1000 km² to ensure a stable relationship between R^* and Ln A (Demoulin et al., 2017). The extraction of each index can be conducted on the ArcGIS 10.2 software platform, and the extraction threshold of the river network was set to 0.5 km².

3.2.2 Reconstruction of a paleo-profile in a stable state

The local slope of the river channel in an eroded basin can be considered to be mainly controlled by the tectonic uplift rate and the catchment area under natural conditions. Therefore, the following relationship exists:

where k_s and θ are steepness index and concavity index respectively (Flint, 1974; Snyder et al., 2000). For regions with consistent bedrock-climatic characteristics and uplift rates, the θ value of the longitudinal profile in a steady state is generally between 0.4 and 0.6 and is determined by the downstream rates of increase in the discharge and channel width (Kirby et al., 2012). As the term “knickpoint” represents the position where the slope of the longitudinal profile changes abruptly, and the traceable movement of knickpoints is usually considered a response to erosion base-level decline (Whipple, 2001), we consider that the current position of a knickpoint is due to an incision wave associated with the base-level fall by comparing the current longitudinal profile with the reconstructed longitudinal paleo-profile and assume that the original paleo-profile was in an equilibrium state in our analysis (Pastor et al., 2015). Based on these reasonable assumptions, the steepness index and concavity index of upper reaches of the knickpoint (Flint, 1974), together with catchment areas along the mainstream and the downstream distance, x (Hack, 1957), the paleo-longitudinal profile can be reconstructed by the following equation (Whipple et al., 1999; Whipple, 2001; Pastor et al., 2015):

where z is the elevation at a distance x; L is the distance from the outlet of river to the drainage divide; x_c is the distance from the divide to the point at which fluvial processes become dominant over gravitational erosion processes; k_a and h are Hack’s coefficient exponent respectively, and the other parameters have the same meaning as they did above.

3.2.3 Drainage reorganization

Drainage networks are adjusted to changing tectonic activity and climatic conditions until a balance is achieved between the catchment areas rather than maintaining static topographic features (Sean et al., 2017). In a basin, rivers on both sides of the divide have the same erosion datum and elevation difference. In this case, a river reach with a larger Chi (χ) value usually has a lower k_s value, and a low k_s value usually implies a low erosion rate. If there are no significant differences in lithology on both sides of the divide, the divide will then move to the side with a low erosion rate (Willett et al., 2014). The χ can be calculated using the following equations:

where U, K, A and A₀ are uplift rate, erosion rate, catchment area and a reference area respectively; z and z_b are the elevation of a river profile at an observation point indicated by x, and the elevation of the profile at a base level indicated by x_b, respectively. When A₀ = 1 m², the slope of χ-z is k_s (Perron et al., 2013). The Chi-plot calculation is finished by code based on Eq. (7) in the ArcGIS 10.2 software platform. According to the Strahler system, we divide the basin into several sub-basins and appropriately merge the broken areas with very small sizes to facilitate the exploration of river reorganization in the basin.

A drainage network is a dynamic system that adjusts to perturbations in the landscape, such as tectonic deformations and climate changes, to reach an equilibrium state between fluvial erosion and base level variations (Struth et al., 2019). Knickpoints represent a transient response to external perturbations. The longitudinal profile of the mainstream shows that there are two obvious knickpoints 42 km and 144 km away from the river source, reflecting the disturbance in the longitudinal profile caused by two large-scale base-level descending events (Figure 2a). Since knickpoint-I is located near the river source, we can conclude that the basin has almost completed the main response phase to the first decline in base level. Utilizing dating based on morphology methods related to estimating the latest decrease in data time points in conjunction with the empirical relationship of t-S_R (Demoulin et al., 2017), we can date knickpoint-II to around 46 ka BP. (Figure 2b).

The channel upstream of knickpoint-I has undergone a long period of erosion evolution and has been in a steady state. The steep index of the steady-state profile is extracted via Chi-plot (Figure 3a), and then the downstream steady-state profile was reconstructed according to Eq. (5), which could reflect the steady-state bed elevation distribution before the first phase of base-level decline (Figure 3b). Additionally, it also can find another reach close to the steady state in knickpoint-II that is located 60-131 km downstream of the river source. According to the extracted steady-state parameters, k_s can also be used to reconstruct the downstream steady-state paleo-longitudinal profile (Figure 3b). At the estuarine location, the elevation differences between the two-stage steady-state longitudinal profiles and the modern riverbed are 299 m and 151 m.

River capture is an inevitable phenomenon in the process of transient geomorphology adjustment, which will lead to changes in the spatial pattern of a drainage system in the basin. Depending on Chi-maps, the degree of disequilibrium of the drainage divides and the trend of horizontal motion needed to reach an equilibrium state can be visually reflected (Struth et al., 2019). We calculated the Chi-value distribution in the basin by taking the river outlet as x_b so that the height differences between the two sides of the sub-basin and the river outlet are the same. In this way, the Chi-value can indirectly reflect the difference in k_s. The results indicate that the two sides of the divide line have significant differences in Chi-distributions in the northern part of the basin, especially at the intersecting location of 3-4 sub-basins (Figure 4). If the influence of lithological differences on the headstream erosion rate is not taken into account, it can be assumed that the headstream erosion rates in different directions in these areas are unbalanced. Besides, Chi-values on both sides of the divide are roughly equal, and the divides without blue arrows in Figure 4a are basically in equilibrium.

The longitudinal profile gives more accurate results compared to the traditional field obser-vations and is the most sensitive geomorphic feature, which maintains evidence of terrain deformation and structural changes in geologic past (Das, 2018). A base level decline can drive a knickpoint migration traceability (Struth et al., 2019). According to the longitudinal profile characteristics, we infer that, after the first phase of base-level decline in the basin, another base-level decline event occurred in the downstream riverbed when the riverbed reached a steady-state during the receding process of knickpoint-II, around 46 ka BP. Since a comparison between the equilibrium longitudinal profile and the current longitudinal profile can reflect the base level decline (Shi, 2019), combining this comparison with the elevation difference between the reconstructed equilibrium longitudinal profile and the current longitudinal profile can allow us to estimate that the incision rate of Buyuan River since the Late Pleistocene is about 0-3.28 mm/yr, during which the position of the maximum incision rate occurred in the river outlet. Given the background of the landform evolution in the Lancang River basin, new dating results of OSL, (U-Th)/He (TL) and ¹⁰Be (CN) on river terrace sediments indicate that the Lancang River has experienced four aggradation events at > 120-100 ka, 90-70 ka, 25-15 ka, and < 9 ka, with each event followed by a rapid incision at 100 ka BP, 45 ka BP, 15 ka BP, and 6 ka BP (Zhang et al., 2018).

As a large tributary of the Lancang River, the Buyuan River is likely to be affected by these significant environmental changes at the corresponding stage. The estimated time of the most recent erosion base-level decline event (46 ka BP) in the Buyuan River basin is close to that of the second rapid cutting event (45 ka BP) of the Lancang River in the Late Pleistocene. Due to knickpoint-I being located upstream of knickpoint-I, it is more likely to be affected by an earlier erosion base-level decline. Therefore, there is a possibility that the rapid cutting event of Lancang River at 100 ka BP has a corresponding time relationship with knickpoint-II. It can be inferred that, during 100-46 ka BP, the incision rate of Buyuan River might have been 0-2.74 mm/yr according to the ratio of elevation difference to time, and the average incision rate during the Late Pleistocene would have been 0-2.99 mm/yr. These results are similar to the incision rate (0.5-2.8 mm/yr) in the Lancang reaches upstream from the intersection of Buyuan River and Lancang River (Figure 5) (Zhang et al., 2018). Note that the two rapid cutting events of Lancang River in 15 ka BP and 6 ka BP mentioned by Zhang et al. (2018) were not interpreted from morphological dating and the longitudinal profile in this analysis. This may be because these two cutting events are localized, but it is more likely that the geomorphologic characteristics of the Buyuan River basin have not yet demonstrated obvious responses to morphological change in the short time since the events, because the time scale of the response of the geomorphologic evolution of the basin to tectonic activity and climate change is generally between 10⁴ and 10⁶ years (Demoulin et al., 2017); this may also be due to the vertical accuracy limitations of the DEM.

Surface processes in the unstable state are complex, but they retain the topographic characteristics before and after the change in boundary conditions and contain information to explain the past and future geomorphic evolution trend (Whipple, 2001; Struth et al., 2019). We use Chi-plots to explore the different degrees of erosion in the Buyuan River basin and the drainage reorganization. Observing the Chi-plots, differences across the divides suggest that the spatial range of the disequilibrium region is relatively limited. In many areas, similar Chi-values at the source of river on opposite sides of a divide indicate that the drainage network is close to a geometric equilibrium state. Along the divide line of the sub-basins, the uneven parts of the Chi-value mainly appear in the intersecting area of the river sources of 3-4 sub-basins, while the lithology differences in these areas are not large (Figures 1 and 4a). Therefore, we believe that there is no severe recombination of water systems in the basin and that the imbalance of traceable erosion only exists in local areas. These analysis results show that river reorganization is not a universal phenomenon in the current geomorphic evolution in the orogenic belt areas around the Qinghai-Tibet Plateau. Geomorphologic evolution dynamics at the basin scale are derived from tectonic uplift and climate change (Whipple et al., 1999; Demoulin et al., 2017); however, recent studies have shown that there is a relatively obvious decoupling phenomenon in the formation and development of the Lancang River basin. The hydrologic driven incision of a large fluvial system in southeast Tibet is mostly likely related to local climate conditions rather than tectonically driven erosion, especially at a time scale of 10-100 ka (Nie et al., 2018; Zhang et al., 2018). Besides, the intensity of neotectonic activity in the basin is relatively weak, according to an evaluation of the tectonic geomorphology index (Gu et al., 2019), which indicates that the basin is less likely to have significant tectonic uplift in the Late Pleistocene. It seems that, since the Late Pleistocene, the change in the geomorphic evolution characteristics in the Buyuan River basin, including river cutting, base-level decline, and river reorganization, is more likely driven by precipitation.

This study presents an analysis of incision and drainage reorganization of the Buyuan River basin in the upper Lancang-Mekong River of China. The application of morphological dating, paleo-profile reconstruction, and Chi-plots has helped us to demonstrate contrasts in the equilibrium state of rivers across divides and then to infer their migration. The results show that: (1) the mainstream of the Buyuan River was affected by at least two large-scale base-level decline events and retained two large knickpoints in the Late Pleistocene. (2) The earlier stage might have occurred in 100 ka BP, and the later stage likely occurred around 46 ka BP. (3) In comparison of the equilibrium profile of the two main channels with the modern river profile, it can be inferred that the average incision rate of the mainstream since the Late Pleistocene is 0-2.99 mm/yr since 100 ka BP and 0-3.28 mm/yr since 46 ka BP. (4) Besides, characteristics in the Chi-values across the divides suggest that the spatial range of the disequilibrium region is relatively limited and that there is no severe reorganization of the river network in the basin, and the imbalance of traceable erosion in the river network exists in local areas. (5) Since the Late Pleistocene, the main driving force of the geomorphological evolution in the Buyuan River basin is likely climate fluctuations.