Israel Cañamón

Rachid Ababou

Adrien Poutrel

Ángel Udías

Ouverte Upscaling 3D coupled hydro mechanical properties of fractured porous rocks

Keywords:

Biot coefficient and modulus Upscaled equivalent continuum

A methodology for upscaling 3D coupled Hydro-Mechanical (HM) properties of fractured porous rocks is developed theoretically, tested on synthetic fractured rock samples, and applied to a real site. The upscaled HM equations take into account the HM coupling in the dual matrix/fracture medium, comprising the cracks system as well as the intact porous matrix, yielding the equivalent stiffness tensor and two different sets of equivalent tensorial Biot coefficients B ij (I) and B ij (II) and moduli M (I) and M (II) for the upscaled system (I and II become identical only under certain hypotheses). We provide detailed theoretical expressions in the general tensorial case,and in the particular case of statistically homogeneous and isotropic cracks. The real site application is performed in a damaged and fractured claystone around the GMR gallery (Meuse/Haute-Marne Underground Research Laboratory). The geometric structure of the fracture set around the "Excavation Damage Zone" of the GMR gallery is described by a hybrid model comprising: (i) a set of 10 000 fissures with radially inhomogeneous statistics (size, thickness and density increasing towards the wall) and; (ii) a deterministic set of large curved fractures, periodically spaced along the axis of the gallery according to a 3D chevron pattern. Both "3D" and "2D transverse" distributions of the upscaled coefficients are calculated, and displayed using spheres or ellipsoids. Global tensorial coefficients are also obtained by upscaling the entire annular fractured zone. Equivalent isotropic coefficients are extracted from these tensors: Young modulus E, bulk modulus K, Lamé shear modulus μ, Poisson ratio ν, and HM coupling coefficients: the Biot coefficient B and the Biot modulus M. In all cases considered, we discuss the impact of the degree of fissuring and fracturing on the upscaled stiffness and hydromechanical coefficients. * C ijkl c ( ) in terms of the fracture normal vector components n i and the Kronecker delta δ ij , we obtain:

Introduction

Understanding the properties of fractured rock masses is fundamental for many industrial applications and in particular for radioactive waste disposal in deep underground facilities. The objective of this paper is to develop a methodology for calculating the upscaled or homogenized 3D coupled Hydro-Mechanical ("HM") properties of a damaged and fractured porous rock. This study will focus on the case of water-saturated elastic rock, with compressible water, deformable porous matrix, and deformable cracks (or fractures).

We aim to demonstrate the methodology by applying it, first, to a 3D isotropic synthetic sample of fractured porous rock, and then, to the radially inhomogeneous damaged zone around a deep gallery excavated in claystone. For the real case application, we consider a stretch of the experimental "GMR" gallery, located in the Meuse/Haute-Marne Underground Research Laboratory (M/HM URL), operated by ANDRA, the French national agency for the management of radioactive waste. This gallery is located in the 155 million years old Callovo-Oxfordian formation, a thick 130 m claystone layer located between depths 400 m and 600 m.

The hydro-mechanical upscaling method used here is basically a superposition method, enhanced to take into account pressure/stress coupling and porous matrix/cracks coupling. In this superposition method, he local strains are superimposed under a frozen "effective" stress field. Similarly, for hydraulics upscaling of the permeable, porous and fractured rock (previously treated in Ref. 1), the local fluxes are Coussy and others is useful for extending poro-elastic models to variably saturated porous media, where suction induces cohesive forces which depend on the shape of the saturation/suction curve. [START_REF] Coussy | [END_REF][START_REF] Dangla | A micromechanical approach to the behaviour of unsaturated porous media[END_REF][START_REF] Grgic | Effect of suction on the mechanical behaviour of iron ore rock[END_REF][START_REF] Lassabatère | Couplages hydromécaniques en milieu poreux non saturé avec changement de phase[END_REF] An overview on the relations between variables that can be measured in quasi-static laboratory experiments, base dona volume averaging approach, is presented in Ref. 23. In the anisotropic case, Pride and Berryman discuss the emergence of a volume averaged tensorial « strain measure » that depends on a weighted areal porosity which yields to the classical Biot and Willis relations.

The remainder of this paper is organized as follows. The 3D geometric structure of the damaged zone around a gallery, used for the field application, is described in Section 2. Section 3 develops the tensorial HM upscaling equations through the superposition approach, starting from a micro-scale description of strains for each component of the fractured rock. Section 4 describes the special case of isotropically fractured porous rock, and implements upscaling tests on synthetic fractured samples. Section 5 shows the various subdomain averaging procedures used to present the results, and then describes the in situ application and analyses the inhomogeneous spatial distribution of the upscaled coefficients thus obtained. The conclusive Section 6 recapitulates and discusses the results, and describes ongoing work towards extensions of those results.

Structure of fractured claystone around the GMR gallery

We consider a 3D stretch of the experimental "GMR" gallery, excavated, at depth ∼500 m, in the Meuse/Haute-Marne Underground Research Laboratory (M/HM URL). The GMR gallery is sub-horizontal, aligned with the direction of the in situ minor horizontal principal stress "σ h ". We assume a cylindrical shape for the drift, and we adopt the following parameters for the geometry of the domain: (i) axial length of the stretch: L = 20 m; (ii) drift diameter: 4 m; (iii) approximate thickness of the annular damaged zone or "EDZ" 1 around the drift: 4 m. The domain of investigation is a 3D rectangular box of size 20 m × 13 m × 13 m, which encloses the drift itself (centered and parallel to the first direction) and its Excavation damaged Zone ("EDZ"). The transverse scales 13 m × 13 m extend beyond the observed "EDZ" at the URL site.

To represent quantitatively the geometric structure of the damaged & fractured clay stone around the cylindrical excavation, we use a hybrid geometric model comprising two sets (Fig. 1): (i) a statistical set of 10000 planar disc-shaped isotropically oriented fissures, with a Poisson location of the centers, but with radially inhomogeneous crack density, radii and apertures (see Fig. 1a andc). Some of the statistical and geometric parameters were calibrated by comparing the upscaled permeabilities with measured permeabilities along radial boreholes; (ii) a deterministic set of large curved "chevron" fractures, periodically distributed along the axis of the gallery according to a 3D chevron or herringbone pattern (see Fig. 1b andc). We developed a modified conoïdal surface that resembles, under various projections, the observed traces of the large "chevron" fractures on gallery walls. For the purposes of upscaling, each curved fracture surface was further discretized into approximately 900 planar triangular pieces (patches). Thus, considering a cylindrical stretch of length 20 m, with a 0.50 m inter-spacing between chevron fractures, 46 curved fractures are needed to account for the fracturing along the gallery stretch, and 41400 planar objects were finally used to represent the chevron fractures in the upscaling procedure. For details on the complete set of geometric parameters of both fissures and fractures see Ref. 1. superimposed under a frozen pressure gradient field. The input parameters are the geometry of the fracture sets, the elastic stiffness moduli and apertures of the joints, and the elastic stiffness moduli and porosity of the rock matrix. The output parameters are the result of a system of upscaled continuum laws with equivalent continuum tensorial coefficients of two types: (i) mechanical stiffnesses and compliances (generalized 4th order tensorial Hooke/Terzaghi); and (ii) HM coupling coefficients and moduli (generalized Biot type laws). In fact, the nonisotropic non-orthotropic formulation yields two different coupling Biot coefficients, which become identical to each other only under some simplifying hypotheses. Results of the upscaling are presented using various averaging criteria and spatial supports (3D distributions, 2D crosssection distributions, radial profiles and global upscaled parameters).

Many of the available works in the literature deal with a granular porous microstructure, isotropic or not, with pore pressure coupling. The case of isotropic granular media is covered by the classical theory of Biot, [START_REF] Biot | General theory of three-dimensional consolidation[END_REF][START_REF] Biot | General solutions of the equations of elasticity and consolidation for a porous material[END_REF][START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF] which generalizes the empirical "effective stress" concept of Terzaghi [START_REF] Terzaghi | The shearing resistance of saturated soils and the angle between the planes of shear[END_REF] (pressure/stress coupling). There have been also various generalizations for anisotropic micro-structures, among other cases. However, we are interested here in the case of fractured porous materials, where the fine grained porous matrix is assumed isotropic but the fractures or "cracks" are geometrically anisotropic. Moreover, the anisotropic nature of cracks occurs both at the scale of individual cracks (strong local anisotropy due to the cracks aspect ratio), and at the larger scale of a rock sample comprising many cracks (due to anisotropic distribution of fractures orientations). In that respect, this paper will present a generalized version of the Terzaghi/Biot coupling, using a three-dimensional upscaling approach based on strain superposition, tailored to deal with the distinct geometric structures of the fine grained saturated porous matrix, and of the fluid-filled crack system.

Superposition methods have been used in the literature for obtaining macro-scale coefficients such as permeability, mechanical stiffnesses, and hydro-mechanical coefficients in fractured media under various restrictive hypotheses. Previous works [START_REF] Oda | An equivalent continuum model for coupled stress and fluid flow analysis in fractured rock masses[END_REF][START_REF] Ababou | Continuum modeling of coupled thermo-hydro-mechanical processes in fractured rock[END_REF][START_REF] Ababou | Coupled thermo-hydro-mechanical modeling for the near field benchmark test 3 (BMT3) of DECOVALEX phase 2 -progress report #[END_REF][START_REF] Stietel | Continuum representation of coupled hydro-mechanical processes of fractured media: homogenisation and parameter identification[END_REF] used a flux and strain superposition approach to develop upscaled hydraulic and HM properties of fractured rock assuming that the deformable matrix is impervious to flow and does not contribute to hydro-mechanical pressure/ stress coupling. Saevik et al. [START_REF] Saevik | A 3D computational study of effective medium methods applied to fractured media[END_REF] studies hydraulic upscaling for fractured media also with impervious matrix, where they consider either superposition methods (also called "effective medium approximation"), or various extensions of the Self Consistent approach. Ababou et al. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF] introduced a permeable porous matrix for permeability upscaling in 2D, and later Cañamón et al. [START_REF] Valera | Coupled Phenomena in 3D Fractured Media. Analysis and Modeling of Thermo-Hydro-Mechanical Processes[END_REF][START_REF] Cañamón | Numerical modeling and upscaling of a 3D fractured rock with Thermo-Hydro-Mechanical coupling[END_REF][START_REF] Cañamón | 3D fracture networks: optimal identification and reconstruction[END_REF] generalized that method for 3D hydromechanical upscaling, in which only the fluid-filled cracks contributed to H-M pressure/stress coupling. Later, [START_REF] Ababou | Equivalent upscaled hydro-mechanical properties of a damaged and fractured claystone around a gallery (Meuse/haute-marne underground Research laboratory)[END_REF] extended the previous method so that both the fluid-filled cracks and the saturated porous matrix contribute to Hydro-Mechanical pressure/stress coupling (and to fluid production/pressure coupling). The present paper develops in detail and validates the full mathematical theory and tensorial equations of the upscaling technique used in Ref. 14.

Some tensorial forms of the hydro-mechanical laws for equivalent fractured porous media have already been used in the literature, albeit with restrictive hypotheses. Thus, Cheng [START_REF] Cheng | Material coefficients of anisotropic poroelasticity[END_REF]16 formulates the set of fully tensorial non-isotropic non-orthotropic laws and coefficients for stressstrain elasticity, with pressure coupling. These laws are typically applicable to materials having a non-isotropic microstructure. These works also presents specializations of these laws for isotropic materials, analyzes the relations between scalar coefficients, and presents convenient formulations for experimental identification of these coefficients under undrained conditions. In particular, the relation between drained & undrained formulations of poro-elasticity is clearly presented. Previously, Coussy [START_REF] Coussy | Mécanique des Milieux Poreux[END_REF][START_REF] Coussy | Mechanics of Porous Continua[END_REF][START_REF] Coussy | [END_REF] have developed a formal thermodynamic approach to poro-elasticity for granular media, which helps explain various coupling terms, e.g., the fluid production caused by H-M coupling. Finally, we note that the thermodynamic approach of 3. Upscaling theory for hydro-mechanics: strain superposition method and upscaled coefficients 3.1. Equivalent "macro-scale" continuum equations, tensorial coefficients, and state variables

The upscaled behavior of the claystone is described by a system of coupled equivalent continuum hydro-mechanical PDE's which enforce mass conservation and (quasi-static) momentum conservation. The state variables are total stress (σ ij or Δσ ij )(normal stresses are taken negative under compression in this paper), total strain (ε ij ), and water pressure (p or Δp). Other intermediate variables are the water flux density vector (q i ) expressed by Darcy's law, [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF] and the increment in fluid content Δξ (also interpreted as the volumetric density of incremental fluid production), which is defined further below. The macroscale continuum equations to be further obtained in this paper involve the following tensorial laws and coefficients:

Firstly, hydro-mechanical Hooke/Biot stress-strain-pressure law, called here the "first Biot constitutive equation":

= - σ R ∆ε B ∆p ij ijkl kl ij I ( ) (1)
The reciprocal law expresses strain vs. stress and water pressure as follows:

= + ε C ∆σ B ∆p īj ijkl kl ij (2)
where:

R ijkl [Pa] is the stiffness tensor. B ij (I) is a dimensionless tensorial "Biot" coupling coefficient, called here "first" Biot coefficient.

C ijkl [Pa -1 ] is the compliance tensor, and. B ¯ij [Pa -1 ] is a "reciprocal" Biot coefficient representing strain/ pressure coupling (known as the poro-elastic expansion coefficient in the geotechnical literature).

Secondly, hydro-mechanical equation, which relates water production via the increment in fluid content Δξ [m 3 of water/m 3 of clay rock], to water pressure variation (Δp) and to the global strain (ε ij ) of the poroelastic continuum. This law is called here the "second Biot constitutive equation":

= + ∆ξ B ε M ∆p 1 ij II ij II ( ) ( ) (3)
where:

B ij (II) is a dimensionless tensorial "Biot" coupling coefficient called here "second" Biot coefficient (see further below), and.

M (II) is the dimensionless Biot modulus, called here "second" Biot modulus (see further below).

And finally, the mass conservation equation for the moving fluid:

∂ ∂ = - ∂ ∂ m t x ρ Q { } w i w i (4)
where.

= ∆m ρ ∆ξ W W
is the incremental variation (or production) of water mass.

ρ W [kg/m 3 ] is the water density, and.

= Q K J i ij j [(m 3 /s)/m 2 ]
is the water flux density vector, computed via a macroscale matrix/fractures Darcy's law, i.e., multiplying the upscaled equivalent hydraulic conductivity tensor, K ij , by the global hydraulic gradient J j .

In the remainder of this work, we focus on mechanics and coupled hydro-mechanics upscaling.

Underlying hypotheses of the upscaling

The medium to be homogenized comprises an isotropic elastic porous matrix and a set of elastic planar cracks with a given constant aperture "a". The crack aperture is assumed much smaller than its other length scales, but the cracks can have otherwise various geometries. In total, all the cracks and pieces of fractures taken into account in this work are planar convex sets, like discs, or triangles, or their intersections with the convex homogenization domain. The diameter or characteristic length "ℓ" of a planar convex set -like a polygonal or an elliptic shaped crack -can be defined in different ways. In mathematical morphology, the diameter of a planar convex set is the largest distance between pairs of parallel lines tangent to its boundary, but it can also be defined as some equivalent disc diameter, e.g., d = 2R = 4 × Area/ Perimeter, or alternatively, d = 2R = 2(Area/π) 1/2 . The latter expression represents the equivalent diameter of a circle having the same area as the convex planar crack, and it is our choice in this work.

We take into account the "internal" H-M coupling within the porous matrix, as well as the H-M coupling due to the cracks, and the nontrivial H-M interaction between pressure couplings in the cracks and in the porous matrix. Therefore, the developments presented here can be interpreted as a theory of "dual continuum upscaling", provided the assumption of local hydro-mechanical equilibrium between the porous matrix and the cracks. Their relative contributions to the H-M coupling will be analyzed by decomposing formally the Biot coefficient B ij and the Biot modulus M of the global matrix/cracks system. This analysis based on crack/matrix decomposition of coupling coefficients will be easier to understand in the special case of isotropic systems of cracks, a case which will be developed in detail for interpretation.

The main physical and mathematical hypotheses considered in the local scale laws are: (i) low Reynolds number for water flow (Darcy & Poiseuille); (ii) quasi-elastic deformations (here, calculations for elastic case); (iii) quasi-static mechanics (no shocks and no wave dynamics); (iv) isotropic homogeneous poro-elastic rock matrix (Biot coefficient B M ) and; (v) non-isotropic elastic cracks: specific normal stiffness k N * ≠ specific shear stiffness k S * [Pa/m], with Terzaghi's law assumed locally for each crack (B C = 1). The main hypotheses of the upscaling methodology itself are: (i) linear superposition of fluxes q i for hydraulic upscaling; (ii) linear superposition of strains ε ij for mechanical upscaling; (iii) two different cases are considered, (iii.a) no H-M coupling within the porous matrix (B M = 0, M M = ∞, θ M = 0) and (iii.b) full contribution of porous matrix & cracks to HM coupling (B M ≠0, M M ≠∞, θ M ≠0).

3.3. Expressions of local 4th rank compliance tensors for the porous matrix, and for each individual crack

Porous matrix compliance tensor

The intact porous matrix is assumed to follow an isotropic elastic behaviour. The intrinsic compliance tensor of the isotropic matrix can be expressed in terms of the matrix intrinsic Young modulus, E M *, and Poisson ratio, ν M *, as follows:

= + - * * * * * C ν E δ ν E δ δ 1 ijkl M M M ijkl M M ij kl (5)
with contractions

= - = * * * * C δ ν E δ K (1 2 )/ /(3 ) ijkk M ij M M ij M , = * * C K 1/ iikk M M , where = + δ δ δ δ δ ( )/2 ijkl ik jl jk il
is the 4th order identitity tensor and K M * is the intrinsic bulk modulus of the matrix. However, in the sequel we will keep the general tensorial notation of * C ijkl M for the sake of generality, and it will only be substituted by its isotropic form at the end. The same isotropic behaviour is assumed for the solid grains of the matrix, and an equivalent definition holds for the solid grains intrinsic compliance, * C ijkl S .

Single fracture compliance tensor

Fig. 2 shows the local e k ' and global e k coordinated systems criteria letting to express the fracture orientation in 2D (Fig. 2a, k = 1, 2) or 3D (Fig. 2b, k = 1, 2, 3), where the crack normal vector is defined as n.

The intrinsic compliance tensor of a single planar crack "c", expressed in its local 3D coordinated system e k ', is as follows:

= + + * * * C δ δ δ δ C δ δ C ( ) 4
where. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF] is the shear compliance of the crack "c" (in the plane e 1 '-e 2 ').

= * * - C a k ( ) N c N ( ) 1 is the normal compliance of the crack "c" (along axis e 3 '). = * * - C a k (2 ) S c S ( )
This can be written using non-equilibrated indicial Kelvin notation, which represents 4th rank symmetric tensors via 6 × 6 matrices:

= ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ * * * * C C C C ( ) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 IJ c LOCAL N S S ( ) (7)
In order to apply a superposition approach for the strains of the fractured medium, we need to rotate the * C ijkl c ( ) of every fracture from the local coordinated system e k ' to a fixed global coordinated system e k with the aid of the rotation matrix A ij relating both vector bases. After some algebraic simplifications, and expressing the compliance tensor

= - + * * * * C C C F C G ( 4 ) 4 ijkl c N S ijkl c S ijkl c ( ) ( ) ( ) (8)
where F ijkl (c) and G ijkl (c) are geometric or fabric tensors associated to crack (c):

⎧ ⎨ ⎩ = = + + + F n n n n G δ F δ F δ F δ F ( ) ijkl c i c j c k c l c ijkl c ik jl c il jk c jk il c jl ik c ( ) ( ) ( ) ( ) ( ) ( ) 1 4 ( ) ( ) ( ) ( ) (9)
with contractions The superposition approach is applied to the strains produced in the porous fractured system, within any shape for the homogenization subdomain. This is accomplished by adding the strain contribution of each individual fracture inside the subdomain to a "total strain due to cracks", and adding to it the local strain produced in the intact porous matrix, weighted by their respective volumetric fractions. To compute the local strains of cracks and matrix we work under the hypothesis of a "frozen" global stress tensor Δσ ij and a frozen global fluid pressure Δp. They can be combined to form a frozen effective stress Δσ ij (EFF) that is formulated differently for matrix and cracks, using the local coupled hydro-mechanical laws for each medium. Although we assume Terzaghi behaviour for all cracks (B C * = B C = 1), we will keep track of the corresponding isotropic Biot coefficient B C δ ij in the formulation below for the sake of completeness. The superposition approach yields:

= = = G F F n n ijkk c ijkk c ij c i c j c ( ) ( ) ( ) ( ) ( ) , = = G F 1 iikk c iikk c ( ) ( ) 3
= + ∑ = - + ∑ = = - -- + ∑ -- = = - + ∑ + - + ∑ = * = * * = * * = * * = * ε ε ε Φ C ∆σ φ C ∆σ Φ C ∆σ B δ ∆p φ C ∆σ B δ ∆p Φ C φ C ∆ σ B Φ C B φ C ∆p (1 ) (1 ) [ ( )] [ ( )]
(

)

ij ij M c N ij c ijkl M kl M EFF c N c ijkl c kl C EFF ijkl M kl M kl c N c ijkl c kl C kl ijkl M c N c ijkl c kl M ijkk M C c N c ijkk c 1 ( ) ( ) 1 ( ) ( ) ( ) 1 ( ) ( ) 1 ( ) ( ) 1 ( ) ( ) (10)
where:

ϕ (c)
is the volumetric fraction of fracture or crack (c) in the homogenization domain;

= ∑ = Φ φ c N c 1 ( ) is the volumetric fraction of N fractures in the homogenization domain;
This equation (reciprocal form of the first Biot constitutive relation) can be compared with eq. ( 2) to identify the coupling coefficients involved on it.

The equivalent compliance tensor of the matrix/cracks system at the scale of the homogenization subdomain is the coefficient multiplying Δσ kl in eq. ( 10):

∑ = - + = + * = * C Φ C φ C C C (1
)

ijkl ijkl M c N c ijkl c ijkl M ijkl C 1 ( ) ( ) (11)
where we can distinguish the equivalent compliance tensor associated to the matrix, C ijkl M , and the equivalent compliance tensor associated to the set of cracks, C ijkl C . The last one can be further developed using eq. ( 8) and the definition of the specific surface area [m 2 /m 3 ] contributed by each individual crack (c),

= σ φ a 2 / c c c ( ) ( ) ( ) : ∑ ∑ ⎜ ⎟ ⎜ ⎟ = = ⎧ ⎨ ⎩ ⎛ ⎝ - ⎞ ⎠ + ⎫ ⎬ ⎭ = ⎛ ⎝ - ⎞ ⎠ + = * = * * * * * * C φ C σ k k F k G k k F k G 2 1 2 2 1 2 2 ijkl C c N c ijkl c c N c N S ijkl c S ijkl c N S ijkl S ijkl 1 ( ) ( ) 1
( ) ( ) ( ) (12)
where we have defined, similarly to Ref. 6 and others, two "geometric" or "fabric" tensors, characterizing the geometry of the fracture network:

⎧ ⎨ ⎩ = ∑ = ∑ = ∑ = + + + = = = F F n n n n G G δ F δ F δ F δ F ( ) ijkl c N σ ijkl c c N σ i c j c k c l c ijkl c N σ ijkl c ik jl il jk jk il jl ik 1 2 ( ) 1 2
( ) ( ) ( ) ( ) 1 2 ( ) 1 4 c c c ( ) ( ) ( ) (13)
with contractions = = G F F ijkk ijkk ij (compared to Oda's equation ( 21), there appears here a factor of 2 in the shear compliance terms of C ijkl C ). The equivalent stiffness tensor (T ijkl or R ijkl ) can be obtained by inverting the equivalent upscaled compliance tensor given in eq. ( 11) via equilibrated 6x6 matrix notation for tensorts as defined in Ref. 25. The equilibrated version of the 6x6 compliance matrix preserves the tensorial properties of the inversion.

The equivalent poro-elastic expansion coefficient, denoted here B ¯ij , of the matrix/cracks system at the scale of the homogenization subdomain is the coefficient that multiplies Δp in eq. ( 10):

∑ = - + = + * = * B B Φ C B φ C B B ¯(1 ) ¯īj M ijkk M C c N c ijkk c ij M ij C 1 ( ) ( ) (14)

Definition of the 1st Biot coefficient B ij (I)

The inversion of the strain-stress law of eq. ( 2) yields a relation for the "first" tensorial Biot coefficient which is

= B R B īj I ijkl kl ( )
. We can now develop this expression by inserting eq. ( 14), obtaining:

⎟ = + = ⎛ ⎝ + - ⎞ ⎠ * * B R B C B C R B k F Φ B K δ ( ) (1 ) 3 ij I ijmn C mnkk C M mnkk M ijmn C N mn M M mn ( ) (15)
The first Biot coefficient B ij (I) only depends on the specific normal stiffness k N * of the cracks and on the specific bulk stiffness modulus of the matrix (K M *).

If no Biot coupling is considered for the matrix, i.e., B M = 0, then the first Biot coefficient (labelled in this case B ij (0)

) represents solely the contribution of cracks to the stress/pressure coupling. Furthermore, its expression can be compared to the expression obtained by [START_REF] Oda | An equivalent continuum model for coupled stress and fluid flow analysis in fractured rock masses[END_REF] in the case of 2D fractured media with non-porous elastic matrix, given Terzaghi's hypothesis (B C = 1). Therefore, eq. ( 15) completes Oda's definition of the Biot coefficient, taking into account the matrix stress/pressure coupling. Let us consider a material control volume V of equivalent fractured poro-elastic continuum medium, of initial volume V 0 . Fig. 3 shows a schema of the conceptual model presented in two different ways: both pores and cracks of Fig. 3a constitute the so-called "pore volume", which is completely filled with water at any time. The pore volume at the initial state is represented in white in Fig. 3b, and the water volume filling that pore space is shown separately in blue. Water undergoes pressure changes and hydro-mechanical coupling, not only through the cracks but also through the pores of the matrix. The total porosity Θ of the porous matrix + cracks system incorporates the contributions of matrix porosity θ M and cracks volumetric fraction Φ:

= + - Θ Φ θ Φ (1 ) M ( 16
)
The control volume V can be decomposed in pore volume and solid volume, and the pore volume itself is formed by the matrix porosity volume and the cracks volume:

= + = + + = + - + - - V V V V V V ΦV θ Φ V θ Φ V (1 ) (1 ) (1 )
≡ - ∆ξ ∆V V / ξ ( 18
)
where ∆V ξ is the volume of fluid entering (negative) or exiting (positive) the control volume due to pore pressure variations. Alternatively, Δξ could be defined as = ∆ξ ∆m ρ / W W , where Δm W represents the variation of the specific fluid mass entering in the control volume and ρ w ≈ρ w0 due to small deformation hypothesis. It was firstly defined by Biot & Willis 4 as the "increment in fluid content", and it is positive when water flows into the control volume and negative when water flows out of the control volume. We will formulate the fluid production term by volume variations, following a volume decomposition approach similar to that presented by, [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF][START_REF] Pride | Connecting theory to experiments in poroelasticity[END_REF] by analyzing volume variations for the control volume (ΔV), pore volume (ΔV Θ ), and sample fluid volume (ΔV W ) due to increments of confining pressure and/or fluid pressure with respect to an initial state {t = 0: V 0 , V Θ0 , V W0 }.

Term ΔV W represents the variation of the volume occupied by the total mass of fluid initially in the pore space (matrix pores and cracks) of the sample. The variation of fluid volume ΔV W will be at any time the sum of two terms:

= + ∆V ∆V ∆V W Θ ξ ( 19
)
We can manipulate equation (19) to obtain, equivalently.

- = ⎛ ⎝ ⎜ - ⎞ ⎠ ⎟ ∆V V Θ ∆V V ∆V V ξ Θ Θ p W VΘ ( 20
)
This equation is equivalent to equation #28 of, [START_REF] Pride | Connecting theory to experiments in poroelasticity[END_REF] but provides a precise interpretation of fluid production in terms of the relative variations at constant p or constant V Θ . The first term corresponds to the pore volume variation (at constant pressure). The total volume of pores is V Θ = V Θ, and the relative variation of pore volume is, under the assumption of small deformations:

= + ⇒ = + ∆V Θ∆V V∆Θ ∆V V ∆V V ∆Θ Θ Θ Θ Θ ( 21
)
Variation of total porosity ΔΘ can be further developed, taking increments in eq. ( 16):

= - + - ∆Θ Φ ∆θ θ ∆Φ (1 ) (1 )
M M (22) leading finally to:

= + - + - ∆V V ∆V V Φ ∆θ Θ θ ∆Φ Θ (1 ) (1 )
Θ Θ M M (23)
The second term of equation ( 20) corresponds to the variation of the fluid volume initially in the sample (at constant pore volume), which is obtained as follows: )

⎜ ⎟ = + ≈ ⎛ ⎝ - ⎞ ⎠ ⇒ ≈ ∆V V ∆V V ∆V ∆V V ∆V V ∆V V Θ ∆V V 1 W W W W W W 0 0 0 0 ( 24
)
where we have used:

+ = -+ ≈ - ε ε O ε ε 1/(1 ) 1 ( ) 1 2 if ε « 1, with ε
≡ ΔV/V 0 , we have neglected mixed products of variations such as (ΔV W /V 0 ) × (ΔV/V 0 ) and we have used, for constant porosity volume, the relation V 0 = V Θ0 /Θ = V W0 /Θ. Introducing expressions ( 23) and ( 24) into equation ( 20), we finally have:

- = = + - + - - ∆V V ∆ξ Θ ∆V V Φ ∆θ θ ∆Φ Θ ∆V V (1 ) (1 )
ξ p M M W W V 0 Θ (25)
Equation ( 25) can be rewritten in a shorter and more intuitive way, by identifying each term of the sum with a particular contribution of the medium to the fluid production:

= + + + ∆ξ ∆ξ ∆ξ ∆ξ ∆ξ SK M C W ( 26
)
To illustrate the previous expressions and the sign of the fluid production term, Fig. 4 shows some examples of stress "tests" (thought experiments) over a generic sample. By definition, the volume fraction ∆V V / represents the mean vo- lumetric strain of the poro-elastic medium or poro-elastic skeleton (to be distinguished from the strain of liquid water), which includes changes of porosity as well as changes of solid density (grains density). This can be computed by introducing the volumetric (spherical) part of the global strain tensor ε ij :

= = ∆V V trace ε ε ( ) kk (27)
Therefore, we can finally write:

= ∆ξ Θε SK kk ( 28
)
FLUID PRODUCTION TERM Δξ W DUE TO FLUID DEFORMATION.

For the fluid mass initially inside the control volume, the following equivalency holds:

= - = + = * ∆ξ Θ ∆V V Θ ∆ρ ρ Θ K ∆p 1 W W W W W W 0 ( 29
)
where Θ 0 ≈Θ and ρ w ≈ρ w0 because of the small deformations hypothesis and an isotropic elastic deformation law is adopted for the fluid, characterizing the fluid compressibility by an intrinsic spherical stiffness modulus K W *.

FLUID PRODUCTION TERM Δξ M DUE TO MATRIX POROSITY VARIATION.

By definition, the matrix porosity variation is equal to:

⎜ ⎟ = ⎛ ⎝ ⎞ ⎠ ∆θ ∆ V V M θM M ( 30
)
Developing the increments, and

introducing identities = - ∆V ∆V ∆V θM M S and, = - V V θ / 1 S M M yields: = - - ⇒ = - ⎧ ⎨ ⎩ - ⎫ ⎬ ⎭ ∆θ V ∆V ∆V V V ∆V V ∆θ θ ∆V V ∆V V 1 ( ) (1 )
M M M S θ M M M M M M M M S S ( 31
)
The two fractions in brackets can be identified with the intrinsic volumetric strains of the matrix and the solid grains, respectively, de-

= -

- + - = - - - + * * * * * * ∆ξ Φ θ C ∆σ B C ∆p C ∆σ Φ θ C C ∆σ B C ∆p (1 )(1 ){ } (1 )(1 ){( ) } M M kkmn M mn M kkll M kkmn S mn M kkmn M kkmn S mn M kkll M (34)
We can express the global σ ij in terms of the global strain ε ij and pressure p with the help of eq. ( 1), and using the relation

- - = - Φ θ Θ (1 ) (1 ) 1
M and regrouping terms yields finally:

= - - + - + * * * * ∆ξ Θ C C R ε B δ B C B C ∆p (1 ) {( ) [( ) ] } M kkmn M kkmn S mnij ij M mn mn I kkmn M mn I kkmn S ( ) ( ) (35)
FLUID PRODUCTION TERM Δξ C DUE TO VARIATION OF CRACKS VOLUMETRIC FRACTION.

By analogy with the matrix porosity variations developed previously, we can define the variation of cracks volumetric fraction ΔΦ as follows:

= ⎛ ⎝ ⎞ ⎠ ∆Φ ∆ V V ΘC (36)
And, after further developing the increments, we obtain:

= - = - ∆Φ ∆V V V V ∆V V ∆V V Φ ∆V V ΘC ΘC ΘC (37)
The second fraction of the right-hand side is, according to (27), the volumetric strain ε kk of the poro-elastic continuum. The first fraction corresponds to the global (not intrinsic) volumetric strain ε kk C of the cracks system:

= + ε C ∆σ B ∆p īj C ijmn C mn ij C ( 38
)
where C ijmn C and B ¯ij C are defined, respectively, by eqs. ( 12) and ( 14), from the superposition approach.

Inserting now the equivalent volumetric strain law of the cracks system in eq. (37), and multiplying by (1-θ M ) we obtain:

= - + - ∆ξ θ C ∆σ B C ∆p Φε (1 ){ } C M kkmn C mn C kkll C kk
where the global stress σ mn can be expressed in terms of the Hooke-Biot law for the cracks.

Finally, by regrouping terms, we get the fluid production term due to the variation of cracks volumetric fraction:

= - - + - ∆ξ θ C R ε Φε B δ B C ∆p (1 ) { ( ) } C M kkmn C mnij ij kk C mn mn I kkmn C ( ) (39)
FINAL EXPRESSION OF FLUID PRODUCTION Δξ. Substituting the four terms developed earlier in equation ( 26), applying the relation (1-θ M ) Φ = (Θ-θ M ), and rearranging deformation and pressure terms separately, we obtain finally the expression for fluid production taking all phenomena into account:

= - - + - + + + - - + + + - - * * * * * { } ∆ξ Θ C C R θ C R θ δ ε Θ B δ B C C B Θ θ B δ B C ∆p {(1 )( ) (1 ) } (1 )[( ) ] (1 )
( ) The coefficient multiplying the deformation term ε ij of eq. ( 40) has been identified as the "second Biot coefficient" B ij (II) = (∂ξ/∂ε ij ) p=0 , as denoted in equation (3). It expresses the fluid production per unit expansion of the dual matrix/cracks poro-elastic medium under drained conditions (p = 0 or more generally p = p 0 ). Its inverse {B (II) } -1 can be interpreted as a drained coefficient of poro-elastic expansion per unit injection of water {B (II) } -1 ij = (∂ε ij /∂ξ) p=0 . We can rearrange the terms by separating the factors multiplied by matrix porosity θ M :

kkmn M kkmn S mnij M kkmn C mnij M ij ij M mn mn I kkmn M kkmn S mn I K M C mn mn I kkmn C ( ) ( ) 1
= - - + + -- - - * * * * B Φ C C R C R θ δ Φ C C R C R (1 )( ) { (1 )( ) } ij II kkmn M kkmn S mnij kkmn C mnij M ij kkmn M kkmn S mnij kkmn C mnij ( ) (41)
Particularized for the case of isotropic matrix and isotropic solid grains, and using the classical expression for the scalar Biot coefficient B M coupling stress/pressure inside the isotropic porous matrix, B M = 1-K M */K S *, yields:

= - + + ⎧ ⎨ ⎩ -- - ⎫ ⎬ ⎭ * * B Φ B K R R C θ δ Φ B K R R C (1 ) 3 (1 ) 3 ij II M M kkij ijmn mnkk C M ij M M kkij ijmn mnkk C ( ) (42)
This finally yields a relation between the second and first Biot coefficients when Terzagui's law (B C = 1) holds for the cracks:

= + - B B θ δ B { } ij II ij I M ij ij I ( ) ( ) ( ) (43)
The first and second Biot coefficients become equal under the two hypotheses B C = 1 and θ M = 0. Note: if matrix porosity is neglected, the sole contribution to total porosity is the volumetric fraction of cracks.

Definition of 1st and 2nd Biot moduli M (I) and M (II)

The coefficient multiplying the pressure term Δp in eq. ( 40) is defined as the inverse of the generalized Biot modulus, 1/M (II) = (∂ξ/ ∂p) εij=0 , whence M (II) = (∂p/∂ξ) εij=0 . The latter may be interpreted as the "stiffness" of the pressure/production coupling, while its inverse 1/ M (II) is the "compliance" [Pa -1 ] of the production/pressure coupling. Fluid production is expressed in m 3 of water per m 3 of sample, and the sample is mechanically confined. The expression of M (II) can be rearranged by separating the factors multiplied by matrix porosity θ M (as was done for the second Biot coefficient). This leads to:

= - - + + - + + -- - + - - - * * * * * * * { } { } Φ B δ B C C B B δ B C Φ θ Φ B δ B C C B B δ B C Φ (1 ) [( ) ] ( ) (1 )
This expression can be cast equivalently in the form:

⎜ ⎟ = + ⎛ ⎝ - ⎞ ⎠ * M M θ K M 1 1 1 1 II I M W I ( ) ( ) ( ) (45
where the "first Biot modulus" M (I) has been defined as:

= - + - + - + * * M B δ B C B δ B C Φ C B Φ K 1 ( ) ( ) (1 )
1

I M mn mn I kkmn M C mn mn I kkmn C kkmn S mn I W ( ) ( ) ( ) ( ) (46)
The two Biot moduli, M (I) and M (II) , coïncide only in the special case where θ M = 0. Clearly, the 2nd Biot modulus M (II) appears as a generalization of the 1 rst modulus M (I) as the first modulus does not incorporate fluid production coupling that takes place within the pores of the rock matrix, and therefore, it does not fully represent coupling unless the rock matrix porosity can be neglected.

Specialization of H-M coupling coefficients for locally isotropic porous matrix and solid grains

The formal expressions of Hydro-Mechanical coupling coefficients can be specialized for the case of locally isotropic porous matrix and locally isotropic solid grains, which is the hypothesis considered in the applications of this study. The first Biot coefficient B (I) , expressed in eq. (15), can be further developed by inserting the isotropic compliance tensor of the porous matrix (from eq. ( 5)), This yields:

⎜ ⎟ = - ⎛ ⎝ - ⎞ ⎠ + * * B Φ K K R B R C (1 ) 1 3 1 1 ij I M S ijkk C ijmn mnkk C ( ) (47)
Equation ( 46) yields for the case of isotropic matrix and solid grains:

= - - + - + * * M Φ B K δ δ B B C δ B B Φ K 1 (1 ) 3 ( ) ( / ) 1 I M M mn mn mn I C kkmn C mn mn I C W ( ) ( ) ( ) (48)
which can then be introduced in expression (45) to obtain 1/M(II).

To summarize the various upscaled laws and coefficients obtained so far, we present in Fig. 5 a flow chart of the full set of HM equations and how are they related to each other.

Specialization for statistically isotropic cracks (and isotropic matrix)

This numerical exercise is aimed at validating the upscaling superposition approach for perfectly isotropic fractured rock, with isotropic matrix and solid grains but also with fractures randomly oriented in 3D space and spatially homogeneous (in a statistical sense). Note that the case of isotropic matrix-cracks systems can arise in applications, but it is also interesting for its own sake, as the resulting H-M coefficients become easier to interpret and analyse. Note, also, that the isotropic theory in this section is applicable either to a truly isotropic fractured rock, or else to an "isotropized" fractured medium deduced from the real fractured medium.

Several tests of numerical upscaling on a perfectly isotropic system (cartesian network) and on quasi-isotropic matrix-crack systems are analyzed. The purpose is to apply the fully tensorial numerical upscaling procedure to those isotropically fractured porous rock, and to compare the numerically upscaled coefficients with the theoretical isotropic expressions to be developed in the next two subsections. Finally, we also want to verify that asymptotic values of the stiffness ratio r = K C /K M lead to physically correct values of the upscaled Biot coefficient B and Biot modulus M. Fig. 5. Flowchart of equivalent continuum constitutive laws and coefficients, as obtained from the Hydro-Mechanical upscaling procedure developed in the text, for locally isotropic matrix and Terzaghi behavior in the cracks.

= + B B K K B K K I C bulk C M bulk M ( ) (49) = + - B B θ B (1
)

II I M I ( ) ( ) ( ) (50)
where

= - K C ( ) C iikk C
1 is the isotropic stiffness modulus of the crack system,

= - * K K Φ /(1 ) M M
is the isotropic stiffness modulus of the matrix and = + -

K K K (1/ 1/ ) bulk C M
1 is the bulk stiffness modulus of the fractured medium. Equation ( 43) would give the scalar second Biot coefficient B (II) .

The original expression of Oda [START_REF] Oda | An equivalent continuum model for coupled stress and fluid flow analysis in fractured rock masses[END_REF] for the equivalent Biot coefficient is obtained as a particular case of (49) by considering the following hypotheses: (i) Terzaghi behaviour on all cracks (B C = 1); (ii) absence of stress-pressure coupling in the matrix (B M = 0); and (iii) constant fracture diameter (setting = ℓ ℓ f in Oda's expressions):

= = * B B K K ρ K k C bulk C bulk N (0) 32 ( 51
)
where the relation

= * K k ρ / C
N 32 has been used, and ρ 32 is the crack density in m 2 /m 3 (specific area, defined as the sum of fracture areas per unit volume of the fractured medium).

To check the asymptotic behavior of the scalar Biot coefficients B (0) , B (I) and B (II) we express them in terms of a stiffness ratio defined as = r K K / C M , and we examine the results for r ≫ 1 and r « 1:

= = = = > > ⇒ → < < ⇒ → = ′ + + + B B K K B K K B B if 0 if 1 (Terzaghi s hypothesis) C K K B K K K B K K B r C M C M C (0) (1 / 1 / ) 1 / 1 (0) (0) bulk C C C C M C C M C (52) = + = + = + > > ⇒ → < < ⇒ → = + + + + - B B B K K B B K K B B if if 1 I C K K M K K B K K B K K B r B r C M I M C M I C ( ) 1 / 1 / 1 1 ( ) ( ) bulk C bulk M C C M M M C C M 1 (53) = + - > > ⇒ → + - < < ⇒ → + - = B B θ B K K B B θ B K K B B θ B (1 ) if (1 ) if (1 ) 1
II I M I C M II M M M C M II C M C ( ) ( ) ( ) ( ) ( ) (54)
4.1.2. Biot moduli M (0) , M (I) , M (II) for isotropic matrix/cracks systems A simplified expression can also be obtained from equations ( 45) and ( 46) for a homogeneous isotropic fractured medium, in which tensors can be reduced to simple scalars by tensor contraction. This leads finally to:

⎜ ⎟ = - + - ⎡ ⎣ ⎢ -⎛ ⎝ - ⎞ ⎠ ⎤ ⎦ ⎥ + * * * * M B B K Φ B K K K B Φ K 1 (1 ) 1 1 1 I C I C M M M S I W ( ) ( ) ( ) (55) = + - = = - - + -⎡ ⎣ - - ⎤ ⎦ + * * * * * ( ) ( ) θ θ B B Θ B Θ (1 ) ( ) (1 ) M M M K M M K C I B K K K I K 1 1 1 1 1 ( ) 1 1 ( ) 1
II I W I C M M M S W ( ) ( ) ( ) (56)
The first term in the right hand side of the latter equation ( 56) is the contribution of cracks to the second Biot modulus, the second term is the contribution of matrix, and the third term is the contribution of water (compressibility). For the case of cracks of constant aperture a, we could use the relation = ρ Φ a / 32 . Oda's 6 original expression for the equivalent Biot modulus is a particular case of equation ( 55) by considering the following hypotheses and simplifications: (i) Terzaghi behaviour for all cracks (B C = 1); (ii) absence of stress-pressure coupling within the matrix (B M = 0); (iii) no porosity within the matrix (θ M = 0); (iv) incompressible water (K W = ∞); and (v) constant fracture diameter (setting = ℓ ℓ f in Oda's expressions):

= - = - * M K B B ρ k B 1 1 ( ) ( 1
) C C N (0) (0) 23 (0) (57)
Contributions of matrix compressibility and of water compressibility are not taken into account in Ref. 6, and neither in Refs. 7-9. Thus, this specialization offers an excellent verification of the validity of our more general expressions (55) and (56).

As for Biot coefficients, it is also interesting to check the asymptotic behavior of the "isotropized" Biot moduli M (0) , M (I) and M (II) , with respect to the stiffness ratio r defined above. This is developed in the equations below:

= - = - = > > ⇒ → ⇒ → = < < ⇒ → ⇒ → = + + + * ( ) ( ) B B B B K K M K K K r M K ( ) if 1 if M ρ k C K C C K K B K r r C M r r K B C C M r r K B M 1 (0) 1 1 1 (0) 1 (0) N C bulk C C C C C M C (0) 23 (58) = + - + - + + = - + + > > ⇒ → - + < < ⇒ → - + + + + + - - - * - * * * * * ( ) ( ) ( )( ) ( ) Φ Φ Φ K K M Φ Φ K K M Φ Φ (1 ) if (1 ) if (1 )
M B K B K K K K B r B r K K B r B r K C M I B K K C M I B K K 1 1 1 1 1 1 1 1 1 1 1 ( ) 1 1 ( ) 1 1 I
C C M M C M S C M W S C M W M S W C S W ( ) 1 1 (59)
= - - + + + - = - + + > > ⇒ → - + < < ⇒ → - + + + + + - - * - * * - * * * * * ( ) ( ) ( ) ( ) θ Φ Φ θ Φ Θ Θ K K M Θ Θ K K M Θ Θ (1 ) (1 ) ( (1 )) (1 ) if (1 ) if (1 )
M M K B r B r M K K B r B r K C M II B K K C M II B K K 1 1 1 1 1 1 1 1 1 ( ) 1 1 ( ) 1 1
II S C M W S C M W M S W C S W ( ) 1 1 (60)
In the forthcoming numerical tests, we will compare the theoretical isotropic expressions given by the above equations with the numerically upscaled values obtained by the general numerical upscaling procedure implemented with our UP3D Matlab based software. The numerical tests parameters are set based on the values of the case study presented in section 6, corresponding to 10000 statistical fissures in the Excavated Damaged Zone around the GMR drift at the Meuse/Haute Marne URL. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF]

Isotropic test #1: symmetric cartesian fracture network

This test provides a perfectly isotropic behaviour for a porous fractured medium. We generate 3 sets of regularly spaced fractures orthogonal to the coordinated axes and completely traversing the 8x8x8m 3 generation domain (the upscaling domain is a sphere with the same volume). The number of fractures (2) and their spacing (5.39 m) for each set is selected such that we obtain a m 2 /m 3 fracture density ρ 32 close to that of the Excavation Damaged Zone micro-fissured system (see original drift dataset in Table 7): , B ij (II) for isotropic matrix/cracks systems A simplified expression can be obtained from equations ( 43) and (47) for a homogeneous isotropic fractured medium, in which tensors can be reduced to simple scalars by tensor contraction; this leads finally to:

= ∑ = = ∑ = ⋅ ≈ ⋅ = = - ρ A m m Φ V a ρ m m 0.6394 / 3.1968 10 / V f N f V f N f 32 1 1 2 3 1 1 32 5 3 3 (61)

Test#1/case#1: fractures stiffer than matrix

In this case (Test#1/Case#1), the fractures are stiffer than the matrix: K M =4.17 MPa and K C = 15.6 MPa (using k N * = 10 MPa/m, k S * = 1 MPa/m). The results are presented in Table 1.

Test#1/case#2: fractures more compliant than matrix

In this case (Test#1/Case#2), the cartesian fractures are more compliant than the matrix: K M =4.17 MPa and K C =1.56 MPa (with

k N * = 1 MPa/m, k S * = 0.1 MPa/m).
The results are presented in Table 2.

Test#1: asymptotic analysis with respect to stiffness ratio r = K C /K M

Finally, asymptotic analyses for the upscaled coupled coefficients of the first numerical test (Test#1) are summarized in Table 3.

We appreciate from both Tables 1 and2 a perfect agreement between the results of the isotropized equations and the direct upscaling for all the Biot coefficients and all the Biot moduli. Table 3 also shows that the asymptotic behaviour of all the coefficients is coherent with the isotropized equations expectations.

Isotropic Test#2: statistically isotropic 3D fractured rock

This test is performed on a cubic upscaling domain of (8 × 8 × 8) m 3 (in particular, in the equivalent volume sphere), where a statistically isotropic fracture network is generated. However, the generation domain has to be bigger than the upscaling domain in order to avoid decreasing fracture density near the borders ("border effect"). Therefore, a generation domain of (10x10x10) m 3 is selected. We set the following parameter values for the fractured network based on criteria also explained below: (i) the fracture centers (X fc , Y fc , Z fc ) are randomly generated following a homogeneous 3D Poisson process; (ii) orientations are statistically isotropic in 3D space, where the dip direction is a uniform random variable in [0,2π], and the cosine of the dip angle is a uniform random variable in [-1,+1]; (iii) fracture size is a constant value equal to the average radius of the 10000 statistical fissures of the GMR case study of Section 5:

∑ 〈 〉 = = ≈ = R N R m m 1 0.2005 0.20 i N i 1 (62)
(iv) In this test, the fractures have the same deterministic constant aperture of the GMR case study, = ⋅ - a m 5 10 5 and (v) the number of fractures (N = 5100) is set to fit approximately the areal fractured density (specific area) of the GMR case. The values of ρ 32 and Φ obtained for this test case are:

= ∑ ≈ = ∑ ≈ ⋅ = = - ρ A m m Φ V m m 0.6394 / 3.1971 10 / V f N f V f N f 32 1 1 2 3 1 1 5 3 3 (63)
These values are almost identical to the ones obtained in Test#1 with three orthogonal sets of fractures, which is very convenient for comparison purposes between Test#1 and Test#2.

The 5100 randomly generated 3D isotropic Poisson cracks are visualized in Fig. 6. The upscaling results are analyzed in the sequel for the same two cases of Test#1: Case#1 with cracks stiffer than the matrix, and Case#2 with cracks more compliant than the matrix. In addition, sensitivity analysis with respect to cracks/matrix stiffness ratio is also performed.

Test#2/case#1: cracks stiffer than matrix

In this case (Test#2/Case#1), the fractures are stiffer than the matrix: K M =4.17

Table 2

Upscaling results of the orthogonal fracture set of the isotropic Test #1/Case#2.

Biot coefficient

Oda's hyp. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF] General case [START_REF] Biot | General theory of three-dimensional consolidation[END_REF] B Oda's hyp. [START_REF] Biot | General solutions of the equations of elasticity and consolidation for a porous material[END_REF] General case [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF] M (0) M I M II) Table 4 Upscaling results for the set of 5100 isotropic cracks of Test #2/Case#1.

Biot coefficient Oda's hyp. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF] General case

Biot modulus [Mpa]

Oda's hyp. [START_REF] Biot | General solutions of the equations of elasticity and consolidation for a porous material[END_REF] General case [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF] M (0) M (I) Table 5 Upscaling results for the set of 5100 isotropic cracks of Test #2/Case#2.

Biot coefficient Oda's hyp. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF] General case [START_REF] Biot | General theory of three-dimensional consolidation[END_REF] B (0) Oda's hyp. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF] General case [START_REF] Biot | General theory of three-dimensional consolidation[END_REF] M (0) M (I) M (II) Isotropized equations 5.7307 9.6490 6.5460 Direct upscaling 5.7712 9.7056 6.5683

B (I) B (II)
Table 3 Asymptotic results of the isotropic Test #1.

described (partially) via two 3x3 matriceseach of which can be represented by ellipsoids-: (i) "Normal stiffness" matrix N ij containing stiffness coefficients R 1111 , …,R 3333, R 1122 , …,etc.; and "Shear stiffness" matrix S ij containing stiffness coefficients R 2323 , …,R 2312 , R 1313 , …,etc. Note that the shear matrix ellipsoïd is easier to interpret, as it reduces to a sphere in the case of isotropic shear stiffnesses, whereas the normal ellipsoïd does not reduce to a sphere in that case. In order to remedy some possible interpretation problems with these mechanical ellipsoids, we have complementarily extracted a set of more "classical" mechanical coefficients (E, K, λ, μ, ν) from the full 4th rank compliance and stiffness tensors. Similarly for the HM coupling coef-

ficients (B ij (II)
, M (II) ), which become (B, M) after isotropization. These calculations were performed using adequate contractions of 2nd rank and 4th rank tensors, based on the definition of spherical quantities similar to the definition of a spherical bulk stiffness K in isotropic elasticity.

Upscaling results around the excavated drift

3D spatial distributions of upscaled coefficients

The 3D upscaled tensorial coefficients around the drift are represented, as explained before, in the form of spatially distributed ellipsoïds on a Cartesian grid of voxels of 0.50x0.50 × 0.50 m 3 . As displaying the spatial distribution over the entire domain of 20 x 13 × 13 m 3 could be cumbersome, we show only the first transverse layer of the 3D distribution, i.e., the (Y,Z) plane with X ∈ [-10.0 m, -9.5 m]. This gives an idea of the "3D heterogeneity" of upscaled coefficients inside the damaged/fractured zone around the drift, as displayed in Fig. 8 for several types of coefficients.

Axial upscaling: 2D distribution transverse to the drift

A sequential upscaling averaging along the drift axis from the 3D distribution of upscaled coefficients is performed. The 2D spatial distribution of axially upscaled tensorial coefficients is shown in Fig. 9. In comparison with the 3D spatial distributions shown earlier, it can be seen that the 2D upscaled coefficients appear less anisotropic, and less variable in space (particularly near the wall). This effect is due to the additional smoothing induced by averaging along the axis of the drift.

Global upscaling: equivalent coefficients for the entire EDZ volume

The "global" upscaling domain, in this case, is the 20 m long cylindrical annular region 4 m thick around the drift known as EDZ, recalling that the radius of the drift itself is R = 2 m. Again, a sequential upscaling has been performed and averaged values have been In particular, Biot modulus can be interpreted as the stiffness of the coupling between pressure variations and fluid production. This global "stiffness", as expected, is reduced as more cracks are "added".

Similarly, as expected, the Biot coefficient B, measuring the interaction between fluid pressure and stresses, increases as fractured network becomes denser, as there is more fluid, considering fully saturated conditions, inside the poroelastic medium.

Radially upscaled coefficients on annular zones around the drift

Upscaled HM coefficients over annular shells 0.5 m thick around the drift are shown in Fig. 10. Isotripized Biot coefficient and Biot modulus are presented in Fig. 10a. Near the wall, the Biot coefficient is much higher (B ≈ 0.90) than the global value B = 0.641 given earlier at the global scale of the entire 4 m thick damaged/fractured rock around the drift. It should also be compared to the value B = B M = 0.50 assumed for the intact rock matrix. Mechanical coefficients (Fig. 10b) remain mostly constant in the first 2 m of damage zone near the wall, and then they suffer a dramatic change and quickly reach the values of the intact rock matrix. Finally, some tensorial coefficients (shear compliance S ij ) are presented in Fig. 10c, in order to examine its degree of anisotropy vs. radial distance: clearly the claystone is not only more compliant but also more anisotropic in the near wall region.

Conclusions

The present work present a hydro-mechanical characterization of the fractured porous rock in terms of equivalent continuum, macroscale coefficients, with generally anisotropic geometric structure of the crack system, comprising the following parameters: fabric tensors of the fractured medium (F ijkl , G ijkl ), the equivalent tensorial compliance C ijkl , the expansion coefficient B ¯ij (or reciprocal Biot coefficient) and a "first" Biot coefficient B ij I ( ) . The fluid production term ∆ξ is defined and de- composed into several hydro-mechanical coupling terms. We show that one of the fluid production coupling terms involves a "second" Biot coefficient, B ij II ( ) , which relates fluid production to strain under drained conditions. We also show that another fluid production coupling term involves two distinct Biot moduli, M (I) and M (II) . Both moduli convey the coupling between fluid production and pressure for the upscaled matrix/cracks system, but the first Biot modulus M (I) does not take into account this coupling within the porous matrix itself, while the second Biot modulus M (II) corrects for this and represents more generally the full coupling for the matrix/cracks system.

In order to gain more insights, the theoretical tensorial expressions were then specialized to the case of isotropic crack systems, and the corresponding Hydro-Mechanical coupling coefficients were analyzed in terms of matrix/crack stiffness ratio. These expressions were then used for synthetic upscaling tests indicating that the numerically computed tensorial values of Biot coefficient (B) and modulus (M) are close to the theoretical isotropic values obtained for isotropic cracks. The validation is satisfactory, as far as the upscaled tensorial coefficients are coherent with the corresponding isotropic coefficients. In the isotropic test #1 (3 orthogonal fractures), the upscaled Biot tensor and Biot modulus coincide in value with the ones obtained by the isotropic Fig. 9. Transverse 2D distributions (averaged along the X axis) of upscaled coefficients, from left to right and up to down, N ij , S ij , B ij (0) and 1/M (0) .

Table 8

Global equivalent coefficients: hydro-mechanical coefficients of the fissured + fractured claystone, upscaled at the global scale of the 4 m thick damaged zone and then isotropized by tensor contractions. equations. In the isotropic test #2 (statistically isotropic fractured network), the slight deviation of the isotropic behavior of the Biot tensor and the slight differences in the values compared to the isotropic ones can be explained by the stochastic nature of the fractured network. The tensorial upscaling theory was then applied, numerically, to a field site case in the damaged and fractured claystone (Excavation Damaged Zone) around the 500 m deep GMR gallery at the Meuse/ Haute Marne URL (France). At the global scale of the 4 m thick damaged zone, the Young and Shear moduli are reduced by a factor 4 compared to the intact rock matrix. The global Poisson ratio is slightly higher, near the walls of the gallery compared to the intact rock value. On the other hand, at the local scale, the stiffness moduli can be reduced by one order of magnitude or more in the near wall region. Furthermore, looking at the "3D" coefficients, locally upscaled on a grid of cubic voxels, it is clear that material compliances like the shear compliance matrix S ij , and the inverse Biot modulus 1/M, increase significantly near the gallery wall. The Biot coefficient B ij and the Biot modulus M are calculated first by taking into account the sole effect of the cracks (hypothesis "iii.a", B ij

), and then by taking into account the complete coupling effects due to porous matrix and cracks (hypothesis "iii.b", B ij (I) and B ij (II) ). Looking at the global values shown in Table 8 under the more general hypothesis "iii.b", and comparing these values with those of the intact rock matrix, it is clear that the upscaled Biot coefficient (B) increases with damage, and that the Biot modulus (M) decreases with damage. For all HM coefficients analyzed here, the "farwall" values obtained beyond R = 4.5 m coincide with the intact rock matrix properties. Since the radius of the gallery is R = 2.0 m, it can be concluded that, according to these upscaling calculations, the hydromechanically damaged rock lies in an annular region of thickness 2.5 m around the gallery wall.

From a theoretical point of view, the proposed upscaling method yields generally non-isotropic/non-orthotropic tensorial stress-strainpressure relations, represented by tensorial coefficients of rank up to 4. A detailed description of the theoretical upscaling method based on strain superposition has been provided. Additionally, we have exploited an algebraic theoretical formulation, which has allowed us to define consistently the equivalent isotropic coefficients for these laws. For a 2nd rank tensor like B ij , one way to obtain the isotropic equivalent is to take the spherical part: B = Tr(B)/3 = B kk /3. The upscaled coefficients can be potentially used as inputs for numerical simulations based on "standard" hydro-mechanical codes (simulators), provided that the numerical simulators have the ability to deal with fully tensorial coefficients in the mechanical and hydro-mechanical coupling laws.

In future, we plan to complete the results of this paper by implementing the upscaling method on a variety of other "spatial supports" … For instance, we plan to formulate the upscaled coefficients as continuously distributed parameters using filtering techniques and moving windows (rather than fixed sub-domains or partitions). See for instance the upscaled radial permeability profiles K(r) obtained by. [START_REF] Ababou | Macro-permeability distribution and anisotropy in a 3D fissured and fractured clay rock: 'excavation damaged zone' around a cylindrical drift in callovo-oxfordian argilite (bure). Special issue on "clays in natural & engineered barriers for radioactive waste confinement[END_REF] Similarly, it is expected that the current hydro-mechanical upscaling will allow us to obtain properties such as K (bulk stiffness) and B (spherical Biot coefficient) along "numerical boreholes". Several other extensions of our upscaling method are currently being considered for hydro-mechanical processes to account for: (i) weakly non elastic/non reversible deformation behavior of the discrete cracks, and (ii) retroactive 'feedback' effects of the mechanical variables on the upscaled properties themselves, such as hydraulic conductivity K ij (very sensitive to crack apertures), but also other hydro-mechanical properties. Our current theoretical developments and extensions of the present work indicate that the Biot modulus "M" is most sensitive to crack apertures and to crack normal stiffness; therefore it is expected that "M" should be quite sensitive to mechanical feedback effects due to aperture variations. Finally, we are continuing the extension of our theoretical frame for obtaining equivalent continuum Hydro-Mechanical (HM) laws and coefficients for more general types of dual porosity media like karstic regions with volumetric cavities and holes.

Fig. 1 .

 1 Fig. 1. (a) 2D transverse view of the statistical set of 10 000 small disc shaped fractures; (b) 3D view of the curved surface of a single curved "chevron" fracture; (c) 2D longitudinal section of the 10 000 small fractures (blue) + 46 periodically spaced "chevron" fractures (green) on a 20 m stretch of cylindrical gallery. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

. 4 .

 4 Equivalent compliance, and development of the first Biot constitutive relation 3.4.1. Total strain vs. stress and pressure (superposition of strains)

Fig. 2 .

 2 Fig. 2. Schematic views of an individual crack and its reference frame: a) 2D; b) 3D.

  ⎜

3. 5 .

 5 "Second" Biot equation & fluid production Δξ 3.5.1. Fluid production term Δξ: decomposition by volume variations

  Fluid production term ξ or Δξ represents the specific production of fluid volume inside the control volume, in [m 3 /m 3 ],

Fig. 3 .

 3 Fig. 3. Hydro-mechanical conceptual model: (a) conceptual model of the fractured poro-elastic rock: micro-scale schematic of the solid grains (orange), matrix pores (light blue), and planar cracks (dark blue); (b) control volume of the equivalent poro-elastic continuum formed by the porous matrix and cracks: solid volume (orange), pore volume (cracks+pores)(white), and fluid volume (filling the cracks+pores)(blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

3. 5 . 2 .

 52 Fluid production Δξ in terms of strains and pressure FLUID PRODUCTION TERM Δξ SK DUE TO PORO-ELASTIC SKELE-TON DEFORMATION.

  eqs. (32) and (33) into eq. (31) and multiplying by (1-Φ), we obtain the matrix fluid production in terms of global stress and global pressure:

Fig. 4 .

 4 Fig. 4. Thought experiments: a) fluid compression test (δp > 0) without bulk deformation. b) compression test with incompressible fluid.

3 .

 3 Definition of the 2nd Biot coefficient B ij (II)

4. 1 .

 1 Isotropized expressions and asymptotic behavior of the HM coefficients 4

a

  MPa and K C = 15.6 MPa (using k N * = 10 MPa/m, Oda's hypotheses for B ij : B C = 1; B M = 0 (and = ℓ ℓ f in Oda's expression). b General case for B ij : B C = 1; B M = 0.5. c Oda's hypotheses for M: B C = 1; B M = 0; θ M = 0; K W * = ∞ (and = ℓ ℓ f in Oda's expression). d General case for M: B C = 1; B M = 0.5; θ M = 0.14; K W * = 2.2 GPa.

Fig. 6 .

 6 Fig. 6. 3D and 2D views of the statistically homogeneous and isotropic fractured network used in the "isotropic" homogenization test #2. Generation domain: red lines, (10x10x10) m3; homogenization domain: yellow lines, (8 × 8 × 8) m3; equivalent spherical domain with same volume: dashed blue circle, R ∼ 5 m. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .

 7 Fig. 7. Schematics illustrating several averaging and upscaling options, leading to spatially distributed (3D, 2D transverse, radial or "global") equivalent continuum coefficients.

Fig. 8 .

 8 Fig. 8. Frontal views of the 3D upscaled coefficients, from left to right and up to down, N ij , S ij , B ij (0) (ellipsoïds) and 1/M (0) (spheres) in the first transverse layer (Y,Z) with X ∈ [ 10.0 m, 9.5 m].

Fig. 10 .

 10 Fig. 10. Radial distributions of upscaled tensorial mechanical & hydro-mechanical coefficients of the damaged claystone over annular cylindrical shells: (a) Hydromechanical coupling coefficients (b) Mechanical coefficients (obtained from C ijkl ) (c) Shear compliance tensor S ij (principal components).

Table 1

 1 Upscaling results of the orthogonal fracture set of the isotropic Test #1/Case#1.k S * = 1 MPa/m). The results are presented in Table4.4.3.2. Test#2/case#2: cracks more compliant than matrixIn this case (Test#2/Case#2), the cartesian fractures are more compliant than the matrix: K M =4.17 MPa and K C =1.56 MPa (with k N * = 1 MPa/m, k S * = 0.1 MPa/m). The results are presented in

	Isotropized equations	5.7309	9.6492	6.5461
	Direct upscaling	5.7309	9.6492	6.5461

Table 5 .

 5 4.3.3. Test#2: asymptotic analysis with respect to stiffness ratio r = K C /K M Finally, asymptotic analyses for the upscaled coefficients are summarized in Table6for different values of crack and matrix stiffnesses, or equivalently, with respect to the stiffness ratio r = K C /K M :

	Case	Crack stiff. (MPa)		Upscaled Biot coefficients		Upscaled Biot moduli (MPa)
	r = K C /K M	k N *	k S *	B (0)	B (I)	B (II)	M (0)	M (I)	M (II)
	K C ≫K M 37.4	≫ 10 13	≫ 10 12	0 0.0003	0.5 0.5001	0.57 0.5701	∞ 15645	16.663 16.659	8.6770 8.6760
	3.74	10 10	10 9	0.2104	0.6052	0.6605	19.808	13.768	7.9301
	0.374	10 9	10 8	0.7271	0.8635	0.8826	5.7309	9.6492	6.5461
	0.0374	10 6	10 5	0.9996	0.9998	0.9998	4.1684	8.3342	5.9943
	K C K M	«	«	1	1	1	4.1667	8.3326	5.9936

Note: the term "EDZ" (Excavation Damaged Zone) is a convenient qualitative term that designates, loosely, the perturbed or damaged zone of rock around an excavation (it is usually micro-fissured and fractured).

Acknowledgements

The first two authors wish to acknowledge the support of the Agence Nationale pour la gestion de Déchets Radioactifs (ANDRA), who provided the main scope and the field data to perform the simulations.

Even if we appreciate slight differences between the isotropized values and the upscaled values, and also in the tensorial isotropic behaviour of the Biot coefficients, which may be due to the stochastic nature of the fractured network, they are both in good agreement. The asymptotic behaviour of the results is also in agreement with the expected values. All the results match almost perfectly the values found in the isotropic Test#1.

Site application: equivalent upscaled "HM" coefficients in fractured claystone around a drift

The tensorial upscaling theory is now being applied, numerically, to a field site that was described earlier in Section 2: case of damaged and fractured claystone (Excavation Damaged Zone) around a 500 m deep drift (the "GMR" gallery in the Meuse/Haute Marne Underground Research Laboratory, France).

Input parameters values

Table 7 shows the list of numerical values of the basic geometric, mechanical and hydro-mechanical input parameters used in this case study.

Some of the parameters of Table 7 (particularly for the rock matrix) were selected by comparison with the internal databank of ANDRA, 26 and by comparison with other works in the literature. The mean value of the Biot coefficient B over different layers in the Callovo-Oxfordian (between depths 400-500 m roughly) were calculated by [START_REF] Cosenza | Theoretical prediction of poroelastic properties of argillaceous rocks from in situ specific storage coefficient[END_REF] in their Table 1. They obtained the value B ≈ 0.4 (range: B = 0.37 to 0.42) and observed that this value is close to that obtained directly by hydromechanical tests in other works. They quote in particular the tri-axial tests of [START_REF] Coste | Evolution des propriétés hydrauliques et hydromécaniques lors de l'endommagement des argilites de l'Est. Journées Scientif[END_REF] , who obtained B = 0.36 "for the same claystone", while other authors like [START_REF] Vincké | Etude du comportement poromécanique d'argilites dans les domaines élastiques et post-élastiques[END_REF] obtained higher values from B = 0.4 up to 0.8, also "for the same claystone". Other data obtained in-situ at the Meuse/Haute-Marne URL seem to point to values B ≈ 0.6 for the "intact" rock matrix (although small fissures cannot be avoided on rock matrix samples extracted in-situ). For these reasons, we have taken B M ≈ 0.5 for the intact rock matrix in this study. Concerning the Biot modulus "M ", based on, [START_REF] Cosenza | Theoretical prediction of poroelastic properties of argillaceous rocks from in situ specific storage coefficient[END_REF] their calculated Biot modulus is about M ≈ 8 GPa (range: M = 7.69-8.36 GPa in their Table 1). There are few other more direct in situ measurements of the Biot modulus of this claystone in the literature. We have verified that this value is close to that which can be obtained theoretically (from Biot's isotropic theory of granular materials) for the intact rock matrix using reasonable micro-scale parameters, so we take this computation as the final value: M M = 8.68 GPa.

Post-treatment to obtain upscaled hydro-mechanical coefficients

Spatial averaging procedures for upscaling the HM coefficients

The theoretical upscaling method described in this work can be implemented, technically, for various types of spatial averaging "supports", and for various multidimensional scales. Thus, coefficients can be obtained as 3D fields distributed in (x, y, z), 2D fields distributed in a transverse cross-section (y,z), annular shells at various radial distances from the wall of the gallery, or even as "global" values by taking the annular domain to be as thick as the entire fractured zone (4 m thick). All these distributions of the upscaling results can be obtained whether by direct upscaling (upscaling domain of the corresponding shape) or by sequential upscaling, obtaining first the 3D distribution of the coefficients and then averaging into the other forms, as shown schematically in Fig. 7.

Ellipsoïdal representation of upscaled tensorial coefficients and isotropization

The tensorial version of Darcy's law (involving 2nd rank tensor K ij ) can be used to define a directional conductivity along the gradient direction (K grad ) and a directional conductivity along the flux direction (K flux ); it is then possible to show that the polar plots of K 1/ grad and of K flux both describe an ellipse in 2D or an ellipsoïd in 3D. However, for simplicity, we have chosen here to represent the matrix quantities A ij (that is K ij , B ij , etc.) with ellipsoïds such that the 3 principal axes of the ellipsoïd indicate the directions of the 3 eigenvectors of matrix A ij , and the 3 principal radii of the ellipsoïd are equal to the corresponding 3 eigenvalues of matrix A ij [START_REF] Cañamón | Numerical modeling and upscaling of a 3D fractured rock with Thermo-Hydro-Mechanical coupling[END_REF] . Thus, strongly anisotropic tensors are indicated by strongly elongated ellipsoïds (prolate, like a rugby ball), or strongly flat ellipsoïds (oblate, like a saucepan). The Biot modulus M, or any other scalar quantity, is plotted as a sphere, whose radius indicates the magnitude of the plotted coefficient.

Mechanical compliance (C ijkl ) and stiffness (R ijkl ) are higher order 4th rank tensors, which make it more difficult to represent them in space, and to visualize their spatial distribution. In the present work, we have chosen to reduce the information and to display only some parts of the C ijkl and R ijkl tensors. For example, the stiffness tensor R ijkl can be