#### 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.

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 aswhere 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, yieldwhere 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 yieldswhere , 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 yieldswherewhere 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, andwhere is the external force of the node.

Introducing the lumped mass matrix yieldswhere 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)**

The experiment used Yamuna sand. The soil density, internal frictional angle, and elastic modulus are 1568 kg/m^{3}, 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

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.

##### 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.

##### 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 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.

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.

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)**

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 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)**

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)**

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.

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)**

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).