Advances in Civil Engineering
Volume 2020 (2020), Article ID 8859344, 13 pages
https://doi.org/10.1155/2020/8859344

Research Article

Forecasting Disastrous Characteristics of Highway Landslides Using the Material Point Method: A Surcharge-Induced Perspective

Jian Wang, Wei Zhang, Chao-yang Xie, Qian Xiao, Yu-chen Zhu, and Bo Pan

School of Earth Sciences and Engineering, Nanjing University, Nanjing 210046, China

Correspondence should be addressed to Wei Zhang; wzhang@nju.edu.cn

Received 7 May 2020; Revised 14 July 2020; Accepted 21 July 2020; Published 7 August 2020

Academic Editor: Qiang Tang

Copyright © 2020 Jian Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Precisely forecasting the disastrous characteristics of landslides can address the optimized preliminary design of the highway landslide stabilization. In this paper, taking the surcharge-induced landslides as an example, we developed the models of the material point method (MPM), a novel meshless numerical method, to calculate the large deformation of the highway landslides, with the affordable computing power. After convincing verification of the MPM model to simulate a surcharge-induced soil landslide, 32 typical postfailure scenarios were analyzed to obtain the highway landslide run-out processes and disastrous characteristics, such as the sliding distance and speed of the leading edge and sliding body morphology. Moreover, linear regression equations of the maximal sliding distance and speed were deduced after verification, to forecast the reasonable avoiding distance. The maximal sliding distance and speed were found to be negatively linearly correlated with the internal frictional angle and cohesion of the soil, and positively linearly correlated with the surcharge and the slope angle, respectively. Optimized preliminary design, of the highways in the mountainous and hilly areas, can be performed based on those insights.

1. Introduction

Recent decades have seen a boom in transportation infrastructures in developing countries like China [1–4]. Given the highway planning, the engineering geological conditions are among the decisive factors when selecting the routes [5]. Compared with the traditional route selection procedure, geologic route selection, performing the engineering exploration before the preliminary design, can avoid major adverse geological areas, reduce the hazardous impact on highways, and improve the safety and economy of highway construction [6].

Landslides constitute one of the most common geological disasters when building highways in mountainous and hilly areas. Once the potential landslide hazard is not fully appreciated, the disastrous geotechnical actions may cause huge economic losses and casualties, both in the construction and operation stages. Given that the disastrous characteristics, such as the sliding distance, sliding speed, and sliding body morphology, can be precisely predicted, the distance of avoiding and winding can be reasonably arranged and the design of highway landslide stabilization can be optimized [7].

However, traditional slope stability analysis methods, such as the limit equilibrium method [8] and strength reduction FEM [9], cannot represent the disastrous characteristics of highway landslides, due to their computational incapability of the large deformation. Various numerical methods have been used to simulate large deformation problems. Zhou et al. [10] used the two-dimensional particle flow code model to study the failure mechanism of Yangjiagou landslide, analyzing the typical dynamic process of the formation of landslides and barrier lakes. Wu et al. [11] simulated the kinematic behavior of sliding rock blocks of the Tsaoling landslide induced by the earthquake, based on seismic discontinuous deformation analysis (DDA). Zhu et al. [12] reproduced the flow process of the landslide at Hongao landfill in China and estimated the arriving area using a full three-dimensional smoothed particle hydrodynamics (SPH) method. Although the above methods can be used for large deformation calculation, they cost such huge computing power that makes them not adequate to be used for the large-scale practical projects [13].

The material point method (MPM) developed by Sulsky et al. [14] is one of the promising meshless methods for overcoming such limitation. The MPM uses a double Lagrangian–Eulerian discretization, in which the material domain is divided into a series of particles. The momentum equations are solved on the background grid so that the grid distortion is avoided. The MPM is capable of solving such large deformation problems as landslides and debris flows [15–24]. Compared with other meshless methods, the MPM has higher computational efficiency and better computational stability [20, 21]. Moreover, developed calculation platforms, such as Anura3D [25] and MPM3D [26], have afforded adequate computing powers for practitioners.

In this paper, taking the surcharge-induced landslides as an example, we developed the MPM model to calculate the large deformation of the highway landslides. Typical postfailure scenarios were analyzed to obtain the landslide run-out process and disastrous characteristics, such as the sliding distance, sliding speed, and sliding body morphology. Moreover, linear regression equations, of the maximal sliding distance and speed, were deduced to forecast the disastrous characteristics of similar highway landslides.

2. Methodology

2.1. Basic Principle of MPM

As shown in Figure 1, the MPM is described by a Lagrange particle system and Euler grid [14]. The discrete particle system carries all of the material information of the medium. Computational grids are used only for solving momentum equations and spatial derivatives.

Figure 1: Schematic diagram of MPM.

Without considering heat exchange, the governing equations of the point-of-matter method include the virtual work equation, momentum equation, geometric equation, constitutive equation, and boundary condition equation, according to the principal parts of literature [14] and its derivation [27].

2.2. Discretization of the Continuum

Using a background grid with linear functions to discretize the continuum, the density of the continuum can be expressed as

where is the total number of particles, is the mass of particle , is the Dirac Delta function, and is the coordinate of particle . According to literature [20, 21], the weak forms of the coupling equations, with boundary conditions, yield

where is the displacement, is the acceleration, is the density at the current moment, is the cauchy stress tensor, is the volume force, is the pointing surface force, A is the area, V is the volume, and i and j are the space coordinates, respectively.

By introducing (1) into (2), the virtual work equation (2) is simplified as the summation form:

where the physical quantity with subscript p means that it is related to the material point p, and h is the imaginary boundary layer thickness.

The layer thickness of the imaginary boundary, h, is introduced to transform the last part of the virtual work equation to the body integral. Since no relative movement exists between the particle and the background mesh, mapping information between the particle and the nodes of the background grid is performed by the finite element shape function, based on the background mesh node [28], which yields

where , and are the coordinates, displacement, displacement derivatives, and virtual displacement of particle , respectively. The variables with subscript represent the variables of the grid node.

Given that the node virtual displacement is zero on the essential boundary and arbitrary at all the other nodes, combining (5) and (7) with (3), the motion equation of the background grid node yields

where

where is the momentum in direction of node number , and is the velocity in direction of node .

where is the mass matrix of background mesh, and and are the values of the shape function of nodes I and J at the particle point p, respectively.

where is the internal force of node, and

where is the external force of the node.

Introducing the lumped mass matrix yields

where is the lumped mass matrix.

We simplify (9) as

The motion equation of the background grid node yields

3. Verification

Surcharge-induced landslides are common disastrous scenarios occurring along the highways in mountainous and hilly areas because buildings on hills and mountains exert loads [16]. To verify the MPM algorithms in simulating the surcharge-induced landslide, the numerical example of Manna et al. is used. Manna et al. [29] performed a model-box test using Yamuna sand to construct the soil slopes with angles of 30°, 45°, and 60°, respectively. As shown in Figure 2, the dimensions of the model box are 0.75 m, 0.3 m, and 0.4 m, in length, width, and height, respectively. A gradually increasing load was applied on the slope crest to observe its instability process. The slope, with an angle of 45°, failed under the load of 4.28 kN, with a settlement of 11.40 mm; the surface bulged and moved about 50 mm forward at the toe.

(a)

(b)

Figure 2: Failure morphology of the test slope [29]: (a) top view, (b) profile.

The experiment used Yamuna sand. The soil density, internal frictional angle, and elastic modulus are 1568 kg/m³, 37°, and 5 MPa, respectively. Note that the cohesion is negligible.

The instability process of the slope model, with an angle of 45°, is simulated using the MPM. The soil parameters of the simulation are consistent with those used in the experiment. As shown in Figure 3, an MPM slope model, with a height of 0.3 m, a length of 0.46 m, a slope angle of 45°, and a foundation thickness of 0.1 m, was established. The length of the foundation is 0.75 m, 0.29 m longer than the bottom of the landslide. The width of this MPM model is 0.3 m. At the left and right boundaries, the grid nodes are fixed along the x-direction. At the bottom boundary, the grid nodes are fixed along the x-, y-, and z-directions. At the boundaries of z = 0 m and z = 0.03 m, the grid nodes are fixed along the z-direction. A uniform load of 93 kPa is initially applied on top of the slope. The Mohr–Coulomb yield criterion is used to calculate the plastic strain. The computational grid size is set as 0.03 m; the total computation time is 5 s. The plastic strains of the slope on typical calculation moments are shown in Figure 4.

Figure 3: Schematic of a slope subjected to surcharge.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Plastic strain evolution process: (a) time = 0 s, (b) time = 1 s, (c) time = 1.5 s, (d) time = 2 s, (e) time = 4 s, (f) time = 5 s.

As seen from Figure 4, under the surcharge, compression deformation occurs at the top of the slope; the soil extrudes outwards along the slope surface; the damage occurs at the slope toe, forming a plastic strain zone. Note that, in Figure 4(f), the red dotted line sketches the initial form and position of the slope; the calculated crest settlement and sliding distance are 12.9 and 56.8 mm, respectively, both having errors less than 15%, suggesting that the MPM model is reasonable and convincing to simulate the run-out of surcharge-induced soil landslides.

4. Postfailure Simulation

4.1. Model Parameters

A typical highway landslide scenario is depicted in Figure 5. As shown, an ideal homogeneous soil slope, with a height of 20 m, a length of 67 m, and a slope angle of 35°, is at the side of the new highway to be built. An MPM model is developed to simulate the postfailure stage of the slope to forecast the disastrous characteristics for the optimized preliminary design. As seen from Figure 5, the length of the foundation is set as 130 m, 63 m longer than the bottom of the landslide, and the thickness as 10 m. The bottom of the model is restricted along all the three directions: x, y, and z; the left and right sides of the model are restricted along the x-direction; the top side and the surface of the slope are not restricted; and the rear and front of the model along the z-direction are restricted to simulate the plane strain state. Three measuring points, A, B, and C, are set on the shoulder, middle, and toe of the slope, respectively. The slope soil is silty clay, whose physical and mechanical parameters are listed in Table 1.

Figure 5: Slope MPM model (total material point number: 30588).

Table 1: Physical and mechanical parameters of the soil.

Density (kg/m³)
Modulus of elasticity (MPa)
Poisson’s ratio
Cohesion (kPa)
Internal friction angle (°)

2000
20
0.3
20
25

4.2. Initial Stress Field

We set up a load with a height of 5 m on top of the slope, ranging from the slope shoulder point A to 20 m inward. Note that the surcharge denotes the buildings on top of the slope [30]; various surcharges can be exerted by changing the density. After applying the gravity, the initial stress field of the slope is obtained, as shown in Figure 6, in which the coordinate axis denotes the slope dimensions and the legend exhibits the magnitude of the effective stress. As can be seen, the effective stress increases gradually, under gravity, from top to bottom of the slope.

Figure 6: Initial stress field.

4.3. Run-Out Simulation

An extra intensity, 4 times of the soil density, is applied to conduct the surcharge for destruction. Note the load on top of the slope is equivalent to 400 kPa. The Mohr–Coulomb yield criterion is used to analyze the run-out process of the landslide. Figure 7 exhibits the snapshots of the run-out process of the slope failure; the development of the rupture surface can be observed in the plastic strain zone. As seen, the plastic strain occurs on the 2 s computational time. At 3 s, the sliding zone occurs within the slope, and the load on top starts inclining. At 4 s, the sliding zone expands, and the integral sliding starts. Subsequently, from 6 s to 10 s, the plastic deformation zone extends further to the slope interior; the slope slides forward along the sliding zone. At 12 s, the landslide deformation ceases, and the plastic strain zone tends to be stable.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7: Evolution of displacement field and plastic strain zone: (a) time = 1 s, (b) time = 2 s, (c) time = 3 s, (d) time = 4 s, (e) time = 6 s, (f) time = 8 s, (g) time = 10 s, (h) time = 12 s.

Figure 8 depicts the horizontal speed curves of the 3 measuring points, A, B, and C, during the run-out process. As shown in the figure, the 3 measuring points accelerate immediately after the loading. At 4 s, the speed tends to be uniform and fluctuates in a small range. At 7 s, the speeds start decreasing, and the sliding stops after 9 s.

Figure 8: Evolution of sliding speed.

Among the 3 measuring points, the point C has the maximal speed, reaching 2.68 m/s, during the acceleration stage. Points B and A both have a lower maximal speed of about 2.31 m/s.

Figure 9 depicts the horizontal displacement curves of A, B, and C. As seen from the figure, point C has the longest sliding distance of 14.23 m; point A has the shortest one of 11.06 m. Therefore, the avoiding distance, of the newly built highway from the toe of this slope, should be greater than 14.23 m.

Figure 9: Evolution of sliding distance.

According to Figures 7-9, the run-out process can be divided into four stages: the first is the local creep stage. As seen from Figure 7(b), after applying the surcharge, the plastic strain occurs at the slope toe and gradually develops towards the interior. The leading edge of the slope begins to bulge, and the soil of the slope toe starts creeping forward as the shear stress exceeds the shear strength.

The second is the acceleration stage. As seen from Figures 7(c) and 7(d), the plastic strain zone penetrates from the slope toe to the interior of the slope, with the sliding zone forming. The slope shears out from the slope toe along the sliding zone, resulting in an integral sliding.

The third is the uniform sliding stage. As seen from Figures 7(e), the plastic strain zone continues to expand, with a high plastic strain zone occurring. Meanwhile, the sliding mass moves forward rapidly as a whole.

The last is the deceleration state. As seen from Figure 8, from 6 s to 7 s, the landslide varies from the uniform sliding stage to the deceleration one, with the sliding mass continuing moving forward and the maximal sliding distance at the leading edge reaching 14.23 m, as shown in Figure 7(h).

To sum up, the slope failure exhibits progressive failure characteristics [31]. At the initial loading stage, the slope experiences the creep deformation. As the plastic strain accumulates and gradually develops towards the slope surface, the slope slides along the arc-shaped sliding zone. The deformation and failure characteristics of the case are similar to those of the slumping slide, which exhibits compressive and shear failure [31], consistent with the experimental phenomena observed in literature [32].

5. Disastrous Characteristic Analysis

A landslide with a high speed has a great impact and wide damage scope. The disaster-inducing scope and damage capability of a landslide can be approximately characterized by the maximal sliding distance and speed, respectively. By analysis with different slope morphologies, surcharges, and physical and mechanical parameters, the insight of the impacts of various factors on landslide-inducing disastrous behaviors can be gained.

Figure 10 exhibits the MPM simulated results of the scenarios with different surcharges on top of the slope. Still, the soil cohesion is 20 kPa; the internal friction angle is 25°; and the slope angle is 35°.

(a)

(b)

Figure 10: Surcharge impact: (a) evolution of sliding speed, (b) evolution of sliding distance.

As shown in Figure 10(a), the maximal speed increases with the increase of surcharge. Compared with the maximal speed under the surcharge of 300 kPa, maximal under 400 kPa increases significantly, while the maximal speed varies a little from 400 to 600 kPa. In Figure 10(b), with the surcharge increase, the sliding distance increases gradually, while the amplitude decreases gradually.

Figure 11 exhibits the results of the scenarios with different soil cohesion. The surcharge remains 400 kPa; the internal friction angle is 25°; and the slope angle is 45°. As seen from the figure, with the increase of the soil cohesion, the maximal sliding speed and distance decrease. Compared with the surcharge of 15–20 kPa, when the cohesion varies from 20–25 kPa to 25–30 kPa, both sliding speed and distance decrease significantly.

(a)

(b)

Figure 11: Cohesion impact: (a) evolution of sliding speed; (b) evolution of sliding distance.

Figure 12 exhibits the results of the scenarios with different internal friction angles. The surcharge is 400 kPa; the cohesion is 20 kPa; and the slope angle is 35°. As shown in the figure, with the increase of the internal friction angle, the maximal sliding speed and distance decrease approximately linearly, suggesting that the impact of the internal friction angle on landslide speed and distance is uniform.

(a)

(b)

Figure 12: Internal friction angle impact: (a) evolution of sliding speed, (b) evolution of sliding distance.

In order to study the impact of different slope shapes on landslide movement, the slope angle is set as 30°, 35°, 40°, and 45°, respectively. The surcharge remains 400 kPa; the friction angle is 25°; and the cohesion is 20 kPa. The results are shown in Figure 13.

(a)

(b)

Figure 13: Slope angle impact: (a) evolution of sliding speed; (b) evolution of sliding distance.

As shown in Figure 13, with the increase of the slope angle, both maximal sliding speed and distance decrease. However, when the slope angle increases to a certain extent, the impact of the surcharge on the sliding speed and distance reduces, the gradient of the deceleration phase decreases gradually, and the deceleration phase extends.

In general, with the increase of the surcharge and slope angle and the decrease of the soil strength parameters, the sliding speed and distance increase significantly. Note that the characteristics of landslide movement, such as sliding distance, speed, and acceleration, are significant for evaluating the scope and degree of a landslide disaster. 32 slope MPM models are established with variations of surcharges, internal friction angles, cohesion, and slope angles. Taking the sliding distance and speed of measuring point C into account, parametric sensitivity analysis is performed and the results are listed in Table 2.

Table 2: Parametric sensitivity analysis.

Scenario number
Surcharge (kPa)
Cohesion (kPa)
Internal friction angle (°)
Slope angle (°)
Sliding distance (m)
Sliding speed (m/s)

1
300
20
25
35
8.492
1.5626

2
350
20
25
35
11.415
1.9783

3
400
20
25
35
14.227
2.6862

4
450
20
25
35
15.569
2.7328

5
500
20
25
35
16.911
2.9698

6
550
20
25
35
18.318
3.3343

7
600
20
25
35
19.185
3.5408

8
400
15
25
35
15.681
2.832

9
400
17
25
35
15.174
2.7772

10
400
19
25
35
14.251
2.5446

11
400
20
25
35
14.227
2.6862

12
400
21
25
35
13.027
2.3332

13
400
23
25
35
11.954
2.0658

14
400
25
25
35
10.759
1.8943

15
400
27
25
35
9.4391
1.716

16
400
29
25
35
8.2382
1.4373

17
400
30
25
35
7.6099
1.2967

18
400
20
20
35
20.886
3.8407

19
400
20
22
35
18.159
3.3991

20
400
20
24
35
15.475
2.9001

21
400
20
25
35
14.227
2.6862

22
400
20
26
35
12.912
2.4793

23
400
20
28
35
10.175
1.9355

24
400
20
30
35
7.61
1.3646

25
400
20
32
35
5.46
0.83115

26
400
20
34
35
3.5008
0.41482

27
400
20
35
35
2.7073
0.25283

28
400
20
25
25
3.9192
0.8591

29
400
20
25
30
8.1405
1.4788

30
400
20
25
35
14.227
2.6862

31
400
20
25
40
18.994
3.7907

32
400
20
25
45
20.833
4.2865

In terms of the 32 scenarios in Table 2, the correlation curves of surcharge, cohesion, internal friction angle, and slope angle to the maximal sliding distance and speed are shown in Figure 14(a)-14(d), respectively. As noted from those figures, the maximal sliding distance and speed vary approximately linearly with the variation of each parameter;

(a)

(b)

(c)

(d)

Figure 14: (a) Correlation curves of MSD and MSS to surcharge. (b) Correlation curves of MSD and MSS to cohesion. (c) Correlation curves of MSD and MSS to internal friction angle. (d) Correlation curves of MSD and MSS to slope angle. (MSD: maximal sliding distance; MSS: maximal sliding speed).

The linear regression equations of the maximal sliding distance and speed, of measuring point C, were deduced, respectively, in terms of surcharge, cohesion and internal friction angle of the soil, and slope angle, as follows:

where p, c, , and are the surcharge, cohesion, internal friction angle, and slope angle of the soil, respectively. s, denote the maximal sliding distance and speed of measuring point C, respectively.

Note that only the first 31 scenarios are used for regression analysis. Equations (16) and (17) are used to forecast the maximal sliding distance and speed of scenario 32, respectively. The predicted results are compared with those computed by MPM. Note that the errors were only 8% and 2%, for maximal sliding distance and speed, respectively, indicating a convincing efficiency of the proposed forecasting models. Furthermore, both equations can be used for optimized preliminary design of the highways in mountainous and hilly areas. Note that using nonlinear regression functions, such as power or logarithmic function, cannot obtain satisfying forecasting results.

6. Discussion

According to Figure 14, the cohesion and the internal friction angle are negatively linearly correlated with sliding distance and speed. The slope can be enhanced by improving the soil properties, such as grouting, electroosmotic drainage, electrochemical reinforcement, and vegetation recovery. Note that the surcharge and slope angle are positively linearly correlated with the sliding distance and speed. Also, the slope can be enhanced by lowering the slope height, decreasing the slope angle, or, reducing the surcharge. If necessary, engineering prevention measures such as retaining walls, anti-slide pile, or anchorage could be deployed at the slope toe.

Note that DEM and DDA can be used to compute the large deformation of landslides; however, they are not capable of addressing fluid-solid coupling issues, which are among the main causes of landslides. Although SPH can simulate the run-out of landslides, its approximate function does not always meet the interpolation condition, which does not guarantee the single-value characteristic of the velocity field, because the material interface may be penetrated or mixed, so that the definition of displacement boundary is vague. As a meshless numerical method, MPM can introduce various soil constitutive models and address fluid-solid coupling. Compared with other numerical methods, MPM holds the potential of gaining higher computational efficiency and numerical stability [33, 34].

7. Conclusions

The typical postfailure scenarios of the surcharge-induced soil landslides are simulated by establishing a series of MPM models; the disastrous characteristics are analyzed and linear regression equations are deduced for forecasting. The following conclusions can be drawn:

(1)MPM is a meshless numerical method, which can effectively simulate large deformation of landslides with affordable computing power. It can reproduce the entire run-out process and obtain disastrous characteristics of landslides, such as sliding distance, sliding speed, and sliding body morphology.(2)The surcharge-induced landslide is a thrust-type landslide; its mechanism lies in the compression-shear failure, which is mainly based on the overall sliding along the arc-shaped sliding zone. The whole run-out process can be divided into four stages: local creep, accelerated sliding, uniform sliding, and decelerated sliding.(3)Optimized preliminary design, concerning avoiding distance and other engineering measures, of the highways in the mountainous and hilly areas can be performed based on those insights deduced from the MPM simulation and subsequent linear regression analysis.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This study was financially supported by the National Key R & D Program of China (Grant no. 2018YFC1505100) and the Natural Science Foundation of Jiangsu Province (Grant no. BK20181182).

References

J. He, S. Wang, Y. Liu et al., “The development of composite bridges with corrugated steel webs in China,” Proceedings of the Institution of Civil Engineers-Bridge Engineering, vol. 172, pp. 1–44, 2019.

N. Hu, G.-L. Dai, B. Yan, and K. Liu, “Recent development of design and construction of medium and long span high-speed railway bridges in China,” Engineering Structures, vol. 74, pp. 233–241, 2014.

C. F. Lu, J. F. Liu, Y. H. Liu, and Y. M. Liu, “Intelligent construction technology of railway engineering in China,” Frontiers of Engineering Management, vol. 6, no. 2, pp. 503–516, 2019.

H. Han, J. Guo, J. Zhang, and Y. Sun, “Technical advances of temporary facilities for the failure prevention of large-span cantilever casting construction of mountainous concrete box-type Arch bridges,” Advances in Civil Engineering, vol. 2020, Article ID 6412613, 20 pages, 2020.

K. Makanae and N. Dawood, “Development and evaluation of a tangible terrain representation system for highway route planning,” Computer-Aided Civil and Infrastructure Engineering, vol. 24, no. 3, pp. 225–235, 2009.

C. C. Osadebe, B. M. Salami, S. E. Obrike, and M. O. Ajala, “The Kano-Kazaure highway, north central Nigeria: the significance of the engineering geology in construction,” Bulletin of Engineering Geology and the Environment, vol. 70, no. 1, pp. 33–40, 2011.

Ministry of Transport of the People’s Republic of China, JTG∕T 3334-2018 “Code for Design of Highway Landslide Prevention, China Communications Press, Beijing, China, 2018.

X. An, N. Li, P. Zhang, and W. Sun, “Error source analysis and precision assessment of limit equilibrium methods for rock slopes,” Advances in Civil Engineering, vol. 2018, Article ID 3280734, 13 pages, 2018.

M. C. M. Nasvi and S. Krishnya, “Stability analysis of colombo-katunayake expressway (CKE) using finite element and limit equilibrium methods,” Indian Geotechnical Journal, vol. 49, no. 6, pp. 620–634, 2019.

J.-w. Zhou, P. Cui, H. Fang et al., “Dynamic process analysis for the formation of Yangjiagou landslide-dammed lake triggered by the Wenchuan earthquake, China,” Landslides, vol. 10, no. 3, pp. 331–342, 2013.

J.-H. Wu and C.-H. Chen, “Application of DDA to simulate characteristics of the Tsaoling landslide,” Computers and Geotechnics, vol. 38, no. 5, pp. 741–750, 2011.

C. Zhu, Y. Huang, L.-t. Zhan et al., “SPH-based simulation of flow process of a landslide at Hongao landfill in China,” Natural Hazards, vol. 93, no. 3, pp. 1113–1126, 2018.

J. Meng, J. Huang, H. Lin, J. Laue, and K. Li, “A static discrete element method with discontinuous deformation analysis,” International Journal for Numerical Methods in Engineering, vol. 120, no. 7, pp. 918–935, 2019.

D. Sulsky, Z. Chen, H. L. Schreyer et al., “A particle method for history-dependent materials,” Computer Methods in Applied Mechanics and Engineering, vol. 118, no. 1-2, pp. 179–196, 1994.

X. P. Li, Q. W. Yan, S. X. Zhao et al., “Investigation of influence of baffles on landslide debris mobility by 3D material point method,” Landslides, vol. 17, no. 5, pp. 1–15, 2020.

S. Andersen and L. Andersen, “Modelling of landslides with the material-point method,” Computational Geosciences, vol. 14, no. 1, pp. 137–147, 2010.

L. Mirada and X. Francesc, “The material point method in slope stability analysis,” Ground Engineering, vol. 48, pp. 31–36, 2015.

T. Belytschko, Y. Krongauz, D. Organ et al., “Meshless methods: an overview and recent developments,” Computer Methods in Applied Mechanics and Engineering, vol. 139, no. 1-4, pp. 3–47, 1996.

S. F. Li and W. K. Liu, “Meshfree and particle methods and their applications,” Applied Mechanics Reviews, vol. 55, no. 1, pp. 1–34, 2020.

A. R. York, D. Sulsky, and H. L. Schreyer, “The material point method for simulation of thin membranes,” International Journal for Numerical Methods in Engineering, vol. 44, no. 10, pp. 1429–1456, 1999.

A. R. York, D. Sulsky, and H. L. Schreyer, “Fluid-membrane interaction based on the material point method,” International Journal for Numerical Methods in Engineering, vol. 48, no. 6, pp. 901–924, 2000.

Y. Dong, D. Wang, M. F. Randolph et al., “A GPU parallel computing strategy for the material point method,” Computers and Geotechnics, vol. 66, pp. 31–38, 2015.

S. Bandara and K. Soga, “Coupling of soil deformation and pore fluid flow using material point method,” Computers and Geotechnics, vol. 63, pp. 199–214, 2015.

B. Shi, Y. Zhang, and W. Zhang, “Run-out of the 2015 Shenzhen landslide using the material point method with the softening model,” Bulletin of Engineering Geology and the Environment, vol. 78, no. 2, pp. 1225–1236, 2019.

E. Conte, L. Pugliese, A. Troncone et al., “Post-failure stage simulation of a landslide using the material point method,” Engineering Geology, vol. 253, pp. 149–159, 2019.

X. Zhang, C. Zheng, and Y. Liu, The Material Point Method-A Continuum-Based Particle Method for Extreme Loading Cases, Academic Press, Cambridge, MA, USA, 2016.

W. Zhang, B. T. Shi, B. Shi, Y. Zhang, and F. D. Wu, “Run-out process simulation of soil landslide using material point method and its application,” Journal of Engineering Geology, vol. 25, no. 3, pp. 815–823, 2017, (in Chinese).

X. Zhang, Y. P. Lian, Y. Liu, and X. Zhou, Material Point Method, Tsinghua University Press, Beijing, China, 2013, (in Chinese).

B. Manna, S. Rawat, R. Zodinpuii et al., “Effect of surcharge load on stability of slopes - testing and analysis,” Electronic Journal of Geotechnical Engineering, vol. 19, pp. 3397–3410, 2014.

R. M. Mendes, M. R. M. de Andrade, C. A. Graminha, C. C. Prieto, F. F. de Ávila, and P. I. M. Camarinha, “Stability analysis on urban slopes: case study of an anthropogenic-induced landslide in são josé dos campos, Brazil,” Geotechnical and Geological Engineering, vol. 36, no. 1, pp. 599–610, 2018.

B. L. Wang, S. E. Peng, and H. K. Chen, “Mechanical deduction of formation mechanism for the thrust load caused landslide,” Journal of Geological Hazards and Environment Preservation, vol. 21, no. 2, pp. 74–77, 2010.

T. F. Hu, “Model test study of deformation mechanism of landslide induced by loading,” Journal of Shijiazhuang Tiedao University, vol. 29, no. 1, pp. 38–42, 2016, (in Chinese).

B. T. Shi, Y. Zhang, and W. Zhang, “Analysis of the entire failure process of the rotational slide using the material point method,” International Journal of Geomechanics, vol. 18, no. 8, 2018.

J. J. Shi, W. Zhang, B. Wang, C. Y. Li, and B. Pan, “Simulation of a submarine landslide using the coupled material point method,” Mathematical Problems in Engineering, vol. 2020, Article ID 4392581, 14 pages, 2020.