Application of Smoothed Particle Hydrodynamics (SPH) for the Simulation of Flow-like Landslides on 3D Terrains


 Flow-type landslide is one type of landslide that generally exhibits characteristics of high flow velocities, long jump distances, and poor predictability. Simulation of it facilitates propagation analysis and provides solutions for risk assessment and mitigation design. The smoothed particle hydrodynamics (SPH) method has been successfully applied to the simulation of two-dimensional (2D) and three-dimensional (3D) flow-like landslides. However, the influence of boundary resistance on the whole process of landslide failure is rarely discussed. In this study, a boundary algorithm considering the friction is proposed, and integrated into the boundary condition of the SPH method, and its accuracy is verified. Moreover, the Navier-Stokes equation combined with the non-Newtonian fluid rheology model was utilized to solve the dynamic behavior of the flow-like landslide. To verify its performance, the Shuicheng landslide event, which occurred in Guizhou, China, was taken as a case study. In the 2D simulation, a sensitivity analysis was conducted, and the results showed that the shearing strength parameters have more influence on the computation accuracy in comparison with the coefficient of viscosity. Afterwards, the dynamic characteristics of the landslide, such as the velocity and the impact area, were analyzed in the 3D simulation. The simulation results are in good agreement with the field investigations. The simulation results demonstrate that the SPH method performs well in reproducing the landslide process, and facilitates the analysis of landslide characteristics as well as the affected areas, which provides a scientific basis for conducting the risk assessment and disaster mitigation design.


30
The occurrence of landslides is usually introduced by earthquakes or heavy rainfall, and are 31 always accompanied by a large number of casualties and extensive damage (Cascini et al. 2012). 32 According to the classification of landslides, there is a term "flow-like landslides", which has 33 the characteristics of fast flow velocity and long run-out distance, that causes severe damage 34 (Mcdougall and Hungr 2004;Huang et al. 2012a; Pastor et al. 2014). Therefore, great attention 35 should be paid to these flow-like landslides to study the flow mechanism and the sliding 36 characteristics. Due to the vast scale, unpredictability, short duration, and destructive character 37 of landslide disasters, it is quite difficult and dangerous to record the process in real time.

38
Numerical methods can simulate the flow behavior of such landslides, and can also be applied 39 to the investigation and evaluation of landslides. Most numerical studies on the landslides are 40 based on solid mechanics analysis methods, such as the limit equilibrium, finite element method 41 (FEM) (Chatterjee and Krishna 2018) and the fracture mechanics method (Zhang 2020). These 42 methods focus on the slope stability analysis, which analyzes the mechanism of landslide 43 damage from the stress-strain perspective. However, these methods are not able to simulate the 44 large deformation state of the landslide, and thus can not reproduce the whole landslide process. through appropriate force-displacement relationships that relate to the overlap and contact 53 forces that can occur between particles. However, sufficiently small time steps are needed to 54 avoid the numerical instability, leading to high costs when simulating the problems containing 55 a large number of elements. To reduce the computational effort, the grain size in a typical DEM 56 simulation is typically significantly less than that in a real problem. Another viable approach is 57 discontinuous deformation analysis (DDA) ( which can be applied in the complex geometries. However, the current simulation of non-88 Newtonian fluids lacks the treatment of fluid-solid boundaries (Dai et al. 2017). The literature 89 suggested that the application of free-slip boundary conditions to treat the fluid-solid boundary 90 will make the non-Newtonian fluid unable to obtain resistance from the solid wall boundary in 91 the simulation, thus making the simulation and experiment inconsistent (Cao et al. 2019).

92
In this study, a boundary treatment method is proposed that considers the friction based on 93 the SPH method. For the SPH method, the entire domain is represented as a set of randomly distributed particles 106 with no connection between them. Two main steps are required to obtain the fluid governing 107 equations in SPH form, namely, the kernel approximation and the particle approximation. 108 Firstly, approximation of the field function and its derivative are derived from an integral 109 representation method by using a smoothing kernel function. Secondly, these approximations 110 are replaced by the sum of all neighboring particles in the support domain. It is the kernel 111 approximation and particle discretization that make the SPH method work well for free surface 112 flows and large deformations (Liu and Liu 2010; Violeau and Rogers 2016).

113
The integral representation of a function () f x may be given as the following kernel 114 approximation using the smoothing technique: 115 where () f x is a field function of the three-dimensional position vector x , and ( the smooth function. The smoothing length h that determines support domain  . The discrete 118 form of a function at particle i can eventually be described as: 119 (2) 120 where j m and j  are the mass and density of particle j, respectively, and indicates the smoothing function. In the same manner, the particle approximate form for the 122 spatial derivative of the function can be formulated as: 123 There are several kernel functions that are widely utilized. Typically, higher-order kernel 125 functions often enhance the computing accuracy, and substantial time is required. The cubic 126 spline kernel, which belongs to the B-spline family, is one of the kernels that is commonly 127 employed in SPH method (Liu and Liu 2010) and is given by: 128 equations can be written as follows: 132 − ν ν ν indicates the velocity difference between particle i and particle j; σ represents the 135 total stress tensor and g is the gravity acceleration, respectively. It is assumed that the 136 temperature is constant throughout the simulation, so the energy equation is not presented. 137 Due to the density approximation, the pressure field in the SPH formulation for weakly 138 compressible fluids generally exhibits instability and numerical noise. Particularly, in the 139 present case, the noise of the pressure field may be critical. In this study, the delta-SPH method 140 (Marrone et al. 2011) was adopted, which implemented by adding a diffusive term in the 141 continuity Eq.(5) and，can be written as: 142 where  is the diffusive coefficient, which is usually adopted as 0.1; 0 c is the speed of sound at 144 reference density; and ij  is written as: where h  is a small constant and is usually adopted as 0.01.

147
Regarding the momentum Eq. (6), the total stress tensor σ can be written as the summation 148 of an isotropic component pressure p and a viscous shear stress component τ . To remove the 149 numerical oscillations, the Monaghan type artificial viscosity ij  was adopted in this study, 150 which is most widely used in the current SPH simulation. After considering the artificial 151 viscosity, the momentum equation can be written as: 152 where  and   are two constants and are commonly set as 0.1 and 0.01, respectively; Detailed

155
parameter definitions can be found in reference (Monaghan 1992).

156
In this work, the temperature of landslides is considered as a constant. Therefore, the 157 energy equation is not considered. In the typical SPH method for resolving the compressible 158 flows, the particle motion is governed by the pressure gradient, whereas the particle pressure is 159 determined using an appropriate Equation of State (EoS). 160 whereis commonly set to be 7.0; 0 c is the speed of sound and is commonly set at least 10 times 162 the maximum fluid velocity for the sake of keeping the density fluctuation within 1% which 163 can be calculated by ( )  is the Bingham viscosity and y τ is the yield stress, respectively. 184 Assuming that under the condition of shear rate, approaches an infinite value, therefore, the direct application is not available. To solve this 186 problem, existing studies employed a variety of strategies to avoid numerical divergence, i.e., that the numerical solution's accuracy is satisfactory as the viscosity at a low shear rate is 1000 206 times greater than that at a high shear rate.

211
Notice that the advantage of the Cross model over the Bingham model is that the effective 212 viscosity is a continuous variable, and the numerical instability is avoided, as shown in Fig. 1.

214
2.3 Boundary condition 215 The boundary treatment has been a main challenge in SPH simulation and affects the calculation 216 accuracy as well as the efficiency. It is vital to select an adequate boundary condition that 217 represents the effect of solid boundary in order to closely describe the dynamic mechanism of 218 landslides, as shown in Fig. 2.

219
A simple and convenient boundary algorithm was proposed that can consider the boundary 220 friction. In this work, two types of boundary conditions need to be set. Firstly, when a particle 221 approaches to a solid boundary, the kernel function will be truncated, and an error in the solution 222 result will occur. To address this issue, a generalized wall boundary condition by Adami ( Through the above two steps, the no-slip boundary conditions considering friction on the 233 solid boundary was successfully set. When calculating the viscous forces on the boundary, both 234 the dynamic viscosity coefficient can be obtained by interpolation from the fluid, or a fixed 235 viscosity coefficient can be set.  vertical section of L-box has a capacity of 2.5 L. The sliding door was lifted up quickly in 2-3 257 seconds. After the door was opened, the flowing distances of the mortar sample were recorded 258 with a video camera at a certain time interval. One of the cases was chosen for the simulation, 259 and the parameters are given in Table 1. 260 Table 1 The input parameters of numerical simulations The simulation process is presented in Fig. 3. After the simulation begins, the fresh mortar 262 column collapses under gravity and moves along the bottom of the box. Due to the boundary 263 resistance and its viscosity, the granular flow stopped flowing at the time of 6 s. Moreover, the 264 evolution of the collapsed fresh mortars front before stopped was quantitatively studied and 265 compared with the experimental data. Comparing the simulation results with the test data ( Fig.  266 3), the peak distance and residual deformation of the flow are well captured by the SPH 267 simulation. Fig. 4 shows the numerical and experimental results of the change in L-flow 268 distance over time. To clarify the effect of slippage resistance, the numerical results without 269 considering the slippage resistance are also shown in Fig. 4. The numerical results are not in 270 good agreement with the experimental results when slippage resistance was not considered.

271
This is because that the use of free-slip boundary conditions will make the non-Newtonian fluid 272 unable to obtain sufficient resistance and reaction force from the solid boundary during the 273 process. This section shows that the proposed boundary treatment can simulate a fluid-solid 274 boundary with friction effectively. The Shuicheng flow-like landslide event was triggered by intense rainfall in July 2019 (Fig. 5  283  a). A total of three heavy rainfall events occurred in the week before the occurrence of landslides, 284 with a cumulative rainfall of 288.9 mm. Under the influence of continuous heavy rainfall, the 285 landslide begun to lose stability and rapidly slid along the slope at high speed. After roughly 286 700 meters of landslide movement, two diversions were formed after hitting a small ridge and 287 subsequently impacting and burying the residential buildings on the lower part of the slope. The where v C is the volume concentration. According to the investigation report, the best-fitted v C is 306 44 %, resulting in a dynamic viscosity of 255 Pa s  . 307 According to the typical cross-section in Fig. 5c, a two-dimensional SPH model was 314 established to simulate the flow behavior of Shuicheng flow-like landslide. 1318 particles were 315 used to represent the landslide and 2678 particles were used to represent the boundary. The 316 initial velocity of the particles was set to zero. After the slope failure, the landslide flow mass 317 particles start to slide downward under the action of gravity. The proposed boundary treatment 318 method considering the friction were utilized, while the boundary particles remained stationary 319 throughout the entire simulation. 320 According to the simulation results, the whole process duration was about 90s. Fig. 6  321 presents the slope configurations at different time points, reflecting the movement of landslides. 322 The particles slipped down from the top of the landslide and then moved along the steep slope 323 by gravity. Finally, these particles reached an equilibrium state and were deposited in the 324 lowlands. The overall performance of the sliding was accelerated motion in the period 0-30 s 325 and decelerated motion after 30 s. Landslides usually have a long distance, and the measurement and selection of various 329 parameters may vary greatly in different locations. Therefore, it is meaningful to conduct a 330 sensitivity analysis to determine how various rheological parameters affect the simulation 331 results. Table 2 lists the eight calculation cases with different rheological parameters. To 332 explicitly quantify the differences of the sensitivity of different parameters v  , was evaluated 333 using the following equation: 334 where j v  is the velocity difference between case j and case 1, and N=18 represents the number

336
of velocity at an interval of 5 s. The results of all the eight cases are presented in Table 2 and 337 Fig.7, which shows that the internal friction angle has more influence on the computation 338 accuracy in comparison with the cohesion of landslides. 339 Since the 2D simulation cannot reflect the diffusion and lateral contraction, 3D simulation on 345 the real complex terrain is necessary. In this section, the 3D terrain is generated from the original 346 topographic map at a scale of 1/2000 and homogenizes the mesh. The previously created mesh 347 is then transformed into a sequence of boundary particle particles. Similarly, the landslide is 348 discretized into a succession of particles each with its own set of properties. In this simulation, 349 the particle distance is set to 7.5 m, resulting in 12825 boundary particles and 7010 fluid 350 particles. The strength parameters adopted in the 3D simulation are the same as in the 2D 351 simulation (Table 3). Based on this model, Shuicheng landslide was numerically simulated in 352 3D terrain, and the results are shown in Fig. 8  The results suggest that the front flow takes about 28 s to travel nearly 520 m to reach the 358 isolated island. The front flow velocity reaches a maximum with an average velocity of 29.8 359 m/s. Due to the good stability of the soil near the isolated island, it diverts the debris to the sides 360 and forms a divergence. The diversion of the stable slope protects the lower residential houses 361 and creates an "island of safety". The existence of isolated islands has dramatically slowed 362 down the landslides. Therefore, the front flow remains stable at about 90 s and is deposited in 363 the lowlands. Eventually, the estimated total volume of deposition is 1.110 6 m 3 , which is less 364 than the field investigation value 1.9110 6 m 3 . The reason might be that the simulation process 365 does not take into account the entrapment process. In addition, the simulated landslide paths 366 and deposited areas matched well with the results of field investigations. The topographic data, 367 in this case, were obtained from the 1/2000 scale topographic map, and it is considered that 368 more accurate topographic data would help to obtain better simulation results. 369 The displacement and velocity time histories of the landslide front and rear were 370 investigated to gain a better understanding of the landslide characteristics in the present 371 simulation. combination of natural and artificial reasons, which makes the simulation more challenging. In 384 addition, it can be seen that the flow's pressure field is formed smoothly. Fig. 11 depicts the 385 evolution of the flow-like landslide pressure field over time. As previously stated, delta-SPH 386 successfully corrects pressure field noise, particularly for the high viscosity fluids. The landslide material is discretized into a succession of SPH particles of varying diameters 396 in this model. The greater the number of particles and the narrower the particle spacing, the 397 greater the computing time required. It should be mentioned that, in the current code, the SPH 398 module is accelerated based on Open Multi-Processing (OpenMP). In theory, the higher the 399 number of particles in the SPH, the more accurate the simulation will be, but a balance needs 400 to be reached with the needed calculation time. Fig. 12 depicts the relationship between the 401 program runtime and thread count. The results show that the computational efficiency of the 402 proposed SPH model improves as the number of threads increases. Although the OpenMP code 403 is more computationally efficient than the serial code on the CPU, the computation time 404 consumption is large when more particles exist. shearing strength parameters has more influence on the computing accuracy in comparison 419 with the coefficient of viscosity. 420 3. A 3D model of Shuicheng landslide was constructed and computed that corresponded to the 421 site conditions. In terms of the solid volume and deposition area, the computed results are 422 generally consistent with the field investigations. In addition, the characteristics of the front 423 and the rear flow velocity and zone of influence of the flow-like landslide were analyzed, 424 which will help to the mitigation or countermeasure design work and offer evidence for the 425 hazard assessment. 426 4. The approach provides an effective tool for studying the dynamic behaviors of landslides. 427 However, the current study does not take into account the entrainment process of flow-like 428 landslides. Large-scale practical problems usually require massive particles, and future 429 research should be directed to GPU-accelerated modules.

431
Acknowledgments 432 This numerical code has been developed based on PySPH, which is an open framework for 433 Smoothed Particle Hydrodynamics (SPH) simulations. The authors are also appreciate for the 434 valuable discussions with Zhang Weijie (College of Civil and Transportation Engineering, 435 Hohai University, China) and Shi Wenbin (College of Civil Enginnering, Guizhou University, 436 China) on topics of landslides. 437