Laboratory and 3-D distinct element analysis of the failure mechanism of a slope under external surcharge

Landslide is a major disaster resulting in considerable loss of human lives and property damages in hilly terrain in Hong Kong, China and many other countries. The factor of safety and the critical slip surface for slope stabilization are the main considerations for slope stability analysis in the past, while the detailed post-failure conditions of the slopes have not been considered in sufficient detail. There is however increasing interest in the consequences after the initiation of failure that includes the development and propagation of the failure surfaces, the amount of failed mass and runoff and the affected region. To assess the development of slope failure in more detail and to consider the potential danger of slopes after failure has initiated, the slope stability problem under external surcharge is analyzed by the distinct element method (DEM) and a laboratory model test in the present research. A more refined study about the development of failure, microcosmic failure mechanisms and the post-failure mechanisms of slopes will be carried out. The numerical modeling method and the various findings from the present work can provide an alternate method of analysis of slope failure, which can give additional information not available from the classical methods of analysis.


Introduction
Slope stability problem is always a main concern in geotechnical engineering.Natural or cut slopes have greatly influenced humans' properties and lives, especially in a hilly terrain like Hong Kong and many other countries where slope failures can be fatal.Landslides have caused loss of life and a significant amount of property damage in Hong Kong, and for the 50 years after 1947, more than 470 people died as a result of failures associated with man-made cut slopes, fill slopes and retaining walls.The government has put forward a Landslip Preventive Measure with an annual expenditure of about HKD 870 million to upgrade sub-standard slopes and HKD 600 million for the maintenance of existing slopes.Even with the continuous upgrading works since 1977, about 300 incidents affecting man-made slopes, walls and natural hillsides are reported to the government every year.Besides Hong Kong, many major cities in China are facing similar situations, and slope stabilization becomes a very important activity, with different types of materials and construction techniques being adopted for soil nail slope stabilization.
Many slope stability analysis methods and tools have also been developed over the last 50 years.Each stability analysis method differs from the others in some of the basic assumptions, but most of these methods will give similar factors of safety/slip surfaces for normal cases, which are sufficient for normal engineering uses (Cheng and Lau, 2008).Besides the factor of safety, the failure mechanism and postfailure development are also important in some cases.For slope stability analysis, the limit equilibrium method (LEM) is widely used by engineers and researchers, and this is a traditional and well-established method.The finite element method (in the form of the strength reduction method, SRM) has become popular recently for the analysis of slopes, but there are also many practical problems to the SRM for some complicated problems, which have been discussed by Cheng et al. (2007).
The distinct element method (DEM) that is initiated by Cundall and Strack (1979) is now available as a commercial program (Itasca, 2008), and it is an alternate numerical method that is more suitable for large-scale ground movement with separation, so that it will be more convenient to  2009), Cheng et al. (2009Cheng et al. ( , 2010)), Jenck et al. (2009), andLu andMcDowell (2010).It should be noted that there are very limited applications of the DEM for slope stability analysis, especially for the progressive failure mechanism of slopes.Cheng and Lau (2008) have considered some slope failure by DEM for simple slopes under gravity, but it appears that there is a lack of detailed study about slope stability with DEM.DEM is also more suitable for qualitative instead of quantitative assessment of the stability of slopes at present, as there are difficulties in defining precisely the microscopic parameters and the contact model, and most of the DEM models cannot reflect completely the grain distribution, initial stress or drainage condition.
There have been many advances in the distinct element method over the last 20 years, and complex mechanical interactions of a discontinuous system can now be analyzed in three dimensions (though this is still very tedious and timeconsuming at present).There is only limited previous 3-D particle flow modeling on real problems (most works are devoted to the fundamental study of soil behavior) due to the various constraints in image capturing, model generation in three dimensions and a large volume of data for manipula- tion.In the present study, PFC3D is used for the analysis of slope laboratory test results where the interaction of the particles is based on the DEM as described by Cundall and Strack (1979).In this paper, the progressive failure mechanism of the model slope under local surcharge is determined both from the model test and numerical analysis for comparison, and the procedures to set up an appropriate 3-D distinct element model is discussed.In the present work, the authors have managed to carry out a tedious and very timeconsuming 3-D DEM analysis together with the laboratory test, aiming at quantitative slope failure analysis, which is an important work in distinct element analysis of the complete failure process of a slope.The present work has illustrated that the results from the numerical modeling and the laboratory test can also agree reasonably well, both qualitatively and quantitatively, which is not commonly found in the literature.

Physical modeling of slopes under external surcharge
A physical test was carried out in the laboratory to investigate the failure process of slopes under external local loading.A physical soil tank was built with the layout as shown in Fig. 1, and the tank was about 1.5 m deep × 1.84 m wide × 1.2 m high.The height of the model slope was 0.7 m, with a slope angle of 45 • .The sectional view of the soil slope model was illustrated in Fig. 2. Sandy soil with 5 % moisture content and particle size distribution as shown in Fig. 3 was compacted to form a slope model in seven layers using an electric vibrator.
The average bulk density of the soil model was determined as 1672 kg m −3 .A direct shear test was conducted to determine the shear strength parameters of the river sand used in the laboratory test, and the shear-normal stress relation was illustrated in Fig. 4, which gave a friction angle of 58.6 • to the compacted sand with a cohesion of 0.6 kPa.The parameters of the river sand used for the physical model were given in Table 1.
In this test, the surcharge loading was applied from the hydraulic jack, and the loading was transferred to the steel plate to simulate a local distributed load as shown in Fig. 1.The intermediate and final failure states of the model slope were illustrated in Fig. 6, and the results from the five transducers were given in Fig. 7. From Fig. 6a and e, we could see that cracks developed firstly around the steel loading plate and extended towards each corner side of the slope crest with an angle of 45 • , which was basically in accordance with the classical theory.During the test, flags with different colors were used to locate the time of appearance and the location of cracks during loading, which could also be seen in Fig. 6.From Fig. 6d, it was noticed that the failure surface was approximately a triangular zone at the top plus a log-  spiral zone similar to the Prandtl mechanism for a bearing capacity problem with an inclined surface below the top triangular zone (Cheng and Au, 2005).The inclination of the failure surface to the horizontal direction in Fig. 6d was also close to 74.3 • (45 • + ϕ/2), which was a good illustration of the general shear failure for the present test.The soil mass actually failed gradually instead of having a sudden failure in the loading process.As the steel plate kept going down, the inclined slope face was covered with vertical cracks and diagonal cracks that were directed towards the steel plate (see Fig. 6e).The soil mass at the slope surface was drawn down layer by layer to the slope toe by the action of the loading steel plate, and this was the typical face failure phenomenon for cohesionless soil.After 4 h of loading, the sand at the middle top of the slope was highly compressed, and a cavity was formed underneath the steel plate, which could be noticed after the removal of the steel plate, while the largest cracks were generated at the critical failure surface within the slope body for the global failure (Fig. 6d).The global soil mass was pushed down to the toe of the slope gradually, and the complete physical slope model collapsed eventually, which was shown in Fig. 6.
An ultimate load of 35 kN had been attained at a displacement of about 6 mm as shown in Fig. 7.For the slope surface, the corresponding displacement at the maximum pressure was about 2 and 1 mm at the top and bottom of the slope, respectively, as shown in Fig. 7, which was much lower than the corresponding displacement at the loading plate.Beyond the peak load, the applied load decreased with increasing jack displacement.From the displacements for the left and right LDTV at the upper and lower levels, it was clear that the displacements of the slope were basically symmetrical.In Fig. 7, after the maximum load had been achieved, the loading force decreased with increasing displacement until a displacement of about 13 mm for which the load maintained constant for about 9 mm.At the beginning of this constant load stage, the local triangular failure zone was fully developed, while the failure zones at the two ends of the plate were still not clearly formed.At the end of this constant load stage, the failure zones at the two ends became visually apparent.When the displacements were further increased, the applied load decreased further and the failure zone propagated towards the slope surface until the failure surface as shown in Fig. 6c and d is obtained, which is the typical global failure as considered in classical stability analysis.Finally, the residual load of the test was around 5 kN.
For this test, there were several interesting phenomena worth discussing.The failure profile and cracks first initiated beneath the footing as shown in Fig. 6e, which was a typical bearing capacity failure with a triangular failure zone.This could also be observed from the upper part of the failure profile as shown in Fig. 6c and d.As the load increased, the failure zone extended and propagated towards the toe of the slope and the final failure surface was shown in Fig. 6d.It was observed that the failure mechanism of the physical model test is a local triangular failure beneath the bearing plate, and the failure surface propagates towards the slope surface with a curved surface similar to a log-spiral curve until a failure mechanism is formed.This type of problem could be considered a bearing capacity problem as well as a slope stability problem, which were demonstrated to be equivalent by Cheng et al. (2013).

Three-dimensional numerical modeling of slopes under local surcharge
To assess the laboratory test results further, the authors have adopted the DEM in the numerical analysis, because the development of cracks, face failure and the final collapse are difficult to assess by finite element analysis.In DEM, there are several methods of model generation.For the present problem, which is relatively simple in geometry and layout, the desired porosity is obtained by the radius expansion method.By using numerical biaxial tests, the micromechanical properties of the assembled material in the numerical models are calibrated in order to match with the macroscopic response of the real material in the physical test.Numerical simulations to reproduce the stress-strain and normal-shear stress relations (Fig. 4) similar to that by Cheng et al. (2009Cheng et al. ( , 2010) ) are carried out under the same conditions as the physical experiments such as porosity, boundary conditions and loading.The micro-properties of the river sand as shown in Table 2 are determined by varying the micro-properties until the macro-properties obtained numerically match with the experimental results (angle of repose and stress-strain relation).The diameters of the particles in the DEM model are kept the same as those as given in Fig. 3.The frictional coefficient of sands is set to 1.638 (corresponding to a friction angle of 58.61 • in Table 1).The  bond strength is fixed at the value of 6 N as referred to by Cheng et al. (2003).The particle density of the sandy soil is 2650 kg m −3 , while the bulk density for the sandy soil is 1650 kg m −3 .For the 3-D DEM numerical simulation as shown in Fig. 8, the dimensions of the numerical model in the particle flow simulation are exactly the same as the physical model.A full-scale distinct element analysis is seldom carried out in the literature due to the great computer resource requirement, and a significant amount of computer time has been used for the analysis in the present work.In this study, two different loading patterns are modeled in order to assess the results of analysis using PFC3D: (1) applying the force to the raft footing, and (2) adding velocity to the loading wall, which are illustrated in Fig. 8.In general, the results from the two methods of simulations are similar, so only the results for simulation 1 are given in the present paper.

Raft footing loading pattern
In this loading pattern, the footing raft of balls is created by bonding the particle with larger contact bonds, fixing velocity constraints on the balls and modifying their stiffness and friction properties.The raft footing is added to the top of the model, and the footing load is applied in the vertical direction (x direction) with a magnitude of 3.5 × 10 4 N, as shown in Fig. 8a.The parameters of the sands are the same as those in Table 2, where both the normal and shear contact bond strengths (n_bond, s_bond) of the particles are 6 N. The numerical failure process development under loading is illustrated in Figs. 9 and 10.It is noticed that, on top of the slope crest, the region below the loaded surface has deformed to form a depression zone from the DEM modeling, which is also observed from the test as shown in Fig. 9e.Sand particles are triggered to move down the slope, dragging more and more sands downward.The depression zone develops larger and deeper, accompanied by considerable settlement at the top of the slope, and the inclined slope face moves forward, with an upheaval at the slope toe as shown from the xy direction view in Fig. 10d.
A slope with soil particles under higher bond strength is also simulated, where the normal and shear contact bond strengths of particle values are increased to 60 N, which is ten times larger than that in the former case.The results of the comparison are shown in Fig. 11.It is observed that noticeable collapse has taken place in case 1, with forward movement of the slope body and an extruding slope toe.On the contrary, the slope remains stable if the bond strength is 60 N under the loading of 500 N, as shown in Fig. 11b.For applied load equal to 2000 N, both cases of slopes experience failure, as the bond strength is completely destroyed by the applied load.It is demonstrated that larger bond strength between soils followed by slower and smaller failure can withstand relatively larger loads, so that the soil is more stable under such an external load.
Soil particles at the left corner, right corner and middle part of the slope crest are selected to monitor the position and velocity change in the z direction (vertical direction), re-   spectively, which are illustrated in Figs.12-13.For the right corner of the slope crest as shown in Fig. 12, the sand at the crest moved down with a high velocity at the beginning, so the displacement increased at a rapid rate within 2 × 10 5 time steps.Such a high initial velocity is normal even for the laboratory test, as the plate is not fully in contact with the soil at the beginning of the test.A small amount of plate movement is required so that the plate is in full contact with all the soil particles at the top of the slope.The increase in the vertical displacement became smooth afterwards with several fluctuation points (a typical phenomenon from DEM calculation).
For the left corner of the slope crest, the sand experiences the same trend as that for the right-hand side.In theory, the problem is symmetric, as the loading plate is located at the center.
In the microscopic view, there will be minor deviation from symmetry, as it is virtually impossible to keep all the soil particles truly symmetrical in the generation of the model or for the actual physical model.Such minor deviation from symmetry will however become negligible as the movement of the plate increases.
For the slope crest corner as shown in Fig. 12, it can be noticed that the particle moved down firstly at the corner of the crest where a crack is generated at the very beginning of the slope failure process.With further increase in the displacement of the loading plate and when a critical failure surface is generated, an instant collapse is observed from the great fluctuation of these history records.
At the middle of the crest of the slope as shown in Fig. 13, the situation is different in that a more rapid and less moderate trend can be noticed.The displacement reaches to more than 0.2 m at the middle of the crest compared with a displacement of almost 0.17 m at the crest corner.However, the particles at the corners move faster than those in the middle in the beginning, which are demonstrated by the sharp increase in the settlement in Fig. 12. Figure 13 shows that the loading at the middle of the slope top has pushed the sands at the middle of the slope crest downward and accumulated at the toe of the slope.
The front and side views for the final failure of the laboratory model test and DEM modeling are shown in Figs. 14 and 15.A steeply inclined failure surface (triangular failure zone underneath the loading plate), which is followed by a curved failure surface, can be observed from Fig. 15b and c.The failure mechanisms are marked by the equivalent slip surface, similar to that by LEM and SRM.The slip surface is determined from the velocity/displacement for all the particles.It is noticed that there is a narrow band of particles where the movement is significantly more than the adjoining particles, and the locus of this band is drawn as the slip surface.The failure profile of the physical test matched reasonably well with that as predicted from the present simulation as shown in Figs. 14 and 15a.At the beginning of the laboratory test, cracks first appear along each crest corner in the laboratory test, which is the effect due to the side plate of the model.In the particle flow simulation, the loading plate is applied with a small velocity of 5 × 10 −5 m per step.Besides the triangular failure zone underneath the loading plate, there is an obvious failure zone at the middle of the slope top, and the soil moved downward and accumulated at the bottom of the slope, which was also observed from the physical test as well as the DEM modeling.Based on the observation and comparison above, it can be concluded that DEM can give a very good qualitative description of the complete failure process, which was not possible with the strength reduction method (SRM) or the limit equilibrium method (LEM).
The relations between the applied load and the displacement curves of the slope at different measuring points from laboratory and numerical results are illustrated in Fig. 16.For the right-hand side, three groups of curves of DEM results were generated from the measuring points of lower position, upper position of slope surface and middle crest, respectively.For the soils at the lower position of the slope surface, there is a rapid increase in the loading applied before the displacement has reached 10 mm, which can be noticed in green and light blue curves in the figure.The applied load reaches the peak value of 37 kN at about 10 mm settlement, and then the loading first drops to around 25 kN, followed by the second drop down to around 10 kN till eventual failure.For the soil at the upper position of the slope surface, the trend is similar to that of lower surface points, but the settlement largely develops to approximately 14 mm under maximum loading.For the soil at the middle of the slope crest, the maximum load is achieved at a larger displacement of 26 mm as compared with the slope surface measuring points, and soil particles experience longer-time settlement with significant displacement values in the whole failure process.
In conclusion, the whole slope undergoes two major failures/stress redistributions in the test and numerical modeling when the bearing loading is 25 kN at three-quarters of the displacement, as well as 10 kN at the end of the test.Under these two conditions, the loading force remains unchanged, while displacement increases, and the reason is that the slope reaches a critical point and the failure surface is generated, with major stress redistribution taking place within the soil mass.The settlement at the lower parts of the slope surface appears earlier than that at the crest; however, the maximum final settlement of 40 mm is located at the slope's middle crest, which is underneath the loading plate and develops larger than the lower positions at the slope surface.In comparison with the laboratory test results in three groups of curves from the left-hand side of Fig. 16, also as shown in Fig. 7, the trends of the three curves are qualitatively similar, and the final failure loading of 35 kN from the laboratory test was also close to the value of 37 kN from the DEM analysis.DEM can hence give a relatively good qualitative assessment of the slope failure process up to the initiation of the failures.
For the post-failure results, the matching between the DEM model and the laboratory test is less satisfactory quantitatively.Gardiner and Tordesillas (2005) have given some discussion about the treatment for the loss of contact, which is not considered in the present study due to the lack of suitable material parameters.In theory, better matching may be achieved by choosing other types of contact rules or other micro-parameters, and the micro-parameters may even be changed after the maximum load has been reached.The authors have not carried out such numerical calibration for this study, and such post-failure parameters/contact rules are difficult to assess unless the corresponding laboratory tests are carried out, which are practically not possible for normal engineering problems.Even though the matching is not good quantitatively, the load displacement relation from DEM analysis and the model test results are similar qualitatively.Based on the comparison as mentioned above, the results of numerical modeling are considered effective and acceptable for engineering application.

Conclusion and discussion
In this paper, a laboratory slope model test is carried out with the corresponding 3-D numerical simulation using the discrete element method.The complete progressive failure mechanism of a slope under a local surcharge is determined both from the model test and numerical analysis for comparison.Furthermore, both qualitative and quantitative studies for full-scale 3-D numerical modeling have been carried out in the present study, for which it is virtually not found in the literature.
For the physical modeling of soil slope under locally external surcharge, the main observations from the failure process are as follows.The first crack was found in the surroundings of the loading plate, and it extended towards each boundary with an angle of approximately 45 • .After 4 h hydraulic loading, the sands at the middle top of the slope is pressed to form a big hollow zone, while the largest cracks are generated at the critical failure surface within the slope body.The soil mass contained within the critical failure surface is then pushed down with a complete failure, but the failure process is a gradual process with major stress redistribution and crack propagation.The failure surface after removal of the failed mass as well as the numerical modeling are found to fit nicely with the classical general shear failure mechanism as discussed by Cheng and Au (2005).
For the 3-D analysis of the particle flow model under local loading, the failure process of the particle flow slope model is generally similar to that of the physical model in most respects.It is observed that an obvious collapse is noticed at the middle of slope crest where the loading plate falls down gradually, forming a depression zone at the top.At the same time, sands at the two crest corners fall with a high velocity at the beginning and slow down to achieve equilibrium at the later stage.The failed soil mass moves down and accumulates at the slope toe.Many surface cracks directed towards the loading plate are found during the loading process.Finally, a critical failure surface basically conforming to the classical theory has formed, followed by a major collapse of the soil mass.
In comparing the physical modeling and numerical simulation, deformation at the center of loading plate develops more than the lower parts of the slope surface in laboratory tests, while the settlement underneath happens earlier than that in the middle, both in the laboratory test and the nu-merical simulation.The relations between the force against displacement curves are basically similar between the physical model test and the numerical model, and the final failure loadings are also very close, which is not accomplished in other previous studies.Although there are some discrepancy between the load-displacement results for the DEM analysis and laboratory test after the initiation of failure, DEM can still produce a reasonable qualitative assessment of the postfailure response of the slope for the considerations of the engineers.These results have useful contributions to the better understanding of the complete slope failure mechanism and the effect and spread of the slope failure, which is not possible for other classical methods.
Through the present work, the method to define a 3-D slope problem, the selection of the micro-parameters, the various precautions required for a good 3-D numerical quantitative modeling and the suitability of DEM in the assessment of the complete failure mechanism of slope have been assessed.If the numerical modeling is continued (under a different lengthy article), the analysis will cover the spread of the debris flow, and the effect of the slope failure can be assessed in greater detail.Currently, the post-failure analysis is not carried out in Hong Kong, except for some critical locations where there is a high risk to human life, but all the post-failure analyses carried out in Hong Kong are limited to 2-D analyses only (nearly the same situation in other countries).The present laboratory and 3-D numerical modeling can also complement the inadequacies for the typical slope hazard assessment, which will be useful for critical slopes.

FigureFigure 1 .
Figure 3 Particle size distribution Figure 4 Shear-Normal stress relation for sand

Figure 3 FigureFigure 2 .
Figure 3 Particle size distribution Figure 4 Shear-Normal stress relation for sand

FigureFigure 3 .
Figure 3 Particle size distribution Figure 4 Shear-Normal stress relation for sand

Figure5Figure 5 .
Figure5 LVDT at top and sloping face of the model test

Figure 6 .
Figure 6.Failure process of the soil slope under increasing loading.

Figure 7 .
Figure 7. Loading force against the displacement of the slope surface.

Figure 8
Figure 8 Two loading patterns of simulation models: (a) applying the force on the raft footing; (b) adding velocity on the loading wall.

Figure 8 .
Figure 8. Two loading patterns of simulation models: (a) applying the force on the raft footing; (b) adding velocity on the loading wall.

Figure 9 .
Figure 9. Failure process of the numerical slope model under the loading raft.

Figure 10 Figure 10 .
Figure 10 Failure process of numerical model under loading raft in XY direction Figure 10.Failure process of the numerical model under the loading raft in the xy direction.

Figure 11 Figure 11 .
Figure 11 Eventual failure of two modeling cases under loading raft in XY direction(with external load 500N)

Figure 12 .
Figure 12.Vertical position history for balls at the right corner of the crest.

Figure 13 .
Figure 13.Vertical position history for the balls at the crest middle.

Figure 13 Figure 14 .
Figure 13 Vertical position history for the balls at the crest middle

Figure 15 .
Figure 15 Side view for the final failure from laboratory test and DEM modelling Figure 15.Side view for the final failure from the laboratory test and DEM modeling.

Figure 16 Figure 16 .
Figure 16 Loading force against displacement curves of the slope at different measuring points from Laboratory and Numerical results Figure 16.Loading force against displacement curves of the slope at different measuring points.

Table 1 .
Shear strength parameters of the river sand.

Table 2 .
Microscopic parameters of the sands for particle flow analysis.