Overhanging rock slope by design:An integrated approach using rock mass strength characterisation,large-scale numerical modelling and limit equilibrium methods

Pul Schlotfeldt,Dvide Elmo,Brd Pnton

aGolder Associates Ltd.,Vancouver,Canada

bNBK Institute of Mining,The University of British Columbia,Vancouver,Canada

1.Introduction

1.1.Methodology

The purpose of this paper is to present an efficient,yet geologically representative engineering approach to design a large complex overhanging rock slope.The slope in question is at a scale that exceeds the spacing and persistence of joints present in the rock mass,such that intact rock bridging contributes to the stability of the slope.For this problem,it is not the failure of a single discrete block that governs the stability of the overhang,but the complex interaction of stresses that develop due to the eccentric loading from the cantilevered rock section combined with stress rotations between non-persistent joints in the rock mass.Numerical modelling suggests that the overhangs may fail with a combined shear-tensile failure mechanism.That is to say,the shear failure tends to dominate in the lower one third of the overhangs as rock bridge strength is exceeded,while the tensile failure dominates the upper two thirds of the slope as a result of rotation and overturning.The limit equilibrium(LE)model presented herein captures this unique failure mechanism.

This simplified LE method developed allows for multiple support design iterations to be completed in a time-efficient manner for the proposed overhanging rock slopes.The steps used to develop this engineering tool include:

(1)Detailed field mapping of existing slopes combined with a detailed geo-referenced survey(X,Y,Z),which allows for the exact locations of persistent fractures to be spatially represented.

Fig.1.(a)Example of a naturally occurring overhanging rock bluff with basal collapse and an associated blocky basalt overhang;and(b)Critical section through overhanging cliff(after Tsesarsky and Hatzor,2009).

(2)Statistical characterisation of mapping data also allows for geotechnical domains to be delineated for different zones in the rock mass.

(3)Development of a discrete fracture network(DFN)specific to each mapped domain.Since field mapping includes detailed survey of fracture locations,the DFN model is generated with both deterministic fractures(at exact locations of persistent or stability controlling fractures),and stochastically generated fractures to represent joints in three dimensions that are not captured deterministically(calibrated to mapping statistics).

(4)Coupling of the DFN with large-scale explicit numerical strength test simulations to understand possible shear and tensile fracturing processes and to derive synthetic rock mass strength properties that accommodate scale effects.

(5)Full-scale numerical simulations of the slopes to verify the failure mechanism as excavation proceeds.

(6)Development of an LE method to calculate factors of safety(with and without support)that provides a simplification of more complex failure mechanisms and allows for the rapid quantification of rock support as a tool to optimise slope geometries where possible.

It is contended that,provided that the failure mechanism is understood sufficiently,the presented approach allows the designer to make simplified assumptions about the potential shape and location of the failure surface and to apply either shear or tensile strength estimates to discrete parts of the assumed failure surface area.The results of the current study have importance with respect to analysing engineering problems with slopes where confining stresses are very low or absent and failure is not induced due to the development of high compressive stresses or just kinematic block release,but a combination of shear and tensile failures that dictates rock mass behaviour.

This study comes from the work undertaken for infrastructure development where overhanging slopes with complex geometries and varying slope orientations were required and support levels needed to be minimised.The minimised support design criteria required optimising the volume of overhanging rock,which was achieved by leaving a rock buttress at the base of the overhanging slopes.Multiple design iterations were required to achieve this optimisation.While the location of these slopes is not shared,it is contended that,provided that the potential failure mechanism is understood for overhanging slopes,there is(in theory)no limitation on the size of overhanging slopes that can be developed.Rather,the size and length of the anchorage required and the difficulty of construction and associated costs become the limiting factors in the design.

1.2.Current state of practice

Naturally occurring overhanging rock slopes can be found on river banks or sea cliffs(scour),at geological contacts where rock mass strength contrasts are present,at the base of flexural toppling slopes and in the areas where key blocks have kinematic potential to failure from the base of a slope(as shown in Fig.1).In general,large rock overhangs are inherently unstable,partly due to weathering processes and partly due to fatigue and stress concentrations that develop locally across rock bridges that can cause brittle failure.On rare occasions,overhanging rock slopes are required by design,but typically at a scale where the tensile strength of individual blocks becomes the limiting factor rather than at a scale where the slopes are significantly larger than the scale of fracturing.The case history presented here describes the latter case.

The challenge with modelling of rock overhangs that are three dimensional(3D)in nature is that most of the existing commercially available modelling codes do not allow for breakage or brittle fracture of rock to be simulated between non-persistent fractures over large distances.Several commercially available codes used to analyse rock overhangs are discussed in this section,and the limitations of each code are presented for completeness.

Tsesarsky et al.(2005)and Tsesarsky and Hatzor(2009)used the discontinuous deformation analysis(DDA)method to assess the stability of a 34 m high overhanging cliff with tightly spaced bedding planes and three sets of vertical joints.The authors illustrated using the numerical method that the stability of the eccentrically loaded overhanging cliff is determined by the depth of face parallel tension cracks,which were observed behind the crest of the cliff,and the critical length of these cracks was assessed using the DDA method.They also illustrated that the main advantage of the DDA method over the available LE methods is that kinematics can be accounted for in the analysis and the failure mode is not assumed prior to running the model.At the time,the DDA method is superior to the LE analysis and available finite element codes since rock mass structure could be input directly into the numerical simulation.However,breakage of intact rock between fractures is not explicitly considered using this method.Simulation of brittle fracture of intact rock was considered by Carter et al.(2010),who used two-dimensional(2D)discrete element method(DEM)code UDEC(Itasca,2017)to assess recession of overhanging coastal cliffs using the Voronoi tessellation algorithm within the commercially available codes.The advantage of adding Voronoi tessellations to the model is that these randomly oriented,multisided elements can be inserted stochastically between non-persistent fractures.Failure along the boundary of these elements can be used to simulate breakage of intact rock present between joint ends and spacings.Characteristic values for intact rock friction angle,cohesion and tensile strength are used as the constitutive relationship to define the strength along the Voronoi“grain”or element interfaces in the DEM models.These strength values are applied independently to the joint strength estimates,which allows the strength of intact rock and joint frictional strength to be more accurately captured.Strain softening can also be considered using userdefined relationships.However,the number of model elements and processing times provide limitations on Voronoi tessellations for complex joint networks for large-scale simulations.

In the study of fracture initiation and propagation in intact rock,Hoek and Martin(2014)noted that the improvements in computing power in the discrete element codes,specifically bonded particle or grain-based codes,have been useful in studying the tensile cracking behaviour at the laboratory scale.However,Hoek and Martin(2014)noted that“while this approach holds much promise,the current grain-based models are still limited to two dimensions,and much research needs to be carried out before this approach becomes state of practice”.The 3D Voronoi and 3D tetrahedral trigon methods have also been used to model tensile failure at the laboratory scale using the commercially available 3D DEM model 3DEC(Gao et al.,2014;Itasca,2017).If the tensile failure behaviour can be effectively simulated at the laboratory scale in three dimensions,and numerical codes are sophisticated enough to input multiple DFN realisations without significant manual processing and large computational times,it is contended that 2D DEM analysis provides a valid alternative for the modelling of larger-scale brittle fracture problems.

For the reasons outlined by Hoek and Martin(2014),large-scale modelling of overhanging rock sections is typically undertaken using 3D finite element methods.Paronuzzi and Serafini(2009)studied the influence of stress state changes of a historic manmade cantilevered rock cut developed for the Cellina Valley state road in Italy in the early 1900’s.The study required a 3D elastic stress state analysis of the cantilevered overhanging rock slab,and concluded that the rock slab above a small road was subjected to flexural failure when the induced tensile stress exceeded the maximum tensile strength of the intact rock in the manmade cut.Estimates for the tensile strength of rock mass were also provided by Paronuzzi and Serafini(2009)based on the critical stress state observed in the numerical models correlated with point load tests over an estimated area of rock bridging.These estimates were discussed by Zhang et al.(2011),who revisited the critical stress state estimated in the models and suggested that the tensile strengths noted in the initial study may have been overestimated.In both cases,the authors used a 3D continuum code to understand the complex stress redistributions within relatively small overhang,however,the authors noted that the model could not be used to reproduce time-dependant complexities such as rock fatigue,thermal-induced stresses,growth of micro-fractures,and intact rock weakening,nor does the scale of considered overhanging slope take into account the scale effects evident in rock mass behaviour.These types of studies,however,provide insight into the importance of rock bridges between joints that play a role in the stability of general and overhanging slopes.Other continuum codes allow for joint orientation to be considered indirectly in analyses using the ubiquitous joint model.The model proposed by Zhang et al.(2011)does not capture the brittle fracture behaviour between non-persistent fractures,however,it can consider structural geology in continuum models.

For the overhanging rock project presented in this paper,the hybrid finite-discrete element method (FDEM)code ELFEN(Rockfield,2017)was chosen for detailed numerical analysis.In the literature,the FDEM approach has been used by Elmo(2012)and more recently by Mahabadi et al.(2016)to analyse the stability of overhanging rock slopes.In the FDEM code ELFEN,an explicit finite element solution governed by a Mohr-Coulomb elasto-plastic model with a rotating crack(Rankine)tensile cut-off is coupled with a discrete element analysis,thus allowing to model the transition from continuum to discontinuum,typical of rock brittle fracture.Additionally,the ELFEN code allows to embed fracture networks within the continuum mesh and assign specific shear strength properties(Mohr-Coulomb or Barton-Bandis model)to individual fractures or collectively to fracture sets.In this paper,2D FDEM simulations were used to understand the failure mechanisms that could potentially affect the rock slope under consideration,and to estimate rock mass strength parameters to be used for the LE model.

Notwithstanding the type of the geomechanical model used in the analysis(e.g.FDEM,DEM or DDA),the creation of a realistic 3D fracture model is fundamental for the generation of fracture traces used in numerical modelling for slope stability analyses(Spreafico et al.,2015).For the project under consideration,access challenges had to be considered as well,due to the nature of the rock slope,which in parts had to be mapped using a combination of scanline mapping from ropes and terrestrial remote sensing techniques.

The purpose of this paper is not to highlight the use of specific numerical codes,but rather to present a toolbox approach for the systematic design of large overhanging rock cuts,while considering the most likely failure mechanism and the fracture characteristics within the modelled failure path.

2.Representation of fractures in numerical modelling

The insertion of fractures in geomechanical simulations allows for stress rotations and anisotropy to be incorporated in the failure mechanism.It is important to use a numerical code capable of simulating breakage of intact rock,since it is believed that the existence of intact rock bridges between the pre-existing fractures provides the primary strength for the rock mass at the project site.For this reason,it is also important to properly map and characterise the persistence and spacing characteristics of various joint sets present at the project site.Very good 3D rock exposures were available for mapping,and detailed(X,Y,Z)coordinates were taken at key locations with the purpose of:(i)serving as input statistical parameters for the development of a realistic DFN model;and(ii)providing a large data set of deterministic fractures to be included in the DFN model.

2.1.Mapping and statistical characterisation of the rock mass

Geological mapping and rigorous statistical characterisation of both the discontinuity sets and intact rock mass are the essential components of the rock engineering design process.In addition,high-quality LiDAR(light detection and ranging)generated surfaces with orthophotograph overlays are integrated with traditional engineering geological mapping methods during an intense field campaign.Photos and images of the actual mapping process are not included in this paper since that would make it possible to identify the site.Despite the confidentiality terms within which the project was completed,the authors believe that it is important to discuss the basic principles of the mapping approach and share those among the scientific community.

Prior to the field mapping campaign,orthophotograph overlays combined with a surface generated from LiDAR were used to trace the location of persistent discontinuities and target mapping locations.The process of field mapping entailed two entities:data collector(DC)and data recorder(DR).The role of the DC was to select points on given joint traces previously identified on the LiDAR surfaces(only traces longer than 1 m were mapped in the field),mark those points on the rock surface,add a survey number and fracture ID,take orientation reading,assess fracture conditions,and then convey those to the DR using a local area Wi-Fi network connection for direct transferring to the DR laptop.Thus,the DC and DR formed a mapping pair that enabled an almost real-time recoding of the mapped data and the development of a comprehensive database.Despite access constraints,a total of nearly 1600 fractures were mapped both on the ropes and at the crest of the natural rock bluff.The mapping reference sheets used by DR included drop-down boxes with qualitative observations(descriptors and descriptive ranges)which included:(i)rock type identification;(ii)field estimates of intact rock strength;(iii)degree of weathering;(iv)fracture type identification;(v)fracture length(persistence);(vi)fracture spacing;(vii)fracture aperture;(viii)fracture orientation;(ix)fracture surface roughness and waviness;(x)fracture infill;and(xi)fracture terminations.

All the parameters above are included in the International Society for Rock Mechanics(ISRM)suggested methods for the quantitative description of discontinuities in rock masses(Brown,1981);however,specific considerations had to be given with respect to the extent at which fractures could be sampled for DFN purposes and which limitations might have been introduced in the analysis by the chosen sampling methods(Elmo et al.,2015).For instance,the ISRM methods use fracture spacing,whereas the fracture intensity parameter in the DFN method can refer to linear,areal and volumetric intensity.Furthermore,the fracture radius(DFN input)does not necessarily correspond to fracture length or the trace length mapped on the rock surface.The definition of fracture radius also requires geotechnical engineers to understand the concepts of truncation and censoring biases(Elmo et al.,2015).Values below a certain fracture length are omitted(truncation)or relatively larger values cannot be measured because of the limited extent of the rock exposure(censoring).Truncation bias plays a major role in defining the correct fracture intensity for DFN analysis(Rogers et al.,2009).Fracture terminations are seldom collected by geotechnical engineers in the field,but provide a hierarchal structure to relate discontinuity sets.To set a common reference with respect to terminology and avoid ambiguity,the term“fracture”is herein used to define structural features independent of their genesis.

The rock mass at the site was divided into three different geotechnical domains,using fracture intensity and surface conditions(shear strength)as the controlling parameters.Note that the quality of the mapped data and the adopted digitisation process made it possible to include a relatively large number of mapped fractures(deterministic features)within the stochastic DFN model,thus resulting in a hybrid deterministic-stochastic DFN model.The process is described in detail in Section 2.2.

2.2.Discrete fracture network(DFN)analysis

The code FracMan(Dershowitz and Einstein,1988;Golder Associates Inc,2017)is used in the current analysis for DFN data synthesis.As discussed by Elmo et al.(2014),to build a“volumetrically simple”DFN model,only fractures orientation,size,and intensity are required(Fig.2).The spatial distribution of fracturing would assume a constant volumetric intensity within the bounding volume,and the use of an enhanced Baecher model would yield an exponential distribution function for fracture spacing along a sampling line in the model corresponding to the scanline mapped in the field.Whereas this approach may be suited for small-scale rock mass model(e.g.synthetic rock mass models),it is argued that for large-scale problems(e.g.open pits and block cave mines),the characterisation of the rock mass would require a certain degree of spatial ordering.As discussed in Rogers et al.(2009),the process of generating a mine-scale DFN model is much more complex than the one shown in Fig.2.

Fig.2.Workflow used for generation of a “volumetrically simple”DFN model(Zhang and Einstein,2000).

Validation of the DFN model is then achieved using the field data and associated statistics as initial inputs and later comparing the orientation,intensity and pattern of the simulated fracture traces with those measured in the field.

In this work,a disaggregate approach was used to generate the DFN model;accordingly,every mapped fracture set(Table 1)was treated independently with respect to the statistical analyses for fracture orientation,intensity and size.Fracture networks were generated separately for each fracture set and then combined to obtain the overall representation of the fracture network.Fracture orientation and size distributions did not change across the three geotechnical domains.The only exception was the fracture frequencyP10,with geotechnical domain II being characterised by a higher degree of fracturing compared to domains I and III.These domains are spatially linked but discrete from each other.Proper domain identification is important in the case of overhanging rock slopes,since everyset may have specific strength and deformability characteristics.Since geotechnical domain I includes a large number of deterministic(mapped)fractures,there is a need to adjust the input fracture intensityP10in the model,by creating a subdomain in which a hybrid DFN modelling approach is applied(Fig.3).Table 2 indicates the percentage of stochastic fractures generated in geotechnical domain I and subdomain I.

AP10conditional approach is then used when generating the DFN model,such that stochastic fractures are generated using a conditional technique that allows specifyingP10intervals along selected artificial scanlines in the model(Fig.3).Because thefracture spacing measured in the field was corrected to true spacings measured for every fracture set,a true(orientation bias removed)linear fracture intensityP10was used to constrain the DFN model.For each fracture set,multiple scanlines were used and inserted within the bounding volume regions representing the three different geotechnical domains.The software DIPS(Rocscience,2017)was initially used to process the mapped fracture data to identify the different fracture sets.The orientation data were then exported into the DFN code for analysis and simulations.Fig.4 compares the stereonet of the mapped fracture sets with that of a simulated DFN realisation.

Table 1Summary of fracture parameters for DFN analysis.

Fig.3.Region boxes used in the DFN model for generating the fracture networks for the three geotechnical domains and simulated scanlines used in the DFN models to condition each fracture set to the assumed fracture intensities defined in Tables 1 and 2.

Table 2Percentage of stochastic fractures generated in geotechnical domain I and subdomain I.

Fractures mapped in the field were assumed to be discs with a diameter equivalent to the measured fracture trace.It is accepted that the measured trace lengths are typically biased,therefore,there is a need to correct the sampling bias to derive the true trace lengths and corresponding fracture radii.Mapped fracture traces were compared with simulated fracture traces in an iterative process,allowing the selection of the correct fracture radius on a setby-set basis(Fig.5).A visual comparison of mapped and simulated fracture traces served as an additional validation of the DFN model(Fig.6).

The process for the generation of complete DFN model is described in Fig.7 and summarised below:

(1)Initially,volume regions are defined to represent the three geotechnical domains(I,II and III)and subdomain I.These are the regions inside which the stochastic fractures are generated.

(2)Deterministic fractures are imported.Note that these fractures may have an undulated profile.

(3)Stochastic fractures are then generated for subdomain I.

(4)For each fracture set,stochastic fractures are then generated independently within domain I(minus the volume corresponding to subdomain I).The process is repeated for each fracture set in domains II and III.

(5)All fractures are then combined into a final,complete DFN model.

3.Synthetic rock mass strength estimates and overhanging rock failure mechanism verification

Fig.4.Comparison between mapped(left)and simulated(right,one generic DFN realisation)fracture orientation(all geotechnical domains).

The synthetic rock mass(SRM)is a numerically based approach that can be used to establish representative rock mass properties;the approach fully accounts for anisotropic effects and rock mass scale effects.Rock mass classification systems are traditionally used to derive properties for numerical analysis of rock engineering problems.However,those are limited byan inherent assumption of isotropic rock mass behaviour,and by the lack of consideration for intact rock fracturing and failure occurring along the natural fractures.According to Elmo et al.(2016),three main components converge into the SRM approach:(i)data collection and characterisation,(ii)DFN modelling,and(iii)the geomechanical model used for simulating rock mass behaviour.The use of SRM modelling allows for the definition of equivalent Mohr-Coulomb and/or Hoek-Brown strength envelopes,with better consideration for rock mass behaviour at relatively low confinement.Those synthetic strength envelopes could then be incorporated into an LE,a continuum finite element or a finite difference model.SRM models could also be used to define equivalent continuum rock mass properties for a portion of the rock mass up to the representative elementary volume(REV).

Because of their underlying stochastic nature,there are an infinite number of possible realisations of 3D fracture system based on the mapped data.In principles,DFN codes allow to generate a very large number of fracture network realisations while minimising the computing time.However,it is not possible to use a very large number of realisations for the subsequent SRM modelling.The limitation does not impact the validity of the DFN modelling approach as a powerful tool to generate realistic fracture networks;rather this limitation stems from the computing time required to run numerical models that explicitly simulate fracturing mechanisms.In the literature(e.g.Pierce et al.,2007;Mas Ivars,2010;Elmo,2012;Elmo et al.,2016),the authors have typically used a very limited number of DFN realisations to determine a representative range of rock mass properties for the area/volume under consideration.Even for research based projects,using more than 10 DFN realisations to determine SRM properties may not be a viable option due to the computing time required to run large-scale SRM models.Furthermore,the determination of the scale at which the SRM properties become representatives of the rock mass may require running multiple DFN realisations for SRM models at different scales.For real-life engineering projects,the number of DFN realisations that can be considered in the analysis may be further constrained by budget considerations.It is recognised that those aspects may be a limitation of the SRM approach;in this context,Elmo et al.(2016)provided a general workflow for planning an SRM analysis,including main variables and controlling factors(Fig.8).

For the case under consideration,it was decided to simplify the 3D model to a 2D scenario.A cross-section was cut for several 3D DFN realisations and the resulting 2D traces were inserted into the FDEM code adopted for the geomechanical analysis(Fig.9).The FDEM code ELFEN was selected for this purpose.ELFEN is a FDEM code(Rockfield,2017),which has the ability to accommodate a continuum to discrete transition and model the post-failure interaction typical of brittle rock material(Owen et al.,2004).In FDEM models,the finite element-based analysis of continua is merged with discrete element-based transient dynamics,contact detection and contact interaction solutions(Munjiza et al.,1995).Use of fracture mechanics principles in combinationwith a brittle fracture driven continuum-discontinuum transition allows for the simulation of new fractures,creation of discrete blocks and full consideration of the failure kinematics(Elmo et al.,2016).The ELFENcomputational methodology has been extensively tested and validated fully against controlled laboratory tests by Yu(1999),Cai and Kaiser(2004),Yan(2008)and Klerck(2010).The rock mass can be represented in the ELFEN using a variety of constitutive criteria including a Mohr-Coulomb capped with a Rankine rotating crack material model,in which fracturing is controlled by tensile strength and fracture energy parameters.Detailed descriptions of the material constitutive models and fracture mechanics criteria implemented in ELFEN can be found in Owen et al.(2004),Pine et al.(2007)and Klerck(2010).

3.1.In fluence of jointing on tensile and compressive failures of rock overhangs

Fig.5.Comparison between mapped and simulated(one generic DFN realisation)data for fracture traces(all geotechnical domains).

This section describes the analysis undertaken to characterise the potential failure mechanisms of a large overhanging rock mass subjected to gravitational loading only.In the models,the strength of the rock mass is explicitly defined as a combination of intact rock strength and shear properties of the embedded fracture network.For SRM models using FDEM codes such as ELFEN,a calibration process of the material properties is generally not required if a mesh discretization size of mm-scale is used,and intact rock cohesion and friction could be estimated using the generalised Hoek-Brown failure criterion for intact rock material(Elmo et al.,2016).However,careful consideration should be given to upscaling the modelled material properties to large-scale models as the intact rock parameters will be dependent on the assumed mesh discretization size(Fig.10).The approach(Elmo et al.,2016)assumes that the strength of the intact rock matrix would decrease as the mesh size increases,and the reduction in strength is expressed by the formulation proposed by Hoek and Brown(1980).Using this approach,Elmo and Stead(2017)successfully characterised rock mass strength for pillars with increasing width.They presented a chart that clearly illustrated the link between model scale,mesh size and minimum size of induced fracture that can be modelled without using a very computationally expensive remeshing technique.Using those results,and considering that the underlying DFN model includes a very detailed representation of the fracture pattern,including both deterministic and stochastic fractures,a mesh size of 100 mm was used in the current FDEM models.Accordingly,the assumed intact rock strength in the model was calibrated to approximately 88%of the laboratory uniaxial compressive strength(UCS)(50 MPa scaled value).

Fig.6.Visual comparison between mapped and simulated(one generic DFN realisation)traces exposed on the slope surface.Note that the DFN model includes a combination of deterministic and stochastic fractures.

Fig.7.Process for the generation of the complete DFN model.

Fig.8.General workflow for the development of a SRM model(after Elmo et al.,2016).

Fig.9.Section used to define trace maps for the numerical analysis in ELFEN.

Estimated properties of intact rock material and fracture for domains I,II and III used for the simulation are listed in Table 3.Shear strength properties of each fracture set were estimated using field data for joint roughness coefficient(JRC)and joint compressive strength(JCS)incorporated into Barton(1971)and Band is(1993)equations.These equations relate scaled joint roughness and joint wall strengths and the normal stress acting on the sliding surfaces to peak friction angles.Cohesion is typically ignored because shear strength of the joints is determined from wall-towall rock contact along the fracture surfaces(based on field observations).Some cohesion is assumed for set J0 in domain II because of field evidence of healing and mineralisation,which would result in the presence of in-plane rock bridges.

The model in Fig.11 can be considered as a field scale analogue of a large-scale vertical shear box test with a large overhang,in which the shear loading is provided by the gravitational force acting on the right-hand side of the model,while the vertical and horizontal boundaries on the left-hand side are kept fixed.The results clearly show that,for this scale of simulation,shear failure initially develops at the base of the model(stage B).Subsequently,the rock mass becomes progressively damaged(stage C)and the resulting shear failure of key blocks in the lower part of the model induces tensile failure of the rock mass across the upper part of the model,which in turn results in a rotational movement of the entire block(stage D).In stage D,an almost sub-vertical failure plane develops in the model that delimits rock damage on the right-hand side of the model at the beginning at the base of the overhang and the left side of the model with little or no failure.Shear failure is predominant in the lower 1/3 of the model(downwards movement with shearing through intact rock and along pre-existing fractures),while tensile failure(intact rock fails in tension and rock mass rotates outwards)is predominant in the upper 2/3 of the model.

These results form the basis of the LE assessment,and in the authors’opinion provide a reliable two-fold failure mechanisms for use in LE analyses,namely the lower one third of the slope failed initially and principally in shear due to the weight of the overhanging rock mass,followed by the tensile failure of the upper two thirds of the rock mass failing primarily duetodilation and rotation.For the overhanging slope under consideration,set J1 is the subvertical fracture set that runs perpendicular to the section plane.

Fig.10.Example of upscaling approach for intact rock matrix properties in SRM models using FDEM approach(after Elmo et al.,2016).

3.2.Verification of SRM strength estimates

The larger-scale numerical simulation shown in Fig.11 allows identifying the two main types of failure mechanisms for the overhanging rock slope under consideration,i.e.shear and tensile failures.With the mechanism understood,the next step is to provide estimates of rock mass strength under shear and tensile loading conditions to be used as inputs in the LE analysis.

Table 3Estimated properties of intact rock material and fracture for domains I,II and III used for FDEM simulations.

3.2.1.SRM modelling under shear loading condition

Rock mass samples representing the lower portion of the rock block(simplified overhanging slope)shown in Fig.11 were tested to simulate shearing with increasing normal stresses across the near vertical failure surface to predict shear strength envelopes for use in the LE analysis.As discussed earlier,SRM models should include testing several DFN realisations to account for the variability related to the simulated fracture pattern.Note that the influence of the underlying DFN model on the simulated SRM results would also be a function of both the size and the shape of the SRM model(Elmo et al.,2016).To account for that,5 DFN realisations were selected to provide a reasonable range of rock mass shear strength(Fig.12).The modelling results are shown in Fig.13.It is apparent that the fracture network embedded in the model has a strong influence on the estimated shear strength.For example,a difference of approximately 1 MPa in shear strength is observed at zero normal stress between models for DFN realisations R3 and R4.In the context of LE analysis,the results presented in Fig.13 potentially provide a lower bound(considered conservative)and upper bound shear strength estimates for the simulated rock mass.An attempt has also been made in Fig.13 to determine an average shear strength envelope for the rock mass.This could be applied,with judgement,in the LE analysis to estimate the increase in shear strength across the simplified lower portion of the failure plane shown in Fig.11 as the normal stress is increased assuming the introduction of rock anchorage that crosses the near vertical failure plane.

Similarly,SRM modelling of shear mechanisms at different scales(2 m,5 m,10 m and 20 m)was conducted to validate the approach and to demonstrate how the proposed SRM approach could effectively capture the reduction of rock mass strength and equivalent GSI rating as the model size is increased(Fig.14).The Hoek-Brown curves in Fig.14(dashed lines)assume a varying GSI(indicated),miparameter of 28,and UCS of 50 MPa.It was noted that the field estimation of GSI was in the range of 60-65(blocky rock mass,fair to good joint conditions).Despite the limitation imposed using a single DFN realisation,the results clearly show how the adopted SRM model provides reasonable estimates of rock mass strength.

3.2.2.SRM modelling under tensile loading conditions

Rock may have reasonable tensile strength at the intact rock scale.However,at the scale considered in this project,it is reasonable to assume that the natural fractures and the tensile strength provided by potential rock bridges would largely control the behaviour of the rock mass under tensile loading.In most rock engineering problems,the tensile strength of the rock mass is typically neglected.However,given the likely failure mechanism of the overhanging slopediscussed in the previous section,it is argued that the rock mass should not be considered as a totally tensionless material.This is particularly relevant to the upper portion of the simplified failure plane shown in Fig.11,where rock bridges fail in tension and the rock mass rotates outwards.If the rock mass tensile strength was set to zero in the model,it would fail in all cases and this is not considered realistic.

The measurement of the tensile strength for an intact rock specimen is generally carried out using either indirect(Brazilian)or direct(pulling)tension tests.The results of Brazilian tests on core samples are not reported here for confidentiality reason.In the simulated Brazilian test,a cylindrical specimen is loaded diametrically across a circular cross-section of rock,with the applied compressive loading to the disk causing tensile deformation(i.e.splitting open)perpendicular to the loading direction.In this study,simulations of large-scale Brazilian tests were carried out for the fractured rock masses at different scales.Synthetic large-scale direct tension tests were also undertaken and the results were compared to those of the indirect Brazilian tests. The same DFN fracture networks as shown in the top two thirds of the model in Fig.11were considered in the tensile SRM models. All the embedded pre-existing fractures,except for the very persistent sub-horizontal features of set J0, do not have any tensile strength pre-applied to them, and sliding conditions (sliding on pre-existing fractures) are accounted for by using the Mohr-Coulomb strength criterion. A slightly pronounced boneshaped geometrical configuration was adopted for the direct tension SRM models to minimise the possibility of failure being localised at the boundary of the models (Fig. 15).

Fig.11.FDEM results for the simulation of failure behaviour of a large overhanging block under gravity loading.The figure shows the progression from shear failure to tensile failure.20 m×20 m box regions are superimposed to show the relative location of different assumed geotechnical domains.

Fig.16 illustrates the typical modes of failure for the indirect and direct tension SRM models.The results show that failure can occur as a combination of tensile fracturing of the intact rock material(shown in red)and sliding along the existing discontinuities.Brazilian tests typically provide estimates of tensile strength that are generally greater than the equivalent measurement obtained in direct tension tests.The modelled direct tensile strengths are only 2.5%-8.6%of the equivalent indirect tensile strengths estimated in the simulated Brazilian test.It is evident that the modelled rock mass would have limited tensile strength at a scale of 20 m.

An additional model was considered to investigate tensile strength(direct tension)at a much larger scale(40 m),as shown in Fig.17.Clearly,at a much large scale,the tensile strength(direct tension)could be as low as 0.01 MPa.However,this simulation neglects the failure mechanisms simulated in Fig.11 and the fact that direct tension cannot be fully developed if the basal one third of the rock mass is stable or stabilised,i.e.shear failure is not present or is arrested.

3.2.3.Summary of modelling results

Fig.12.Shear model setup and DFN geometries used to study possible uncertainties related to different embedded fracture networks.

With reference to Fig.11,while the potential failure path may be complex,it could be approximated to a near vertical line(i.e.a near vertical surface in three dimensions).It is contended that treating the failure as a discrete plane that intersects the base of the overhanging rock mass and propagates upward is a reasonable simplification for LE work.It is also considered a reasonable simplification to apply the simulated strength estimates(shear strength applies to the lower one third of the failure plane and tensile strength to the upper two thirds of the failure surface)to the vertical failure plane area to establish the resisting force provided by the rock mass along the failure plane.

This concept is developed further below,but in summary the following results were adopted in the LE analysis for determining anchorage support levels for the overhanging slopes:

(1)Vertical rock mass shear strength parameter(τrm)in the lower one third part of the slope was taken as 0.25 MPa.This value is lower than the lowest shear strength estimate(rock mass cohesion)simulated for the large-scale sheartests(R5 model)and is considered suitably conservative.

Fig.13.SRM modelling under shear loading conditions,results for different fracture networks shown in Fig.12.

Fig.14.Variation of SRM strength and estimated correlation between strength and sample size(modified from Elmo,2012).

Fig.15.Synthetic SRM models.Brazilian(20 m in diameter)and direct tension test(20 m in width)models for the three different geotechnical domains.

Fig.16.Results for SRM models(Brazilian and direct tension tests).Zone of fracturing of the intact rock material is shown in red.The modelled tensile strength values are also indicated.

Fig.17.Larger scale(40 m)SRM model for direct tension.Model setup(left)and modelling results(right),in which the zone of tensile fracturing of the intact rock material is shown in red.The modelled tensile strength is also indicated.

(2)Rock mass tensile strength for the upper two thirds of the slope is taken as 0.01 MPa.This value is close to the lowest results for direct tension and it also reflects a balance that includes the scale effects apparent in the results provided in Fig.15.The value is considered to be suitably conservative,given the fact that tensile stresses are induced in the rock mass only as a consequence of shear failure in the lower one third of the slope.Provided that the shear mechanism is mitigated in the design,then tensile stresses should be kept to a minimum and a tensile strength of 0.01 MPa is likely conservative because direct tension cannot be fully developed if the basal one third of the rock mass is stable or stabilised,i.e.conceptually the lower one third of the rock provides a stable pedestal,albeit the fact that it is cantilevered,for the upper two thirds of the overhanging rock mass to sit on it.

However,the modelled SRM results should be viewed in light of the limitations below(Elmo and Stead,2017):

(1)The computational time and cost associated with SRM construction and testing limit the number of DFN realisations that may be considered in the analysis.

(2)The parameterisation of joint properties used within the SRM model.As SRM models require the knowledge of joint properties,the problem of establishing scale dependent relationship for rock mass properties is associated with the extrapolation of large-scale shear strength properties from laboratory tests(i.e.mm-scale).

(3)The inherent variability within a geotechnical domain makes the number of fracture network realisations needed to derive a reliable SRM a mostly undefined problem.Except for either a very closely fractured or an almost massive rock mass,the mechanical response is generally non-uniform due to the orientation,spacing and persistence of the discontinuities(Pine and Harrison,2003).

4.Support design using a graphical 3D limit equilibrium method

In the design of overhanging slopes for the civil infrastructure project in question,it is important that the level of support installed was sufficient to ensure that no major failure occurs within the project design life.It is contended that as numerical models become more sophisticated,the ability to complete sensitivity analysis using a range of input parameters becomes less realistic in the design process(given the length of time it takes to run multiple iterations)and it is therefore important to separate the design process from the numerical modelling process.In the case of this project,the construction process and the design process were run in parallel and answers were required to be faster than the length of time it takes to run multiple iterations.

In the case of the project under discussion,the separation from geomechanical modelling to LE method occurred at the junction where strength parameters were established(described above).LE methods were used to undertake detailed design.Not only is the LE method time-efficient(simplified analysis undertaken with spreadsheets),it uses strength parameters determined by numerical simulations,which account for scale effects in the rock mass,and a simplified failure mechanism,thereby allowing for an optimisation of the engineering design.LE method is further simplified in two distinct calculations;given the failure mechanism,it is considered to be reasonable for developing the force equilibrium and moment equations into two parts instead of considering strength and eccentricity in its entirety as one equation that evaluates both shear and tensile failures at the same time(work flow shown in Fig.18).

The method allows for the calculation of driving and resisting forces and moments based on the changing geometries and provides the designer a method to rapidly calculate the support forces to obtain the required factors of safety(FOSs),which in turn can be devolved into rock support designs for changing dimensions and geometries.The method,which analyses representative slides(Fig.19),also allows for multiple iterations in order that optimisations can be introduced to balance support levels and the possible extent of buttressing required at the base of the slope to ensure stability and construct ability of each slice.In addition,this method can be used to determine the global stability of the entire overhanging slope by summing the driving and resisting moments for the tensile block portion of individual slices and separately summing the driving and resisting forces for the shear block portion of individual slices.

The dimensions of the individual LE slice in the model were defined based on the overall failure path,which in turn is controlled by orientation of joint sets observed in the field.The graphical LE method breaks the overhanging slope into discrete slices or “blocks”(at the scale of the slope)that consider the following changes:slope orientation and angles,dimensions of the slope,and likely location of simplified failure surfaces.The potential failure plane,while vertical in nature,is a curved surface located behind the overhanging slopes.The dip directions of the overhanging slopes almost rotate through a full 360°.Spreadsheets that consider all these variables were developed to facilitate analyses of support forces requirements and hence options for rock support.

4.1.A geologically representative unit block(or slice)to capture the failure mechanism and facilitate the LE analysis

Depending on the orientation of the slope,failure was observed in the numerical models to initially develop along traces of the subvertical joint sets by shear failure along joints and through rock bridges in the lower one third of the failure area of the overhanging slope(Fig.11).In the upper two thirds of the slope,tensile failure dominated as the vertical joints opened up and sliding occurred along the less steeply inclined joint set(J0)and near vertical rock bridges failed in tension between fractures.

Fig.18.Workflow to determine support levels for changing geometries of overhanging rock slopes.

Fig.19.Workflow for analysis of composite limit equilibrium slices for a simple overhanging slope(schematic representation of the slope shown)with a single dip direction.The green part of the sub-vertical failure plane primarily provides shear resistance and the red part of the plane primarily provides tensile resistance.

In general,the sub-vertical joint sets(J1 or J2)are either roughly orthogonal or perpendicular to the overhanging slope depending on the final orientation of the cut face.As described above,a simplified failure surface is assumed to develop in the rock mass that is subparallel to the vertical joint sets as shown in Fig.20.Since the designed slope orientation rotates by roughly 90°,joint sets roughly parallel to the overhanging slope in one location become roughly perpendicular in another location,and the failure plane changes orientation and becomes a near vertical but curved surface.

Slices of 3 m in width were also delineated for the LE method(method of slices)at the scale of the slopes.These slices were cut perpendicular to the overhanging cut slopes.This width(rather than a unit width)was chosen since it approximated to the average spacing between the two sub-vertical joint sets.This is not to say that each slice coincides with the exact sub-vertical planes.Rather it acknowledges the fact that two near orthogonal sub-vertical sets exist in the rock mass with an average spacing of around 3 m and one set is parallel to the back-release failure plane,and the other sets are sub-parallel to the proposed side-release planes for the purposes of LE analysis.Also given the large scale and the complex geometry of the slopes,slices of 3 m in width provided a simplification for the LE analysis to ensure that individual slices have equal or similar dimensions on either edge of the slice but reduced number of slices required in the LE analysis.As shown in Fig.20,planes that delineate individual slices can be considered as potential side-release surfaces that intersect the failure surface or the vertical back-release plane which initiates at the base of the overhang.

After the geometry of the slice has been defined,the area of the shear failure plane(As)and the area of the tensile failure plane(At)are calculated in AutoCAD along the simplified failure surface for individual blocks.The block height(H),block volume(V)and lever arm distance(d)to the centroid of each block are calculated or measured by the same graphical means.

4.2.Force and moment limit equilibrium equations for shear and tensile failures

Fig.20.An example of a single geologically representative limit equilibrium slice divided into a shear block(green)and a tensile block(red).Each slice has a width of 3 m based on a simplification using the knowledge that the spacing of the sub-vertical joint sets(either J1 or J2-depending on face orientation)is approximately 3 m.A is the area.

Fig.21.A cross-section through an overhanging section of a limit equilibrium slice where rock buttressing was used to reduce the size of an overhang.Green area indicates subblock that may fail primarily in shear.Red sub-block indicates area that may fail primarily in tension.

For each LE slice(before support is considered,left-hand diagram of Fig.21),the failure process is simplified into separate force and moment equilibrium equations(Eqs.(1)and(2)).Although computed independently of each other,the failure processes are considered to develop progressively over time with shear developing first and tensile failure later or occurring nearly concurrently.The LE method presented deliberately separates shear and tensile failure modes.The notion is that if shear failure does not occur or cannot occur,the likelihoods of tensile failure diminishes or certainly the volume of rock prone to tensile failure reduces substantially.It is not necessary to develop one equation that considers both eccentricity and compressive and tensile stresses at the same time.This is because,in terms of assessing support level,if the shear sub-block(green area in Fig.21)is stable(self-supporting or with anchorage support),it provides a stable pedestal that the tensile sub-block sitting on the top of overhanging part of the slope diminishes and the potential tensile failure plane moves to the outside edge of the shear sub-block(overhanging rock mass volume is reduced substantially).It is therefore more critical to assess the stability of the shear block and determine the amount of support(bolting and/or increasing the extent of the buttress located below the overhang)required to ensure that the shear sub-block has an adequate FOS so that shear failure does not initiate.Once this is achieved,it provides a safe pedestal to stop or limit the rotation of the upper part of the overhanging slice.

whereFOStensileis the factor of safety for moment equilibrium about the base of the LE block,σtrmis the tensile strength of the rock mass(determined from SRM modelling),Atis the area of the tensile portion of failure plane(from the AutoCAD LE model),His the total height of the LE block at the failure plane(from the AutoCAD LE model),Vis the volume of the overhanging section of rock(from the AutoCAD LE model),γis the unit weight of the rock,anddis the lever arm distance from the block centroid to the pivot point at the base of the overhang(from the AutoCAD LE model).

whereFOSshearis the factor of safety for vertical force equilibrium of the LE block,τrmis the shear strength of the rock mass determined from SRM modelling,andAsis the area of the shear portion of failure plane(from the AutoCAD LE model).

The tensile failure FOS(Eq.(1))is checked first,and it is assumed that the shear failure has already been initiated(the lower one third of the failure plane provides no shear or tensile resistance)and the total volume of the slice produces a driving moment.Thus,Eq.(1)only considers the upper tensile failure when calculating resisting moments.The rationale for calculating the tensile failure FOS first,without considering eccentricity and shear and tensile stresses in its entirety,i.e.one equation,goes tothe excavation sequence of the overhanging rock slopes.

The slopes are to be excavated in a top-down sequence and the tensile sub-block will be exposed first;this will be the part of the overhanging slope that would require sequential support installed as each bench is exposed and access is available.The shear subblock of each slice is the last part of the overhanging slope to be exposed in the excavation sequence and will also be the part of the slopes where the most critical but final support is installed.The goal for the upper part of the slope is to install minimal support and typically only fully grouted dowels are needed to provide confinement to stop relaxation of the rock mass rather than active support.The results of Eq.(1)provide a first-order check to determine if there is a rotational issue.Depending on the height and volume of overhanging rock slope,the moment FOS is typically inadequate and this provides a warning that support in the form of tensioned bolts will likely be required in the shear sub-block portion of the slope to provide a stable pedestal.Since the excavation of the base of the overhang occurs in the final,it then becomes an iterative process to check both Eqs.(1)and(2)and subsequently using Eqs.(3)and(4)to determine a reasonable support level for the shear sub-block once it is finally exposed and allows for the minimisation of the support for the upper tensile sub-block.Again,depending on the results of this analysis,it may be required to shift the rock buttress left in place to reduce the size of the overhanging portion of the rock slope and rationalise the final support configuration.The level and specified type of supports need to be minimised to ensure construct ability,e.g.the length and size of individual bolts need to be optimised such that they can be drilled and installed on narrow bench with limited access.

whereFOSis the factor of safety for vertical force equilibrium of the LE block,Tis the tension force in rock bolts installed at the specified installation angle(minimum to achieve the desiredFOS),αis the installation angle for the post tensioned rock bolts(30°upward from the horizontal),and φ′is the friction angle across the shear failure plane.

whereSis the spacing of rock bolts of a specified type(minimum to achieve the desired FOS),Nis the number of rock bolts of a specified type(minimum to achieve the desired FOS),andPuis the ultimate bolt capacity(manufacturer specification).

In summary,the LE method assumes a lower pedestal that has tensioned rock bolt support and an upper portion of the slope that only requires confinement,i.e.the dowels(only go into tension if movement occurs)provide additional tensile strength to the tensile failure plane(Fig.22)and the rock mass forward of the sub-vertical failure plane.Excavation in a top-down fashion can then proceed safely.The tensile portion of the failure plane provides resistance to tensile failure and the rock mass tensile strength(σtrm)is determined from SRM modelling as described above.The tensile failure plane area(At)is factored into the moment equilibrium as shown in Eq.(1).For the shear portion of the failure plane(lower one third),the shear strength(τrm)is also determined from SRM modelling as described above.This shear strength considers both the frictional resistance of joints and cohesion from intact rock bridges.This part of the surface resists the full block weight-unit weight(γ)and block volume(V)-and is factored into the force(shear)equilibrium as shown in Eq.(2).

As described above,the design criteria adopted demanded that both force and moment equilibriums should be satisfied separately,i.e.FOS≥1.5 is required for both conditions for individual slices,and then support requirements can be quantified for each slice.The global stability of the entire overhanging slope is ensured because the stability of each individual slice is considered first(Eqs.(1)-(4)).This approach is considered to be reasonable because the slice method of analysis ignores the contribution of side-release planes(see Fig.20).The strength of the side-release planes(sub-parallel to set J1)will be similar to the shear strength of the primary failure plane(τrmsub-parallel to set J2)located behind the overhanging slope.However,because individual slices are confined in a direction perpendicular to the side-release planes,the entire surface area of the side-release planes would provide shear resistance,rather than just the lower one third as is the case with the primary failure plane.The analogy is that multiple blocks can be picked up with the block edges aligned vertically without individual blocks dropping to the floor,provided that there is enough clamping applied to the two outside blocks on either end of the suspended stack.In the case of the overhanging slopes,there is enough resistance to shear failure of individual slices along the side release planes for this is not to be considered as a failure mechanism,provided that both ends of the overhanging slopes are bolted in place.

Fig.22.A cross-section through an overhanging block with fully grouted dowels shown conceptually in the tensile failure part of the block(red)and tensioned bolts with a free stressing length installed through the shear part of the same block(green).See Fig.11 to visualise the failure path.

For the project in question,the LE method allowed for the design and installation of support in a time-efficient and sequential manner and ensured that confinement can be maintained and the effects of dilation and shear and tensile failure potentials can be significantly reduced.Tensioned rock supports in the shear sub-block area were designed to be installed at an angle of 30°upward,as shown in Fig.22.The amount of tension force required from bolting of each slice can be calculated using Eq.(3)and the spacing of bolts for the shear sub-block can be calculated using Eq.(4).

5.Conclusions

The explicit consideration of fractures in geomechanical simulations allows for stress rotations and anisotropy to be incorporated in the failure mechanism of large overhangs.Since very good rock exposures were available at the project site,and a detailed geological and geotechnical survey was completed during mapping,a statically representative DFN model with deterministic and stochastic fractures could be generated for the project site.Within the constraints of the limited number of DFN realisations used in the simulations,the authors believed that the SRM models provided reasonable estimates of rock mass strength,and,more importantly,allowed to study the failure mechanisms in detail.

While the sophisticated FDEM models presented in this paper were important for understanding the rock mass strength and failure behaviours,it can be difficult to incorporate rock support in these simulations in a timely manner for multiple design iterations.For this reason,the results from modelling presented herein were used to develop a simplified LE method,which distilled the failure mechanism and pathway through the rock mass,and was used to optimise support requirements.Note that the modelled failure mode would be site-specific,the authors are not suggesting that the modelled failure mechanisms would apply to every overhanging slope.

Because the scale of the slopes and structural geology are unique to this project,and play such a crucial role in determining how the rock mass will behave,the equations presented should not just be applied to other overhanging slope designs without due consideration of alternative failure mechanism.The design methodology and steps presented,however,are considered to be applicable to other projects where overhanging slopes are being considered.In the design of overhanging rock slopes,at an intuitive level anyway,it is important to understand the likely failure through observation and hypothesis.In this study,full-scale numerical simulations were proved useful in developing a refined understanding of the failure mechanism of overhanging cut slopes.In the case presented in this paper,shear failure dominated at the base of the overhanging slopes in question while tensile failure dominated due to eccentric loading.The LE equations presented here reflect this,but will need to be adapted if other failure mechanisms are recognised for geological environments.

Conflict of interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its o utcome.

Bandis SC.Engineering properties and characterization of rock discontinuities.In:Hudson JA,editor.Comprehensive rock engineering,vol.1.Pergamon Press;1993.p.155-83.

Barton N.A relationship between joint roughness and joint shear strength.Nancy.In:Proceedings of the symposium of international society of rock mechanics on rock fracture;1971.Paper I.8.

Brown ET.The ISRM suggested methods for rock characterization testing and monitoring.Oxford:Pergamon;1981.

Cai M,Kaiser PK.Numerical simulation of the Brazilian test and the tensile strength of anisotropic rocks and rocks with pre-existing cracks.International Journal of Rock Mechanics and Mining Sciences 2004;41(Suppl.1):478-783.

Carter T,Cottrell B,Brewster L,Carvalho J,Poschmann A.Cliff recession process modelling as a basis for setback definition.In:Coasts,marine structures and breakwaters:adapting to changes.Thomas Telford Limited;2010.https://doi.org/10.1680/cmsb.41301.0044.

Dershowitz WS,Einstein HH.Characterizing rock joint geometry with joint system models.Rock Mechanics and Rock Engineering 1988;21(1):21-51.

Elmo D,Liu Y,Rogers S.Principles of discrete fracture network modelling for geotechnical applications.In:Proceedings of the 1st conference on discrete fracture network engineering.Vancouver,Canada:DFNE;2014.

Elmo D,Moffitt K,Carvalho J.Synthetic rock mass modelling:experience gained and lessons learned.In:Proceedings of the 50th US Rock Mechanics Symposium,vol.3.American Rock Mechanics Association;2016.p.2208-16.

Elmo D,Stead D,Rogers S.Guidelines for the quantitative description of discontinuities for use in discrete fracture network engineering.Montreal.In:Proceedings of the 13th ISRM congress;2015.Paper 587.

Elmo D,Stead D.Applications of fracture mechanics to rock slopes.In:Feng X-T,editor.Rock mechanics and engineering,Vol.3:analysis,modelling&design.CRC Press;2017.

Elmo D.FDEM&DFN modelling and applications to rock engineering problems.Faculty of Engineering,Turin University;2012.

Gao F,Stead D,Coggan J.Evaluation of coal longwall caving characteristics using an innovative UDEC Trigon approach.Computers and Geotechnics 2014;55:448-60.

Golder Associates Inc.Frac Man.Redmond,USA:Golder Associates Inc.;2017.

Hoek E,Brown ET.Underground excavations in rock.London:Institution of Mining and Metallurgy;1980.

Hoek E,Martin CD.Fracture initiation and propagation in intact rock-a review.Journal of Rock Mechanics and Geotechnical Engineering 2014;6(4):287-300.

Itasca.UDEC,3DEC,PFC user manuals.Itasca Consulting Group Inc.;2017.

Klerck PA.The finite element modelling of discrete fracture in quasi-brittle materials.PhD Thesis.University of Swansea;2010.

Mahabadi OK,Lisjak A,He L,Tatone BSA,Kaifosh P,Grasselli G.Development of a new fully-parallel finite-discrete element code:Irazu.In:Proceedings of the 50th US rock mechanics symposium.American Rock Mechanics Association;2016.Paper ARMA 16-516.

Mas Ivars D.Bonded particle model for jointed rock mass.PhD Thesis.KTH-Engineering Geology and Geophysics Research Group,Department of Land and Water Resources Engineering,Royal Institute of Technology(KTH);2010.

Munjiza A,Owen DRJ,Bicanic N.A combined finite-discrete element method in transient dynamics of fracturing solids.Engineering Computations 1995;12(2):145-74.

Owen DRJ,Feng YT,de Souza Neto EA,Cottrell MG,Wang F,Andrade Pires FM,Yu J.The modelling of multifracturing solids and particulate media.International Journal for Numerical Methods in Engineering 2004;60(1):317-39.

Paronuzzi P,Serafini W.Stress state analysis of a collapsed overhanging rock slab:a case study.Journal of Engineering Geology 2009;108(1-2):65-75.

Pierce M,Cundall P,Potyondy P,Mas Ivars D.A synthetic rock mass model for jointed rock.In:Eberhardt E,Stead D,Morrison T,editors.Rock mechanics:meeting Society’s challenges and demands,proceedings of the 1st Canada-US rock mechanics symposium.London:Taylors&Francis Group;2007.p.341-9.

Pine RJ,Harrison JP.Rock mass properties for engineering design.Quarterly Journal of Engineering Geology and Hydrogeology 2003;36:5-16.

Pine RJ,Owen DRJ,Coggan JS,Rance JM.A new discrete modelling approach for rock masses.Geotechnique 2007;57(9):757-66.

Rockfield.Elfen version 4.7.1.Swansea,UK:Rockfield Software Ltd.;2017.http://www.rockfield.co.uk.

Rocscience.DIPS.Toronto,Canada:Rocscience Ltd.;2017.http://www.rocscience.com.

Rogers S,Elmo D,Beddoes R,Dershowitz B.Mine scale DFN modelling and rapid upscaling in geomechanical simulations of large open pits.Santiago,Chile.In:International symposium on rock slope stability in open pit mining and civil engineering;2009.

Spreafico M,Perotti L,Cervi F,Bacenetti M,Bitelli G,Girelli VA,Mandanici M,Tini MA,Borgatti L.Terrestrial Remote Sensing techniques to complement conventional geomechanical surveys for the assessment of landslide hazard:the San Leo case study(Italy).European Journal of Remote Sensing 2015;48:639-60.

Tsesarsky M,Hatzor YH,Leviathan I,Saltzman U,Sokolowsky M.Structural control on the stability of overhanging,discontinuous rock slopes.In:Alaska rocks 2005,the 40th US symposium on rock mechanics.American Rock Mechanics Association;2005.Paper 771.

Tsesarsky M,Hatzor YH.Kinematics of overhanging slopes in discontinuous rock.Journal of Geotechnical and Geoenvironmental Engineering 2009;135(8):1122-9.

Yan M.Numerical modelling of brittle fracture and step-path failure:from laboratory to rock slope scale.PhD Thesis.Vancouver,Canada:Simon Fraser University;2008.

Yu J.A contact interaction framework for numerical simulation of multi-body problems and aspects of damage and fracture for brittle materials.PhD Thesis.Swansea,UK:University of Wales;1999.

Zhang LQ,Zhou J,Wang XL.Discussion of the paper“Stress state analysis of a collapsed overhanging rock slab:a case study by Paranuzzi and Serafini 2009”.Journal of Engineering Geology 2011;119:120-30.

Zhang L,Einstein HH.Estimating the intensity of rock discontinuities.International Journal of Rock Mechanics and Mining Sciences 2000;37:819-37.