Failure mechanism of boulder-embedded slope under excavation disturbance and rainfall

Due to the uneven weathering of rocks, boulders can exist inside a slope, making the deformation and failure mechanism of the slope very complex. By analyzing the failure characteristics of a boulder-embedded slope under alternating excavation and rainfall, two classical instability modes are proposed, i.e., boulder instability and soil instability. For the soil instability, three failure processes may occur, including the sliding surface above the boulder, the sliding surface below the boulder, and the sliding surface intersecting the boulder. Meanwhile, the interaction between soil and boulder can also vary during different failure phases. Furthermore, the slope sliding deformation, failure mechanism, and soil-boulder interaction are investigated by on-site monitoring and numerical simulation. The results show that the boulders play an anti-sliding role and block the formation of deep sliding surface, resulting in the shallow soil instability and local large deformation of the studied slope. Besides, during the slope sliding deformation, three failure processes of the soils appear one after another, and the soils may slide along the bottom or top of the boulders because of the hindering effect of the boulders.


Introduction
Boulder is a particular geomorphic element formed by the physical and chemical weathering of rocks. It has a roughly spherical or elliptical shape (Ollier 1971). As rainfall contributes to the rock weathering process, boulders exist widely in humid climates, such as the Penedo da Sobreira boulders in Spain (Alejano et al. 2010), the boulder clusters of the Serra da Estrela mountain in Portugal (Migoń and Vieira 2014), the half dome of Yellowstone National Park in the USA (Chapin et al. 2014), and a large number of boulders distributed in Malaysia (Alavi Nezhad Khalil Abad et al. 2016). The uncertain distribution and concealment of boulders in slope can make the engineering works with foundation pit and slope construction complex. Sometimes, they can even cause natural disasters such as rockfall and landslide, which seriously threaten human lives and infrastructures (Liang et al. 2014;Messenzehl and Dikau 2017;Gong and Tang 2017;Vick et al. 2019). Therefore, it is of great practical significance to assess the stability and failure mechanism of boulder-embedded slopes.
According to the position inside slope, boulders can be divided into three type, i.e., wholly exposed, partially exposed, and embedded (Liu et al. 2018). For the first type, the boulder is exposed on the slope surface with only a tiny amount of contacts with soil or other boulders. Due to the small constraints, the boulder is easy to roll or slide along the slope surface under the external forces and causes rockfall disasters (Alejano et al. 2015;Perez et al. 2019). For the second type, its lower part is buried in residual soil, while its upper part is exposed on the surface. Its stability depends on the buried depth and the strength of residual soil, and the failure mode is similar to the wholly exposed boulder. For the third type, the boulder is surrounded by soil, which means that the slope original stratum structure is changed to the soil-rock binary structure. Clearly, the large difference between soil and boulder can make the slope deformation and failure process very complex. The boulder may block the sliding surface of the soil mass, and the soil mass may also move the boulder. The interaction between them can lead to the soil failure, boulder instability, or mud-rock flow sliding, which brings great challenges to slope design and construction. However, there is little research on boulder-embedded slope. Liu et al. (2019) and Li et al. (2021) investigated the influence of the rotation angle of a single boulder and the different sizes and shapes of boulders on slope stability numerically. Their results show that the grain size of the boulder is positively related to the factor of slope safety. Although its rotation angle and shape have little influence on the stability, the shape and position of the sliding surface will be changed. Furthermore, they pointed out that the sliding surface can expand around the boulder. Liu et al. (2018) made a mechanical analysis on the buried boulder and suggested that the boulder can only lose its stability when the surrounding soil was damaged.
At present, the appropriate theoretical model and systematic mechanical analysis on boulder-embedded slopes are still lacking. If the relative slipping and detaching between soil and boulder cannot be considered appropriately, the hindering effect of the boulder could be only considered as the coordinated deformation between soil and boulder, leading to the misjudgment of slope safety. Meanwhile, the actual boulder slopes can be affected by many external factors, and the failure mechanism of boulder-embedded slopes has not been fully understood under the engineering and natural disturbances. Some researchers evaluated slope safety by low-discrepancy sampling (Hu et al. 2022) and supplementary sampling ) and proposed a probabilistic strategy for analyzing slope instability . In fact, the excavation and rainfall are the two main factors causing slope instability (Sun et al. 2019;Zhao et al. 2019;Li et al. 2019;Gong 2021;Han et al. 2021). For a boulder-embedded slope, excavation breaks the initial mechanical balance of the slope (Wang et al. 2016), resulting in the mass difference between the soils on the upper and lower sides of the boulder and thus causing the slope instability. Besides, rainfall can reduce the effective stress and shear strength of soils (Bishop and Blight 1963;Zhu et al. 2022), drive the expansion of sliding surface (Yao 2021), and cause soil failure or debris flow. Under the combined influence of excavation and rainfall, the failure modes of the slope will vary. In this study, the calculation formulas to determine the stability of boulder-embedded slopes are deduced, and two classical failure modes under excavation and rainfall, i.e., boulder instability and soil instability are proposed. Then, the real slope containing multiple boulders is analyzed to explore the mechanical interaction between soil and boulder, failure mechanism, instability pattern under alternating excavation, and rainfall.

Analysis of slope failure
In general, the influence of excavation and rainfall can be modeled by setting the proper boundary conditions and physico-mechanical properties. The excavation can trigger an unloading rebound effect, and the rainfall can reduce the strength and increase the weight of soils. For boulder-embedded slopes, the two disturbances will lead to the sliding deformation of soil and influence the force equilibrium of boulders. Under some extreme conditions, the boulder can become unstable. Because the large difference of the physical and mechanical parameters exists between soil and boulder, their interaction becomes nonnegligible on slope stability. To investigate the potential local and global instabilities affected by their interaction, it is necessary to further establish the effective mechanical models and discuss their roles in different failure modes. Figure 1 shows the mechanical model used in this study. In the analysis, several assumptions have been adopted: 1. A boulder exists in the slope overburden layer, and the soil mass remains stable initially. 2. The inclination angles θ of the overburden and the slope are the same. 3. The local coordinate system is defined as the X-axis aligns with the slope inclination and direct upwards, and the Y-axis is normal to the slope. 4. As a rigid body, the boulder will not deform in surrounding soils.

Mechanical model of boulder instability
Before excavation, the boulder is subjected to the gravity force G, and the surrounding soil pressure contains the decomposed force components N 1 along the Y direction and N 2 along the X direction. The friction force S acts at the interface between the soil and boulder, as shown in Fig. 1.
Since the boulder is firmly embedded in soil, it is assumed that the boulder is a rigid body and cannot deform or rotate within the overburden soil at the critical state of failure. According to Fig. 1a, the force equilibrium of the boulder before excavation can be expressed as follows: After excavation, the earth pressure is partly released on the right of the boulder, resulting in the reduction of resistant force along the slope inclination to N ′ 2 , as shown in Fig. 1b. If the overall resistant force is less than the downslope component of gravity, N ′ 2 , may turn into a downward thrust, and the boulder will slide downwards along the bedrock. Meanwhile, when the rainfall is active, it will increase the soil weight and reduce the soil strength and the soil-boulder interface friction to S' which further intensifies the boulder (1) instability. The calculation of total forces acting on the boulder along the X and Y directions under excavation and subsequent rainfall is expressed as follows: In Eq. 2, S' < S and | N ′ 2 | ≤ |N 2 |. N ′ 2 is preceded by a minus sign for a supporting force and a plus sign for a thrust force.
It can be seen from the analysis that the possible failure mode of the boulder is sliding along the bedrock layer. Under rainfall, the sliding resistant force of the boulder is the lowest, and it is prone to slide. Before excavation, the soil plays an anti-sliding effect by exerting embedded earth pressure and basal friction on the boulder. After excavation, the original force balance cannot be held, and part of the earth pressure is transformed into a downward active loading, leading to the initiation of boulder sliding.
(2) Figure 2a shows the forces acting on the overburden soil in the slope. The gravities of the excavation zone and upper slope soil are W 2 and W 1 , respectively. The stability overburden is affected by the gravity G and the anti-slide force F R of the boulder. Before excavation, it is assumed that the boulder remains stable, which will move only after the surrounding soil loses stability and slides. F R is equal to the sum of N 2 and the interface friction force S in the mechanical model of boulder instability. The straight-line fracture surface method has been employed to analyze the overburden soil stability, which assumes that the sliding surface is approximately straight, as shown in Fig. 2a. In this research, the fracture surface  is assumed to have the length L and occurs at the interface between the overburden and bedrock layers. The cohesion and internal friction angles of soil in natural state are c and φ, respectively. The sliding and resistant forces in the initial stable state can be calculated as follows:

Mechanical model of soil instability
After excavation, the upper overburden soil loses the effective support and can slide along the weak surface. The subsequent rainfall can reduce the shear strength of the soil and increase its bulk density, aggravating the failure process, as shown in Fig. 2b. Furthermore, the failure of the overburden soil can be divided into three processes as follows: (1) the sliding surface occurring above the boulder, (2) the sliding surface intersecting the boulder, and (3) the sliding surface occurring below the boulder.
For the first failure process, since the boulder is not involved, it is the same as the general soil slope failure. The sliding and resistant forces of the slope under excavation and rainfall can be expressed as follows: where W and W' are the natural and saturated unit weight of the sliding soil, L' is the effective length of the fracture surface, and c' and φ' are the effective cohesion and internal friction angle of the saturated soil. As the sliding surface does not connect with the boulder, only the shallow surface layer and the frontal part of the soil slip. This local failure is the most likely to occur than the other modes.
For the second failure process with the sliding surface intersecting the boulder, the sliding and resistant forces can be expressed as follows: where N ′ 2 is the boulder's anti-sliding force. When the surrounding soil is saturated, F ′ R < F R . As the friction between the saturated soil and the boulder is weakened, this force is also reduced. In this mode, since the sliding surface intersects the boulder, the soil mass needs to overcome the strong anti-slide force of the boulder.
For the third failure process where the slip surface is under the boulder, the slope stability calculation formulas can be expressed as follows: It can be seen from the above formulas, the boulder's anti-sliding force F R disappears completely when the sliding surface is above and below the boulder. The damage caused by the infiltration of heavy rainfall or the rise of groundwater level can lead to the weakening of the soil from inside and form a wide range of deep-seated sliding surface. The potential slope failure and landslide hazards show a great potential to cause serious damages.

Numerical method
According to the analysis of slope failure mode, the boulderembedded slope can fail as either the boulder instability or soil instability. Depending on the external loading condition and the location of sliding surface, the boulder and the soil can play diverse roles. In the slope analysis, rainfall is considered after the excavation. The elements without seepage properties, such as boulders and lower bedrock, cannot be directly analyzed in the rainfall infiltration simulation. Therefore, the separate simulation and the superposition analysis method (Yu et al. 2021) have been used to solve this problem. Besides, the relative deformation between boulder and soil can be modeled by the interface elements. The calculation process is shown in Fig. 3, and the specific steps are described as follows: 1. The numerical model is established according to the real size and stratum distribution of the slope. Then, the corresponding material parameters are assigned to each layer and boulder, and the contact elements are set at the interface between boulder and soil. 2. The saturated-unsaturated seepage theory is adopted to analyze the water seepage in slope under the rainfall condition. Clearly, the non-seepage characteristic mesh groups such as boulders, interfaces, and bedrock are deactivated in seepage analysis. The rainfall intensity is applied to the excavated surface and slope top to obtain the pore water pressure and transient saturation zone distribution. 3. According to time sequence and pore water pressure, the specific seepage results are selected and extracted to create the pressure boundary conditions, and the distribution of transient saturated zone will be recorded. The physical test, parameter inversion, and other methods are adopted to obtain the physical and mechanical parameters of the saturated soil. 4. The excavation-induced stress field is analyzed by the construction stage and strength reduction method during which the boulder, interface, and bedrock grid group are reactivated sequentially. After excavation, the additional stages are set based on the previous selected seepage calculation results. In these stages, the corresponding boundary conditions of pore water pressure are activated, and the element parameters in the transient saturated zone are replaced with saturated parameters. Meanwhile, when applying new pore water pressure and modificated boundaries in each stage, the boundary of the previous stage needs to be passivated to avoid repeatedly applying pore water pressure and changed parameters.
In this way, the mesh groups without seepage characteristics such as boulders, bedrocks, and interfaces can participate in the stress calculation at rainfall stage. Meanwhile, as the deformation analysis is not needed in the seepage simulation, the change of soil porosity and the lower support loss caused by the deactivated boulders and bedrock will not affect the seepage calculation. This numerical approach separates the analysis of stress and seepage fields. The rainfall infiltration is simulated by applying pore water pressure to the elements and changing their mechanical parameters, which solve the problem that the materials without seepage characteristics cannot participate in the rainfall analysis under the alternating excavation-rainfall conditions.

Establishment of slope model under excavation and rainfall
To verify the accuracy of the theoretical analysis and the feasibility of the developed simulation method, a cutting slope containing multiple boulders caused by excavation and rainfall is analyzed. The numerical results are compared with the theoretical analysis and on-site monitoring data. Furthermore, the failure mode, instability mechanism, and interaction between soil and boulder are discussed.

Geological survey
The highway slope is located in Guangdong Province, China. It contains three main sections, i.e., K288 + 850 ~ 980, K288 + 980 ~ K289 + 180, and K289 + 180 ~ 420 as shown in Fig. 4a, which are situated at the junction of the Guangdong-Hunan fold belt with a typical geological structure of syncline. A normal fault with an occurrence of 177°∠87° passes through the southeast part. According to the measurement of the exposed bedrock, the attitude of the stratum is 73~80°∠15~37°, which belongs to the gently anti-dip stratum and is conducive to slope stability. The annual rainfall in this area is 1400-2000 mm and mainly occurs in spring and summer. No groundwater is found in the slope during drilling exploration. The geological data shows that boulders are deposited in the slope area. To determine the location and dimension of boulders, 5-geophysical prospecting lines containing 22-logging wells (BZK01 ~ 22) are set along the slope strike direction. Meanwhile, the high-density electrical method is used in these wells for exploration. The results show that the slope is covered entirely by silty clay, sandy clay, and fully to moderately weathered granites from top to bottom. The boulders are mainly located in sandy clay and fully weathered granite strata. Their shapes are basically ellipsoid or sphere, and their sizes vary from 0.8 to 15 m, as shown in Fig. 4b, c.
Since sliding failures occurred in the middle of the K288 + 980 ~ K289 + 180 slope section under excavation and rainfall and the boulders are mainly distributed near the K289 + 100 section, the K289 + 100 section is selected as the slope case for the study. The slope has a height of 54 m, an inclination of 265°, and a natural slope ratio of 1:2.5. The initial design involved the four-level excavations. The ratios of the lower first and second slopes are 1:1 and 1:0.8 with the height of 10 m and 4.7 m, respectively. The ratios of the upper third and fourth slopes are 1:1.2 and 1:1.3 with the designed height of 4.7 m and 7.2 m. The platforms are gradually narrowed from low to high with the widths of 9.5 m, 6.2 m, and 3.5 m. The slope morphological and geological conditions are shown in Fig. 5, in which the weathering degree of granite is classified according to the standard for engineering classification of rock mass (GB/T50218) published by the Ministry of Housing and Urban-Rural Development, of the People's Republic of China (2014). Two large elliptic boulders exist in the fully weathered layer. The upper one is smaller and close to the interface between strongly weathered granite and fully one. The lower one is twice times larger than the upper boulder and locates near the toe of the second slope.
Deformation, failure process, and displacement characteristics

Slope deformation and failure
The K288 + 980 ~ K289 + 180 slope section exhibited a typical progressive failure induced by rainfall infiltration on the exposed slope after the excavation, which delayed the progress of slope supporting and protection. The slope construction started in December 2015, and all four slope excavations were completed at the end of February 2016. In mid-March of the same year, the continuous rainfall occurred. After 5 days of continuous rainfall, the slight tension cracks appeared on the first platform, as shown in Fig. 6a. On March 27, 2016, rainwater began to seep out at the slope toe and formed a small ditch, as shown in Fig. 6b.
When April comes, the rock and soil masses were softened and deformed due to the heavy rainfall infiltration. The water accumulated gradually at the slope toe, and the soils on both sides of the ditch were in the saturated state, as shown in Fig. 6c. On April 18, a shallow landslide occurred in the first slope. The sliding area extended from the toe to the top platform, but did not reach the lower boulder, as shown in Fig. 6d. On April 20, the lower part of the second slope was dragged and damaged by the sliding bodies, with the cracks appearing on the second platform. By the end of the rainfall on April 30, the sliding area developed to the upper part of the second slope.
Landslides were also observed in the fourth slope. During rainfall, the surface runoff from the slope top continuously scoured the fourth slope. The field survey on March 25 showed that part of the soils had been washed away, as shown in Fig. 7a. Affected by continuous rainfall, many tension cracks appeared on the slope crest in April and continued to expand and develop. By April 21, the fourth slope slid as a whole. The sliding area started from the tension crack at the slope crest and extended to the toe of the fourth slope. The edge dislocation formed at the original tension crack, as shown in Fig. 7b. Simultaneously, the deformation and extrusion of the landslide body promoted the initiation of shear cracks on both sides of the slope.

Characteristics of deep slope deformation
The slope monitoring was conducted from March 18 to April 29, 2016. The inclinometers were buried with the logging wells to investigate the stability of boulders and soil masses by measuring the horizontal displacement along the slope dip direction. Four measuring points were arranged at the K289 + 100 section, among which ZK1-1, ZK1-3, and ZK1-4 were located at the slope toe, secondary platform, and slope crest, corresponding to the BZK10, BZK11, and BZK12 logging wells, respectively. The lower boulder is situated at the south of the ZK1-1 pipe, with the depth of 6 ~ 16 m. The upper boulder is located at the north of the ZK1-4 pipe, with the depth of 10 ~ 15 m. ZK1-2 is an additional measuring point, and the inclinometer pipe was buried at the middle of the first platform. The four measuring points are used to monitor the displacement of slope toe, first slope, second slope, and fourth slope. The layout of the measuring points and the relative positions of boulders and pipes are shown in Fig. 8.
To ensure the accuracy of monitoring, the data was recorded at every 0.5 m in depth, and the initial position of each pipe was set as the zero displacement. Each pipe was partially embedded in strongly weathered bedrock. Particularly, ZK1-3 and ZK1-4 went through the boulder area. To facilitate the data analysis, the corresponding areas of each pipe at different depths have been marked on the displacement curves as illustrated in Fig. 9. In Fig. 9, the lower and upper boulders point out the measured displacements of the inclinometers at the same depth as the boulders. Fig. 4 The topography and geology of the slope before excavation a terrain and section division, b the vertical geological profile I, and c the vertical geological profile II (the gray filling represents granite boulders) ◂ As illustrated in the displacement curves, although the rainfall started on March 15, 2016, apparent slope deformations only occurred after 10 days on March 25, 2016. The slope toe deformed significantly, and the largest deformation reached 32.88 mm on the surface. The maximum displacement at the first, second, and fourth slopes were 22.11 mm, 19.32 mm, and 21.62 mm. The soil around the boulder also began to deform, and its displacement gradually decreased with the depth, which indicates that the slide and deformation were activated. When April came, the heavy rainfall aggravated the slope deformation. The monitoring data on April 1, 2016, show that the average displacement at the slope toe and the middle-upper part of the first slope were 40.89 mm and 29.15 mm, respectively. The curves of the second and the fourth slopes change significantly. The displacement of the second slope increases gradually with the depth, demonstrating the traction deformation of the first slope. A data gap occurred at the middle fourth slope, indicating the shallow area of the fourth slope began to slide. By April 8, 2016, the traction sliding in the middle-lower part of the second slope was intensified, two data gaps appeared at the middle and the junction between the lower part and the boulder area. Meanwhile, the data gap of the fourth slope began to move downward, indicating that the failure of the slope began to expand from shallow to deep.
On April 18, 2016, the first slope collapsed. On April 20 and 21, 2016, the lower parts of the second and fourth slopes slid. The monitoring results on April 22 show that the average displacements of the collapsing area on the slope toe and the first, second, and fourth slopes were 83.9, 78.78, and 78.97 mm, respectively. The damage range of the fourth slope was extended from 0 ~ 5m to 0 ~ 8m, and the second data gap also appeared at the junction between the lower part and the boulder area. The displacement at the depth of the two boulders is significantly lower than the upper area, shows that the boulders play an anti-sliding effect. By the end of the monitoring, the deformation area of the second slope expanded to the upper part with the average displacement of 74.84 mm.

Numerical configuration
Based on the site and deep displacement characteristics, the first, second, and fourth slopes were instable under rainfall, exhibiting shallow sliding failure patterns. The deformation of the soil mass on the side of the boulder was significantly lower than the upper areas, indicating that the boulder had not slipped but might have an anti-sliding effect. Following the theoretical analysis, two-slope failure modes are possible.
In the first mode, the sliding surface was above the boulder, and the failures of the first, second, and fourth slopes were all local sliding. The plastic and the stress concentration zone were close to or coincide with the damaging zone. The anti-slide effect of boulders might be caused by the significant displacement difference between the upper and lower soil mass. In the second failure mode, the sliding surface intersected with the boulder. Namely, the plastic zones were located inside the slope. Because the boulder has fully exerted its anti-sliding capacity, the plastic zones cannot connect at boulder's position. The soil mass was damaged in the shallow area due to the initiation of landslides. On the basis of the developed simulation approach, the above cutting slope containing multiple landslides caused by excavation and rainfall has been analyzed. The numerical results are compared with the on-site monitoring and the theoretical failure modes. The failure mode of cutting slope and the soil-boulder interaction mechanism are further discussed.

The slope analysis model
In this study, the two-dimensional finite element software Midas GTX (MIDAS Information Technology Co. Ltd. 2020) was used to analyze the slope failure. The numerical model consisting of 11,067 elements was configured according to Fig. 5. The vertical sliding bearings were applied to the left and right boundaries, and the fixed bearings were applied to the bottom boundary. The slope safety factor was calculated by the strength reduction method. Since the traditional failure criterion of numerical convergence or plastic zone penetration may cause infinite displacement or plastic zone discontinuity to the boulder-embedded slope, the displacement mutation was used as the failure criterion in this research. It is challenging to calculate the seepage analysis under rainfall conditions when affected by materials with no seepage characteristics such as the bedrock and boulders. Therefore, the separate simulation and superposition approach were used to analyze the excavation and rainfall separately, and the pore water pressure was applied to the corresponding rainfall stage of the stress field. The input parameters of physical mechanics and seepage characteristics are discussed as follows:

Geomechanical parameters
Since the relative sliding and detaching between the boulder and surrounding soils may exist, it is necessary to set a suitable contact type at the interface to avoid the compatibility problem. Based on the experimental study of soil-rock interface shear behavior (Marinho and do Amaral Vargas 2020;Xu et al. 2013), the Coulomb friction elements were set on the soil-boulder interface. The cohesion force and internal friction angle of the elements were assumed to be 5 kPa and 30°. The normal and tangential stiffness moduli setting as 5 × 10 6 kPa and 1 × 10 6 kPa, respectively. The natural and saturated state parameters of each layer of the slope were selected according to the site investigation report provided by the construction unit. Τhe parameters of the boulders were assumed based on a large number of excavated boulder data (Du et al. 2013), as shown in Table 1.

Seepage characteristics
The parameters of soil seepage characteristics include the coefficient of unsaturated permeability and the soilwater characteristic curve. The former determines soil permeability, and the latter determines its hydraulic gradient distribution in an unsaturated state (Fredlund  . 1998). These are all highly nonlinear functions and are difficult to calculate directly. Therefore, the Van Genuchten model (Parker et al. 1985) was used to calculate the water content function and then combined with the coefficient of saturated permeability to obtain the permeability coefficient function in different states. The governing equation of the V-G water content function model is as follows:  where θ w is the volumetric water content; θ r is the residual volumetric water content; θ s is the saturated volumetric water content; ψ is the negative pore water pressure; and a, n, and m are the curve fitting parameters.
Because the strongly to moderately weathered granites are relatively deep and have little influence on the slope stability, their permeability characteristics are not considered. In accordance with the geological survey data and related references (Lin et al. 2013;Tan et al. 2017;Yu et al. 2020;Zhao et al. 2017), the VG model parameters and saturated permeability coefficients K s of the silty clay, sandy clay, and fully weathered granite are listed in Table 2.

Excavation and rainfall analysis conditions
According to the steps of the developed simulation approach, the Midas GTX transient seepage module was used to analyze the rainfall process. The daily precipitation of the simulated rainfall was input into the software. Since no groundwater was found on site, the initial groundwater level was set at the bottom of the slope model to meet the transient seepage analysis conditions without affecting the results. The rainfall lasted for 46 days, and the accumulated (7) w = r + s − r 1 + a n m rainfall was 478.3 mm. The daily precipitation of the rainfall is shown in Fig. 10.
To analyze the influence of rainfall infiltration on the slope stress variation, the pore water pressure of slope elements, extracted from the rainfall seepage analysis results on March 25, April 6, 18, and 30, was created as four pressure boundary conditions. The transient saturated zone of these four days' rainfall seepage analysis results was also recorded. As the volume of excavation was small, only one stage in the construction phase was set by deactivating the mesh group of the soil masses in the excavation area for stress analysis. After excavation, four rainfall stages were set to simulate the influence of pore water pressure and shear strength parameter weakening of saturated rock and soil on slope stress variation. The pressure boundary conditions were applied to four stages in chronological order of the seepage analysis results. Meanwhile, the natural parameters of the soils in the transient saturation zone were replaced by saturated parameters. In the analysis of each rainfall stage, the pore water pressure and parameter modification boundaries of the previous stage were deactivated.

Seepage analysis
As the rainfall has lasted for a long time and the simulated rainfall affects the stress analysis through the pore water pressure and parameter weakening boundary, this section will mainly discuss the four seepage analysis results applied to the stress field in a chronological order. The distributions of pore water pressure on different days of rainfall are shown in Fig. 11. The white dotted lines represent the infiltration lines, and the dusty blue area is the lower bedrock. Figure 11 shows the distribution of pore water pressure in the early, middle, late, and final stages of the rainfall. Due to the low rainfall intensity, the infiltration of rainwater only formed a transient saturated zone on the surface in the early stage (March 25, 2016). Below this zone, the pore water pressure of the soil was all negative. In the middle stage (April 6, 2016), the rainfall intensity gradually increased, and the loss of slope matrix suction caused the transient saturation zone to extend downward. When entering the late stage (April 18, 2016), the rain began to subside, but the previously infiltrated rainwater gathered in the silty clay layer and formed a saturated zone on the slope top. Affected by the steep slope inclination, most rainwater flowed down to the middle and lower part of the slope through the shallow area, which is also consistent with the phenomenon that the upper fourth grade slope was scoured by surface runoff, as shown in Fig. 7a. This rainwater infiltrated into the fully weathered granite through the high permeability silty and sandy clay, making the expansion of the transient saturation zone in this area much higher than the others. In the final stage, the rainfall intensity increased again. The lower boulder blocked the infiltration of rainfall at the slope front, so that the soil below the first and second platforms was all saturated. This can be confirmed by the phenomenon that a large amount of rainwater seeped out at the slope toe, as shown in Fig. 6d. Influenced by the surface runoff and infiltration capacity of rock and soil, the matric suction of the soil mass at the slope front rapidly dissipated and turned to a positive value during the entire rainfall process. However, the middle and upper soils did not change after forming a transient saturated zone in the shallow area.

Slope deformation and failure
Slope stability analysis Figure 12 shows the evolution of the slope effective plastic strain under the excavation and the subsequent rainfall of March 25 to April 30, 2016. For clarity of the presentation, the areas without obvious plastic strain below the fully weathered layer are not displayed in the figure. It can be seen from the figure that the potential failure mode of the slope is deep-seated sliding. After the excavation, two unconnected plastic zones from the slope toe to the lower boulder and to the upper boulder formed inside the slope. The sliding surface shows a trend to bypass the boulders and extend along the interface between fully and strongly weathered granite, indicating the hindering effect.
In the early period of rainfall (March 25, 2016), a plastic zone formed at the slope top and extended to the upper boulder, and the lower two plastic zones also expanded. In the middle stage of rainfall (April 6, 2016), the deep plastic zone started to extend towards the excavated surface, and the stress was mainly concentrated on the both sides of the upper boulder and the slope toe. When entering the late rainfall stage (April 18, 2016), the plastic zone of the slope grown significantly; the upper two sliding surfaces were connected along the bottom of the boulder, while the lower sliding surfaces were still blocked. The plastic zone extended upward to the first and second platforms and the slope crest. A stress concentration belt formed between two boulders and the rest stress concentration area of the slope toe; the sides of boulders also show expansion. At the end of rainfall (April 30, 2016), the slip surface still did not penetrate along the bottom of the lower boulder, and the plastic zone no longer changed. Figure 13 shows the variation of the factor of slope safety in each modeling stage. The slope stability gradually decreased under the influence of rainfalls. According to Wang et al. (2017), the slope deformation can be divided into the creeping phase (1.05 < F s < 1.1), the extrusion phase (1.02 < F s < 1.05), the sliding phase (0.98 < F s < 1.02), and the sudden slip phase (0.95 < F s < 0.98). Under the initial and excavation conditions, the slope was in the creeping deformation phase. Because of the rainfall infiltration, the slope transformed into extrusion deformation and finally slid.

Characteristic of creeping deformation phase
As shown in Fig. 13, before the rainfall, the slope was in the creeping deformation phase. The slope was stable at the initial stage with a safety factor of 1.112. After excavation, the factor of slope safety dropped to 1.068. The horizontal deformation concentrated in the middle-lower part of the first and second slopes, with the maximum deformation of 16.4 mm at the slope toe, as shown in Fig. 14a. As the soil masses were relatively soft and the excavated mass was small, the slope did not show any unloading signs. Instead, the high stresses concentrated at the shallow areas due to the loss of effective support, as shown in Fig. 14b. At the bottom of the fully weathered granite layer, an intermittent plastic zone extended from the slope toe to the upper boulder and the maximum plastic strain were concentrated at the slope toe, showing the characteristics of the traction landslide.

Characteristics of extrusion deformation phase
After the rainfall infiltration, the slope reached the extrusion deformation stage, characterized by the rapid extension of the plastic zone in soils. On the March 25, 2016, the factor of safety was reduced to 1.044. Three deep sliding surfaces formed and were separated by the two boulders, as shown in Fig. 12b. Apart from the first and second slopes, large displacements were at the shallow area from the slope crest to the upper fourth slope. At this point, three deformation areas initially formed: the middle-lower parts of the first and second slopes and the slope crest to the upper fourth slope, as shown in Fig. 15a. The maximum displacement values of the three areas were 31.5 mm, 20.1 mm, and 21.8 mm, respectively, which are basically close to the maximum values of 32.9 mm, 22.11 mm, and 19.32 mm of the monitoring points at the slope toe, the second platform and the slope crest on the same day, indicating that the simulation captured the main features of the slope failure effectively.
On April 6, 2016, the factor of safety was reduced to 1.025, and the slope approached the sliding deformation phase. The deformation area of the second slope extended to the middle-upper part, and the area of the crest extended to the middle fourth slope as shown in Fig. 15b. The maximum horizontal displacements of the three areas were 59.4 mm, 46.9 mm, and 40.5 mm, which are twice as high as the previous stage. The significant differences appeared between the displacement values of the deformation area and the rest. This phenomenon can also be seen from the deep displacement curve of the second and fourth slopes on April 8, 2016. The data gaps at the depth of the middle second, fourth slopes, and the junction between the lower second slope and the boulder area are coincide with the boundary locations of the deformation area. Combined with the distribution of the slope plastic zone at this moment, the plastic zone expanded, and the sliding surface on both sides of the upper boulder gradually connected; the sliding surface tried to bypass the boulder and went through the soil at the bottom, meaning the third failure process of the soil instability where the sliding surface occurs under the boulder.

Characteristics of sliding deformation phase
As the rainfall continued, the slope eventually slid. The factor of slope safety on April 18, 2016, was 0.997, which means the sliding deformation phase. The deformation area of the first slope extended to its top and connected with the second slope, and a landslide occurred from the lower second slope to the slope toe. Meanwhile, the area of the slope crest expanded downward to the upper boulder, and soil masses slide along the fourth slope toe, as shown in Fig. 16a. The comparison between the simulated and actual failures shows that the simulated failure at this stage is basically consistent with the collapse of the first slope on April 18, 2016, the failure at the lower second slope on April 20, 2016, and the landslide of the fourth slope on April 21, 2016. The plastic zone at this stage further expanded significantly, as shown in Fig. 12d. The fracture surfaces connected along the bottom of the upper boulder, showing the third failure process of the soil instability where the sliding surface occurs under the boulder. However, the large-scale deep landslide did not occur because the sliding surface was still blocked by the lower boulder. The stress concentration belt shows that the soil masses have been destroyed and attempted to move the boulder when the bypass is impossible. At this moment, the slope shows the second failure process of the soil instability, i.e., the sliding surface intersects the boulder. Due to the large difference of the mechanical parameters, the boulder cannot be pushed. The plastic zone expanded to the slope surface, which eventually led to the damage of soil masses above boulder, and the failure process is similar to the first one for soil instability.
By the end of the rainfall (April 30, 2016), the horizontal displacement did not change. The damage area of the second slope extended to the upper part, as shown in Fig. 16b. The deformation of the fully weathered granite on the sides of boulders is 3.71 ~ 38.9 mm, which is much lower than the 59.9 ~ 97.3 mm of the three damage areas. The soil displacement in these ranges gradually decreased with the depth, which is consistent with the measured displacement curves at the depth of the boulders. The boulders effectively restrain the deformation of deep soil mass.

The interaction between soil and boulder
According to the simulated failure process of the boulder-embedded slope, the underlying failure mechanism and interaction between soil and boulder are discussed in this section. Along with the construction stage, the slope shows the progressive failure characteristics. Based on the evolution of the displacement field, the local failure occurred on the slope under the alternating excavation and rainfall. In the excavation stage, the loss of the front soil Fig. 13 The factor of slope safety at each phase masses makes the slope form a free face and lose effective mechanical support, causing the stress concentration at the excavated slope surface and the exposure of the weak zone at the slope toe and further resulting in the deformation of the middle-lower first and second slopes. At the early stage of rainfall, the rainwater infiltrates into the shallow areas and forms a transient saturation zone, which increases the unit weight of soil mass and reduces its effective stress and shear strength and therefore causes the deformation of the slope crest to the upper fourth slope. In the middle of rainfall, the enhanced rainfall intensity makes the soil matrix suction lose rapidly. With the expansion of the saturation zone, the deformation area of the first and second grade slopes also extends to the middle, showing typical tensile failure characteristics. In the late to the final rainfall stage, the runoff accumulation and the blocking of rainwater flow by the lower boulder completely saturate the soil from the toe to the second slope and induce the collapse of the first slope and the traction failure of the second slope. At the same time, the runoff continues to scour the fourth slope and erodes the soil on the surface, making the deformation area extend downward and a landslide occurs from the crest to toe.
According to the evolution of the plastic zone, the shallow failure occurs on the slope since the deep sliding surface is blocked by the boulders. With the process of excavation and rainfall, three failure processes of soil instability appear one after another. In the excavation stage, the slope mechanical balance is broken; the plastic zone appears at the slope toe and expands upward to cause two sliding surfaces from the toe to the lower boulder and to the upper boulder, showing the characteristics of the traction landslide. In the early stage of rainfall, the rainwater drives the plastic zone from the upper boulder to the top of the slope to connect, thus resulting in an intermittent sliding surface from the toe and top of the slope. In the middle to late rainfall stage, the plastic zones on both sides of the upper boulder gradually penetrate along its bottom, showing the third failure process that the sliding surface bypasses the boulder. Meanwhile, the stress concentrates at the lower two sliding surfaces, presenting the second failure process of the soil moving the boulder. When the lower boulder is difficult to be pushed and the deep sliding surface cannot form, the plastic zones extend to the excavation slope surface and the soil masses above the boulder side, which can be manifested as the first failure process.
In accordance of the instability mode, the soil mass in the slope is a sliding body, which tries to bypass the boulder and slip as a whole, showing the typical hindering effect of the boulder. The boulders mainly played an anti-sliding role. However, they also had adverse effects on stability. The negative effects are reflected in the following aspects: (1) in the excavation stage, the boulders increase the difficulty of anchor/bolt construction and slow the progress of slope support and protection so that rainfall can infiltrate into the excavated bare surface; (2) during the rainfall infiltration, the lower boulder can block the infiltration rainwater flow at the slope front and make the soil masses of the first and second slope become saturated quickly. However, during excavation and rainfall, the two boulders can effectively block the formation of the deep sliding surface, reducing the deep deformation and forcing the soil masses to slide along the shallow area. Namely, the scale of sliding failure can be effectively controlled.

Conclusions
To reveal the nonlinear deformation and failure mechanism of a boulder-embedded slope under alternating excavation and rainfall, the related mechanical model has been developed. The main conclusions can be drawn as follows: 1. To assess the stability of a boulder-embedded slope under excavation and rainfall, a systematic theoretical analysis has been carried out through mechanical calculation. Two classical failure modes of boulder and soil were proposed with four representative failure processes. By taking an actual slope with two embedded boulders as a study case, the real failure process was reproduced numerically and compared with field investigation and measured deformation. The failure mechanism and interaction between boulder and soil were analyzed. The results show that the boulders block the formation of the deep sliding surface, causing the shallow soil instability. With the process of excavation and rainfall, three failure processes during soil instability appear one after another.
Funding This work was financially supported by the UK Engineering and Physical Sciences Research Council (EPSRC) New Investigator Award (Grant No. EP/V028723/1), for which the authors are grateful.

Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding author upon request.

Competing interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.