Calibration of coupled hydro-mechanical properties of grain-based model for simulating fracture process and associated pore pressure evolution in excavation damage zone around deep tunnels

Kirsh Frhmnd, Mrk S. Diederichs

a Department of Geological Sciences and Geological Engineering, Queen’s University, Kingston, Ontario, K7L 3N6, Canada

b Golder Associates, Vancouver, British Columbia, Canada

Keywords: Coupled hydro-mechanical properties Excavation damage zone (EDZ)Grain-based model (GBM) calibration Stress-fracturing of rock Cohesive crack model Stress-dependent permeability

ABSTRACT The objective of this paper is to develop a methodology for calibration of a discrete element grain-based model (GBM) to replicate the hydro-mechanical properties of a brittle rock measured in the laboratory,and to apply the calibrated model to simulating the formation of excavation damage zone (EDZ) around underground excavations.Firstly,a new cohesive crack model is implemented into the universal distinct element code(UDEC)to control the fracturing behaviour of materials under various loading modes.Next,a methodology for calibration of the components of the UDEC-Voronoi model is discussed. The role of connectivity of induced microcracks on increasing the permeability of laboratory-scale samples is investigated. The calibrated samples are used to investigate the influence of pore fluid pressure on weakening the drained strength of the laboratory-scale rock.The validity of the Terzaghi’s effective stress law for the drained peak strength of low-porosity rock is tested by performing a series of biaxial compression test simulations. Finally, the evolution of damage and pore pressure around two unsupported circular tunnels in crystalline granitic rock is studied.

1. Introduction

The process of brittle rock fracture during compression and tension is the result of initiation, growth, and coalescence of multiple individual microcracks that eventually lead to the formation of multiple clustered regions of macro-fractures in rock. As the compressive or tensile stress is applied across the boundaries of a rock sample, a complex heterogeneous stress system will be distributed through the rock in which the tensile and shear stresses will be concentrated on pre-existing flaws, which are also termed microcracks, grain boundaries, cavities, and cleavages (Kranz,1983). If the localised tensile stress exceeds the local strength of the microstructure,some microcracks start to form at the point on the boundary of a pre-existing flaw where the maximum tensile stress occurs. These axially aligned extensional microcracks develop during the early loading stages of laboratory compression tests.As the applied deviatoric stresses increase in the sample,the density of compression-induced tensile cracks increases, and eventual interaction and coalescence of these cracks result in the formation of multiple localised and macroscopic damage zones in the material.It is the presence and creation of such micro-fractures that cause the axial stress-strain curve of rock to deviate from true elastic (linearity) one in the pre-failure region (Hazzard et al.,2000).

In addition to pure mechanical effects involved in the deformation of rock,the hydro-mechanical couplings between the solid phase of rock and the pore fluid are also playing a role in the change of stress state and eventual deformation of the material. In this case,if the external load is applied rapidly,the compression of the pore volume causes a rise in pore pressure(P)as the pore fluid does not have time to be dissipated from the pore space. In this undrained hydro-mechanical condition,the rise in pore pressure leads to a reduction in effective stress(σ′=σ-P),and in turn,may cause hydraulic fracturing. However, in the drained condition, the slow loading of the material gives enough time for pore fluid to escape the compressing volume.As a result,the escaping pore fluid allows the pore volume to be consolidated, and in turn, reduces the rock permeability (Rutqvist and Stephansson, 2003).

Various numerical simulation methodologies, such as discrete element method (DEM), hybrid finite element-discrete element method (FEM-DEM), damage mechanics methods, and grain-based DEMs (e.g. particle flow code (PFC) (Itasca, 2008), and universal distinct element code (UDEC)-Voronoi code (Itasca, 2014)), have been utilised successfully to simulate the brittle failure mechanisms of rocks (Bruno, 1994; Diederichs, 2000; Potyondy and Cundall,2004; Kemeny, 2005; Lan et al., 2010; Ivars et al., 2011; Kazerani et al., 2012; Bewick et al., 2014; Farahmand and Diederichs, 2014;Farahmand et al.,2018;Vazaios et al.,2018).The analyses presented in this study focus on simulating the mechanical responses of the rock to account for the role of interaction between the pore fluid and the solid phase of the rock.Therefore,the ability of the grain-based model(GBM)to realistically replicate the influence of pore pressure change on the deformation and weakening the strength of the rock is the main topic of this investigation.

The microcracking-induced permeability variation and poromechanical coupling effects on change of stress state in rock are two important scenarios that the DEM-Voronoi model should be able to model with sufficient accuracy. This study examines the capability of the model to successfully simulate the rock fracture process and its associated poromechanical multiphysics, and proposes a systematic methodology to calibrate the hydro-mechanical parameters of the model. In the DEM-Voronoi model, the rock material is represented by an assembly of polygonal grains bonded through their cohesive contacts. Here, each bond formed between the grains represents an existing plane of weakness defined by the surrounding solid matrix,and is assigned with hydraulic properties that are described by the cubic law (Witherspoon et al.,1980). In this way,the fluid flow in the model is governed by the aperture of the Voronoi joints,and in return,the joint aperture is a function of the stress state and fluid pressure acting on the joint boundary. In GBMs,it is generally the micro-mechanical properties of the model,which include those of the grains and contacts that produce the emergent macro-mechanical properties of the bulk rock. As a result,when modelling the fracture development behaviour of the porous media filled with fluid, the first step is to calibrate the hydro-mechanical micro-parameters of the model to match the macro-mechanical and hydraulic material properties.

In order to establish a holistic modelling technique to simulate the evolution of cracked-modified rock properties and associated pore fluid behaviour during rock deformation, the proposed technique should be able to reproduce a number of key mechanisms.Firstly,the model(e.g.uniaxial or triaxial compression test models)should be able to correctly simulate the mechanical response in the pre- and post-failure states of the material. Hence, the model should be able to successfully simulate the pre-failure elastic response, strength levels such as crack initiation threshold, crack damage threshold,and peak compressive strength(σy),and finally dilation angle (ψ) evolution in the post-failure state (Walton and Diederichs, 2015). In terms of hydro-mechanical coupling processes,four poromechanical effects are required to be simulated to mimic the complex interplay occurring between the pore fluid and solid components of the rock. These effects include: (i) the deformation- and damage-induced permeability changes (Zoback and Byerlee, 1975); (ii) the effective stress change due to pore pressure evolution (Biot, 1941); (iii) the pore pressure change due to change in stress state of the porous medium(Skempton,1954);and(iv) the effect of pore pressure on weakening or strengthening of the rock strengths (crack initiation and damage thresholds, and peak compressive strength) (Brace and Martin,1968).

The applied model should be calibrated in order to consider the permeability change of the bulk rock in two fundamental mechanisms:(i) the initiation and growth of newly developed cracks that lead to permeability increase of the material, and (ii) the pore volume contracts or dilations as a function of stress,which therefore lead to a change of rock permeability.Correct numerical modelling of rock permeability change is of importance when evaluating the permeability evolution adjacent to an underground excavation in both the undamaged and damaged regions, which are termed the excavation influence zone(EIZ)and excavation damage zone(EDZ),respectively (Tsang et al., 2005; Perras and Diederichs, 2016). The relationship between the degree of rock deformation and associated permeability evolution in both intact rock sample and rock mass surrounding a hypothetical tunnel is illustrated in Fig.1.

The second poromechanical effect that must be addressed is the effect of pore pressure change on the magnitude of effective stress(σ′) experienced by the rock. By modelling this effect, the pore pressure influence on changing stress will be simulated accordingly. Based on Biot’s theory of consolidation (Biot,1941), when a porous material is saturated with fluid, the stress states that the material experiences are governed by

Fig. 1. Schematic of the excavation damage zones (EDZs) formed around a tunnel section and corresponding state of stress damage-permeability for each sub-zone of the EDZ. Regions I, II, III and IV represent the elastic deformation (EIZ), isolated stable cracking (outer EDZ, EDZouter), strain-hardening unstable cracking (inner EDZ, EDZin-ner), and post-peak (highly damage zone, HDZ), respectively (modified after ANDRA(2005)). The photograph shows a cut section of the EDZ fractures at the tunnel wall at the ?SPO Hard Rock Laboratory in Sweden (Emsley et al.,1997).

where α is a poroelastic material constant known as Biot’s coefficient(Biot,1941). The numerical modelling in this study examines whether the reaction forces occurring between pore fluid and grains when subjected to loading can reproduce the macroscopic stress state in a large volume of the rock according to the Biot’s effective stress law. In this case, the macroscopic response of the model during elastic deformation (non-damage state) should therefore follow the effective stress law.

In terms of continuum mechanics theory, it is Skempton’s coefficient(B),defined as the ratio of the induced pore pressure to the change in applied stress for undrained conditions, which controls the relationship between the applied confining stress and the pore fluid pressure induced in the sample (Skempton,1954). This pore pressure change in rock in an undrained condition has a profound influence on the deformation and state of damage in rocks. This undrained response and the associated pore pressure evolution are especially important in cases where the rock porosity is very low such that the pore fluid does not have enough time to escape from the sample.As a result of undrained loading,rock fracturing may be initiated due to the sudden pore pressure increase induced by the change in loading condition. Hence, the micro-components of GBMs should be calibrated in a way that the pore pressure response of the bulk model during rock elastic deformation is replicated correctly. The ability of the model to simulate such fluid pressure changes may significantly affect the shape and extent of the EDZ developed around excavations.

Another important topic is how effectively the GBM can account for the influence of pore pressure on weakening the rock strength.There is a considerable body7 of experimental evidence showing that the failure strength of a wide variety of rocks is a unique function of σ3-P,where σ3is the total confining stress.As a result,to a fairly good approximation, Terzaghi’s effective stress law(σ′3= σ3-P) (Terzaghi, 1943) governs the failure strength of porous rocks (Serdengecti and Boozer, 1961; Handin et al., 1963;Murrell, 1965; Robinson and Holland, 1969; Dunn et al., 1973;Byerlee,1975;Gowd and Rummel,1997;Paterson and Wong,2005;Jaeger et al., 2009). An example is given in Fig. 2. In general, the differential stress for the fracture strength of a particular rock is approximately the same as the effective confining pressure, when the latter is taken as the total confining pressure minus the pore pressure.According to this,the capability of the model needs to be tested based on Terzaghi’s effective stress notion that the failure of the rock would be controlled by Terzaghi’s effective stress law.

The remainder of this study is organised as follows.In Section 3,the implementation of a constitutive model,referred to as“cohesive crack model”, to control the contact behaviour is presented. The constitutive model defines the fracture behaviours of the contacts formed between grain boundaries in tension(mode-I),shear(mode-II), and mixed mode (mode I-II). The model determines the contact forces based on the relative displacement that the contact has undergone according to a set of stress-displacement laws,and accounts for the effects of material softening in front of the crack tip known as the fracture process zone (Dugdale, 1960) after the bond between two contacts breaks. In Section 4, a systematic process to calibrate the hydro-mechanical parameters is described when the grains making up the model are distributed according to the real mineral composition of the Lac du Bonnet(LdB)granite.The appropriateness of the DEM-Voronoi model to reproduce the coupling effects between fluid flow and rock deformation during fracturing of a laboratory-scale sample is investigated in Section 4.3.1.The ability of the calibrated model to mimic the rock permeability evolution coupled with the microcracking phenomenon is examined by simulating a biaxial compression test.Later,the ability of the model to replicate correct values of Biot’s and Skempton’s coefficients measured experimentally through laboratory testing is evaluated in Section 4.3.2. The effectiveness of the DEM-Voronoi model to replicate the effects of pore pressure on weakening the undrained strength of the rock according to Terzaghi’s effective stress law is given in Section 5.Finally,the calibrated model is utilised to simulate the formation of EDZ around an excavation scale numerical simulation of a drift in Canada’s Underground Research Laboratory(URL).The pore pressure evolution due to the re-distribution of in situ stresses and accompanied material fracturing is compared to the pore pressure measured at the tunnel site.

Fig.2. Effect of pore pressure on peak strength in triaxial compression tests at various confining stresses for Darley Dale sandstone with 21% porosity, showing approximate conformity with the Terzaghi’s effective stress law. Dashed lines are the constant“effective stresses” (after Murrell,1965).

2. Generation of the grain-based model

The Voronoi tessellation algorithm implemented in RS2 FEM software (RocScience Inc, 2015a) is used to generate a twodimensional (2D) rectangular intact rock sample with a height of 150 mm and width of 60 mm (Fig. 3a), according to the recommended height-to-width aspect ratio of 2.5 for a laboratory unconfined compressive strength (UCS) test (Ulusay and Hudson,2007). A circular sample with a radius of 80 mm is used to simulate a Brazilian indirect tension test,as shown in Fig.3b.The grain size is selected based on the average grain size of the LdB granite,which is 3.2 mm in diameter(Martin,1993).Hence,an average edge length of 2 mm is used to generate the Voronoi grain-based body of the model. Note that the generated Voronoi network used in this study is not based on the actual geometries of mineral grains obtained from digital images of the rock.Hence,the numerical models and the initial geometries used in this study are not trying to reproduce the exact rock structure observed in the actual LdB granite. In order to import the geometrical data of the Voronoi network into the UDEC software (Itasca, 2014), the embedded language FISH is used. The Voronoi tessellation algorithm implemented in UDEC and RS2 codes, as explained in Jing and Stephansson (2007), is presented as follows:

(1) Step 1 - Expand the tessellation area slightly over the boundary of the domain of interest to reduce the boundary effect.

(2) Step 2 - Distribute nodes along the domain boundaries according to the desired edge length of the elements.

(3) Step 3-Generate interior nodes randomly until the desired density of interior nodes is achieved.

(4) Step 4 - Smooth the position of nodes so that a more uniform grid of nodes can be generated by an iteration process(the more the iteration,the more uniform the grid).

(5) Step 5-Connect the interior nodes to form the intermediate triangular element mesh.

(6) Step 6 - Form the Voronoi polygonal elements by constructing perpendicular bisectors of all triangles which share a common edge.The polygonal elements are truncated at the boundaries of the initial tessellation region.

As discussed in Lan et al. (2010), there are generally three options available to simulate the grain structure:(i) polygonal grains commonly used in UDEC-Voronoi code (Itasca, 2014); (ii) diskshaped grains used in discrete element codes,such as PFC;and(iii)square-shaped elements used in some finite element codes,such as RFPA2D (Tang,1997). Although Voronoi models are computationally more expensive than particle-based and finite element squareshaped ones, from a geometrical point of view, it appears that the polygonal structures such as those used in UDEC-Voronoi code may be more representative of the mineral structure observed in crystalline rock(Lan et al.,2010).Additionally,owing to the full contact between grains and better interlocking offered by the random polygonal shapes,the Voronoi-shaped grain assemblage provides a closer resemblance in shape of the unit block to the actual lowporosity rocks such as LdB granite for hydro-mechanical modelling purposes.

Note that the numerical analyses presented in this study are carried out using 2D UDEC code. The selection of 2D over threedimensional (3D) analysis is due to the constraints imposed by computational capacity which is limited in simulating the hydromechanical behaviour of a large number of small (centimeter to millimeter range) 3D elements as it would be required to fully reproduce brittle failure processes.As a result,2D analysis is more efficient than 3D one in terms of model resolution (minimum element size) to replicate rock fracture and associated hydromechanical coupling processes.

In this study, the rock mechanics sign convention is used in which the compressive stress is positive, whereas the tensile stress is considered negative.The convention here is that the fluid pressure is always pointing toward the boundary of blocks or elements.

Fig. 3. Model geometry for (a) uniaxial compression and (b) Brazilian testing. Grains with purple, light red, light brown, and dark red represent potassium feldspar (Kfeldspar), plagioclase feldspar (P-feldspar), quartz, and biotite minerals, respectively.(c)Thin section of LdB granite(Lim et al.,2012).(d)Histogram showing the abundance of different mineral phases in real LdB granite and the synthetic rock mass (SRM)model.

The calibrations and simulations of these uniaxial compression and Brazilian tests are based on the short-term physical laboratory test properties of LdB granite sampled from the 240 m level at the URL in Pinawa, Manitoba, Canada (Martin, 1993). The grain-scale material heterogeneity is introduced into the model based on the mineral composition of the LdB granite.Four material grain types are defined,including potassium feldspar(K-feldspar), plagioclase feldspar(P-feldspar),quartz,and biotite grains,which are distributed in the model based on the abundance of each mineral type reported in Table 1.The histogram in Fig.3d compares the distribution of various mineral phases in both actual granite and the UCS model.

3. Cohesive crack model defining the behaviour of the grain contact

Experimental observations demonstrate that during fracturing of quasi-brittle materials such as concrete and rock, there is intermediate space between cracked and intact portions of the material(Labuz et al.,1985).The region formed at the tip of an opening crack(mode-I) is defined as the fracture process zone. The fracture process zone is a zone of partially damaged and interlocked material which is still able to withstand the applied stress and transfer the load from one surface to the other(Fig.4).The material outside the fracture process zone is assumed to be linear elastic(Dugdale,1960;Hoagland et al.,1973; Labuz et al.,1985). As the crack propagates,the microcracks in the fracture process zone merge and become a single structure to give continuity to the pre-existing crack. Actually, fracture process zone acts as a bridging zone between the cracked and intact regions.

In an attempt to idealise the effect of the fracture process zone on the fracturing of quasi-brittle materials, several cohesive crack models for both mode-I(crack opening)and mode-II(crack sliding)have been proposed. The most important cohesive crack models(for idealised mode-I fracturing) and slip-weakening models (for idealised mode-II) were reported in Dugdale (1960), Hillerborg et al. (1976) and Munjiza (2004), and in Ida (1972) and Lisjak et al. (2014), respectively.

In this study, the initiation and growth of the explicit microcracks are simulated using DEM software by means of breakage of the contacts between the polygonal grains. Thus, composite fracture trajectories are free to develop within the constraints imposed by the geometry of the Voronoi grain network. Depending on the local stress and relative contact wall displacements, the contact may yield under mode-I(tensile mode),mode-II(sliding mode)or mixed mode(mode I-II)conditions.A cohesive crack model shown in Fig.5 defines the behaviour of the grain interfaces under variousloading conditions. In this model, the stress-displacement curves for mode-I and mode-II conditions have a nonlinear pre-peak branch and a softening branch (Fig. 5). The nonlinear pre-yield branch is intended to represent the decay of stiffness in the preyield state due to the progression of damage development.

Table 1 Micro-material properties for different minerals of LdB granite.

According to this cohesive model, the contact boundary between grains yields when the relative contact displacement reaches critical values of crack opening, up(mode-I), and shear displacement,sp(mode-II).Displacements at peak strength are evaluated as

Fig. 4. Phenomenological description of a quasi-brittle fracturing process in (a) mode-I (modified from Huespe and Oliver, 2011) and (b) mode-II fracturing.

where e=exp(1)=2.718281 is the base of the natural logarithm;kn0,cand ksh0,crepresent the initial normal and shear stiffnesses of the contacts, respectively; and Tcand fshare the strengths of the contact in tension (mode-I) and shearing (mode-II), respectively.

The peak shear strength of the grain contact(fsh)is determined based on the following bilinear criterion (as demonstrated in Fig. 6):

where ccand cc，2ndare the primary and secondary cohesions of the contact, respectively; φcand φc，2ndare the primary and secondary friction angles of the contact, respectively; σnis the normal stress acting across the contact length;and σn，Transitionalpointis the normal stress at which the shear failure envelope transitions from a steeper to a shallower slope.

The application of the bilinear shear strength envelope is supported by studies on the influence of confining pressure on the mode-II fracture toughness(KII)and shear strength of various intact rocks(Lundborg,1968;Backers,2005).For example,Backers(2005)observed that the determined KIIincreases in a bilinear fashion with increasing confining stress and shows a change in slope at about 30 MPa for the tested rock types (Fig. 7). For the case of the ?SPO diorite with approximately similar mechanical properties(UCS=219 MPa,Ti=15 MPa,Ei=68 GPa,and νi=0.24) with LdB granite, the transition from a steep to a shallow slope occurs at approximately 30 MPa with the friction angle of φc= 70°for σn<30 MPa and φc,2nd=18°for σn>30 MPa.The estimated friction angles for these two regimes suggest a ratio of φc/φc,2ndequal to approximately 4.

In the post-yield state (softening branch), the following stressdisplacement (traction-separation) relationship that links the cohesive stress transmitted by the contact to the contact displacement in tensile and shear states is used:

where τ is the shear stresses acting across the contact length, and χ（Di） is the softening function which defines the decay of contact strength in the post-yield regions.

Fig.6. Bilinear failure criterion with residual strength and tension cutoff to define the strength properties of the contacts in shear and tension (Farahmand, 2017). The bilinear peak shear strength envelope is approximated based on the experimental data obtained from the shear strength of the punch-through rock samples in triaxial compression as reported by Backers (2005). The purple bilinear envelope defines the peak shear strength of the contact, while the red envelope defines the residual strength of the contact after yielding.

Fig. 7. Influence of confining stress on shear strength of the punch-through shear samples. Two regimes of different slopes are separated, and a transitional zone between the two regimes is indicated.Experimental data are taken from Backers(2005).

The shape of the softening branch of the stress-displacement curves in mode-I and mode-II fracturing (see Fig. 5a and b) can be approximated using the following empirical softening equation proposed by Evans and Marathe (1968):

where A,B,and C are the empirical curve fitting parameters equal to 0.63,1.8,and 6,respectively;and Di(i=I,II and I-II)is the damage variable with a value between 0 and 1. This equation defines the shape of the softening branch of the stress-displacement curves in both mode-I and mode-II fracturing.Eq.(6)is derived based on the curve fitting to the experimental stress-strain data obtained from concrete in direct tension (Munjiza, 2004).

The damage variable for contacts subjected to mode-I(Munjiza,2004), mode-II (Tatone, 2014), and mode I-II (Tatone, 2014) displacements are respectively defined as follows:

where u and s are the relative displacements in normal and shear directions of contact,respectively;and uresand sresare the residual opening and residual slip, respectively, which correspond to the relative displacements that a cohesive contact requires to reach the residual normal and shear bonding stresses, respectively.Substituting the damage variables into Eq. (6), the shape of the softening branch of the stress-displacement curves in both mode-I and mode-II can be obtained.

The values of residual opening (ures) and residual slip (sres) in Eqs. (7)-(9) can be calculated by

where GIcand GIIcare defined as the energies required to extend the crack surface of a unit area in tension and shear modes, respectively.GIcin the plane strain condition can be estimated according to the following equation proposed by Griffith (1921):

In this study, the values of GIIcare considered as 2GIc(Lisjak et al., 2014; Farahmand and Diederichs, 2015). The behaviour of the contacts under normal compressive loading is represented by the hyperbolic function suggested by Bandis et al. (1983):

where u0is the contact closure under compression,and umaxis the maximum possible contact closure.

To determine if the cohesive model is able to replicate the defined stress-displacement behaviours of the contacts under tensile and compressive loadings, a simple numerical test is performed,as shown in Fig.8.Here,two blocks,each measuring 30 cm in width and 10 cm in height,represent the quartz material and are glued together by their cohesive contacts. At the first stage, the upper block is subjected to the tensile forces acting on the top external boundary (Fig. 8, top right) to induce tension in the cohesive contact while the bottom external boundary of the lower block is fixed in x-and y-direction in zero displacement conditions(the red arrows in Fig.8,left,illustrate the stress paths for a contact in the first tension cycle). Following the breakage of the contacts,compressive forces are imposed on the upper block to force the contacts to close (stress paths illustrated by the blue arrows in Fig. 8, left). This loading condition aims to validate if the model is able to replicate the normal stress-normal displacement behaviour of the contact under compressive loading as defined by Eq. (13).Finally,tensile stresses are applied to the top external boundary of the upper block again to check the stress-displacement response of the contact in tension (yellow arrows in Fig. 8, left). The result of this numerical test validates the successful implementation of the constitutive model into the UDEC software using FISH commands.

Fig. 8. Normal stress versus normal displacement of contact for a series of pull and push tests, which is used to validate the implementation of the cohesive constitutive contact model in UDEC for this study. Tensile stresses have a positive sign, in accordance with the general sign convention for internal stresses in UDEC, and positive normal displacement indicates extension (opening).

4. Model calibration process

In continuum models,the input parameters are directly derived from the hydro-mechanical properties measured using standard laboratory tests.For the case of GBMs,however,the input values for the micro-components of the model cannot be determined directly from the macro-properties measured from physical laboratory tests. Instead, the values of the micro-mechanical parameters assigned to the grains and their interfaces must be determined through a numerical calibration process, in which the emergent macro-properties of the model are compared to the relevant measured response of the material (Potyondy and Cundall, 2004;Tatone, 2014).

In this section, a new systematic calibration procedure for obtaining the hydro-mechanical input parameters of the DEMVoronoi model is developed. The calibration process involves two main steps, as illustrated in Fig. 9. In the first step, the micromechanical deformation and strength properties of the grains and interfaces are required to be calibrated to reproduce the exact macro-mechanical properties tested in the laboratory. In the next step,the values of initial and residual apertures are determined by a series of numerical experiments, thus the model replicates the permeability of the material at various amounts of deformation and damage. Finally, the initial and residual contact apertures are refined until the model reproduces poroelastic properties similar to those derived by laboratory tests. The poroelastic properties that must be calibrated in the model are the Biot’s coefficient, α, and Skempton’s coefficient, B. The step-by-step procedures for calibration of the hydro-mechanical properties are discussed in the following sections.

Note that the main goal of the calibration procedure presented here is to find out a set of input parameters that simultaneously yield the correct stress-strain behaviour for triaxial and Brazilian simulations and also to provide a representative permeability and pore pressure response of the intact rock under study. However,through an unprecedentedly large number of parameter sensitivity analyses, it was determined that a range of parameter values yielded the target hydro-mechanical behaviour.To constrain this overdetermined problem, additional information regarding typical damage patterns, failure modes, permeability evolution with damage progression and poromechanical response in tested samples needed to be considered. After doing so, a narrow range of potential input values could be defined which yielded the desired stress-strain behaviours with emergent damage patterns, permeability and pore pressure response that agreed with laboratory observations (Tatone,2014).

4.1. Calibration of the model to obtain the mechanical behaviour of the rock

To calibrate parameters controlling the deformability and strength properties of the DEM-Voronoi model,a series of uniaxial compression and Brazilian tension tests is simulated in UDEC. The laboratory experimental data of LdB granite are used as a case to calibrate DEM-Voronoi parameters in UDEC, as presented in Table 2.In the following subsections,the detailed procedure used in each step of the mechanical parameter calibration process and relationships between the model micro-parameters and macroproperties of intact rock samples are discussed.

Fig. 9. Flowchart illustrating the micro-mechanical and hydraulic property calibration procedure for grain-based models.

Table 2 Experimental mechanical properties of LdB granite (laboratory experimental data are from Martin (1993)).

4.1.1. Selection of deformability parameters

(1) Finding a representative grain size

The modelled grain size should be of the same order as the average grain size in the real rock fabric.In the case of LdB granite,an average edge length of 1 mm is selected for the Voronoi grains to result in an approximately 2.5 mm diameter grain size.However,in the case of medium-and large-scale problems such as underground excavations, satisfying this condition requires impractical computational power to handle a huge number of grains, which is currently unrealistic for geotechnical modelling practitioners to achieve (Gao and Stead, 2014). As a rule of thumb, modelled grain size can be selected for calibration of large-scale models following the ISRM suggestion (Ulusay and Hudson,2007) for a UCS test,for which the grain size should be at least 10-20 times smaller than the width of the sample. In addition, the concept of characteristic length of material,Ich,can be used as the maximum allowable grain size for a proper fracture process zone representation in a quasibrittle rock (Hillerborg et al.,1976):

where Eiis the Young’s modulus of the intact rock.

(2) Calibration of contact stiffness in normal and shear directions

The values of contact normal stiffness(kn0,c)and shear stiffness(ksh0,c) control the macroscopic Young’s modulus (Ei), Poisson’s ratio(νi),crack initiation threshold,and brittle to ductile transition of the sample in the post-failure region. Diederichs (2000) and Kazerani et al. (2012) suggested that the ratio of the contact shear stiffness to normal stiffness(ksh0,c/kn0,c)should be equal to the ratio of Gi/Ei,where Giis the shear modulus of the undamaged rock, i.e.ksh0,c/kn0,c= Gi/Ei. This ratio is normally between 0.35 and 0.5.

Kazerani et al. (2012) also suggested closed-form expressions based on the width of the cohesive zone to calculate bond stiffnesses in plane strain problems such that:

where β is a constant multiplier, and Tcis the contact tensile strength.

However, the following relationship as suggested in the UDEC documentation(Itasca,2014)is used in this study to determine the normal stiffness of the contacts:

where n is a multiplier factor,Kiis the bulk modulus of intact rock,and ΔZminis the smallest width of an adjoining zone in the normal direction. The shear stiffness of the contacts can be determined according to the ratio of the contact shear to normal stiffness(ksh0,c/kn0,c),which is equal to the ratio of Gi/Ei(Diederichs,2000;Kazerani et al., 2012).

(3) Calibration of macroscopic Young’s modulus (Ei)

For calibrating the macroscopic Young’s modulus of the sample,the density(ρgr),Young’s modulus(Egr),and Poisson’s ratio(νgr)of all grain mineral phases are assigned based on elastic properties of each individual mineral phase obtained from laboratory test data.In this calibration example,the elastic values for all types of mineral grains are given in Table 1.After assigning the grain properties,the value of kn0,cshould be adjusted first until the macroscopic Young’s modulus of the modelled uniaxial compression sample matches the target laboratory test value.It should be noted that decreasing the ksh0,c/kn0,cratio increases the macroscopic Young’s modulus and also decreases the ductility of the material.

(4) Calibration of macroscopic Poisson’s ratio (νi)

Diederichs (2000) showed that the macroscopic Poisson’s ratio of material in the bonded particle models (BPMs) is related to the ratio of contact shear to normal stiffness (ksh0,c/kn0,c). As this stiffness ratio decreases, Poisson’s ratio increases, and the sample becomes more dilatant.Results reported by Kazerani et al.(2012)and Gao and Stead(2014)support this idea for the case of DEM-Voronoi models. Thus, the contact stiffness ratio is calibrated to fit the material’s Poisson’s ratio.

4.1.2. Selection of strength parameters

In this section, the procedure for calibration of microscopic strength parameters, including contact tensile strength, cohesion,and friction angle is discussed.

(1) Calibration of contact tensile strength (Ti)

Tensile strength assigned to contacts significantly affects the macroscopic tensile strength and crack initiation stress of intact rock. As a result, accurate calibration of this parameter is of paramount importance for reproducing the correct crack initiation threshold and tensile strength of rock. The values for the tensile strength of different mineral interfaces in this study are taken from the results of laboratory testing on different mineral samples reported in Savanick and Johnson (1974). As the Brazilian and UCS tests result in lower macroscopic tensile strength and crack initiation threshold reported for the LdB granite, the contact tensile strength for each mineral phase is increased until the correct Brazilian indirect tensile strength and crack initiation stress are obtained.The reason for this discrepancy is attributed to the source of tensile strength data which are from laboratory experimentation on Rockville granite (Savanick and Johnson, 1974). As a result, some modifications are needed to be made to the given values of contact properties to match the macroscopic behaviour of LdB granite.

In order to calibrate the tensile strength of the models,a suite of indirect tensile(Brazilian)and UCS physical laboratory tests should be performed to provide necessary corrections to the values of contact tensile strength to obtain the target values of macroscopic tensile strength and crack initiation stress.

(2) Calibration of contact cohesion (ci)

Laboratory and in situ observations reveal that the process of rock damage consists of two dominant mechanisms, i.e. extensile cracking and shear cracking,where the extensile cracking is known as a primary mechanism of damage and the shear cracking becomes important only after a sufficient number of extensile cracks(after crack damage threshold is reached) extend through the sample and start to interact(Diederichs,2007).In the DEM models,the contact tensile strength and cohesion define the tensile(mode-I) and shear (mode-II) behaviours of the cracks, respectively. As a result,applying a high ratio of contact cohesion to tensile strength(cc/Tc)would force the bonds to dominantly yield in a tensile state.

The Griffith criterion(Griffith,1921)predicts a cohesion value of two times the tensile strength.Fracture studies by Okubo and Fukui(1996) and Laqueche et al. (1986) show that for very small cracks(smaller than grain size), the ratio of crack cohesion to tensile strength is approximately equal to 4 (Diederichs, 2000).

For calibration purposes, the initial value of contact cohesion should be selected in order to reproduce the target peak UCS.In this study, the initial cohesion value for each type of contact is set to approximately four times the contact tensile strength (cc/Tc= 4)based on the results reported in Okubo and Fukui (1996) and Laqueche et al. (1986). Afterwards, a series of 2D biaxial compression test models with varying confining stresses is executed to improve the accuracy of the macroscopic results.For each series of simulated tests, the confining stress and the corresponding strength values are required to be plotted in σ1-σ3space to calculate the emergent material cohesion (ci) and friction angle (φi). In the present study, after a series of iterative tests, cc/Tc= 2.7 was determined to be suitable to replicate the desired macroscopic contact cohesion properties.

(3) Calibration of contact friction angle (φi)

As mentioned previously, the shear failure mechanism becomes important only after sufficient tensile cracks accumulate and start to interact.Thus,it is expected that the contact friction angle affects the behaviour of rock at applied axial stresses greater than the crack damage threshold of the material.The contact residual friction angle(φres,c) is considered to be the most important parameter that controls the post-peak behaviour of the material, and it can also affect the degree of the brittleness of the rock. However, determining the micro-frictional properties of the grains from experimental data is difficult.As a result,the friction angles of contacts must be calibrated by conducting a series of 2D biaxial compression model tests with different confining stresses. For the LdB granite case, the most appropriate value of the microscopic contact friction angle (φc) is approximately equal to the macroscopic friction angle (φi) of the rock.

4.2. Results of the calibration to obtain macroscopic mechanical properties

Since the models are composed of four different mineral phases,ten contact types are required to be calibrated. Calibrated contact properties for the defined interfaces of the model are listed in Table 3. It should be noted that the value for the fracture energy release rate (GIc) for each contact type can be estimated using Eq.(12).The fracture toughness(KIc)values listed in Table 1 are used to calculate GIcfor most interfaces. However, there are no available data for GIcvalues of the interface between quartz and feldspar,and the interface between biotite and other minerals. Therefore,following Mahabadi et al.(2012),the fracture energy release rate of the contact is assumed to be proportional to its tensile strength(i.e.This assumption is based on the fact that both GIcand KIcof the contact are directly proportional to their corresponding values of Tc.

Since the solution algorithm of the UDEC software is dynamic and based on a time step solution, the loading rate applied to the platens must be defined.High loading rates or platen velocities(e.g.>0.1 m/s) result in numerical oscillations within the sample;therefore, it is important that the velocity of converging platens must be sufficiently slow and also applied damping ought to be high enough to ensure that the sample remains in quasi-static equilibrium and perturbations are dispersed throughout the sample before new loads and displacements are applied(Kazerani et al.,2012). In general, the effect of loading rate is most pronounced on the post-peak response.Faster platen velocities result in less brittle behaviour,while slower loading results in localised rupture zones.This is expected since, for slower loading rates, fewer contacts break at each time step, which allows for load re-distribution and consequent concentration of subsequent damage around previous contact rupture points (Diederichs, 2000). For faster loading rates,more distributed cracks can occur (in a heterogeneous material with a given strength distribution)in each time step,resulting in a more distributed rupture network (Diederichs, 2000).

The axial stress versus axial strain behaviour is monitored during testing by FISH functions.Compressive stresses are measured in the upper platen. The axial strain is calculated based on the displacement of the two points at the very top and bottom of the samples for both cases of compression and Brazilian tension tests.In addition,14 sets of measuring points along the external lateral sides of the compression sample are selected to monitor the lateral strain,and the average of the lateral strain of the monitoring points represents the overall lateral strain of the rock under compression.In Brazilian tests,the tension is calculated from the forces acting on the upper platen according to the following equation:

Table 3 Calibrated model parameters for LdB granite.

where F is the load at failure,b is the thickness of the disk,and d is the diameter of the disk sample.

The central straight-through crack Brazilian disk (CSCBD) samples (Fowell,1995) are used to measure the emergent mode-I and mode-II fracture toughness values of the Brazilian tensile strength models. As illustrated in Fig.10, the CSCBD samples are 80 mm in diameter,and have a 0.024 mm long vertical notch or a 0.024 mm long notch with a 27.2°inclination to the loading direction, for mode-I and mode-II fracture toughness tests, respectively. The geometries of the CSCBD samples are chosen according to Fowell(1995) (Fig. 10b). The Brazilian tensile strength (Ti) is calculated by Eq.(17),and mode-I(KIc)and mode-II(KIIc)fracture toughnesses are deduced as follows (Atkinson et al.,1982):

where R is the radius of the Brazilian disk,and anis the half-length of the notch. Since the analysis performed herein is in two dimensions, b is assumed to be 1 m. The two dimensionless coefficients NIand NIIin Eqs.(18)and(19)can be calculated according to the following equations:

where θ is the notch orientation with respect to the loading direction.It should be noted that when θ=0°and an/R=0.3,the pure mode-I condition will be achieved,as NI=1 and NII= 0.However,for the sample with the notch oriented at θ=27.2°and an/R=0.3,the sample will fail in pure mode-II condition. In this case, the values of NIand NIIare equal to approximately 0 and 1.72342,respectively(Fowell,1995).

Figs. 10a and 11a show the resultant stress-strain curves for Brazilian and unconfined compression tests, respectively. The Brazilian tensile strength Ti= 10.8 MPa is achieved for the calibrated model, while KIc,i= 1.44 MPa m1/2, and KIIc,i= 3.5 MPa m1/2.

The axial stress-axial strain response of the model in compression deviates from linearity at σ1=181 MPa(87%of UCS),which corresponds to the reversal point in the volumetric strainaxial strain curve (Fig. 11). The crack damage threshold corresponds to the stress at which this reversal occurs.The onset of crack initiation occurs at 86 MPa,which is determined based on the onset of tensile cracking in the model (red lines in Fig. 11a). The crack initiation threshold occurs approximately at 41% of UCS. These results are comparable to the experimental results, as listed in Table 2,in which the crack initiation and damage occur at 43%and 83% of UCS, respectively.

The number of induced tensile and shear cracks demonstrates that the primary mechanism of damage in the UCS model is extensile cracking, while the shear (cohesive) cracking starts to become the dominant behaviour only after the axial stress applied to the sample exceeds the crack damage threshold. A summary of the experimental and simulated results is given in Table 4.Strength properties obtained from the model results demonstrate a very good agreement with the experimental data reported by Martin(1993).

Fig.12 shows the comparison between the crack initiation and damage thresholds and peak strength obtained for LdB granite from experiments and simulations.The results demonstrate a good match between experimental and numerical simulation findings.

Fig.10. (a)Peak stress-axial strain curves obtained from calibrated Brazilian and CSCBD samples for measuring mode-I and mode-II fracture toughnesses.Damage in the samples is shown. (b) Geometrical configuration of the CSCBD samples as suggested by Fowell (1995).

Fig.11. Model response during uniaxial compression test:(a) Axial stress-axial strain curves and associated incremental accumulation of tensile (extension) cracks and shear (cohesive) cracks; and (b) Volumetric strain-axial strain curve.

4.3. Calibration of the model to determine the hydraulic and poromechanical behaviours of the rock

In this section, the procedure to calibrate the hydraulic parameters of biaxial test simulations is discussed. In addition, the influence of contact apertures on the emergent poroelastic properties of the model is examined.

The fluid flow simulation scheme in UDEC is restricted to modelling the flow through joint elements, but not in the intact material grains. As a result, in the context of the GBM, the fluid is only able to flow within the joints between adjacent mineral grains.In UDEC,the numerical implementation for fluid flow is based on a“domain” structure. The domain is a 2D area defined to represent the void space created at joint intersections by contacts around that spot and the void space along the joints (Jing and Stephansson,2007). A joint formed between two blocks is separated into a number of domains by the contact points.The location of a contact point is determined based on the internal discretisation by the finite difference mesh elements.The domain is assumed to be filled with fluid at uniform pressure and it communicates with its neighbours through the contacts (Itasca, 2014).

It is assumed that the fluid flow occurring in the joint space can be generally simplified into flow between two parallel plates with smooth surfaces, where the pressure gradient between adjacent joints governs the fluid flow.The flow rate for each contact can then be determined based on the cubic law(Witherspoon et al.,1980)as follows:

where μ is the dynamic viscosity of water, a is the contact (joint)aperture, Lcis the assigned contact length between two adjacent domains, and ΔPdis the pressure difference between the adjacent domains.The dynamic viscosity of the water is equal to 0.0001 Pa s at about 25°C (Detournay and Cheng,1993).

Table 4 Comparison between experimental and simulation results for the LdB granite (experimental data are from Martin (1993)).

Fig. 12. Strength envelopes from experiments and simulations (maximum principal stress, σ1, versus minimum principal stress, σ3).

In UDEC,the deformation of material changes the fluid pressure in joints by changing the joint aperture. As a result, the state of stress in the adjacent grains will change according to the magnitude of the fluid pressure acting on them. This simulates the twoway hydro-mechanical coupling effects between pore pressure and rock deformation. During the development of cracks along grain boundaries, the aperture evolution as a result of opening(mode-I)and shear dilation(mode-II)is recorded,and according to the concurrent magnitude of the aperture, the crack fluid rate is calculated.

Note that the UDEC implements a transient hydro-mechanical interactive analysis algorithm using the dynamic relation technique based on the domain structure(Jing and Stephansson,2007).The algorithm is implemented by performing a series of iterations(termed cycles in the UDEC) for fluid flow with the time step defined by the user but controlled by a criterion to maintain numerical stability of the explicit flow analysis.For a typical cycle i,the following tasks are carried out in the given order:

(1) Computing the flow rate Qi/i+1from domain i to the domain i+1 as a function of unbalanced fluid pressures Piand Pi+1in each domain.

(2) Computing the initial volume of each domain.

(3) Performing a series of mechanical relaxation steps until a continuous flow in each domain is reached. Assuming incompressibility of the fluid,the net volume of fluid flowing into a domain during a hydraulic time step must equal the total volume change of the domain. The unbalanced fluid volume being the difference between the two is gradually reduced during the relaxation process.

(4) The domain pressure is changed in accordance with the unbalanced volume of each domain following an adaptive change of domain volume until a stable pressure is obtained(Jing and Stephansson, 2007).

4.3.1. Permeability calibration

It is well documented in geomechanics literature that the permeability of rock is very dependent on the stress condition(Zoback and Byerlee, 1975; Souley et al., 2001; Mitchell and Faulkner, 2008; Popp et al., 2008; Jiang et al., 2010). When rock is subjected to mechanical loading, and therefore changes in stress condition, the rock microstructure closes, opens, extends and induces new cracks, which in turn affect the mechanical and hydraulic properties of the rock. As a result, the stress-dependent permeability is a function of two fundamental mechanisms:

(i) microcracking-induced permeability change which is related to the growth of the microcracks, and

(ii) stress-induced pore volume change where the pore volume is a function of the interaction between fluid pressure and mechanical stress.

In the first case, damage-modified permeability evolution is related to the heterogeneous nature of rock material and the presence of internal microstructure with different sizes and properties. Thus, when the rock is subjected to a mechanical compressive load, local tensile stress concentrations induced on these internal microstructures initiates a complex process of fracture propagation in the rock. Generation of these crack corridors(channels) modifies the overall permeability of the rock to an extent that is dependent on the opening aperture, density, and connectivity of such cracks.

In this section,a series of biaxial compression test simulations is presented. Throughout the loading stages of the simulations, the permeability evolution in the model as a result of elastic deformation and microcracking are monitored. For each test, the initial and residual apertures are modified manually in order to find a set of values for contact apertures that results in a stress-permeability curve similar to those obtained experimentally. It should be noted that in the DEM-Voronoi model, it is the assigned initial and residual apertures that control the overall hydraulic properties of the samples. Hence, to calibrate the model, the values for the initial aperture of the contacts are adjusted such that the permeability of the sample matches that of the real rock under a specific hydrostatic loading condition. Then, the values of the residual apertures are found by performing a series of biaxial hydro-mechanical tests such that the model permeability in the elastic portion of the axial stress-strain curve matches the permeability of LdB granite at the corresponding stress state.

During the biaxial test, fluid flow through the rock body is simulated by applying fixed hydraulic heads to the top and bottom of the model while the two lateral boundaries of the model are considered to be impermeable.Next,the mass continuity equation is solved at joint (domains) intersections through an iterative scheme under the prescribed hydraulic conditions to obtain the flow field.During the test,the sample should be loaded in a drained loading condition. To do so, the fluid is assumed to be incompressible; therefore, the deformation of grains and change of the joint aperture will not be able to induce any pressure change in the fluid.

The axial permeability of the compression model (kyy) is obtained from Darcy’s law for flow in anisotropic and homogeneous porous media (Bear,1972) as follows:

where Qyis the flow rate in y-direction,A is the cross-sectional area of the compression model,ΔLsis the height of the rock sample,and ΔPsis the fluid pressure difference between the top and bottom of the rock sample.

The sensitivity of the emergent permeability of a biaxial sample to the assigned values of the initial and residual apertures is shown in Fig. 13, when the sample is subjected to a confining stress of 10 MPa, and the fluid pressure difference between the top and bottom of the rock sample is set to 5 MPa. When the values of the initial and residual apertures are set to 0.2 μm and 0.05 μm,respectively, the resultant stress-permeability curve mimics the behaviour obtained from a laboratory-tested sample. The experimental results for the permeability of LdB granite are taken from Souley et al.(2001).Note that under compressive loading,a crack or fracture will never be closed even under extremely high normal stresses due to the presence of contact asperities on the surfaces of the discontinuity which provide channel-like connected void spaces between asperities,therefore leading to a residual hydraulic aperture.This is the reason for the small but undiminished residual flow observed under high normal stresses(Bandis et al.,1983;Jing and Stephansson, 2007). In the fluid flow simulation scheme in UDEC, a minimum value is assumed for hydraulic aperture (ares),below which mechanical closure does not affect the contact permeability.

As illustrated in Fig.14, a permeability increase from its initial intact value by approximately six orders of magnitude is obtained for the model under the confining pressure of 10 MPa, which is consistent with the reported experimental data. The reduction in the permeability of the rock at the early stage of the deviatoric loading is achieved by assigning a lower value to the residual aperture than the initial aperture. As a result, as the deviatoric stress increases in the model, the hydraulic conductivity of the cracks decreases due to their closure under compression.After this stage, the axial permeability remains almost constant because the failed cracks in the models are mostly isolated and do not contribute to the fluid conduction in the model.However,as failed cracks start to accumulate, the axial permeability gradually increases by more than three orders of magnitude.Through the postpeak region, a continuous increase in the number of connected damaged microcracks leads to a dramatic increase in hydraulic conductivity by approximately six orders of magnitude from its initial non-damaged permeability. In summary, the permeability change can be sorted into three different stages.At the first stage(at stresses less than the crack initiation threshold), the permeability reduction is caused by the closure of microstructures during early stages of deviatoric stress.In this range of deviatoric stress,no crack nucleation or crack propagation occurs. However, when the stress exceeds the crack initiation stress, microstructures start to propagate.The stable crack growth in the deviatoric stress state between crack initiation and damage thresholds does not affect the overall permeability of the sample.By increasing the axial stress above the crack damage threshold, the unstable crack propagation initiates which in return results in the dramatic increase in permeability.The reason for the red curve in Fig.13 returning to lower abscissa is that the rock loses its bearing capacity(post-peak stress drop)after reaching the yield stress, and therefore, the deviatoric stress induced in the sample is gradually reduced until a residual strength is reached. This gradual reduction in the post-peak stress is accompanied by the increase in the sample permeability due to post-peak dilation of the sample(i.e.formation of new cracks,and shear dilatancy of the induced cracks).

Fig. 13. Permeability evolution of the biaxial samples under the confining stress of 10 MPa when various sets of initial and residual apertures(a0 and ares)are assigned to the contacts. The experimental results for the permeability of LdB granite are taken from Souley et al. (2001).

4.3.2. Calibration of poromechanical parameters

This section focuses on the calibration of Biot’s coefficient (α),bulk moduli (K, Ksand Ku), and Skempton’s coefficient (B) for the UDEC models of LdB granite biaxial test simulations. Hydrostatic compression tests are simulated in order to determine the compressibility of the granite and its poroelastic parameters.Three hydraulic boundary conditions,as shown in Fig.15,are considered.While the applied confining stress deforms the sample,the changes in sample volume (ΔV) and pore pressure (ΔP) are monitored during the simulations (Fig.16).

In the first test simulation,the hydrostatic compressive stresses are applied to the four external boundaries of a square sample in an undrained condition(Fig.15a).The undrained loading characterises the condition where the fluid is trapped in the porous solid such that the variation of fluid content,ξ,in the sample is equal to 0.The goal is to calculate the emergent Skempton’s coefficient,B,as well as the undrained bulk modulus (Ku) of the model. Under such testing conditions,any variation of confining stress(Δσ)would lead to volumetric strain (Δεv) as

It is apparent that pore pressure (P) will be induced under the undrained condition when ξ=0.Skempton’s coefficient(B)can be calculated as (Detournay and Cheng,1993):

Fig.14. Comparison between the laboratory testing (Souley et al., 2001) and numerical simulation results from this study for a biaxial sample of LdB granite under the confining stress of 10 MPa.

In the second test,the saturated sample is gradually loaded.The hydrostatic stress is increased in a drained condition while maintaining constant pore pressure(ΔP=0).It should be noted that this testing condition is known as“drained loading condition”,since the confining stress is increased by Δσ, but the pore pressure (P) is maintained at its initial value (P0) (Fig.15b). As explained by Lau and Chandler (2004), the compressibility modulus of the model in the drained condition can be derived as

Fig.15. Sketch of(a)undrained,(b)drained,and(c)unjacketed boundary conditions with principal compressive stresses(σ1,σ2 and σ3)and pore fluid pressure P(Makhnenko and Labuz, 2016).

In the third test, the bulk modulus of the solid matrix (Ks) is calculated according to Eq. (27), when the hydrostatic stress and pore pressure are simultaneously increased with the same increment (Δσ =ΔP) (Fig.15c):

Fig. 16. Pore pressure-confining stress response of the grain-based models under undrained loading condition: (a) Effect of assigned apertures to the contacts on the induced pore pressure in the samples; and (b) Effect of changing contact initial aperture on the emergent Skempton’s coefficient B of the models. The experimental data for the Skempton’s coefficient of LdB granite are taken from Lau and Chandler(2004).

Fig. 16 shows the change in the calculated Skempton’s coefficient (B) when different values of apertures are assigned to the contacts formed between grains. As expected, the pore pressure induced by the confining stress increases with decreasing contact apertures. Skempton’s coefficient of LdB granite measured in the laboratory(Lau and Chandler,2004)is used to calibrate the contact apertures. The bulk modulus of water is assumed to be 2.2 GPa(Detournay and Cheng,1993).

As shown in Fig.16a, the model with initial and residual apertures of a0= 0.2 μm and ares= 0.05 μm, respectively, yields approximately the same value of B as determined in the laboratory.The results of this sensitivity analysis,as shown in Fig.16b,suggest the following exponential relationship between the contact aperture and the resultant Skempton’s coefficient:

In Fig. 17, the derived hydrostatic stress-volumetric strain response of the models under the undrained hydrostatic loading is compared to that of the LdB granite observed in the laboratory.Similar to the induced pore pressure response,the calibrated model with the values of a0= 0.2 μm and ares= 0.05 μm successfully reproduces the rock response that was measured experimentally.

The estimated values of K and Ksobtained from numerical modelling(Fig.18a and b)are used to calculate the emergent Biot’s coefficient,α,of the GBM.Biot’s coefficient(α)is defined as the ratio of the volume change of fluid fill porosity to the volume change of the solid material when the fluid is free to move out of the solid material (Detournay and Cheng, 1993). This parameter can be calculated as follows (Biot,1941):

As mentioned previously, the capability of the calibrated DEMVoronoi model to replicate the impact of pore pressure change on the effective stress, according to Biot’s theory, is essential for simulating the hydro-mechanical behaviour of large-scale rock. It should be noted that during mechanical loading, in order to maintain a constant value for the fluid pressure through the model,the fluid should be assumed to be incompressible. As a result, the fluid pressure does not change with the grain deformation,whereas the fluid pressure changes the state of stress acting on the grains. Hence, a drained loading condition will be simulated.Calculating α using Eq. (29) shows that, while laboratory tests on core samples by Lau and Chandler(2004)indicate α=0.73,a value of α = 0.69 is estimated for the calibrated model. Table 5 lists the calculated values for compressibility moduli and Biot’s coefficient of the LdB granite.

The consistency of the emergent poroelastic parameters of the model can be checked (Detournay and Berchenko, 2001) by calculating K using the calculated values of α, B and Kuand the following relation (Fig.17):

The K value of 13.36 GPa calculated using Eq. (30) is relatively consistent with the value of 16.72 GPa obtained from numerical modelling (Fig.18).

Table 5 Comparison between the poroelastic parameters obtained from the numerical simulations and experimental data (Lau and Chandler,2004).

5. Influence of pore pressure on the strength thresholds of the rock

As mentioned in Section 1, experimental data show that the failure strength of a wide variety of rocks is a unique function of σ3- P (Byerlee, 1975; Gowd and Rummel, 1997; Paterson and Wong, 2005). As a result, to a fairly good approximation, Terzaghi’s effective stress law （σ′3= σ3-P） (see Fig. 19a) governs the failure strength of porous rocks (Serdengecti and Boozer, 1961;Handin et al., 1963; Murrell, 1965; Dunn et al., 1973; Byerlee,1975; Gowd and Rummel, 1997; Paterson and Wong, 2005).That is, the differential stress for fracture strength of a particular rock is approximately the same at the same “effective confining pressure” when the latter is taken to be the total confining pressure minus the pore pressure (see Fig. 19b). In addition, according to a series of triaxial compression tests,Brace and Martin(1968) concluded that the law of effective stress holds for peak strength of crystalline rocks of low porosity (0.1%-3%), such as Westerly granite, as long as loading rates are kept below certain critical values in order to maintain the drained loading condition during the test. For example,the critical strain rate was about 0.1 s-1for cylindrical samples of Westerly granite with the length of several centimetres (Brace and Martin,1968). As stated by Jaeger et al. (2009), “most experiments support the conclusion that the Biot’s coefficient (α) for failure of rocks is unity. And that, the Biot’s coefficient for failure processes has no particular connection to the effective stress coefficient (Biot’s coefficient) that appears in the theory of linear poroelasticity”. As a result, the concept of the effective stress coefficient in which the effective stress is equal to the difference between the total stress (σ3) and effective pore pressure (αP) only applies to linear elastic deformation of the porous material,but not to the nonlinear part of the axial stressaxial strain curve. However, the peak strength of materials is governed by Terzaghi’s effective stress law in which Biot’s coefficient is considered to be unity.

In this section, the impacts of pore fluid pressure on the emergent crack initiation and damage thresholds and peak strength of the models under the drained loading condition are explored. The question of interest here is to what extent the DEM-Voronoi model is able to reproduce the effect of pore pressure on weakening the strength of the rock. For this purpose, a series of biaxial compression tests is simulated to evaluate if the crack damage stress and peak strength of the samples follow Terzaghi’s effective stress law.To do so, saturated samples are subjected to various confining stresses when different pore pressures are initialised along the grain interfaces. During each test, the axial load is incrementally applied to the top and bottom of the sample, while the pore pressure is maintained at a constant value to create a drained loading condition.

The peak strength(σy),crack damage stress and crack initiation strength obtained by applying different sets of pore pressures and confining stresses (σ3) to the biaxial samples are plotted in Fig.20a-c.The points with the same colours in the plots show the strength of the samples that experience the same effective stress loading conditions (the difference between applied σ3and P is the same). The coloured dashed lines in Fig. 20a demonstrate that for the rocks under the same effective stress condition, the resultant peak strengths are similar. As a result, the numerical modelling results demonstrate that σyfor brittle rocks follows Terzaghi’s effective stress law. For the case of crack initiation strength, Terzaghi’s law is shown to be valid as well (Fig. 20b).

Fig.17. Plot of the confining stress-volumetric strain obtained from numerical simulations and laboratory testing of LdB granite under undrained hydrostatic loading.The experimental data for the undrained bulk modulus (Ku) of LdB granite are taken from Lau and Chandler (2004).

Fig. 18. Response of the models in hydrostatic compression tests with two different loading conditions for the measurement of (a) drained bulk modulus, and (b) solid matrix bulk modulus of the samples. The experimental data for the drained bulk modulus (K) and the bulk modulus of solid matrix (Ks) of LdB granite are taken from Lau and Chandler (2004).

As shown in Fig. 20c, the simulations are able to replicate the effects of pore pressure on the crack damage stress of the rock based on the notion of effective stress.Similar to Fig.20a and b,the differential stress for the crack damage of a particular rock is approximately the same at the same effective confining pressure.As mentioned earlier,the crack damage threshold can be identified by examining the plot of volumetric strain and axial strain,in which the reversal point in the plot refers to the crack damage stress.This reversal point marks the stress value in which the rock starts to experience dilation (transition from compaction behaviour to dilational deformation).The volumetric strain-axial strain curves in Fig.20d are used to identify the crack damage threshold for each of the loading conditions.Final fracture patterns obtained from biaxial testing of material under the same effective stress loading are illustrated in Fig. 21.

Fig.19. Concept of Terzaghi’s effective stress law: (a) Resultant effective stress acting on porous medium; and (b) Representation of impact of pore pressure on modifying the Mohr’s strength envelope.Note that effect of fluid pressure is to horizontally move the Mohr circle to the left by an amount equal to the magnitude of the fluid pressure.

6. Application of the model to simulating the formation of breakout and fluid-pressure response observed in the URL

The goal of this section is to assess if the short-term response of the undisturbed EIZ and EDZ around two case-history tunnels at a URL in Pinawa, Manitoba, Canada, can be predicted using the calibrated GBM.The URL is located approximately 120 km northeast of Winnipeg, Manitoba, Canada, within the massive LdB granite near the western edge of Canadian Shield.It consists of two main levels at depths of 240 m(240 level)and 420 m(420 level)below ground surface. The Mine-by Experiment. Tunnel and the Tunnel Sealing Experiment(TSX),which are excavated at a depth of 420 m,are two of a series of field studies conducted in order to characterise the hydro-mechanical properties of the EDZ.The in situ stresses at the tunnel sites are determined to be a sub-vertical minimum principal stress of σ3=11 MPa,a sub-horizontal maximum principal stress of σ1=60 MPa,inclined at 11°to the horizontal direction,and a subhorizontal intermediate principal stress of σ2= 43 MPa (Martin,1997). The Mine-by Experiment Tunnel with a diameter of 3.5 m was excavated in the LdB rock mass using a smooth non-damaging drill method.To maximise the EDZ,the tunnel axis was oriented to be parallel to the direction of maximum principal stress(σ1)while the minimum principal stress (σ3) was sub-vertical. The experiment involved installing instrumentation in an undisturbed volume of rock before excavating the tunnel. Extensometers are installed in the boreholes drilled from the adjacent underground galleries. The installed extensometers recorded rock mass deformation as a result of tunnel excavation.In addition,development of the EDZ as a result of the microcracking process and associated permeability changes were investigated at the tunnel site.The 40 m diameter TSX tunnel was oriented parallel to the direction of the maximum principal stress (σ1) in order to minimise the extent of the EDZ by reducing the concentration of tangential stresses induced adjacent to the tunnel.During the excavation process,the microseismic events, changes in microvelocities and rock transmissivity were monitored. Furthermore, pore pressure evolution for a period of one year was measured at various distances from the excavation (Martino and Chandler, 2004). The field measurements show that the general pore pressure response can be described as gradual rise over a time period of ～40 d, followed by gradual dissipation of pressure due to the pore fluid diffusion (desaturation) in the rock mass.

For the Mine-by Experiment,the excavation-induced maximum tangential stresses at the roof and floor were estimated to be approximately 160 MPa, while a stress magnitude of 100 MPa was determined at the same locations for the case of the TSX tunnel(Martino and Chandler, 2004). The dimensions of the developed EDZ around the TSX tunnel are illustrated in Fig. 22.

To simulate the formation of the EDZs around the Mine-by Experiment and the TSX tunnels, two DEM-Voronoi models with the average grain size of 3 cm were constructed.The validity of the DEM-Voronoi model to simulate the observed damage progression and its associated hydro-mechanical responses at field scale are investigated. The hydro-mechanical properties of the models are calibrated(as listed inTable 6)to the long-term properties of the LdB granite samples obtained from the URL.The studies by Potyondyand Cundall (2004) and Radakovic-Guzina et al. (2015) showed that in order to successfully simulate the notch development around tunnels, the model parameters must be calibrated to match the longterm strength of the rock (～0.4UCS). The short-term UCS of the LdB granite obtained from 420 level is approximately 150 MPa(Lau and Chandler,2004).The in situ(initial pre-excavation)pore pressure is approximatedto be around 3 MPa(Detournayand Berchenko,2001;Rutqvist et al.,2009).To avoid numerical instability,the materials within the tunnel boundary are excavated by gradually decreasing the Young’s modulus of grains to zero.

Fig.23a-e compares the simulated re-distribution of stress due to excavation around the Mine-by Experiment tunnel. The simulated stresses are obtained from the boundary element code Examine2D (RocScience Inc, 2015b) and the GBM in UDEC. To evaluate the elastic response of the UDEC model to the excavation of the tunnel,the values of the strength properties assigned to the grain contact are set to be 1000 times greater than the contact strength property values as listed in Table 6.As a result,the contacts will not yield under the applied loads.The maximum compressive stress that occurs in the roof and floor of the tunnel is approximately 160 MPa for elastic analyses(Fig.23a and b),while the same stress magnitude is estimated at the tip of developing notch(Fig. 23c). The simulated stress concentration corresponds very well to the measured value of 160 MPa at the tunnel site(Martino and Chandler,2004).For the case of the TSX model,the maximum compressive stress at the tip of the notch, which is approximately 95-100 MPa, corresponds to the reported stress value of 100 MPa measured in the field (Read, 2004). For both of the Mine-by Experiment and TSX models, the induced stresses create a zone of compressive stress in the roof,whereas materials adjacent to the sidewall of the tunnel experience tensile stress. The stress states and accompanied fracturing at the roof and sidewalls of the tunnel are responsible for the measured several orders of magnitude increase in permeability of those regions (Souley et al., 2001).Compared to the field observations and measurements(see Fig.22),the models are able to successfully simulate the observed dimension and depth of failure due to the spalling mechanism.The notch formation due to the progressive slabbing of material(Germanovich and Dyskin, 2000; Diederichs, 2007) is captured well. The developed notches at the roof and floor progressively develop until the induced compressive zone at the tip of the notch(process zone (Martin et al.,1999)) exceeds the tensile strength of the grain contacts. The presence of unloaded fractures in the EDZ creates a zone of increased permeability in the roof and floor of the models,while the permeability increase at the side of the tunnel is attributed to extensional microcracking of material and associated fracture opening.

Fig.20. (a)Peak strength,(b)crack initiation,and(c)crack damage thresholds for the DEM-Voronoi samples as a function of pore pressure,for several different values of effective stresses.(d)Volumetric strain-axial strain curves for biaxial compression simulations.Note that the reversal points in volumetric strain-axial strain curves are used to identify the crack damage thresholds for each biaxial simulation.

Fig. 21. Final fracture patterns obtained for the biaxial samples under the same effective stress loading condition. The white lines in the biaxial compression samples are the stress-induced fractures.

Fig. 24 shows the predicted induced pore pressure in the rock mass immediately after excavation operation is finished (shortterm response). The modelling results do not account for the effect of time on changing the pore pressure as a result of diffusion/recovery (desaturation/saturation) processes. It can be seen that the simulated pore pressures adequately match the field data, including the pore pressure spike at the tip of the developed notch. The pore pressure decrease in the tunnel sidewalls, due to the tensile cracking, was captured successfully.

As illustrated in Fig.25a-d,both the brittle and elastic models of the Mine-by Experiment are mostly able to predict the rock mass deformation to match the field deformation recorded by extensometers.The exception is the monitoring line located in the floor (Fig. 25c) for which the brittle model over-predicts the deformation. The recorded deformation along the monitoring lines at sidewalls and roof,as illustrated in Fig.25a,b and d,shows a good fit to the deformation recorded in situ. The discrepancy between the roof deformations obtained from the brittle model and the extensometer data (Fig. 25 d) is due to the existence of a detached block located at the roof of the brittle model which causes the recorded deformation to be larger than the in situ measurement.

Fig. 22. Development of damage around the Mine-by Experiment and Tunnel Sealing Experiment tunnels. The depths of failure observed in situ(Read,2004)are compared with the DEM-Voronoi simulation results. The induced fractures in the models are displayed with black dots.

Table 6Calibrated parameters for simulating damage progression around the URL tunnel.

7. Discussion

In this study,a GBM was first calibrated to the hydro-mechanical properties of tested laboratory samples of LdB granite and then was used to investigate the formation of damaged zones around deep tunnels.The model is capable of simulating the contribution of fully hydro-mechanical coupling processes between stress, flow and damage in the formation of fractured zones in low-porosity crystalline rocks. A cohesive crack constitutive model based on the theories of nonlinear fracture mechanics (NLEFM) (Dugdale,1960)is implemented in the distinct element numerical code UDEC to define the constitutive behaviour of the grain interfaces before and after yielding.

The new calibration procedure aims to find suitable values of micro-parameters of the model that can reproduce the macroproperties of the rock observed by physical laboratory tests. The calibration procedure presented herein is an iterative process that involves developing a series of Brazilian, UCS, and drained/undrained biaxial simulations to find a set of input parameters that yields correct hydro-mechanical response of the tested rock. The procedure involves calibration of firstly the macro-mechanical properties, and secondly the macro-hydraulic/poromechanical properties of the rock. This process begins with assigning the exact mineral deformability properties to the model grains based on reported data in the literature (e.g. Lama and Vutukuri, 1978;Mavko et al., 2009). Later, the adjustment of the contact stiffness values(kn0,cand ksh0,c)continues until the desired Young’s modulus and Poisson’s ratio for the material are achieved. The microstrength properties assigned to contacts need to be adjusted until the stress-strain curves for the simulated Brazilian and UCS tests resemble the laboratory-derived curves. By doing so, the model simultaneously yields the correct values for the following parameters: Ei, νI, Ti, UCS, and crack initiation and damage thresholds.Simulated fracture patterns and failure modes must also be inspected to ensure consistency with those observed on laboratory samples. After that, a series of biaxial test simulations under varying confining stresses needs to be modelled to adjust the contact friction angle,φc.The iterative refinement for values of cc,φcand Tcmust be carried out until the simulated biaxial strength envelope and the resultant macro-mechanical ciand φimatch those obtained from biaxial laboratory tests. Afterwards, the emergent strength and deformability of the model from UCS and Brazilian simulations must be checked again to confirm that they remain in agreement with the experimental values.

After the calibration of mechanical properties, a series of biaxial compression tests under a drained loading condition with different input initial and residual apertures(a0and ares)must be performed to adjust the emergent hydraulic properties of the model.The adjustment of the apertures should continue until the permeability-stress curve for the simulated biaxial test displays consistency with that obtained from laboratory measurements.In addition, the assigned values for the model components must be able to successfully reproduce the pore fluid pressure change as a result of deformation of the solid material.To assure that the pore pressure increase during the undrained compression of the model can be monitored, it is necessary to check if the resultant Skempton’s coefficient(B)confirms the response measured in the laboratory. Ultimately, a series of drained and unjacketed hydrostatic simulations, as explained in Section 4.3.2, is required to be carried out to derive the compressibility behaviour of the model(K and Ks). The values of K and Ksare used to calculate the Biot’s coefficient (α) of the model. By determination of the parameters controlling the poromechanical properties,the calibration process is completed.

Fig. 23. Comparison between the maximum principal stress (σ1) distributed around the Mine-by Experiment and TSX models. Stress re-distribution around Mine-by Experiment tunnel obtained from (a) elastic analysis using DEM-Voronoi model, (b) damage analysis using DEM-Voronoi model, and (c) elastic analysis using boundary element code(RocScience Inc, 2015b). Stress re-distribution around TSX tunnel obtained from (d) damage analysis using DEM-Voronoi model, and (e) elastic analysis using boundary element code. (f) Population of induced tensile and shear cracks for the Mine-by Experiment and TSX models.

Numerical experimentations demonstrate that the model has the ability to capture micro-mechanisms involved in modifying rock permeability during compressive loading. The biaxial samples accurately reproduce the permeability decrease in the elastic region of the stress-strain curve. For this range of deviatoric stress, permeability decreases from its initial value to approximately two orders of magnitude lower.The onset of microcracking(crack initiation) and the following crack growth do not significantly change the axial permeability. The reason for such hydraulic behaviour is attributed to the isolated nature of propagating cracks for this range of deviatoric stress(stress levels between crack initiation threshold and damage threshold).However, the accumulation of the microcracks that corresponds to an unstable crack propagation state results in a permeability increase of three orders of magnitude. Further progression of damage and accompanied dilation of the material lead to approximately six orders of magnitude of permeability increase for the post-failure state,compared to its non-damaged state.It is noteworthy to mention that the hydraulic properties of rock in the post-failure state are similar to the permeability of fractured rock within the highly damage zone(HDZ)around deep tunnels,where the fluid flow is governed by the fracture permeability of the broken material. These modelling results are in reasonably good agreement with the stress-dependent behaviour of the brittle rocks observed in laboratory experiments.Careful examination of the relationship between the occurrence of fracturing and the resultant damage geometries with permeability evolution in samples indicated that only the major family of the grown microcracks could produce the changes in permeability. This is consistent with the field measurements in which the maximum hydraulic conductivity is recorded in the zones where the borehole camera detected the majority of stress-induced damage along the borehole wall(Martino and Chandler,2004).In general,it can be said that the hydraulic behaviour of the laboratory-scale rock samples subjected to biaxial loading is analogous to the permeability of the sub-regions of the EDZ. In this sense, the hydraulic characteristics of EIZ, EDZouter, EDZinnerand HDZ are similar to the permeability of the rock in the elastic deformation,isolated stable cracking, strain-hardening unstable cracking, and post-failure states, respectively.

It has been demonstrated that the calibrated GBM is able to accurately produce the depth and shape of the brittle spalling failure that occurs in the roof and floor of the TSX and Mine-by Experiment tunnels in a Canada’s URL when the models are calibrated to the long-term strength of the granite.The re-distribution of in situ stresses after excavation results in the development of zones of compressive stresses at the top and floor of the tunnels,while tensile stresses occur in the sidewalls. In these zones, the heterogeneity of stress distribution due to the microstructure heterogeneity leads to the initiation and growth of the tensile and shear cracks. The microcracking mechanisms in forms of tensile and shear failures eventually result in creation of the V-shaped notches at the roof and floor of the tunnels (see Fig. 22). As expected,the observed permeability increases were achieved for the regions around the tunnels where the maximum degree of fracturing occurred(roof and floor).Moreover,in the case of the Mineby Experiment,the fracturing of rock initiated at the regions around the tunnels where the materials were subjected to tangential stress of 160 MPa. This value of tangential stress is consistent with the stress level at which crack initiation has been observed in situ. In general, good agreement was obtained between the deformation that occurred in the field and the modelling results.It has also been shown that the observed stress-induced pore pressure changes in the EIZ and EDZ could be explained and captured in the numerical modelling.

Finally,it should be noted that the ability of the model presented herein is restricted to predicting the fluid pressure response of a rock mass under quasi-static mechanical loading (relatively low loading rate). As a result, the calibrated model does not guarantee the reproduction of the fracture behaviour of the rock mass surrounding an excavation under dynamic loading conditions (e.g.earthquake).In particular,a sudden spike in the pressure of trapped pore fluid and the potential for hydrofracturing of material as a result of an explosion or earthquake may not be correctly simulated.

Fig.24. Excavation-induced pore pressure calculated from the Mine-by Experiment model.The field pore pressure responses due to the excavation are taken from Chandler et al.(2002) and Rutqvist et al. (2009). The contours show the magnitude of the pore pressure at the contacts.

8. Conclusions

The primary conclusions of this study are summarised as follows:

(1) The capability of the grain-based DEM (GDEM) models to reproduce the mechanical impacts of fluid pore pressure on the weakening of the fracture strength was examined. The numerical experimentations showed that Terzaghi’s effective stress law（σ′3= σ3-P）governs the peak strength(σy)of the porous rock under a drained loading condition. It has been shown that Terzaghi’s effective stress law also holds in predicting the crack initiation and damage thresholds of the low-permeable hard rocks.

(2) The discrete modelling approach was able to reproduce the pore pressure response of the rock under undrained loading conditions for both the laboratory and excavation-scale problems. The fluid pressure increase as a result of material deformation was correctly simulated quantitatively when the hydraulic parameters of the model were calibrated. A very good agreement was found between the predicted values of permeability during biaxial compression and those measured in the laboratory.Numerical results demonstrated that interaction of multiple cracks and creation of localised macro-fractures increased the permeability of the material significantly. Spontaneous formation of these fluid channels that occurs at crack damage stress and the resultant fluid flow increase are qualitatively and quantitatively similar to the permeability evolution mechanism observed in real lowporosity rocks. This ability is especially important when simulating the hydraulic conductivity increase in the EDZin-nerand creation of macro-corridors for the migration of particles dissolved in fluid.

(3) The GBM was able to accurately predict the magnitude of wall displacements occurring in the Mine-by Experiment at a Canada’s URL in Pinawa, Manitoba. The displacements were obtained from the installed extensometers at the tunnel site.Similar to the brittle spalling mechanism observed experimentally, the model successfully simulated the progressive formation of slabs driven by extension cracking. The calibrated model, which was validated against laboratory and field measurements,can serve as a reliable predictive tool for designing excavations in highly stressed grounds.

Fig.25. Comparison between the displacement profiles predicted by the model and recorded displacements in the Mine-by Experiment tunnel.The monitoring points at which the displacement data are collected are displayed in the plots.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors wish to thank the Nuclear Waste Management Organization (NWMO) of Canada, and the Natural Sciences and Engineering Research Council of Canada (NSERC) for funding this research work. The discussion with Drs. Jennifer J. Day, Giovanni Grasselli, Matthew Perras, Ioannis Vazaios, Chrysothemis Paraskevopoulou,Alireza Baghbanan,and Omid Mahabadi and Mr.Tom Lam helped with the preparation of this paper.