Soil-Rock Slope Stability Analysis under Top Loading considering the Nonuniformity of Rocks

Soil-rock slopes are widely distributed in central or western China. With the development of transportation, many subgrades are being built on mountainsides and therefore, slope stability has to be estimated under high loadings. To obtain better estimation results, a new rock contour establishing algorithm was developed, capable of considering interlock effect between rocks. /en, computed tomography (CT) and unconfined triaxial tests with ring top loadings were conducted. Based on rock distribution characteristics (obtained by CTphotos) and the appearance of shear failure surfaces in slopes under ring top loadings, four rock skeleton status and five shear failure surface developing models were introduced. Based on the developed rock contour establishing algorithm, ten groups (twelve models per group) were established and calculated by finite element method (FEM). After this, normalized ultimate loading increasing multiple N, which was the ultimate loading ratio of rock-containing slope to uniform soil slope, was introduced to evaluate the influence of rock distributions on slope stability. /e value of N was increased with the increase of rock content due to rock skeleton status. /e values of N in slopes with angular rocks were about three times higher than those with round rocks which was due to complex geometric shape and distribution characteristics of angular rocks. /en, considering different slope angles (50°–60°), rock contents (0%–60%), and rock shapes (round and angular), the ultimate loading increasing multiple N of soil-rock slopes under high loadings was calculated and suggested for engineering designs. Finally, based on the failure surfaces of numerical modes, three typical failure modes were developed, which could be reference for designers to deal with slopes.


Introduction
Soil-rock slopes consisted of low-strength soils and highstrength rocks [1][2][3] and are widely distributed around the world, especially in central and western China [4,5]. With the development of transportation in China, many highways or railways have been built on mountainsides (shown in Figure 1), among which some mountains contained soilrock mixtures. In designing subgrades on soil-rock slopes, designers usually assume slope as being composed of uniform materials ignoring the effect of rocks among slopes.
Strength reduce method (SRM) has been usually used in estimating slope stability. In estimations, where plastic belts are developed from the bottom to the top of slope, obtaining specific strength reduction coefficient of soil with SRM seems to be reasonable. However, it is different from slope instability due to additional loading, where plastic belt (or failure position) possibly starts from loading position, slope bottom, or both. As shown in Figure 2, failures appeared at the top of the slope (near subgrade) and the bottom of slope was intact. Hence, it is necessary to study soil-rock slope stability under top loading.
Developing suitable soil-rock slope analytical models is necessary to improve the stability analytical results of soilrock slopes. Hence, many researchers [6][7][8][9][10] have focused on the establishment of soil-rock slope models considering rock contents, shapes, distributions, sizes, and so on. According to the above rock characteristics, several random establishment algorithms have been developed for soil-rock slopes. Liu et al. [7] applied random media method (changing material properties in meshed elements) to develop a soil-rock slope model capable of taking into account the effect of rock content, but rock integrity could not be controlled in analyses. en, Liu et al. [4] applied random rock establishment algorithm which could consider the contents and shapes of rocks to study the stability of soilrock slopes with different rock contents, sizes, shapes, and relative distributions. Emad et al. [6] employed random rock establishment algorithm to study the stability of soil-rock slopes with different volumetric block proportions (VBP) and conducted low scale model tests to verify the obtained analytical results. ey concluded that VBP played an important role in the stability of soil-rock slopes and applied three-step conservative approach to estimating soil-rock slope stability. Based on different VBP values, Khorasani et al. [9] established 180 theoretical models to study the stability of soil-rock slopes considering the inclination of rocks in slopes. eir results showed that the safety factor of soil-rock slopes tended to maximize when greater differences appeared between the dip angle of rocks and slope angle. e above algorithms were not able to establish soilrock slope models capable of taking into account interlock effects [8] between rocks which had a great effect on the stability analysis of soil-rock slopes.
Slope stability [11][12][13][14] is an international problem and has been studied for decades. With the development of theoretical sciences and computer technologies, finite element method (FEM), extended finite element method (XFEM), limited equilibrium method (LEM), discrete element method (DEM), general particle dynamics (GPD), and peridynamics (PD) have been developed. All above methods have been employed in the investigation of heterogeneous slope stability and have obtained great achievements. FEM is a common method for heterogeneous slope stability analyses. Liu et al. [4] applied this method to study the stability of soil-rock slopes, through changing material properties to reflect material nonuniformity. XFEN is an extended version of FEM in fracture mechanics. Zhou and Yang [15,16] used this method to study the stability of rock slopes containing cracks and holes and concluded the development model of shear belts in slopes. LEM is a slope stability analysis method based on equilibrium equations. In analyses, the minimum safety factor of slope is calculated by adjusting the position and shape of sliding belt. Zhou and Cheng [17,18] applied this method to study the sliding surface and stability of Xiangjiashan landslide body and achieved good results. DEM is a discontinuous medium analysis method, which was applied by Ye et al. [19] to investigate slope stability under seismic loading, but its rationality was not verified. PD and GPD are often used to analyze the fracture development characteristics of rock-like materials. Bi et al. [20][21][22][23] applied these two methods to study rock-like material fracture development modes containing preexisting embedded flaw. Based on the calculation characteristics of above analysis methods and existing research results on soil-rock slope stability, considering compatibility with developed software for soil-rock slope models, FEM was finally used to analyze soil-rock slope stability in this paper.
Comprehensive understanding of soil-rock slope instability modes, which is significant in slope reinforcement, is of great importance for designers [19,[23][24][25]. However, most studies on slope instability modes have focused on rock slopes, slopes with weak slide belts, and those under extreme conditions. Han et al. [26] performed experimental tests to study rock slope failure modes under vibration and proposed three typical failure modes including crack development, slope spalling, and collapse sliding. Masoudian et al. [27] employed numerical methods to study the slip surface of slopes considering nonhomogeneity and rainfall and developed 6 slope instability modes for 0, 1, 2, 4, 7, and 10 days of constant raining. At shallow landslides near Changbai Mountain in China, Zhan et al. [28] deeply investigated slope failure mechanism and belt failure extent. However, it is blurred about instability modes of soil-rock slopes, especially under the effect of top loadings.
In this study, soil-rock stability has been analyzed under top loadings considering the effects of the contents  and shapes of rocks using FEM. Firstly, a new soil-rock slope establishing algorithm was developed which was able to consider the effects of single rock geometric shapes and rock distribution characteristics (especially interlock effect between rocks). en, CT and unconfined triaxial tests under ring top loading were conducted to analyze rock skeleton status and failure surfaces of soil-rock mixtures.
en, using the proposed soil-rock slope establishing algorithm, several studies with different slope angles considering different rock contents and shapes were conducted on soil-rock slope stability under top loading. Finally, based on statistical results, increasing multiple N of ultimate loading was suggested considering rock angles, shapes, and contents and three typical failure modes were developed for soil-rock slopes under top loading.

Shape Characteristic of Rocks.
In soil-rock slope stability analyses, the integrity and shape characteristics of rocks are of high importance [4,29]; therefore, it is necessary to well understand the shapes characteristics of rocks.
As shown in Figure 3, round rocks have high granule integrity and no defects appeared in rocks. Soil-rock mixtures with round rocks are usually found in strata exposed to water and wind transportations and sedimentation. During water transportation, sharp edges of rocks are washed and rubbed resulting in the rock shapes to be streamline with high integrity. During wind transportation, the streamline and integrity of rocks are mainly due to the effect of particle friction arisen from wind.
As shown in Figure 4, sharp edges and initial cracks appeared in angular rocks. Soil-rock mixtures with angular rocks usually existed in soil layers formed by mountain rock collapse, where large rocks break into small pieces along initial cracks or rock beddings. If complete breaking does not happen, some initial cracks are included in rock. Hence, angular rocks have sharp edges and sometimes include initial cracks.

Size Characteristic of Rocks.
e definition of rock size in soil-rock mixtures is different for different purposes [1,10,30]. Xu et al. [8,30,31] found that rocks with smaller size cannot affect the mechanical properties of soil-rock mixtures and rocks with greater size lead to serious size effect [32]. Hence, suitable rock size is important in the analysis of soil-rock mixtures. According to the sample preparation standard (GB/T 50123-1999) and existed references [1,33,34], Kalender et al. [2,8,35,36] concluded that, in soil-rock mixtures, minimum thresholds for soil and rock are suggested to be the minimum value of 2 mm and 0.05 Lc (where Lc is the characteristic length of model), and maximum thresholds are advised to be 0.2 Lc. Based on the advice, Zhang et al. [1] conducted triaxial tests with soil-rock mixtures containing 2-10 mm rocks (samples are 50 mm in diameter and 100 mm in height) and obtained good experimental results. e definition of rock size in soil-rock slopes is different from soil-rock mixtures, based on numerical studies on soilrock slopes [3,4,37], the threshold for rock and soil is usually limited to 0.1 Lc∼0.4 Lc (where Lc is the width of slope). When rock size is smaller than 0.1 Lc, rock distribution has less effect on slope stability. When rock size is greater than 0.4 Lc, rock distribution has a decisive influence on slope stability and failure modes and will lead to serious size effect. Hence, the suitable rock size in soil-rock slopes is in 0.1 Lc-0.4 Lc. Liu et al. [4] adopted the suggestion to study the stability of soil-rock slope and obtained good analysis results.

Establishment of Soil-Rock Slope Model
3.1. Establishment of Single Rock. Liu et al. [4] developed a random rock contour generation algorithm which was able to consider the shapes and sizes of rocks. In Liu's algorithm shown in Figure 5, a random circle was first established which was used to check existing rock intrusions. en, the new rock A was established in circle where the maximum diagonal points (P A2 and P A4 ) of rock were located on circle circumference and other points (P A1 and P A3 ) were inside the circle. e algorithm could well consider rock shape and size but was not able to take into account interlock effect between rocks, like rock B in Figure 6. Hence, a new random rock contour generation algorithm was developed (as shown in Fig Rock establishment algorithm considering rock shape and size is shown in Figure 7. (1) A random center point P 1 (x 1 , y 1 ) was created which was limited among known boundaries (slope contour). (2) Symmetrical random points P 2 and P 3 were created (shown in equations (1) and (2), resp.) such that lines P 2 −P 3 were the maximum diagonal of rock.
where d max is the maximum diagonal of rock; θ max is random angle corresponding to maximum diagonal ranging from 0 to π. (3) Minimum diameter point P 4 was created as follows: x 1 + 0.5d min cos θ min , y 1 + 0.5d min sin θ min , (3)

Advances in Civil Engineering
where d min is the minimum diagonal of rock (larger than 0) and defined as d min � μd max + k(2× rand − 1); μ is distortion coefficient (between 0 and 1); k represents the variation threshold value of distortion coefficient (between 0 and 0.5); θ min is random angle corresponding to minimum diagonal (between 0 and 2π) which is not equal to θ max . (4) Random rock points P 5 and P n were created as shown in equations (4) and (5).
x 1 + 0.5d n cos θ n , y 1 + 0.5d n sin θ n , where d 5 and d n are the random diameter of rock and are between d min and d max ; θ 5 and θ n are random angles corresponding to random diameter (between 0 and 2π); and θ n is not equal to θ n−1 , θ n−2 , . . . , θ 1 (n is 1, 2, 3, . . .). e algorithm of rock intrusion checking is described as follows (shown in Figure 6): (1) A random point such as point P B1 was established and was set as the center of circle and diameter was gradually increased until contacting with existing points (P 1 in rock R 1 ). (2) Rock and contact point were labeled as R 1 and P 1 , respectively, and other points were labeled as P 2 ·P 3 in a clockwise order. (3) e area S R1 of rock R 1 was calculated based on coordinates.
where the coordinates of P 1 , P 2 , and P 3 points were (x 1 , y 1 ), (x 2 , y 2 ), and (x 3 , y 3 ), respectively. (4) e area S R1 of rock R 1 considering the effect of point P B1 was also calculated as

Advances in Civil Engineering
where S 1 is the area of triangle P B1 -P 1 -P 2 , S 2 is the area of triangle P B1 -P 2 -P 3 , and S 3 is the area of triangle P B1 -P 3 -P 1 . (5) e area S R1 of rock R 1 was checked: if S R1 ′ was equal to S R1 , P B1 was in R 1; otherwise, it was out.
rough checking rock intrusions after establishing every random point, the established rock could interlock other rocks, as shown in Figure 6, which was more realistic than the rock distribution characteristic shown in Figure 5 (Liu et al.'s method [4]).

Establishment of Soil-Rock Slope.
Slope rock establishment process is shown in Figure 8. Firstly, corresponding   [4].  Advances in Civil Engineering rock boundary points were established according to rock characteristics such as graduation, content, and shape. Rock intrusion checking was conducted after each boundary point being established. Finally, corresponding soil-rock slope analysis model was developed and operation interface for soil-rock slope establishing software was shown in Figure 9.
For better understanding of the parameters in soil-rock slope establishing software (shown in Figure 9), a case study was performed on a soil-rock slope with 10% rock content and round rocks. Specific parameters for slope contour were "0, 0, 40, 0, 40, 15, 25, 15, 15, 30, 10, 30, 0, 35, 0, 0," which filled in "known boundary." "Mini-sides" and "Maxi-sides" were both set at 15, which controlled the side number of established polygons (15 polygons were defined as round); "distortion parameters" and "control distortion parameters" were set at 1 and 0, respectively, and "known boundary distribution" filled string "angle, average," which kept the length of established polygons constant; "relative diameter/ %" of "Mini-sides" and "Control-side1" were set at 0.199 and 0.201, respectively, to control rock size to be 0.2 Lc, and corresponding "Contents/%" was set at 10%. e established results of soil-rock slope with 10% rock content and round rocks were shown in Figure 10(a).

Typical Analytical Models.
Initial analytical models were obtained from engineering experiences of subgrades on mountainsides in western China. Original slope height was 42 m, which consisted of soil-rock mixtures. For building a subgrade on a mountainside, the top of the slope was excavated ( Figure 11) and a rigid subgrade layer was built with 5 m width and 0.5 m height. Final slope height and width were 15 m and 10 m, respectively, and the numerical analytical model is shown in Figure 11. In order to take into account the influence of rocks, combining the discussion in Section 2.2 and the precalculation results with different rock sizes, 2 m (0.2Lc, where Lc is the width of slope) was finally considered as the threshold for soil and rock in soil-rock slope model. By setting different parameters, such as rock sides and rock contents, 10 groups (12 models per group) of soil-rock slope models were established with different rock contents (10%, 20%, 30%, 40%, 50%, and 60%) and shapes (angular and round). Based on the Monte-Carlo stochastic method, ten groups (twelve models per group) stochastic models for soilrock slopes were established. One representative group was shown in Figure 10 where the slope model with high rock content was generated on the basis of low rock content for comparing the results with different rock contents.

Rock Distribution Characteristics in Soil-Rock Mixtures.
For better analysis of rock distribution characteristics in soil-rock slopes, 6 groups of soil-rock mixture samples were made consisting of different rock shapes (angular and round) and contents (10%, 30%, and 60%). e specimens were cylindrical in shape with diameter 50 mm and height 100 mm and each layer was compacted 20   [1,34,38], rock size in specimen was defined to be between 2 and 10 mm. CT [39] is a good method to observe the internal structures of objects without any disturbance and has been widely utilized in medicine and engineering. As shown in Figure 12, 6 section photos were obtained by CT and rock contours were described by black lines. It was seen in Figures 12(a) and 12(d) that some rocks were surrounded by soil and did not contact each other and rock-float skeleton was observed as shown in Figure 13(a). By increasing rock content, the probability of rocks contacting was also increasing, as shown in Figures 12(b) and 12(e). Some rocks contacted with each other and formed rock-part skeleton structures, as shown in Figure 13(b). When rock content reached 60%, they contacted with each other as shown in Figures 12(c) and 12(f ), forming rock-complete skeleton structures as shown in Figure 13(c). When rock content exceeded 60%, soil particles could not fill the voids between rocks and rock-void skeleton structure was formed, as shown in Figure 13(d), which is not much likely to appear in soil-rock slopes because slopes cannot keep their original shapes under this situation.

e Effect of Rock Distribution on Shear Failure Surface.
In order to better explain the shear failure surfaces of soil-rock slopes under top loading, unconfined triaxial tests under ring top loadings were conducted according to Wang's studies [4,6,38]. In their studies, the rationality of applying the obtained developing law of shear failure surface in soil-rock mixtures with small size to explain the shear failure surface of soil-rock slope with large scale had been proved. Unconfined triaxial tests were shown in Figure 14 and ring shaped loading was attached on the top of specimens prepared as described in Section 4.1. en, as loading was increased, the specimen was destroyed and a shear failure surface was formed. Based on failure surface, four typical failure appearances were concluded including rock loosen, rock no loosen, rock drop, and rock friction, as shown in Figure 14. en, five typical developing models of shear failure surfaces under top loading were determined which included "bypass, diversion, inclusion, penetration, and combination" modes, as shown in Table 1.   Advances in Civil Engineering slope but can also reflect its failure surface. However, it is not directly applicable to achieve the ultimate loading of slope. For obtaining ultimate loading on the mountainside shown in Figure 10, additional loading was applied to subgrade and the safety factor of slope was calculated based on strength reduce method. Ultimate loading was defined as the value when the obtained safety factor was between 0.99 and 1.01. A calculation example for ultimate load of slope is shown in Figure 15 where the geometric and material parameters of slope are shown in Figure 10 and Table 2, respectively. en, fixed constraint was set on the bottom of slope and horizontal constraint was set on left and right sides. e model was discretized using a triangular mesh, with maximum unit size of 0.02 Lc and the total number of meshed elements was 11239. As shown in Figure 15(a), the initial safety factor of slope was 1.18 which was then decreased with the increase of loading and the final ultimate loading value was obtained to be 56 kPa (safety factor was 1.00) as shown in Figure 15(b).

Analysis of Results considering Different Rock Contents.
Typical calculation results of one group (12 models) with FEM considering different rock contents and shapes are summarized in Table 3. For better studying the effect of rock distribution on slope ultimate loading, dimensionless results were utilized and the increasing multiple N of ultimate loading was put forward (the ratio between obtained ultimate loading value on slope with rocks and that on uniform slope).
As shown in Table 3, soil-rock slopes showed different ultimate loads under different rock contents and shapes, which can be explained by the developing modes of shear failure surface. As shown in Figure A, rocks were not distributed on shear failure surface; hence, ultimate load was still 56 kPa (N � 1.0) which was the same to uniform slope. As shown in Figure B, due to the distribution of rocks, shear failure surface showed "bypass" mode which caused shear failure surfaces to be not smooth; hence, the increasing multiple N increased to 2.6. As shown in Figure C, under the effect of "diversion" mode, the slope was hard to form throughout continuous shear failure surface; hence, increasing multiple N was 3.0. According to Figure D, under the effects of "bypass" and "penetration" modes shown in slope bottom and top, respectively, the developing location of shear failure surface was changed and antisliding force was increased; hence, increasing multiple N increased to 3.1.
As shown in Figure E, under the effect of "bypass" mode, the developing location of shear failure surface was seriously changed, which tended to appear as part failure at the top of slope; hence, ultimate load increased less (N � 1.5). As shown in Figure F, under the effects of "bypass," "diversion," and "inclusion" modes, the developing location of shear failure surface was changed, which tended to appear as part failure at the surface of slope and increasing multiple N was 4.0. As shown in Figure G, under the effects of "bypass," "diversion," and "inclusion" modes, the slope was hard to form smooth main shear failure surface and increasing multiple N increased to 6.1. As shown in Figure H, under the effect of "bypass," "diversion," and "combination" modes, the slope was hard to form throughout main shear failure surface and increasing multiple N increased to 5.6.
As shown in Figure I, under the effects of "penetration" and "combination" modes, the developing location of shear failure surface was changed, which tended to appear as part failure at the surface of slope. At the same time, shear failure surface between rocks showed greater antisliding force and Small rocks integrate into larger rocks and shear failure surface is seriously affected therefore, increasing multiple N increased to 9.8. As shown in Figure J, under the effects of "bypass" and "inclusion" modes, the developing location of shear failure surface was changed, which tended to appear as part failure at the surface of slope and increasing multiple N increased to 10.4. As shown in Figure K, under the effects of "bypass," "diversion," and "inclusion" modes, the slope was hard to form throughout and smooth main shear failure surface and increasing multiple N was 5.9. As shown in Figure L, under the effects of "bypass" and "inclusion" modes, the developing location of shear failure surface was changed, which tended to appear as part failure at the surface of slope and increasing multiple N increased to 10.6. Combining the analytical results of Figure 15 and Tables 2 and 3, it was concluded that the influence of developing modes of shear failure surface on slope included ① adjusting the developing location of shear failure, which tends to cause different failure modes of slope, such as part failure shown in Figure E; ② changing the shape of shear failure surface, which makes shear failure surface be not throughout, continuous, and smooth, such as the failure surface shown in Figure H; ③ improving the antisliding force of slope, such as the failure surface shown in Figure I.
For better estimating the influence of rock contents on the ultimate loads of slope and reducing the impact of the specificity distribution of rocks on analytical results, statistical results between increasing multiple N and rock content were obtained as shown in Figure 16 and Table 4.
As shown in Figure 16, the analytical results obtained from FEM revealed that increasing multiple N of ultimate loading was increased with the increase of rock content but with different amplitudes. Comparing the differences in increasing amplitudes, two sections were obtained including rapid (0%-30%) and slow (30%-60%) increasing periods. In rapid increasing period, the value of N reached 17.7 times and changing amplitude was 16.7, while in slow increasing period, N value reached 23.65 times and changing amplitude was 5.95. By comparing discrete regions between different rock contents, it was observed that, in total, discrete regions were gradually increased and corresponding results are shown in Table 4.

Mechanical Analysis considering Different Rock Contents.
In analytical models, the volume of each rock was set to be the same; therefore, the amount of rock was increased with the increase of rock content. By increasing rock content from 0% to 60%, more rocks appeared in slopes, which increased the amount of rocks distributed in shear failure surface. Hence, with the increase of rock content, failure surface was more seriously affected by rocks including adjusting the developing location of shear failure, changing the shape of shear failure surface, and improving the antisliding force of slope. e above mentioned influences improved slope stability and increased ultimate loading, and the amount of improvement was increased with the increase of rock amount due to rock combinations. At the same time, with the increase of rock content, rock distribution was more diversified in the premise of random algorithm. Hence, the average value of N and discrete region were both increased by increasing rock content.
Different increasing amplitudes of the average value of N could be explained by rock skeleton status. When rock content exceeded 30%, soil-rock mixtures were in rockcomplete skeleton status. Under this condition, based on experimental results described in section 4.1, rocks fully filled   10 Advances in Civil Engineering slopes and increasing rock content played more important role in changing relative distances between rocks that had less improving effect on soil-rock slope stability. When relative distance between rocks was small, they could be recognized as one combined rock (Table 2), which reduced the amount of rocks in slopes and led to small increases in N. Hence, less increasing amplitude of N appeared and discrete region was slowly increased when rock content was above 30%.

Analytical Results considering Different Rock Shapes.
As shown in Figure 17 and Table 5, slopes containing angular rocks had greater average ultimate loading multiple N and discrete region than those containing round rocks. In FEM results, the maximum average value of N (60% rock content) obtained for slopes with angular rocks was 25.6, which was about 2 times greater than 11.9 obtained for slopes with round rocks. Discrete regions in slopes with angular rocks were at least 3 times greater than those with round rocks when rock content was in the range of 10% to 60%.

Mechanical Analyses considering Different Rock Shapes.
By comparing results obtained for angular and round rocks, it was concluded that the main difference was in rock shape. Different from round rocks which had streamlined shapes, angular rocks had obvious sharp edges and irregular shapes. Based on Liu et al.'s study [4] and experimental findings in Table 1, it was concluded that angular rocks tended to guide failure surface development along edge tangent. In other        words, failure surface was easier to be affected by angular rocks than round rocks, which greatly improved slope stability under top loading. Hence, soil-rock slopes with angular rocks had higher average N values than those with round rocks. In addition, in terms of random algorithm, angular rocks had higher distribution characteristics than round rocks due to geometric properties such as rock inner angles, the angle between the long axis of rock and slope, and number of sides. It means that higher distribution characteristics and more effective models of rocks for shear failure surface appeared in soil-rock slopes with angular rocks leading to higher discrete regions.

Stability Analytical Results of Soil-Rock
Slopes under Topping Loading 6 Figure 18). Meanwhile, lower boundary was also determined which could be considered as reference in estimating soil-rock slope stability.

ree Typical Failure
Modes. Based on FEM calculation results, three typical failure surfaces were concluded (Figure 19), including deep failure surface (surface-1), top partial failure surface (surface-2), and shallow failure surface (surface-3). Surface-1 is a common slope failure surface, especially in soil-rock slopes with low rock contents. At lower rock contents, the influence of rock on shear failure surface is mainly increasing the length of shear failure surface and antishear force of slope, which do not change the developing tendency of shear failure surfaces. Surface-2 usually happens when a single rock or a pile of rocks is distributed in the middle of slope preventing the development of shear failure surface deep in the slope resulting in the appearance of shear failure surfaces toward the surface of slopes.
Surface-3 usually happens when a single rock or a pile of rocks is distributed under subgrade. Due to the accumulation of rocks, the soil under the subgrade is strengthened and therefore, shear failure surface appears on the surface of slope. In this failure model, slope has greater ultimate loading than other models.
In engineering, if the slope needs to be treated, it is advised for designers to consider the different failure modes of soil-rock slope for obtaining the best treated effect.

Conclusion
In this study, for obtaining better estimation of limited loadings on slopes under topping loading, a new random rock contour establishment algorithm was introduced and corresponding soil-rock slope models were developed considering different rock contents and shapes. en, based on experimental findings, four rock skeleton models and five shear failure modes were introduced. Finally, the influence of rock distribution on the ultimate loading of slopes was discussed through FEM and three typical failure modes were obtained in soil-rock slopes.
(1) Traditional random rock contour establishment algorithm cannot take into account interlock effect between rocks. Hence, a new rock counter establishment algorithm was proposed, which was able to take into account the interlock effect of rocks through first establishing random rock boundary points and then checking for rock intrusions. (2) According to CT section images in soil-rock mixtures with different rock contents and shapes, four typical rock skeleton status were developed based on rock distribution characteristics in specimens, including rock-float, rock-part, and rock-complete skeleton status. en, unconfined triaxial tests with ring top loadings were performed on specimens. Based on the developing paths of shear failure surfaces and rock distribution characteristics, five shear failure surface developing modes, including bypass, diversion, inclusion, penetration, and combination modes, were put forward and corresponding appearances in tests were given.
(3) By the proposed rock contour establishment algorithms, ten groups (twelve models per group) were established and ultimate loading on slope top considering different rock contents and shapes was calculated using FEM. en, normalized increasing multiple N was utilized to describe the influence of rocks on slope ultimate loading. It was concluded that ultimate loading was nonlinearly increased with increasing rock content, which could be explained in terms of rock skeleton status. Soil-rock slopes with angular rocks had about three times higher ultimate loading than those with round rocks, which could be explained by more complex distribution characteristics of angular rocks. Finally, it was shown that assuming soil-rock slopes as uniform slopes underestimated ultimate loading on slope top and rock contents and shapes had to be considered. (4) In order to improve the applicable range of the calculated results in this paper, more studies were conducted and the reference values of ultimate loading increasing multiple N of soil-rock slopes under top loading were suggested considering that different slope angles (50°-60), rock contents (0%-60%), and rock shapes (round and angular) were obtained. (5) Based on shear failure surfaces of soil-rock slopes under topping loading, three typical failure surfaces, including deep failure surface, top partial failure surface, and shallow failure surface, were developed which could be reference for designers to treat slopes.
In this paper, through the combination of stochastic modeling and finite element model (FEM), the influence of rock contents and shapes on the ultimate load of slope is studied. In order to better serve the engineering about soil-rock slope, the effect of soil-rock contacting surface and cracks on the slope will be studied and general particle dynamics (GPD) and extended finite element method (XFEM) will be used to consider material discontinuity in the next research.

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

Conflicts of Interest
e authors of this article confirm that the received funding that has been mentioned did not lead to any conflicts of interest regarding the publication of this manuscript.