Mapping of topographic changes
Pre-flood topography data (SRTM, 1 arc-second, flown onboard the space shuttle Endeavour in February 200031,43) was subtracted from post-flood topography data (ALOS, 1 arc-second, launched in January 200631,43), to generate a map of elevation difference referred to below as DOD.
Before calculation and interpretation of the elevation difference, the vertical accuracy of the two different data products was compared44. The official elevations of three international airports (Qamdo Bangda Airport, Gongga Airport, and Nyinchi Airport, https://www.fsaerodata.com/) around Yigong Lake were used as control points to calibrate the elevation data (Supplementary Table 1). The results show that the vertical precision of both DEMs is within an error range of ~ 5 m45 after calibration (Supplementary Figs. 3, 10). A positive DOD value indicates accumulation, and a negative value indicates erosion by the flood (Fig. 2; Supplementary Fig. 1). We assumed that, compared with effects induced by the 2000 flood, topographic changes were negligible between February 2000 and 10 June 2000 (the date of the outburst flood) and between 10 June 2000 and the acquisition time of ALOS data in 2006 (see Supplementary Text). We used the difference map, Landsat and Google Earth images, and a field expedition to delineate topographic changes and quantify the erosional and depositional effects of the 2000 Yigong outburst flood.
The material, bedrock or alluvial sediments, eroded and excavated by the outburst flood, cannot readily be distinguished and quantified limited by the lack of real-time monitoring. Based on the continuous high shear stress (e.g. >1 kPa) during the flood and the post-flood landforms, we provisionally assume that freshly eroded areas of bedrock and alluvium visible after the flood imply that both types of material were eroded.
Flood Hydraulics
The dam-failure hydrodynamic routing with unsteady flow was modeled using the HEC-RAS 2D hydraulic model46. The simulations are configured with a nominal mesh resolution of 30 m using a time step of 1s and a simulation time of 17 hours. The input data are as follows. (1) The SRTM DEM before the flood was used as terrain input for exploring geomorphic changes induced by the flood. (2) The location of Yigong Lake and landslide dam was determined according to the remote sensing images of May 12, 2000 (Supplementary Table 2). The physical shape parameters of the landslide dam before and after the flood, and the storage capacity curve of Yigong lake (Supplementary Fig. 11d), were used to calculate the discharge after the dam breach. (3) The Manning coefficient (n) is uniformly assigned as 0.04 for the whole domain. Previous studies47–49 of the choice of n for calculating very large discharges in narrow, deep bedrock canyons show low sensitivity. (4) The inflow discharge to the lake is set to 0. The impact of this background flow on the flood depth can be ignored22 because the discharge of the outburst flood was at least 100 times larger. The HEC-RAS simulations yield hydraulic parameters that include flow depth, velocity, shear stress and discharge for each cross-section.
Shear stress was calculated using:
$$\begin{array}{c}{\tau }_{b}=\gamma RS\#\left(1\right)\end{array}$$
where \({\tau }_{b}\) is bed shear stress; \(\gamma\) is the unit weight of water (i.e., 9800 N\(\bullet\)m−3); \(R\) is hydraulic radius; \(S\) is energy slope.
According to the simulation results (Supplementary Figs. 5, 11), the flood propagation was as follows. The outburst occurred at 19:00 hr on June 10th, 2000, and 10 minutes later, the flow at the breach site reached ~ 14×104 m3/s. After 20 minutes, the flood wave reached Tongmai Bridge. At 19:40 hr on June 10, the peak flow at Tongmai Bridge was ~ 12×104 m3/s (Supplementary Fig. 11). Ten hours after the outburst, Yigong Lake was almost completely drained, and the discharge has been stable since then (Supplementary Figs. 5, 11). The calculated peak discharges at the breach and Tongmai Bridge are slightly less than those calculated by Turzewski30 and Morey22 using the GeoClaw model, which were 17.3×104 m3/s and 13.9×104 m3/s (n = 0.04), respectively. In addition, Delaney and Evans29 estimated peak discharge of 11×104 m3/s at Tongmai Bridge from a synthetic hydrograph scaled to their estimated peak breach discharge of 61,461 m3/s, using the FLO-2D model. Although the terrain data selected for simulation and methods to estimate the peak discharge are different, our calculated results are among the range of the results from previous work22,29−30. In addition, the peak discharge gradually decreases downstream (Supplementary Fig. 11), conforming to the basic law of outburst floods, supporting the validity of our hydraulic simulation.
The magnitudes of discharge, flow velocity, and other flood parameters were large compared to base flow. The calculated peak discharge of the flood was ~ 12×104 m3/s, ~ 300 times the mean annual discharge, and more than 50 times the annual peak flood during the monsoon period30,39. Downstream from Tsangpo Gorge, the peak discharge of the flood was ~ 3–6×104 m3/s (Supplementary Fig. 5), which is 15–30 times the mean annual discharge30,39. The monsoon baseflow in this area can generate water depths of ~ 5–10 m30,39, whereas the Yigong flood could have been 50 m deep, and possibly 100 m deep, through the Tsangpo Gorge region (Supplementary Fig. 11). Usually, river flow velocities rarely exceed 5 m/s, but the flood attained a velocity of around 20 m/s (Supplementary Fig. 11).
Thresholds Of Shear Stress
Our simulation produces a spatial assessment of the bed shear stress. Here, we estimated thresholds of shear stress within Po-Tsangpo Gorge for incipient motion, plucking a protruding block via sliding, and suspension using the calculation formula in previous work24,30,50. Results show the relationship between the given particle size and the thresholds of shear stress for different erosion processes (Fig. 6).
The threshold shear stress for incipient motion was calculated using:
$$\begin{array}{c}{\tau }_{c}^{*}=0.15{S}^{0.25}\#\left(2\right)\end{array}$$
$$\begin{array}{c}{\tau }_{b}=D\times {\tau }_{c}^{*}\times g\times ({\rho }_{s}-\rho )\#\left(3\right)\end{array}$$
where \({\tau }_{c}^{*}\) is the dimensionless critical shear stress for incipient motion; \(S\) is the river gradient (set to 0.01); \({\tau }_{b}\) is the threshold shear stress for movable blocks; \(D\) is the median block size; \(g\) is gravitational acceleration; \({\rho }_{s}\) is the density of granite (2700 kg/m3); and \(\rho\) is the density of water (1000 kg/m3).
The threshold shear stress for suspension was calculated using:
$$\begin{array}{c}R=\frac{({\rho }_{s}-\rho )}{\rho }\#\left(4\right)\end{array}$$
$$\begin{array}{c}{w}_{s}=\frac{Rg{D}^{2}}{{{c}_{1}v+\left(0.75{c}_{2}Rg{D}^{3}\right)}^{0.5}}\#\left(5\right)\end{array}$$
$$\begin{array}{c}{\tau }_{b}={\rho \left(0.8{w}_{s}\right)}^{2}\#\left(6\right)\end{array}$$
where \({w}_{s}\) is the setting velocity; \({c}_{1}\) is the grain shape, an empirical constant of 20; \({c}_{2}\)is roughness, an empirical constant of 1.1; \(v\) is the kinematic viscosity, an empirical constant of 10− 6 m2/s; \({\tau }_{b}\) is the threshold shear stress for suspended deposits.
The threshold shear stress for block plucking was calculated using:
$$\begin{array}{c}{\tau }_{pc}^{*}=\frac{\text{cos}\theta \left(\text{tan}\phi -\text{tan}\theta \right)+2{\tau }_{w}^{*}}{\left[1+0.5{C}_{d}{\left(\frac{u}{{u}^{*}}\right)}^{2}\frac{P}{L}\right]\left[1+{F}_{L}^{*}\text{tan}\phi \right]}\#\left(7\right)\end{array}$$
$$\begin{array}{c}{\tau }_{pc}=D\times {\tau }_{pc}^{*}\times g\times ({\rho }_{s}-\rho )\#\left(8\right)\end{array}$$
where \({\tau }_{pc}^{*}\) is the dimensionless critical shear stress for plucking; \(\theta\) is the bed angle (set to 1°); \(\phi\) is the bed friction angle (set to 35°); \({\tau }_{w}^{*}\) is dimensionless block sidewall stress (set to 0); \({C}_{d}\) is the local drag coefficient (set to 1); \(P\) is the block protrusion height (0.1 or 0.2); \(L\) is the block length; \({F}_{L}^{*}\) is the dimensionless hydraulic lift force (set to 0.85); and \(\frac{u}{{u}^{*}}\) is 8.3, the same as in Turzewski125 ; \({\tau }_{b}\)is the threshold shear stress for plucking. Here, we notice that the bed angle \(\theta\) has little influence, with the critical shear stress changing slightly (3%) when the bed angle changed from 0.5° to 1.5°. We measured several bed angle values from DEM, and finally set the median at 1°.
Concurrent Landslide Mapping And Volume Estimation
Landslide mapping. Landslide mapping along ~ 80 km downstream from the breach in 1990–2020 was based on Landsat remote sensing images under conditions of relatively clear sky, mostly from April to September. During this period, the vegetation grows well and the mountain shadow is weak, making the images more accurate and reliable compared to winter conditions for the recognition and classification of ground features. The selection of Landsat data is shown in Table S2. Field checking of landslide inventories validated the remote sensing techniques along a ~ 30 km reach downstream from the landslide-dam, where roads allow investigation.
The following steps were used to identify and map concurrent landslides51: (1) Extraction of the landslide connected to the valley floor employed supervised classification of features using ENVI software. (2) Landslides, floodplains and roads show some similarity in remote sensing images. Independent work showed that areas with a gradient of < 20° have very low landslide densities52, and the results were re-extracted for areas with slopes > 20°. (3) Finally, visual interpretation was used to confirm or modify the landslide boundaries.
Landslide volume estimation. The landslides were identified and numbered as L1-L22 (Fig. 2). The DOD value was used to calculate the height/erosion depth of the landslide. The total volume of a landslide body is the sum of the DOD values of all grid cells multiplied by the unit grid cell area. Thus, we calculated the volume of each landslide individually as:
$$\begin{array}{c}{V}_{Lj}={\sum }_{i=1}^{n}{H}_{i}\times {A}_{cell}\#\left(9\right)\end{array}$$
where \({V}_{Lj}\) is the volume of the landslide numbered \(Lj\), with \(j\) = 1–22; \(n\) is the total number of grid cells for the landslide; \({H}_{i}\) is the elevation difference of each grid cell in the landslide area, that is, the DOD value of each grid cell; \({A}_{cell}\) represents the area of the unit grid cell (i.e., 30 m × 30 m). The vertical accuracy of the DEM data on the hillside is less than that of the valley floor (see Supplementary Text). This may result in a specific error in calculating the landslide volume using the DOD value. Empirical formulas18,32, which were summarized from a dataset comprising 4,231 bedrock and soil landslides that included data from the Namche Barwa massif, were employed to verify the landslide area-volume (Fig. 5c).
Valley Floor Width Mapping
We manually mapped the regional valley floor width using Google Earth images and Landsat remote sensing images. For each image, we manually mapped the margins of the valley bottom, defined as the break in slope at the upper edge of the (near) vertical walls of the gorge, mostly marked by the absence of vegetation. The valley floor was divided into 53 segments with an average length of 1.5 km. The total area of the valley floor divided by the channel centerline length before and after the flood is the valley floor width of each segment.
Erosional Impact Of 2000 Yigong Outburst Flood Compared With Long-term Denudation
The study area lies in the central part of the Namche Barwa–Gyala Peri (NBGP) massif, where a close coupling between tectonics and climate governs surface processes, leading to erosion through river incision and landslides. Assessment of the 2000 Yigong outburst flood is based on the volume and area of alluvial sediment and bedrock eroded across the valley floor during the 2000 outburst flood, combined with the volume and area of concurrent landslides on the adjoining hillslopes. For the valley floor, we calculated the difference between digital elevation maps before and after the flood (DOD) along the ~ 80 km reach downstream from the dam breach. The assessment takes into account erosion and deposition across the whole valley floor affected by the flood, yielding a mean erosion depth. For concurrent landslides above the valley floor, we used the landslide erosion volume divided by the landslide area, yielding a mean erosion depth. Based on their respective areas along the reach, the depths were integrated to calculate the relevant average erosion depth for the 2000 outburst flood in the study area. The study reach lies directly below the breached landslide dam. Consequently, we infer that the bulk of the deposited material originated in the study reach, allowing a robust calculation of erosion and deposition in the reach. However, some may have been swept into the reach from the lake bed and margin upstream.
Several methods have yielded denudation rates over different timescales for the NBGP massif, and they are generally in agreement with each other, despite some discrepancies. Thermochronological methods13,15,16 provide estimates averaged s over ~ 106 years and have yielded a denudation rate of 7 ± 2 mm/yr (5–9 mm/yr13,15,16). Cosmogenic radionuclide methods14 provide estimates averaged over ~ 103 years and have yielded a denudation rate of 4–28 mm/yr14. Suspended sediment load of rivers within the catchment12,13 provides estimates for annual to decadal timescales and has yielded an average annual denudation rate for the 1971–1979 period of 10 mm/yr (5–17 mm/yr)12,13. The denudation rates obtained by the different methods have limitations. Modern suspended sediment loads may be affected by enhanced sediment flux due to anthropogenic effects, the stochastic nature of precipitation, and by the fact that only the suspended sediment load is used and not the total sediment load. Cosmogenic radionuclide data have a large error range due to active landslides along the Tsangpo Gorge14, and the time scale of the thermochronological methods is long. Nevertheless, the values accord broadly.
The present study considers erosion over the duration of the Quaternary. Therefore, in order to calculate how many years of regular meteorological processes might be needed to generate an equivalent amount of erosion, we apply denudation rates based on thermochronological and cosmogenic data. The relevant average erosion depth for the 2000 flood was divided by a range of 5 mm/yr to10 mm/yr.