Role of natural fractures in damage evolution around tunnel excavation in fractured rocks

This paper studies the role of pre-existing fractures in the damage evolution around tunnel excavation in fractured rocks. The length distribution of natural fractures can be described by a power law model, whose exponent a defines the relative proportion of large and small fractures in the system. The larger a is, the higher proportion of small fractures is. A series of two-dimensional discrete fracture networks (DFNs) associated with different length exponent a and fracture intensity P21 is generated to represent various scenarios of distributed preexisting fractures in the rock. The geomechanical behaviour of the fractured rock embedded with DFN geometry in response to isotropic/anisotropic in-situ stress conditions and excavation-induced perturbations is simulated using the hybrid finite-discrete element method (FEMDEM), which can capture the deformation of intact rocks, the interaction of matrix blocks, the displacement of natural fractures, and the propagation of new cracks. An excavation damaged zone (EDZ) develops around the man-made opening as a result of reactivation of preexisting fractures and propagation of wing cracks. The simulation results show that when a is small, the system which is dominated by large fractures can remain stable after excavation given that P21 is not very high; however, intensive structurally-governed kinematic instability can occur if P21 is sufficiently high and the fracture spacing is much smaller than the tunnel size. With the increase of a, the system becomes more dominated by small fractures, and the EDZ is mainly created by the coalescence of small fractures near the tunnel boundary. The results of this study have important implications for designing stable underground openings for radioactive waste repositories as well as other engineering facilities that are intended to generate minimal damage in the host rock mass.


Introduction
Discontinuities are ubiquitous in crustal rocks in the form of faults, bedding planes, joints and veins. Such geological structures associated with mechanical weakness often dominate various properties of the host rock including strength, deformability, permeability and anisotropy (Barton and Quadros, 2014;Lei et al., 2017a). Subsurface rocks embedded with naturally occurring fractures are often encountered in engineering excavations for tunnel construction, hydrocarbon extraction, mining operation and geological disposal of radioactive waste. Underground excavations that perturb the rock mass from an originally equilibrated state can engender stress redistribution with tension, compression and shear occurring in different parts around the opening (Read, 2004). Such perturbations are expected to trigger the formation of an excavation disturbed zone (EdZ) and/or excavation damaged zone (EDZ) (Tsang et al., 2005). For crystalline rocks, EdZ corresponds to the region where only reversible elastic deformation has occurred, whereas EDZ refers to the region where irreversible deformation involving new crack propagation has developed (Hudson et al., 2009;Tsang et al., 2005). Progressive failure in the rock can also lead to a failed zone featured with wedge failure, spalling and even rockbursting at the periphery of an underground cavity (Hoek and Martin, 2014;Liu et al., 2017;Martin and Christiansson, 2009). In the context of radioactive waste repositories in crystalline rocks, the predication of failure depth and position is critical for safety management during the excavation stage and/or the open-drift stage, while the understanding of EDZ properties is important for assessing long-term performance involved with radionuclide migration and mineral dissolution (Cho et al., 2013;Dao et al., 2015;Kwon and Cho, 2008;Tsang et al., 2015Tsang et al., , 2005Walton et al., 2015). Thus, the non-trivial interaction between man-made openings, pre-existing fractures and stress-driven new cracks is worthy of more study.
Damage and failure around excavations in intact rocks has been extensively studied based on various numerical methods including the damage mechanics-based rock failure model (Tang et al., 1998), the bonded-particle discrete element method (DEM) (Potyondy and Cundall, 2004), the elasto-plastic cellular automaton model (Feng et al., 2006) and the hybrid finite-discrete element method (FEMDEM) . In reviewing the literature, it seems that the role of preexisting fractures in EdZ/EDZ formation, which was recognised as a bottleneck issue (Tsang et al., 2005), has not been adequately investigated in the past decade. In this paper, we use the FEMDEM method to simulate the damage evolution around excavation in fractured crystalline rocks embedded with pre-existing fractures following different power law length distributions. We focus on the excavation stage and model the progressive rock mass collapse in an extreme case where no artificial support is introduced. The crystalline rock is assumed to have incompressible grains, i.e. the bulk modulus of grains is much larger than that of the rock, which leads to a Biot-Willis coefficient of 1.0. Geomechanical response of the fractured rock is then modelled based on the Terzaghi's principle, whereas the transient dissipation of pore fluid pressure as well as the dynamic poroelastic coupling is beyond the current scope.
The remainder of the paper is organised as follows. In Section 2, the generation of natural fracture networks is described. In Section 3, the basic principles of the FEMDEM model are presented. The model setup including material properties, boundary conditions and simulation procedures are described in Section 4, while the simulation results are reported in Section 5. Finally, a brief discussion is given and conclusions are drawn.

Fracture networks
Natural fracture systems often exhibit a broad range of fracture sizes, which can be described by a power law statistical model (Bonnet et al., 2001;Davy, 1993;Lei and Wang, 2016): where n(l)dl is the number of fractures with sizes l belonging to the interval [l, l + dl] (dl < < l) in an elementary volume of characteristic size L, a is the power law length exponent, and α is the density term. The only intrinsic characteristic length scales in this model are the smallest and largest fracture lengths, i.e. l min and l max , respectively. In numerical simulations, L is the scale of the modelling domain and should be in between the fracture length bounds, i.e. l min < < L < < l max . The length exponent a defines the manner that the fracture frequency decreases with the fracture size, and extensive outcrop data suggest that a of 2D fracture patterns generally falls in [1.3, 3.5] (Bonnet et al., 2001). A small a value corresponds to a system dominated by large fractures, while a large a value relates to a network of small fractures with their size close to l min . The fracture intensity P 21 , i.e. total length of fractures per unit area, in the domain of size L is related to the first moment of the density distribution of fracture lengths by (Darcel et al., 2003;de Dreuzy et al., 2004): where l′ denotes the fracture length included in the domain of size L (i.e. after some being censored by the domain boundary).
In this study, we choose the modelling domain size L = 20 m. The fracture length bounds are defined as l min = L/100 = 0.2 m, and l max = 100 × L = 2000 m. The spatial and orientation distribution of fractures are assumed to be uniform (i.e. nominally homogeneous and isotropic). We explore three different length exponent scenarios, i.e. a = 1.5, 2.5 and 3.5, and three different fracture intensity cases, i.e. P 21 = 1.0, 2.0, and 3.0 m − 1 . Fig. 1 shows the discrete fracture networks, DFN1-9, generated using different a and P 21 values. It can be seen that when a = 1.5, the system is dominated by very large fractures and has higher connectivity; when a = 3.5, the system mainly consists of small fractures and exhibits lower connectivity; when a = 2.5, both large and small fractures can be important.

Finite-discrete element method
During the 1990s, many algorithmic solutions to discontinuum problems, which later became known as the combined finite-discrete element method (FEMDEM), have been developed (Munjiza et al., 1999;Munjiza et al., 1995;Munjiza and Andrews, 2000). Extensive development and application of the FEMDEM method have been conducted after the release of the open source Y-Code (Munjiza, 2004), with different versions emerged including the code collaboratively developed by Queen Mary University of London in the UK and Los Alamos National Laboratory in the USA (Munjiza et al., 2011), the Y-Geo by the University of Toronto in Canada Mahabadi et al., 2012), and the Solidity platform by Imperial College London in the UK (Latham et al., 2013;Lei, 2016;Xiang et al., 2009). The FEMDEM model that accommodates the finite strain elasticity coupled with a smeared crack model is able to capture the complex behaviour of discontinuum systems involving deformation, rotation, interaction, fracturing and fragmentation (Lei et al., 2017b;Lei et al., 2015bLei et al., , 2014. The key features of the FEMDEM method relevant to the current study are briefly described as follows, while for more details of the method, the reader is referred to (Munjiza, 2004;Munjiza et al., 2011).
The FEMDEM model represents a two-dimensional (2D) fractured rock using a fully discontinuous mesh of three-noded triangular finite elements, which are linked by four-noded joint elements (Fig. 2). There are two types of joint elements: unbroken joint elements inside the matrix and broken joint elements along fractures . The deformation of the bulk material is captured by the linear-elastic constant-strain triangular finite elements with the impenetrability enforced by a penalty function and the continuity constrained by bonding forces of unbroken joint elements (Munjiza et al., 1999). The interaction of discrete matrix bodies is calculated based on the penetration of finite elements via broken joint elements (Munjiza and Andrews, 2000). The joint elements (either broken or unbroken) are created and embedded between the edges of triangular element pairs before the numerical simulation and no remeshing is performed during later computation. The propagation of new cracks is captured by the transition of unbroken joint elements to broken ones, governed by the tensile and shear failure criteria, in an unstructured grid.
The motion of finite elements is governed by the forces acting on the elemental nodes (Munjiza, 2004): where M is the lumped nodal mass matrix, x is the vector of nodal displacements, f int are the internal nodal forces induced by the deformation of triangular elements, f ext are the external nodal forces including external loads f l contributed by boundary conditions and body forces, cohesive bonding forces f b caused by the deformation of unbroken joint elements, and contact forces f c generated by the contact interaction via broken joint elements. The equations of motion of the FEMDEM system are solved by an explicit time integration scheme based on the forward Euler method. The contact interaction between two triangular elements (one is named the contactor and another the target) along fractures is computed based on the penalty function method (Munjiza and Andrews, 2000) as: where n is the outward unit normal to the penetration boundary Γ c , while φ c and φ t are the potential functions for the contactor and target elements, respectively. The normal contact force f n and tangential friction force f t exerted by a contactor onto a target edge are given by , Munjiza (2004): where L p is the penetration length, φ is the potential function along the target edge, v r is the relative velocity (at the Gauss point) between the contactor and the target edge, p is the penalty term, and μ mob is the mobilised friction coefficient which varies during the shearing process for rough fractures (Barton, 1973;Barton et al., 1985): where JRC mob is the mobilised joint roughness coefficient that varies with the displacement δ s , JCS is the joint compressive strength (MPa), σ n is the normal compressive stress (MPa) derived from the Cauchy stress tensor of adjacent finite elements, and φ r is the residual friction angle. Both JRC and JCS are scale-dependent parameters with their Fig. 1. Generated random discrete fracture networks (domain size is 20 m × 20 m) associated with different length exponent a and fracture intensity P 21 .

Fig. 2.
Representation of a 2D fractured rock using an unstructured, fully-discontinuous mesh of three-noded triangular finite elements linked by four-noded joint elements.
Q. Lei et al. Engineering Geology 231 (2017) [100][101][102][103][104][105][106][107][108][109][110][111][112][113] field scale values JRC n and JCS n given by Bandis et al. (1981): where JRC 0 and JCS 0 are the values measured based on a laboratory sample of length l 0 , and l n is the effective joint length (i.e. between fracture intersections) that can be calculated through a binary-tree search of connected broken joint element chains in a complex fracture network . The mobilised JRC mob can be calculated using the dimensionless model (Barton et al., 1985) based on the ratio of the current shear displacement δ s to the peak shear displacement δ peak , which is given by Barton et al. (1985): This roughness-governed non-linear friction law (also called Barton-Bandis criterion) has been implemented into the FEMDEM framework  and validated against the well-established empirical solutions for the direct shear test of different sized fracture samples (Fig. 3). It can be seen that the rough fractures exhibit significant nonlinear shear strength behaviour with the maximum value reached at the peak shear displacement, beyond which the strength gradually decreases to the residual value. As the fracture size increases, a transition from a "brittle" to "plastic" shear failure mode occurs with an increased peak shear displacement and a reduced peak shear strength. Thus, the implemented non-linear friction model permits to capture the important geomechanical behaviour of large fractures, e.g. lower frictional resistance to shearing.

Model setup
The material properties of the natural fractures and intact rock are based on the typical ranges for crystalline rocks (Barton and Choubey, 1977;Lama and Vutukuri, 1978;Zoback, 2007), as presented in Table 1. The energy release rates are estimated based on empirical correlations (Jin et al., 2011;Zhang, 2002). The fractured rock is considered to be at a depth of 500 m, which is the typical depth for geological disposal. The pore fluid pressure is assumed to be hydrostatically distributed in the subsurface with the water table located at the ground surface. Thus, the effective overburden stress is calculated to be S′ v = 8.0 MPa. Two different stress ratio scenarios are considered: S′ H /S′ v = 1.0 and 2.0, corresponding to S′ H = 8.0 and 16.0 MPa, respectively.
The 20 m × 20 m modelling domain embedded with distributed fractures is discretised by an unstructured mesh with the pre-existing fracture walls set as the internal boundaries. An average element size of h ≈ 5.0 cm is used based on the suggestion in the literature . A circular tunnel with a diameter of d = 2 m is placed at the centre of the model domain. The distance from the tunnel centre to the domain boundary is five times the tunnel diameter, for which the boundary effect is considered minor. The bottom of the model is constrained by a roller condition (i.e. no displacement in the y direction) to accommodate the gravitational forces. The effective overburden stress S′ v is applied to the top of the domain while the effective horizontal stress S′ H is imposed uniformly (i.e. gravity-induced gradient is neglected) to the left and right sides. Based on the recommendations in the literature (Mahabadi, 2012), the penalty term p is chosen to be 10 times Young's modulus, i.e. p = 600 GPa, and the damping coefficient η is assigned to be 0.5 times the theoretical critical value, i.e. η = 6.34 × 10 5 kg/m·s.
The plane strain numerical experiment is designed with multiple sequential deformation-solving phases (Fig. 4): (i) the fractured rock is deformed from an unstressed state to an equilibrated state under the insitu stresses loaded through a ramping stage (for avoiding sudden violent failure), (ii) the circular excavation core is relaxed by a gradually reduced deformation modulus to mimic the unloading process as the tunnelling face advances in the out-of-plane direction, (iii) rock materials inside the tunnel are removed, after which an interior free surface is created with no support introduced, and an EDZ progressively develops around the opening with the emergence of new cracks that can link pre-existing discontinuities. Here, the relaxation phase is introduced mainly for the purpose of avoiding unrealistic artificial shocks that can arise by instant removal of excavation materials.

Stress distribution
Figs. 5 and 6 show the distribution of local maximum principal stresses in the fractured rocks embedded with different DFNs under the isotropic in-situ stress condition before and after the tunnel excavation, respectively. Slight stress heterogeneity emerge in the fractured rock before excavation, due to the presence of natural fractures (Fig. 5). After excavation, a heterogeneous stress field can be generated in the  rock, especially in the region close to the tunnel (Fig. 6). In the fracture networks of a = 1.5, the high stress zone around the tunnel tends to be bounded by the neighbouring large fractures and exhibits an anisotropic shape. Very high stresses can occur in the "slab" region between the tunnel boundary and a very close, bypassing large fracture (see DFN4 and DFN7 in Fig. 6). With the increase of P 21 , the effect of preexisting fractures on excavation-induced stress redistribution becomes more significant. In DFN7, a large interior stress loosening zone is formed around the tunnel boundary, where a failure zone involving intensive rock mass failure develops as a result of structurally-governed kinematic instability (e.g. key blocks) and stress-driven breakage (e.g. wing cracks). The shape of the failure zone is strongly affected by the spatial distribution of pre-existing large fractures close to the tunnel. As a is increased, the effect of pre-existing fractures is weakened because of the reduced number of large fractures. In the networks of a = 3.5, the small fractures seem to have minor impacts on the stress redistribution around the tunnel, and only slight rock mass failure occurs at the tunnel periphery.
In the anisotropic in-situ stress field, more stress heterogeneity has been accommodated in the fractured rocks both before (Fig. 7) and after ( Fig. 8) the excavation, compared to the isotropic stress case. In the networks of a = 1.5, stress concentration mainly occurs at the deadends or intersections of large fractures, as a result of sliding of fracture walls and rotation of interlocked blocks under differential stresses, while the local stress fields around isolated, very small fractures seem to be relatively homogeneous (Fig. 7). After excavation, a high stress zone is formed around the tunnel and bounded by large fractures nearby, the effect of which becomes more pronounced as P 21 increases (Fig. 8). Note that the stress loosening zone around the tunnel in DFN7 is smaller than that of the isotropic stress case, which may be attributed to the enhanced shear resistance of some sub-vertical fractures under higher horizontal confinement and the stress arching effect which compacts the blocks around the tunnel. In the networks of a = 2.5, stress concentration mainly emerges at the dead ends of large fractures (Figs. 7 and 8). The high stress zone around the tunnel can also be affected by neighbouring large fractures (Fig. 8). In the networks of a = 3.5, stress concentration occurs at the tips of only a few relatively large fractures, while the local stress fields around most small fractures are very homogeneous (Fig. 7). The high stress zone around the tunnel also tends to exhibit a shape similar to that for excavation in isotropic, Q. Lei et al. Engineering Geology 231 (2017) [100][101][102][103][104][105][106][107][108][109][110][111][112][113] homogeneous media (Fig. 8). The simulation results suggest that the effect of fractures on stress distribution in the fractured rock tends to be enhanced as a decreases or P 21 increases.

Shear displacement
Under the isotropic stress condition, most fractures are quite suppressed for shearing before the excavation (Fig. 9), while only a few large fractures (e.g. in the networks of a = 1.5) intersecting with or very close to the tunnel exhibit certain amount of shearing after the excavation (Fig. 10). In the anisotropic stress field, before excavation, the fracture networks with a = 1.5 accommodate substantial shearing along some preferentially-oriented large fractures (i.e. oblique to S′ H ) (Fig. 11). Some large fractures in the networks of a = 2.5 are also slightly reactivated for shearing. However, most fractures in the networks of a = 3.5 are suppressed for shearing. Hence, the fracture network with smaller a and higher P 21 can accommodate more shear displacement. Furthermore, it seems that the total shear displacement of fractures in the rock after excavation is controlled by the anisotropic  Q. Lei et al. Engineering Geology 231 (2017) [100][101][102][103][104][105][106][107][108][109][110][111][112][113] stress effect, while the excavation-induced perturbation has a minor contribution (Fig. 12). Fig. 13 shows the near-field fracture pattern around the tunnel excavation under the isotropic in-situ stress condition. If P 21 is not high (e.g. P 21 = 1.0 or 2.0 m − 1 ), in the networks of a = 1.5 (i.e. DFN1 and DFN4), the spacing of pre-existing large fractures is greater than or comparable to the tunnel diameter, such that the rock mass made up of interlocking large blocks remains very stable after the excavation; however, in the networks of a = 2.5 and 3.5, more new cracks can propagate from the tips of pre-existing fractures driven by tensile failure and link with other small fractures nearby. If P 21 is high (i.e. P 21 = 3.0 m − 1 ), in the network of a = 1.5 (i.e. DFN7), the spacing of pre-existing large fractures is much smaller than the tunnel size such that intensive failure occurs as a result of structurally-governed kinematic instability; in the network of a = 2.5 (i.e. DFN8), the spacing of large fractures is still greater than the tunnel diameter so that the rock Fig. 7. Stress distribution in the fractured rocks under the anisotropic in-situ stress condition before the tunnel excavation. Q. Lei et al. Engineering Geology 231 (2017) 100-113 mass failure is still manifested by the coalescence of small fractures, which becomes more significant in DFN9 having a higher proportion of small fractures. Generally, under the isotropic stress condition, larger P 21 leads to more new crack propagation given that a is fixed. The EDZ distribution in the fractured rock is characterised using an ellipse that has a minimal volume to cover, say 95%, the area with excavation-induced new cracks (Fig. 14). It shows consistency with the observation based on Fig. 13: if P 21 is small, the damage increases as a increases; if P 21 is large, more damage occurs in the case of a = 1.5. The directional variation of the damage inside each ellipse is further illustrated using a graphic rose with the distance to the tunnel centre also indicated (Fig. 15). The longest spoke of the rose corresponds to the direction with the greatest damage, which is significantly affected by distribution of pre-existing fractures near the tunnel (compare Figs. 13  and 15). Fig. 16 shows the local view of fracture development around the tunnel under the anisotropic stress condition. Figs. 17 and 18 give the spatial and directional distribution of damage in the fractured rocks, Fig. 9. Shear displacement in the fracture networks before excavation under the isotropic in-situ stress condition. Q. Lei et al. Engineering Geology 231 (2017) 100-113 respectively. The fracture networks accommodate more new cracks under the anisotropic stress condition compared to those under the isotropic stress condition. The propagation of wing cracks tends to follow the direction of the maximum principal stress S′ H , but is also affected by the variation of local stress field, especially in regions close to the tunnel. These new cracks can be generated by interaction of small fractures near the tunnel boundary where high stresses are localised, or shearing of large fractures which can initiate new cracks at their tips located far from the tunnel. Thus, when P 21 = 1.0 or 2.0 m − 1 , the damage increases as a increases, as a result of more short fractures coalescing around the tunnel, whereas the role of large fractures is not significant; when P 21 = 3.0 m − 1 , the network of a = 1.5 (i.e. DFN7) has the most damage driven by structurally-governed kinematic instability. A very large EDZ envelope is formed in DFN7, because the propagation of new cracks induced by shearing of large fractures can occur in the rock far from the excavation. Generally, in the anisotropic stress field, for a fixed a value, higher P 21 also leads to more new crack propagation. Fig. 11. Shear displacement in the fracture networks before excavation under the anisotropic in-situ stress condition. Q. Lei et al. Engineering Geology 231 (2017) 100-113

Effects of excavation position
The relative position of excavation to major large fractures may also influence the stress redistribution and damage evolution in the fractured rock. We choose DFN5 as an example and compare the response of the fractured rock in the anisotropic stress field after excavation at three different locations (Fig. 19). If the tunnel cuts in between the two dead-ends of the nearby large fracture (Fig. 19a), the localised high stresses above the tunnel seem to be dissipated by the fracture. When the tunnel just bypasses the large fracture (Fig. 19b), the high stress distribution above the tunnel is obstructed by the fracture, and intensive new cracks are generated in the "slab" region experiencing Fig. 13. Fracture development at the near-field to the excavation boundary in different fracture networks under the isotropic in-situ stress condition.
Fig. 14. The EDZ distribution under the isotropic stress condition is illustrated using damage envelopes covering 95% excavation-induced broken joint elements.

Discussion and conclusions
In this paper, we studied the stress distribution and damage evolution around an unsupported circular tunnel in fractured crystalline rocks based on the FEMDEM simulation. The numerical model realistically captured the EDZ development due to shearing of pre-existing fractures and propagation of new cracks, and the failed zone as a result of structurally-governed kinematic instability and/or stressdriven brittle failure. The EDZ characteristics were found significantly affected by the presence of natural fractures as well as their distribution in the subsurface. If the length exponent a < 2, the system is dominated by the behaviour of large fractures; if a > 3, the damage is Fig. 15. The directional frequency of excavationinduced broken joint elements under the isotropic stress condition is illustrated using the rose diagram. The colour bands show the ranges of distance from the centroids of broken joint elements to the tunnel centre and N total gives the total number of broken joint elements. Fig. 16. Fracture development at the near-field to the excavation boundary in different fracture networks under the anisotropic in-situ stress condition.
Q. Lei et al. Engineering Geology 231 (2017) [100][101][102][103][104][105][106][107][108][109][110][111][112][113] controlled by the interaction and linkage of small fractures; if 2 < a < 3, both large and small fractures play important roles. Thus, the classical effective medium theory that treats fractured rocks using homogeneous representations with equivalent material properties may only be applicable to the scenario of a > 3, while the effects of large fractures may have to be considered explicitly if a < 3, for which a representative elementary volume (REV) may not exist (Bonnet et al., 2001). The observed different failure regimes depending on the length exponent a shows similarity with previous studies of connectivity and permeability regimes of fracture networks, for which the relative proportion of large and small fractures is critical (Bour and Davy, 1997;de Dreuzy et al., 2002, Dreuzy et al., 2001a, 2001bLei et al., 2015a). It is Fig. 17. The EDZ distribution under the anisotropic stress condition is illustrated using damage envelopes covering 95% excavationinduced broken joint elements. Fig. 18. The directional frequency of excavationinduced broken joint elements under the anisotropic stress condition is illustrated using the rose diagram. The colour bands show the ranges of distance from the centroids of broken joint elements to the tunnel centre and N total gives the total number of broken joint elements.
Q. Lei et al. Engineering Geology 231 (2017) [100][101][102][103][104][105][106][107][108][109][110][111][112][113] worth mentioning that the spatial and orientation distribution of natural fractures may not be uniform, but fractal and anisotropic. The clustering and anisotropy of fracture networks were found having important impacts on the connectivity, permeability and strength of fractured rocks (Barton and Quadros, 2014;Darcel et al., 2003;Davy et al., 2006;de Dreuzy et al., 2004;Figueiredo et al., 2015;Ghosh and Daemen, 1993;Harthong et al., 2012;Lei et al., 2014;Wang et al., 2017;Zhang and Sanderson, 1995). The investigation of the fractal spatial distribution of natural fractures with distinguishable orientation sets on damage evolution around excavation will be a focus of our future study. The simulation results in this paper also suggest that the insitu stress condition can also have important effects on reactivation of natural fractures, propagation of wing cracks and stress/damage distribution around the excavation. The research findings of this study have important implications for designing stable underground openings for radioactive waste repositories as well as other engineering facilities that are intended to generate minimal damage in the host rock mass. To draw more definitive conclusions, Monte Carlo simulations may be needed, which, however, requires excessive runtime based on the current processing power (each single DFN model takes~300 h on a desktop computer equipped with an Intel(R) Xeon(R) CPU E5-2697@ 2.30 GHz). Such an issue may be resolved in the future with the rapid development of computational technologies. Fig. 19. The stress distribution, spatial and directional distribution of damage in the fractured rock with a = 2.5 and P 21 = 2.0 m − 1 under the anisotropic stress condition for three different excavation positions relative to the nearby large fracture: (a) the tunnel cuts in between the two dead-ends of the large fracture, (b) the tunnel just bypasses the large fracture, and (c) the tunnel cuts one dead-end of the large fracture.