RETRACTED: Forecasting the Landslide Blocking River Process and Cascading Dam Breach Flood Propagation by an Integrated Numerical Approach: A Reservoir Area Case Study

Abstract

This paper aims to introduce a numerical technique for forecasting the hazard caused by the disaster chain of landslide blocking river-dam breach floods through an integration of the distinct element method (DEM) and a well-balanced finite volume type shallow water model (SFLOW). A toppling slope in a reservoir area, the southeastern Tibetan Plateau, was chosen for the study. Creep has been observed in the potential instability area, and a possible sliding surface was identified based on the data collected from adits and boreholes. Catastrophic rock avalanches may be triggered after reservoir impoundment, and the associated landslide disaster chain needed to be predicted. First, the landslide blocking river process was modeled by the DEM using the three-dimensional particle flow code (PFC 3D). The landslide duration, runout distance, and kinematic characteristics were obtained. In addition, the landslide dam and barrier lake were constructed. Then, the cascading dam breach flood propagation was simulated using the self-developed SFLOW. The flow velocity, inundation depth, and area were obtained. The hazard maps derived from the combined numerical technique provided a quantitative reference for risk mitigation. The influences of two involved parameters on the final hazard-affected area are discussed herein. It is expected that the presented model will be applied in more prediction cases.

1. Introduction

Many great rivers flow through the southeastern margin of the Tibetan Plateau, and their rapid downward incision caused by sustained bedrock uplift shapes the current deep valley landscape [1,2]. This region is, therefore, the most promising area for hydropower development in China. Meanwhile, numerous catastrophic, high-speed landslides have occurred in this alpine canyon region, such as the Yigong long runout rock slide-debris avalanche, the Tangjiashan short runout landslide, and the Baige rock slide avalanche. They blocked narrow rivers and caused a series of cascading hazards, including landslide-generated impulse waves, backwater flooding, and dam breach flooding [3,4,5]. As a result, forecasting the hazard caused by the disaster chain of a landslide-blocking river-cascading dam breach flood in the reservoir area is required for risk mitigation.

Research methods used in landslide runout analysis mainly include empirical methods, experimental methods, and physically-based numerical methods [6,7]. Historical case-based empirical methods that rely on statistical geometric correlations can be used in runout prediction on a regional scale [8]. However, these methods can hardly be applied to local scale hazard mapping because the comprehensive information regarding the dynamic process of moving mass and on the final deposit morphology of landslide dams cannot be accurately evaluated [9]. Physical modeling experiments use chute, flume, or actual rugged basal to investigate the transport mechanism, dynamic process, and deposition area of landslides [10,11,12]. These phenomenological studies can provide insights on runout prediction by analyzing the relevant physical parameters. However, the size effect is always the main limitation of laboratory-scale tests [7,13,14]. To overcome these limitations, physically-based numerical methods have become the most attractive technique with which to investigate landslide-blocking rivers in the past decade [15,16,17,18]. Of these, two types of approaches were well-established, namely, continuum-based methods and discontinuum-based methods. Previous studies have reported that discontinuum models can more effectively simulate mass fragmentation and large displaced fractures, which are crucial factors in an actual event [13,14,19,20,21,22,23]. For this reason, the distinct element method (DEM), as one of the most representative discontinuum-based methods, is used in this study to forecast the dynamic process of landslide-blocking rivers.

The dam breach flood analysis methods can be grouped into two broad categories: empirical methods and numerical methods [24]. Based on a database of documented dam failure cases, the empirical methods predict breaching parameters using statistically derived regression equations [25,26,27,28,29]. Because the dam breach process and the flood propagation cannot both be obtained from empirical models, the physically based numerical models have been favored by most researchers in the past decade [30]. They can be further classified into simplified and detailed models according to the degree of generalization of the physical process [24]. Of these, the simplified, physically based numerical models can simulate the entire dam breach process and produce the associated flow discharge hydrograph by incorporating breach morphology, sediment transport, and side slope stability [30,31,32,33]. Multidimensional, physically based models can simulate the breach flood characteristics and predict flood propagation along the downstream river valley in more detail [4,34,35]. Among these, two-dimensional (2D) shallow water equations (SWEs) have been developed to describe the flow process resulting from dam breaches [36,37,38]. They were derived by simplifying the Navier–Stokes equations with depth integration based on the boundary conditions of bottom and surface water. To solve them numerically, the finite volume method (FVM) was used, as it has advantages in an irregular computational domain such as easy implementation, mass and momentum conservation, and high efficiency [39,40,41,42]. For this reason, a self-developed, well-balanced finite-volume-type shallow water model (SFLOW) was used in this study to predict the cascading dam breach flood propagation.

Through simulating a disaster chain, the landslide-affected area, barrier lake storage capacity, and breach flood inundation can be quantitatively assessed. These parameters are crucial targets for landslide hazard mapping and risk mitigation. Previous studies have mainly focused on back analysis of past catastrophic events [9], and several works have been related to single hazard predictions [20,43,44]. However, few attempts have been made to forecast a disaster chain by combining different numerical models [5]. The difficulties include the identification of potential sliding surfaces and the construction of landslide dams.

With the aforementioned motivation, this study aims to present a numerical approach for landslide hazard chain prediction by integrating DEM-based landslide runout simulation and SFLOW-based flood propagation simulation. A toppling slope in a reservoir area, southeastern Tibetan Plateau, was chosen as the study area. Its potential area of instability has been identified and may slide along an existing fracture surface after reservoir impoundment. The landslide motion process and deposition morphology were simulated by the DEM using the three-dimensional particle flow code (PFC 3D). Also, the landslide dam topography and associated dammed lake were constructed. Following this, the cascading dam breach flood propagation in the reservoir area was modeled using SFLOW, and flood inundation maps were obtained. Finally, the influences of two parameters from the combined model on the hazard prediction results were investigated.

2. Methodologies

2.1. Distinct Element Method

The original version of the distinct element method (DEM) was first proposed by Cundall (1971) [45] for the analysis of rock mechanic problems. It was later applied to granular material by Cundall and Strack (1979) [46]. The three-dimensional particle flow code (PFC 3D) developed by Itasca Consulting Group, Inc. is a simplified implementation of the DEM. It treats the sliding mass as an assembly of solid elements and describes the motion process from the perspective of particle interactions. Its dynamic behavior is represented numerically using a time-stepping algorithm, and the solution scheme is based on the explicit finite-difference method. The calculation cycles performed in the PFC alternate between the application of Newton’s second law to the particles and a force–displacement law at the contacts (Figure 1). From this, the velocity and displacement of each particle can be obtained in real time at a time step scale. Because PFC 3D performed well when simulating past events [13,20,21,37,47], it was used for forecasting the landslide blocking river process. Following this, the landside runout behavior and deposit morphology could be obtained.

Ball and wall are the two basic contact objects in the PFC. They are defined as non-deformable rigid elements, but allow for interpenetration. Using different elements, both ball–wall and ball–ball models are available for constructing a landslide numerical model [48]. Regarding the ball–wall model (also called the smooth bottom model), the sliding mass and slip surface are modeled as ball assembly and walls, respectively. This model is appropriate for a landslide scenario in which a slip surface is known and mass erosion or entrainment is not considered [13,19,21,47,48,49]. Otherwise, the ball–ball model is preferable, although its computation is expensive.

A well-recognized contact model is the preferred option for forecasting the landslide motion process. Therefore, the parallel bond contact model proposed by Potyondy and Cundall (2004) [50] was adopted, as it has many successful applications in landslide runout simulations [14,20,47]. It is a linear-based bonded particle model and contains three groups. Of these, the linear group can be reactivated when the surface gap is less than or equal to zero. It can simulate the mechanical behavior of recontact among multiple separated blocks. The dashpot group, as a kind of viscous damping, uses a spring–dashpot system and adds normal and shear dashpots at each contact. It can illustrate the energy dissipation caused by particle collisions. The parallel bond group is removed when the bond is broken, and broken parallel bonds cannot be activated again. This group can mimic mass separation and unfragmented blocks. The properties of each group can be described using multiple micro-parameters:

$$\left\{\begin{array}{c}〈k_n,k^{\mathrm{*}},{\mu }_b〉, \mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}\hfill \\ 〈{\beta }_n,{\beta }_s〉, \mathrm{d}\mathrm{a}\mathrm{s}\mathrm{h}\mathrm{p}\mathrm{o}\mathrm{t} \mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}\hfill \\ 〈\overline{\lambda },{\overline{k}}_n,{\overline{k}}^{\mathrm{*}},{\overline{\sigma }}_c,{\overline{\tau }}_c〉, \mathrm{p}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{l} \mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d} \mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}\hfill \end{array}\right.$$

(1)

where $k_n$ and ${\overline{k}}_n$ are the normal stiffness of the ball element and bond, respectively; $k^{\mathrm{*}}$ and ${\overline{k}}^{\mathrm{*}}$ are the normal-to-shear stiffness ratios of the ball element and bond, respectively; ${\mu }_b$ is the ball friction coefficient; ${\beta }_n$ and ${\beta }_s$ are the normal and shear critical damping ratios, respectively; $\overline{\lambda }$ is the radius multiplier used to set the parallel-bond radii; and ${\overline{\sigma }}_c$ and ${\overline{\tau }}_c$ are the normal and shear strengths of the bonds, respectively. Notably, the $k_n$ and ${\overline{k}}_n$ are radius-dependent micro-parameters which can be expressed as:

$$\left\{\begin{array}{c}{k}_n=4R{E}^{*}, \mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}\hfill \\ {\overline{k}}_n={\overline{E}}^{*}/(R_1+R_2), \mathrm{p}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{l}\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}\hfill \end{array}\right.$$

(2)

where $E^{*}$ and ${\overline{E}}^{*}$ are the effective moduli of the ball element and bond, respectively, and $R$ is the ball radius.

PFC derives macro-scale material properties from the interactions of micro-scale properties [20,50]. However, the universal analytical solution between the macro- and micro-parameters has not been established. The most frequently used method to determine the values of micro-parameters is the trial and error method [19,21]. It utilizes a series of numerical tests to calibrate micro-parameters based on the macroscopic responses of samples, which results in low efficiency. Alternatively, some recently proposed parameter calibration methods use machine learning algorithms, such as support vector machines [47] and deep neural networks [51]. However, these approaches require a large amount of sample data from laboratory tests for training sets, which also requires time and effort.

2.2. SFLOW Model

The self-developed SFLOW model, introduced by Han et al. (2017) [41], was used for simulating the cascading dam breach flood propagation. It was constructed on the 2D shallow water equations (SWEs) which are a type of simplified implementation which depth-integrate the 3D Navier–Stokes equations [52]. In the SFLOW, the governing equations for free-surface shallow flow were corrected by setting the basal flow frictional resistance coefficients with the quadratic rheologic friction model (Equation (3)) [36]. Using this friction model, the frictional, viscous, and turbulent effects on the shallow water flow could be incorporated [53].

$$S_f=S_1+S_2+S_3=\frac{\tau }{{\rho }_m}+\frac{k\eta U}{8{\rho }_m{h}_w}+\frac{g{n}_{td}^2{U}^2}{{h_w}^{1/3}}$$

(3)

where $S_f$ is the flow frictional resistance coefficient; $S_1$, $S_2,$ and $S_3$ are the fluid yield factor, fluid viscous factor, and fluid turbulent factor, respectively; $\tau$ is the yield stress; ${\rho }_m$ is the sediment density; $k$ is the resistance coefficient; $\eta$ is the viscosity of the flood; $U$ is the flow velocity; $h_w$ is the water depth; $g$ is the gravitational acceleration; and $n_{td}$ is the equivalent Manning resistance coefficient that incorporates both turbulent boundary friction and internal collisional stresses, which is expressed as [53]:

$$n_{td}=0.0538n\mathrm{e}\mathrm{x}\mathrm{p}\left(6.0896C_v\right)$$

(4)

where $n$ is the Manning resistance coefficient and $C_v$ is the sediment concentration by volume.

In the SFLOW, the governing equations were solved by a well-balanced numerical scheme based on non-staggered rectangular grids. First, a Godunov-type finite volume method (FVM) was applied to update the flow variables to guarantee first-order accuracy in time. Then, a linear slope-limited reconstruction method was used to obtain second-order accuracy in space [54]. Finally, the bed slope and friction terms were approximated by a central differentiation approach and a full implicit scheme, respectively, to ensure the stability of the numerical scheme [55]. Notably, the performance of such a well-balanced numerical scheme has been verified in some existing studies [37,39,40,42].

To quantitatively validate whether the SFLOW model can adapt to complex topography, three types of irregular terrain were used to perform simulation tests of flood propagation. In the first scenario, the flood flowed through the riverbed without obstruction (Figure 2a). In the second scenario, the front flood bypassed the bulge terrain and was divided into two currents (Figure 2b). In the third scenario, the flood flowed into the marsh terrain and filled it (Figure 2c). Consequently, the simulation results indicated that the SFLOW model is capable of modeling dam breach flood propagation on irregular topography. Following this, the flow velocity, inundation depth, and area could be obtained.

The authors developed software with a graphics user interface (GUI) based on the SFLOW model (Figure 3). It includes two menu bars for uploading elevation and parameter files and outputting and visualizing simulation results; a dark blue rectangular window for recording flow depth and velocity along with elapsed time; a chart for displaying critical parameters; and a circular window for updating the calculation progress in real time.

In a simulation, the computational domain is a user-defined complex terrain consisting of discrete elevation points with a fixed spacing (Figure 4a). Benefiting from this, the breach geometry (i.e., breach width and depth) can be constructed by modifying the present terrain. The inflow boundary conditions were set by a flood hydrograph at breach (e.g., a generalized flood hydrograph in Figure 4b), while the outflow and control boundaries were adapted to the actual terrain.

3. Application to a Reservoir Area Case

3.1. Overview of the Mogu Toppling Slope

3.1.1. Geographical Setting and Engineering Background

The hydropower station is geographically located in the confluence of the YL River, XS River, and QD River, southeastern margin of the Tibetan Plateau, China (Figure 5a). The Mogu toppling slope is on the left bank of the XS River, 3 km upstream of the dam site. The elevation of this deformation body varies from 2960 to 3200 m, and the surface area is 1.35 × 10⁶ m² with an average gradient of 32.5°. The slope area susceptible to failure should first be determined, and the disaster chain of the landslide-blocking river-dam breach flood also needs to be predicted for risk mitigation purposes.

3.1.2. Geological Setting and Slope Deformation Mechanism

The reservoir area is tectonically located in the northeast of the Yajiang–Daocheng block (Figure 6a). In this region, large-scale faults have not developed, except for some early Pleistocene minor active faults. The geologic strata exposed at the Mogu slope are mainly composed of Triassic silty slate (Figure 6b). Under such a geological setting, the slope failure mode is mainly associated with the angular relationship between the slope surface and minor structures such as slaty cleavage.

The attitude of the slaty cleavage in the fresh rock is 60°N~70°W/NE∠60°~80°. The strike of the strata is subparallel to the slope face. Thus, the Mogu slope is a bedded anti-dip rock slope, and it satisfies the basic kinematic conditions of toppling failure mode. For this slope structure, the surficial steep-inclined rock slabs are prone to bending towards the deep valley under self-gravity [56,57]. Accordingly, the closer the distance from the slope surface, the stronger the degree of toppling deformation and the smaller the dip angle of slaty cleavage.

3.1.3. Toppling Zonation and Potential Instability Area

Two deep gullies divide the entire slope into three zones along the river flow direction: Ⅰ, Ⅱ and Ⅲ (Figure 5b). Geophysical methods, including adits prospecting and borehole drilling, were conducted to investigate the internal deformation characteristics of each slope zone (Figure 7a).

For slope zones Ⅰ and Ⅲ, explorations demonstrated that the maximum horizontal and vertical toppling depths were 97 and 42 m, respectively. The surficial rock mass was strongly toppled, but not fractured or broken. Failure precursors such as large-scale tension cracks, slip zones, and creep were not observed either on the slope surface or in the internal rock mass. This indicates that the two toppling zones were stable as a whole; thus, they were not investigated further.

For slope zone Ⅱ, significant differences in the deformation characteristics were observed compared with the aforementioned two zones. Two exploration adits revealed that the rock mass in zone Ⅱ could be separated into four sections according to the degree of toppling deformation (Figure 7b,c), namely, a fracture zone, strongly toppled zone, slightly toppled zone, or un-toppled zone. Of these, the average horizontal depth of the first section was approximately 20 m, and the vertical depth increased with the elevation. The rock mass within it was broken with tension cracks, which were filled by debris and angular gravels (Figure 7(b-1)). The structure of the rock mass was catalastic. A slip zone located between the fracture zone and the strongly toppled zone was observed in PD02 (Figure 7(b-2)). It had a thickness varying from 40 to 60 cm and was mainly composed of off-white debris mingled with fine gravel and mud. In addition, a toppling fracture surface at the boundary of the first two sections was formed in PD04 (Figure 7(c-1)). Outer rock creep along the fracture surface was observed in both adits, and rock collapse occurred at the toe of the slope (Figure 5c). For the other sections, no failure precursors were founded (Figure 7(b-3,b-4,c-2,c-3)).

Xu et al. (2022) [58] reported that the toppling fracture surface in slope zone Ⅱ was a possible sliding surface, and the upper fracture zone was the most likely potential instability area (Figure 8). The sliding surface may currently be developing and extending to the entire slope after reservoir impoundment. A retrogressive rock avalanche (hereafter referred to as a Mogu landslide) may be triggered by the decrease in the friction coefficient at the sliding surface [59]. For this reason, only the identified potential instability area was studied.

3.2. Landslide Blocking River Process Modeling

3.2.1. Construction of the Landslide Numerical Model

For the Mogu landslide, the potential instability area was identified; erosion was not considered due to the short runout distance. Therefore, the ball–wall model was adopted to construct a landslide numerical model. The modeling procedures were carried out using Rhino 6 software and can be described as follows. First, the present topography was constructed based on a 12.5 m digital elevation model (Figure 9a). It was divided into three independent parts (left bank, river, and right bank) to facilitate the assignment of different parameters. Second, the landslide boundary was determined by importing the boundary of the identified potential instability area in Section 3.1.3. Third, the depth of the sliding surface was extrapolated based on the records from exploration adits and drilling boreholes; the derived profile AA’ is shown in Figure 8. From this, the basal slip plane was approximated by projecting and interpolating the records along the profile AA’ (Figure 9b). The topography, including the sliding surface, was then represented by 67,265 triangle meshes used for DEM simulation (Figure 9c). Finally, the sliding mass area was generated by converting the enclosed space bounded by the present topography and the sliding surface into a geometric entity. Its volume was 9.1 × 10⁶ m³ (volume is obtained using the Rhino 6 software), and it was filled by 76,143 balls with radii between 2.0 and 3.32 m (Figure 9c).

3.2.2. Macro-Parameters and Micro-Parameters

The macro-parameters of physical samples, including peak shear strength (${\tau }_p$), cohesion ($c$), and internal friction angle ($\phi$), were first obtained from a direct shear test of laboratory rock (Figure 10). Next, a same-scale numerical shear box filled with 8572 balls with radii between 0.002 and 0.00332 m was constructed (Figure 11a). Then, direct numerical shear tests were repeated by adjusting the micro-parameters until the recorded shear stress–shear–displacement curves were close to the laboratory curves (Figure 11b). Specifically, trial and error was adopted due to the limited laboratory data. Figure 11c presents the contact distribution of a numerical sample after 20 mm displacement under constant normal stress of 800 kPa. The figure shows that almost all bonds in the middle of the sample were broken. These broken parallel bonds cannot be reactivated, thus can simulate rock mass fragmentation as described above. Finally, a comparison of the results obtained from numerical and laboratory tests is plotted in Figure 12. As indicated in the figure, the simulation curve and the derived macro-parameters were close to those of laboratory tests. Therefore, the calibration procedure was finished, and the obtained micro-parameters are listed in

Table 1. Micro-parameters of the numerical direct shear test and numerical landslide model.

Micro-Parameter	Symbol	Numerical Direct Shear Test	Numerical Landslide Model
Linear group
Number of ball elements	$n_b$	8572	76,143
Minimum ball radius (m)	$R_{min}$	0.002	2
Ball radius ratio	$R_{max}$/$R_{min}$	1.66	1.66
Ball density (kg/m3)	${\rho }_b$	2400	2400
Ball friction coefficient	${\mu }_b$	0.3	0.3
Normal stiffness (N/m)	$k_n$	8 × 105~1.33 × 106	8 × 108~1.33 × 109
Normal-to-shear stiffness ratio	$k^{*}$	3	3
Parallel bond group
Bond normal stiffness (N/m3)	${\overline{k}}_n$	1.51 × 1010~2.5 × 1010	1.51 × 107~2.5 × 107
Bond normal-to-shear stiffness ratio	${\overline{k}}^{*}$	3	3
Normal strength (Pa)	${\overline{\sigma }}_c$	2.1 × 107	2.1 × 107
Shear strength (Pa)	${\overline{\tau }}_c$	2.2 × 107	2.2 × 107
Dashpot group
Normal critical damping ratio	${\beta }_n$	0.1	0.1
Shear critical damping ratio	${\beta }_s$	0.1	0.1

.

Studies have shown that when the index RES (Equation (5)) is larger than 10, the number and radii of particles have a weak effect on the macro-parameters [61]. In this study, the RES was equal to 18.8, much larger than the threshold value. Therefore, the number and radii of the particles were adapted to the size of the numerical model.

$$\mathrm{R}\mathrm{E}\mathrm{S}=\left(\frac{L_{min}}{R_{min}}\right)\times \left(\frac{1}{1+R_{max}/R_{min}}\right)$$

(5)

where $\mathrm{R}\mathrm{E}\mathrm{S}$ is the number of ball elements in the minimum side length of the numerical model, and $L_{min}$ is the side length of the shear box in this study.

The ball radius range used for landslide simulation was enlarged by 1000 times due to the limitations of computing efficiency and different scales between the laboratory test and the actual event. Due to this, some radius-dependent micro-parameters, such as $k_n$ and ${\overline{k}}_n,$ were readjusted according to Equation (2) [50]. The final micro-parameters of each group used for landslide simulation are listed in Table 1.

For rock avalanches, low residual friction was found between sliding mass and slip surface [14,62,63,64]. Assuming the self-lubrication mechanism, a low friction coefficient (0.05) was assigned for the sliding surface [20,21,22,49]. In contrast, a slightly higher friction coefficient (0.2) was assigned for the riverbed and the topography across the river valley. The increase in the friction coefficient was attributed to the resistance of the river bed and the opposite slope to the landslide debris.

3.2.3. Runout Behavior and Landslide Dam Predictions

The Mogu landslide was initiated by reducing the friction coefficient of the slip surface [21,22] to model the triggering mechanism mentioned in Section 3.1.3. The calculation cycles were terminated when the average velocity of all particles was less than 0.1 m/s. According to the history curves of average velocity and displacement (Figure 13), the whole sliding process lasted 45 s and could be partitioned into three stages, namely, the acceleration stage (0~10 s), the extremely high-speed motion stage (10~20 s), and the accumulation stage (20~45 s).

In the initial stage, the lower particles lost stability first (Figure 14a). Next, the middle and upper particles were initiated in sequence due to the absence of bottom support (Figure 14b). Then, the sliding mass maintained an extremely high speed for 10 s in the second stage (Figure 14c,d). Its average velocity reached its maximum (29.4 m/s) at approximately 15 s (Figure 13a). In this stage, the front particles reached the river valley and their velocity began to decrease due to an increase in the wall friction coefficient and flat terrain. In the third stage, the front particles across the river valley hit the opposite bank. They were deposited at the bottom of the valley, while the rear particles continued to move, but were obstructed by the stopped particles in the front (Figure 14e). Therefore, the average velocity began to gradually decrease during this stage. Nearly all of the particles stopped moving after 45 s (Figure 14f). The elevation of the sliding mass dropped by 213 m, and the released energy was 3 × 10¹³ J. The final horizontal runout distance was only 462.5 m, because the sliding process was severely restrained by the deep and narrow valley. Generally, the Mogu landslide is a high-speed, short-runout rock avalanche, according to the classification scheme updated by Hungr et al. (2014) [65].

To investigate the variation characteristics of the landslide velocity and displacement fields, 6 groups consisting of 22 ball were monitored (Figure 15a). Their runout paths are plotted in Figure 15b. The velocity field of each path showed three stages, which were identical to the aforementioned descriptions. The velocity and displacement profiles of all monitoring groups indicated that the main sliding process lasted 30 s (Figure 15c,d). Regarding velocity profiles, three points can be summarized. First, the monitoring ball velocities fluctuated (the velocities were partially regained). This can probably be attributed to particle collisions on the rugged terrain during sliding process [20]. Second, the peak velocities of the lower groups were smaller than the upper groups (Figure 15c). This is because the lower particles were first obstructed by the opposite slope and, therefore, experienced a shorter acceleration stage. Third, the lower monitoring groups reached their peak velocity earlier than the upper groups did. In other words, the rear particles showed hysteresis during the whole sliding process which was also observed in other rock avalanches [22,47]. Regarding displacement profiles, particle runout distance is positively correlated with elevation (Figure 15d). The maximum and minimum displacements were 691 and 216 m, respectively.

The Mogu landslide blocked the river and formed a dam 1240 m long, 390 m wide, and 56 m high (Figure 14f and Figure 16a). The associated barrier lake was constructed using the ArcGIS software (Figure 16b). Its maximum storage capacity was 1.34 × 10⁷ m³, with a surface area of 0.55 km². The average discharge of the Yalong River was 1860 m³/s, and the net flux of the landslide damming area was assumed to be fifty percent of the average discharge [66]. In this scenario, the water storage of the barrier lake would reach its maximum after four hours. Because most landslide dams are short-lived in the southeastern Tibetan Plateau [67,68], the next section focuses on forecasting the cascading dam breach flood hazard.

3.3. Dam Breach Flood Propagation Modelling

3.3.1. Parameters and Data Preprocessing

Based on the descriptions in Section 2.2, the necessary input data of the SFLOW includes rheological parameters, a raster map and a flood hydrograph.

The rheological parameters used for simulating flood propagation are listed in

Table 2. The parameters used in the simulation of dam breach flood propagation.

Reheological Parameter	Symbal	Value
Sediment density (kg/m3)	${\rho }_m$	1500
Sediment concentration by volume	$C_v$	0.1
Yield stress (Pa)	$\tau$	500
Viscosity of flood ( Pa∙s)	$\eta$	0.0325
Resistance coefficient	$k$	100
Manning resistance coefficient (s/m1/3)	$n$	0.025

. Of these, the parameter $C_v$ is set to 0.1, as suggested by O’Brien (2006) [69], because of the low content of solid materials in the floodwater. The parameter $\tau$ controls the inundation area and $\eta$ affects the flow velocity [41]. The physical and mechanical properties of a dam breach flood should fall in between diluted debris flow and clean water. Therefore, the parameters $\tau$ and $\eta$ of the flood take the middle value of that used in diluted debris flow [70] and clean water. They are defined as 500 Pa and 0.0325 Pa∙s, respectively. The parameter $k$ is related to the surface material, and is defined as 100, as recommended by O’Brien (2006) [69]. The parameter $n$ is dependent on surface roughness, which refers to the value (0.025) in the Tangjiashan and Baige outburst flood cases [33].

A reservoir area (13.5 km²), including the landslide damming area and dam site, is defined as the computational domain (Figure 16b). Then, a 5 m resolution raster map of the simulation domain is converted to a ASCⅡ file, which can be recognized by the developed software. The worst scenario in flood inundation assessment, namely, instantaneous complete dam-break, occurred after the barrier lake reached its storage capacity [38]. The breach depth and top width in this scenario are 56 and 168 m, respectively, based on the cross-section of the landslide dam revealed by the DEM simulation (Figure 16a). For the failure mode of instantaneous complete dam-break, the flood hydrograph at breach is usually generalized using a quartic parabola model (Figure 4b) [71,72]. In this model, the peak discharge occurred at the initial moment, and then the discharge dropped rapidly with the constant decrease in water level in the barrier lake. The area under the flood hydrograph was equal to the volume of water released ($V_w$), namely, the storage capacity of the landslide-dammed lake in this scenario.

The peak discharge ($Q_p$) and duration ($T_{max}$) in the flood hydrograph are two important parameters associated with flood inundation. However, neither of them can be obtained from field measurements for a future event. Alternatively, the parameter $Q_p$ was estimated using five well-recognized empirical models, as shown in

Table 3. Dam breach peak discharge, predicted by different empirical models.

Equation	Investigator	Predicted Value (m3/s)
$Q_p=0.981{\left(h_d{V}_l\right)}^{0.42}$	Costa (1985) [25]	5239
$Q_p=0.72{\left(V_w\right)}^{0.53}$	Evans (1986) [26]	4311
$Q_p=0.0158{\left({\rho }_{w}g{h}_d{V}_l\right)}^{0.41}$	Costa and Schuster (1988) [27]	2979
$Q_p=0.99{\left(d{V}_w\right)}^{0.40}$	Walder and O’Connor (1997) [28]	3513
$Q_p=\sqrt{g{h_d}^5}{{\left(\frac{h_d}{h_r}\right)}^{-1.371}\left(\frac{{V_l}^{1/3}}{h_d}\right)}^{1.536}{e}^{\alpha }$	Peng and Zhang (2012) [29]	1837

Note: $Q_p$ is the peak discharge (m³/s); $h_d$ is the dam height (m); $V_l$ is the lake volume (m³); $V_w$ is the volume of water released (m³); ${\rho }_w$ is the water density (kg/m³); $g$ is the gravitational acceleration (m/s²); $d$ is the drop in water level (m), $h_r$ is referential dam height (1 m); and $\alpha$ equals 1.276, −0.336, and −1.532 for high, medium, and low dam erodibility.

. It was determined by taking the average value of the middle three predicted values, and was equal to 3601 m³/s. Following this, the $T_{max}$ was derived by Equation (6), and it was equal to 4.65 h. Finally, the flood hydrograph used for this scenario was obtained by magnifying the generalized flood hydrograph according to $Q_p$ and $T_{max}$.

$$V_w={\int }_0^{T_{max}}{Q}_{p}Q\left(T\right)\mathrm{d}T$$

(6)

where $Q\left(T\right)$ is the function of the flood hydrograph.

3.3.2. Flood Propagation Predictions

The simulation results of dam breach flood propagation are presented in Figure 17. After the collapse of the landslide dam, the outburst flood first flowed downstream along the river bed of the XS River (to facilitate reading, the locations of the XS River, YL River, and QD river are shown in the top right corner of Figure 17a). At 0.25 h, the front flood reached the confluence of the XS River and the YL River (Figure 17a). Subsequently, the flood diverged at the confluence (Figure 17b). A small portion flowed into the YL River, while the rest of the water continued forward. At 1 h, the front flood hit the dam site (Figure 17c) and began to flow into the QD River. Backwater flooding caused by the obstruction of the dam site can be observed in Figure 17d–f, and therefore, the front water level of the dam site continued to increase in this stage. Finally, the maximum inundation depth and area were 48 m and 0.79 km², respectively. On the one hand, the inundation depth was significantly smaller than the dam height (294 m). On the other hand, the inhabitants in the potential inundation area were evacuated. Therefore, the risk caused by the possible Mogu landslide blocking river event can be controlled. Even so, monitoring the deformation of the Mogu slope is also necessary.

The maximum inundation depth and flow velocity in the flood propagation were monitored, and the results are shown in Figure 18. Regarding the maximum inundation depth, it increased rapidly in the first 0.25 h, and then experienced a short deceleration process as the flood flowed forward. After the front flood hit the dam site (Figure 17c), it gradually increased due to the obstruction of the dam site (Figure 18a). Regarding the maximum flow velocity, it reached its peak (15.6 m/s) within one minute. Then, it decreased rapidly because the energy dissipated due to basal flow friction. Subsequently, the flow velocity was regained, in part, after two strong impacts between the water and the slope/dam site (Figure 18b). Especially for the first impact, the flow direction changed significantly (Figure 18a,b), and the flow velocity returned to 15.4 m/s. Finally, the velocity gradually decreased along with the time elapsed in the stage of backwater flooding.

4. Discussion

4.1. Effect of the Friction Coefficient of the Sliding Surface

The PFC 3D model fails to consider all mechanisms in a rock avalanche. This study focuses on a possible triggering mechanism of the Mogu landslide, namely, the decrease in the friction coefficient at the sliding surface (${\mu }_s$) caused by reservoir impoundment. To investigate the influence of this mechanism on the final deposit shape, a sensitivity analysis of ${\mu }_s,$ varying from 0.05 to 0.4 with an increment of 0.05, was conducted (Figure 19). With increasing ${\mu }_s$ ($\le$0.3), the peak velocity decreased to 14.4 m/s, and more sliding masses were deposited on the left bank slope. As a result, the landslide damming area became smaller. When the ${\mu }_s$ was larger than 0.3, the final landslide dam shape was not sensitive to ${\mu }_s$ (Figure 19b). It was suggested that 0.3 was the threshold of this landslide numerical model. Of course, the universality of the threshold value needs to be further discussed by examining additional cases. The scenario with a ${\mu }_s$ of 0.05 was adopted for predictions mainly due to self-lubrication between the sliding mass and the slip surface during runout [20,21,22,49].

4.2. Effect of the Breach Depth

In Section 3.3, the worst scenario (complete dam-break) was adopted for flood hazard predictions. However, a landslide dam may not be completed destroyed in an actual event, since the dam length along the longitudinal river profile becomes longer as the dam elevation decreases. To investigate the influence of breach depth ($h_b$) on flood hazard mapping, different scenarios were established by modifying it. Assuming that the breach morphology was an inverted isosceles trapezoid (Figure 20a) [24,29,30,31,32], the breaching parameters of each scenario are listed in

Table 4. Breaching parameters of six scenarios used for simulation.

Scenario	Dam Failure Mode	Breach Morphology (Cross-Section)
		$\mathbf{Depth}/{\mathit{h}}_{\mathit{b}}$ (m)	$\mathbf{Top} \mathbf{Width}/{\mathit{W}}_{\mathit{t}}$ (m)	$\mathbf{Bottom} \mathbf{Width}/{\mathit{W}}_{\mathit{b}}$ (m)
No. 1	Complete dam-break	56	186	---
No. 2	Partial dam-break	42	157	73
No. 3		28	159	103
No. 4		21	165	123
No. 5		14	166	138
No. 6		7	165	151

Note: For partial dam-break failure mode, the breach morphology was assumed to be an isosceles inverted trapezoid (the acute angle between leg and base was 45°: $W_t=W_b+2h_b$).

. Then, the flood hydrographs of the inflow points were recalculated according to Table 3 and Equation (6) (Figure 20b), and the breach morphology of each scenario was modified in the digital elevation model.

Figure 20c reveals that the parameter $h_b$ was positively correlated with the maximum inundation depth/area. However, with increasing $h_b$, the curve slope became smaller. This is because that the shape of a landslide-dammed lake is usually a wedge with an inverted trapezoid cross-section in the alpine canyon region. As a result, the increment of water released gradually reduced with the increase in breach depth, which decreased the increment of inundation depth/area.

The simulated results of scenarios No. 3~No. 6 are presented in Figure 21. After the front flood hit the dam site, the front water level of the dam site did not keep rising in scenario No. 5 or No. 6 (Figure 21a,b) due to a limited volume of water released, which was different from the worst scenario (Figure 17). In contrast, the trends of flood propagation in scenarios No. 3 and No. 4 were consistent with the worst scenario, except for the predicted inundation depth and area (Figure 21c,d).

4.3. Limitations of the PFC 3D Model and the SFLOW Model

In this paper, the combinations of PFC 3D and SFLOW models demonstrated enormous potential in forecasting chain hazards. However, there are some limitations of the integrated model. First, the geometry of the potential sliding surface can hardly be accurately constructed. Nevertheless, on the basis of the available data collected from geophysical methods, a potential instability area was determined, and the extrapolated slip surface seems to be a reasonable approximation. Second, the PFC 3D model cannot simulate all failure mechanisms in detail. For example, the loss of partial shear strength along the slip surface may have occurred due to an increase in the porewater pressure after reservoir impoundment. Alternatively, this mechanism was reflected by decreasing the friction coefficient between the sliding mass and the slip surface. Lastly, the SFLOW model has limitations in terms of simulating the evolution of breach morphology. Therefore, an instantaneous dam-break failure mode can be assumed, and different scenarios of breach morphology were established.

5. Conclusions

This study mainly introduces an integrated numerical technique for landslide disaster chain predictions by a reservoir area case in the southeastern Tibetan Plateau, China. Of these, PFC 3D was used for modeling the blocking river process, and the cascading dam breach flood propagation was simulated by the SFLOW. The major findings which were derived are: (1) The toppling fracture surface in zone Ⅱ of the Mogu slope may develop into a sliding surface, and the upper creep zone is the potential sliding mass. A high-speed, short-runout rock avalanche with a volume of 9.1 × 10⁶ m³ may be triggered by the reservoir impoundment. (2) DEM simulation revealed that a landslide with a peak velocity of 29.4 m/s can block the river in 45 s, forming a landslide dam 1240 m long, 390 m wide, and 56 m high. The water storage of the associated barrier lake is 1.34 × 10⁷ m³, which may be filled in four hours. (3) The SFLOW simulation revealed that the flood reached the dam site in one hour, with a flow velocity of 15.4 m/s. The final inundation depth and area were 48 m and 0.79 km², respectively. (4) Sensitivity analyses demonstrated that the friction coefficient at the sliding surface and breach depth are two critical parameters associated with the hazard-affected area. (5) Finally, benefiting from hydraulic engineering, the risks caused by outburst flood can be controlled. It is expected that the method of combined PFC 3D and SFLOW can be improved and applied in more prediction cases to mitigate risks in the chain disaster of landslide blocking river-cascading dam breach floods.