Some considerations on the role of joint spacing in rock slope stability analyses using a 3D discrete element approach

Stability assessment of potential wedge failures in rock slopes is usually carried out based on stereographic projection techniques and limit equilibrium analyses. Nevertheless, this methodology presents several limitations (i.e., lack of information on joint mechanical properties, oversimplification of the structure, etc.) that could lead to poorly accurate conclusions. In line with this idea, the main objective of this work is to evaluate the effect of joint spacing in the estimation of the factor of safety of rock slopes affected by potential wedge failures. For this purpose, a 3D numerical distinct element code (3DEC) was selected to carry out a substantial number of simulations in which the factor of safety of two types of slopes, affected by different discontinuity sets, was studied by using the shear strength reduction (SSR) method. Different values of joint spacing, cohesion, and friction angles were considered, combined with two angles of the slope face under study. A discrete fracture network (DFN) model was also implemented, with the aim of benchmarking results with those obtained from the original models. The joint spacing has been found to relevantly affect the values of the factor of safety of the slope, which showed variations of up to 40% in comparison with those obtained from limit equilibrium methods. This work provides an insight to a more realistic interpretation of rock slope analyses against wedge failures, particularly to better estimations of the factor of safety.


Introduction
Rock-mass excavation is a common practice in the fields of civil and mining engineering when developing a variety of structures like foundations, tunnels, or slopes. All these activities necessarily involve an interaction with the rock mass, which originally presents characteristics that are determined by its formation and geological history, and that will be disturbed when the excavation is being developed.

3
A well-known feature of any rock mass is that its mechanical behavior and, hence, the stability of any structure placed in or on it is greatly determined by the presence and characteristics of discontinuities. According to Priest (1993), the discontinuity spacing can be considered as one of the most relevant indicators when determining the so-called quality of a rock mass, seeming reasonable to think that the magnitude of this parameter will be related to the general stability of slopes. This fact can be illustrated by the photographs shown in Fig. 1, in which two rock masses, affected by several families of discontinuities but with very different spacings, are shown.
A preliminary stability analysis of the slopes like those shown in Fig. 1 would often require plotting the main joint sets by means of the stereographic projection, in such a way that potential instability mechanisms can be identified and kinematic feasibility of failures determined (Wyllie 2018). This can be facilitated by resorting to commercial software like Dips (Rocscience 2020a), even though other freeware options like Stereographic Projection (Sakellariou and Kozanis 2006), InnStereo (Schönberg et al. 2018), or Stereonet (Allmendinger 2020) are currently available. After that, a widely applied analytical technique like the limit equilibrium method (LEM) is often used for rock-slope design (Stead et al. 2001;Wyllie 2018), through the estimation of the socalled factor of safety (FS). This FS is expressed in terms of the ratio of the resisting forces or moments, to the driving ones. The deformability of materials is not accounted for in this approach. The stability threshold is defined by FS = 1. In a similar way as presented for the stereographic projection technique, LEM is nowadays implemented in programs like SWedge (Rocscience 2020b) or RockStability (GEO5 2021).
Even though remarkable advances were made to offer user-friendly interfaces, simplicity, and good graphical outputs on the already referred programs, they still present inherent limitations. In the case of LEM, the oversimplification of complex phenomena and the non-consideration of strain/displacement compatibility or, more particularly, deficiencies in the case of wedge stability analyses like a correct definition of the shear forces on the slip surfaces, as pointed out by Deng (2021). The stereographic-based solutions solely consider the dip and dip direction of the principal joint sets and slope for generating the wedge geometry to be analyzed, disregarding the effect of spacing, as illustrated in Fig. 2, where two slopes affected by two families of discontinuities (J1 and J2) are represented. As can be seen, the dip and dip directions of the joint sets are the same in both examples, which is evident when analyzing their representation by stereographic projection.
In this line, the existence of smaller, probably less-stable wedges is disregarded and so, the overall safety factor of the slope. Therefore, the main goal of the present study was to assess the effect of joint spacing on the safetyfactor values of rock slopes, highlighting the importance of considering this feature in order to obtain more reliable stability analyses.

The wedge failure: a background
Wedge failure can be considered one of the most common problems within the instabilities observed in rock masses (Goodman and Kieffer 2000), in addition to being able to Fig. 1 a Photograph of a rock slope, where a typical wedge failure can be observed (Wyllie 2018); b photograph of a slope excavated in a highly fractured rock mass (relatively low joint spacing), where multiple wedge failures occur (photo taken by the authors in Benavente, province of Zamora, Spain) occur in a much wider range of geological and geometric conditions than, for example, planar failure (Wyllie 2018).
The first studies which consider the problem of the stability of a rock block in the face of a translation movement (landslide) and which inspire the analysis of wedge failures were developed by various authors around late 1960s (Goodman 1964;Londe 1965;Wittke 1965;Londe et al. 1969), at a time when the discontinuous nature of rock masses and the great influence of discontinuities in their mechanical behavior was starting to be understood. Most of these referred analyses were based on the limit equilibrium method (LEM), where only the sliding mechanism of the rock block is considered.
The first approaches to the resolution of the particular problem of the stability of a wedge were developed by Hoek et al. (1973), using graphic methods and stereographic projection, techniques which were already well established (Goodman 1964) and later on often resorted to (Hendron et al. 1980). Hoek and Bray (1981) published the book "Rock Slope Engineering" in which the methodology for the analysis of the safety factor of a wedge by using vector calculation is included. This method was later implemented into successive versions of the SWedge software (Rocscience 2020b) for the resolution of the wedge stability problem, although this application is also based on the "block theory" development, as presented by Goodman and Shi (1985), when determining the size of the wedge to carry out the calculations. Some limitations were found when applying the LEM for wedge stability analyses so far, like the need to assume that the block slides in a direction parallel to that of the intersection line. These features became evident for several authors, leading them to develop alternative methods (Chan and Einstein 1981;Chen et al. 1999), but also to the introduction of several parameters that affected the behavior of the joints, such as dilatancy (Wang and Yin 2002) or dynamic effects (Kumsar et al. 2000).
Another disadvantage of the LEM implementation into software packages like SWedge is that they mainly consider, as structural characteristics, the effect of the orientation of the planes for the determination of the maximum wedge but any other parameters that may significantly affect the general stability of the slope, such as the joint spacing. This fact has not been disregarded in recent SWedge versions (Rocscience 2019), where the influence of the joint spacing on the FS can be implemented by means of two levels of spacing ("large spacing" or "low spacing-ubiquitous joint model") in the "probabilistic analysis" type. On the contrary, SWedge has the advantages of allowing a good variety of graphic outputs through a user-friendly interface, besides coupling with other Rocscience softwares, i.e., Dips (Rocscience 2020a). At the end of the 1980s, the effect of the uncertainty, inherent to rock masses, began to be taken into account in the analysis of their stability. This was introduced through the application of statistical techniques in the field of slope engineering (Frangopol and Hong 1987;Tamimi et al. 1989), more or less recently applied to the study of rock wedges (Jiménez-Rodríguez and Sitar 2007;Li et al. 2009;Jiang et al. 2013;Ma et al. 2019). Since 1990, with the development of computing and the improvement of the speed of processors, numerical analyses began to be applied more and more frequently to the resolution of rock engineering problems and, in particular, to the study of wedge failures (Wang et al. 2004). These techniques allowed to evaluate in the last years, based on 3D numerical models, the effect of, i.e., intact rock bridges on the stability of this type of blocks (Paronuzzi et al. 2014(Paronuzzi et al. , 2016.

Methodology
In the present work, stability calculations for each given model of a rock slope affected by two joint sets prone to cause wedge failures were carried out by means of two methodologies. A deterministic method as that implemented in SWedge (Rocscience 2020b) was first selected to perform stability analyses of single wedges. These models were performed for different cohesion and friction angle values, to be able to calibrate 3DEC models involving a single wedge. The 3D discrete element analysis implemented in 3DEC (Itasca 2016a) was used to perform the numerical simulations so as to estimate the factors of safety for all the slope models affected by the same joint-set orientations, but considering different spacing values. Additionally, some numerical models were also studied by generating discrete fracture networks prone to cause wedge failures. The estimation of all factors of safety carried out with 3DEC was obtained based on the shear strength reduction (SSR) technique, by modeling the shear strength of discontinuities with the Mohr-Coulomb criterion.

Definition of the simplified models
Two slope models were created, each from an initial block measuring 220 × 160 × 60 m 3 . These models were defined by the general slope angle (ψ f,1 = 65° and ψ f,2 = 65°), height (H = 50 m), and length (L = 220 m), as presented in Fig. 3. Dip direction of the slope face was set to α = 180°. "Intact rock," that is, those parts of the rock mass not affected by discontinuities, was modeled as rigid (negligible deformation) blocks, with density, ρ = 2700 kg·m −3 .
Each slope model was affected by two joint sets prone to cause wedge failures. Six values of uniformly distributed joint spacing were selected. For these models, joint persistence was considered "infinite." The shear strength of all joint sets was modeled with the Mohr-Coulomb failure criterion (Labuz and Zang 2012), considering six cohesion values and two friction angles. The Mohr-Coulomb criterion was selected to be able to apply, in a simpler way, the shear strength reduction technique for estimating the factor of safety. All geometrical features and mechanical properties of the joint sets defined for these models are presented in Table 1.

Definition of discrete fracture network models
According to Itasca (2016b), a discrete fracture network (DFN) corresponds to a statistical description of joints. A   Fig. 3 a 3DEC model for a slope with ψ f,1 = 65° ("model 1") and b 3DEC model for a slope with ψ f,2 = 90° ("model 2"). Fixed parts of the slope are indicated in green single DFN model can have multiple realizations, being each realization a set of fractures whose characteristics obey the statistical description set in the DFN model. A DFN model considers the fracture population within a rock mass as a set of discrete, planar, and finite size (disk-shaped) fractures. A DFN modeling approach is stochastic, in such a way that the geometrical characteristics are defined through independent statistical distributions of the fracture size (diameter), orientation, and position.
Two DFN models, whose combination in the numerical simulations generated potential wedge failures, were implemented by following similar structural characteristics as selected for the rock slope models presented in the "Definition of the simplified models" section, as described in Table 1. The idea of implementing DFN models was to check if similar trends regarding the effect of spacing on the FS could be observed by studying numerical results for rock slopes presenting a variability closer to that observed in actual rock masses. The characteristics of the implemented DFN models are presented in Table 2.
The mechanical properties of the generated DFNs were assigned in a similar way as for the conventional (uniformly distributed spacings) models, that is, cohesion (c j ) values ranging from 0 to 0.2 MPa (0, 0.01, 0.05, 0.1, 0.15, and 0.2 MPa) and friction angle, ϕ j = 40°.

Calculation of the factor of safety
The factor of safety was firstly estimated through a deterministic analysis by means of SWedge, for a rock slope affected by two joint sets prone to cause wedge failures. Different values of cohesion and two values of friction angle were considered for these calculations, according to Table 1. The factors of safety obtained for each combination of cohesion and friction angle with SWedge (FSi) were used to calibrate the corresponding ones obtained with 3DEC, by considering the highest joint spacing (S ≥ 80 m) in order to generate a single, maximum-size wedge. The estimation of the factor of safety with 3DEC required the application of the SSR method (Zienkiewicz et al. 1975;Shen and Karakus 2013;You et al. 2018), in such a way that an increasing trial factor of safety (f) is used to reduce the shear strength properties of the joints (c j and ϕ j , in the current case) in the way shown by Eqs. 1 and 2 until the slope fails. Failure was considered to be attained when the displacement of any block of the model reached a value of approximately 10 −2 m. At failure, the factor of safety equals the trial factor of safety (f) (Wyllie 2018).
The calibration process is sketched in Fig. 4. After carrying simulations for a single wedge with both programs and considering different cohesion values and a constant friction angle, the factors of safety were found to be very similar, as shown by the results presented in Fig. 5.

Results
In the current section, the results obtained after performing 126 numerical models with 3DEC by combining different values of the shear-strength parameters of the joints (c j and ϕ j ) with the selected spacing values are presented. Complementarily, some numerical models carried out by implementing a DFN approach are provided.
(1)  Represented in Fig. 6 are the obtained values of the FS normalized with respect to FSi after carrying out numerical simulations with 3DEC. In this case, six cohesion values and a constant friction angle ϕ j = 40° were considered, for each value of spacing of J1 and J2. The rock-slope face dip was set as ψ f = 90°.
As it can be seen, for spacing values greater than the height of the slope (S ≥ 80 m), there is hardly any variation in the FS with respect to that obtained with the limit equilibrium approach (FSi). However, as spacing decreases below S/H values = 0.8 (that is, spacing values approximately below the height of the slope), a pronounced decay of the FS with respect to that obtained for the assumed simple wedge model with SWedge is observed. This decay is more pronounced as the initial cohesion of the joints is increased, causing the FS to decrease to minimums of about 20% of its initial value; this occurs for spacings 10 times smaller than the height of the slope and cohesion values of 0.20 MPa. Some examples of the numerical models carried out for a rock slope with ψ f = 90° are presented in Fig. 7.

Rock-slope model with ψ f = 65°, ϕ j = 40° and two joint sets (J1 and J2)
In the case of the rock slope model with a face dip, ψ f = 65°, a similar trend can be observed as in the case of ψ f = 90° (Fig. 8).
For values of the normalized spacing (FS/FSi) greater than 0.8, a slight variation of the FS is observed with respect to that of the instability produced for a simple wedge (S ≥ 80 m). However, the observed decays in FS/FSi are more pronounced than in the case of the slope with ψ f = 90°, reaching a reduction in FS of around 25% of its initial value as the spacing of the discontinuities begins to be of lower order than the height of the slope. In Fig. 9, some examples of the numerical models for a rock slope with ψ f = 65° are presented, for different levels of spacing and cohesion equal to 0.2 MPa.

Rock-slope model with ψ f = 65°, ϕ j = 30° and two joint sets (J1 and J2)
The effect of considering a lower joint friction angle ϕ j = 30° as that considered in the case presented in the "Rock-slope model with ψ f = 65°, friction angle, ϕ = 40° and two joint sets (J1 and J2)" section was also studied. As can be appreciated in the curves representing the normalized FS against the normalized spacing for different cohesion values, a similar decaying trend could also be observed (Fig. 10).
Although the trends are similar to the case presented in Fig. 8, in the present scenario it can be seen that the FS decay is somewhat greater when a smaller friction angle is considered, especially for those models with cohesion. It should be noted that cohesion values lower than 0.05 MPa have not been considered for this case, since they generated models with FS < 1.

Rock-slope model with ψ f = 65°, ϕ j = 40° and two joint sets defined through a DFN
In a similar manner as performed for the conventional (uniformly spaced joint sets) rock slope models, a study of the influence of spacing (indirectly introduced trough the parameter P 10 in 3DEC models) on the factor of safety was carried out. Figure 11 shows normalized results of the factors of safety obtained for rock slope models affected by different levels of joint spacing.
The fact of considering a DFN model prevented the generation of fully persistent joint sets. Based on this modeling limitation and for the sake of brevity, a maximum joint spacing, S = 15 m, had to be selected for the present analysis, since greater S values did not generate wedges within the models unless a large number of numerical models was performed. Based on results presented in Fig. 11, it can be observed a general decaying trend of the normalized FS when the spacing is reduced. It should be pointed out that, regarding the stochastic nature of these models, results will depend, in some manner, on the number of wedges generated in the model, which will depend also on the number of intersections between the joint sets. This fact may explain the scarce dependence of FS results in relation to the cohesion values for those larger values of spacing, whereas a clear dependence between these parameters can be observed for S/H values below 0.1.

Fig. 7
Displacement results after applying SSR method, obtained with 3DEC for the slope "model 2" (ψ f = 90°) affected by joint sets J1 and J2, with six values of decreasing spacing (a-f, correspondingly). For the current case, c j = 0.20 MPa and ϕ j = 40° were used. Overlay-ing the slope model is presented a graph showing the evolution of the trial factor of safety (f) against the number of steps, and the obtained FS Figure 12 provides an example of the numerical models carried out by implementing the DFN approach for generating the joint sets. The value of the FS is also provided for a case with joint sets cohesion, c j = 0.2 MPa and a friction angle, ϕ j = 40°.

Discussion
Even though the main objectives established for the present work have been achieved, it has been considered necessary to carry out a critical review of the obtained results as well as some remarks on the proposed numerical models. It has to be highlighted that this work represents a parametric study aiming at evaluating, from a theoretical point of view and based on 3D numerical models, the effect of spacing on the interpretation of the stability of rock slopes.
Regarding the shear strength behavior of discontinuities, it was implemented through the Mohr-Coulomb model with a non-associated flow rule (zero dilatancy) and flat joints. This representation of the shear strength of discontinuities is far from reality and does not capture the characteristic non-linear shear strength behavior of discontinuities nor does it represent the non-linear behavior of a rock mass. The reason behind the selection of this criterion was, on the one hand, to facilitate the implementation of the SSR technique into the simulations-it is relevant to note that non-linear/Barton-Bandis-type shear strength criteria have barely been implemented in 3D codes, nor SSR technique can be easily applied in such scenarios (Koçal 2008). On the other hand, although involving simplified approaches, to help in the understanding of the phenomena under study from a conservative point of view.
It should also be noted that a full persistence of discontinuities has been considered for the simplified models, dismissing a more realistic scenario including the existence of rock bridges. This aspect has been brought to light by Paronuzzi et al. (2014Paronuzzi et al. ( , 2016 and discussed recently by Bu et al. (2022), who stressed the necessity of not neglecting rock bridges when performing rock slope stability analyses. Nevertheless, the selection of completely persistent discontinuities allowed the comparison of numerical results with those obtained through the LEM approach. Nevertheless, and trying to improve the significance of numerical results in this regard, some DFN models were implemented for some rock-slope scenarios.
The mechanical properties of the discontinuities involved the selection of a realistic set of cohesion and friction angle values ranging from 0 to 0.2 MPa and 30° and 40°, respectively. These properties were assigned to all discontinuities generated in the rock mass and reduced at a time for each considered scenario when applying the SSR technique. The fact of considering uniform and constant cohesion and friction angle values instead of probability density functions, although far from what one can expect in an actual rock mass, helped us in reducing the number of parameters to be controlled when carrying out the numerical models, also facilitating the interpretation of results. In a similar way, the fact of considering uniformly distributed spacing values leaves room for improvement of our work, by considering probability distribution functions as inputs instead of constant values for this parameter, as developed by Priest (1993). Additionally, and regarding the selected values for modeling joint normal and shear stiffnesses, they have been found not to relevantly affect results.
"Intact rock" was modeled as rigid blocks in the numerical simulations carried out with 3DEC. This modeling is appropriate when considering hard rock masses like those that are composed of, i.e., unweathered granite, limestone, or quartzite, where deformation mainly occurs due to block movements. LEM models also do not consider block deformability, so including it in the numerical model would complicate comparisons.
The types of instabilities observed in the simulations generally correspond to wedge failures, although in some cases, especially those involving DFN, the configuration of the structural characteristics of the model causes instabilities that differ from the type of failure under study (blocks different than tetrahedrons and sliding surfaces associated Fig. 9 Displacement results after applying SSR method, obtained with 3DEC for the slope "model 1" (ψ f = 65°) affected by joint sets J1 and J2, with six values of decreasing spacing (a-f, correspondingly). For the current case, c j = 0.20 MPa and ϕ j = 30° were used. Overlay-ing the slope model is presented a graph showing the evolution of the trial factor of safety (f) against the number of steps, and the obtained FS with planar failures). This may affect the general decaying trend of the FS for those models where DFN were implemented. Additionally, it has to be highlighted that the statistical nature of the DFN approach, which prevents the generation of fully persistent joints, adds variability on the obtained results, in such a way that the more spaced ig. 11 Normalized factors of safety (FS/FSi) plotted against five normalized spacings as obtained for the numerical models carried out by implementing two wedge-prone joint sets through a DFN approach, for a rock slope with ψ f = 65°, six cohesion values and ϕ j = 40°P age 11 of 14 1639 the joint sets are, the less the probability of intersections thus the formation of rock wedges. This was particularly reflected on the results obtained for the largest joint set spacings, where a low dependence between joint cohesion and FS could be observed.
Obtaining the factor of safety using the shear strength reduction method can respond to very different criteria (such as the choice of an appropriate convergence criterion or the boundary conditions of the model), which is not the case with the limit equilibrium methods. In the present case, an attempt to find trial factors of safety with a 1% accuracy was made, by following recommendations provided by Itasca (2019).

Conclusions
In the present work, an analysis of the effect of spacing and strength parameters of the discontinuities on the estimation of the factor of safety of rock was carried out.
For this purpose, 126 simulations were carried out using a 3D distinct element code (3DEC) by combining, for two rock slope models, different values of cohesion, friction angles, and spacing of the involved joint sets.
For all the cases studied, a general trend of decreasing factor of safety has been observed as the spacing of the potentially wedge-forming joint sets decreases, an effect that becomes more relevant when higher strength parameters Fig. 12 Displacement results after applying SSR method, obtained with 3DEC for the slope "model 1" (ψ f = 65°) affected by DFN models, with six values of decreasing spacing (a-e, correspondingly). For the current case, c j = 0.20 MPa and ϕ j = 40° were used. Overlaying the slope model is presented a graph showing the evolution of the trial factor of safety (f) against the number of steps, and the obtained FS (specially cohesion) are introduced. These trends were also found in rock slope models where joint sets were generated through a stochastic DFN approach, mimicking the original joint sets produced for the simplified models.
The obtained factors of safety can reach values that represent a 40% decrease with respect to that obtained using models based on the limit equilibrium method (LEM), reducing the analyses only to the orientation of joint sets and slope together with the mechanical properties of discontinuities.
This work shows that the extrapolation of the stability analysis of a simple wedge to the whole rock slope is not totally realistic and may yield results that are on the side of uncertainty, depending on the structural characteristics of the rock mass. Based on this, a special effort should be made on the use of more sophisticated models, especially for the assessment of those rock masses with a relevant degree of fracturing.
The conclusions obtained in this work aim at being useful in the field of rock engineering when analyzing the stability of slopes but also when obtaining more realistic conclusions when making retrospective failure analyses in rock masses.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature, through Universidade de Vigo/CISUG.

Conflict of interest 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/.