Three-dimensional modeling of the impact behavior of debris flows in areas affected by earthquakes

The 2008 Wenchuan earthquake triggered many landslides and caused high volumes of loose solid material to collect in gully drainage areas; over time, this material eroded and became the source material for multiple debris flows. On August 13, 2010, a heavy rainstorm struck the areas that were impacted by the Wenchuan earthquake. The resulting debris flows caused more than 2000 deaths and destroyed numerous buildings, bridges, and traffic facilities. To mitigate and/or prevent the negative consequences of such disasters in the future, it is important to have reliable and effective methods for predicting the impact behavior of these debris flows. In this study, we investigated the impact behavior of debris flows in post-earthquake areas using a three-dimensional model that is based on the particle flow code PFC3D. We verified our model by simulating a granular flow model test in a laboratory. We quantified the flow process and the impact force of the granular material; by comparing our model data to the test data, we found that the PFC3D model was capable of accurately simulating the dynamic behavior of the granular flows. We then used our model to estimate the propagation processes and the impact forces exerted on check dams by two large debris flows that occurred in the Wenchuan earthquake area. The results agreed with data collected during field surveys, and the PFC3D model can reproduce the dynamic behavior of debris flows in areas that had recently experienced significant seismic activity.


Introduction
Post-earthquake debris flows are a major hazard in mountainous areas affected by mega earthquakes (Lin et al. 2002;Lin et al. 2004;Chen et al. 2009;Li et al. 2016;Dahlquist et al. 2019;Fan et al. 2019). The 2008 Wenchuan earthquake in China triggered approximately 56,000 landslides (Parker et al. 2011;Gorum et al. 2011). Much of the material in these landslides was left behind in river channels or on steep slopes with grades exceeding 1 3 40°. These sediments are quasi-stable during the dry season but can act as source material for debris flows during the rainy season (Zhang et al. 2012a). In the three years following the Wenchuan earthquake, multiple heavy rainfall events resulted in many large-scale debris flows that caused the deaths of ~ 2000 people and impeded restoration and reconstruction efforts in the areas affected by the earthquake. For example, from August 13 to 14, 2010, a heavy rainstorm near Yingxiu Town, the epicenter of the Wenchuan earthquake, triggered multiple hillside and channel debris flows. Among these debris flows was the Hongchun gully debris flow, which blocked and ultimately diverted the Minjiang River so that it flooded the newly rebuilt Yingxiu Town (Tang et al. 2009;Xu et al. 2012). Unsurprisingly, both academicians and government organizations are invested heavily in finding ways to accurately predict the dynamic impacts of debris flows.
At present, there are several effective methods to gain insight into the dynamic behavior of debris flow disasters: field surveys, the interpretation of remote sensing data (Yin et al. 2010;Adam et al. 2015;Ma et al. 2018a,b), theoretical calculations (Ingles et al. 2006;Pastor et al. 2014;Du 2018;Ji et al. 2020;Perna et al. 2022), and laboratory models (Wang et al. 2011;Marco et al. 2016;Zhang et al. 2017).
Collecting data in the field and analyzing remote sensing images provide direct observations of the dynamic characteristics and the impact area of debris flows. For example, Crowley et al. (2003) used remote sensing observations to characterize altered rock areas that represented potential debris flow sources. Zhuang et al. (2010) assessed the probability of the Minjiang River becoming blocked by debris flows in response to the Wenchuan Earthquake. Lu et al. (2012) presented a fast debris flow disasters areas detection method based on remote sensed images. Despite obtaining useful results, this qualitative method requires a lot of manpower, financial resources, and time. As such, this method should be used in conjunction with other debris flow tools to reveal the kinetic characteristics of debris flows.
Laboratory simulations can be used to explore the possible deformation and damage mechanisms of debris flows. For example, Koo et al. (2017) investigated the debris-barrier interaction using flume tests and observed that, without accounting for debris-barrier interactions, the projected maximum momentum of a debris flow during impact would be 30% lower than the testing results. Ng et al. (2017) investigated the effect of the debris flow on rigid barriers using geotechnical centrifuge testing. Tan et al. (2020) introduced a new fast door-opening method implemented in a large-scale physical model built in Hong Kong to study the impacts of debris flows on a flexible barrier. Wendeler et al. (2019) proposed a debris flow load model for flexible barriers in torrent channels based on the results from laboratory tests. Unfortunately, laboratory models are often hindered by space limitations because it is difficult to build a large-scale three-dimensional model that accurately simulates natural processes. In addition, these laboratory tests are time-consuming and labor-intensive.
Numerical methods such as the finite element method (FEM) and the finite difference method (FDM) have been extensively used in the analysis of debris flow behavior as an alternative to or in conjunction with the limit equilibrium method. For example, Zhou et al. (2013a, b) and Zhang et al. (2021) applied the finite difference method (FDM) to numerical simulations of debris flows. Medina et al. (2008) developed a 2D finite volume method (FVM) known as the FLAT model to calculate the dynamic properties of the flows. In addition, Chen and Lee (2000) used the finite element method (FEM) to simulate the sliding process of debris flows in the context of a 3D dynamic model with a Lagrangian framework. It is worth noting that these grid-based methods are mostly used to conduct dynamic analyses of the debris flows. Because these grid-based methods are usually limited to the 1 3 analysis of relatively small displacements due to mesh distortion, these studies mainly focus on analyzing the stability of their models. Recently, some meshfree methods were proposed to simulate the dynamic behavior of debris flows (Cuomo 2020). For example, Dai et al. (2017) proposed a fluid-structure coupling model based on Smoothed Particle Hydrodynamics (SPH) method to estimate the impact force of debris flow on check dams. Li et al. (2018) developed a soft-rigid contact model based on Material Point Method (MPM) for granular flow impact on retaining structures. The discrete element method (DEM) (Cundall 1971) is geared towards analyzing the dynamic characteristics of discontinuous media. The particle flow method, an extension of the discrete element method, has recently been applied to the field of debris flow modeling. Zhang et al. (2012b) used the three-dimensional particle flow discrete element code PFC3D to establish a three-dimensional debris flow model of the Jiweishan Expressway remote debris flows to determine the trajectory, maximum displacement, rock fragmentation, and accumulation morphology of these debris flows. In this study, we use the discrete element method to construct a PFC3D model that will allow us to investigate the movement characteristics and the impact forces of debris flows. Our results indicate that this numerical model accurately predicts the trajectory and impact forces of debris flows.

Simulation method
Discrete element particle flow software is often employed in the fields of geotechnical engineering and mining engineering to analyze the discontinuous deformation and failure caused by slope instability. There are many advantages to using the PFC3D particle flow software to simulate disaster processes such as slope collapse and debris flows. This software package monitors changes in the energy and force chain of the model, which provides valuable insight into the deformation mechanisms of a given landslide. The particle flow discrete element program treats the area in question as a combination of disc particles (ball) and wall boundaries (wall) and then generates a suitable contact model. During the loading process, the particles fall away from the parent body, which allows for more accurate modeling of large deformations such as slip and the separation between rock and soil particles. This software is an iterative solution. During the iteration process, the force-displacement law and Newton's second law are applied alternately. The force-displacement law reflects the changes in the contact force between particles and particles and between particles and walls. Newton's second law reflects particle movements such as translation and rotation. In this study, we assume linear parallel bonding between particles and linear contact between particles and walls.

Parameters calibration and model verification
Before applying the PFC model, it is important to calibrate the micro-parameters used in the simulation. The commonly-used approach is trial-and-error method (Wang and Tian 2018;Huang et al. 2019). In this study, a cylinder specimen (Φ 250 mm × 500 mm) of granular material under axial compression is simulated, as shown in Fig. 1. The granular material is simulated by a series of spherical particles with a radius of 0.006 ~ 0.014 m and a density of 2020 kg/m 3 . The confining pressures exerted on the specimen are 200 kPa, 400 kPa and 600 kPa, respectively. Through the trial-and-error method, a shear strength envelope curve (as shown in Fig. 2) can be obtained by using the micro-parameters listed in Table 1. Then, it is easy to obtain that the cohesion and internal frictional angle of the granular material are 0 kPa and 30° which match the macro-strength of the granular material used in the granular flow model test in Moriguchi et al. (2009). Therefore, the microparameters in Table 1 are suitable to simulate the mechanical behavior of the granular material in the model test.
To verify that our PFC3D model accurately simulates the debris flow and the resulting impact forces, we simulated a granular flow model test described in Moriguchi et al.  about 2.65. The weight of the sand was fixed at 50 kg (Moriguchi et al. 2009). Our threedimensional PFC numerical model uses 11,926 particles with diameter of 0.008 m to simulate the granular material. The granular material contact model assumes linear parallel contact. The wall acts as both boundary and measuring plate. The parameters used in the simulation are calibrated and listed in Table 1. Snapshots of the sliding of the granular material are shown in Fig. 4. At t = 0.2 s, the granular material began to slip with a velocity of ~ 1 m/s. At t = 0.4 s, the front of the granular material is close to the right boundary of the flume and the velocity increases to ~ 2 m/s. At t = 0.8 s, the granular material impacts the measuring plate at a velocity of ~ 3.0 m/s. When t = 1.2 s, most of the granular material has reached the right boundary of the flume. When t = 1.6 s, the velocity of the granular material fell to ~ 1 m/s. By t = 2.0 s, all the granular material had come to rest at the right boundary of the flume.
We simulated changing process of the impact force for different flume inclinations. As shown in Fig. 5, the PFC3D model accurately simulates the evolution of the impact force over time. We simulated five granular flow tests for each slope angle and investigated how the peak impact force varied with the flume inclination ( Fig. 6). In Fig. 6, the purple triangles represent the experimental peak impact forces from Moriguchi et al. (2009) and the red triangles represent the results generated using our PFC model. As expected, the peak impact forces generally increase with the flume inclination. Furthermore, the numerical results agreed well with the experimental data. As we have demonstrated, our PFC3D model accurately simulates the propagation of granular material and provides reasonable estimates of the impact forces that arise during the propagation of the debris flow.
In the PFC simulation, the particle size used in the simulation may affect the numerical results. To investigate the effect of particle size on the impact force, the granular material was discreted into different particle sizes, and then the impact force on the measuring plate was calculated. The results are shown in Table 2. With the increase of the particle size, the peak impact force of the granular material on the plate became a little larger. Therefore, in the PFC simulation, it is better to use the same particle size with the granular material. However, the computational efficiency of the PFC3D model is closely related to the particle numbers used in the simulation. As shown in Table 2, the running time increased sharply with the particle number used in the simulation. Therefore, the determination of the particle size in the PFC simulation should balance the computational accuracy and efficiency.

Study area
The 2008 Wenchuan earthquake triggered thousands of landslides and provided large amounts of loose soil sources in debris flow valleys. After the earthquake, a large number of debris flows triggered by subsequent rainstorms in the earthquake-hit area, such as the "9.24" debris flow event in 2008, the "August 13" debris flow event in 2010, the "July 10" debris flow event in 2013, and the "August 20" debris flow event in 2019, and thus resulted in serious casualties, huge economic loss and long-term impact (Tang et al. 2011a, b;Xu et al. 2012;Yang et al. 2021;Chen et al. 2021;Wang et al. 2022). Yingxiu Town, one of the most severely devastated sites by the Wenchuan earthquake, is prone to rainstorms and experienced significant damage from catastrophic debris flows after the earthquake. On 14 August 2010, in total of 21 debris flows were triggered by the rainstorm around Yingxiu Town (Tang et al. 2011a, b). Among those events, the debris flow of Shaofang gully and Hongchun gully destroyed the National Road G213, Dujiangyan-Yingxiu Expressway, and the S303 provincial road. In addition, the debris flows blocked the Minjing River and  formed barrier dams, and resulted in flooding of the newly reconstructed Yingxiu Town (Tang et al. 2011a, b;2012). According to the in situ tests, the density of the source material in main channel and three upstream branches of the Hongchun gully is 1790-2050 kg/ m 3 , the internal friction angle is about 35°, and the cohesion is 2900 kPa (Ouyang et al. 2015).

Shaofang gully debris flow modeling
Shaofang gully (Fig. 7), which has an area of 0.61 km 2 , is located on the left bank of the Minjiang River near the Longmen Mountain Central Fault (i.e., the Yingxiu-Beichuan Fault). The mouth of the ditch is the intersection of the Dujiangyan-Yingxiu Expressway and Wenchuan National Road G213. This deep-cut structure has eroded through the mid-mountain landform, creating a relative height difference of 1014 m. The gully area is characterized as having a willow leaf geometry. The middle and upper valleys, which have asymmetrical V-shaped profiles, are 100-320 m wide. The mouth of the ditch is narrow with a width at the valley bottom of ~ 35 m. The length of the main ditch is ~ 1.58 km, with an average vertical slope drop of 46.5%. With eight waterfall segments, this ditch follows the "steep-slow-steep" ladder pattern. The mouth of the ditch has an average depth of 10 m and a maximum depth of 21 m. The Shaofang gully watershed mainly consists of Quaternary residual slope deposits, earthquake-triggered landslide deposits, and outcrops of the Proterozoic Chenjiang Jinning Formation. Unsurprisingly, landslide deposits are the primary source of debris flows. After the 2008 Wenchuan earthquake, the volume of loosely deposited material in the Shaofang gully ditch reached ~ 2.52 × 10 6 m 3 . The "8.14" debris flow, which had a debris volume of ~ 2.5 × 10 5 m 3 , only activated 10% of the available source material. Satellite images (with the resolution of 12.5 m) both before and after the debris flow event are used to generate the topography model of the slope and identify the source material area. Using the discrete element method, we simulated the Shaofang gully debris flow and characterized its trajectory and impact force with our PFC3D model. We devised a three-dimensional terrain model with rigid walls and a simulated debris flow sliding body that consists of 7575 spherical particles with a radius of 2.0 m. In this scenario, when particles encounter other particles, we apply the linear and parallel  Table 2. An axial compression test similar to the one in Sect. 3 is simulated and the strength envelope curve is obtained in Fig. 8. The cohesion and the internal frictional angle are 2957 Pa and 35.60°, which are close to the shear strength of the source material in Shaofang gully (Ouyang et al. 2015). Therefore, the parameters in Table 2 are suitable to be used in the PFC model to simulate the dynamic behavior of Shaofang gully.
The simulated trajectories, velocities, and displacements of the Shaofang gully debris flow are shown in Figs. 9, 10, and 11, respectively. The front of the debris flow reached a maximum velocity of 17.5 m/s about 45 s after the initial slope failure while the rear part of the debris flow reached a maximum velocity of 12.5 m/s ~ 58 s after the initial failure. The overall average velocity reached a maximum value of 7.5 m/s after 46 s had elapsed. The total displacements of the front and rear parts of the debris flow were 640 m and 750 m, respectively (average value of 590 m). Eventually, the debris washed down toward the opposite slope and the velocity of the body gradually slowed until it came to rest. To verify our model, we compared the simulated and measured final deposition area of the Shaofang gully debris flow. As shown in Fig. 12, our predicted accumulation is in good agreement with the field survey data. As such, we conclude that our PFC3D model effectively reproduces the three-dimensional motion of a debris flow and that the numerical parameters used in the simulation are reasonable.
To study the impact force of the Shaofang gully debris flow, we added three check dams with a maximum height of 15 m into our model of the debris flow gully (Fig. 13). As shown in Fig. 14

Hongchun gully debris flow modeling
Hongchun gully (Fig. 15) is located in Yingxiu Town, which resides on the left bank of the Minjiang River. The fan-shaped drainage basin of the Hongchun gully has an area of 5.39 km 2 and its altitude varies from 875 to 2140 m. The bedrock is composed of weathered granite, Sinian igneous rock, Carboniferous limestone, and Triassic sandstone. Because the 2008 Wenchuan earthquake occurred on the Yingxiu Beichuan fault, which passes through the Hongchun gully, there is a high volume of loose deposits located on the surrounding slopes (Li et al. 2014;Domènech et al. 2019). At 3:00 am on August 14, 2010, slope failure caused a significant volume of debris (~ 711,000 m 3 ) to wash out at the mouth of the Hongchun gully (Tang et al. 2011a, b). This debris flow blocked the Du-wen highway and the Minjiang River; as a result, a natural debris dam with dimensions of 100 m × 150 m × 10 m (L × W × H) formed in the Minjiang River . The dam caused river diversions and flooding, resulting in 32 deaths and the forced evacuation of more than 8,000 residents. To simulate the motion and impact force changes over time, we used the discrete Fig. 9 Snapshots of simulated propagation of the Shaofang gully debris flow element method to establish a geological model and simulate the motion and impact force of a typical Hongchun debris flow. The three-dimensional terrain was constrained by rigid walls and the debris flow was simulated with 21,545 particles. In this scenario, when particles encounter other particles, we apply the linear and parallel bonding contact model. Because the Hongchun gully is very close to the Shaofang gully and the source materials in these two gullies are similar, we used the same microscopic parameters in both test cases (Table 2). Representative snapshots of the simulated propagation, velocities, and displacements of the Hongchun gully debris flow are shown in Figs. 16, 17, and 18, respectively. The front of the debris flow reached a maximum velocity of 15 m/s about 15 s after the failure, which is consistent with the numerical results based on the SPH model (Cheng et al. 2022). After that time, the flow rushed toward the opposite slope, slowed, and then began to flow in the opposite direction. The rear of the debris flow reached a maximum velocity of 12.5 m/s A comparison of the simulated and measured accumulation area of the Hongchun gully debris flow is shown in Fig. 19. Our results indicate that our PFC3D model results are in good agreement with the site survey data. To examine the impact behavior of the Hongchun gully debris flow, we added three artificial check dams (maximum height of 15 m) into the gully of our model and analyzed the resulting debris flow motion (Fig. 20) and impact forces exerted on the check dams (Fig. 21). The debris flow encountered the check dams ~ 10 s after the initial slope failure; after this time, the impact forces experienced by the three dams increased significantly; the maximum impact forces recorded on dam 1, dam 2, and dam 3 were 200 kPa, 850 kPa, and 370 kPa, respectively. The kinetic energy of the debris flow mainly depends on the mass and velocity. The terrain near the dam 2 is steep which results in a large flowing velocity of the debris flow when reaching the dam 2, as shown in Fig. 20b. Besides, the mass of debris flowing impacting the dam 2 is larger than those impacting to the dam 1 and 3. Therefore, the impact pressure of dam 2 was much larger than that of dam 1 and dam 3.
By applying our PFC3D model to the Shaofang gully and Hongchun gully debris flows, we demonstrated that our model accurately simulates the dynamic behavior of these debris flows such that the data produced in the model is consistent with actual site data collected during field surveys. Furthermore, the impact forces experienced by artificial check dams are reasonable compared to the range of impact force values reported in similar studies.

Discussion
The source material of a debris flow is usually very complicated which usually includes liquid and solid phases. The interaction of different phases may influence the macromechanical property of the geo-material. Therefore, it is better to simulate the mechanical behavior of debris flows using a multi-phase coupled model. However, for the large-scale events such as debris flows, it is difficult to quantitatively consider the microinteraction among different phases when modeling their propagation behavior. Therefore, the debris flow is usually considered as a single-phase viscous fluid to simulate the propagation behavior in the literature. For example, Huang et al. (2015) proposed an SPH model to predict the runout distance of debris flows. The debris flows were simulated as a single-phase Bingham fluid, and the accuracy of the modeling results was verified. Zhu et al. (2021) assumed a debris flow as a viscous fluid to analyze the large deformation behavior. Ray et al. (2022) set up a physical model in the laboratory to calibrate the material contact properties during numerical simulation using PFC software. The efficiency of PFC solver for debris flow modeling was verified. Similar research work can be found in Staron (2008), Tang et al. (2013), Wei et al. (2019), Shen et al. (2020), and Li et al. (2021). Besides, according to the experimental and numerical results presented in the literature Ng et al. 2020;Shi et al. 2022), the impact force of a debris flow mainly depends on the density, velocity, and debris The effect of the pore water pressure on the impact force of debris flows is not significant. Therefore, in the presented study, the debris flow is assumed as a single-phase material for simplified calculation and the material's multi-phase nature is ignored. Hydrostatic approach is widely used to estimate the dynamic impact force (F d ) acting on a check dam due to debris impact (Armanini 1997;Koo et al. 2017). The impact force can be calculated as where ρ is debris density, g is gravitational acceleration, w is the barrier width, and h is debris flow depth. β is a dimensionless empirical factor to account for the dynamic effect of flow impact. In this study, the value of β is taken to be 7 according to Armanini (1997). Based on this empirical equation, the residual pressures on the check dams in the Shaofang gully and Hongchun gully can be estimated. Table 3  where the coefficient k = 1 and α = 2.5 for rigid barriers according to Kwan (2012). Figure 22 shows the comparison of the peak impact forces obtained by Eq. 2 and PFC simulations in this study. It shows that the simulated peak impact forces can match the empirical  (Table 4).

Limitations of the PFC model
As a numerical software based on DEM method, the PFC3D introduced in this study has obvious advantages in modeling extremely large deformation behavior of geo-material. Therefore, it has been widely applied in the simulation of geological disaster propagation, such as flow-like landslide (Feng et al. 2017;Li et al. 2021), rock avalanche (Nichol et al. 2002;Zhang et al. 2019), and debris flow (Zou et al. 2017;Ray et al. 2022) with both dry and saturated source materials. The validity and reliability of the PFC models have been extensively verified. However, the limitations of the PFC model come to light during the application to the debris flow modeling. First, the relationship between micro-parameters of the PFC particles with the macro-mechanical properties of the granular assembly remains unknown. The determination of the micro-parameters in the model still needs theoretical support. Second, some important phenomena involved in debris flow propagation, such as multi-phase interaction (Scaringi et al. 2018), initiation and erosion processes (Zhou et al. 2013a, b), are observed in field investigation, but cannot be considered in the presented PFC model. Therefore, the performance of the model should be improved and further verified. Third, it is difficult to reflect the particle size distribution of natural debris flows in PFC simulation. As a result, the localized impulse load cannot be considered due to the impractical particle size used in the simulation. In addition, the computation efficiency of the PFC3D is quite low when a large number of PFC particles are generated for the three-dimensional modeling of large-scale geological disasters. Therefore, although the simulation results in this study are acceptable, much more attentions and effects should be paid to break through the above limitations and improve the model performance.

Conclusions
Catastrophic debris flows frequently occur in the parts of Wenchuan that were heavily impacted by the 2008 earthquake. In this study, we demonstrated that our PFC3D model could be used to accurately predict the trajectories and impact forces of these debris flows; as such, this model would be a useful tool in seismic hazard analysis and mitigation design. Our conclusions were as follows: A) To verify our PFC3D model, we simulated a granular flow model test and compared the simulated impact force results to the measured data in the model tests. The parameters used in the simulation were determined by trial-and-error method. The general agree- Fig. 19 Comparison of the simulated and measured damage scope for the deposition zone of the Hongchun gully debris flow. ment between the lab data and the model results indicated that the PFC3D model can be used to accurately simulate debris flows. B) We also applied our PFC3D model to the Shaofang gully and Hongchun gully debris flows in post-earthquake areas. Our model accurately reproduced the propagation processes of the debris flows. We also calculated the velocity and displacement time series; with these values, we were able to determine the peak impact force values experienced by hypothetical check dams. These numerical results could provide valuable insight into the engineering and design of future infrastructure that prevents or mitigates the negative effects of debris flows.   Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.