Effective Elastic and Hydraulic Properties of Fractured Rock Masses with High Contrast of Permeability : Numerical Calculation by an Embedded Fracture Continuum Approach

In this work, the hydromechanical modeling of the fractured rock masses was conducted based on a new numerical simulation method named as embedded fracture continuum (EFC) approach. As the principal advantage, this approach allows to simplify the meshing procedure by using the simple Cartesian meshes to model the fractures that can be explicitly introduced in the porous medium based on the notion of fracture cells. +ese last elements represent the grid cells intersected by at least one fracture in the medium. Each fracture cell in the EFC approach present a continuum porous medium whose hydromechanical properties are calculated from ones of the matrix and ones of the intersected fractures, thanks for using the well-known solution of the joint model. +e determination of the hydromechanical properties of the fracture cells as presented in this work allows to provide the theoretical base and to complete some simple approximations introduced in the literature. +rough different verification tests, the capability of the developed EFC approach to model the hydromechanical behavior of fractured rock was highlighted. An analysis of different parameters notably the influence of the fracture cell size on the precision of the proposed approach was also conducted. +is novel approach was then applied to investigate the effective permeability and elastic compliance tensor of a fractured rock masses taken from a real field, the Sellafield site. +e comparison of the results calculated from this approach with ones conducted in the literature based on the distinct element code (UDEC) presents a good agreement. However, unlike the previous studies using UDEC, which limits only in the case of fractured rock masses without dead-end fractures, our approach allows accounting for this kind of fractures in the medium.+e numerical simulations show that the dead-end fractures could have a considerable contribution on the effective compliance moduli, while their effect can be neglected to calculate the overall permeability of the of fractured rock masses.


Introduction
Evaluation of effective mechanical and hydraulic properties of a geological reservoir, represented as porous rock masses with complex fracture networks, is a major issue in many engineering applications such as oil and gas field recovery, CO 2 storage, and geothermal exploitation to just mention a few.
e porosity of the matrix constitutes the primary porosity and the essential storage capacity of the reservoir, while the capacity storage of fractures represents only some fraction of the total porosity.Likewise, the hydraulic conductivity of such fractured reservoirs is composed by the conductivity of the rock matrix (pores and their connectivity) and the conductivity of cracks/fractures; the later one is usually some orders of magnitude greater than the former.us, the presence of fractures induces in the reservoir a quite complex fluid flow and hence significantly affects the production rate of oil recovery or geothermal energy.e idea of dissociation of flow inside the fracture network from that of the matrix is first proposed in [1,2].Since then, the double permeability system is still used as a common nomination for fractured reservoirs.
Albeit a huge number of studies on fractured rock masses in the past, their modeling is yet a challenging and dynamic research topic, with two principal trends grouping continuum and discontinuum approaches [3][4][5][6][7][8][9][10]. In the former approach, the fractured medium is considered as an equivalent continuum medium whose properties are derived from an appropriate upscaling procedure by using a homogenization technique.For this approach, the effect of fractures is implicitly accounted for in the equivalent constitutive model that can be described by the principles of continuum mechanics: the properties of the equivalent continuum medium are homogenous and defined at any point of the domain.is approach has the advantage for solving problems of large scales with the high densities of fractures, but its results may be sensitive to the domain's size especially when the studied domain is smaller than the representative elementary volume [10].Otherwise, the homogenization-based continuum model may not take into account the connectivity effect of fractures, the effect of clustering and spatial distribution of fractures, or the individual characteristics of a fracture since these characteristics are, at best, taken into account through some statistical features of all crack sets.
For the discontinuous approach, the medium with an embedded discrete fracture network (DFN) is represented as an assemblage of blocks (discrete elements) bounded by a number of intersecting discontinuities [3-7, 10, 11].e overall behavior of the fractured medium is mainly affected by the interactions between the matrix blocks and fractures (represented by interfaces between blocks).Although this approach can better describe the behavior of the fractured medium on a small scale, it can be very expensive on computer memory and time simulation particularly for a large scale and/or longtime transient problems.e main issue for this approach is the reproduction by numerical models of the complexity of in field fractures.In fact, because of the complexity and variation of fracture properties at a given site, the statistics of geometric characteristics and the properties of fractures/cracks are issued from limited and potentially biased field measurements (e.g., one-dimensional (1D) borehole imaging or two-dimensional (2D) outcrop mapping).e challenge here is not only how to input a great number of fractures into the model but also the conceptual model to express the fracture network and its connectivity.
To overcome limitations, as well as to take advantage of these two-mentioned approaches, some hybrid approaches borrowing the concept of continuum models incorporating the effect of fractures explicitly were introduced under different names: "fracture continuum approach" [9,[12][13][14][15][16][17], "fracture-cell model" [18], or "embedded discrete fracture model" [19][20][21][22][23][24]. e common point and the principal idea of these approaches lies on the concept of "fracture cell" which represents the grid cell of the discretized medium intersected by one or many fractures in the medium.Each fracture cell presents an equivalent porous medium with its own properties calculated from the contributed properties of intact matrix and fractures intersecting the cell.us, the fracture cell properties are different from that of the intact matrix and vary from cell to cell in respect with the geometry and properties of intersecting fractures.is approach was successfully applied to the study flow problem in fractured rocks [9,[12][13][14][15][16]18] or the coupled flow and heat transfer process [17] in the fractured reservoir.Some applications were recently conducted for the coupled hydromechanical problem [21][22][23][24].
In this work, a hybrid method, called from now on as the embedded fracture continuum (EFC) approach, is chosen to study the coupled hydromechanical behavior of the fractured porous media.e proposed approach can be applied for all fractured media owing a weak to a high density of fracture.But for the demonstration purpose, in the present work, we focused notably in the case of high density with an application in a real fractured rock of the well-known Sellafield site.
e work is limited in the case of single phase flow.For some approaches on multiphase flow, the interested reader could consult [25,26].e structure of this paper is organized as follows.Following this introduction, the EFC approach based on the notion of the fracture cell with its equivalent (poroelastic and hydraulic) properties will be detailed.e implementation of this approach in the open source code deal.II [27] is validated through different numerical verification tests.As the first application, this developed approach will be utilized to investigate the effective (hydraulic and hydromechanical) properties of the rock mass taken from the Sellafield site [28,29].Finally, the paper will be finished with some conclusions.

Fracture Modeling by the Embedded Fracture Continuum (EFC) Approach
Shearing the same idea with some other approaches proposed in the literature such as the "fracture continuum approach" [9,[12][13][14][15][16][17], the "fracture-cell model" [18], or the "embedded discrete fracture model" [19][20][21][22], the embedded fracture continuum (EFC) approach is developed in this work to model explicitly the fracture network on the coupled hydromechanical behavior of fractured rocks.e common point of this approach with abovementioned models is the use of the fracture cell notion as the fundamental basic brick.Follow that, each fracture is modeled explicitly in the discretized porous medium by a group of grid cells (fracture cells) intersected by this fracture (Figure 1) [30].us, we distinguish a fracture cell which represents the grid cell of the discretized media intersected at least by one fracture [9,[12][13][14][15][16][17][18] with the grid cells that are not intersected by any fracture (called the matrix cell).If the matrix cells represent the continuum porous matrix owing the physical properties of the intact rock, the fracture cell properties represent those of an equivalent continuum porous medium counting for hydraulic and mechanical properties of intact rock and ones of fractures intersecting the cell.Note also that, in the EFC approach, we distinguish two types of fracture cell: type I represents the cell intersected by only one fracture, while type II gathers all fracture cells intersected by multiple fractures (more than one fracture), as indicated in Figure 2.
Based on the fracture-cell concept, the meshing procedure of fractured media in the EFC approach will be an easy task because the mesh is not necessary to be conformed to fracture geometry, and fracture could intersect whatever element and hence a simple uniform Cartesian mesh can be used.At this point, the present approach shares the same advantage as the well-known extended finite element method (XFEM) largely developed in the last decade [31][32][33].However, if this latter approach uses the lower dimensional interface elements to model the fractures, in our proposed approach, they are described by the equidimensional elements.

Advances in Civil Engineering
As highlighted in some previous works [9, 12-18, 21, 23, 24], it is required that the fracture cell size could be quite small in order to satisfy the computational error tolerance and to bring the geometrical description of fractures by a mesh to an acceptable level.is requirement could lead to a very heavy numerical calculation if the uniform Cartesian mesh (called hereafter as the global refinement to distinguish to the local refinement) is utilized in the whole studied domain.In order to satisfy the exigencies of the tolerated numerical error and minimize the number of degree of freedom and hence the computational time, the local mesh refining strategies may be one of the most appropriated choice.In this work, the local refinement technique based on the hanging node method [27,28,[34][35][36] is used.e efficiency of this technique in combination with the EFC approach will be demonstrated in the next section through some verification tests.

Equivalent Poroelastic Based on Micromechanical Solution.
As previously mentioned, each fracture cell in the model represents an equivalent continuum porous medium that their hydraulic and hydromechanical properties will be calculated from the properties of the intact rock and ones of the embedded fracture in the cell.In this subsection, the determination procedure of these equivalent properties will be presented.
In the literature, some scholars attempt to determine the equivalent properties of fracture cell by using the micromechanical solution.As for instance, in [10], the authors divided the fractured domain into different subregions, and they used Oda's crack tensor to calculate the upscaled properties for each element representing each subregion.By comparing the results with the alternative discrete fracture network (DFN) models, a reasonably good agreement was observed by these authors.However, the authors pointed out that some difficulties can be encountered in this approach based on Oda's crack tensor where the size of element becomes small or/and only one or a few fractures intersect the element.Recently, Figueiredo and his collaborators [23,24] propose to determine Young's modulus of the elements intersected by a fracture from the following equation: 1 where E m , k n , and h indicate, respectively, Young's modulus of the intact matrix, the normal stiffness of fracture, and the  element size.Note that the same formula was used in [37] to determine Young's modulus of the elements containing fracture that are defined as weak elements.
In effect, this last formula is extracted from the solution of the joint model as detailed below.Moreover, it is important to note that, in their works, Figueiredo and coauthors [23,24] assumed that the equivalent properties are isotropic and homogeneous for all fracture cells.However, the determination of other equivalent poroelastic properties like the Biot coefficient and Biot's modulus was not discussed in the abovementioned references.To validate the proposed solution, Figueiredo and coauthors [23,24] compared their numerical results with the analytical solution in a simple case of one vertical and completely open fracture.A relative difference smaller than 5% was reported by these authors in terms of the stress concentration result at the tip of fracture.
Despites that the use of a micromechanical approach to determine properties of the fracture cell could not be justifiable in all cases, various authors state its efficiency [10,23,24].e main concern for the use of such approach in the effective properties of the fracture cell is the size of the fracture cell often smaller than that of the representative elementary volume (REV).While the micromechanical approach based on Oda's crack tensor can be a good choice for the medium with high density of fractures embedded in a quite large fracture cell, the joint model seems to be more appropriate for small cell with only one or few fractures intersecting the element.is observation is coherent with the discussion of [38] who studied the poroelasticbehavior of jointed rock.Following these authors, the homogenization techniques based on the wellknown Eshelby's theory [39] can be used for short joints as a limiting situation of embedded ellipsoidal inclusions (known as the cracked medium) while for the long joints cross-cutting the REV, the joint model is suggested.
In this work, our goal is to model the rock masses weakened by a quite long fracture network, and thus the joint model could be appropriate to determine the equivalent properties of fracture cell.e detail about this type of model was largely discussed in the literature ( [38] and different references cited therein), and hence in what it follows, only some important results of this model will be captured.Furthermore, in our study, we assume that all the joints have the same Biot coefficient (which is equal to 1), while the Biot modulus of fracture is very large and hence, the overall poroelastic properties of the jointed medium can be written as follows [38]: where S hom , S m , and S i represent, respectively, the compliance tensor of the overall medium, intact matrix, and i th joint family.k i n and k i t designate the normal and shear stiffness, while n i , t i , and t i ′ indicate the unit normal and tangential vectors, and the parameter d i is the spacing of the i th joint set.
Consider now the simplest case of a fracture cell of type I that is intersected by one horizontal fracture (the normal vector of fracture coincides with the vertical direction n � e 3 ).For an isotropic intact matrix characterized by elastic properties (E m , ] m ), one can deduce straightforwardly from Equation (2) the equivalent transversely isotropic poroelastic properties of the fracture cell: where E fc 2 , E fc 3 , ] fc 12 , ] fc 31 , and ] fc 23 represent the five elastic modulus of the fracture cell; b fc 2 , b fc 3 , and M fc are, respectively, the Biot coefficients and Biot modulus of the fracture cell, while the parameter G m � E m /(2(1 + ] m )) indicates the shear modulus of intact isotropic rock.
For the general case of a fracture cell of type I intersected by an arbitrary orientated fracture (Figure 3(a)), the equivalent properties calculated in the global coordinates system (e 1 , e 2 , and e 3 ) can be determined from ones calculated in the local coordinate system associated with the fracture (n, t, and t ′ ) through a rotation with an angle θ(θ ∈ [0, π]) as shown in Figure 3 around the axis e 1 .In the local reference (n, t, and t ′ ), the equivalent poroelastic properties of the fracture cell are merely those previously obtained in Equation (3) with n coinciding with e 3 and d in that case is calculated as follows: For that case, the general form of the equivalent poroelastic properties (compliance tensor and Biot tensor) of the fracture cell in the global reference (e 1 , e 2 , and e 3 ) after rotation becomes as follows: where For the fracture cell intersected by multiple inclined fractures (fracture cell type II), its equivalent poroelastic properties by using the solution of the joint model (Equation ( 2)) have the same matrix form as presented in Equation (5) and their components depend also on the elastic modulus (E m and ] m ) of the intact matrix, the mechanical properties (k i n and k i t ), and the corresponding inclination angle θ i of each fracture i that intersects the cell.e expression of the diagonal components S 11 , S 22 , and S 33 written in this case are as follows: e other nonzero components of the compliance tensor and Biot tensor are not presented here just to simplify the presentation.

Isotropic Approximation of the Equivalent Poroelastic
Properties of Fracture Cell.As previously mentioned, in their work, Figueiredo and coauthors [23,24] considered that the elastic properties of all fracture cells generated in the whole medium are isotropic and are characterized by the same Young's modulus of the elements calculated from the fracture cell intersected by a horizontal fracture.Otherwise, Poisson's ratio is equal to one of the intact rock: where E fc 3 is defined in Equation (3).More precisely, in these last works, one does not distinguish fracture cells intersected by one or many fractures.It means also that this approximation is only applicable for fractures owing the same mechanical properties.
In this work, another isotropic approximation of the equivalent poroelastic properties of the fracture cell is proposed.Following this, the isotropic compliance tensor of this proposed approximation has the in-plane compliances coefficients equal to the weakest ones calculated from the joint model.us, in difference with the previous approximation (Equation ( 8)), the isotropic approximation introduced here allows to model the fractures having different mechanical properties by accounting for the contribution of each fracture in the isotropic properties of the fracture cell.More precisely, for the fracture cell of type I intersected by an inclined fracture, their approximated isotropic moduli are proposed as follows: Advances in Civil Engineering Similarly, for the fracture cells of type II, intersected by more than one fracture, the effective properties are obtained using the hypothesis of noninteraction leading to the summation of all fracture contributions: To distinguish with the isotropic approximation proposed in this work, we call hereafter the approximation proposed in Equation ( 8) as Figueiredo isotropic approximation.
e corresponding isotropic Biot coefficient and Biot modulus of the fracture cell can be deduced as function of the isotropic elastic modulus of the intact matrix and ones of fracture cell (Equation ( 2)) as follows: anks for its simplicity, the isotropic approximation of the equivalent poroelastic properties of the fracture cell will be used in our proposed EFC approach.e accuracy of the two approximations as detailed here will be investigated in the next section of the paper.

Equivalent Permeability of Fracture Cell.
e determination of the equivalent permeability of the fractured rock masses has been discussed in a lot of work since several decades [9,15,16,18,40,41].It is largely accepted that, in fractured rocks where the permeability of the intact matrix is negligible as compared to the permeability of fracture, the flow is principally conducted through the fracture inside the domain.e permeability of the fracture K f is calculated from the fracture's hydraulic aperture that, in the 2D case, assuming smooth fractures coincides with the width w of the fracture.Combined with laminar flow hypothesis, the permeability K f of the fracture is finally obtained as follows [42,43]: e effective permeability of the fracture cell of type I is obtained by considering the general case of an arbitrary fracture orientated to an angle θ to the direction of pressure gradient as shown in Figure 3(b).In Figure 3(b), two constant heads p 1 and p 2 are applied, respectively, on the left and on the right of a fracture cell, while the top and bottom sides of the cell are considered to be impermeable.Considering that the flow in the fracture cell is exclusively through the fracture, permeability of the matrix is small with respect to the permeability of fracture, the equivalent permeability of the fracture cell in the x direction [40,41]: 6

Advances in Civil Engineering
Similar result is obtained for the permeability of the fracture cell in the vertical direction as follows: ese equations show that the effective permeability of the fracture cell equals that of the fracture intersecting the cell, corrected for its orientation and multiplied by the ratio between the cell size and hydraulic aperture.Clearly, the orientation of the anisotropic permeability of the fracture cell depends on the orientation of the fracture with the highest value in the direction of the plane of the fracture.e permeability tensor in the global coordinate could be obtained again by calculating firstly the permeability in the local fracture coordinate system and then making a rotation to bring results on the global coordinate system.
Similar to the mechanical case, for simplicity, we propose to use for a fracture cell (type I) an isotropic hydraulic permeability that is defined as the maximum between K x and K y : Correspondingly, for the fracture cell of type II, the effective permeability is obtained by adding the contribution of all cracks intersecting the fracture cell: e isotropic approximation of the permeability in this last case is obtained by generalizing Equation ( 15):

Verification Tests
In this section, through some numerical verification tests, the accuracy of the developed EFC approach will be verified.
In the first stage, two examples which represent correspondingly the purely mechanical and purely hydraulic problems are considered with the focus on the convergence solution of EFC with respect to the fracture-cell size and on the correctness of the isotropic approximations of the equivalents (hydraulic and hydromechanical) properties of the fracture cell.For the first purpose, the size of fracture cell will be gradually decreased by using both the global refinement with uniform meshes and local refinement technique utilizing the hanging nodes.e efficiency of this local refinement technique in combination with the EFC approach will be highlighted.Finally, the same methodology will also be applied in the third example to verify the validity of the proposed EFC approach in the coupled hydromechanical modeling of fractured rocks.

Purely Mechanical Test.
As the first test case, we consider an elastic medium with dimension L � 1.0 m of width and H � 1.25 m in height which is intersected by a fracture with 1 mm of width as shown in Figure 4 Different numerical simulations can be conducted to study this problem.In the first step, the commercial code Flac3D (Itasca Consulting Group Inc.) is used in which the fracture is explicitly modeled as the elastic interface elements (Figure 5(a)).e results obtained from this simulation will serve as the referent results for the comparison purpose of the results calculated from the EFC approach.In addition, for this last approach, two meshing configurations of the fractured medium are considered.
e first configuration uses a uniform mesh, resulting from the global refinement technique as shown in Figure 5(b).For this configuration, the commercial code Flac3D is also chosen in which we implemented the EFC approach.is simulation demonstrates the easy use of the EFC approach presented in this work in an existence commercial code when only a simple algorithmic to search the mesh cells intersected by the fracture is required.e other configuration corresponding to the case uses the nonuniform meshing based on the local refinement technique (Figure 5(d)) to describe the fracture in the EFC approach.
e numerical simulation of this second configuration is performed using open source library deal.II.e reasons for the choice of this code are not only on its capacity to deal with different problems, as demonstrated in numerous works ( [44][45][46][47] and different references cited therein) but also on its capacity for local refinement of meshes using the hanging nodes technique [27,30,[34][35][36].Finally, as shown in Figure 5(c), the fracture can also be modeled with its real geometry by using the conformed mesh model which is also simulated in deal.II.e grid cells in this case coincide with the boundaries of the fracture, i.e., the size of elements is exactly equal to the fracture aperture (h � w � 1 mm) and mechanical properties of elements are equal to the elastic moduli of the fractures (E f , ] f , and G f ).ese latter parameters can be determined from the Advances in Civil Engineering normal and shear sti ness of fracture (k n and k t ) through the following relationships [48]: In Figure 6 are depicted the results of the vertical displacement following a vertical pro le at the middle of the domain.ese results are calculated with di erent models of fracture inclined at 45 °to the horizontal direction.As the rst statement, a superposition of the results determined from the models based on the interface elements and the conformed mesh can be noted as expected which demonstrated the correctness of the numerical simulation by deal.II and justi es Equation (18) deriving the relationships of the mechanical properties of fracture.Our attention is then focused to the results evaluated from the EFC approach, which were conducted on the two other models using the uniform and nonuniform meshes and the isotropic approximations of the fracture cell properties.For the clarity purpose, we detail here that, in the uniform mesh model 128 × 156, elements are used which is similar to the interface element model.e fracture cell size of the uniform mesh model (which is about h 7.8125 mm) is similar as one chosen in the simulation conducted on the nonuniform mesh model.However, in this last model, thanks to the local re nement technique, the matrix cell size can be taken much more important (about 8 times greater) with respect to the one of the uniform mesh model.Consequently, both models use the same size of the fracture cell, but only 920 elements are necessary in the nonuniform mesh model as compared to 19968 elements for the uniform mesh model.As illustrated in Figure 6, the results achieved from these two models are similar which demonstrated the e ciency of the local re nement based on the hanging node technique in combination with the EFC approach to model the fractured medium.
Note that in all these simulations the isotropic approximations of e ective properties of the fracture cell are used.e results exhibit a quite signi cant di erence between the two isotropic approximations of the equivalent mechanical properties of the fracture cell.In comparison with the referent results provided from the interface element model, a maximum di erence of about 2.86% is noted for the isotropic approximation proposed in this work, while the di erence of up to 12.5% is evaluated for the approximation presented by Figuereido et al. [23,24].e gap between the two isotropic approximations can be explained by the fact that, in the Figuereido isotropic approximation, Young's modulus is assumed constant and independent on the inclination angle of fracture.In addition, as shown in Figure 7, the oriented fracture (represented by an inclination angle θ) can present an important e ect on the equivalent elastic properties of the fracture cell, particularly in case of large fracture cell size (with respect to the fracture aperture) and high inclination angle (around 45 °).For a chosen fracture cell size h, the elastic Young modulus of the fracture cell is smallest in case of θ 45 °and highest in case of θ 0 °(or θ 90 °).us, the isotropic elastic Young's modulus proposed by Figuereido et al. [23,24] is overestimated in case of inclined fracture.As consequence, the approximation of these last authors seems appropriated only in cases of quasi-horizontal or quasi-vertical fractures.
As indicated by di erent scholars [9, 12-18, 21, 22], the correctness of the approaches based on the fracture cell concept could depend on the size of this element.is discussion can be veri ed now when we investigate the in uence of this parameter on the convergence of the EFC approach.Note that from the results presented in the previous paragraph, in the last part of this paper, we discuss only on the EFC approach using our isotropic approximation and the local re nement technique just to simplify   Advances in Civil Engineering presentation.In deal.II this could be done easily by changing the local refinement level around the fracture while keeping a constant prechosen matrix cell size.In addition, our investigation will be carried on with different inclination angles of fracture ranging from 0 °to 45 °.
In Figure 8(a), we present the vertical displacement determined at the controlled point "A" located on the middle of the top boundary.ese results are evaluated from different values of fracture cell size and inclination angle of fracture ranging from 0 °to 45 °.In these simulations, the matrix cell size is fixed at 62.5 mm, while the fracture cell size is gradually decreased (from h � 62.5 mm to h � 1.9531 mm) by changing the local refinement level around the fracture in deal.II.In general, we can state that the vertical displacement at this controlled point decreases by decreasing the size of the fracture cell. is displacement tends to an asymptotic value which approaches to the results calculated with the conformed mesh (case h � w).However, from a fracture cell of size h � 7.81 mm, it can be noted that the relative difference of results is smaller than 1% for the whole considered range of inclination angle of fracture.
Figure 8(b) highlights the displacement of the controlled point calculated from the EFC approach in comparison with the referent results based on the elastic interface element model (which is exactly the results conducted with the conformed mesh model).It shows that the EFC approach could underestimate the result for the inclination angle of fracture ranging from 0 °to about 30 °, while it could overestimate for the case of fracture owning the inclination angle from 30 °to 45 °.However, for all considered cases of inclined fracture, the relative difference is smaller than 3%. is difference is small enough to validate the EFC approach using the isotropic approximation proposed in this work in modeling the mechanical behavior of the fractured medium.Advances in Civil Engineering

Purely Hydraulic Test.
In this second test, the same geometric model as the previous example (Figure 4(b)) is considered to simulate the hydraulic behavior in the fractured medium.e validation of the EFC approach is also done by comparing with the referent results provided from the calculation conducted on the conformed mesh model.
More precisely, for the referent case, the fracture is modeled explicitly with his real aperture (w �1 mm), and the permeability determined from Equation ( 12) is equal to ). e permeability of the porous intact matrix is taken equal to 2.4 × 10 −15 (m 2 ), a typical value for a Sellafield granitic rock (UK) as reported in [49].us, the   Advances in Civil Engineering permeability of fracture is seven times higher than the matrix permeability.e viscosity and compressibility of fluid are 0.001(Pa.s)and 5.0 × 10 −10 (Pa −1 ), respectively.Note here that the permeability of the fracture cell is supposed to be isotropic and calculated from the permeability of fracture following Equation (15).As shown in Figure 4(b), the top and bottom boundaries of the model are considered to be impermeable, while a constant gradient pressure is imposed on two lateral boundaries by applying the pressure p 1 � 10 (MPa) on the left and p 2 � 0 (MPa) on the right of the model.Similar to the previous mechanical test, in the first stage, we will investigate the convergence of the EFC approach by modeling the fracture with different local refinement levels corresponding to different values of fracture-cell size ranging from h � 62.5 mm to h � 1.95 mm. e results depicted in Figure 9(a) present the horizontal fluid flux across the fractured medium, which tends to decrease with respect to the fracture cell size.However, the obtained results seem to converge at a constant value from a local refinement of mesh corresponding to the fracture cell size about h � 15.625 mm with relative differences of about 1% for all considered inclination angles of fracture.
In the second stage, the results calculated with the EFC approach (with the fracture cell size h � 7.8125 mm) are  Advances in Civil Engineering compared with the referent results conducted with conformed mesh configuration as illustrated in Figure 9(b).A quite similar tendency as in the purely mechanical problem can be observed.An overestimation of the fluid flux for the inclination angle between 0 °and 40 °is noted, while for the inclined fracture of about 45 °, it seems that the EFC approach proposed here may slightly underestimate the result.In general, a maximum difference of about 5% (case of fracture with inclination angle 25 °) between the two models is observed which seems acceptable to validate the capability of the EFC approach to solve the hydraulic problem of fractured media.

Coupled Hydromechanical (HM) Test.
As the last verification test, in this part, we use the EFC approach to model the coupled hydromechanical behavior of the fractured medium.To this end, we retain the outer geometry (L � 1.0 m and H � 1.25 m) and the mechanical properties and hydraulic properties of the porous matrix as the two previous examples.
e two other necessary coupling parameters of the porous matrix, namely, the Biot coefficient and Biot modulus are taken, respectively, equal to b � 0.8 and M � 20.8 (GPa).Two configurations are considered: the first one consists of a short fracture (l � 0.6 m) placed in the center of the domain and orients at 45 °(Figure 10(a)) in respect with the horizontal direction, and the second case accounts for two fractures (with the same length l � 0.6 m) with two (arbitrary) inclination angles of 25 °and 90 °(vertical fracture), respectively, as illustrated in Figure 10(b).All the parameters of the fractures like the aperture, mechanical, and hydraulic properties are taken as in the previous examples, i.e., w � 1 mm, k n � 434 (GPa/m), k t � 86.8 (GPa/m), and K f � 8.33 × 10 −8 (m 2 ).e boundary conditions for this test are also highlighted in Figure 11 where the left, right, and bottom boundaries of the saturated medium are fixed and impermeable, while on the top boundary, it is permeable, and a compressive normal stress σ 1 � 24.46 (MPa) is applied.Again, the validation of the EFC approach will be conducted by comparing with the referent results calculated from the conformed mesh model (Figure 11) in which the real geometry and poroelastic properties of fractures are used.For this purpose, two controlled points "A" and "B" located, respectively, at (0.5 m and 1.25 m) and (1.0 m and 0.0 m) are selected to observe the results during the transient calculations (Figure 10).By using the same size, the hydraulic and mechanical properties of the fracture cells (case of the nonuniform mesh model) are similar as ones of the previous tests while their coupling HM parameters are determined from Equation (11).
e comparison of the results determined from the two models is presented in Figure 12, where the vertical displacement at the controlled point A and excess pore pressure at the controlled point B are traced against the elapsed time.As expected, during the consolidation, the decrease of the excess pore pressure induces an increase of effective stress which explains the increase of the vertical displacement in the medium versus elapsed time.A very good agreement of the results was calculated from two models (nonuniform mesh and conformed mesh) for both configurations of embedded fractures in the porous medium.
e relative error seems more significant in the case of two fractures which is about 4.5% for the displacement and 1.5% for the result of pore pressure.Again, these moderate differences allow confirming the validity of our proposed EFC approach to model the HM behavior of the fractured media.

Application of the EFC Approach to Investigate the Effective Permeability and Elastic Compliance Tensor of Fractured Rock Masses in Real Field
e capacity of the EFC approach will be also demonstrated in this section through a numerical investigation of the effective permeability and elastic compliance tensor of the real rock mass taken from the Sellafield site.e abundancy of information about this site, thanks to the investigation program performed by United Kingdom Nirex Limited [28,29] as part of the second Bench Mark Test (BMT2) of the international collaborative program DECOVALEX III, is the main reason of this choice.Additionally, the information of the site is also well-synthetized in various works [4][5][6][7][49][50][51]. e necessary information about the spatial distribution of fractures is summarized in Table 1 and is sufficient to generate without difficulty the fracture network in the rock mass.ese parameters consist of input data about the fractures network such as the distribution of fracture's length, orientation, location, and fractures' aperture.Following this, the fracture trace lengths can be characterized by the power-law distribution [4][5][6][7]49]: where N F is the number of fractures (with fracture length greater than the length L) per unit area, C is the constant density, and D is the fractal dimension.
Concerning the orientations of fractures, the Fisher distribution [4-7, 28, 29, 50] is largely adopted.It was also pointed out that, in this site, there are four principal sets of fractures (Table 1) and the probability of a fracture of angle θ deviating from the mean orientation angle of a fracture set is given by Baghbanan and Jing [51]: where Fisher's constant K is assigned for each fracture set [28,29].Besides, in our investigation, the fractures with uniform aperture w � 65 μm as indicated in the works of Min and his collaborators [4,5] will be chosen, while the spatial location of fractures' center is represented by the Poisson distribution as largely accepted in the literature [4][5][6][7]51].
In the upscaling procedure to determine the effective properties of the fractured rock masses, the investigation of REV is an important step.e REV size is defined as the minimum size beyond which the upscaled values are stable [52].For the fractured rock masses in the Sellafield site, Min and his collaborators [4,5] by using the UDEC code demonstrated the existence of the REV whose size can be chosen from 2 m to 6 m with the coefficient of variation 12 Advances in Civil Engineering varying from 10% to 5%, respectively [4], for the purely mechanical problem.A quite similar REV size was also proposed for the hydraulic problem [5], which can range from 5 m to 8 m with the corresponding 20% and 10% acceptable variations, respectively.Based on these discussions, in this work, the REV size of 5 m is chosen to calculate the effective permeability and mechanical properties of this fractured rock masses using the EFC approach.

Methodology for the Determination of the Effective Elastic
Properties. e effective elastic properties of fractured rock masses in the Sellafield site will be calculated in this work based on the 2D plane strain hypothesis.e same methodology to determine the effective compliance matrix in Equation (21) as proposed in [4,53] will be adopted here.For plane strain hypothesis, the linear elastic behavior of a general anisotropic media written in the matrix form is as follows [4,53]:  Advances in Civil Engineering 13 Note here that the coordinate Oxyz presents the global coordinate (z � 1, x � 2, y � 3) as shown in Figure 3(a).
Assuming that the Young modulus and Poisson ratios in the z direction are equal to the ones of the intact matrix (] zx � ] zy � ] m , E z � E m ), it follows S 12 � S 21 � S 13 � S 31 � −] m /E m , S 11 � 1/E m .e study in [4] also revealed that the values of S 14 , S 24 , S 41 , and S 42 are very small in comparison with the other components, while S 34 and S 43 are equal to zero due to the fact that shear stress σ xy does not affect the z direction deformation.us, the compliance tensor in Eq. ( 21) can be rewritten in an explicit form as follows:

Advances in Civil Engineering
Hence, calculating the effective elastic properties of the fractured rock masses reduces to the determination of only five components S 22 , S 23 , S 32 , S 33 , and S 44 which can be conducted through three numerical experiments (independent loads), as illustrated in Figure 13 (also see the details in [52,53]).
e first loading set (Figure 13(a)) represents a uniaxial compressive test with the normal stress f a x applied in the x direction.Two components S 22 and S 32 can be calculated from this test taking into account the condition of ε (a)  zz � 0 (hence, σ a zz � ] m f a x ): x f (a) x . ( Similarly, the two other components S 23 and S 33 can be calculated from the second loading set (Figure 13(b)): e last component S 44 can be determined from the third loading set (Figure 13(c)): Note that ε (a) xx , ε (a) yy , and ε (b) xx , ε (b) yy designate the averaged normal and lateral strain under uniaxial stress in x and y directions correspondingly, while ε (c)  xy indicates the average shear strain under pure stress loading.
From these components of the effective compliance tensor, we can deduce without difficulty the in-plane effective Young modulus, Poisson ratio, and shear modulus of the equivalent medium: (26)

Methodology for the Determination of the Effective
Permeability.e evaluation of the effective permeability of the fractured rock masses in each direction x or y is calculated with a constant hydraulic pressure gradient imposed on the boundaries of the corresponding direction (x or y) and no-flow boundaries of the other direction (Figure 14).With the assumption that the permeability tensor is diagonal in the global coordinates system and by adopting Darcy's law for the upscaled fractured media [5], their effective permeability K x and K y following the corresponding the x and y directions are determined from the following equation: where Q x and Q y are the steady flow rates at the boundaries of the corresponding x and y directions, A is the crosssectional area, and μ is the dynamic viscosity of fluid.
In the numerical simulations, all terms on the right-hand side of Equation ( 25), except for K x and K y , are specified, while the steady flow rate Q x and Q y are calculated numerically.us, K x and K y can be back calculated.

Monte Carlo Simulation of the Effective Hydraulic and Mechanical Properties of Sellafield Fractured Rock Masses.
To obtain the representative results of the effective hydraulic and mechanical properties, the Monte Carlo simulations were carried out to obtain the effective mechanical and hydraulic properties.For that, ten realizations of the REV of Sellafield fractured rock mass were generated using the synthetized parameters in Table 1.As in the previous verification tests, an investigation of the fracture cell size on the results of the effective properties was also carried out.However, to not complicate the presentation, this study will not be detailed here.e results of this investigation indicate that a fracture cell size of about h � 9.77 mm can be used.
is fracture cell size presents an acceptable convergence of the effective properties of fractures rock masses when a maximum relative difference of about 6% can be stated in comparison with the case of the fracture cell having two times smaller size (h � 4.89 mm).
It is worth to note that, in the literature, the effective permeability and elastic compliance tensor of fractured rock masses are usually numerically determined by utilizing the universal distinct element code (UDEC) for only sample in which the dead ends and isolated fractures were eliminated [4][5][6][7][54][55][56][57].In this work, the EFC approach as proposed in work allows to simulate both two types of samples with and without dead-end fractures.e ten realizations of the REV of Sellafield fractured rock masses accounting for the dead-end and without dead-end fractures are depicted in Figure 15.
In Table 2, the mean value, standard deviation, and coefficient of variation of the effective permeability and elastic properties of the Sellafield fractured rock masses calculated with these ten realizations of REV are summarized.
e obtained results show that this fractured rock mass can be considered as isotropic.More precisely, the mean values of Young's moduli and permeability following the two directions x and y calculated with the dead-end fractures model are correspondingly equal to 16 Advances in Civil Engineering  e corresponding results of this rock calculated in the model without dead-end fractures are E � 25.24(GPa), ] � 0.225, and K � 1.23 ×10 −13 (m 2 ).ese last results seem consistent with the ones presented in [4,5,51].More precisely, by using the distinct element method (UDEC), these authors showed that the normalized elastic moduli in x and y directions are about 31% and 34% than that in intact rock (E m � 84.6 GPa) meaning that the averaged Young's modulus of this fractured rock masses is about E � 27.5(GPa).Concerning Poisson's ratio, a mean value ] � 0.35 presented by these authors seems higher than our result.With a 95% confidence level, the permeability in x and y directions determined by [5] are, respectively, K x � (1.00 ± 0.099) × 10 −13 (m 2 ) and K y � (1.29 ± 0.09) × 10 −13 (m 2 ) which present an averaged value of K � 1.15 × 10 −13 (m 2 ).18

Advances in Civil Engineering
A comparison of the mean results calculated from the dead-end and without dead-end fractures models highlighted a quite significant difference of about 17% of the effective Young and shear moduli, while for the effective permeability, a more moderate difference of about 4% is observed.is difference confirmed the discussion of Yang et al. [53] who concluded that the results calculated from the model without dead-end fractures overestimate the deformation modulus of the fractured rock masses.

Conclusions
We introduce in this work a new numerical simulation method to model the coupled hydromechanical behavior of fractured porous media.is method called the embedded fracture continuum (EFC) approach allows to model explicitly the embedded fractures network in the continuum porous media by using the fracture-cell concept.Each fracture cell in this approach presents an equivalent continuum porous medium that has its own properties calculated from the contributed properties of intact matrix and fractures.
e details of the determination procedure of equivalent poroelastic properties of the fracture cell are derived in this paper using the well-known joint model.rough different verification tests, the capability of the EFC approach is demonstrated, while the isotropic approximation of both the hydraulic and poroelastic properties of the fracture cell seems an appropriate choice.is EFC approach which can be easily implemented in a commercial or open-source code introduces an efficient numerical tool to model the hydromechanical behavior of the fractured media.To overcome the required small size of fracture cell in the EFC approach, the local refinement technique based on the hanging nodes is considered and presents its efficiency.As the first application, in this paper, we use this approach to investigate the effective permeability and elastic properties of the Sellafield fractured rock mass.e Monte Carlo simulation conducted with different generations of the fractures network in this rock presents a good consistency with the available results in the literature calculated from the distinct element method.

Figure 1 :
Figure 1: Sketch of fracture cell and matrix cell concept in the EFC approach using the uniform Cartesian mesh obtained from the global refinement technique (a) and nonuniform Cartesian mesh based on the local refinement technique with the hanging nodes (b).

Figure 2 :
Figure 2: Two types of fracture cells in the EFC approach: type I represents the grid cell intersected by only one fracture (a) and type II indicates the grid cell intersected by more than one fracture (b).

Figure 3 :
Figure 3: Inclined fracture cell with its local coordinate system is characterized by an inclination angle θ(θ ∈ [0, π]) with respect to the global reference Oxyz (case 2D plane strain model with t′ � e 1 ) (a); the fracture cell considered in the hydraulic problem (b).

Figure 4 :
Figure 4: Geometric model with boundary conditions of the purely (a) mechanical and (b) hydraulic veri cation tests.

Figure 5 :
Figure 5: Modeling of the oriented fracture with inclination angle θ � 45 °in the medium using the interface element model (a), uniform mesh model with the global refinement technique (b), the conformed mesh model (c), and the nonuniform mesh obtained from the local refinement technique with hanging nodes (d).

Figure 7 :Figure 6 :
Figure 7: Young's modulus of the fracture cell E iso calculated from the two isotropic approximations as function of the fracture cell size h (in mm) and inclination angle (in degrees) of fracture: (a) E iso versus h; (b) E iso versus inclination angle of fracture.

FIGURE 8 :
FIGURE 8: Vertical displacement at the controlled point "A" with respect to the size of the fracture cell (a) and in comparison with the referent results (using the interface element model or conformed mesh model) and the results calculated from the EFC using the nonuniform model based on the local refinement technique (b).

Figure 9 :
Figure 9: Horizontal flux going through the fractured medium and its convergence as function of the fracture cell size and inclination angle of fracture (a); comparison of the horizontal flux determined from the conformed mesh and nonuniform mesh models (b).

Figure 11 :Figure 10 :FIGURE 12 :
Figure 11: Hydromechanical simulation of the fractured medium with two fractures: case of the nonuniform mesh model using the local refinement technique (h � 7.81 mm) (a); case of the conformed mesh model (b).

Figure 14 :
Figure 14: Boundary conditions for the effective permeability calculation: (a) x direction; (b) y direction.

Figure 13 :
Figure 13: ree loading sets: uniaxial compression in the x direction (case a) (a), uniaxial compression in the y direction (case b) (b), and pure shear test (case c) (c) to determine five homogenized components S 22 , S 23 , S 32 , S 33 , and S 44 of fractured rock masses.

Figure 15 :
Figure 15: Ten realizations of the REV with 5 m size of the Sellafield fractured rock by considering or not the dead-end fractures.
S 11 S 12 S 13 S 14 0 0 S 12 S 22 S 23 S 24 0 0 S 13 S 23 S 33 S 34 0 0 S 14 S 24 S 34 S 44 0 0 0 0 0 0 S 55 S 56 0 0 0 0 S 56 S 66 [28,29]fracture is inclined with respect to the horizontal axis by an angle ϕ.Due to the symmetry of the considered problem, we consider only the inclination angle ϕ ranging from 0 °to 45 °(with interval 5 °for each case of study).eelasticproperties of the isotropic elastic medium consist of elements of Young's modulus E m � 84.6 (GPa) and Poisson ratio ] m � 0.24, while the normal and shear stiffness of the fracture are, respectively, k n � 434 (GPa/m) and k t � 86.8 (GPa/m).esearetypical values published for a granitic fractured rock at the Sellafield site in UK[28,29].On the top boundary of the considered model, a normal compressive stress 24.46 MPa is applied, while on the bottom boundary, all the displacements are fixed.

Table 1 :
Selected parameters for the generation of the fractures network in the Sellafield fractured rock mass.(x max − x min ) and y i � y min + R y,i (y max − y min ) where (x min , x max ,) and (y min , y max ,) are coordinate ranges.R x,i and R y,i are number in the range [0, 1] i and y i ) of every fracturex i � x min + R x,i

Table 2 :
Effective permeability and elastic properties of Sellafield fractured rock masses calculated with two models with and without deadend fractures.