The Numerical Simulation of Hard Rocks for Tunnelling Purposes at Great Depths : A Comparison between the Hybrid FDEM Method and Continuous Techniques

Tunnelling processes lead to stress changes surrounding an underground opening resulting in the disturbance and potential damage of the surrounding ground. Especially, when it comes to hard rocks at great depths, the rockmass is more likely to respond in a brittle manner during the excavation. Continuum numerical modelling and discontinuum techniques have been employed in order to capture the complex nature of fracture initiation and propagation at low-confinement conditions surrounding an underground opening. In the present study, the hybrid finite-discrete element method (FDEM) is used and compared to techniques using the finite element method (FEM), in order to investigate the efficiency of these methods in simulating brittle fracturing. +e numerical models are calibrated based on data and observations from the Underground Research Laboratory (URL) Test Tunnel, located inManitoba, Canada. Following the comparison of these models, additional analyses are performed by integrating discrete fracture network (DFN) geometries in order to examine the effect of the explicit simulation of joints in brittle rockmasses. +e results show that in both cases, the FDEMmethod is more capable of capturing the highly damaged zone (HDZ) and the excavation damaged zone (EDZ) compared to results of continuum numerical techniques in such excavations.


Introduction
Significant changes in the stress regime and material properties of a rockmass are the result of the construction of underground openings [1].Induced stresses because of an excavation and subsequent stress redistributions result in damage, fracturing, desaturation, and so on, hence leading to the creation of a zone around the underground opening that the rockmass is disturbed with properties (strength, deformability, permeability, etc.) that differ from the original material properties [2].Understanding the geomaterial response under such conditions is critical to the design and construction of an underground project, especially at great depths.erefore, assessment of the damaged zone around an excavation is of great significance in evaluating the effects on the support requirements and the underground opening stability.
e damage zone surrounding an underground opening is comprised of different subzones depending on the intensity of the induced damage, but these zones are usually referred to collectively as the excavation damage zone (EDZ).
e intensity of the damage usually depends on the distance of an examined point of the rockmass from the excavation, with the increase in distance resulting in the decrease in the influence of the excavation to the surrounding rock [3].By disregarding rockmass damage associated with the construction method employed, as it can be eliminated by taking necessary precautions [4,5], damage associated with stress changes, the excavation geometrical features, preexistent joints within the rockmass, and so on can be evaluated based on fracture coalescence [6] and divided into the highly damaged zone (HDZ), the excavation damage zone (EDZ), and the excavation influence zone (EIZ) [3].
Given the specific requirements of a project, numerical modelling can be utilized in order to provide a better insight of the rockmass response during the excavation and assist in the geotechnical and geological design.Especially, regarding the numerical simulation of rock fracturing, this has been the focus of interest for various researchers [7][8][9][10][11][12][13][14].However, the conducted numerical modelling is only as good as the applied input parameters and assumptions that are integrated into the model, which is required to be able to simulate the expected rockmass behaviour under specific conditions.Within this paper, the short-term mechanical response of an underground excavation within a hard, highly interlocked rockmass is examined by employing two different numerical techniques, the finite element method (FEM) and the hybrid finite-discrete element method (FDEM).For the developed FEM models, the constitutive assumptions adopted include the use of the Hoek-Brown criterion [15] and the damage initiation-spalling limit (DISL) model as proposed by Diederichs [16], while for the developed FDEM model, the finite-discrete element method as proposed by Munjiza [17] is adopted.Furthermore, two different scenarios including a "fracture-free" and a fractured model are utilized by integrating a discrete fracture network (DFN) within the numerical models.For the properties of the intact rock, the well-established case of the Atomic Energy of Canada Limited's (AECL) Underground Research Laboratory (URL) [18] is used in order to examine the mechanical response of the geomaterial during the excavation.Based on the obtained numerical results, it is shown that continuum techniques can capture brittle failure in hard rockmasses at some extent, by making appropriate modifications in the constitutive assumptions.However, certain aspects such as the clear distinction between the highly damaged zone (HDZ) (collapsed material) and the excavation damage zone (EDZ) (fractured material that maintains its structural integrity and does not collapse) are not captured by continuum approaches adequately.On the contrary, the FDEM method appears to be a better fit for simulating brittle fracturing, as it is capable of capturing the complex phenomena associated with brittle failure in hard rock excavations in low-confinement environments.

Short-Term Mechanical Response of Hard Rockmasses and Constitutive Assumptions for Continuum Codes
Instability of underground openings occurs as a result of gravity-driven fallouts which are controlled by the rockmass structure or are stress driven as the rockmass strength is exceeded.In both cases, the rockmass behaviour is controlled by two major factors including the in situ stress regime and the rockmass degree of fracturing [19].Common assumptions within the engineering design include full persistence of joints.However, at great depths, nonpersistent jointing environments away from zones where tectonic processes take place (faulting, shear zones, folding, etc.) are more likely to be encountered.erefore, the presence of rock bridges and the existent joints are the two contributing factors controlling the behaviour of massive and/or rockmasses with nonpersistent joints.For cases in which stress driven rockmass failure is expected to occur, constitutive models that are based on the shear strength of the examined materials are commonly applied.However, their ability to capture the behaviour of massive or moderately fractured rockmasses around underground excavations has been proven to be limited.As documented in the literature, damage within hard rocks is the result of extensile fracturing parallel to the direction of the maximum principal stress σ 1 , as a result of exceeding the tensile strength of the rock [16,[19][20][21]. is results in a lower overall rockmass strength observed in situ, which cannot be predicted by shear failure-based criteria which consistently overestimate the rockmass strength in the numerical models; hence, shear failure based-criteria are not appropriate to use in such cases and other techniques need to be employed.
where σ 1 and σ 3 are the major and minor principal stresses, respectively, σ ci is the unconfined compressive strength of the intact rock, and m b is the reduced value of the material constant for the intact rock m i according to the following: where GSI is the Geological Strength Index [23] and D is a factor which depends on the degree of the ground disturbance to which the rockmass has been subjected to blast damage and stress relaxation.s and α are constants for the rockmass given by the following equations: For underground excavations, such conditions arise for shallow openings, and therefore low confinement conditions (as long as structurally driven failures such as wedge failures, unravelling, etc. are not expected), or deep openings for which the rockmass strength is relatively lower to the in situ stress regime (e.g., squeezing ground conditions).However, as previously mentioned, for brittle rockmasses with nonpersistent joints and of high strength under high stresses, the Hoek-Brown criterion is not appropriate to use.
Brittle damage as a result of induced stresses during an excavation is commonly simulated by employing a cohesion weakening and friction mobilization approach [19,[24][25][26][27], with this approach being used for massive and moderately jointed rocks.Diederichs [16] developed the DISL method in order to capture the response of brittle rockmasses by using 2 Advances in Civil Engineering the generalized Hoek-Brown criterion [22] as a base by using (1). is modi ed approach utilizes the peak and residual strength envelopes of the Hoek-Brown criterion, based on the parameters listed in Table 1, in order to simulate brittle fracturing within conventional and commonly available numerical packages.e method involves the damage initiation as representation of an elevated cohesion and low friction, transitioning to the spalling limit which is represented by cohesion loss and friction mobilization, in order to capture the brittle behaviour of the rockmass at low con nement environments around an underground opening as the tunnel advances.In Figure 1, the Hoek-Brown and DISL strength envelopes are demonstrated.
In this study, the Hoek-Brown criterion and the DISL method are used within the FEM numerical software RS2 [28], in order to be compared to the developed FDEM model as described in the following sections.

The Finite-Discrete Element Method
e nite-discrete element method (FDEM) [17,[29][30][31] is a numerical technique that combines the nite element method with a smeared crack model in order to capture the behaviour of systems that involve discontinuum mechanics under complex deformation, rotation, interaction, and fracturing conditions [32].
e method is based on the discretization of the medium domain into 3-noded, triangular elements that form the mesh.e constant strain triangular elements that comprise the mesh are connected by 4-node cohesive elements which are employed in order to simulate the nonlinear material behaviour ahead of the crack tip due to interlocking and microcracking [33].
e assignment of these cohesive elements between the edges of the adjacent triangular elastic elements allows for the initiation and propagation of fractures once their tensile strength (opening-Mode I), their shear strength (sliding-Mode II), or both (mixed mode-Mode I-II) are exceeded, depending on the stress and deformation state of the medium.e main advantage of this type of simulation is that the potential fracturing paths do not need to be determined a priori.On the contrary, the trajectories of the fractures are controlled by the paths of the induced stresses imposed by the loading conditions.e fracture trajectories though depend on the mesh topology, a limitation that can be overcome as long as the selected element size is small enough so that the fracture patterns can be considered independent of the mesh con guration [34][35][36].
Prior to the occurrence of any fracturing, and as long as the strength of the cohesive elements is not exceeded, in order to maintain the elastic response of the material, an arti cial sti ness for the cohesive elements is introduced in the simulated system by using normal, tangential, and fracture penalty coe cients [13,17].Regarding the strength of the cohesive elements, this is de ned by their tensile and shear strengths as previously mentioned.e tensile strength of the cohesive elements is controlled by the peak tensile strength f t and the fracture energy in tension G I .In a similar fashion, the shear strength is controlled by the peak shear strength f s and the fracture energy in shear G II .e peak shear strength is expressed as the friction coe cient μ, the cohesion c, and the normal stress σ n , as show in (4) [37]: Once the cohesive elements reach their peak strength in Mode I or Mode II, they yield and energy is dissipated through the assigned fracture energy values depending on the failure mode of the cohesive element.Once the fracture energy is depleted, the cohesive element breaks and is removed from the simulation.e forces that act on the nodes of the triangular elements govern the motion of the nite elements according to (5) [17]: Table 1: e equations for determining the DISL model input parameters after [16].e parameters α, s, and m shown here are de ned as material constants based on the crack initiation (CI) stress, the uncon ned compressive strength (UCS), and the tensile strength (T) of the rock.e subscripts s and r stand for the peak and residual values, respectively (adopted from [3]).(modi ed after [17]).

Modelling method
Advances in Civil Engineering 3 where M is the nodal mass matrix, € u is the nodal displacement vector, f int are the internal nodal forces as a result of the deformation of the constant strain elements, and f ext are the external nodal forces that are generated by the imposed external loads, cohesive bonding forces because of the deformation of the unbroken cohesive elements, and contact forces due to the interaction of the broken cohesive elements.e motion equations of the system are solved by an explicit time integration scheme.Once a cohesive element breaks, the two previously connected triangular elements, referred as the contactor and target elements, respectively, interact with one another along the fracture based on the penalty function method [31]: where f c is the contact force, n is the outward unit normal to the penetration boundary Γ c , and φ c and φ t are the potential functions for the contactor and target elements, respectively.As described in the previous section, the primary form of damage of hard rock materials is that of extensile fracturing because of defects and flaws within the geomaterial [21,38,39].e progressive failure of hard rock materials under the low confinement conditions that are encountered around an excavation boundary can be captured from an FDEM model with crack initiation and propagation explicitly according to the nonlinear-elastic fracture mechanics principles [40,41].
For the purposes of this study, in order to capture the brittle behaviour of hard rock excavations at great depths, the FDEM model was developed in the 2D FDEM numerical code Irazu [42], as described in the following sections.

Geological Conditions
In order to study the brittle behaviour of a hard rockmass during the excavation of an underground opening, the material properties of the Lac du Bonnet (LdB) granite and information collected from the case study of the URL Test Tunnel, located in Pinawa, Manitoba, Canada were used [18].e URL Test Tunnel was excavated at a depth of 420 m comprising of a circular cross section of a 3.5 m diameter within the LdB granite in a virtually fracture-free environment [43].Overall, the LdB granite at that depth and within the vicinity of the tunnel can be assumed as a hard, massive, brittle rockmass that is homogeneous and isotropic.Based on data obtained from laboratory testing on intact specimens [44], the mechanical properties of the intact LdB Granite are listed in Table 2. Additionally, the stresses measured in situ at the URL Test Tunnel [43], as shown in Table 3, were used for the FEM and FDEM models in this study.

Geometry, Mesh Configuration, and Excavation
Sequence.
e geometry of both the FEM and FDEM models that were developed follow the geometrical characteristics of the URL Test Tunnel.In Figure 2, it can be seen that both models are comprised of a 60 m × 60 m master domain, with the FEM model having a smaller 15 m × 15 m domain in which the employed mesh is finer.e size of the outer domain was selected as such so that potential boundary effects would be avoided [45].In a similar fashion, the FDEM model is also comprised of smaller subdomains with elements varying in size per domain in order to optimize the computational cost without compromising the accuracy of the obtained results.
e subdomain close to the vicinity of the excavation consists of approximately 0.03 m elements in order to secure that the potential fracture patterns would be independent of the mesh topology, as previously mentioned.
e specifics of the models are listed in Table 4. e equations of motion for the discretized system were integrated with a time step of 6.1 × 10 −8 s for the intact model and 4.2 × 10 −8 s for the fractured model in order to ensure numerical stability for the explicit solver of the code.
Regarding the simulation of the tunnel advancement, the induced three-dimensional (3D) effects on the rockmass are simulated by applying the excavation-induced stresses as a distributed load on the excavation boundary for the FEM models, and the face replacement method [45] for the FDEM models.For the FEM models, after the initialization of the geostatic stresses, the excavation material is removed and substituted with a uniform load applied on the tunnel boundary.is load is gradually reduced to zero in order to cause the destressing effect of the tunnel excavation.For the FDEM models, the advancement of the tunnel involves the replacement of the material within the excavation boundary with unstressed, elastic material during each step, hence leading to the "softening" of the tunnel core and the tunnel circumference to converge.e stability of the numerical analyses conducted both with the FEM and the FDEM models was secured by employing a relatively large number of steps.For the FEM models, the simulation process was comprised of twenty stages in order to simulate the Advances in Civil Engineering advancement of the excavation face, and at each stage, the induced stresses were decreased by 5% until no load was applied on the excavation boundary.For the FDEM models, a large number of steps was employed in order to secure the equilibrium at the establishment of the geostatic stress stage (600,000 steps), the elimination of dynamic e ects during the advancing excavation face (2,200,000 steps), and the subsequent complete removal of the material within the excavation boundary (700,000 steps).

Field Stresses and Boundary Conditions.
e eld stresses assigned to the tunnel models replicate the in situ stress conditions encountered during the excavation of the URL Test Tunnel [43].For the purposes of this study, only the inplane stresses were used for the FDEM models since the adopted cohesive elements do not account for the in uence of the out-of-plane stress on fracture nucleation and growth.On the contrary, the FEM models were assigned an out-ofplane stress of the same magnitude as the recorded intermediate principal stress σ 2 in the URL Test Tunnel, which was 45 MPa. e e ect of gravitational forces and gravityinduced stress gradients were not taken into consideration in the numerical model, and a uniform constant stress eld was assigned instead (deep tunnel assumption).Regarding the far-eld boundary conditions, displacements are xed in the horizontal and vertical directions for both the FDEM and FEM models.Additionally, for the case of the FDEM model, an absorbing boundary condition was employed in order to minimize dynamic oscillations.

Mechanical Input Parameters, Model Calibration, and
Discrete Fracture Networks 5.3.1.Deformability Parameters.In order to determine the deformability parameters that are required for the numerical analysis, the properties of the LdB granite and the in situ rockmass conditions at the URL were assessed.Since the rockmass within the tunnel was virtually free of fractures, it can be assumed that deformability properties (Young's modulus E rm and Poisson's ratio ] rm ) of the rockmass may be approximately the same as these of the intact rock.For the conducted FEM analyses, the modulus and Poisson's ratio values of the intact LdB granite, as listed in Table 1, were used directly into RS2.However, for the FDEM model these values cannot be used directly, as the overall sti ness of the system is a ected by the employed penalty coe cients [36].erefore, the selected modulus and Poisson's ratio values need to be adjusted properly.In order to overcome this, a numerical model of an uncon ned compressive strength (UCS) test is used in order to calibrate the required deformability parameters and penalty coe cients so that the modulus and Poisson's ratio of the LdB granite can be obtained.In Figures 3 and 4, the UCS model con guration and the obtained stress-strain curves are illustrated, respectively.A viscous damping factor is introduced into the model in order to minimize numerical oscillations in the linear part of the stress-strain curves.As observed, before any fracturing occurs, the specimen has a linear-elastic response with the model behaviour governed by the 3-node elastic triangular elements, as long as the strength of the cohesive elements is not exceeded.e established parameters are listed in Table 5.

Strength Parameters.
For the FEM models, there are two major assumptions in order to simulate the strength of   1. e second assumption is that the rockmass is expected to behave according to the DISL model, as it was described in Section 2, with the UCS of the intact rock being in this case as well 200 MPa, and the rest of the parameters of the modi ed Hoek-Brown adopted from [16].e complete set of parameters used is shown in Table 5.
For the FDEM models, once the deformability parameters have been established, the strength properties of the cohesive elements need to be determined in order to replicate the eld conditions and failure mechanisms observed in eld.It has to be noted that for the FDEM method, no constitutive model is required, and the fracturing occurring depends solely on the strength parameters of the cohesive elements.In order to achieve an agreement between eld and model observations, a trial-and-error calibration process was employed.In Figure 5, the formation of the "v-shaped" notch observed in situ and the spalling processes simulated in the FDEM model are illustrated.e established strength parameters of the cohesive elements are listed in Table 5.In Figure 6, a complete owchart of the calibration methodology employed is demonstrated based on [46].

Discrete Fracture Networks.
For the purposes of this study, two types of models were developed within RS2 and Irazu.
e rst type of model was created in order to simulate the fracture conditions encountered at the URL Test Tunnel; hence, no fractures are present [43].e second type of model was assigned a discrete fracture network comprised of one subvertical (mean dip 80 °) and one subhorizontal (mean dip 10 °) joint sets.
e number of the generated fractures is determined by determining two additional parameters: (a) the areal fracture intensity P 21 , de ned as the sum of fracture trace lengths divided by the mapping area [47] and (b) the mean trace length of the fracture traces.e input parameters for the generation of the employed DFN geometries used for both the FEM and FDEM models are listed in Table 6.
e generated DFN geometries were installed within the high-resolution areas of all the models, as described in Section 5.1.Additionally, the fractures are assumed to have a purely frictional behaviour in order to simplify the conducted analyses and the interpretation of the obtained results.It has to be noted that these values have been selected arbitrarily in order to investigate the e ect of the rockmass structure within the FEM and FDEM models, on the development of the fracturing mechanisms, and the extent of the damage based on each numerical method.

Results and Discussion
6.1.Intact Models.In order to assess the impact of the selection of the numerical method on the rockmass response and whether or not it captures the in situ rockmass behaviour, based on the eld observations from the URL Test Tunnel, initially the extent of the potential damage on the   Advances in Civil Engineering rockmass due to the excavation is examined for the models that do not have a DFN assigned to them.In the FEM models, the extent of the damage is initially assessed by identifying the yielded elements around the excavation circumference.Regarding the FDEM model, the numerical results are assessed based on the material that collapses as a result of the fracturing processes taking place.e yielded elements for the FEM models are illustrated in Figures 7 and  8. e collapsed material within the FDEM model is highlighted in Figure 9.
In Figures (7a), 8(a), and 9(a), results obtained for the same amount of deconfinement (40% decrease of the initial stress state) show that all numerical models still have a linear-elastic behaviour, as no fracturing has occurred yet.Based on the major principal stress contours, the magnitude of the monitored stresses is approximately the same, as expected, as the models are still in a prepeak state.However, as further destressing is taking place with the advancement of the face, yielding and fracturing occur.By recording the extent of damage based on the criteria above, in Figure 10, the damage profiles are compared for the fracture-free models between one another and the actual damage profile from the URL Test Tunnel.As it can be observed, the damage profile resulting from the FEM model using the DISL approach resembles the damage profile obtained from the FDEM model.Furthermore, both of them show that they are able to provide a reasonable estimate of the highly damaged area.However, the model using the conventional Hoek-Brown criterion is not able to capture the failure mechanism in a realistic manner.
As discussed in the previous sections, the excavation response of a brittle rockmass cannot be captured by Advances in Civil Engineering conventional shear failure criteria since the failure mechanism is not controlled by the shear strength of the rockmass.More speci cally, in Figure 10, it is shown that at the locations in which eld observations suggest material collapse, the numerical model using the Hoek-Brown criterion shows the lowest extent of damage based on the recorded yielded elements.Furthermore, the shape of the damage extent is completely di erent from the estimations of the other two numerical models and the collapsed material monitored in situ.More speci cally, it resembles more a "butter y" shape than the characteristic "v-shaped notch." is is attributed to the application of the Hoek-Brown criterion which essentially predicts a higher rockmass strength than that of the actual rockmass strength in situ. is higher strength allows the rockmass to withstand the high compressive stresses that manifest at the crown and oor of the excavation due to the .For more information, the reader is referred to [46].

8
Advances in Civil Engineering tunnel advancement.is results in zero stress release at these locations, and it promotes failure at the sidewalls of the excavation which are under a tensile regime.On the contrary, the DISL model efficiently captures the shape of the highly damaged zone (HDZ) and is able to simulate the brittle response of the rockmass based on the modified parameters of the Hoek-Brown criterion.Furthermore, the model predicts tensile failure at the sidewalls of the excavation.While in situ, no tensile cracks were observed, acoustic emission events [49] suggest tensile microcracking around the area that the elements are yielding (Figures 8 and 10).A similar observation can be done within the FDEM model as shown in Figures 9 and 10.
e DISL approach is able to capture the brittle response of the rockmass based on the instantaneous cohesionweakening frictional strengthening model with the Hoek-Brown parameters. is method is relatively simpler than other approaches such as the one proposed by Hajiabdolmajid et al. [27] in which the cohesive strength reduction and frictional strength mobilization reach their post-peak values depending on the magnitude of the accumulated plastic strain.However, the DISL approach does not come without some limitations.As observed in Figure 10, while the FEM model using the DISL approach captures the highly damaged zone (i.e., the collapsed material as a result of the brittle fracturing), the rockmass beyond that area behaves elastically as no more elements are yielding past that extent.Microseismic events (MS) [49], however, suggest that stressinduced damage because of the excavation is extending beyond the highly damaged zone.In FEM models, the elements are assigned a constitutive model in order to simulate the material failure.As a result of this, the HDZ is captured quite well with the DISL approach which provides a good constitutive assumption for simulating brittle fracturing based on the yielding of the elements.However, beyond that yielded material area, the yielded elements cannot be used as a precursor of the damage, and other measured quantities (volumetric strain and principal stress concentrations) need to be assessed in order to evaluate the damage extent [3] and the captured fracturing processes.As observed in Figure 11, a distinction has been made for the FEM-DISL model between the HDZ and EDZ based on the major principal stress.By comparing the results with the recorded MS events, it is shown that the potential EDZ extent is underestimated as the proper mechanics of stress redistribution due to fracturing are not adequately captured.
e FDEM approach, on the other hand, allows for the direct distinction between HDZ and EDZ, and it is more capable of capturing fracturing microprocesses as it allows for fracture propagation and interaction while the overall material maintains its structural integrity.As shown in Figures 9 and  10, this excavation-damaged zone (EDZ) can be captured by the calibrated FDEM model, hence overcoming the limitations of the FEM-DISL model.Furthermore, as observed in Figure 9, the collapsed material fails as a result of extensile fracturing (red lines) due to exceeding the tensile strength of the rock close to the excavation boundary where the confinement is low.erefore, the developed FDEM model and the FDEM method show that are mechanisms capable of simulating the actual failure in order to replicate the in situ conditions.
e accurate simulation of the physical phenomena taking place and resulting in brittle fracturing under high stresses is of great importance, as they affect the design process and the selection of the appropriate temporary support measures.In Figure 12, a monitoring line starting at the crown of the excavation and going into the rockmass for a length of 2.25 m is shown.In Figures 11 and 13, stress measurements taken along this line highlight the difference between each model and field stress measurements conducted at the URL Test Tunnel (SM-5 stress cell) [50] and their subsequent implications in the design process.It can be observed that despite the fact that the in situ measured stresses are not in a very good agreement with the numerical results in terms of the absolute values and location, the FDEM model is capturing the possible stress conditions at the URL Test Tunnel.More specifically, in Figures 11 and 13, the principal stresses σ 1 and σ 3 which were recorded from the numerical models are higher than the stresses recorded at the SM-5 stress cell of the URL Test Tunnel.However, that stress measurement point was located in the disturbed zone around the excavation.
e numerical results from the FDEM model show that the highly damaged zone along the monitored line extends up to 25 cm from the excavation boundary (blue dashed line).Within that range, the stresses monitored are zero, since the material is collapsing.Beyond that distance (blue box), however, a great fluctuation in the recorded stresses can be observed as the monitoring line now lies within the excavation damaged zone.
e occurred fracturing results in a nonsmooth stress distribution along the line, with some points close to the SM-5 stress cell recording similar stress magnitude values.Of course, the exact same location is not the one monitored, a factor that also contributes in not obtaining the exact same value.On the contrary, from the FEM models, a distinction between the HDZ and the EDZ cannot be made instantaneously based on the extent of the yielded elements.In this case, as previously mentioned, the principal stress components were examined in order to determine the extent of the HDZ and EDZ, as suggested in [3].e model using the Hoek-Brown criterion predicts a significantly smaller area affected by the excavation, along the monitoring line.Furthermore, the stresses beyond that distance correspond to these obtained from the material with a linear-elastic behaviour; from acoustic

Advances in Civil Engineering
emission events [49], it is clear that the rockmass is damaged even beyond the collapsed material.For the FEM model in which the DISL approach was applied, according to the major principal stress measurements (change in the slope of the curve in Figure 11), it is shown that the estimated HDZ extent (red dashed line) is approximately the same as the one measured in the FDEM model (blue dashed line).However, the estimated extent of the EDZ is less than the one predicted from the FDEM model, which is more in agreement with the eld observations.Furthermore, it fails to replicate the stress conditions which are monitored by the SM-5 stress cell.
erefore, the conducted stress measurements clearly show the merits of the creation of a more advanced numerical model in order to simulate such rockmass and stress conditions and promote a more e cient support design.

Fractured Models.
In the previous section, the FEM and FDEM models were used in order to replicate the conditions of the URL Test Tunnel for the virtually "fracture-free" LdB granite, based on the use of di erent constitutive assumptions for the FEM models and the employment of the FDEM method.By assuming that the material parameters used in the previous section were representative of the intact rockmass, the same models were modi ed in order to include structure explicitly simulated within them.In order to   Advances in Civil Engineering achieve this, DFN geometry was assigned to each of the numerical models.Because the internal DFN generators incorporated within the Irazu and RS2 codes were used, while the initial DFN parameters were the same for both (Table 6), the exact same geometry was not able to be created for both codes.However, after generating a number of different geometries, two were selected, one in Irazu and one in RS2, with similar joint patterns in order to obtain comparable results.
In Figures 14-16 the high-resolution area of the numerical models containing the DFN geometries are illustrated.In the FEM model cases (Figures 14 and 15), the yielding elements, as in the previous section, are used in order to assess the stress-induced damage due to the tunnel advancement.Within the FDEM model (Figure 16), the excavation damage zone is divided into two separate groups including (a) the material collapsing because of the spalling processes (HDZ) and (b) the damaged rockmass that still maintains its structural integrity due to the con ning stresses (EDZ).In all three di erent cases, it becomes evident that the explicit simulation of joints within the numerical models a ects the stress-induced damage.ese discontinuity elements act both as stress barriers and stress concentrators, and the length and orientation of the rockmass  Numerical results are consistent with the results by [16].
Advances in Civil Engineering joints control both the shape and extent of the stress-induced damage.
More speci cally, in Figure 17, the di erent damage pro les are compared to one another, and they are also compared to the excavation damage pro les obtained from the intact models, as discussed in the previous section.As it can be seen, the fractured models are strongly in uenced by the explicit integration of the DFN geometries within the numerical models in all three cases.In the case of the FEM models (Figures 17(a) and 17(b)) which have the same fracture pattern, it can be observed that similar damage pro les are obtained, both in terms of magnitude and shape.e joints highlighted in Figures 14 and 15 control the in icted damage to the rockmass (yielded elements) around the excavation.However, the use of a di erent constitutive approach in these two models results in a larger damage extent at the crown and oor of the excavation when the DISL approach is employed.On the contrary, when the conventional Hoek-Brown criterion is employed, the crown and oor damage appear to be minimal.In Figures 14 and  15, it can be seen that at the crown and oor of the tunnel, two relatively big areas are fracture free.
erefore, the rockmass behaviour in these two particular regions is controlled by the presence of these rock bridges (intact parts of rock), and the selection of the constitutive model for the rockmass simulation a ects how the material is going to fail.More particularly, the DISL approach is promoting the brittle rockmass failure, as in the case of the intact model.On the contrary, the Hoek-Brown criterion overestimates the rockmass strength resulting in minimal damage at the crown and oor of the excavation.
As discussed in the previous section, the FDEM method can capture the behaviour of brittle rockmasses by simulating the extensile fracturing that is taking place, once the tensile 12 Advances in Civil Engineering strength of the intact rockmass is exceeded, and hence forming fractures that propagate along the direction of the major principal stresses σ 1 .As a result of this, the damage pro les between the FDEM and FEM-DISL models have similarities due to their capability of simulating the brittle response of the excavation as the face advances.As observed in Figures 17(b) and 17(c), the damage in icted on the rockmass at the crown and the oor of the excavation due to brittle failure, where intact parts of rock are present, is captured in both models.However, in the FDEM model, a clearer distinction between the HDZ and the EDZ can be made as result of the complete simulation of fracture initiation and growth.Following the comparison between the fractured models with one another, each intact model is also compared to its corresponding fractured model in order to investigate both the impact of the constitutive model and the employed fracture pattern.In sparsely to moderately jointed rockmasses, the material behaviour is governed by the existence of intact rock parts between the joints present.In this study, joint structure was explicitly simulated, and therefore the behaviour of the medium and the damage extent due to the excavation are controlled by both the constitutive assumptions of the intact rock (rock bridges) and the overall geometry of the joint system.In Figure 17(a), it can be observed that the applied DFN geometry in uences the damage pro le at some extent, especially in areas that the fracture intensity is higher.However, in areas where the in uence of the structure is not that signi cant (excavation crown and oor), the damage extent is similar in both models due to the higher strength predicted by the Hoek-Brown criterion.e overall shape of the damage pro les though is not signi cantly di erent between the two analyses as the overall rockmass behaviour is dictated by the employed constitutive model that promotes yielding towards speci c directions because of the rockmass shear strength, hence decreasing the impact of the preexistent joints.
In the FEM-DISL model, however, that is not the case.As observed in Figure 17(b), by taking into account the brittle response of the intact material, through the application of the DISL approach, results in signi cantly di erent damage pro le for the intact and fractured model.Since the constitutive model of the intact rock promotes brittle failure, this leads to the accumulation of damage at the intact rock parts located at the roof and oor of the tunnel in both cases.However, the presence of structure within the rockmass in the fractured model governs the propagation of the yielded elements.e subhorizontal and subvertical joints close to  In Figure 8, it can be observed that the collapsed material is detaching from the rockmass.On the contrary, the damaged material fractures but maintains its integrity and does not collapse.e results are compared to the actual damaged pro le as observed at the URL Test Tunnel [48], and acoustic emission (AE) and microseismic (MS) events are also plotted for comparison [49].Advances in Civil Engineering 13 the excavation act as stress barriers and stress elevators at speci c locations, hence shaping the formed damaged zone. is becomes more evident at the sidewalls of the excavation.In the intact model, the high horizontal stress leads to a high bending moment at the sidewalls, hence resulting in tensile failure due to bending.However, in the fractured model, that is not happening.e rockmass structure leads to a redistribution of the induced stresses in a way that this bending e ect is not in uencing signi cantly the excavation.On the contrary, it results in accumulation of damage at the upper and lower half of the tunnel.
A similar observation can be made for the FDEM model (Figure 17(c)).e e ect of the preexistent joints in this case becomes evident on both the collapsed material (HDZ) and the damaged material (EDZ).As seen, the collapsed material in the case of the intact model (blue dotted line) is dictated by the applied anisotropic stress regime, with the formation of the notch following the general direction of the geostatic minor principal stress σ 3 .On the contrary, for the fractured model, the notch formation (black continuous line) is in uenced by the occurring extensile fracturing but at the crown is also controlled by a strong subhorizontal joint. is joint redirects the induced stresses and makes the notch to form to the right.at is not the case though for the oor of the excavation (black continuous line).In both models, a large area of intact rock can be seen.erefore, the stressinduced fracturing is mainly the result of spalling, hence leading to similar HDZ pro les for the excavation oor.Regarding the EDZ, the in uence of the DFN geometry is stronger.In the intact model, the EDZ follows the general direction of the HDZ.However, the damaged material is contained at the upper and lower parts around the excavation (blue dash-dot line), with minor tensile fracturing occurring at the sidewalls (blue dotted line) because of bending.In the fractured model, the preexistent joints result in a wider area around the excavation to be damaged (black dashed line).In Figure 16, the rockmass structure controlling the EDZ (white lines) creates speci c blocks of intact material that restrain the stress-induced cracks.More speci cally, strong subhorizontal joints at the crown guide fracture propagation to the right of the excavation, and the  Advances in Civil Engineering subvertical joints at this location act as constraints for the propagating cracks.For the excavation oor, fracturing is intensi ed due to the presence of some joints, but the damage of the material is mainly controlled by spalling due to the presence of rock bridges.Finally, the sidewalls are not imposed to signi cant damage due to the stress redistribution.However, it is evident that preexistent subvertical joints interact with one another resulting in the coalescence of adjacent discontinuities.

Conclusions
e numerical modelling of hard rockmasses has been attempted by employing di erent numerical techniques, including the use of continuum-and discontinuum-based approaches.In the present study, the numerical codes RS2 and Irazu were used in order to perform the numerical simulations using the FEM and the FDEM methods, respectively.
e FEM models simulating the fracture-free rockmass were based on data obtained from the URL Test Tunnel for the LdB granite and relevant research that has been done on that speci c site.Furthermore, the intact FDEM model was calibrated in order to replicate the failure mechanism observed at the URL Test Tunnel.
Once the intact numerical models were established, DFN geometries were integrated into them in order to investigate the impact of joints on the rockmass behaviour during an excavation when di erent constitutive assumptions are made Advances in Civil Engineering and the joints are explicitly simulated.For the FEM models using the conventional Hoek-Brown criterion, for a GSI 80, it was shown that the estimated damage pro le does not correspond to the eld observations for the intact model.Furthermore, the addition of joints does change the shape and the extent of the damage pro le; however, it does not result in signi cant changes to the overall rockmass behaviour.e rockmass constitutive model in this case appears to be the main contributing factor in the rockmass response.However, this does not lead to a realistic damage pro le. e FEM models using the DISL approach have the capability of capturing the brittle response of the rockmass when a continuum technique is employed.For the intact model, the URL eld observations are replicated to an extent, and the simulation of the physical processes taking place is accurate than using the conventional Hoek-Brown criterion.However, HDZ and EDZ cannot be distinguished directly and additional processing is required.On the other hand, the intact FDEM model is capable of simulating the extensile fracturing occurring during spalling, and a clear distinction between the HDZ and the EDZ can be made.While both the FEM-DISL model and the FDEM model provide similar predictions of the HDZ, the EDZ predicted by the FDEM model is in a better agreement with eld observations, as the numerical results are more consistent with AE and MS data from the URL Test Tunnel.Additionally, in situ stress measurements correspond better with stress measurements from the FDEM model.Once joints are added in the FEM-DISL and the FDEM models, it can be observed that similar damage profiles are obtained, since both the effect of the intact rock and the rockmass structure can be adequately captured.However, again the distinction between the HDZ and the EDZ is clearer in the FDEM model where the effect of preexistent joints on the rockmass response is more profound.is highlights the merits of such an advanced numerical technique which has the potential of assisting greatly in the design of an efficient support system for such rockmasses, as the stress and damage state of the material surrounding an excavation can be more realistically captured.

Figure 1 :
Figure1: Schematic of the Hoek-Brown criterion strength envelope and the composite strength envelope of the DISL model[16] (modi ed after[17]).

Figure 2 :
Figure 2: Tunnel model con guration created in Irazu (left) and in RS2 (right).e geometrical characteristics of the excavation and the applied geostatic stresses correspond to the ones recorded and monitored at the URL Test Tunnel.

Figure 3 :
Figure 3: Uncon ned compressive strength (UCS) test for the calibration of the FDEM elastic microparameters based on the deformability properties of the LdB granite.(a) e initial UCS con guration and (b) the UCS specimen at its postpeak condition.Black fractures indicate failure in tension (Mode I), and yellow fractures indicate failure in shear (Mode II).

Figure 4 :
Figure4: Axial stress-axial strain and axial stress-lateral strain curves obtained from UCS testing of the calibrated FDEM model.Recording of stress-strain results ceased once the peak strength was achieved, as the UCS numerical model is used only for the calibration of the elastic parameters of the model.

Figure 5 :Figure 6 :
Figure 5: (a) Photograph of URL Test Tunnel (after [16] modi ed from [24]) showing the damage pro le observed in situ.(b) Damage pro le from the FDEM model (highlighted black) after the completion of the numerical analysis.

Figure 7 :
Figure 7: Major principal stress σ 1 contours of the intact FEM model using the Hoek-Brown criterion (GSI 80).(a) Elastic response of the model at 60% of the induced stresses (tunnel advancing) and (b) nal model state after the completion of the excavation.Yielded elements in shear (×) and in tension (○) are plotted.

Figure 8 :
Figure 8: Major principal stress σ 1 contours of the intact FEM model using the DISL approach [16].(a) Elastic response of the model at 60% of the induced stresses (tunnel advancing), and (b) nal model state after the completion of the excavation.Yielded elements in shear (×) and in tension (○) are plotted.Yielded elements in tension at the crown and oor of the excavation are the result of extensile fracturing.Numerical results are consistent with the results by[16].

Figure 9 :
Figure 9: Major principal stress σ 1 contours of the intact FDEM model.(a) Elastic response of the model at 60% of the induced stresses (tunnel advancing), and (b) nal model state after the completion of the excavation.Cohesive elements failing in tension (Mode I fractures) are coloured red, and cohesive elements failing in shear (Mode II fractures) are coloured yellow.

Figure 10 :
Figure 10: Damage pro les obtained from each intact numerical model (DFN geometries are not integrated in the models): purple continuous line-FEM Hoek-Brown model, red dashed line-FEM-DISL model, blue dotted line-collapsed material of the FDEM model (HDZ), blue dash-dotted line-damaged (fractured) material of the FDEM model (EDZ).In Figure8, it can be observed that the collapsed material is detaching from the rockmass.On the contrary, the damaged material fractures but maintains its integrity and does not collapse.e results are compared to the actual damaged pro le as observed at the URL Test Tunnel[48], and acoustic emission (AE) and microseismic (MS) events are also plotted for comparison[49].

Figure 11 :
Figure 11: Major principal stresses σ 1 monitored along the line (AB) (Figure 12) for the intact FEM Hoek-Brown model (purple line), the intact FEM-DISL model (red line), and the intact FDEM model (blue line).In situ stress measurements (stress cell SM-5) are plotted for comparison [50].e purple hatched area indicates the extent of the damaged zone (yielded elements) for the FEM Hoek-Brown model along line (AB).e red hatched area indicates the extent of the damaged zone (yielded elements) for the FEM-DISL model along line (AB).e blue hatched area indicates the extent of the damaged zone (failure of cohesive elements) for the FDEM model along line (AB).e blue dashed line indicates the ending point of the collapsed material and the starting point of the damaged zone in the FDEM model.e red dashed line indicates the ending point of the HDZ and the starting point of the EDZ in the FEM-DISL model.e red continuous line indicates the ending point of the EDZ in the FEM-DISL model.

1 Figure 12 :Figure 13 :
Figure 12: Monitoring line along the Y axis placed at the crown of the tunnel within the intact FDEM (left) and FEM (right) models for monitoring the principal stresses.e symbols (A) and (B) indicate the starting and ending point of the monitoring line, respectively.e obtained stress results are shown in Figures 11 and 13 for the major and minor principal stresses, respectively.(Left) Cohesive elements failing in tension (Mode I fractures) are coloured red, and cohesive elements failing in shear (Mode II fractures) are coloured yellow.(Right) Yielded elements in shear (×) and in tension (○) are plotted.

R = 1 Figure 14 :Figure 15 :
Figure 14: Major principal stress σ 1 contours of the fractured FEM model using the Hoek-Brown criterion (GSI 80).Yielded elements in shear (×) and in tension (○) are plotted.e dashed black line indicates the damage pro le based on the recorded yielded elements.White lines indicate the contributing joints to the obtained damage pro le.

Figure 16 :
Figure 16: Major principal stress σ 1 contours of the fractured FDEM model.Cohesive elements failing in tension (Mode I fractures) are coloured red, and cohesive elements failing in shear (Mode II fractures) are coloured yellow.e dashed black line indicates the damage pro le based on the monitored fractures.e continuous black line indicates the collapsed material.White lines indicate the contributing joints to the obtained damage pro le.

Figure
Figure Comparison between damage pro les for (a) the FEM Hoek-Brown intact (purple continuous line) and fractured (black dashed line) models, (b) the FEM-DISL intact (red dashed line) and fractured (black dashed line) models, and (c) the FDEM intact (blue dotted line-collapsed material and blue dash-dotted line-damaged material) and fractured (black continuous line-collapsed material and black dashed line-damaged material) models.

Table 2 :
[44]rimental values (average value and range) of the mechanical properties of the intact Lac du Bonnet granite[44].

Table 3 :
Stress field conditions applied in the FDEM model in Irazu (negative values denote compression).

Table 4 :
Model con guration speci cations.or sparse fractures and moderate to good quality discontinuity surfaces) and a UCS value equal to 200 MPa, which falls within the range of the UCS of the LdB granite, as listed in Table [23]tact model con gurations (FEM and FDEM) and the fractured model congurations (FEM-DFN and FDEM-DFN) are listed.therockmass.erstassumption is that the rockmass is going to behave according to the Hoek-Brown criterion for a geological strength index (GSI)[23]value of 80 (massive rockmass with no

Table 5 :
Input material parameters for all numerical models.

Table 6 :
Input parameters for the generation of the discrete fracture network (DFN) geometries.