Amplified Risks of the Yarlung Zangbo–Brahmaputra River to Glacier Hazard Chains due to Multi-Hazard Transformation

Ruochen Jiang; Limin Zhang; Ming Peng; Wenjun Lu; Dalei Peng; Shihao Xiao; Xin He

doi:10.1016/j.eng.2025.07.002

Engineering ›› 2025, Vol. 53 ›› Issue (10) :187 -202. DOI: 10.1016/j.eng.2025.07.002

Research

Research-article

Amplified Risks of the Yarlung Zangbo–Brahmaputra River to Glacier Hazard Chains due to Multi-Hazard Transformation

Author information +

History +

PDF (15668KB)

Abstract

A glacier hazard chain can form a long-runout mass flow and generate a large flood, affecting downstream areas hundreds of kilometers away from the initiating hazard site. This study focuses on the Yarlung Zangbo Daxiagu. The objective is to address two key unresolved issues: the evolution of detached glacier materials into debris flows or debris floods and the amplification of the impact range and threats. A comprehensive framework is developed that considers the impacts of near-field and far-field hazards. Numerical modeling, remote sensing, and field investigations were integrated to understand the interactions, transformations, and amplifications of hazards in the glacier hazard chain. The results indicate that extensive, nearly saturated sediments on the glacier valley floor, when entrained, amplify the magnitude of the mass flow. The topography plays a crucial role. When the valley outlet is perpendicular to the river course, topographic obstacles cause immediate halting, resulting in the formation of high barrier dams. Conversely, when the glacier valley aligns nearly parallel to the river course, the mass flow can travel a much longer distance upon entering the river, causing an enlarged affected area. The barrier dams can breach rapidly, causing breaching floods that amplify the downstream impact from several kilometers to hundreds of kilometers. Our analysis reveals that the overall impacts remain spatially limited. Specifically, downstream areas along the Yarlung Zangbo–Brahmaputra River are unlikely to face greater threats from the upstream floods than local monsoon floods. Our findings provide the foundation for the management of glacier hazard chains.

Graphical abstract

Keywords

Glacier hazard chain / Multi-hazard transformation / Risk amplification / Mass flow / River damming / Flood risk / Yarlung Zangbo–Brahmaputra River

Cite this article

Download citation ▾

Ruochen Jiang, Limin Zhang, Ming Peng, Wenjun Lu, Dalei Peng, Shihao Xiao, Xin He. Amplified Risks of the Yarlung Zangbo–Brahmaputra River to Glacier Hazard Chains due to Multi-Hazard Transformation. Engineering, 2025, 53(10): 187-202 DOI:10.1016/j.eng.2025.07.002

登录浏览全文

4963

注册一个新账户忘记密码

1. Introduction

High Mountain Asia (HMA) has extensive glacier coverage and exhibits the greatest topographic variations on Earth [1,2]. These distinct characteristics provide rich water resources for neighboring nations and stimulate active hydropower developments, where more than 650 hydropower projects are under construction or planned [[3], [4], [5]]. However, with a temperature increase rate (0.3 °C per decade) that is twice the global average (0.16 °C per decade) over the past few decades [[6], [7], [8]], HMA is highly sensitive to hazard-triggering climate conditions. Consequently, frequent glacier hazards have occurred in recent years [9], such as glacier detachments in Sedongpu Glacier [10,11] and Aru Glacier [11,12], each releasing more than 1 × 10⁸ m³ of sediment. Such hazards occurred in a chain of events and destroyed infrastructure in the form of either rock–ice avalanches, as in the two Aru events in 2016 [11,12], or river damming and long-runout outburst debris floods (> 400 km), as in the 2000 Yigong event and the 2018 Sedongpu hazards [12,13] (Fig. 1).

Catastrophic consequences have occurred, such as tragic loss of life. For example, a devastating incident caused by intense rainfall and snowmelt, the 2013 Kedarnath hazard in Uttarakhand of India, resulted in the loss of more than 6000 lives and the destruction of at least ten hydropower stations [9]. A rock–ice hazard chain in India’s Chamoli district in 2021 resulted in extensive damage to two hydropower stations and claimed over 200 lives [14,15]. These examples serve as stark reminders of immense threats and highlight the urgent need for proactive measures to mitigate glacier hazards in HMA. In particular, there is growing concern about the safety of river systems and hydropower development due to the immense magnitude of these hazards and their cascading effects [9,12]. A critical scientific challenge is to understand how complex transformations of multiple hazard types (e.g., debris flows and breaching floods), known as multi-hazard transformations, occur and how they amplify the exposure risk within a hazard chain process [2,9,16].

Data-driven approaches, including statistical and machine learning models, have been utilized to analyze individual hazards by leveraging historical hazard inventories. These approaches utilize large datasets to identify patterns, correlations, and trends of hazardous events. Statistical models provide insights into the factors influencing hazard occurrence, whereas machine learning algorithms are capable of extracting complex relationships to predict future hazards [[17], [18], [19]]. However, this strategy tends to underestimate the overall impacts because of insufficient consideration of any cascading effects. To achieve a more accurate assessment, it is imperative to consider the hazard chain process and physical interactions [17,[20], [21], [22]].

Several physical models, such as EDDA [23,24], TOCHNOG [25], DAN3D [26], and RAMMS [27], have been developed and employed to reconstruct the dynamic process of mass flows and assess potential influencing factors. Mergili et al [28] developed an advanced, comprehensive, and efficient computational tool (r.avaflow), which has been used in real case analyses [29]. The propagation and impact amplification of such mass flows are associated with moving geomaterials and the basal bed [30]. When large-area calculations are conducted under multiple hazard scenarios, mass flows are commonly treated as an equivalent single-phase fluid mixture comprising solid and water. This approach allows the analysis of the dynamic characteristics of mass flow, as well as the distribution of erosion and deposition, ultimately providing insights into the extent of the disaster impact range. This strategy has been applied in case studies in the mountain cryosphere [15,30]. Pudasaini [31] and Pudasaini and Mergili [32] constructed a multi-phase simulation framework for mass flow movements to conduct a more comprehensive analysis of the interactions between different phases of materials during mass flow movement, the interactions between mass flow and channels, and the consequential impact on the evolution of dynamic processes [32,33]. This multi-phase flow analysis method has already been applied to the analysis of cascading hazards in the Piz Cengalo–Bondo, Huascarán, Zelongnong, and Chamoli events [14,29,[34], [35], [36]].

However, there is also a notable research gap in understanding the relationship between the transformation of hazard types and the topographic characteristics of glacier valleys and mountainous rivers. There is an urgent need for a comprehensive framework to assess the impacts of hazards and the amplification of risks. To bridge this gap, we combine numerical modeling, remote sensing, and field investigations to investigate the interactions, transformations, and amplifications of hazards within a glacier hazard chain and propose a framework to assess the amplified risks of mountain river systems exposed to glacier cascading hazards. This study focuses on the lower reach of the Yarlung Zangbo–Brahmaputra River in eastern Himalaya (Figs. 1(a)–(d)), a region characterized by a widespread glacier distribution [37,38]. The outcomes of this study will contribute to the development of effective risk mitigation strategies for ensuring water security at the Asia Water Tower.

2. Research area

The research area is located in the lower reach of the Yarlung Zangbo–Brahmaputra River, eastern Himalaya (Fig. 1). The river has annual mean runoff rates of 1.654 × 10¹¹ m³ at the outlet of Yarlung Zangbo Daxiagu and 5.372 × 10¹¹ m³ at Bahadurabad [39]. This river sharply bends in Yarlung Zangbo Daxiagu, passing between the peaks of Gyala Peri (7294 m a.s.l.) and Namcha Barwa (7782 m a.s.l.), with an average depth of 2268 m and a maximum depth of 6009 m [10,39]. Polythermal and temperate glaciers dominate the landscape, covering 7400 km² [40,41]. Historically, frequent glacier cascading hazards struck this region, such as the 1950 and 1964 Zelongnong disasters (Glacier No. 9) [11,42], the 2000 Yigong event [13], the 2018 Sedongpu event (Glacier No. 11) [10,11], and several hazards in the Peilong glaciers during 1980–2013 (Glacier No. 14) [43,44].

3. Methodology

We integrate numerical simulation, remote sensing, and field investigations to evaluate glacier hazard chains and assess the potential exposure risks in this area (Fig. 2). Table 1 summarizes the data sources used in this study. Multi-source satellite images were taken to delineate geomorphological characteristics (Fig. S1 in Appendix A) and detect glacier movements using an offset tracking technique (Figs. S2 and S3 in Appendix A). Stereo satellite data is used to produce high-resolution ground elevation data and analyze ground elevation changes before and after selected historical hazards (Fig. S1 and Section S1 in Appendix A). The geographical information, including the terrain and river system, was processed in ArcGIS Pro software (Esri Company, USA) to represent regional topographic features. Long-term regional climate and meteorological data were collected at Nyingchi (29.57°N, 94.47°E, 3001 m a.s.l.) and Bomi (29.87°N, 95.77°E, 2735 m a.s.l.) (Fig. 1(d)). These historical cases serve as benchmarks to help refine the hazard chain analysis in this region (Fig. S4 in Appendix A). We have also made several field trips to collect samples, measure in situ physical properties and investigate historical hazards (Fig. 3, Fig. 4, and Figs. S5 and S6 in Appendix A). The details are presented as follows.

3.1. Remote sensing analysis of glaciers and historical hazards

The majority of valley glaciers in the study area are located in Yarlung Zangbo Daxiagu. Owing to their remote and inaccessible locations, comprehensive and efficient monitoring of glaciers and their associated hazards relies heavily on remote sensing techniques. By conducting a comparative analysis of the study region at different time intervals via multi-source satellite images, we quantified the surface deformation of the ground surface using glacier image velocimetry (GIV) [45]. The fundamental idea of this strategy is to measure the movement velocity of glaciers by tracking the displacement characteristics of image pixels of satellite images over the research area. The analysis is conducted using co-registered satellite images, and the properties of two images are compared to determine the most suitable alignment between them. If any motion has occurred, the corresponding pixel features exhibit displacement relative to their original locations within the reference satellite image. This displacement signifies the movement or shifts that has taken place in the period between the two captured images. By analyzing and quantifying these pixel displacements, the magnitude and direction of motion can be obtained.

On the basis of the obtained glacier surface motion and slope maps of the research area (Fig. S2 and Section S2 in Appendix A), the thickness of the glacier sediments was calculated using a glacier flow model. The model employed in this study adopts the shallow-ice approximation, assuming that glacier motion is predominantly governed by basal sliding and sediment rheology. The rheological behavior is governed by Glen’s law [1]. The glacier thickness (H) is expressed as Eq. (1):

(1)

H = {(\frac{v_{s} (1 - β) (n_{g} + 1)}{2 C_{g} {(ρ_{g} g)}^{n_{g}} {‖ \nabla H_{s} ‖}^{n_{g}}})}^{n_{g} + 1}

where

\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y})

is the Hamiltonian operator; v_s is the glacier surface velocity and

β

is a glacier flow parameter determined following Ref. [1]; n_g is an exponent and taken as 3; H_s is the surface elevation;

ρ_{g}

is the valley sediment density; g is the gravitational acceleration; and C_g is a creep parameter taken as 2.4 × 10⁻²⁴ Pa·s⁻¹.

3.2. Field investigations

We conducted several field trips and in situ tests in southeastern Xizang to collect samples, measure in situ physical properties and investigate historical hazards. Figs. 3, S3, and S5 show several historical hazards and sampling locations, respectively. The mass flows formed by detached sediments cause channel erosion, where nearly saturated deposits are distributed (Figs. S5 and S6). The particle size distribution curves and physical properties of the samples are summarized in Fig. 4 and Table 2, respectively.

3.3. Physically based simulation program

The dynamics of the detached mass flow and debris flood are analyzed and quantified via the EDDA simulation program [23,24]. The flow dynamics are described via the following equations:

(2)

\frac{\partial h}{\partial t} + \frac{\partial (h v_{x})}{\partial x} + \frac{\partial (h v_{y})}{\partial y} = i [C_{V^{*}} + (1 - C_{V^{*}}) s_{b}]

(3)

\frac{\partial (C_{V} h)}{\partial t} + \frac{\partial (C_{V} h v_{x})}{\partial x} + \frac{\partial (C_{V} h v_{y})}{\partial y} = i C_{V^{*}}

(4)

\frac{\partial (h v_{x})}{\partial t} + \frac{\partial (h {v_{x}}^{2} + 0.5 g h^{2})}{\partial x} + \frac{\partial (h v_{x} v_{y})}{\partial y} = g h [- sgn (v_{x}) S_{f x} - \frac{\partial H_{b}}{\partial x}] - v_{x} \{i [C_{V^{*}} + (1 - C_{V^{*}}) s_{b}]\}

(5)

\frac{\partial (h v_{y})}{\partial t} + \frac{\partial (h v_{x} v_{y})}{\partial x} + \frac{\partial (h {v_{y}}^{2} + 0.5 g h^{2})}{\partial y} = g h [- sgn (v_{y}) S_{f y} - \frac{\partial H_{b}}{\partial y}] - v_{y} \{i [C_{V^{*}} + (1 - C_{V^{*}}) s_{b}]\}

where t is time; h denotes the flow depth; v_x and v_y represent the x- and y-flow velocities, respectively;

i = \partial H_{b} / \partial t

signifies the erosion or deposition rate; H_b is the bed elevation; s_b is the saturation degree in the erodible bed; sgn(·) represents a sign function that follows the conditions: sgn(x) = 1 when x > 0, sgn(x) = 0 when x = 0, and sgn(x) = −1 when x < 0; S_f_x and S_f_y indicate the flow resistance slopes in the x- and y-directions; C_V* is the volume fraction of solid in the erodible bed; and

C_{V}

is the volumetric sediment concentration of the mixture. Eqs. (4), (5) simplify the erosion-induced momentum change by representing it as a reduced production term [34].

The flow mixture can be considered a water flow (

C_{V} < 0.2

) with the same dynamic viscosity as water, a hyper-concentrated flow (

0.2 {< C}_{V} < 0.45

) with increased yield stress and dynamic viscosity, or a full debris flow or mass flow (

C_{V} > 0.45

) with substantially increased yield stress and dynamic viscosity. For

C_{V} < 0.2

S_{f}

is the total flow resistance slope and can be obtained using the Manning equation, and the following quadratic rheological model is employed when

C_{V} > 0.2

[46]:

(6)

S_{f} = \frac{τ_{y}}{ρ g h} + \frac{K μ v}{8 ρ g h^{2}} + \frac{n_{t d}^{2} v^{2}}{h^{4 / 3}}

where

τ_{y}

denotes the yield stress of the mass flow; v is the absolute value of the depth-averaged flow velocity;

ρ

is the density of the mass flow; K and

μ

are the laminar flow resistance and the dynamic viscosity, respectively; and n_td represents the equivalent Manning coefficient of the mass flow:

(7)

n_{t d} = 0.0538 n_{m} e^{6.0896 C_{V}}

where

n_{m}

is the Manning coefficient.

τ_{y}

and

μ

are estimated as [47]:

(8)

τ_{y} = α_{1} e^{β_{1} C_{V}}

(9)

μ = α_{2} e^{β_{2} C_{V}}

where

α_{1}

α_{2}

β_{1}

, and

β_{2}

are empirical parameters. Eqs. (8), (9) are taken when

C_{V}

> 0.1. When

C_{V}

is no more than 0.1,

τ_{y}

is determined using:

(10)

τ_{y} = (1 - C_{s}) C_{V} (ρ_{s} - ρ_{w}) g h \cos θ^{2} \tan ϕ

where

C_{s}

is the coefficient of suspension of solid particles;

θ

is the slope angle;

ρ_{s}

and

ρ_{w}

are the densities of solids and water, respectively; and

ϕ

is the friction angle of the solid particles.

Erosion takes place when the bed shear stress exceeds the critical erosive shear stress and when the volumetric sediment concentration (i.e.,

C_{V}

) is lower than a specific equilibrium value,

C_{V \infty}

[48]:

(11)

C_{V \infty} = \frac{ρ_{w} \tan θ}{(ρ_{s} - ρ_{w}) (\tan ϕ_{b e d} - \tan θ)}

where

ϕ_{b e d}

is the internal friction angle of the erodible bed. The erosion rate is calculated as follows [49]:

(12)

i = K_{e} (τ - τ_{c})

where

K_{e}

and

τ_{c}

are the erosion coefficient and the critical erosive shear stress, respectively; the shear stress (

τ

) can be expressed as [50]:

(13)

τ = ρ g h S_{f}

The erosion coefficient

K_{e}

and the critical erosive shear stress

τ_{c}

are calculated using the following equations [51]:

(14)

K_{e} = 20 075 {e_{b}}^{4.77} C_{u}^{- 0.76}

(15)

τ_{c} = \frac{2}{3} g d_{50} (ρ_{s} - ρ_{w}) \tan ϕ

where

e_{b}

denotes the void ratio of the bed materials and C_u = d₅₀/d₁₀ is the coefficient of uniformity, in which d₅₀ and d₁₀ represent the particle diameters at 50% and 10% finer by weight, respectively.

The deposition process of solid material occurs when

C_{V}

C_{V \infty}

and v < v_e, and the deposition rate can be expressed as [52]:

(16)

i = C_{d e p} (1 - \frac{v}{p v_{e}}) \frac{C_{V \infty} - C_{V}}{C_{V *}} \times v

(17)

v_{e} = \frac{2}{5 d_{50}} (\frac{g \sin θ_{e} ρ}{0.02 ρ_{s}}) λ^{- 1} h^{1.5}

(18)

θ_{e} = \arctan (\frac{C_{V} (ρ_{s} - ρ_{w}) \tan ϕ_{b e d}}{C_{V} (ρ_{s} - ρ_{w}) + ρ_{w}})

(19)

λ^{- 1} = {(\frac{C_{V^{*}}}{C_{V}})}^{1 / 3} - 1

where C_dep denotes the deposition rate coefficient and is calibrated via benchmark analysis (Section 5.1);

λ

is a coefficient related to values of

C_{V}

and

C_{V *}

, which can be calculated using Eq. (19).

θ_{e}

is the critical angle;

v_{e}

is the critical flow velocity; and p is an empirical coefficient, and taken as 0.67 according to Ref. [48].

When river damming occurs, the dam breaching process is simulated as an erosion process of the dam materials by the overtopping river flow [52]. The erosion rate of dam materials in water flow is quantified using Eqs. (12), (13), (14), (15), and the shear stress

τ_{d}

can be calculated as:

(20)

τ_{d} = γ_{w} R_{h} S

where

γ_{w}

is the unit weight of water;

R_{h}

is the hydraulic radius; and S is the energy slope and is equal to the slope gradient of the breaching channel. The cross-section of the barrier dam is assumed to be a trapezoid, so

R_{h}

can be calculated as:

(21)

R_{h} = \frac{(H_{w} - H_{bb}) \cos φ + W_{b} \sin φ}{2 (H_{w} - H_{bb}) + W_{b} \sin φ} (H_{w} - H_{bb})

where

H_{w}

is the elevation of the water surface; H_bb and

W_{b}

are the elevation and width of the breach bottom, respectively; and

φ

is the angle of the side slope. The breaching discharge

Q_{b}

is calculated as [53]:

(22)

Q_{b} = 1.7 [W_{b} + (H_{w} - H_{bb}) \tan φ] {(H_{w} - H_{bb})}^{3 / 2}

The hydrodynamic process of dam breaching is obtained by solving the following equation:

(23)

A_{s} \frac{d H}{d t} = Q_{i n} - Q_{o u t}

where

A_{s}

is the surface area of the barrier lake and

Q_{i n}

and

Q_{o u t}

are the water inflow and outflow through the dam crest, respectively. More details about the dam breaching model are available in the Ref. [52].

3.4. Assessment framework of multi-hazard exposure risks

A comprehensive evaluation of multi-hazard exposure risk is imperative, considering that a single river stretch may experience landslides, debris flows, or outburst floods triggered by the instability of glaciers. In this study, we propose a framework to assess exposure risk (Fig. 5) using three indicators, including the equivalent dam height (EDH), equivalent peak discharge (EPD), and equivalent sediment thickness (EST), on the basis of multi-scenario analysis:

(24)

E D H = \sum_{k = 1}^{4} \sum_{j = 1}^{N_{k}} D H_{jk} P_{k}

(25)

E P D = {m a x}_{j = 1, 2, \dots, N_{k}} (\sum_{k = 1}^{4} {P D}_{jk} P_{k})

(26)

E S T = {m a x}_{j = 1, 2, \dots, N_{k}} (\sum_{k = 1}^{4} {S T}_{jk} P_{k})

where P_k and N_k indicate the occurrence frequency of scenario k and the number of hazard events, respectively; DH_jk, PD_jk, and ST_jk are the dam height, peak discharge, and flood sediment thickness caused by hazard event j in scenario k, respectively. The determination of multiple calculation scenarios is based on the observed characteristics of mountain glaciers, which will be explained in detail in Section 4.

4. Glacier distribution and likely detachment magnitudes

Fig. 6 shows the obtained thicknesses of the glacial deposits. The associated error analysis is presented in Section S2. Our remote sensing results show that most valley glaciers are located on gentle valley floors (9°–18°) in the Yarlung Zangbo Daxiagu along the Yarlung Zangbo–Brahmaputra River (i.e., Paizhen–Medog reach) (Fig. 6(a)). Glacial sediments are distributed extensively, with an average thickness of 200 m and a maximum thickness of 550 m (Fig. 6(a)). These valley glaciers have accelerated markedly at velocities four times greater than they did ten years ago (Fig. 6 (b) and Fig. S2). The massive valley glaciers pose significant threats to mountain rivers because of their combination of rapid movement velocities, large volumes, and close distances to the river course (Figs. 6(c)–(e)). The morphometric characteristics of these glaciers, including length and volume distributions, are systematically summarized in Fig. 6(f). Local meteorological observations in Nyingchi (29.57°N, 94.47°E, 3001 m a.s.l.) and Bomi (29.87°N, 95.77°E, 2735 m a.s.l.) indicate a regional warming trend of 0.25–0.35 °C per decade, which is twice the global average (0.16 °C per decade) observed over the past decades (Fig. 7). As presented in Fig. 6, Fig. 7, climatic warming constitutes a primary driver of accelerated glacier movements in the study region.

Previous studies have overlooked a crucial aspect: the determination of hazard scenarios that account for the varying volumes of detached glaciers. This information is indispensable for accurate physical modeling and risk analysis of glacier hazard chains. By incorporating diverse detached glacier volumes into the analysis, a better understanding of potential hazards can be attained, allowing for more robust assessments of the associated risks. Considering that glacier detachment typically occurs at the glacier front, we employed the glacier detachment length (L_D) as an indicator to assess the magnitude of initial hazards. Fig. 6(f) shows that the L_D and glacier detachment volume (V_D) follow an exponential relationship (i.e., L_D and log₁₀(V_D) satisfy a linear relationship with R² = 0.78). The length (L_G) and volume (V_G) of the valley glaciers in this region have the same relationship, with R² = 0.80, and V_G is close to the detached volume V_D when the L_D approaches the L_G. The fitting results in Fig. 6(f) provide a valuable tool for predicting and assessing the magnitude of glacial instability based on a measurable parameter.

Furthermore, the occurrence frequency of glacier detachments of various magnitudes under these scenarios can be assumed to be the occurrence frequency of historical events (Fig. 6(f)). We define four hazard scenarios, along with their corresponding magnitude and occurrence frequency of glacier detachments, to assess the impacts of their subsequent cascading hazards: ① Scenario 1, with low-magnitude and L_D = 1 km (occurrence probability: 54%); ② Scenario 2, with moderate-magnitude and L_D = 2 km (occurrence probability: 23%); ③ Scenario 3, with high-magnitude and L_D = 3 km (occurrence probability: 15%); and ④ Scenario 4, an extreme scenario with the destabilization of all sediments from each valley glacier (occurrence probability: 8%).

5. Hazard propagation and transformation

5.1. Benchmark analysis based on historical cases

In the research area, the 2018 Sedonpgu glacier detachment in Yarlung Zangbo Daxiagu (Fig. 8) is taken as a benchmark to validate the parameters for the simulation of all valley glaciers. This hazard involves two glacier detachments, occurring at Sedongpu at 10:48 PM Beijing time on October 16, 2018, and at 8:03 AM Beijing time on October 29, 2018, according to monitoring information from local seismic and hydrological stations. The glacier detachment area and ground elevation changes were derived using ZY-3 stereo satellite data (Fig. 8). The volume of the detached mass is 1.5 × 10⁸ m³, and the average and maximum depths in the detachment zone (P1–P2, 3100–4100 m a.s.l.) are 150 and 190 m, respectively. P2–P4 represents the deposition zone, which includes a valley deposition area (P2–P3, 2750–3100 m a.s.l.) and a river deposition area (P3–P4) (2750 m a.s.l.). The barrier dam, with a volume of V_Dam = 5.5 × 10⁷ m³, blocked the river course and triggered an outburst flood that traveled along the 460 km canyon.

EDDA is applied to quantify the runout process and deposition distribution in the Sedongpu case. The discharge of river flow along the Yarlung Zangbo–Brahmaputra River is set as 2000 m³·s⁻¹ according to hydrologic data [10]. The geo-material properties and physical parameters are determined based on field investigations and geotechnical tests. The deposition rate coefficient, C_dep, is evaluated to ensure that the simulated deposition results are close to the ground elevation derived from the ZY-3 stereo images. As shown in Fig. 8, C_dep = 0.03 provides the best agreement between the simulated and measured results, which is then employed in the multi-scenario simulations for other valley glaciers in the Yarlung Zangbo Daxiagu.

5.2. Multi-scenario hazard chain analysis

Numerous glaciers are situated in remote and inaccessible regions, making human access and soil sampling impractical. Notably, the limited soil samples obtained presented limited variations in their physical properties. Hence, we utilize the parameters listed in Table 3 for analysis. Furthermore, our analysis focuses on terrain factors and their impact on the movement of the formed mass flows and river damming. To ensure consistency, we employed the same set of physical parameters throughout the analysis. Multiple simulation scenarios were determined and introduced in 3.4 Assessment framework of multi-hazard exposure risks, 4 Glacier distribution and likely detachment magnitudes, with an illustrative example provided in Section S3 in Appendix A.

Fig. 9, Fig. 10, Fig. 11, Fig. 12 and Figs. S5–S8 in Appendix A show the simulated multi-scenario dynamics and evolution of hazard chains. As the volume of the unstable initial mass increases, the subsequent mass flow exhibits higher velocities (Fig. 9(a)). Consequently, the mass flow acquires a greater momentum, which intensifies the basal erosion forces acting on the channel deposits. This results in a more significant level of erosion and entrainment (Eqs. (3), (13), and (14), and Fig. 9(b)), potentially causing amplification of the impact range. Our site photographs show that the valley channel beds are often saturated with sediments (Fig. 3). The entrainment of such sediments causes changes in the solid content of the moving mixture. An increased amount of detached mass enhances channel entrainment, with observed entrainment mass ratio values (i.e., the ratio of the entrainment volume to the detached volume) in the range of 3% to 14%. These two values are close to those in the historical cases of the 2013 and 2015 Flat Creek events in Alaska (61.66°N, 141.54°W and 61.65°N, 141.55°W) (17% and 3.7%, respectively) [54] and the two Aru Glacier detachments in western Himalayas (34.04°N, 82.26°E and 34.02°N, 82.28°E) (5.8% and 21.6%, respectively) [12]. Furthermore, at the same solid volumetric ratio, a large-magnitude mass flow is associated with a lower flow resistance slope (S_f) because of its greater flow depth. Consequently, large-magnitude detachments trigger more rapid and longer-runout mass flows (Fig. 9).

With the glacier front closer to the river course (Fig. 6(e)) contributes to more river deposits in the river course (e.g., Glacier IDs 4–9 in the Paizhen–Ganglang reach), even if the initial hazard magnitude is the same, as only a small amount of mass flow material is deposited on the short flow path (Figs. S7–S10 in Appendix A). On the other hand, when the distance between the glacier front and the river channel is large (> 30 km), even in the event of a large-scale glacial collapse, the resulting debris flow cannot reach the Yarlung Zangbo Jiang, such as Glacier ID 1 (Fig. 6(c) and Figs. S7–S10). The low-magnitude detachment hazard (Scenario 1) causes relatively small-scale barrier dams 10⁶–10⁷ m³ in volume (green points in Fig. 10(a)) and 15–20 m in dam height (green points in Fig. 10(b)) when the distance is less than 5 km (L_f < 5 km), whereas moderate- and high-magnitude detachment hazards (Scenarios 2 and 3) can extend the impact range to L_f = 8–10 km. It is likely that barrier dams over 50 m high formed if glacier detachments occurred at the same scale as those in Sedongpu in 2018 (Glacier ID 11) (Fig. 10(b)).

Local topographic characteristics at the confluences of glacier valley outlets and the river course also affect the geometric features of the dam. Figs. 11(a)–(h) show longitudinal profiles of river deposits under a variety of hazard intensities in the Paizhen–Ganglang reach. Specifically, Fig. 11(d) shows the simulated river deposition and changes in riverbed elevation resulting from the Sedongpu event in 2018. The valley outlets of Glaciers IDs 4 and 11 are almost perpendicular to the river course (Figs. 12(a) and (b)). This topographic condition allows instant dissipation of kinetic energy and rapid sediment deposition when mass flows out of the valley and impacts the riverbed. A narrower river course in the river reaches of Glacier ID 4 (110 m riverbed width) results in a longer barrier dam along the river course than that of Glacier ID 11 (390 m riverbed width). When the valley outlet is almost parallel to the river (Figs. 12(c) and (d)), less kinetic energy is dissipated during the impact process, contributing to a longer runout distance along the river course and a lower dam height (Fig. 11). The variations in the entrainment–deposition process and topographic factors directly affect the subsequent breaching floods.

6. Amplified impacts from dam-breaching floods

6.1. Benchmark analysis of historical dam-breaching outburst floods

The deposition of solid materials in a mass flow does not signify the termination of the hazard chain. When massive amounts of sediment enter the river course, a barrier dam can form, obstructing the natural flow of the river. Such a barrier dam is often short-lived and can breach within a few days or one month after its formation. We simulated the benchmark dam breaching floods of the 2018 Sedongpu hazards with a peak discharge of 32 000 m³·s⁻¹ [10] and the 2000 Yigong event with a peak discharge of 120 000 m³·s⁻¹ [13] based on the flow dynamic equations of the EDDA to calibrate and validate the input parameters. The reconstruction of these events based on in situ hydrological records provides a basis for assessing the magnitudes and impacts of outburst floods in downstream areas.

Fig. 13 shows the characteristics of outburst dam-breaching floods. The simulated flood hygrographs at Sedongpu, Yigong, Dexing, Duding, and the outlet of Yarlung Zangbo Daxiagu indicate that the peak flow gradually decreases while the duration of the flood increases during the propagation of the dam-breaching flood. The discharge attenuations during the 2000 Yigong event and the 2018 Sedongpu hazard were similar. The peak discharges at the outlet of Yarlung Zangbo Daxiagu are approximately one-third of the peak flows at the dam breach locations. The in-situ records and simulation results highly agree when the Manning coefficient (

n_{m}

) is set to 0.04 (Fig. 13). The calibrated parameters are utilized in all glacier hazard analysis scenarios.

Furthermore, in 1998, the river experienced its most severe meteorological flood in recent history, caused by heavy monsoon rains. The peak flow rate reached 13 487 m³·s⁻¹ at the Nuxia hydrological station, which is 2.6 times higher than the maximum monthly mean flow rate of 5119 m³·s⁻¹ recorded in August (Fig. 14) [10,38,55]. This flood event serves as a benchmark for evaluating downstream flood threats.

6.2. Multiple-scenario risk assessment

The calculation results of the three quantitative indicators along the Yarlung Zangbo–Brahmaputra River are presented in Fig. 15. Given the detached length (L_D ≈ 3 km) and dam height (∼60 m) in Scenario 3 (Fig. 10, Fig. 11) close to the 2018 Sedongpu hazards, Scenario 3 can serve as a benchmark case for assessing potential hazards with similar impacts to those observed in the 2018 Sedongpu hazards. The results show that the peak flood discharge attenuates more slowly in the steeper sections of the river with a narrower channel (e.g., Sedongpu–Ganglang and Yigong–Ganglang) than in the downstream reaches (e.g., Ganglang–Dexing and Dexing–Beibeng) (Figs. 15(a) and (b)). The flood peak discharge attenuates to 31% ± 3.1% of the peak discharge at the outlet of Yarlung Zangbo Daxiagu, with discharges of 3 900, 6 500, 9 500, and 18 000 m³·s⁻¹ for scenarios 1, 2, 3, and 4, respectively. The large EDH values in the Paizhen–Ganglang reach indicate that this reach is more susceptible to river damming (Fig. 15(c)). In contrast, the downstream area faces significantly less river damming (Fig. 15(c)). The EPD and EST results indicate that the most intensive area of hydropower development is at the greatest risk from breaching floods and flood debris sediments (Figs. 15(c) and (d)). Figs. 15(e) and (f) further demonstrates that these evolutionary processes are coupled with the variations in the width and elevation of the river channel. Moreover, the results also indicate that a risk assessment method that is solely based on one type of hazard is inadequate and may underestimate the overall impact. For example, only the Sedongpu–Ganglang reach is hazard prone if only the EDH index is considered. Similarly, the Sedongpu–Ganglang reach and the Ganglang–Beibeng reach are found to have the same level of exposure risk if only the EPD index is used. The risk level decreases downstream along the river, but the Duding–outlet of Yarlung Zangbo Daxiagu reach is still susceptible to breaching floods (Fig. 15(g)).

Despite the potential threat of upstream outburst floods to downstream areas, a comparative analysis of historical flood records indicates that this threat is relatively less significant than the risk posed by summer monsoon floods. Fig. 14(a) show the analysis area of summer monsoon floods along the Yarlung Zangbo–Brahmaputra River. Seasonal river flow discharges are presented in Fig. 14(b), and hazard consequences caused by summer monsoon floods are shown in Figs. 14(c)–(e). It is noted that dam-breaching floods caused by high-magnitude glacier detachments in Yarlung Zangbo Daxiagu (e.g., the 2018 Sedongpu hazard) result in flow discharges almost equivalent to the monsoon floods in the outlet of Yarlung Zangbo Daxiagu, with a peak discharge of ∼10 000 m³·s⁻¹ (Figs. 14(b) and 15(b)). Even under such circumstances, the river discharge at the outlet of Yarlung Zangbo Daxiagu is less than half of that at Guwahati and Bahadurabad (Fig. 14(b)), which is associated with the summer monsoons and heavy rainfall that occur between April and September (Fig. 14(a)) [38]. These findings further support that the impact of upstream outburst floods downstream is less significant than that of monsoon floods. Monsoon flood events No. 1–3 and No. 7 in this river reach have peak charges of 75 000–80 000 m³·s⁻¹ and 90 000 m³·s⁻¹ at Bahadurabad, respectively [56]. The peak discharge of the 1998 flood at Bahadurabad reached 94 000 m³·s⁻¹ [38]. As a result, monsoon floods can cause greater impacts and affect more cropland areas (Figs. 14(c)–(e)). For example, the 1998 flood event, characterized by large discharge volumes (Fig. S11 in Appendix A), affected millions of people along the Brahmaputra downstream. Moreover, even though the magnitude of the floods at Yarlung Zangbo Daxiagu has the same initial discharges as the 1998 floods at Nuxia, with a maximum flow discharge of approximately 13 487 m³·s⁻¹, and the 2018 Sedongpu event, with a peak discharge of 32 000 m³·s⁻¹, the peak discharges at the outlet of Yarlung Zangbo Daxiagu are in the range of 10 000–20 000 m³·s⁻¹, which are below the flow discharges during the monsoon season (Fig. 14). Therefore, in the lower reaches of the Brahmaputra, such as Bahadurabad and Guwahati, outburst floods caused by glacier detachments do not appear to pose greater threats than local downstream monsoon floods.

7. Discussion

Over the past few decades, climate warming has caused frequent cascading hazards in various glaciers, such as the Kolka–Karmadon glaciers in 2002 [57], the Flat Creek glaciers in 2013 and 2015 [54], the Aru glaciers in 2016 [11,12], and the Sedongpu glaciers in 2018 [11]. Rising atmosphere temperatures undoubtedly play a key role in increasing land surface temperatures, potentially causing rapid ice ablation and weakening of geo-materials. The increasing trend in global temperature has resulted in a noticeable increase in the volume of glacial meltwater, accelerating glacier movements [12]. The deterioration of glacial material strength caused by ice melting [22], as well as the acceleration of glacial movement due to increased meltwater, has driven mountain glaciers to a state of imminent instability. Sudden external disturbances, such as heavy rainfall, intense ice melting, and earthquakes [22,23], can easily trigger large-scale glacial collapses, posing a threat to human communities and infrastructure in affected areas [[10], [11], [12]].

Our comprehensive study shows that such glacier failures trigger cascading hazards and highlights the sequential progression of cascading hazards in a typical glacier hazard chain. These hazards encompass glacier detachments, followed by mass flows/avalanches and river damming, ultimately culminating in outbursts of dam-breaching floods (Fig. 16). Traditionally, the assessment of hazard threats to a river system is often conducted by analyzing one hazard at a time. Our study demonstrates that such an approach can underestimate the inherent risks associated with glacier cascading hazards. This is because the hazard chain not only directly threatens the river reach and infrastructure in the vicinity but also, through multi-hazard transformation, expands the impact range of the initial hazard tens or hundreds of times (e.g., > 400 km) (Fig. 16). For example, the 2018 Sedongpu glacier rupture spanned approximately 3 km and resulted in a mass flow traveling 8 km and a dam-breaching flood spreading over 400 km. Therefore, a multi-hazard analysis approach must be adopted to assess exposure risk.

Our results reveal that the entrainment–deposition process and glacier valley–river terrain characteristics are two key factors in the multi-hazard transformation of glacier hazard chains. The erosion and entrainment of nearly saturated sediments, along with the deposition of a solid mass, contribute to the transformation of the initial avalanche into a debris flow or a debris flood through changes in the volume of solid content. Consequently, the flow resistance decreases, and the mobility increases. Furthermore, a larger detached glacier mass can experience less flow resistance and is therefore likely to cause more rapid and long-runout mass flow. The glacier detachments, avalanches, and debris flows have impact ranges within approximately 20 km, which are also observed in some historical cases as 8.2 and 7.2 km in the two Aru cases and 8 km in the Sedongpu event [[10], [11], [12]]. If the river course is perpendicular to the glacier valley outlet, it introduces another significant threat known as river damming, as seen in the Sedongpu hazard. The propagation of mass flow will be impeded by river damming, but the breach of the barrier dam can trigger an outburst flood, significantly amplifying the extent of hazard impact by up to 10–100 times, with runout distances greater than 400 km in the 2000 Yigong event [13] and the 2018 Sedongpu hazard [10]. Glacier hazards can thus cause serious damage downstream in the form of outburst floods (Fig. 16).

On the other hand, our numerical simulation model drastically simplifies real-world problems. Momentum exchange between the fluid and the particles and phase segregation in the flowing mass are not considered in this depth-averaged simulation model. Particle segregation, the solid-dominated front and the fluid-dominated front must be simulated via multi-phase mass flow models [31,32,58]. Multi-phase flow simulation methods can provide physical processes for analyzing the impacts of hazard chains involving interactions between solid and fluid phases. These advanced physics-based methods require governing equations for the mass and momentum of different phases naturally involved in real field events [33,34]. For multi-phase flow analysis, a corresponding unified erosion model, such as that first proposed in the Ref. [33], is necessary. Pudasaini and Krautblatter [33] pioneered a mechanical model for the energy budget of erosive mass flows, which controls enhanced or reduced mobility. This model provides the mechanical conditions for when, how, and how much energy erosive landslides gain or lose. Complex erosion-induced mobility should be further investigated in the future by employing the unified, multi-phase erosion models in previous studies [33,34].

A comprehensive framework has been proposed to assess exposure risks associated with glacier hazard chains in the Yarlung Zangbo Daxiagu. Our analysis provides an evidence-based reference for hazard risk assessment in the Yarlung Zangbo–Brahmaputra basin. First, our results suggest that the river reach (i.e., Sedongpu–Ganglang) in Yarlung Zangbo Daxiagu is at the highest risk. Second, we clarify a crucial concern regarding whether glacier geohazards in the Yarlung Zangbo Daxiagu could result in catastrophic flooding hazards that threaten downstream regions. Our findings clarify that while outburst floods caused by glacial detachments occur in Yarlung Zangbo Daxiagu, outburst floods attenuate quickly along the river so that their impacts decrease in downstream areas, such as Bahadurabad and Guwahati. The outburst floods resulting from glacial detachments do not pose a greater threat than the local downstream floods during the monsoon season.

8. Conclusions

Systematically assessing the threats of glacier hazard chains is essential for public and infrastructure safety in alpine gorges. Our research reveals significantly amplified impact space and exposure risks of hazard chains. We propose a framework to assess the exposure risks associated with glacier hazard chains. The main findings are summarized as follows.

(1) This paper defines several scenarios of exposure risks of mountain rivers to glacier hazard chains based on glacier distribution, thickness, motion, and historical hazards. We propose the use of potential detached length L_D as a hazard indicator. This indicator can be used to estimate the volume of detached glacier mass using the fitting equation in Fig. 6(f), even without the use of remote sensing inversion or other means to obtain the glacier thickness.

(2) The exposure risks to mountain rivers posed by glacier hazard chains are greatly intensified due to multi-hazard transformations. Two key factors for this amplification are the entrainment–deposition of mass flows and the orientation of the valley outlet to the river course. If the river course is perpendicular to the glacier valley outlet, another notable threat known as river damming may arise. The extent of impact can be significantly amplified due to dam-breaching floods. Furthermore, multiphase simulation models are recommended for a detailed analysis of physical processes such as phase segregation, phase change, and phase interactions, as well as interactions between the flowing mass and infrastructures.

(3) We propose a framework for evaluating the exposure risks associated with glacier hazard chains, which allows us to identify potential threat hotspots along the Yarlung Zangbo Daxiagu. The framework assesses the risk level by considering the impacts of both near-field and far-field disaster events on glacier hazard chains. The analysis in Yarlung Zangbo Daxiagu demonstrated that, despite the presence of damming and flooding threats in the Paizhen–Ganglang reach of Yarlung Zangbo Daxiagu, the downstream areas along the Yarlung Zangbo–Brahmaputra River, such as the outlet of Yarlung Zangbo Daxiagu area, are unlikely to face greater risks than those from local monsoon-driven floods. The findings of our study provide the foundation for the effective management of risks associated with glacier hazard chains and the safety of the river systems in this area.

CRediT authorship contribution statement

Ruochen Jiang: Validation, Investigation, Conceptualization, Writing – original draft, Methodology, Formal analysis. Limin Zhang: Writing – review & editing, Project administration, Funding acquisition, Conceptualization, Supervision, Investigation, Data curation. Ming Peng: Writing – review & editing, Funding acquisition, Investigation, Data curation. Wenjun Lu: Writing – review & editing, Data curation, Investigation, Conceptualization. Dalei Peng: Visualization, Data curation. Shihao Xiao: Writing – review & editing, Data curation. Xin He: Data curation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work obtained support from the National Natural Science Foundation of China (U20A20112, 42061160480, 42377196, and 52479095), the NSFC/RGC Joint Research Scheme (42061160480 and N_HKUST620/20), the Research Grants Council of the Hong Kong SAR Government (16203720, T22-606/23-R, and JRFS2526-6S09) and the Project of Hetao Shenzhen–Hong Kong Science and Technology Innovation Cooperation Zone (HZQB-KCZYB-2020083).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2025.07.002.

References

Publishing order | Descend order by publishing year | Descend order by cited within

[1]	Millan R, Mouginot J, Rabatel A, Morlighem M.Ice velocity and thickness of the world’s glaciers.Nat Geosci 2022; 15(2):124-129.

[2]	Mani P, Allen S, Evans SG, Kargel JS, Mergili M, Petrakov D, et al.Geomorphic process chains in high-mountain regions—a review and classification approach for natural hazards assessment.Rev Geophys 2023; 61:e2022RG000791.

[3]	Marzeion B, Cogley JG, Richer K, Parkes D.Attribution of global glacier mass loss to anthropogenic and natural causes.Science 2014; 345(6199):919-921.

[4]	Jiang R, Dai F, Liu Y, Li A.Fast marching method for microseismic source location in cavern-containing rockmass: performance analysis and engineering application.Engineering 2021; 7(7):1023-1034.

[5]	Bao Z, Zhang J, Lian Y, Wang G, Jin J, Ning Z, et al.Changes in headwater streamflow from impacts of climate change in the Tibetan Plateau.Engineering 2024; 34:133-142.

[6]	Hugonnet R, Mcnabb R, Berthier E, Menounos B, Nuth C, Girod L, et al.Accelerated global glacier mass loss in the early twenty-first century.Nature 2021; 592(7856):726-731.

[7]	Moore JC.Glaciology and global climate change.Engineering 2018; 4(1):6-8.

[8]	Lu C, Cai C.Challenges and countermeasures for construction safety during the Sichuan–Tibet Railway Project.Engineering 2019; 5(5):833-838.

[9]	Li D, Lu X, Walling DE, Zhang T, Steiner JF, Wasson RJ, et al.High Mountain Asia hydropower systems threatened by climate-driven landscape instability.Nat Geosci 2022; 15(7):520-530.

[10]	Chen C, Zhang LM, Xiao T, He J.Barrier lake bursting and flood routing in the Yarlung Tsangpo Grand Canyon in October 2018.J Hydrol 2020; 583:124603.

[11]	Kääb A, Jacquemart M, Gilbert A, Leinss S, Girod L, Huggel C, et al.Sudden large-volume detachments of low-angle mountain glaciers—more frequent than thought?.Cryosphere 2021; 15(4):1751-1785.

[12]	Kääb A, Leinss S, Gilbert A, Bühler Y, Gascoin S, Evans SG, et al.Massive collapse of two glaciers in western Tibet in 2016 after surge-like instability.Nat Geosci 2018; 11(2):114-120.

[13]	Xu Q, Shang YJ, van T Asch, Wang ST, Zhang ZY, Dong XJ.Observations from the large, rapid Yigong rock slide–debris avalanche, southeast Tibet.Can Geotech J 2012; 49(5):589-606.

[14]	Shugar DH, Jacquemart M, Shean D, Bhushan S, Upadhyay K, Sattar A, et al.A massive rock and ice avalanche caused the 2021 disaster at Chamoli, Indian Himalaya.Science 2021; 373(6552):300-306.

[15]	Jiang R, Zhang L, Peng D, He X, He J.The landslide hazard chain in the Tapovan of the Himalayas on 7 February 2021.Geophys Res Lett 2021; 48:e2021GL093723.

[16]	Dubey S, Sattar A, Goyal MK, Allen S, Frey H, Haritashya UK, et al.Mass movement hazard and exposure in the Himalaya.Earth's Future 2023; 11(9):e2022EF003253.

[17]	Dal SF Sasso, Sole A, Pascale S, Sdao F, Bateman A Pinzón, Medina V.Assessment methodology for the prediction of landslide dam hazard.Nat Hazards Earth Syst Sci 2014; 14(3):557-567.

[18]	Tacconi C Stefanelli, Segoni S, Casagli N, Catani F.Geomorphic indexing of landslide dams evolution.Eng Geol 2016; 208:1-10.

[19]	Tacconi C Stefanelli, Vilímek V, Emmer A, Catani F.Morphological analysis and features of the landslide dams in the Cordillera Blanca, Peru.Landslides 2018; 15(3):507-521.

[20]	Sansone S, Zugliani D, Rosatti G.A mathematical framework for modelling rock–ice avalanches.J Fluid Mech 2021; 919:A8.

[21]	Gill JC, Malamud BD.Hazard interactions and interaction networks (cascades) within multihazard methodologies.Earth Syst Dyn 2016; 7(3):659-679.

[22]	Pudasaini SP, Krautblatter M.A two-phase mechanical model for rock–ice avalanches.J Geophys Res Earth Surf 2014; 119(10):2272-2290.

[23]	Chen H, Zhang L.EDDA. 1.0: integrated simulation of debris flow erosion, deposition and property changes.Geosci Model Dev 2015; 8(3):829-844.

[24]	Shen P, Zhang L, Wong HF, Peng D, Zhou S, Zhang S, et al.Debris flow enlargement from entrainment: a case study for comparison of three entrainment models.Eng Geol 2020; 270:105581.

[25]	Crosta GB, Imposimato S, Roddeman DG.Numerical modelling of large landslides stability and runout.Nat Hazards Earth Syst Sci 2003; 3(6):523-538.

[26]	Hungr O, McDougall S.Two numerical models for landslide dynamic analysis.Comput Geosci 2009; 35(5):978-992.

[27]	Bartelt P, Buehler Y, Christen M, Deubelbeiss Y, Graf C, McArdell B, et al.A numerical model for debris flow in research and practice: user manual v1.5 (debris flow). WSL Institute for Snow and Avalanche Research SLF, Davos Dorf (2013)

[28]	Mergili M, Fischer J, Krenn J, Pudasaini SP.r.avaflow v1, an advanced open-source computational framework for the propagation and interaction of two-phase mass flows.Geosci Model Dev 2017; 10(2):553-569.

[29]	Peng DL, Zhang LM, Jiang RC, Zhang S, Shen P, Lu WJ, et al.Initiation mechanisms and dynamics of a debris flow originated from debris–ice mixture slope failure in southeast Tibet, China.Eng Geol 2022; 307:106783.

[30]	Iverson RM, Reid ME, Logan M, LaHusen RG, Godt JW, Griswold JP.Positive feedback and momentum growth during debris-flow entrainment of wet bed sediment.Nat Geosci 2011; 4(2):116-121.

[31]	Pudasaini SP.A general two‐phase debris flow model.J Geophys Res Earth Surf 2012; 117(F3):F03010.

[32]	Pudasaini SP, Mergili M.A multi-phase mass flow model.J Geophys Res Earth Surf 2019; 124:2920-2942.

[33]	Pudasaini SP, Krautblatter M.The mechanics of landslide mobility with erosion.Nat Commun 2021; 12(1):6793.

[34]	Pudasaini SP.Unified mechanical erosion model for multiphase mass flows.2022. ar Xiv: 10880.

[35]	Mergili M, Jaboyedoff M, Pullarello J, Pudasaini SP.Back calculation of the 2017 Piz Cengalo–Bondo landslide cascade with r.avaflow: what we can do and what we can learn.Nat Hazards Earth Syst Sci 2020; 20(2):505-520.

[36]	Mergili M, Frank B, Fischer J, Huggel C, Pudasaini SP.Computational experiments on the 1962 and 1970 landslide events at Huascarán (Peru) with r.avaflow: lessons learned for predictive mass flow simulations.Geomorphology 2018; 322:15-28.

[37]	Doman M, Shatoba K, Palmer A.A mega dam on the Great Bend of China [Internet].New York City: ABC News; 2021 May 25 [cited 2025 Apr 8]. Available from: https://www.abc.net.au/news/2021-05-25/chinas-plan-to-build-mega-dam-on-yarlung-tsangpo-brahmaputra/100146344.

[38]	Rao MP, Cook ER, Cook BI, D RD’Arrigo, Palmer JG, Lall U, et al.Seven centuries of reconstructed Brahmaputra River discharge demonstrate underestimated high discharge and flood hazard frequency.Nat Commun 2020; 11(1):6017.

[39]	Rasul G.Water for growth and development in the Ganges, Brahmaputra, and Meghna basins: an economic perspective.Int J River Basin Manag 2015; 13(3):387-400.

[40]	Nuimura T, Sakai A, Taniguchi K, Nagai H, Lamsal D, Tsutaki S, et al.The GAMDAM glacier inventory: a quality-controlled inventory of Asian glaciers.Cryosphere 2015; 9(3):849-864.

[41]	Zhao FY, Long D, Li XD, Huang Q, Han PF.Rapid glacier mass loss in the southeastern Tibetan Plateau since the year 2000 from satellite observations.Remote Sens Environ 2022; 270:112853.

[42]	Zhang W.Identification of glaciers with surge characteristics on the Tibetan Plateau.Ann Glaciol 1992; 16:168-172.

[43]	Wang Z, Hu K, Ma C, Li Y, Liu S.Landscape change in response to multiperiod glacial debris flows in Peilong catchment, southeastern Tibet.J Mt Sci 2021; 18(3):567-582.

[44]	Li Z, Zhou F, Han X, Chen J, Li Y, Zhai S, et al.Numerical simulation and analysis of a geological disaster chain in the Peilong valley, SE Tibetan Plateau.Bull Eng Geol Environ 2021; 80(4):3405-3422.

[45]	Wyk V, de M Vries, Wickert AD.Glacier Image Velocimetry: an open-source toolbox for easy and rapid calculation of high-resolution glacier velocity fields.Cryosphere 2021; 15(4):2115-2132.

[46]	Julien PY, Lan Y.Rheology of hyper concentrations.J Hydraul Res 1991; 117:346-353.

[47]	O JS’Brien, Julien PY.Laboratory analysis of mudflow properties.J Hydraul Res 1988; 114:877-887.

[48]	Takahashi T, Nakagawa H, Harada T, Yamashiki Y.Routing debris flows with particle segregation.J Hydraul Res 1992; 118:1490-1507.

[49]	Hanson GJ, Simon A.Erodibility of cohesive streambeds in the loess area of the midwestern USA.Hydrol Process 2001; 15(1):23-38.

[50]	Graf WH.Hydraulics of sediment transport. Water Resources Publications, Littleton (1984)

[51]	Annandale GW.Scour technology-mechanics and engineering practice. McGraw-Hill, New York City (2006)

[52]	Chang DS, Zhang LM.Simulation of the erosion process of landslide dams due to overtopping considering variations in soil erodibility along depth.Nat Hazards Earth Syst Sci 2010; 10(4):933-946.

[53]	Singh VP, Scarlatos PD.Analysis of gradual earth-dam failure.J Hydraul Eng 1988; 114(1):21-42.

[54]	Jacquemart M, Welty E, Leopold M, Loso M, Lajoie L, Tiampo K.Geomorphic and sedimentary signatures of catastrophic glacier detachments: a first assessment from Flat Creek, Alaska.Geomorphology 2022; 414:108376.

[55]	Jian J, Webster PJ, Hoyos CD.Large‐scale controls on Ganges and Brahmaputra River discharge on intraseasonal and seasonal time‐scales.Q J R Meteorol Soc 2009; 135:353-370.

[56]	Hossain S, Cloke HL, Turner AG, Stephens E.Hydrometeorological drivers of the 2017 flood in the Brahmaputra basin in Bangladesh.Hydrol Earth Syst Sci Discuss. In press.

[57]	Haeberli W, Huggel C, Kääb A, Zgraggen-Oswald S, Polkvoj A, Galushkin I, et al.The Kolka–Karmadon rock/ice slide of 20 September 2002: an extraordinary event of historical dimensions in north Ossetia, Russian Caucasus.J Glaciol 2004; 50:533-546.