Stability evaluation and potential failure process of rock slopes characterized by non-persistent fractures

.

and structural fractures). Discontinuities in rock slopes were formerly assumed in slope analysis as through-going and fully persistent planes. However, many researchers have suggested that discontinuities are generally not fully persistent (Eberhardt et al., 2004;Scavia, 1995;Terzaghi, 1962) and that failure surface is a combination of preexisting non-persistent discontinuities and newly propagated cracks (Brideau et al., 2009;Einstein et al., 1983;Frayssines and Hantz, 2006;Lajtai, 1969;Zhang et al., 2017). Non-persistent fractures are highly complex and come in various sizes and properties, and their 35 locations and characteristics are difficult to determine. Thus, the influence of non-persistent fractures on slope analysis poses a great challenge (Fan et al., 2015;Wasantha et al., 2014).
However, these non-persistent fractures were artificially specified by researchers as having the same sizes and orderly 40 arranged locations. The complexity and random distribution of non-persistent fractures in rock slopes should be considered in reflecting and modeling the stability and failure process of real rock slopes in nature.
According to a major research work on fractured rock masses, the discrete fracture network (DFN) modeling technique maximizes the use of discontinuity data from exposed surfaces and can become the best option for simulating realistic fractured rock masses (Bonilla-Sierra et al., 2015;Chen, 2001;Elmo et al., 2013;Pine et al., 2006). By coupling the DFN 45 technique with continuum, discontinuum, and hybrid modeling approaches, synthetic rock mass (SRM) models can be set up to investigate the mechanical properties of fractured rock masses (Elmo and Stead, 2009;Mas Ivars, 2010;Pierce et al., 2007). SRM approach has been widely used to determine the representative elementary volume size of a fractured rock mass (Esmaieli et al., 2010), reproduce rock mass properties and behaviours (Mas Ivars et al., 2011), and simulate fracture propagation in a fractured rock mass (Zhang and Stead, 2014). SRM models have been primarily used to simulate failure and 50 deformation of fractured rock slopes (Bonilla-Sierra, 2015;Elmo and Stead, 2013). DFN simulation included in SRM modeling program exists a high degree of variability. The variability of DFN simulation means there are an infinite number of possible realizations of 2D fracture systems exist given specified input parameters (Pine et al., 2006). The results of SRM model analyses that incorporate different DFN models may vary significantly (Elmo et al., 2018;Zhang et al., 2012). Therefore, rock slope analyses based on only one DFN model can lead to erroneous results (Elmouttie and Poropat, 2014;55 Mas Ivars, 2006). A statistical analysis based on a large number of DFN models may reduce the aforementioned errors and provide reasonable results for rock slope analysis and support (Ferrero et al., 2016;Zhang et al., 2013;Zhang et al., 2017).
The current study proposes a comprehensive approach that combines several well-established methods to conduct a stability evaluation and failure process analysis of a fractured rock slope in Tianjin City. First, 100 DFN models are generated on the basis of the fractures collected in the field. Second, slope models are constructed using SRM approach. 60 Third, the improved gravity increase method is employed to determine the stability of a single-slope model. Fourth, the potential failure processes of the fractured rock slope models are simulated. Fifth, the final critical slip surface, factor of safety, and accumulation distance are determined on the basis of the statistical analysis of 100 slope models generated with different DFN models. https://doi.org/10.5194/nhess-2020-58 Preprint. Discussion started: 30 March 2020 c Author(s) 2020. CC BY 4.0 License.

4
The mean value and probability density function (PDF) of the corrected trace lengths for each fracture set can be determined according to the methods proposed by Kulatilake and Wu (1984) (Table 1).
A 2D analysis was performed in this study. The cross section normal to the exposed surface was used to perform the 2D stability analysis and potential failure process of the rock slope. A fracture plane intersected by the cross section with a small angle acted as the surface of separation, which will not influence the results of stability analysis. Therefore, the fractures 100 intersected by the cross section with an angle smaller than 20° were artificially deleted prior to the stability analysis. We derived the characteristics of the fractures (location, orientation, trace length, and density) in the cross section from those features measured from the exposed surface. Subsequently, the 2D fracture traces were generated using Monte Carlo simulation by synthesizing the aforementioned characteristics. The specific processes are clearly described in the research of Zhang et al. (2017) and are briefly introduced as follows. The locations of the fractures were assumed to follow Poisson's 105 distribution. P 21 (the 2D fracture density) in the cross section was proved to be tanη i times that of the exposed surface (where η i is the acute angle between the mean normal vector direction of fracture set i and the NS direction). The 2D orientation of a fracture is reflected by its trace gradient k, which can be expressed as sinα/cotβ (where α and β are the dip direction and dip angle of a 3D fracture, respectively). Empirical distribution was followed to generate the 3D fracture orientations based on the fracture orientation frequency, and then k in the cross section can be determined. The results of Zhang et al. (2017) 110 revealed that the mean [E (ch 2 )] and variance [V (ch 2 )] values of the trace length square are constant. E (ch 2 ) and V (ch 2 ) of the collected fracture traces in the cross section can be determined according to those features on the exposed surface. The rooted numbers result in the trace lengths in the cross section. The cross section of the investigated rock slope was 20 m high, and the smallest mean dip angle of the three fracture sets was 31.7°. Thus, the length of the cross section was 30 m [20 × cot (31.7°)]. Finally, we used Monte Carlo simulation to merge the aforementioned parameters to generate the 2D DFN model in 115 a 30 m × 20 m cross section.
More than one DFN model can be generated with the same fracture data. For example, Fig. 4 exhibits four DFN models with different fracture characteristics. Therefore, the rock slope analyses based on only one DFN model can lead to wrong results (Mas Ivars, 2006;Xu et al., 2014). To solve this problem, we generated and verified the validity of numerous DFN models by applying the following method. Initially, we intersected the fractures in the exposed surface with lines extending 120 along the dip direction. Doing so produced a series of intersection points for an individual line. Subsequently, we intersected the DFN models using the exposed surface and consequently generated a series of intersection points for an individual network. We compared these two sets of intersection points by using their probability density curves. The DFN models generated in the cross section proved reasonable when the results were identical to one another. Finally, we selected 100 reasonable DFN models to construct different slope models and conduct stability analysis and potential failure process 125 simulation.

SRM model for rock slope analysis
SRM approach is used to construct fractured rock slope models. In the present study, the SRM model used was a combination of two well-established models, namely, the bonded-particle model (BPM) in PFC2D and the DFN model (details are available in the research of Pierce et al. (2007). SRM approach is widely used to simulate the mechanical 130 properties of fractured rock masses, perform rock masses stability analysis, and simulate the failure and deformation of fractured rock slopes.

Parameter determination for SRM model
The SRM model in PFC2D is defined by many parameters that cannot be directly identified via laboratory and field experiments. Therefore, the parameters of the model should be predetermined according to the macroscopic characteristics of 135 rock slopes. In particular, several numerical tests should be performed to ascertain and quantify the input parameters of the BPM and DFN model. In PFC2D, a material is considered an assemblage of bonded rigid circular particles, and particle size distribution dramatically influences modeling behaviour (Mas Ivars et al., 2011). Therefore, determining the particle size of a numerical specimen remains a challenge. Generally, the ratio between the maximum and minimum radii of particles is 1.66 (Potyondy and Cundall, 2004). A small particle size is indicative of accurate simulation results. However, the number of 140 particles of rock models increases with a decrease in particle size. A large amount of particles will result in a long calculation time and low computational efficiency. In the present study, particles with radii between 0.05 and 0.083 m were finally selected to fill rock models after repeatedly changing the particle size. In that case, the size of particles is small enough for the comparatively reasonable results and the calculation time is acceptable.
The BPM parameters are calibrated against laboratory uniaxial and biaxial compression tests. The specimens (the 145 height-to-diameter ratio is 2:1) for the numerical uniaxial and biaxial compression tests are set up to reproduce the macroproperties of real rock materials, such as uniaxial compressive strength (σ c ), elastic modulus (E), Poisson's ratio, friction angle, and cohesion obtained in laboratory tests. Different macro-properties are influenced by different parameters.
Specifically, E is mainly controlled by several parameters, including particle contact modulus (E c ), particle normal/shear stiffness ratio (k n /k s ), parallel bond modulus (E b ), and parallel bond normal/shear stiffness ratio (k nb /k sb ); Poisson's ratio is 150 influenced by k n /k s and k nb /k sb ; σ c is influenced by normal strength, internal friction angle and cohesion of parallel bonds; and cohesion and friction angles are influenced by the friction coefficient (Bahaaddini et al., 2013). The values of the aforementioned parameters are empirically assigned in advance, and then the numerical uniaxial and biaxial compression tests are carried out. When the macro-properties of the numerical tests correspond to the results from the laboratory tests, the parameters are considered reasonable. Otherwise, the parameters are adjusted until the rock specimens have the same macro-155 properties as real rock materials (Park and Song, 2009;Yang et al., 2006). The calibrated parameters of particles and bonds are listed in Table 2. The macro-properties of the rock specimens for the numerical tests and those of the real rock materials for the laboratory tests are listed in Table 3.

6
In PFC2D, a smooth joint (SJ) model is commonly adopted to generate fractures. In this model, two particles that lie on the opposite sides of the intended fracture plane can overlap and pass through each other instead of moving along their 160 perimeters, thereby reproducing the real physical and mechanical properties of fractures (Bahaaddini et al., 2013). The SJ parameters are normal stiffness (k nj ), shear stiffness (k sj ), and coefficient of friction (μ j ). To determine these parameters, we built a numerical specimen (with a width-to-height ratio of 1:1) for the direct shear tests and normal deformability tests to simulate macro-properties, including shear stiffness, normal stiffness, and friction angle. The SJ parameters can be obtained as long as the test results approximate those of the laboratory tests. Specifically, k nj was obtained by the numerical normal 165 deformability tests; k sj and μ j were determined by the numerical direct shear test. Table 4 lists the values of the SJ parameters and the results of the normal deformability and direct shear tests of the numerical and laboratory tests. The specific calibration procedures are complex and are thus not introduced in this paper in detail. Additional details about the calibration can be found in the works of Bahaaddini et al. (2013), Cheung et al. (2013), and Duan et al. (2016).

Model generation 170
The investigated rock slope with a relatively low height is located on a ground surface with a relatively flat terrain. The study area is dry, crustal movement is not obvious, and no active fault exists nearby (Sect. 2.1). Thus, the effects of ground stress, water, and earthquake on the slope analysis were not considered in the current work.
The size of the slope section is 20 m (height) × 30 m (length). The bottom and right sides of the slope section were expanded by 10 m as the boundary section (Fan et al., 2004). Ultimately, the size of the SRM model was determined to be 30 175 m (height) × 40 m (length) (Fig. 5). The upper boundary of the model was free. Moreover, the left, right, and bottom boundaries were assumed to be smooth rigid walls.
The particles with the same radii (0.05 to 0.083 m) as those reported in Sect. 3.1 were applied to fill the 20 m × 30 m slope section. Considering the small effect of boundary section on the slope analysis, we filled the boundary section with particles with larger radii (0.1 to 0.15 m) to improve computational efficiency (Fig. 5). A total of 48,947 particles were 180 generated, and the parameters presented in Table 2 were adopted for the bonds and rock particles.
Gravity was applied to the model, and the model was calculated (cycled) until the particle assemblage reached an equilibrium state (i.e., the unbalanced forces reached the required standard of 10 −6 ). Then, the generated DFN model reported in Sect. 2.2 was added to the slope section by defining the FISH functions, thereby forming the SRM model (Fig. 5, with the DFN model in Fig. 4a as an example). The SJ parameters presented in Table 4 were adopted for the fractures. This process 185 ignored the stress concentration at the tips of the structural fractures generated by tectonic stress. Nevertheless, the process was considered reasonable in this study because the stress concentration was intensely reduced after the long-term stability of the rock slope.
The slope was formed within one operation. Therefore, the left boundary (smooth rigid wall) was removed to simulate a one-time excavation of the slope. To trigger instability, an improved gravity increase method is used in this study. The improved gravity increase method was proposed by Meng et al. (2015). This method leads to the failure of a slope in PFC2D by slowly increasing gravity acceleration and reducing the friction coefficient of particles while keeping other parameters constant. The amplitude of 195 reduction of the friction coefficient is the same as that of the increase in gravity acceleration to ensure an approximate invariance of the resisting force. When the slope model is in a limit equilibrium state, i.e., fractures in the slope are propagated and coalesced until a through-going slip surface (i.e., the critical slip surface) is initially formed, the factor of safety F can be defined as the ratio of the gravity acceleration in the limit equilibrium state (g') to that in the initial state (g), i.e., F= g' / g. Taking the model in Fig. 5 as an example, we calculate the factor of safety by using the calculation procedure 200 shown in Fig. 6. The simulation results indicate that the factor of safety of the slope model using the improved gravity increase method is 25. Additional details on the factor of safety can be found in Sect. 7 of this paper.

Initiation and propagation of fractures
When the factor of safety is determined, the critical slip surface is simultaneously obtained. To get an in-depth understanding of the fracture propagation mechanism, the propagation process of fractures and the evolution of force chains during the 205 formation of the critical slip surface are recorded (Fig. 7). The force magnitude is proportional to the thickness of the line in the force chain plots (lower right corner of Fig. 7), in which the blue and green lines denote compressive and tensile forces, respectively. The color of the line segments is obvious in the region where the stress concentration is strong.
Non-homogeneous stresses are distributed throughout the slope owing to the heterogeneity of the slope model. After 2,000 time steps, the compressive stress, which slowly increases from top to bottom under the action of gravity, is initially 210 distributed throughout the slope. The tensile stress exists only in the tips of the fractures (Fig. 7a). Then, the degree of tensile stress concentration increases at the fracture endpoints (Fig. 7b). The fractures (black lines surrounded by a pink circle in  Fig. 7 corresponds to the fracture propagation direction. Finally, a decrease in forces, especially the compressive forces (lighter color) throughout the slope, is observed; furthermore, a through-going surface, i.e., the critical slip surface, is formed when the propagated fractures arrive one after another to the tip of the original fractures (Fig. 7d). 220 https://doi.org/10.5194/nhess-2020-58 Preprint. Discussion started: 30 March 2020 c Author(s) 2020. CC BY 4.0 License.

5 Potential failure and accumulation process
According to the aforementioned results of stability analysis, the factor of safety is extremely high such that the investigated rock slope is highly stable. Nonetheless, the potential failure and the accumulation process are simulated in the present study to obtain good knowledge of the failure mechanism of fractured rock slopes and provide references for similar slope projects.
To maintain coherence of analysis, we use the model in Fig. 5 as an example. Section 4 presents an analysis of the formation 225 process of the critical slip surface of the gravity-increased model. Subsequently, another round of analysis is performed to determine the failure process after forming the critical slip surface. indicates the occurrence of a small deformation in the slip mass above the critical slip surface; the bedrock remains stable without a distinguished gap with the slip mass. The displacements of the slip mass continuously increase, and the largest displacement of the particles (red particles in Fig. 8b) is nearly 0.2 m, which indicates the aggravation of rock failure (Fig.   8b). After 60000 time steps, the slip mass is slowly fractured into many rock blocks along the non-persistent fractures (Fig.   8c). The displacement and size of these blocks vary from one another. The block fragmentation is apparent near the critical 235 slip surface because of the newly propagated fractures, whereas most of the rock blocks far from the critical slip surface remain intact (Fig. 8d). Ultimately, the slip mass is completely separated from the bedrock (Fig. 8d).
As soon as the slip mass is detached from the bedrock, a large displacement and movement of blocks occurs. To reflect the actual failure and accumulation process of rock slopes in nature, the PFC2D procedure should be manually interrupted to make corresponding adjustments (when the time steps are 80,000, which correspond to Fig. 8d). First, the particles in the 240 bedrock are deleted to improve computational efficiency, and the critical slip surface is replaced with a rough rigid wall. The frictional properties of the rigid wall are the same as those of the propagated fractures. Then, the DFN model is removed, and the body force in the model is initialized to avoid splitting caused by the release of stress. Finally, the gravity acceleration is restored to 10.0, which corresponds to natural conditions. Figure 9 presents the computed results. In the figure, the blue sections denote the assemblage of rock particles from a macroscopic level while the red lines represent the bonds (contacts) 245 between particles. For a convenient description, we numbered several rock blocks from 1 to 6.
In Fig. 9a (5×10 4 time steps), except for block 2, the blocks rotate under inertia force and gravity and are separated from one another. The blocks near the slip surface are disintegrated into numerous sub-blocks, and even crushed particles. Rock blocks 1, 3, 4, and 5 rotate counterclockwise by approximately 60°, 30°, 260°, and 40°, respectively; whereas block 6 rotates clockwise by 75° (Fig. 9b) when the time steps are 2×10 5 . Block 6 initially crashes to the ground, and the collision results in 250 the bond breakage between particle blocks (red segment in Fig. 9c). Block 6 is divided into two sub-blocks and some crushed particles. Blocks 4, 5, and 2 successively crash to the ground one after another but their shapes are kept intact, whereas block 3 is split owing to the collision with block 6, as shown in Figs. 9c and 9d. As soon as all crushed particles and https://doi.org/10.5194/nhess-2020-58 Preprint. Discussion started: 30 March 2020 c Author(s) 2020. CC BY 4.0 License. blocks fall to the ground or the critical slip surface, the blocks and particles slide or roll forward as a whole and are accompanied by a fragmentation of the block edges (Fig. 9e). Figure 9f presents the final shape of the rock model. The final 255 deposit is composed of relatively intact rock blocks and crushed particles, and the blocks pile up above the crushed particles, presenting an inverse grading phenomenon. We also record the final accumulation distance, which is represented by the farthest reach distance of intact blocks (bonds exist in particles in these blocks). The final accumulation distance of the rock slope model shown in Fig. 5 is 28 m.

Stability analysis 260
The slope analysis in Sects. 4 and 5 are based on a SRM model with one DFN model. However, as mentioned in Sect. 2, numerous DFN models can be generated because of the variability of DFN simulation, and the stability analysis of only one SRM model may lead to erroneous results. In the present study, 100 SRM models on the basis of the different DFN models (generated in Sect. 2), are built to conduct the aforementioned analysis. In particular, the critical slip surfaces and safety factors of these models are calculated using the improved gravity increase method following the calculation procedures in 265 Figure 6. Then, the method mentioned in Sect. 5 is employed to simulate the potential failure and accumulation process. 12a and 12b present the factors of safety of the 100 models, ranging from 25 to 73.5. Therefore, the variability of the simulation results should be emphasized in the stability evaluation of fractured rock slopes. The outcome that considers numerous model calculations may lead to a rational result. In the present study, a statistical analysis method is applied to solve this problem. The final critical slip surface and factor of safety are attributed to conservative considerations. The 275 potential critical slip surface constitutes a large slip mass (Fig. 11). The final critical slip surface covers over 90 % of the slip mass to guarantee safety in the rockslide analysis and support. Additional information on the critical slip surface can be found in the research of Zhang et al. (2017). Critical slip surfaces are supposed to have an arc-shaped geometry that differs from their actual linear or broken line morphologies. Arc morphology is easily defined and is convenient for practical engineering designs. F s is the final factor of safety of the rock slope when the F values of the other 90 models are greater 280 than F s The result shows that the final factor of safety of the rock slope are 43.5.
Similarly, the results of the potential failure and accumulation process vary. Figure 13 presents the final accumulation states of the models shown in Fig. 10. The plots indicate that the sizes and quantities of the fractured blocks, as well as the rotating degrees, are significantly different. The final accumulation distances of the models are 28, 34, 35, and 40 m.  (Fig. 12c). The potential failure and accumulation process based on one SRM model may obtain erroneous results. A statistical analysis method is also applied to solve this problem. To guarantee safety in the rock slope analysis, we attribute the final accumulation distance to conservative considerations. In particular, when the distance values of the other 90 models are lower than a certain value, then the value is the final accumulation distance (D a ) of the rockslide. Therefore, the final accumulation distance is 87 m. 290

Discussion
Slope failure is a 3D stability problem. Thus, constructing a 3D SRM model for stability evaluation and failure process analysis is highly convincing and accurate. A 3D SRM model also combines rock masses in the BPM with DFN model. However, the quantity of rock particles in a 3D slope is extremely large to conduct an effective calculation in PFC. In addition, the factor of safety obtained by 3D analysis is generally higher than that in 2D analysis. Given that many theories 295 and technologies cannot be established and that their adoption cannot perfectly reflect rock masses at present, the accuracy of analyses may not satisfy engineering project requirements. Moreover, deriving a high factor of safety is sometimes unfavourable. Accordingly, 2D analysis, which is simple and commonly used in reality, is adopted in the present study for the stability evaluation and potential failure analysis of the investigated slope. The 2D analytical result can be regarded a good reference, and it provides a theoretical and practical basis for future initiatives that utilize 3D analysis. 300 The factor of safety of the investigated slope is extremely high. In the field investigation, unlike that for various nonpersistent fractures, weak interlayer and through-going discontinuities are not observed; thus, the non-persistent fractures play a vitally important role in the stability of the investigated slope. The strength of the rock bridges (intact rock) is considerably higher than that of the fractures. Therefore, the obtained factor of safety is extremely high and is thus reasonable. In addition, the effects of water (rainfall) and earthquakes were ignored in the present study. However, the 305 accuracy of the calculation result would increase if rainfall (seepage analysis) and earthquakes (kinetic analysis) are considered. This topic will be the direction of our future research.
The failure process is unlikely to occur in the investigated rock slope unless it is subjected to significant environmental changes, such as earthquakes, rainfall, unloading, or overloading. Nonetheless, the potential failure process is simulated in this study to understand the rockslide mechanism and subsequently provide a good reference for similar slope projects. For 310 example, the size and motion of rock blocks can be utilized to predict risk degree, and the final accumulation of deposits can contribute to hazardous area division. The final arrangement of deposits (a combination of blocks and crushed particles) provides a good explanation for the inverse grading of rock avalanche reported by Cruden and Hungr, (1986), Imre et al. provide conservative results for slope support. Finally, although the final results of the factor of safety, critical slip surface, and accumulation distance can guarantee safety in the rockslide analysis and support, they are not the exact values. Statistical analysis provides a new method for deriving an in-depth understanding of solid earth, where specific locations and characteristics of geological materials and structures remain unknown, such as discontinuities, especially for small-scale 320 non-persistent fractures. Meanwhile, new theories and technologies are required to obtain precise forecasts with respect to the range values characterized by statistical methods.

conclusion
The present study combines several methods, namely, DFN simulation, SRM approach, and statistical analysis, to conduct stability evaluation and potential failure process analysis of a fractured rock slope in Laohuding Quarry in Jixian County, 325 Tianjin. The SRM technique is utilized to generate a slope model with non-persistent fractures in the form of a DFN. An innovative formula for the safety factor is deduced by considering stress concentration on the basis of the improved gravity increase and strength reduction methods. The formation of a critical slip surface is also investigated. The potential failure and accumulation processes are simulated to provide a reference for similar slope projects. Numerous slope models are calculated, and the final results of the safety factor, critical slip surface, and accumulation distance are determined by 330 statistical analysis. The major findings are summarized as follows.
(1) The slope model with non-persistent fractures can be effectively constructed on the basis of SRM technology. The instability of the slope model can be attained by combining the improved gravity increase method or the strength reduction method, thereby obtaining the safety factor and critical slip surface. An innovative formula to calculate the safety factor is proposed by considering stress concentration and the calculation principle of PFC2D. 335 (2) Fracture propagation is closely related to stress concentration. Fractures initially propagate from the tips of the original fractures where the tensile stress is concentrated. Then, the stress is released, and a new stress concentration occurs at the tip of the propagated fractures when the fracture propagates downward to the neighbouring fractures. The critical slip surface is formed by the coalescence of preexisting fractures and newly propagated fractures.
(3) In the initiation of failure, the slip mass is fractured into rock blocks along the preexisting fractures. Then, most 340 blocks rotate and collapse under inertia force and gravity. Several blocks are split into sub-blocks owing to the collision between the blocks and the ground. The final deposit is composed of intact blocks and crushed particles, presenting inverse grading phenomena.
(4) The critical slip surfaces, factors of safety, and accumulation distances of the slope models with different DFN models vary. Therefore, the final outcome is obtained by statistical analysis. It ensures engineering safety for rockslide 345 analysis and support. The factor of safety (reserve) of the studied rock slope is determined to be 43.5. The critical slip surface is confirmed, as shown in Fig. 11. The final accumulation distance is 87 m.