3D Numerical Simulation of Landslides for the Full High Waste Dump Using SPH Method

Waste dump that is generally composed of a large number of loose geotechnical materials is prone to landslides under external loads. In this work, the smoothed particle hydrodynamics (SPH) method combined with the Mohr–Coulomb model is used to study the dynamic characteristics of the landslides that occurred in the waste dump during the failure process. A benchmark test is firstly conducted to verify the effectiveness of the SPH model. Then, taking the Nanfen full high waste dump with a vertical drop of 300 m in Benxi City, China, as an example, the most dangerous section is selected to establish the SPH numerical model for the waste dump landslides, and the overall dynamic process of the landslides is simulated. The simulation results show that the particles in the middle and upper of the slope have larger potential energy, and their sliding distance is larger. On the contrary, the sliding distance of particles in the lower of the slope is smaller. The particles' sliding distance decreases as the depth increases in the vertical direction of both shoulder and middle of the slope. The particles undergo a process of first acceleration and then deceleration. The sliding distance is in good agreement with the field survey result, and the landslides profile is basically consistent with the actual one. The sensitivity analysis of different particle numbers shows that the number of particles has little effect on the numerical results. The SPH method can vividly reproduce the dynamic process of the landslides in the full high waste dump. The evaluation of the sliding characteristics and risk impact range can provide the key parameters and basis for the prevention and control of the landslides in the full high waste dump and ensure the safety of the mine life cycle.


Introduction
As the waste dumps increase in scale and height, their stability can be reduced due to rainfall, earthquake, and human activities. Sometimes, these harmful factors can lead to serious disasters such as landslides and debris flows, which seriously threaten the safety of mining. Currently, there is little work dealt with the mechanism and dynamic process of dump landslides. With the development of mining construction and deep open-pit mines, the design and operation of waste dump require understanding of the potentially dynamic failure mechanism of landslides to ensure the safety of waste dump operation.
At present, the study of waste dump is rare. Most of the studies are focused on the evaluation of the stability of the dump slope, but few studies focus on the landslides or debris flow of the dump. For the evaluation of the safety stability of waste dumps, the mostly used methods are the limit equilibrium analysis techniques [1][2][3][4]. Due to the action by selfweight, external loads, earthquake, water, and environmental effect, the stability of the waste dump can be significantly reduced. Even though the safety factor of dump slope could be obtained by solving the limit equilibrium equations of the dump system, such safety factor cannot be used to describe the evolution of the slope, which usually involves large deformation of the dump materials. Over the years, the finite difference methods were also applied to analyze the stability of waste dump [5][6][7][8], but most of them were used to study the static problems, and only the evolution of potential sliding surface or plastic zone was obtained. Due to the large scale and height of the dumps, researchers usually use physical model tests to study their stability [9]. However, the physical model test has the problem of size effect. Numerical simulation has become the most common method for studying the stability of the dump slope.
In the past few decades, people have developed numerical models with different accuracies. At first, people used the finite element method (FEM) to study structural problems. Later, when the fracture problem was designed, the extended finite element method (XFEM) [10] was proposed. Another method, the efficient remeshing algorithm, has obtained a good result in the study of crack propagation [11]. e discrete element method (DEM) is widely used in geotechnical engineering, geological engineering, and other fields to deal with discontinuous problems. It can simulate the contact and collision process of particles or blocks. Certainly, the discrete element method can be used to deal with large deformation problems [12], but due to the limit of computing scale, it is difficult to simulate the problems of real engineering scale currently. Another method in the numerical field is the meshless method. e cracking particles method (CPM) [13], the peridynamics (PD) method, and the dual-horizon peridynamics (DH-PD) method [14][15][16] are also widely used in the field of fracture mechanics. ey have made important progress in the study of crack growth and propagation. ese meshfree methods may be applied to the study of large deformation problems such as landslides in the future.
Landslides are a hot research topic in the field of geotechnical engineering. Now researchers usually use the computational fluid dynamics (CFD), material point method (MPM), and particle finite element method (PFEM) (see [17][18][19]) to study debris flow, landslides, and other flowtype problems in fluid mechanics, geotechnical mechanics, and other fields, and a series of results have been obtained.
In recent years, the smoothed particle hydrodynamics (SPH) [20,21] method has been employed to simulate landslides. e extensive application of the SPH method in the field of landslides has developed many constitutive models. e common ones are the fluid-based Bingham model [22][23][24]. e SPH constitutive models based on solids are the Mohr-Coulomb model [25], the elastic-plastic soil model based on the Drucker-Prager yield condition [26], and so on. In solving the SPH instability problem, the Lagrangian gradient smoothing method [27] was developed.
One of the most significant advantages of the SPH method is that it is flexible to accommodate the large deformation of landslides and debris flows. As a meshfree and pure Lagrangian method, it can avoid the interface problem between Euler mesh and material in Euler description so that it is especially suitable for solving dynamic large deformation problems. In addition, the SPH method has high computational efficiency and occupies less computer memory, so it can be competent for engineering scale numerical model. In this study, the vertical drop of waste dump is more than 300 m, so the calculation scale is very large. For other numerical methods, it may take more calculation time and computer memory. e SPH method can run efficiently and obtain satisfactory results.
In the present work, we employ the SPH method combined with the Mohr-Coulomb model to investigate the overall dynamic process of the full high waste dump landslides of Nanfen Open-Pit Mine in China. e SPH method and the elastic-plastic constitutive model are both seldom used in waste dump landslides simulation in the literature. Meanwhile, field investigation, satellite data, and geological data are combined to establish the three-dimensional waste dump landslides model with the most dangerous section of the dump slope, which is made in the same size as the real one. A comprehensive understanding of the movement characteristics of the landslides, such as sliding distance, velocity, and profile, has been obtained. Besides, the sensitivity analysis of different particle numbers has been carried out. All this information can help identify the potentially dangerous areas of the dump slope. Additionally, they can provide disaster prevention and mitigation strategies to ensure the long-term safe operation of the waste dump. Our work provides a good case for the engineering research of the landslides in the full high waste dump with a vertical drop of more than 300 m.

The SPH Method
e SPH's essence is to discrete the computational domain into a certain number of particles; each of them occupies a certain mass and volume of the material and carries the physical parameters such as density, velocity, acceleration, and pressure/stress [28].
e SPH method is established on interpolation theory, which involves kernel approximation and particle approximation. A particle is approximated by all the other particles weighted through a kernel function in the influence domain. Figure 1 illustrates the principle of the SPH interpolation. W is the smoothing kernel function, i is the concerned particle, and j is the neighbor particle.Ω defines the influence domain, and κh is the radius of the support domain, which determines the accuracy of the problem being studied, where κ is a constant and associated with the smoothing kernel function. e particles interplay each other within the influence domain and are independent of each other beyond the domain. A larger smoothing length or influence domain generates a smoother or more continuous behavior.
In the SPH method, the kernel approximation of a function f(x) on the domain can be defined as follows: where h is the smoothing length, and W(x − x ′ , h) is the kernel function. e kernel function may be represented by the cubic splines, which satisfy both accuracy and smoothing requirements, as well as minimal control of endpoints and end derivatives. In this study, we choose the popular cubic spline function [29] as the interpolation polynomial of the kernel function.

Advances in Civil Engineering
where z � (|x − x ′ |/h), and |x − x ′ | is the particle distance. D stands for the dimension of the question (1, 2, or 3), and C equals (2/3), (10/7π), and (1/π), respectively, when D � 1, 2, and 3. e governing equations to describe the particle motion and mass conservation of landslides are read: where ρ, v α , t, σ αβ , x β , and f α are the density, the velocity, the time, the stress tensor, the position, and the body force, respectively, and α and β indicate the Cartesian coordinate directions. e discretized governing equations in SPH can be written as follows: where a is the acceleration caused by body force. In order to eliminate unphysical shocks and oscillations and stabilize the numerical calculation, the artificial viscosity π ij was introduced into the momentum equation: and it can be written as [28] In equation (6), α π and β π are constants, which are all typically set around 1.0 [51-52]. c ij � ((c i + c j )/2) is the average speed of sound for geomaterials, andv ij represents the particle velocity vector.ρ ij � ((ρ i + ρ j )/2) is the average density at points i and j.
where φ � 0.1h ij is inserted to prevent numerical divergences when two particles are approaching each other.
To alleviate the tensile instability during SPH computation, the artificial stress method [30] was adopted. en, the momentum equation can be written in the SPH formula: where F n ij R αβ ij is the artificial stress term applied to reduce the tensile instabilities.
For the general elastic-plastic behavior, the waste dump landslides can be modeled according to the Mohr-Coulomb criterion, which can be written in terms of three stress invariants [31]: e flow potential g [32] for the Mohr-Coulomb yield surface can be written as follows: In the equations above, J 2 is the second invariant of deviatoric stress, and J 3 is the third invariant of deviatoric stress. ϕ is the friction angle, c is the cohesion, p is the equivalent pressure stress, and q is the Mises equivalent stress. ε is the meridional eccentricity. c 0 is the initial cohesion yield stress, and ψ is the dilation angle measured in the p − R mw q plane at the high confining pressure. e is referred to as the deviatoric eccentricity. R mc and R mw are related to ϕ and ψ, respectively.
When dealing with large deformation problems, the Jaumann rate of the Cauchy stress that is frame-indifferent is used here and can be written as where _ ω αβ is the rate-of-rotation tensor. en, the final form of the stress-strain relationship for the model in SPH is

Advances in Civil Engineering
where m, r and t denote Einstein's notation. G is the shear modulus, and K is the bulk modulus. _ λ is the rate of the plastic multiplier.
Finally, the SPH formula of _ ω αβ and _ ε αβ can be written as, respectively, In the SPH model, the boundary conditions are defined as described in Boundary Conditions in Abaqus/Explicit. In this work, the interaction between SPH particles (the waste dump soils) and the boundary (the bedrock and the rigid wall) is based on node-to-surface contact. In Abaqus, the improved boundary behavior is based on the ghost particle method [33], and the tangential-type of SPH boundary behavior is defined. e SPH SURFACE BEHAVIOR option in conjunction with the SURFACE INTERACTION option is used to define boundary surface interaction properties between SPH particles and Lagrangian surfaces.

A Benchmark of the SPH Model
Axisymmetric collapse tests on dry grains of sand, salt, couscous, rice, and sugar were conducted by Lube et al. [34]. In these experiments, the granular materials were first placed in a cylindrical container and then moved away from the container quickly, allowing the granular materials to collapse freely under the action of gravity. In the literature [35][36][37][38], the granular material flows have been simulated by the SPH method.
In order to verify the accuracy and effectiveness of the SPH method, two groups of dry sand collapse simulation tests were carried out and compared with reference [38]. e numerical computational parameters are shown in Table 1.
In the present benchmark, the initial halfwidth, the initial height, and the initial column height to halfwidth ratio are d 0 � 0.025 m, H 0 � 0.1 m, and a � H 0 /d 0 � 4, respectively. Figure 2 shows the initial and final schematic configurations of the sand. To examine the sensitivity of the number of particles in the SPH method, the test is divided into two groups, in which the corresponding numbers are 6174 and 133570, respectively. In the two simulations, the computational parameters are completely identified. e simulating results of dry sand collapse are shown in Figures 3 and 4, respectively.
For the axisymmetric experiments, the following equations [38] were proposed to estimate the final runout radius and height of the sand column: In the equations, when a � 4, it yields that (d ∞ /d 0 ) � 4.229, and (h ∞ /d 0 ) � 1.086. Table 2 summarizes the basic information for each benchmark simulation.
It can be seen from the table that the final runout radius to halfwidth ratio and the final height to halfwidth ratio obtained by the numerical simulation of two groups of different particle numbers have no significant difference, and their errors are within the reasonable range shown in the literature [38]. e collapse tests simulated by the SPH method show that the results are basically consistent with those in the literature [38] and are also in good agreement with other similar studies [34][35][36][37]. ese show that the SPH method has high accuracy and effectiveness, and it may be used to conduct the numerical simulation of the landslides in the full high waste dump.

SPH Simulation for the Full High
Waste Dump

Background of the Landslides.
e Nanfen open-pit mine is located in Nanfen District, Benxi City, Liaoning Province,  Figure 2: Axisymmetric collapse of a cylindrical sand body: the yellow region denotes initial sand column, and the dashed curve denotes final configuration (drawn according to reference [38]).
which is the largest black metallurgical mine in China. e full high waste dump (Dump 2, as shown in Figure 5) is perched on the hanging wall of the stope with a top and bottom elevation of about 600 m and 300-320 m, respectively. e vertical drop of waste dump is about 300 m, and the slope is 33°-35°. e waste soils are discharged with a transportation belt and dumping machine. e U-shaped valley where the dump is located is along the approximate east-west direction. e waste dump site belongs to a low mountain valley landform, with exposed gravelly soil and some weeds at the bottom of the valley. e main strata in this area include the gray-white, medium-fine quartz sandstone with thin black shale in Diaoyutai formation, and the egg-blue marl in the South Group, of which the marl is the most widely distributed. e strike of marl is 75°-160°. e dip angle is 15°-30°. e strike of quartz sandstone is 145°. e dip angle is 16°, and the Quaternary Diluvium covers it, with a thickness of 0-6.3 m. e mining area is located in the monsoon climate area of the north temperate zone. e annual average, lowest, extremely low, and high temperature is 7.

Simulation of the Landslides.
Nanfen open-pit mine has been constructed and operated for a long time. According to the existing topographic map of the mining area, section 1-1 passes through the original terrain valley. Section 2-2 has the largest slope angle, while section 3-3 crosses the original terrain ridge. In the present study, section 2-2 is selected for the established geometry model because it has the largest slope angle (see in Figure 8). e detailed geological generalization model of section 2-2 is shown in Figure 9.
Waste dump soils are a mixture of gravel, block stone, ore debris, gravel, fine sand, and clay. e lithology of the dump soils is mainly composed of mixed granite gneiss, quartz-feldspar schist, green mudstone, and amphibole schist. e soils have low viscosity, high particle strength, and good water stability. e content of each grade of grains in the dump soils is shown in Table 3. e soils are well graded. In the field of geotechnical engineering, the Mohr-Coulomb model is widely used because of its few parameters, it is easy to obtain, it has simple concept, and it can reflect the friction characteristics of geotechnical materials. In this study, the mechanical properties of dump soils are in good agreement with the Mohr-Coulomb model. erefore, the SPH method combined with the Mohr-Coulomb model is used to simulate the waste dump landslides. e simulation model size is set as follows (as seen in Figure 10). e slope top has the length of 300 m and the width of 300 m. e height difference from the slope toe to the slope top is 300 m, and the slope angle of the dump is about 33°. e horizontal distance from the toe of the slope to the left boundary of the model is 450 m. e bedrock thickness at the left boundary is 30 m. Combined with the field investigation, and laboratory tests conducted before and after landsliding, the mechanical properties are shown in Table 4.
In Figure 10, the boundary conditions are as follows. Two rigid walls are set on the left and right of the waste dump to limit the displacement in the Z direction, and a rigid wall is set on the rear of the waste dump to limit the displacement in the X direction. All of the displacements on the bottom boundary of the bedrock are fixed, whereas the top surface is without any constraint. e gravity loads are applied to the entire model. During the simulation, the general contact is activated between the dump soils and the ground, as well as the dump soils and the rigid walls. e contact friction coefficient is 0.6. e total simulation time is 100 seconds. e dump soils are transformed into particles when the calculation starts. e total number of particles is 39240. e initial particle distance is 10 m. e scheme of explicit time integration is used to solve the discrete equations of concentrated mass matrix [39].
To observe the sliding process of landslides and analyze the movement characteristics of the dump soils after failure, the displacement contours (U1, the displacement component along the direction of X-axis) and movement shapes of the full high waste dump at t � 0 s, 8 s, 17 s, 30 s, 40 s, and 100 s are given in Figure 11. e landslides reach their maximum displacement of 396 m from the slope toe to the front of the landslides at about 33.5 s after failure, which is close to the actual distance of 420 m by field investigation. Figure 12 shows the sliding distance (U magnitude, total displacement) and movement shapes of the waste dump landslides at t � 0 s, 8 s, 17 s, 30 s, 40 s, and 100 s.  Figure 13. It is shown that the landslides firstly experience a process of acceleration and then decelerate gradually. At about t � 17 s, the velocity reaches its maximum and then gradually decreases to zero. In the sliding process, the bodies of the landslides gradually advance forward and finally cease slowly.
To clearly understand the dynamic process of the dump landslides, we take five points at the different positions of the slope with equal spacing to track their displacements and velocity information. At the same time, the particle displacements with different heights in the vertical axis are presented. e particle positions on the dump slope are shown in Figure 14.
e particle displacements at the different positions on the slope are obviously different, but the displacement (U 1)      (Figures 15(a) and 15(b)). However, the displacements of the two vertical particles show that the higher the elevation of the particle is, the larger the displacement is (Figures 15(c) and 15(d)). Particles at different positions on the slope undergo a process of acceleration and then deceleration (Figure 15(e)). Figure 15(f) shows the energy information with respect to time. e sliding process and profile evolution of the dump landslides simulated by the SPH model at the different times (t � 8 s, 17 s, 30 s, and 100 s) are presented in Figure 16. e final profile of the landslides is basically consistent with the Table 4: e material properties of the SPH model [9].

Materials
Parameters       actual landslides, only except for the failure range of the SPH model at the slope top, which is slightly larger than the actual one. Besides, the landslides shape obtained by the SPH simulation is flatter than that of the field investigation. e distance from the original slope toe to the front of the actual landslides is about 396 m, which is basically consistent with the field measured value of 420 m.

Discussion.
is study shows that the SPH simulation can well reproduce the dynamic process of the full high waste dump landslides, and the numerical predictions are in good agreement with the sliding distance, velocity, and movement shapes on the actual landslides.
It can be seen from Figures 10 to 12 that, due to the limitation of boundary conditions, the displacement and velocity of particles are relatively large near the center of the model and small on both sides. It is difficult for the finite element method to describe the large deformation of materials. Sometimes it will encounter problems such as mesh distortion and nonconvergence. Conversely, the SPH simulation can fully present the entire dynamic process of the dump landslides, and its computational efficiency is very high in dealing with large deformation problems.
In addition, the SPH simulation can track the movement characteristics of particles at any position to understand the overall dynamic process. For example, in Figure 14, we select some representative particles to analyze the movement characteristics of the dump landslides. rough the analysis of Figures 15(a) and 15(b), we find that the displacements of particles near the middle and upper part of the slope are the largest, while, followed by the middle and lower part, the particle displacements near the foot of the slope are the smallest. e reason for these is that the lower particles have the least potential energy. Besides, in the process of sliding, they are first limited by the friction resistance of the ground and could be buried by the particles behind. Furthermore, the particles in the middle and upper part have the larger potential energy. Also, they are pushed by the upper particles in the downward sliding process, which makes the sliding distance of particles in this area possibly the largest. In Figure 15(d), the particles (on the vertical axis of the slope shoulder) with the elevation of 400-500 m have a small sliding distance of 0-100 m, which indicates that the deeper the particle is, or the closer the particle is to the underlying layer, the smaller the sliding distance is. is is consistent with the actual situation. e velocity-time curves of particles at different positions of the slope are similar (Figure 15(e)), and they all experience the process of acceleration and deceleration. However, it is obvious from the curve that the particles at the lower and midlower of the slope experienced a relatively short time of motion. is can also be seen in Figures 15(a) and 15(b). Certainly, there exists some local velocity fluctuation, which can be explained by the interaction of the particles with each other during the sliding process.
On the other hand, in the process of landslides movement, the shape of the sedimentary area is affected by the topography and the change and dissipation of energy. In Figure 16, the SPH simulated landslides profile is close to the actual landslides, but only the deposition height at the original toe is slightly different (Figure 16). e SPH simulated landslides profile is a smooth curve, whereas the actual profile is slightly uneven in some places, which is related to the terrain. e literature [40] has carried out a sensitivity analysis to quantify the influence of all uncertain input parameters. Inspired by it, we carried out a sensitivity analysis of different particle numbers to the results of numerical simulation. Under the same conditions, we set the particle number of the waste dump to 11700 and compare the simulation results with the former with the particle number of 39240. e displacement and velocity contours of the full high waste dump landslides at different times (t � 0 s, 8 s, 17 s, 30 s, 40 s, and 100 s) are shown in Figures 17 and 18, respectively. Compared with Figures 12 and 13, it can be seen that the displacement and velocity contours of the waste dump with different particle numbers are basically the same. e simulated  Figure 19 shows the displacement-time and velocity-time curves of particles at the 5 positions on the slope. Compared with Figure 15, the displacement and velocity distributions of particles at the same position are basically the same. e two groups of dump landslides simulated by the SPH model with different particle numbers are relatively close in shape, sliding distance, and speed. rough the sensitivity analysis of different particle numbers, it can be seen that the number of particles has little effect on the numerical results. Besides, the waste dump soils are composed of gravel soil with a maximum diameter of no more than 300 mm formed after secondary crushing, and the gradation is good, which is particularly suitable for SPH modeling by the meshless particle method. From the numerical simulation of the landslide profile, the slope is a relatively smooth curve, which is relatively close to the actual landslide. e material source of the dump landslides is different from other ordinary slopes. ere is no    obvious sliding surface in the waste dump landslides, and it appears as a flow-type landslides. erefore, the landslides profile is relatively smooth.
Based on the numerical results (as shown in Figure 16), the risk range can be delineated, safety measures can be formulated, and protective fortifications can be constructed, such as interception dams. ese measures can reduce the risk of waste dumps and provide guiding recommendations for the safe operation of mines.

Conclusions
is work first introduced the SPH theory and the elasticplastic Mohr-Coulomb model. en, a benchmark test was carried out to verify the effectiveness and accuracy of the SPH method. After that, the background of the landslides in the full high waste dump was given, and the 3D numerical model was established for the most dangerous sections, and the whole process of landslide calculation was analyzed. According to the calculation results and comparative analysis with the on-site investigation, a comprehensive analysis was carried out in terms of landslides profile, sliding distance, and velocity, as well as the sensitivity analysis of different particle numbers. Specifically, the following conclusions can be drawn: (1) is work has reproduced the entire dynamics process of the full high waste dump landslides by using the SPH method combined with the Mohr-Coulomb model. e simulated results are basically consistent with the actual ones. is is an innovative application of the SPH method in the large deformation problem of the open-pit dump slope.
(2) e sliding distance of particles in the middle and midupper of the dump slope is relatively large, and their maximum speed is also slightly larger. e sliding distance of particles located in the lower and midlower of the slope is relatively small, and their sliding distances and sliding time are relatively short.
(3) e monitoring particles in the vertical direction of both the shoulder of the waste dump and the middle of the slope were traced, and their sliding distances decrease with the increase in depth. (4) Compared with other numerical methods, the SPH method has the advantages of high computational efficiency and low computer memory usage. e SPH model is able to reproduce a postfailure configuration aiming at the sort of landslide problem. Meanwhile, with field investigation and satellite data, the 3D landslide model can be accurately and easily built up. It provides a new way that is interesting and useful to build the full high waste dump landslides model. (5) e sensitivity analysis of different particle numbers shows that different particle numbers have little effect on the numerical simulation results. e results of the two sets of numerical models (the number of particles is 39240 and 11700, respectively) are basically the same.  14 Advances in Civil Engineering simulation will be combined with the in situ monitoring information to optimize the support and reinforcement design scheme of waste dump slope and to improve the ability of disaster prevention and control in the waste dump.

Data Availability
e data used to support the findings of this study are included in the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.