Numerical Investigation of Multiple-Impact Behavior of Granular Flow on a Rigid Barrier

: The debris–barrier interaction issue has gained considerable attention among the engineering community, but most researches have only focused on the single-surge impact condition, with the multiple-surge impact mechanism still lacking clarity. However, multiple-surge impact is more typical in the ﬁeld. Thus, we conduct some numerical simulations based on the discrete element method (DEM) and present a series of results that provide preliminary insights into the multiple-surge impact mechanism. The DEM model is ﬁrstly calibrated using physical experimental results and then used to investigate the ﬂow kinematics, impact dynamics and energy evolution of the successive impact process. The results indicate that compared with single-surge conditions, the barrier is safer under multiple-surge impact as the deposition spreading distance is extended by 6–20% and the impact force is reduced by 6–30%. The dead zone formed by the previous surge behaves as a cushioning layer and a medium for momentum transfer. Three mechanisms of energy dissipation during surge–dead-zone interactions were identiﬁed: friction and penetration at the interaction face between the surge and dead zone, inelastic deformation of the dead zone, and inter-particle interaction within the surge. Each component was analyzed, which shows that inter-particle collision friction accounts for over 60% of the total energy loss during surge–dead-zone interaction. In addition, the performance of granular jump theory in predicting the multiple-surge impact force is assessed, and some possible modiﬁcations are proposed. Finally, some engineering implications from the presented numerical results are discussed.


Introduction
Strengthening the impact resistance of structures to withstand rapid granular flows is of concern to design engineers aiming to prevent the harmful effects of landslide disasters [1,2]. Various types of mitigation structures have been designed to reduce the destructive power of such disasters, for example, check dams, flexible barriers, baffle arrays and slit dams [3][4][5][6][7][8][9][10]. The design of impact-resisting structures relies on good knowledge of the debris-barrier interaction; however, our understanding of this mechanism is still lacking [11][12][13][14][15][16][17].
In nature, granular flows occur where loose overlying deposits are disturbed, and, in specific conditions, multiple small flow events may occur within a short time period [7]. When a flow reaches

Materials and Methods
For decades, researchers have developed various numerical methods to simulate flow-type landslides, such as rock avalanches and debris flow [26]. Among these methods, the DEM has advantages over the continuum mechanics method [3,27] in modeling granular flows, as it can well address the discrete nature of such flows. The DEM is also widely used to investigate debris-barrier interaction [9,28] because it provides micro-scale information that cannot be obtained in physical modelling [29]. Based on this information, we can determine the relation between the micro-mechanism and macro-behavior [5,17,24].
In DEM simulation, the granular flow is represented by an assembly of spherical or non-spherical particles, and the single-particle motion (translation and rotation) is governed by Newton's second law of motion. By calculating the contact force between the particles and solving Newton's second law of motion, the particles' position and velocity can be obtained step by step, providing a reliable model of the granular flow movement, impact and deposition.
In our study, we adopt a commercial software named EDEM to conduct simulations, in which the micro-contact force could be calculated by the Hertz-Mindlin (no-slip) contact model for its computational efficiency. The model calculates the normal force (F n ) by Hertz's theory and the tangential force (F t ) by Mindlin's no-slip model: here, the subscript n is the normal direction, and t is the tangential direction; K denotes the elastic stiffness constant, and D is the damping coefficient; u represents the overlapping or relative displacement between two particles in contact; v rel is the relative velocity; and µ s is the coefficient of the Coulomb friction. Equation (2) shows that the tangential force is limited by Coulomb's law of friction and accounts for the gross sliding movement between two particles in contact. A rolling torque is adopted to address the rolling friction, expressed as: where µ r is the rolling friction, R is the distance between the contact point and the center of mass, and ω is the unit angular velocity of the particle at the contact point. The numerical model ( Figure 1) is based on the physical test conducted by Jiang and Towhata [14]. The flume is inclined at 40 • , and a 0.4 m-high rigid barrier is set at the end of the flume, perpendicular to the flume base. The flume sidewall is 0.35 m high and 0.3 m wide, and the total flume length is 2.19 m. Initially, we used spherical particles with a particle diameter of 10-20 mm to model the granular soil, forming a rectangular deposition body, 0.15 m high, 0.44 m long and 0.3 m wide. The total mass of the particles is 27 kg, and the bulk density of the initial disposition body is 13.6 kN/m 3 , similar to that in the experiment of Jiang and Towhata [14]. The flume and barrier were represented in the model by a wall element. The DEM input parameters are listed in Table 1, which are obtained by experimental measurement and calibration or from other studies.
Water 2020, 12, x FOR PEER REVIEW 3 of 22 here, the subscript n is the normal direction, and t is the tangential direction; K denotes the elastic stiffness constant, and D is the damping coefficient; u represents the overlapping or relative displacement between two particles in contact; is the relative velocity; and is the coefficient of the Coulomb friction. Equation (2) shows that the tangential force is limited by Coulomb's law of friction and accounts for the gross sliding movement between two particles in contact. A rolling torque is adopted to address the rolling friction, expressed as: where is the rolling friction, is the distance between the contact point and the center of mass, and is the unit angular velocity of the particle at the contact point.
The numerical model ( Figure 1) is based on the physical test conducted by Jiang and Towhata [14]. The flume is inclined at 40°, and a 0.4 m-high rigid barrier is set at the end of the flume, perpendicular to the flume base. The flume sidewall is 0.35 m high and 0.3 m wide, and the total flume length is 2.19 m. Initially, we used spherical particles with a particle diameter of 10-20 mm to model the granular soil, forming a rectangular deposition body, 0.15 m high, 0.44 m long and 0.3 m wide. The total mass of the particles is 27 kg, and the bulk density of the initial disposition body is 13.6 kN/m 3 , similar to that in the experiment of Jiang and Towhata [14]. The flume and barrier were represented in the model by a wall element. The DEM input parameters are listed in Table 1, which are obtained by experimental measurement and calibration or from other studies.   The proposed numerical model was validated by comparing the flow kinematics and impact dynamics of granular flow with that observed in laboratory experiments of Jiang and Towhata [14]. The comparison results are presented in Figure 2, which verifies that our DEM model could generally capture the four stages of flow evolution process reported by Jiang and Towhata [14]. (1) The flow starts from a dam break failure of a rectangular debris deposition with obvious particle motion at flow surface (Figure 2a,a1,b,b1). (2) Then, a mature flow is developed: the flow length is increased, while the flow depth becomes smaller, especially at the flow front where the particles collision is significant (Figure 2c,c1). (3) When the flow front reaches the barrier, some discrete particles become agitated because of impact and rebound. Meanwhile, a portion of flow deposits behind barrier, and the dead zone develops with ongoing debris-barrier interaction (Figure 2d,e,d1,e1). (4) Finally, all particles settle behind the barrier that formed a trapezoidal deposition morphology (Figure 2f,f1). We also compared the Froude number (F r = v/ gH, where v and H denote the average velocity and depth of the flow front, respectively) of the numerical and experiment results. Our result, F r = 7.8, was close to the experimental results (6.7). In addition, the comparison of time-dependent impact force is shown in Figure 3. It is noticed the peak value, residual value and evolution trend of impact force registered in physical tests are reasonably reproduced. These comparisons indicate that the proposed numerical model can be used reliably to investigate the multiple-surge impact behavior of granular flow. More details of calibration process of the DEM model can be accessed in Appendix A.
Based on the calibrated model, we first consider two scenarios: a single-surge event and double-surge impact. The latter scenario follows the methodology of Albaba et al. [21]. Initially, surge1, which has half the volume of the single-surge scenario, was formed. This flow reaches the rigid barrier and forms the dead zone, then another surge (surge2) with a volume equal to that of surge1 is released and interacts with the dead zone and the barrier. The total volume of surge and the slope angle varied in the modeling to investigate the effect of the morphology of the dead zone on the multiple-surge impact mechanism. The simulation program is presented in Table 2.    Based on the calibrated model, we first consider two scenarios: a single-surge event and doublesurge impact. The latter scenario follows the methodology of Albaba et al. [21]. Initially, surge1, which has half the volume of the single-surge scenario, was formed. This flow reaches the rigid barrier and forms the dead zone, then another surge (surge2) with a volume equal to that of surge1 is released and interacts with the dead zone and the barrier. The total volume of surge and the slope angle varied in the modeling to investigate the effect of the morphology of the dead zone on the multiple-surge impact mechanism. The simulation program is presented in Table 2.

Flow-Dead-Zone Interaction Mechanism
Here, we present the flow kinematics of multiple-surge impact and discuss the general interaction mechanisms. Figure 4 shows snapshots of the interaction between surge2 and the dead zone formed by surge1 during a simulation with a slope angle of 40° and a total surge volume of 39.6 L. The diagrams are rotated to lie horizontal. At t = 2.6 s (the initial time is set as the moment of the release of surge1), the front of surge2 is approaching the foot of the dead zone; at t = 2.65 s the front impacts the tail of the dead zone with an average flow velocity of 4.3 m/s, followed by a drastic decrease in the flow velocity by about 90% at t = 2.8 s. At the same time, some particles of the tail of the dead zone are entrained upward at a very low velocity (about 0.5 m/s). As more flow material arrives, the front of surge2 moves along the upper surface of the dead zone, causing minor

Flow-Dead-Zone Interaction Mechanism
Here, we present the flow kinematics of multiple-surge impact and discuss the general interaction mechanisms. Figure 4 shows snapshots of the interaction between surge2 and the dead zone formed by surge1 during a simulation with a slope angle of 40 • and a total surge volume of 39.6 L. The diagrams are rotated to lie horizontal. At t = 2.6 s (the initial time is set as the moment of the release of surge1), the front of surge2 is approaching the foot of the dead zone; at t = 2.65 s the front impacts the tail of the dead zone with an average flow velocity of 4.3 m/s, followed by a drastic decrease in the flow velocity by about 90% at t = 2.8 s. At the same time, some particles of the tail of the dead zone are entrained upward at a very low velocity (about 0.5 m/s). As more flow material arrives, the front of surge2 moves along the upper surface of the dead zone, causing minor deformation in the dead zone mass. This deformation is more obvious at t = 2.95 s. Some finer particles of surge2 are embedded in the larger voids of the dead zone by the pressure of upper flow, thus increasing the flow resistance of surge2. This phenomenon, which was also observed in other studies, serves as a mechanism for particle-size segregation during debris-barrier interaction [30]. At t = 3.1 s, a few particles of surge2 impact the barrier; at same time, the dead zone deformation reaches its maximum. Toward the end of the process (t = 3.25 s and t = 3.4 s), more particles impact the barrier, generating a steady impact force on the barrier. At the final stage, a large portion of the particles are deposited at the foot of the dead zone and on top of it, and only a small portion of the particles have interacted with the barrier.
s, a few particles of surge2 impact the barrier; at same time, the dead zone deformation reaches its maximum. Toward the end of the process (t = 3.25 s and t = 3.4 s), more particles impact the barrier, generating a steady impact force on the barrier. At the final stage, a large portion of the particles are deposited at the foot of the dead zone and on top of it, and only a small portion of the particles have interacted with the barrier.  During the interaction process described above, at least three mechanisms of energy consumption can be identified in addition to the energy dissipated by the barrier and by flume friction: (1) friction and penetration at the interaction face between surge2 and the dead zone (surge1), (2) inelastic deformation of the dead zone and (3) inter-particle interaction within surge2. A more quantitative discussion about these three mechanisms is presented in Section 3.3.
At the final stage, because of the effect of the dead zone during multiple-surge impact, many particles are deposited at the foot of the dead zone; thus, the final deposition morphology is more spread out compared with that of the single-surge impact. Figure 5 shows the deposition morphology of granular flows at the final stage, with flume angles of 40 • , 35 • and 30 • and surge volumes V1, V2 and V3, respectively. For each volume V, the simulation was repeated with two flows, each of volume 0.5·V. The diagrams were rotated to lie horizontally. Overall, the final deposition after the double-surge impact is more spread out than that of the single-surge impact; the difference increases when the flume angle is smaller ( Figure 5). For quantitative comparison, we measured the spreading distance, which is defined as the length of the deposition along the flume base from the bottom of the barrier, normalized using Equation (4): where L ∆ is the additional spreading distance of the deposition, and L s and L m represent the spreading distance of the deposition formed by the single-surge and double-surge impact, respectively. The normalized results are presented in Table 3. As the flume angle decreases, the difference between the spreading distance of the single-and double-surge depositions increases, reaching almost 20%. When the slope is relatively gentle, the flow is slower because of a smaller potential energy, hence the granular flow has lower kinetic energy. The dead zone formed by surge1 has a gentler free surface; therefore, less particles are transported to the top of the dead zone because of the longer energy dissipation path, that is, more particles will be deposited at the tail of the dead zone. The combination of these two effects results in a larger deposition spreading distance when the flume slope is gentler ( Figure 5). The volume of the dead zone (block volume) has a small influence on enlargement of the spreading distance of the deposition. For example, when the flume was steeper (flume angle of 40 • ), for a block volume of 19.8 L (V1), the spreading distance was 6.48% greater for the double-surge case, and when the block volume decreased to 15.84 L (V2), the spreading distance increased by an additional 4.59% (a total of 11.074); however, when the block volume decreased further to 9.9 L (V3), the spreading distance only increased by 10%. The results under all the flow conditions considered in this paper show a similar trend, which suggests that there is an optimal block volume where the difference between the single-and double-surge spreading distance is the largest; this should be further investigated in the future. Water 2020, 12, x FOR PEER REVIEW 9 of 22

Impact Dynamics
The impact force on a barrier is a big concern in engineering practice [3,11,14,28]. In this section, the evolution of the impact force is investigated in detail. Figure 6c shows the simulated time-history curve of the impact force on a barrier for single-surge and double-surge impacts with a slope angle of 40 • and total surge volume of 39.6 L. Under single-surge impact, the granular flow front reached the barrier at t = 0.65 s; as the particles collided with the barrier, the impact force increased sharply and reached a peak at t = 1.14 s, then gradually decreased to a level determined by the earth pressure exerted by the dead zone at the end of the approached flow. Under the multiple-surge impact, the curve of surge1 is similar to that of the single-surge impact, but the peak impact force and residual force are lower by 39.44% and 42.87%, respectively. As shown above, surge2 impacts the dead zone formed by surge1 and flows up the free surface of the dead zone ( Figure 4) until reaching the barrier. A large portion of the impact energy is dissipated, leaving only a small impact force acting on the barrier (about 1.29% of the peak force generated by the single-surge impact). However, after surge2 impacts the dead zone (t = 2.65 s; Figure 6c), the force on the barrier also increases sharply, peaking within a short time (about 0.58 s) and then decreasing again. This indicates that the impact energy of surge2 is transformed to an impact force on the barrier mainly through the dead zone. In this process, the dead zone serves as a cushioning layer that protects the barrier from the high-impact force. Figure 6 also gives the influence of slope angle on debris impact force. With a lower slope angle, the time-history curve of impact force is gentler.  In our quantitative analysis, we calculated the total force generated by surge1 and surge2 and compared this with the force exerted by the single-surge impact. The peak force reduction was calculated and normalized by the force of the single-surge impact. The results are presented in Table  4 and reveal the cushioning effect. The highest peak force reduction was 32.65%, which demonstrates that the dead zone serves as an effective cushioning layer and provides the barrier with a stabilizing effect against debris flow. Table 4 also shows the effect of the slope angle and block volume on the peak force reduction, with a similar effect to that shown in Table 3. For example, when the block volume is 19.8 L, as the slope angle decreases from 40° to 30°, the peak force reduction increases from 7.43% to 32.65%. At a slope angle of 40°, when the block volume decreases from 19.8 L to 9.9 L, the In our quantitative analysis, we calculated the total force generated by surge1 and surge2 and compared this with the force exerted by the single-surge impact. The peak force reduction was calculated and normalized by the force of the single-surge impact. The results are presented in Table 4 and reveal the cushioning effect. The highest peak force reduction was 32.65%, which demonstrates that the dead zone serves as an effective cushioning layer and provides the barrier with a stabilizing effect against debris flow. Table 4 also shows the effect of the slope angle and block volume on the peak force reduction, with a similar effect to that shown in Table 3. For example, when the block volume is 19.8 L, as the slope angle decreases from 40 • to 30 • , the peak force reduction increases from 7.43% to 32.65%. At a slope angle of 40 • , when the block volume decreases from 19.8 L to 9.9 L, the peak force reduction first increases from 7.43% to 17.84% and then decreases to 12.51%. All the results exhibit a similar trend, which again suggests an optimal block volume that causes the largest peak force reduction. This aspect should be further investigated in the future.

Energy Evolution
Modeling the energy evolution can offer insights on the multiple-surge impact behavior. The evolution of the normalized total energy of the debris surge was simulated for a slope angle of 40 • and total surge volume of 39.6 L ( Figure 7). As soon as the debris is released, the total energy begins to decrease, initially slowly, then at a higher rate. As the flow motion ends, the energy curve flattens at a low level and maintains that level, which represents the total energy of the dead zone. The single-surge and double-surge impacts have a similar total energy evolution trend. However, during the single-surge impact, about 85.44% of the total energy is dissipated, while for surge1 and surge2, 89.31% and 79.85% of the energy is dissipated, respectively. In the interaction between the surge and the dead zone, most of the debris surge particles are deposited at the foot of the dead zone and above it in the flume (Figure 4). This results in a lower energy loss for surge2 because the dead zone of surge2 has a higher potential energy at the final stage. Therefore, surge2 cannot fully interact with the flume and the barrier, hence, more energy can be dissipated by particle-particle interaction as the surge-dead-zone interaction increases during the double-surge impact. To quantitatively assess this process, we directly measured the energy absorbed by the particle-particle interaction. The energy loss caused by inter-particle interaction could be directly computed in DEM simulation because the particle energy loss through normal damping, tangential damping or shear friction and rolling friction could all be automatically registered during particle contact. For the single-surge impact, particle-particle interaction accounts for 40.75% of the total energy loss, while for the double-surge impact, it accounts for 60.14% of the energy loss. This indicates that the dead zone facilitates the dissipation of almost 20% of additional impact energy.
In Section 3.1, we identified three energy consumption mechanisms for multiple-surge impact during surge-dead-zone interaction: the friction and penetration at the interaction face between surge2 and the dead zone (E ∆S1−2 ), the inelastic deformation of the dead zone (E ∆S1−1 ) and the inter-particle interaction within surge2 (E ∆S2−2 ). In our detailed analysis of the energy consumption, we focused on the interaction process. Thus, the start time (t = 0) was set as the time when surge2 reached the tail of the dead zone formed by surge1 (originally t = 2.65 s), and the time-history curve of the total energy loss (E ∆T ) was recalculated. The amount of energy absorbed by the three identified mechanisms was measured and normalized by the total energy loss (E ∆T ). The results are presented in Figure 8 with slope angles of 30 • , 35 • , 40 • and a total surge volume of 39.6 L. For the slope angle of 40 • , the total energy loss increases with time, reaching 65.78% at the final stage. The inelastic deformation of the dead zone (E ∆S1−1 ) only accounts for 3.59% of the total energy loss, and the friction and penetration at the interaction face between surge2 and the dead zone (E ∆S1−2 ) account for 1.94%. However, the surge2 inter-particle interaction accounts for 62.99% of the total energy loss. The remaining energy is absorbed by the flume and barrier friction. Using the same method, we calculated the energy loss at the other flow conditions (Table 5), and the results show a very small variation. Thus, we conclude that during multiple-surge impact, the enhanced inter-particle interaction within the subsequent surge is the main process of energy consumption.
surge-dead-zone interaction increases during the double-surge impact. To quantitatively assess this process, we directly measured the energy absorbed by the particle-particle interaction. The energy loss caused by inter-particle interaction could be directly computed in DEM simulation because the particle energy loss through normal damping, tangential damping or shear friction and rolling friction could all be automatically registered during particle contact. For the single-surge impact, particle-particle interaction accounts for 40.75% of the total energy loss, while for the double-surge impact, it accounts for 60.14% of the energy loss. This indicates that the dead zone facilitates the dissipation of almost 20% of additional impact energy. In Section 3.1, we identified three energy consumption mechanisms for multiple-surge impact during surge-dead-zone interaction: the friction and penetration at the interaction face between surge2 and the dead zone ( ∆ ), the inelastic deformation of the dead zone ( ∆ ) and the interparticle interaction within surge2 ( ∆ ). In our detailed analysis of the energy consumption, we focused on the interaction process. Thus, the start time (t = 0) was set as the time when surge2 reached the tail of the dead zone formed by surge1 (originally t = 2.65 s), and the time-history curve of the total energy loss ( ∆ ) was recalculated. The amount of energy absorbed by the three identified mechanisms was measured and normalized by the total energy loss ( ∆ ). The results are presented in Figure 8 with slope angles of 30°, 35°, 40° and a total surge volume of 39.6 L. For the slope angle of 40°, the total energy loss increases with time, reaching 65.78% at the final stage. The inelastic deformation of the dead zone ( ∆ ) only accounts for 3.59% of the total energy loss, and the friction and penetration at the interaction face between surge2 and the dead zone ( ∆ ) account for 1.94%. However, the surge2 inter-particle interaction accounts for 62.99% of the total energy loss. The remaining energy is absorbed by the flume and barrier friction. Using the same method, we calculated the energy loss at the other flow conditions (Table 5), and the results show a very small variation. Thus, we conclude that during multiple-surge impact, the enhanced inter-particle interaction within the subsequent surge is the main process of energy consumption.    The initial time is readjusted to t = 2.65 s when the front of surge2 reaches the tail of the dead zone, which means we focus only on the interaction process.

Prediction of Multiple-Surge Impact Forces
Calculating debris flow impact forces is an important part of structure design and risk mapping [3,11,14,28,31]. For decades, many impact force prediction models have been proposed for engineering purposes [11], and the hydraulic model is commonly used because it is simple and easy to use [11,14,31]. Generally, the maximum flow velocity and depth are used to predict the peak force on the barrier [14,31]. However, this strategy ignores the complex debris-barrier interaction mechanism, such as the dead zone formation, which has been shown to be an important factor when addressing the impact forces on the barrier [13,23,31]. Therefore, the impact force exerted on a barrier by dry granular flows may be overestimated, leading to needless construction costs [23]. Koo et al. [23] proposed a simplified load model based on rigid body motion, assuming that the dead zone formed during debris-barrier interaction is triangular shaped with an inclined free surface. When the second flow reaches the tail of the dead zone, the velocity and depth are recorded before the flow moves up along the free surface of the dead zone and finally impacts the barrier. During this process, Koo et al. [23] only considered the kinetic energy loss caused by interface friction and gravity, which is similar to rigid body motion. In our simulation (Section 3), the interface friction only accounts for a small fraction of the energy loss (~2%), and most of the energy is dissipated by the inter-particle interaction within the moving flow; the inelastic deformation of the dead zone also dissipates some of the energy (~5%). Thus, predicting the impact force on a barrier is complex and requires further investigation.
While multiple-surge impact is a more complex process, it is more realistic than single-surge impact; however, very few models account for the multiple-surge impact forces. Tan et al. [7] studied multiple-surge impact on a flexible barrier; however, the impact force was calculated separately and was based on single-surge impact. Albaba et al. [28] proposed an analytic solution based on the granular jump theory and depth-averaged method for calculating the time-history curve of granular flow impact on a rigid barrier. This method may be suitable for predicting the force of multiple-surge impact because it is strongly time-dependent.
In the model of Albaba et al. [28], the total impact force on the barrier is composed of the static force exerted by the dead zone behind the barrier and the inertia force of the flowing mass: where C and K are time-dependent factors calculated using Equations (6) and (7); ρ p is the particle density; φ 1 and u 1 are the depth-averaged solid volume fraction and flow velocity of the subsequent flow, respectively; h 1 is the flow depth; w is the flume width; φ 2 is the solid volume fraction of the dead zone (0.6 is suggested by Albaba et al. [28]); g is gravity; and α is the flume angle.
here, h 2 denotes the granular jump height as in Equation (8); F r is the Froude number; λ φ is the ratio of φ 2 to φ 1 ; γ is a shape coefficient; µ dz is the Coulomb friction coefficient of the dead zone (Albaba et al. [28] used µ dz = 0); and l dz and h dz are the length and height of the dead zone, respectively.
Albaba et al. [28] suggested that when the slope angle (α) is lower than the inter-friction angle (arctan µ s ), the dead zone formed behind the barrier has a rectangular shape; thus, the shape coefficient (γ) equals 1, h dz and h 2 are assumed equal, and l dz can be estimated from the previous time step, as formulated by Equation (9): The impact forces on the barrier in the cases of multiple-surge impact with a slope angle of 30 • , 35 • and 40 • (arctan µ s = 53 • , hence γ = 1) are back-analyzed based on the model of Albaba et al. [28]. Table 6 shows that the model prediction is not satisfactory, especially when the flume is steep. For example, under α = 40 • , the peak force of the first surge impact is underestimated by~31.84%, and the total peak force during the second surge impact is underestimated by~14.71%. However, under α = 30 • , the model performs well with the peak force and main evolution being reasonably well captured. The discrepancy may be explained by the inaccurate representation of the shape of the dead zone. Based on Albaba et al. [28], the dead zone formed by the granular flow was modeled as rectangular; however, our numerical modeling of the dead zone showed an almost trapezoidal shape with a steep free surface (Figures 4 and 5). This misrepresentation of the dead-zone shape effect is the main flaw of this model. To overcome this problem, we used a larger shape coefficient (γ) ( Table 6), which improved the accuracy of the results. Additionally, when the second surge impacts the dead zone formed by the first surge, it causes deformation of the dead zone (Figure 2), which may account for the overestimation of the prediction at the second stage, especially for high γ values (Table 6). In addition, the momentum transfer effect of the dead zone has also not been properly addressed.

Engineering Implications and Limitations
Our multiple-surge simulation presents a more realistic case than the single-surge case when addressing debris flows or debris avalanche barrier design in engineering practice because multiple-surge impact is more commonly encountered in nature [7,21]. As multiple-surge impact is not normally incorporated into barrier design, our results offer preliminarily insights for structural engineers. Overall, for flow masses with the same total volume, a barrier will be safer under multiple-surge impact than under single-surge impact because the deposition spreading distance is extended by 6%~20% and the impact force is reduced by 6%~30%. These findings are useful when considering barrier design that utilizes only single-surge impact theory because they provide assurance that the design has an adequate safety margin in the more common multiple-surge events. The multiple-surge impact mechanism discussed in this paper is essential for establishing the design strategy for barrier structures planned to withstand multiple-surge impact, as these should be carefully addressed in the impact force prediction model. For example, to improve the performance of the original depth-averaged model, two factors should be addressed: the dead-zone shape effect and the momentum transfer effect of the dead zone. As shown in Figure 5 and Table 3, a portion of the material deposited behind the barrier can be retained, especially when the slope is gentle, as the impact force decreases and the capture capacity of the barrier increases. However, one should be extremely careful when using our conclusions when the slope is steeper, because under such a condition, the retaining capacity of the barrier is small and overflow may be more likely to occur. Additionally, our study simulates dry granular flows; thus, the results may be more suitable for arid regions and where small flow events are dominant. It is important to note that our research aims to highlight the importance of the multiple-surge impact process in debris-barrier interaction and to present some preliminarily guidelines for barrier structure design, although a comprehensive design strategy needs to consider multiple aspects. In future research, more experiments, including numerical simulation and physical modeling, are needed to obtain a deeper understanding of multiple-surge impact and establish a reliable and verified design method.
Being a preliminary study, we acknowledge that there are still some open issues that have not been resolved. For example, our analysis showed interesting trends in the effect of the flow-mass volume on the deposition morphology and peak force on the barrier, which requires further study. Adding viscous flows to our model is another research path which should be investigated in detail because the impact behavior would be considerably different from the friction flow pattern in the current model [13]; viscous flows may develop distinct energy distribution routes and will need specific barrier design strategies.

Conclusions
When granular flows encounter a rigid barrier, the debris-barrier interaction is a crucial but challenging issue in landslide mitigation measurements design and vulnerability analysis of buildings. Current research on this issue focuses on single-surge impact conditions while the multiple-surge impact mechanism remains unclear, although multi-surge events are more common in engineering practice. In this study, we present some DEM-based numerical simulations that offer preliminary insights into barrier design against multiple-surge impact. The main conclusions are as follows: (1) For cases in which the total volume of flow material is the same, the barrier is safer under double-surge impact than under single-surge impact, as the deposition spreading distance is extended by 6%~20% and the impact force is reduced by 6%~30%.
(2) Compared with single-surge impact, the energy dissipated by particle-particle interaction during double-surge impact is enhanced by 20% because of the effect of the dead zone. (3) After the first surge is deposited and the dead zone is formed, the second surge arrives, and a large portion of the energy is dissipated when the second surge interacts with the dead zone. The second flow then impacts the barrier, generating a negligible impact force. The momentum of the subsequent surge is transmitted to the barrier mainly through the dead zone. Thus, the dead zone serves as an effective cushioning layer that protects the barrier. (4) Multiple-surge impact is a complex process. During the interaction between subsequent surges and the dead zone, three energy consumption mechanisms were identified: (1) friction and penetration at the interaction surface between surge2 and the dead zone, (2) inelastic deformation of the dead zone and (3) inter-particle interaction within surge2. These three mechanisms account for about 2%, 3% and 63% of the total energy loss, respectively. (5) The depth-averaged model based on granular jump theory may be a promising solution for predicting the multiple-surge impact force on the barrier. However, the performance of the original version is less satisfactory due to the dead zone shape effect, the cushion effect and the dead zone transition, leading to deformation not being reasonably addressed.