Numerical Investigation on the Kinetic Characteristics of the Yigong Debris Flow in Tibet, China

: To analyze the kinetic characteristics of a debris ﬂow that occurred on 9 April 2000 in Tibet, China, a meshfree numerical method named smoothed particle hydrodynamics (SPH) is introduced, and two-dimensional and three-dimensional models are established in this work. Based on the numerical simulation, the motion process of this debris ﬂow is reproduced, and the kinetic characteristics are analyzed combining with the ﬁeld investigation data. In the kinetic analysis, the ﬂow velocity, runout distance, deposition, and energy features are discussed. Simulation results show that the debris ﬂow mass undergoes an acceleration stage after failure, then the kinetic energy gradually dissipates due to the friction and collision during debris ﬂow propagation. Finally, the debris ﬂow mass blocks the Yigong river and forms a huge dam and an extensive barrier lake. The peak velocity is calculated to be about 100 m/s, and the runout distance is approximately 8000 m. The simulation results basically match the data measured in ﬁeld, thus verifying the good performance of the presented SPH model. This approach can predict hazardous areas and estimate the hazard intensity of catastrophic debris ﬂow.


Introduction
Debris flows are a kind of catastrophic geological hazard which can cause very serious economic and human losses [1]. According to Huang [2], about 80% of large-scale debris flows in China occurred in the Tibetan Plateau, of which over 50% were distributed along the Sichuan-Tibet Highway [3]. Therefore, debris flows pose a serious threat to the human engineering activities in southwest China.
Study of the kinetic characteristics of debris flows can contribute to prediction of the impact area of disasters and has recently attracted the extensive attention of scholars around the world [4]. Field survey combined with remote sensing technology is the most direct approach to obtain the basic dynamic characteristics and the impact area of debris flows. For example, Leonardi and Pirulli [5] installed a filter barrier in an experimental site in the Italian Alps and monitored the load exerted by debris flows and the structural response of the barrier. Liu et al. [6] conducted a field investigation along the Sichuan-Tibet railway and systematically analyzed 55 samples of glacial debris flow deposits to determine their grading and rheological properties. Chang et al. [7] interpreted the distribution characteristics of slag in Hou Gully, Shimian, China using the method of remote sensing and field investigation, and conducted the hazard assessment of the catastrophic mine waste debris flow. Ma et al. [8] carried out field investigations to study the triggering conditions and erosion process of a runoff-triggered debris flow in Miyun County, Beijing, China. Xiong et al. [9] studied the activity characteristics and enlightenment of the debris flow triggered by the rainstorm in Wenchuan County, China. However, field investigation simulation results are compared with the survey data in field, which shows that the SPH model can accurately analyze the kinetic characteristics of catastrophic debris flows.

Geological Setting
On 9 April 2000, a rock avalanche occurred at Zhamunong gully in Bomi County, southeastern Tibet, China. After detaching from the parent rock, it transformed to a highspeed and long-distance debris flow. The geographical coordinates of the debris flow are 30 • 12 11" N, 94 • 58 03" E [36]. Along the banks of Yigong river, the mountains are very high and steep, which are covered with thick snow over 4000 m and with dense vegetation below 3500 m. The valleys in this area are very deep under the erosions of glaciers and rivers.
The rock masses in the study area are mainly granitoid rocks, which have experienced strong weathering and have been partially metamorphosed into slate and granitic gneiss [37,38]. The surface of the slope is composed of quaternary loose colluvial deposits. Thick glaciers and snow covered the slope rock, which could decrease the shear strength of the geomaterial after melting and increase the weight of the sliding mass. Due to the collision between the Eurasian Plate and Indian Plate, active faults are well-developed in the Tibetan Plateau. Jiali Fault and Yigong-Lulang Fault, two of the major active strike-slip faults [39], meet at the mouth of Zhamunong gully, as shown in Figure 1. Earthquakes frequently occurred in the study area. For example, 14 moderate earthquakes (Ms = 4.0-5.9) were recorded around the Yigong lake from 1980 to 1996. Therefore, the tectonic activities in this area caused the rock structure fractured, loosened, and weakened, which provided favorable conditions for the Yigong debris flow occurrence. ergy features. The simulation results are compared with the survey data in field shows that the SPH model can accurately analyze the kinetic characteristics strophic debris flows.

Geological Setting
On 9 April 2000, a rock avalanche occurred at Zhamunong gully in Bomi C southeastern Tibet, China. After detaching from the parent rock, it transformed to speed and long-distance debris flow. The geographical coordinates of the debris f 30°12′11″ N, 94°58′03″ E [36]. Along the banks of Yigong river, the mountains a high and steep, which are covered with thick snow over 4000 m and with dense veg below 3500 m. The valleys in this area are very deep under the erosions of glaci rivers.
The rock masses in the study area are mainly granitoid rocks, which have enced strong weathering and have been partially metamorphosed into slate and gneiss [37,38]. The surface of the slope is composed of quaternary loose colluvial d Thick glaciers and snow covered the slope rock, which could decrease the shear s of the geomaterial after melting and increase the weight of the sliding mass. Du collision between the Eurasian Plate and Indian Plate, active faults are well-devel the Tibetan Plateau. Jiali Fault and Yigong-Lulang Fault, two of the major active slip faults [39], meet at the mouth of Zhamunong gully, as shown in Figure 1. Earth frequently occurred in the study area. For example, 14 moderate earthquakes (M 5.9) were recorded around the Yigong lake from 1980 to 1996. Therefore, the tect tivities in this area caused the rock structure fractured, loosened, and weakened provided favorable conditions for the Yigong debris flow occurrence. The study area belongs to the temperate subhumid plateau monsoon clima Influenced by the warm-wet air currents from the Indian Ocean, the weather is with clear four seasons and abundant in rainfall and sunshine. According to the lo teorological station, the annual rainfall averages 876.9 mm, and the cumulative s hours is 1544 h. It was reported that the antecedent precipitation from 1 to 9 Ap was 42.9 mm, which was a main trigger of the debris flow. According to the three w stations around Yigong area, the mean ground surface temperature in this area gr increased before the occurrence of the Yigong debris flow [40]. It might result in the melting in the source area and, thus, may have increased the pore water pressur geomaterials and decreased the shear strength. The study area belongs to the temperate subhumid plateau monsoon climate area. Influenced by the warm-wet air currents from the Indian Ocean, the weather is humid, with clear four seasons and abundant in rainfall and sunshine. According to the local meteorological station, the annual rainfall averages 876.9 mm, and the cumulative sunshine hours is 1544 h. It was reported that the antecedent precipitation from 1 to 9 April 2000 was 42.9 mm, which was a main trigger of the debris flow. According to the three weather stations around Yigong area, the mean ground surface temperature in this area gradually increased before the occurrence of the Yigong debris flow [40]. It might result in the glacier melting in the source area and, thus, may have increased the pore water pressure in the geomaterials and decreased the shear strength. Figure 2 shows an aerial view of the Yigong debris flow occurred in Zhamunong gully (the base map is taken from Google Earth). About 3.0 × 10 8 m 3 geomaterials slid down along the gully for about 3 min [36,37], and the sliding direction is around 225 • . The horizontal runout distance is about 8000 m, and the vertical dropdown is about 3330 m from its source area at 5520 m to its sediment fan at 2190 m. Deduced from seismic surveillance data, the maximum velocity of the debris flow is higher than 100 m/s, and the average velocity is about 40 m/s [40,41].

Debris Flow Features
Water 2021, 13, x 4 of 21 Figure 2 shows an aerial view of the Yigong debris flow occurred in Zhamunong gully (the base map is taken from Google Earth). About 3.0 × 10 8 m 3 geomaterials slid down along the gully for about 3 min [36,37], and the sliding direction is around 225°. The horizontal runout distance is about 8000 m, and the vertical dropdown is about 3330 m from its source area at 5520 m to its sediment fan at 2190 m. Deduced from seismic surveillance data, the maximum velocity of the debris flow is higher than 100 m/s, and the average velocity is about 40 m/s [40,41].   Figure 4 shows the path profile of this debris flow. In this figure, the original slope surface (blue dashed line) and the present slope surface (green solid line) are from [42]. As shown in Figure 4, the debris flow could be identified by three major zones: source zone, propagation zone, and deposit zone. The characteristics of the three zones are described below.   Figure 4 shows the path profile of this debris flow. In this figure, the original slope surface (blue dashed line) and the present slope surface (green solid line) are from [42]. As shown in Figure 4, the debris flow could be identified by three major zones: source zone, propagation zone, and deposit zone. The characteristics of the three zones are described below.

Propagation Zone
The propagation zone of this debris flow covered an area of about 3.46 km 2 . The axial length of this zone is about 3200 m, and the width ranges from 780 to 1500 m. The elevation of this zone ranges from 3790 to 2840 m, with a height difference of 950 m. The average slope of this zone is about 16.0 • , which was much gentler than that of the source zone. Many boulders are distributed in the gully. Most of these are angular with a diameter of over 0.5 m.

Propagation Zone
The propagation zone of this debris flow covered an area of about 3.46 km 2 . The axial length of this zone is about 3200 m, and the width ranges from 780 to 1500 m. The elevation of this zone ranges from 3790 to 2840 m, with a height difference of 950 m. The average slope of this zone is about 16.0°, which was much gentler than that of the source zone. Many boulders are distributed in the gully. Most of these are angular with a diameter of over 0.5 m. Figure 6 is the front view of the deposit zone of the Yigong debris flow. The elevation of this zone ranges from 2200 to 2800 m, and the average slope is about 6.0°. The area of the debris flow deposition is about 5.0 × 10 6 m 2 , and the average depth of sediment is about 50 m [3]. Due to the high motion velocity, the debris flow flushed into the Yigong river and formed a huge dam and an extensive dammed lake. The location of the lake is shown in Figure 2. The length of the trumpet-shaped dam is about 4.6 km, the maximum width is 3.0 km, and the dam height is 60-120 m. The dam sloped at 5° at the upstream side and 8° at the downstream side [35]. After the dam formation, water level of the Yigong lake continuously rose at a rate of about 1 m/day, which flooded the Yigong tea farm, schools, and villages surrounding the barrier lake. On 10 June 2000, the dam failed and resulted in devastating flooding, which destroyed farms, villages, bridges, and highways along its route. In recent years, the loose sediment was eroded by water from the Zhamunong gully and formed a debris fan in the Yigong river channel.  Figure 6 is the front view of the deposit zone of the Yigong debris flow. The elevation of this zone ranges from 2200 to 2800 m, and the average slope is about 6.0 • . The area of the debris flow deposition is about 5.0 × 10 6 m 2 , and the average depth of sediment is about 50 m [3]. Due to the high motion velocity, the debris flow flushed into the Yigong river and formed a huge dam and an extensive dammed lake. The location of the lake is shown in Figure 2. The length of the trumpet-shaped dam is about 4.6 km, the maximum width is 3.0 km, and the dam height is 60-120 m. The dam sloped at 5 • at the upstream side and 8 • at the downstream side [35]. After the dam formation, water level of the Yigong lake continuously rose at a rate of about 1 m/day, which flooded the Yigong tea farm, schools, and villages surrounding the barrier lake. On 10 June 2000, the dam failed and resulted in devastating flooding, which destroyed farms, villages, bridges, and highways along its route. In recent years, the loose sediment was eroded by water from the Zhamunong gully and formed a debris fan in the Yigong river channel.

Numerical Model
To investigate the kinetic characteristics of the Yigong debris flow, a meshfree numerical method named smoothed particle hydrodynamics (SPH) is applied, and two-dimensional and three-dimensional models were established for the rapid debris flow propagation simulation.

Numerical Model
To investigate the kinetic characteristics of the Yigong debris flow, a meshfree numerical method named smoothed particle hydrodynamics (SPH) is applied, and twodimensional and three-dimensional models were established for the rapid debris flow propagation simulation.

SPH Algorithm
The SPH method was proposed in 1997 for astrophysical applications [43]. Recently, this method has been widely applied to a large variety of engineering fields [44][45][46][47]. Compared to the mesh-based method, the major advantage of the SPH method is to bypass the need for numerical meshes and avoid the mesh distortion issue and a great deal of computational work to renew the mesh [48].
In the SPH method, the subject is represented by a set of particles to which the material properties such as velocity, density, and pressure are associated. The properties are updated for each time step of the simulation following the conservation laws of mass and momentum [49].
In this study, the debris flow is assumed as a kind of weakly compressible viscous fluid. Therefore, the continuity and momentum equations are expressed by: where ρ is the particle density, t is the time, and m is the particle mass, v is the velocity, p is the pressure, and F is the body force. The subscripts i and j are the concerning particle and its neighboring particle. δ is the Kronecker delta and Π is an artificial viscosity, which is used to improve the stability of the numerical results [50]. W is a smooth function. In this model, the cubic B-spline function, originally used by Monaghan and Lattanzio [50], is selected as the smooth function. The function formula is where α d is a normalization factor in two-and three-dimensional space, α d = 15/7 πh 2 and 3/2 πh 3 , respectively. R is the normalized distance between particles i and j, defined as R = r/h. Here, r is the distance between particles i and j. The pressure p can be calculated by an equation of state in this study: where ρ is the density calculated by the continuity equation, see Equation (1). ρ 0 is the reference density which can be measured through laboratory tests. c s is the sound speed at the reference density, which can be set equal to ten times the maximum velocity [51]. γ is the exponent of the equation of state and is usually set to 7.0 for a good simulation of geomaterial flow behavior [52].

Material Model
The debris flow mass is a mixture of water, soil, and rock, which is complicated to describe. Hungr [53] proposed a concept of "equivalent fluid", which was intended to simulate the bulk behavior of the debris flow mass. Rickenmann et al. [54] used three different fluid models based on Voellmy fluid rheology, quadratic rheology, and Herschel-Bulkley rheology to simulate the propagation of debris flows across two-dimensional terrain. Recently, some viscous fluid models have been widely used in the numerical modeling of debris flows [29,55,56]. In the presented SPH model, the debris flow mass is assumed as a Bingham fluid, which is widely used to describe the motion behavior of debris flows because of its simple structure [57,58]. The relationship between the shear stress and the shear strain rate in the Bingham fluid model is defined as where τ is the shear stress of the fluid, η is the yield viscosity coefficient in fluid dynamics, and τ y is the yield shear stress, which is commonly defined as the Mohr-Coulomb yield criterion with the cohesion c and frictional angle ϕ [29,59]. p is the pressure which can be obtained by Equation (3). D and D Π are the strain rate and its second invariant.

Boundary Treatment
SPH method is ideal to deal with the free surface boundary. In the presented model, the free surface is identified through the criterion below: where ρ* is present value of the particle density, which equals the initial density ρ 0 plus the density increment dρ. The density increment dρ can be obtained according to the mass conservation equation, as shown in Equation (1). k is the free surface parameter. When the particle is identified as a free boundary particle, then zero pressure is applied. For the solid wall boundary, ghost particles are placed on the boundary lines to exert repulsive forces and avoid the particles crossing the boundary. The velocities of the ghost particles are set to be zero to satisfy the non-slip boundary condition. For a detailed description of the non-slip boundary condition, please refer to [44].

OpenMP Parallelism
To simulate the propagation of a debris flow across complex terrain, it is necessary to develop a three-dimensional numerical model. In the 3D SPH model, however, the computational efficiency is sharply reduced as the particle number increases. To improve the efficiency, it is necessary to parallelize the numerical code without suffering from a loss of precision.
The open multiprocessing (OpenMP) API for shared-memory programming enables loop-level parallelism by the insertion of pragmas within the source code. By adding special directives at the beginning and end of the loop, OpenMP parallel implementation can be easily conducted. The cycles of the loop are then randomly assigned to the available threads. In the present work, the paralleled numerical code was written in FORTRAN 95, and the program was compiled using Microsoft Visual studio 2015 in a PC with the quad-core 8-thread CPU, Intel Core i7-7820HQ, and run at 2.90 GHz clock with 32 GB main memory under the Windows 10 Professional 64-bit operating system.

Time Integration
In a Lagrangian framework, the coordinates of each particles are updated at each time steps. A velocity Verlet scheme is introduced in this SPH model to perform time integration. X n+1 = X n + V n ∆t + 1 2 a n ∆t 2 V n+1/2 = V n + 1 2 a n ∆t where X, V, and a are the displacement, velocity, and acceleration field, respectively.

Two-Dimensional Modeling
According to the debris flow profile in Figure 4, a two-dimensional SPH model was established in this study to simulate the flow behavior of the Yigong debris flow, as shown in Figure 7a. The number of SPH particles used in the numerical model can simultaneously influence the computational efficiency and accuracy [49]. Therefore, to reach an appropriate balance between the computational efficiency and accuracy, 7662 blue particles were used to represent the debris flow mass, and 5906 grey image particles were used to simulate the bed surface. The diameters of those particles are 8 m. The initial velocities of the particles were set to zero. After slope failure, the debris flow mass particles slide down the slope under the action of gravity, while the boundary particles remain stationary throughout the simulation. According to Li et al. [60], the average density of the Yigong debris flow mass was about 2000 kg/m 3 . The strength characteristics of the debris flow mass were studied through a series of high-speed ring shear tests and rotary shear tests in the previous studies [61,62]. According to the test results, the values of the c and ϕ of the geomaterial can be approximately set to be 10 kPa and 20 • , respectively. The selection of dynamic viscosity η is often challenging. In the previous simulation, Bingham model was widely used to simulate debris flows considering a range of dynamic viscosities from 20 to 500 Pa·s [29,63,64]. The sound speed c s is set to be 10 v max (v max is the maximum velocity). The parameter γ in the equation of state is set to be 7.0 for a good simulation of geomaterial flow behavior.
According to the simulation results, the motion process of the debris flow mass takes about 200 s, which is basically consistent with the witnesses' description [36]. Figure 7 presents the slope configurations at different points in time, which indicates the motion process of the debris flow mass after slope failure. The particles slide down from the top of the Zhamunong gully, and then move along the steep slope by gravity. Finally, these particles run to an equilibrium state and accumulate in the Yigong river channel. The color represents the velocity of the particles, which shows that the maximum flow velocity of the debris flow mass is about 100 m/s. To investigate the kinetic characteristics of the debris flow, Figure 8 shows the time history curves of the flow velocity. The blue, red, and green curves represent the velocity evolutions of the rear edge, front edge, and the average value of the debris flow mass. The peak velocities are 102.6 m/s at the front edge and 72.4 m/s at the rear edge. The average velocity of the debris flow mass during the propagation is approximately 39.8 m/s. After slope failure, debris flow mass rapidly slides down and accelerates due to gravity. In this stage, most of the potential energy of the debris flow mass is converted into kinetic energy. After the peak velocity, the kinetic energy is consumed due to the friction, collision, and the breakage of the sliding mass, and the velocity gradually decreases. The overall performance of the sliding mass is accelerated motion in the period 0-50 s and decelerated motion after 50 s.
In Figure 8, the orange and purple dashed lines are the time history curves of the front velocity obtained by the energy models with and without bed entrainment proposed by Kang et al. (2017). The comparison results show that the variation tendencies of the debris flow velocity obtained by the SPH model and the energy models are similar, and the maximum front velocities are also very close. However, the motion time of the debris flow simulated by SPH model is a bit longer than that simulated by the energy models.

Three-Dimensional Modeling
To simulate the debris flow propagation across 3D complex terrain, the 2D SPH model is developed into a 3D version. In the 3D model, the diameter of the SPH particles is set to be 20 m since the resolution of the digital elevation model (DEM) is 20 m × 20 m. The debris flow mass in the source area is discretized into about 11,000 particles so that the total volume of the debris flow is about 3.0 × 10 8 m 3 . The number of particles along the vertical direction varies in different positions according to the depth of the sliding surface at that position. The strength parameters used in 3D simulation are the same as those used in 2D model. Based on this model, the numerical modeling of the Yigong debris flow motion across 3D terrain is conducted, and the results are shown in Figure 10. The color of the particles in the figures represents the sliding velocity. After slope failure, the debris flow mass goes through an acceleration process since the slope is quite steep in the source area. The maximum sliding velocity is about 98.4 m/s, which appears at 47.5 s after the slope failure. Afterwards, the debris flow mass slows down gradually due to the friction and the collision during the propagation. Finally, the debris flow mass crashes into a mountain on the opposite bank of Yigong river and then blocks the river channel. The whole motion process takes about 200 s, and the final depositions of the debris flow mass on the runout path are shown in Figure 10g. Figure 11 shows the Yigong debris flow deposition. The red dashed line is the simulated debris flow deposition, with the area of 4.76 km 2 , which is close to the measured data 5.0 km 2 [37]. The maximum length and width of the deposition belt are 4.62 and 2.84 km, respectively, which are close to the observed  Figure 9 compares the simulated debris flow deposition with the measured data recorded in [42]. It is obvious that the predicted debris flow deposition area is consistent with the measured data. In order to identity the influence degree of each rheological parameter on the numerical results, sensitivity analysis was conducted by varying the rheological parameters. Table 1 lists the seven calculation cases with different rheological parameters. To explicitly quantify the differences between measured and calculated deposits, the L 2 relative error norm in the deposition depth, ε L2 , was evaluated using the following equation: where Y is the measured deposition depth, ∆Y is the deviation of numerical and measured depth, and N is the total number of points at which the depths are compared. As shown in Figure 9, a total of 11 points with a space of 500 m were selected to calculate the error norm, and the results of all the seven cases are listed in Table 1, which shows that the coefficient of viscosity has more influence on the computing accuracy in comparison with the shearing strength parameters.

Three-Dimensional Modeling
To simulate the debris flow propagation across 3D complex terrain, the 2D SPH model is developed into a 3D version. In the 3D model, the diameter of the SPH particles is set to be 20 m since the resolution of the digital elevation model (DEM) is 20 m × 20 m The debris flow mass in the source area is discretized into about 11,000 particles so tha the total volume of the debris flow is about 3.0 × 10 8 m 3 . The number of particles along the

Three-Dimensional Modeling
To simulate the debris flow propagation across 3D complex terrain, the 2D SPH model is developed into a 3D version. In the 3D model, the diameter of the SPH particles is set to be 20 m since the resolution of the digital elevation model (DEM) is 20 m × 20 m. The debris flow mass in the source area is discretized into about 11,000 particles so that the total volume of the debris flow is about 3.0 × 10 8 m 3 . The number of particles along the vertical direction varies in different positions according to the depth of the sliding surface at that position. The strength parameters used in 3D simulation are the same as those used in 2D model. Based on this model, the numerical modeling of the Yigong debris flow motion across 3D terrain is conducted, and the results are shown in Figure 10. The color of the particles in the figures represents the sliding velocity. After slope failure, the debris flow mass goes through an acceleration process since the slope is quite steep in the source area. The maximum sliding velocity is about 98.4 m/s, which appears at 47.5 s after the slope failure. Afterwards, the debris flow mass slows down gradually due to the friction and the collision during the propagation. Finally, the debris flow mass crashes into a mountain on the opposite bank of Yigong river and then blocks the river channel. The whole motion process takes about 200 s, and the final depositions of the debris flow mass on the runout path are shown in Figure 10g. Figure 11 shows the Yigong debris flow deposition. The red dashed line is the simulated debris flow deposition, with the area of 4.76 km 2 , which is close to the measured data 5.0 km 2 [37]. The maximum length and width of the deposition belt are 4.62 and 2.84 km, respectively, which are close to the observed values of 4.60 and 3.0 km, and its shape is basically in agreement with the observed shape (blue solid line in Figure 11).
To verify the performance of parallel computation, the 3D SPH modeling was carried out using different thread numbers (1, 2, 4, 6, and 8). Figure 12 shows the relationship between the average program running time and the thread number. It is obvious that the computation efficiency of the presented SPH model increases with the thread number.

Analysis of Simulation Results
Runout distance, deposition area, and deposition depth are very important for debris flow disaster prediction. About two months after the Yigong debris flow occurrence, the dam broke down, and most of the debris flow deposit was washed away by the flood. Therefore, it is difficult to measure the post-event topography in the field. In this study, the SPH modeling results showed that the total runout distance of the Yigong debris flow was about 8000 m, which are consistent with the measured results [40]. The final deposition area was obtained by 3D SPH modeling, which is in correspondence with the satellite image provided in [65]. The 2D simulated deposition depths along the topographic profile can match the observed results provided in [42] well, while in the 3D modeling, the comparative analysis of deposition depths was not carried out due to lack of measured data.  To verify the performance of parallel computation, the 3D SPH modeling was carried out using different thread numbers (1, 2, 4, 6, and 8). Figure 12 shows the relationship between the average program running time and the thread number. It is obvious that the computation efficiency of the presented SPH model increases with the thread number.   To verify the performance of parallel computation, the 3D SPH modeling was carried out using different thread numbers (1, 2, 4, 6, and 8). Figure 12 shows the relationship between the average program running time and the thread number. It is obvious that the computation efficiency of the presented SPH model increases with the thread number.

Analysis of Simulation Results
Runout distance, deposition area, and deposition depth are very important for debris flow disaster prediction. About two months after the Yigong debris flow occurrence, the dam broke down, and most of the debris flow deposit was washed away by the flood. Therefore, it is difficult to measure the post-event topography in the field. In this study, the SPH modeling results showed that the total runout distance of the Yigong debris flow was about 8000 m, which are consistent with the measured results [40]. The final deposition area was obtained by 3D SPH modeling, which is in correspondence with the satellite image provided in [65]. The 2D simulated deposition depths along the topographic profile can match the observed results provided in [42] well, while in the 3D modeling, the comparative analysis of deposition depths was not carried out due to lack of measured data.
Velocity is one of the key kinetic characteristics during debris flow propagation, which is difficult to measure in field. According to eyewitness accounts, the total sliding time of the Yigong debris flow was about 3 min. The runout distance was about 8000 m. Therefore, the average flow velocity of the debris flow was estimated to be about 40 m/s. According to the dynamic analysis results [40,60], the maximum velocity during the debris flow propagation was more than 100 m/s. Therefore, the velocity-time history predicted by the SPH model in this work fits the literature data well and is reasonable and reliable.
During debris flow propagation, bed entrainment can increase the source volume and then affect the kinematic characteristics. Recently, several numerical models regarding bed entrainment have been proposed with some satisfactory results [36,40]. However, due to many influencing factors, such as particle shape and gradation, pore water pressure, thickness, and compactness of the bed material, it is quite difficult to effectively incorporate the entrainment algorithm into the presented SPH model. In spite of this imperfection, the SPH model presented in this study can simulate the kinetic characteristics of the Yigong debris flow and reach a certain accuracy.

Conclusions
All over the world, debris flows always lead to property loss and human death. This work investigates the kinetic features of the catastrophic Yigong debris flow in the Tibetan Plateau, China. On the basis of the SPH method, 2D and 3D numerical simulations were conducted to reproduce the motion process of the Yigong debris flow. Based on the numerical results and combined with field investigation data and remote-sensing images, Velocity is one of the key kinetic characteristics during debris flow propagation, which is difficult to measure in field. According to eyewitness accounts, the total sliding time of the Yigong debris flow was about 3 min. The runout distance was about 8000 m. Therefore, the average flow velocity of the debris flow was estimated to be about 40 m/s. According to the dynamic analysis results [40,60], the maximum velocity during the debris flow propagation was more than 100 m/s. Therefore, the velocity-time history predicted by the SPH model in this work fits the literature data well and is reasonable and reliable.
During debris flow propagation, bed entrainment can increase the source volume and then affect the kinematic characteristics. Recently, several numerical models regarding bed entrainment have been proposed with some satisfactory results [36,40]. However, due to many influencing factors, such as particle shape and gradation, pore water pressure, thickness, and compactness of the bed material, it is quite difficult to effectively incorporate the entrainment algorithm into the presented SPH model. In spite of this imperfection, the SPH model presented in this study can simulate the kinetic characteristics of the Yigong debris flow and reach a certain accuracy.

Conclusions
All over the world, debris flows always lead to property loss and human death. This work investigates the kinetic features of the catastrophic Yigong debris flow in the Tibetan Plateau, China. On the basis of the SPH method, 2D and 3D numerical simulations were conducted to reproduce the motion process of the Yigong debris flow. Based on the numerical results and combined with field investigation data and remote-sensing images, the kinetic characteristics of the debris flow were analyzed. Some conclusions can be obtained as below.
In the early stage, debris flow mass slides down along the deep slope in the source area. It experiences an acceleration and reaches its peak velocity (about 100 m/s). In this stage, most of the potential energy of the debris flow mass is converted into kinetic energy. During debris flow propagation in the Zhamunong gully, the kinetic energy continuously dissipates due to the friction and the collision. The velocity gradually slows down in this stage. After rushing out of the Zhamunong gully, the debris flow mass crashes into a mountain on the opposite bank of Yigong river and then blocks the river channel. The velocity evolution of the debris flow is obtained based on numerical results, and the final debris flow deposition is predicted, which basically fits the measured data in field.
Although the SPH model presented in this work can reproduce the motion process and analyze the kinetic characteristics of the Yigong debris flow, there are still some problems need to be solved. For example, the bed material entrainment during debris flow propaga-tion has some effect on the debris flow volume and its kinetic characteristics, but was not considered in this work. On the other hand, the disintegration and fragmentation of the rock blocks was not considered in the presented SPH model, which may lead to some error during the simulation of debris flow propagation. Moreover, high-performance parallel computing technology is necessary to improve the calculation efficiency in 3D modeling.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.