J.-F Barthélémy
email: j-francois.barthelemy@ifp.fr

Effective permeability of media with a dense network of long and micro fractures

Keywords: anisotropic effective permeability, self-consistent scheme, fracture permeability, double scale fracture network Nomenclature

This paper presents a new methodology to estimate the effective permeability of random fractured media of any anisotropy containing both micro-fractures and a large number of long fractures cross-cutting the representative volume element. The fractures are replaced by fictitious permeable materials for which the tangential permeability is deduced from a Poiseuille flow. A self-consistent scheme is proposed to derive the macroscopic permeability. On the one hand, the contribution of long fractures to the effective permeability writes by simple superposition of the fracture tangential permeabilities. On the other hand, the contribution of micro-fractures needs to resort to auxiliary problems requiring the computation of second-order Hill (or Eshelby) tensors related to ellipsoids embedded in an anisotropic matrix, for which a complete procedure is detailed. The effect of the micro-fracture normal permeability is put in evidence in the upscaling scheme and analyzed. In particular, it is shown that it must be chosen large enough to allow the connections between families. Examples are finally developed and compared to numerical simulations in the 2D case.

Introduction

A lot of industrial applications, e.g. hydrocarbon production, water supply or subcritical CO 2 storage, are concerned by fluid flows in rocks. The issue of modeling such flows, for example at the scale of a reservoir for the assessment of field productivity and reserves, requires to focus on the porous domain and more specifically on the fracture network. Indeed the latter can be composed of several geological entities of different orders of length magnitude which can all play a very important role in the overall hydraulic behavior of the field. In addition to the pores and fractures at very small scale (decimetric or less) accounting for the intrinsic porosity and permeability of the rock matrix, several fracturation scales are put in evidence in [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models : Implementation and benefits[END_REF]. On the one hand, large scale faults observable by seismic surveys, sub-seismic faults and fracture swarms consisting in clusters of fractures can all be considered as large scale fractures (hectometric or more) taken explicitely into account in flow simulation models. On the other hand, diffuse fractures (metric to decametric) are so numerous that they can not reasonably individually appear in those models; however, their presence may significantly affect the fluid flows. A workflow incorporating the hydraulic effects of fractures at all scales and their possible interactions is detailed in [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models : Implementation and benefits[END_REF] so as to build a complete flow simulation model. An example of model built and calibrated from available static and dynamic data and containing fractures at different scales can be found in [de Galard et al., 2005]. Many real case studies show that, as expected, the presence of large scale fractures of high transmissivity has a strong impact on the field behavior and particularly on well production data, especially when they cross the well ( [START_REF] Cosentino | Integrated study of a fractured middle east reservoir with stratiform super-K intervals -part2: Upscaling and dual media simulation[END_REF], [START_REF] Jenni | History matching of a stochastic model of field-scale fractures: Methodology and case study[END_REF]). Nevertheless, as emphasized in [START_REF] Cosentino | Integrated study of a fractured middle east reservoir with stratiform super-K intervals -part2: Upscaling and dual media simulation[END_REF], the requested accuracy of the simulation used for prediction imposes to build a dual medium model accounting for the interactions between the large draining faults and their neighborhood, namely media of lower permeability but with a high fluid storage capacity as diffuse fracture networks or even the matrix itself. For this purpose, the intrinsic porosity and permeability of the matrix can sometimes be deduced from fracture observations at the core scale (centimetric) as in [START_REF] Howard | Description of natural fracture systems for quantitative use in petroleum geology[END_REF] or from inversion methods based on flowmeter data as in [START_REF] Cosentino | Integrated study of a fractured middle east reservoir with stratiform super-K intervals -part2: Upscaling and dual media simulation[END_REF]. For some reservoirs, dynamic data can be successfully numerically reproduced by a model involving only a permeable matrix and networks of diffuse (metric to decametric) fractures but no large fault [START_REF] Verga | Calibration of fractured reservoirs with dynamic data[END_REF]. The combined influence of diffuse fractures and a permeable matrix suggests resorting to techniques allowing to derive effective properties of dual porosity/permeability media [START_REF] Bourbiaux | A rapid and efficient methodology to convert fractured reservoir images into a dual-porosity model[END_REF]]. In many papers proposing workflows for flow simulation ( [START_REF] Bourbiaux | A rapid and efficient methodology to convert fractured reservoir images into a dual-porosity model[END_REF]], [START_REF] Sabathier | A new approach of fractured reservoirs[END_REF]], [START_REF] Cacas | Nested geological modelling of naturally fractured reservoirs[END_REF], [START_REF] Garcia | Fast and efficient modeling and conditioning of naturally fractured reservoir models using static and dynamic data[END_REF]) or presenting real cases [START_REF] Cosentino | Integrated study of a fractured middle east reservoir with stratiform super-K intervals -part2: Upscaling and dual media simulation[END_REF], [de Galard et al., 2005]), an important step consists in representing the hydraulic effect of diffuse fractures by an effective permeability. In some cases, the proposed upscaling method relies on a numerical resolution on an explicit local discrete fracture network ( [START_REF] Sabathier | A new approach of fractured reservoirs[END_REF]], [Pouya, 2005]). Besides, the issue of upscaling hydraulic properties of fracture networks by analytical methods has also been raised by several authors: [Oda, 1986] presents a method based on a linear superposition of the contributions of different fracture families, [Zimmerman, 1996] implements the differential scheme and Maxwell's method for 2D isotropic fracture networks, [Fokker, 2001] and [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF]] take advantage of the selfconsistent scheme for, respectively, 2D anisotropic networks and 3D isotropic networks. This paper deals with the estimate of the permeability of a fractured medium with a 3D anisotropic fracture network. As the existence of a r.v.e. (representative volume element) is assumed in the sequel, it becomes natural to define the macroscopic scale as that considering the r.v.e. as an elementary particle with effective properties and the microscopic scale as that refering to the different phases composing the r.v.e. (fracture families and solid phase). The diffuse fracture network considered here contains a high density of families of micro-fractures and long fractures with small spacing. The adjective "long" characterizes fractures cross-cutting the r.v.e.. The fracture centers are randomly distributed and the orientations of the families are arbitrary, leading then to a possible anisotropic effective permeability. The notion of r.v.e. implies the separation scale between its characteristic length and those of the fracture families, namely the size of micro-fractures or the average spacing between long fractures cross-cutting the r.v.e.. It should be recalled that the prefix "micro" only refers to a small size compared to that of the r.v.e., whatever the order of magnitude of the latter (decimetric to hectometric according to real case studies [START_REF] Sabathier | A new approach of fractured reservoirs[END_REF]). The framework of random media upscaling, which is also commonly used to estimate mechanical effective properties [Zaoui, 2002], is considered here so as to provide a methodology based on the self-consistent scheme to estimate the effective permeability K hom in the general case of a 3D anisotropic distribution of fractures. The contribution of micro-fractures, requiring the determination of anisotropic Hill polarization tensors, extends to 3D the results obtained in 2D in [Fokker, 2001]. The self-consistent scheme allows to highlight percolation thresholds in the case of a strong contrast between the phases (fracture space and matrix). Considering fractures as fictitious Darcean materials, i.e. characterized by a permeability tensor, the paper focuses on the relevance of the different components of this tensor. In particular, the influence of the normal component on the macroscopic permeability is analyzed in details for both fracture types and illustrated by (3D or 2D) isotropic distributions of fractures allowing to derive closed-form solutions of K hom . In the general case of an anisotropic distribution, numerical solutions of the self-consistent scheme can be obtained and compared to calculations made on explicit discrete fracture networks.

1 Description of the phases

Description of the fracture network

In natural diffuse fractured media as reservoirs, the issue of the existence of a r.v.e. on which upscaling techniques can provide macroscopic properties such as permeability and percolation thresholds has sometimes to be reconsidered. This is for example the case for fracture networks with a fractal correlation ( [START_REF] Bonnet | Scaling of fracture systems in geological media[END_REF], [START_REF] Darcel | Connectivity properties of two-dimensional fracture networks with stochastic fractal correlation[END_REF], [Darcel, 2002]). Nevertheless, aiming at applying upscaling methods, the media which are studied in the present paper are based on the assumption of a random fracture distribution. In this case, some papers ([de Dreuzy et al., 2001a], [START_REF] De Dreuzy | Hydraulic properties of twodimensional random fracture networks following a power law length distribution: 2-Permeability[END_REF]) have shown the influence of the geometrical parameters of the fractures, in particular the length distribution, on the connectivity and the permeability of the medium. A domain of validity of the effective medium theory based on the self-consistent scheme that will be developed hereafter is provided in [de Dreuzy et al., 2001a]. According to the latter, this technique can only be applied for a connectivity beyond the percolation threshold whereas other theories or numerical devices depending on the length distribution should be applied before and at the threshold. Nevertheless, the latter as well as the permeability tensor can be rather well estimated by the self-consistent scheme ( [Fokker, 2001], [START_REF] Pozdniakov | A self-consistent approach for calculating the effective hydraulic conductivity of a bimodal, heterogeneous medium[END_REF]) except in the very neighborhood of the threshold [Sahimi, 1995] which appears as a transition at which the permeability progressively changes of order of magnitude. It is worth mentioning that a percolation threshold can also be obtained by means of Maxwell's method for a 2D isotropic network in [Zimmerman, 1996]. Although thresholds derived by such effective media theories do not exactly coincide with numerical simulations, those theories can anyway provide rough but quick estimates which can be refined by numerical computations on effective fracture networks or which can also be used to assess the sensitivity to the fracture network characteristics. The present study relies on the assumption of existence of a r.v.e. with N fracture families. The i th family (1 ≤ i ≤ N ) contains a large number of fractures all characterized by the same length, aspect ratio, orientation and hydraulic properties. It must be emphasized here that the following developments are based on a finite number of families. This means that the characteristics (length, shape, orientation or conductivity) must take discrete values. Nevertheless, continuous distributions of the latter (defined for example by means of probability density functions) could be accounted for by building an equivalent discrete set of families. For example, in the framework of mechanics, [Zhu, 2006] presents methods either based on numerical integration by gaussian quadrature or on the definition of a representative set of families allowing to replace a continuous orientation distribution by a discrete one. The set of N fracture families is partitioned in two subsets: families of long fractures cross-cutting the r.v.e. (i ∈ C) and families of micro-fractures (i ∈ M).

• Families of fractures cross-cutting the r.v.e. (i ∈ C)

A fracture cross-cutting the r.v.e. is geometrically described here as a plane, assuming then that its curvature can be neglected at the scale of the r.v.e.. The average spacing between fractures of the i th family is denoted by d i . This distance can be considered as the characteristic dimension of the family and must be very small compared to the size of the fractures and to that of the r.v.e. (Fig. 1). The orientation is entirely defined by a normal vector i.e. two Euler angles. Adopting a 3D description of the fracture, the latter can be seen as a volume delimited by two parallel planes. Introducing the average geometrical aperture h i , possibly different from the hydraulical aperture defined later, the infinitesimal volume fraction of the family is:

f i = h i d i (1) L d i
Figure 1: Fractures cross-cutting the r.v.e.

• Families of micro-fractures (i ∈ M) Considering fractures as 2D plane objects, the micro-fractures are modeled as ellipses and the characteristic length of each family (largest radius a i ) is assumed to be very small compared to the characteristic dimension of the r.v.e.. The aspect ratio characterizing the shape in the fracture plane is denoted by η i = b i /a i where b i is the small radius. If N i denotes the number of fractures per volume unit and ǫ i = N i a 3 i the classical density introduced by Budiansky [Budiansky and O'Connell, 1976], the average spacing between two elements of a family is then the inverse of the average number of fractures intersecting a segment of unit length orthogonal to the fracture plane:

d i = 1 N i π a 2 i η i = a i π ǫ i η i (2)
The orientation of all the micro-fractures of a given family is defined by a 3D frame (three orthonormal vectors) i.e. three Euler angles as described later in Fig. 7. A network composed of three families of micro-fractures is shown on Fig. 2. For sake of convenience in the self-consistent scheme based on Eshelby's problem and developed in the following, we also introduce a 3D model of micro-fractures considered as ellipsoids of largest radii a i and b i . The third dimension is then defined by the small radius c i or the ratio ω i = c i /a i ≪ 1 (see Fig. 3). Hence, the volume fraction of the micro-fracture family writes:

f i = 4 3 π ǫ i η i ω i (3)
Figure 2: Network of micro-fractures

a i b i = η i a i c i = ω i a i (ω i ≪ 1) Figure 3: Ellipsoidal fracture
It is assumed that all the characteristic lengths of the families (spacing in the case of fractures cross-cutting the r.v.e. and size in the case of micro-fractures) are of the same order of magnitude, i.e. infinitesimal with respect to the dimension of the r.v.e., in order to incorporate all the families to the same upscaling process. Otherwise several successive steps have to be performed.

Hydraulic properties of the phases

Ene and Sanchez-Palencia have shown, in the framework of periodic homogenization, that the steady-state, laminar, single-phase and incompressible flow of a newtonian viscous fluid in a porous medium obeys to the Darcy law at the upper scale, which means that the filtration velocity is linearly related to the pressure gradient [START_REF] Ene | Équations et phénomènes de surface pour l'écoulement dans un modèle de milieu poreux[END_REF]:

Q = - K hom µ • ∇P (4)
where the second-order tensor K hom has the dimension of a surface (m 2 ) and depends on the sole morphology of the porous space (i.e. not on the fluid viscosity µ). The aim of the upscaling procedure is to determine or at least estimate K hom . It is commonly assumed that the relationship (4) also applies to random media but estimating K hom from Stokes equations at the lower scale proves to be a hard task. That is why some authors ( [Fokker, 2001], [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF], [START_REF] Pozdniakov | A self-consistent approach for calculating the effective hydraulic conductivity of a bimodal, heterogeneous medium[END_REF]) have proposed to replace the different constituents of the microscopic scale (fractures and solid phase) by fictitious Darcy media. The advantage of such a method is that the same type of equation, namely the Darcy law, applies at both the microscopic and macroscopic scale.

The permeability of the solid phase is assumed to be isotropic K s = K s 1. This permeability results from the presence of micropores at a lower scale and is possibly obtained by a previous upscaling procedure. The permeability of the fractures writes with tangential and normal permeabilities which are a priori different:

K i = K i t + K i n n i ⊗ n i with K i t = K i t 1 ni = K i t (1 -n i ⊗ n i ) (5)
Following a classical result of rock hydraulics [START_REF] Guéguen | Introduction to Physics of Rocks[END_REF], K i t accounts for the overall tangential permeability of the fracture and is identified with the ratio between the average velocity of a Poiseuille flow between two idealized planes and the pressure gradient in the direction of the planes. The space e i between those fictitious planes is called the hydraulic aperture. It may differ from its geometrical counterpart because the latter may not be uniform along the fracture due to the face rugosity or simply due to the specific 3D geometry, for instance in the case of an ellipsoid. The equivalent Poiseuille flow then gives:

K i t = e 2 i 12 (6)
It is assumed in the following that the tangential permeability K i t is much larger than K s implying the impermeability of the fracture faces. As the equivalent Poiseuille problem is based on a flow parallel to the fracture plane, it does not give any indication on how to estimate K i n . Recalling that a permeability is the macroscopic evidence of charge loss due to friction at the fluid-solid interface at the microscopic scale, the intuition could lead us to choose a value of K i n much larger than K i t because the latter should be much more impacted by the narrowness of the fracture aperture (6) than the former. Nevertheless, some authors have proposed different other choices: for example, the micro-fracture permeability is taken isotropic [Fokker, 2001] or [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF]. Thus, replacing the real fluid flow in the fracture space by an effective flow in a fictitious Darcy medium somehow makes K i n in (5) a mathematical parameter which does not seem so obvious to choose. It is then necessary to raise the issue of the relevance of any choice of K i n by analyzing its impact on the overall permeability, which is the purpose of section 3.

(K i n = K i t ) in
K i n = K s ≪ K i t in

Determination of the macroscopic permeability

The upscaling problem to solve on the r.v.e. Ω is similar to that allowing to exhibit effective thermal or electrical properties i.e. is the following elliptic problem:

                         div q = 0 (Ω) (a) q = - K s µ grad p (Ω s ) (b) q = - K i µ • grad p (Ω i ) (c) p = ∇P • x (∂Ω) (d) (7)
where Ω s is the space occupied by the matrix, Ω i the space occupied by the i th fracture family and ∇P is a given uniform vector identified as the macroscopic pressure gradient. The linearity of the equations at stake implies the existence of a concentration tensor field such that the solution of ( 7) writes:

∀x ∈ Ω, grad p(x) = A(x) • ∇P (8)
The macroscopic permeability is then given by:

Q =< q > Ω = - K hom µ • ∇P with K hom =< K • A > Ω = 1 - N i=1 f i K s • < A > Ω s + N i=1 f i K i • < A > Ω i (9)
where < • > ω denotes the average operator over the domain ω. The upscaling methods aim at providing estimates of the averages of the concentration tensor or equivalently the relationships between the averages of q in each phase and ∇P .

2.1 Contribution of the fractures cross-cutting the r.v.e.

As they cross-cut the r.v.e., those fractures are subjected to a pressure difference directly related to the boundary condition (7 d). In fact, as pointed out in [Oda, 1986], the fluid flow in those fractures sees an effective pressure gradient obtained by projection of the macroscopic pressure gradient onto the fracture plane. By construction of K i t and recalling that the fracture faces are impermeable, the flow in such a family is purely tangential and writes:

< q > Ω i = - K i t µ 1 n i • ∇P = - K i t µ • ∇P (∀i ∈ C) (10)
The contribution of this family to K hom (9) can then be identified in (10) by:

K i • < A > Ω i = K i t = K i t 1 ni (∀i ∈ C) (11)
which appears as a superposed term in ( 9) without any interaction with other fracture families. K i n is obviously not involved in the explicit contribution of the corresponding family to K hom (11). Focusing on the estimate of the pressure gradient average i.e. < A > Ω i , it can first be noted that (11) leads to

1 ni • < A > Ω i = 1 ni and K i n n i ⊗ n i • < A > Ω i ≪ K i t (∀i ∈ C) (12)
The first relationship of ( 12) allows to determine without ambiguity the projection < A > Ω i onto the fracture plane whereas the projection onto the normal depends on the choice of K i n . If K i n is chosen at least of the same order of magnitude as K i t , then (12) implies that

< A > Ω i = 1 n i (∀i ∈ C) (13)
The choice K i n ≪ K i t still satisfies the consistency with (11) but leads to three undetermined components of < A > Ω i . In fact, the following developments in sections 2.2 and 2.4 will show that those components of < A > Ω i do not actually play a role in the upscaling process provided that they keep bounded values. ( 13) can therefore still be used.

Contribution of the micro-fractures

Unlike the previous case of fractures cross-cutting the r.v.e., the interactions of smaller fractures with their neighborhood must be strong. Their influence on the overall permeability comes from the connections between them, which means that each phase of the heterogeneous medium is in contact with elements of the other phases. According to [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF], the self-consistent scheme is adapted for its capacity to account for a disordered microstructure and to exhibit percolation thresholds. In the framework of the self-consistent scheme, the average of the concentration tensor in each phase is estimated by means of an auxiliary problem of an ellipsoid I embedded in the infinite effective medium itself and subjected to a remote pressure gradient ∇P o . The use of ∇P o different from ∇P corresponds to a classical way to take into account interactions in upscaling schemes (see for example [Zaoui, 2002] in the mechanical framework). The relationship between ∇P o and ∇P is obtained by means of the consistency rule:

< grad p > Ω = ∇P (14)
The auxiliary problem, known as the Eshelby problem [Eshelby, 1957], shows an homogeneous pressure gradient within the ellipsoid, which is used as an estimate of the pressure gradient average in the fractures of the i th family in the r.v.e. problem:

< grad p > Ω i = Λ i • ∇P o with Λ i = 1 + P i hom • K i -K hom -1 (∀i ∈ M) (15)
where the second-order Hill tensor P i hom depends on K hom and on the shape and orientation of the ellipsoid. The ellipsoid representing the i th family is defined thanks to a second-order tensor H i by the equation:

x • t H i • H i -1 • x ≤ 1 (16)
For any ellipsoid, P i hom results from K hom in the general 3D anisotropic case according to the following procedure (see Appendix. for a detailed proof):

1. Diagonalize K hom and introduce the tensor R such that

t Q • Q = 1 and t Q • K hom • Q = 3 j=1 K hom j e j ⊗ e j (17) R = 3 j=1 K hom j Q • e j ⊗ e j (18)
2. Compute the tensor

H ′ i = H i • t R -1 and diagonalize t H ′ i • H ′ i t S i • S i = 1 and t S i • t H ′ i • H ′ i • S i = A 2 i e 1 ⊗ e 1 + B 2 i e 2 ⊗ e 2 + C 2 i e 3 ⊗ e 3 (19) with A i ≥ B i ≥ C i * 3. Compute the tensor ∆ i = I i A e 1 ⊗ e 1 + I i B e 2 ⊗ e 2 + I i C e 3 ⊗ e 3 (20)
with

• if A i > B i > C i I i A = A i B i C i (A 2 i -B 2 i ) A 2 i -C 2 i (F i -E i ) (21)
I i C = A i B i C i (B 2 i -C 2 i ) A 2 i -C 2 i B i A 2 i -C 2 i A i C i -E i (22)
I i B = 1 -I i A -I i C ( 23
)
where

F i = F (θ i , k i
) and E i = E(θ i , k i ) are respectively the elliptic integrals of the first and second kinds of amplitude and parameter:

θ i = arcsin 1 - C 2 i A 2 i ; k i = A 2 i -B 2 i A 2 i -C 2 i (24) • if A i > B i = C i I i B = I i C = 1 2 A i A i A 2 i -C 2 i -C 2 i argch (A i /C i ) (A 2 i -C 2 i ) 3/2
(25)

I i A = 1 -2 I i B ( 26
)
• if A i = B i > C i I i A = I i B = 1 2 C i A 2 i arccos (C i /A i ) -C i A 2 i -C 2 i (A 2 i -C 2 i ) 3/2 (27)
I i C = 1 -2 I i A (28) • if A i = B i = C i I i A = I i B = I i C = 1 3 (
29) * A permutation between the eigenvectors is possibly used so as to respect this order.

Finally compute

P i hom = t R -1 • S i • ∆ i • t S i • R -1 (30)
The determination of the average of the concentration tensor < A > Ω i will be achieved when the relationship between ∇P o and ∇P is found. However, as this relationship comes from the consistency rule ( 14), the average of the pressure gradient within the matrix must also be estimated through an auxiliary problem even if the explicit contribution of the matrix in ( 9) could be neglected because of the low permeability of the matrix.

Contribution of the matrix

In the self-consistent scheme, the matrix phase also requires the resolution of an auxiliary problem to estimate the average of the pressure gradient. Nevertheless, unlike polycrystals in which each phase plays an equivalent role, the matrix here plays a specific role from the geometrical point of view. The latter, which occupies almost all the r.v.e. space, is also represented by an ellipsoid in its auxiliary problem but the shape of the ellipsoid may not be as obvious as for the polycrystals. This ellipsoid is supposed to account for the spatial distribution of the matrix with respect to the fractures crossing it. More precisely, it seems judicious to quantify the anisotropy of the inter-fracture space by means of the family orientations n i and spacings d i as shown on Fig. 4 and build the H s tensor related to the matrix ellipsoid according to the definition:

H s = N i=1 1 d i n i ⊗ n i -1 (31)
Knowing the shape of the ellipsoid, the average of the pressure gradient within the matrix can then be estimated

< grad p > Ω s = Λ s • ∇P o with Λ s = 1 + P s hom • K s -K hom -1 (32)
where P s hom depends on K hom and H s (31).

Macroscopic permeability

Using the estimates ( 15) and (32) in the consistency rule ( 14) provides the relationship between the auxiliary pressure gradient ∇P o and the actual one ∇P :

∇P o = 1 - N i=1 f i Λ s + i∈M f i Λ i -1 • 1 - i∈C f i < A > Ω i • ∇P (33)
Since the fracture volume fractions are infinitesimal, the average pressure gradients i.e. < A > Ω i of the families of long fractures play a role in (33) only if it presents singular values. If < A > Ω i is bounded as in ( 13), (33) gives:

∇P o = Λ s + i∈M f i Λ i -1 • ∇P (34)
A singular component of < A > Ω i , along n i ⊗ n i to remain consistent with (12), would have implied a reduction of the level of ∇P o and consequently of K hom , which would have made the corresponding long fracture family play the unexpected role of a barrier against the fluid flow possibly coming from micro-fracture families. Finally, the average of the concentration tensor in the micro-fracture families (resp. in the matrix) can be written with the help of (15) (resp. ( 32)) and (34):

< A > Ω i = Λ i • Λ s + i∈M f i Λ i -1 (∀i ∈ M) (35)
< A > Ω s = Λ s • Λ s + i∈M f i Λ i -1 (36)
The macroscopic permeability K hom can now be estimated by taking advantage of (10), ( 35) and ( 36) in (9):

K hom = 1 - N i=1 f i K s • Λ s + i∈M f i K i • Λ i • Λ s + i∈M f i Λ i -1 + i∈C f i K i t ( 37
)
It is recalled that, in the framework of the self-consistent scheme, the tensors Λ i (15) and Λ s (32) nonlinearly depend on K hom , which makes (37) an implicit equation characterizing K hom :

K hom = G K hom (38)
Although closed-form solutions can be obtained in simple cases (see section 4), it is necessary to use a numerical method to find a solution of (38). A recursive fixed-point method will be implemented in the following examples:

K 0 = K s and K j = G (K j-1 ) until K j -K j-1 < ε K j-1 (39)
However, the relevance of the upscaling scheme still relies on the choice of the normal component of the microfracture permeability, which is the purpose of the next paragraph.

Influence of the normal component of the fracture permeability

It has been shown earlier that the macroscopic permeability does not depend on K i n of the long fractures, which contributions are obtained by simple superpositions. The case of micro-fractures is different since they are confined inside the r.v.e. and the flow and the pressure gradient established in them strongly depend on their neighborhood. In the self-consistent scheme, the solution in the micro-fractures has been estimated by that of an ellipsoid surrounded by the effective material. Moreover, according to Voigt bound, K hom as a quadratic form is always lower than < K > Ω and is reasonably of the same order of magnitude of the latter when the fractures are connected or in the presence of long fractures cross-cutting the r.v.e. (37). This implies that K hom is necessarily much lower than the tangential permeability K i t of any fracture contributing to the overall permeability since the fracture volume fraction is infinitesimal. Therefore, the fluid flow is almost parallel to the fracture plane in the auxiliary problem without any condition on K i n as the neighboring material is much less permeable than the ellipsoidal fracture. The quantitative influence of K i n on K hom in the case of micro-fractures has then to be examined more in details. To this end, the behavior of Λ i (15) and P i hom for small values of the aspect ratio ω i is first studied. The Taylor expansion of P i hom with respect to ω i can be obtained from another expression of the Hill tensor than the one using the Green function exploited in section 2.2 [START_REF] Dormieux | Microporomechanics[END_REF]:

P i hom = det H i 4 π ξ =1 ξ ⊗ ξ (ξ • K hom • ξ) ξ • t H i • H i • ξ 3/2 dS ξ (40)
Calling (l i , m i , n i ) the eigenvectors of the ellipsoid (see Fig. 3), H i writes:

H i = a i (l i ⊗ l i + η i m i ⊗ m i + ω i n i ⊗ n i ) (41)
Using in (40) the parametrization ξ = √ 1z 2 (cos φ l i + sin φ m i ) + z n i and dS ξ = dz dφ as well as an analogy with a result obtained in the framework of elasticity in [START_REF] Huang | On a flat ellipsoidal inclusion or crack in three-dimensional anisotropic media[END_REF]], P i hom writes:

P i hom = n i ⊗ n i n i • K hom • n i + ω i Π i hom + O ω 2 i K hom (42)
with:

Π i hom = - η i 4 π 1 z=-1 2 π φ=0 z √ 1 -z 2 1 cos 2 φ + η 2 i sin 2 φ 3/2 ∂ ∂z ξ ⊗ ξ ξ • K hom • ξ dφ dz (43)
It is worth noting on (43) that Π i hom = O 1/ K hom . In the case of an oblate spheroid (η i = 1) in an isotropic effective medium K hom = K hom 1, (42) becomes:

P i hom = n i ⊗ n i K hom + ω i Π i hom + O ω 2 i K hom with Π i hom = π 4 K hom (1 n i -2 n i ⊗ n i ) (44)
Substituting ( 42) in ( 15), it comes out that:

Λ i -1 = 1 - n i ⊗ K hom • n i K hom nini + K i n n i ⊗ n i K hom nini + ω i Π i hom • K i -K hom + O ω 2 i K i K hom (45)
where

K hom nini = n i • K hom • n i .
As already previously mentioned, it seems reasonable to assume that, beyond the percolation threshold, K hom is of the same order of magnitude as f i K i t if the i th family significantly contributes to K hom . Keeping the terms of highest order, Λ i -1 can be written by blocks:

Λ i -1 = 1 n i + ω i 1 n i • Π i hom • K i t a i +ω i 1 n i • -Π i hom • K hom • n i ⊗ n i b i + n i ⊗ n i • ω i K i t Π i hom - K hom K hom nini • 1 n i c i + K i n K hom nini -ω i n i • Π i hom • K hom • n i δ i n i ⊗ n i (46)
The tensors a i , b i and c i introduced in ( 46) do not explicitely depend on K i n and are all of the same order of magnitude, namely O (1). As an endomorphism of R 3 , a i is singular whereas, as an endomorphism of the fracture plane, a i can be inverted using the following notation:

a i -1 = 1 n i • 1 + ω i 1 n i • Π i hom • K i t -1 = 1 n i • 1 + ω i K i t 1 n i • Π i hom • 1 n i -1 (47)
In addition, (46) shows that K i n explicitely appears only in δ i (bottom right term i.e. along n i ⊗ n i ), which order of magnitude controls that of Λ i . As regards the order of magnitude of δ i and its dependence on K i n , two cases can be examined:

• K i n ≫ ω i K hom ⇒ δ i = K i n /K hom nini
For this condition, Λ i writes:

Λ i = a i -1 + K hom nini K i n -c i • a i -1 + n i ⊗ n i (48)
On the one hand, recalling that f i is given by ( 3), it comes that f i Λ i can be neglected in the interaction term (34). On the other hand, ( 48) also allows to compute the term K i • Λ i which will be crucial in (37). Indeed, as 1 n i • c i = 0 (46), it comes:

K i • Λ i = K i t • a i -1 + O K hom = K i t • 1 + ω i 1 n i • Π i hom • K i t -1 + O K hom (49)
and, in the case of an oblate spheroid in an isotropic effective medium:

K i • Λ i = K i t 1 + π 4 ωi K i t K hom + O K hom (50)
Taking ( 49) and ( 32) into account as well as the impermeability of the matrix, the implicit expression of K hom (37) can be simplified:

K hom = i∈M f i K i t • 1 + ω i 1 n i • Π i hom • K i t -1 • 1 -P s hom • K hom + i∈C f i K i t (51)
• K i n = O ω i K hom In this case, δ i can be written:

δ i = ω i δ ′ i with δ ′ i = K i n ω i K hom nini -n i • Π i hom • K hom • n i = O (1) (52)
To the first non-zero order, Λ i writes:

Λ i = 1 ni - c i ω i δ ′ i • a i -1 • 1 - b i • c i • a i -1 δ ′ i -1 • 1 ni - b i δ ′ i + n i ⊗ n i ω i δ ′ i (53)
Unlike the previous case, f i Λ i can not be neglected when ω i takes small values. Using (3), f i Λ i writes to the first order:

f i Λ i = 4 3 π ǫ i η i δ ′ i -c i • a i -1 • 1 - b i • c i • a i -1 δ ′ i -1 • 1 ni - b i δ ′ i + n i ⊗ n i (54)
leading, in the case of an oblate spheroid in an isotropic effective medium for which P i hom is given by (44), to:

f i Λ i = 4 3 π ǫ i n i ⊗ n i K i n ωi K hom + π 2 ( 55
)
K i • Λ i can also be obtained:

K i • Λ i = K i t • a i -1 • 1 - b i • c i • a i -1 δ ′ i -1 • 1 n i - b i δ ′ i + O K hom (56)
leading to (50) in the case of an oblate spheroid in an isotropic effective medium. The non-zero term (54) involved in ( 34) and (37) highlights the negative influence of a too small value of K i n on K hom . Indeed, this choice induces an unexpected role of barrier of the fracture in its transverse direction n i contributing to decrease at least the value of K hom nini . Nevertheless, keeping such a small value of K i n could be used to somehow account for cementing on the fracture faces impeding the fluid flow from other fractures.

From the results obtained hereabove, one can now conclude that, in the modeling consisting in replacing a fracture by a fictitious Darcy medium, the most relevant choice of the transverse permeability is specific to the fracture type. Whereas the choice of K i n of the long fractures has no consequence on the macroscopic permeability, in the case of a micro-fracture, K i n must be large enough (K i n ≫ ω 2 i K i t ) to avoid creating an artificial obstacle in the transverse direction except if it corresponds to the aim of the model.

Implementation of the self-consistent scheme on examples

4.1 Isotropic distributions of micro-fractures

• 3D case
In this paragraph, the fracture orientation distribution is assumed to be isotropic and all the fractures are modeled by identical oblate spheroids of aspect ratio ω with the same permeability. The macroscopic permeability is isotropic i.e. K hom = K hom 1 and the ellipsoid related to the matrix is a sphere according to (31). It is worth noting that, given an overall fracture density ǫ, the same permeability tensor K hom can be obtained with a network of three orthogonal families of density ǫ/3 each. Recalling that P i hom is written in (44) for oblate spheroids and P s hom = (3K hom ) -1 1 for a sphere, the macroscopic permeability is obtained by solving (37), which leads to a relationship of the type:

K hom ω K t = L K s K t , K n K t , ǫ, ω (57)
Fig. 5 shows numerical computations of K hom /K t with respect to the density ǫ for ω = 10 -3 and different values of K s /K t and K n /K t , thus illustrating the role of the order of magnitude of the latter which has been highlighted in section 3. When K s tends towards 0, closed-form solutions of K hom can be found.

-If K n ≫ ω 2 K t , a percolation threshold is detected for ǫ = 27/64 and, as already found in [Fokker, 2001], the permeability beyond the threshold writes:

K hom ω K t = π 108 (64 ǫ -27) + O (ω) , ǫ ≥ 27 64 (58)
It is recalled that this percolation threshold ǫ = 27/64 ≈ 0.42 has been compared in [Fokker, 2001] to an estimate ǫ ≈ 0.30 (i.e. a fracture volume fraction f ≈ 1.27ω) coming from numerical computations on networks of randomly placed ellipsoids provided by [START_REF] Garboczi | Geometrical percolation threshold of overlapping ellipsoids[END_REF], showing that the selfconsistent scheme overestimates the threshold in 3D.

-If K n ≪ ω 2 K t , the percolation threshold is reached for a larger value namely ǫ = 9/16 highlighting the negative effect of a too small value of K n compared to the previous case. As already found in [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF], the permeability beyond the threshold is inferior to that of (58) and writes:

K hom ω K t = 3 π 4 16 ǫ -9 16 ǫ + 27 + O (ω) , ǫ ≥ 9 16 (59)
The irrelevance of the model K n ≪ ω 2 K t can also be deduced from the expression (59). Indeed, let the number of fractures per volume unit N and the aperture c be given and let the radius a increase so as to examine its influence on K hom . Recalling that ǫ = N a 3 and K t = c 2 /3, it comes from (59) that K hom varies as 1/a when a is large, meaning that the fracture length would have a negative 13 0 0.25 0.50 0.75 1.00 1.25 1.50

0 0.5 1.0 1.5 2.0 2.5 K hom ω Kt ǫ K s /Kt = 3.10 -4 ; Kn/Kt ≥ 10 -3 K s /Kt = 10 -4 Kn/Kt ≥ 10 -3 K s /Kt = 10 -5 ; Kn/Kt ≥ 10 -3 K s /Kt = 10 -5 ; Kn/Kt = 10 -6 K s /Kt = 10 -5 ; Kn/Kt ≤ 10 -8 K s /Kt = 0 ; Kn/Kt ≥ 10 -3 K s /Kt = 0 ; Kn/Kt ≤ 10 -9
Figure 5: 3D isotropic distribution of spheroidal fractures (ω = 10 -3 ) effect on the permeability, which is unrealistic. This, once again, shows that a small value of K n makes the fractures play the unexpected role of obstacles in their transverse direction.

• 2D case

The same computation can be done on a 2D fractured medium. The geometry is invariant along e 3 and isotropic in the (e 1 , e 2 ) plane i.e. K hom is transversely isotropic:

K hom = K hom // (e 1 ⊗ e 1 + e 2 ⊗ e 2 ) + K hom ⊥ e 3 ⊗ e 3 (60)
The fractures are modeled as flat cylinders of elliptic section. As regards the permeability, one can equivalently consider a random distribution of orientation or two orthogonal orientations. Following the procedure described in section 2.2 in the limit case of a largest axis tending towards infinity, the Hill tensor related to a cylinder of elliptic section (same aspect ratio ω for all the families) embedded in a transversely isotropic medium with the same axis is:

P i hom = 1 K hom // ω 1 + ω m i ⊗ m i + 1 1 + ω n i ⊗ n i (61)
where n i denotes the normal vector i.e. the small axis direction and m i the direction orthogonal to n i in the (e 1 , e 2 ) plane i.e. the large axis direction. It is worth noting that the expression (61) is also found in [Fokker, 2001] in the framework of a complete 2D resolution based on the 2D Green function. The isotropy in the (e 1 , e 2 ) plane leads to the following expression of P s hom

P s hom = 1 2 K hom // (e 1 ⊗ e 1 + e 2 ⊗ e 2 ) (62)
If f denotes the total fracture volume fraction, the fracture density ǫ in the 2D framework is such that f = π ǫ ω. By projection of (37) on e 3 , K hom ⊥ is found, as expected, as an arithmetic average:

K hom ⊥ = (1 -π ǫ ω) K s + π ǫ ω K t (63) 0 0.5 1.0 1.5 2.0 2.5 3.0 0 0.5 1.0 1.5 K hom / / ω Kt ǫ K s /K t = 10 -4 ; Kn /K t ≥ 10 -3 K s /K t = 10 -5 Kn /K t ≥ 10 -3 K s /K t = 10 -6 ; Kn /K t = 10 -5 K s /K t = 10 -6 ; Kn /K t = 10 -6 K s /K t = 0 ; Kn/K t = 5.10 -8 K s /K t = 0 ; Kn/K t ≥ 10 -3 Figure 6: 2D isotropic distribution of elliptic fractures (ω = 10 -3 )
where the K s term disappears for an impermeable matrix. The solution in the (e 1 , e 2 ) plane is shown on Fig. 6 for ω = 10 -3 for different values of K s /K t and K n /K t . Moreover, a closed-form solution is available in the case of an impermeable matrix and K n ≫ ω 2 K t :

K hom // ω K t = π 4 ǫ -1 + O (ω) , ǫ ≥ 4 π (64)
This percolation threshold ǫ = 4/π ≈ 1.27 can be compared to the value ǫ ≈ 1.4 obtained from numerical computations by several authors for the same problem of draining fractures in an impermeable matrix ( [Robinson, 1983], [Berkowitz, 1995], [START_REF] Bour | Connectivity of random fault networks following a power law fault length distribution[END_REF]). Unlike the 3D case, the self-consistent scheme underestimates the percolation threshold but gives a better estimate. It can be mentioned that, in the inverted problem of impermeable needle-like slits in a permeable matrix, the self-consistent scheme as well as Maxwell's method give the same percolation threshold ǫ = 4/π (above which K hom // = 0) and even the same function K hom // (ǫ) which is rather successfully compared to numerical simulations except in the very neighborhood of the threshold [Zimmerman, 1996]. Whereas a percolation threshold is found here for a large enough value of K n , the case K n ≪ ω 2 K t does not allow the macroscopic permeability to reach the order of magnitude of ω K t i.e. does not exhibit any percolation threshold for a finite value of ǫ (see Fig. 6), once more illustrating the effect of K n .

Anisotropic distributions of fractures • 3D case

The algorithm leading to the solution to ( 37) is now implemented in the case of an anisotropic 3D network.

As shown on Fig. 7, the orientation of a given family with respect to a fixed basis (e 1 , e 2 , e 3 ) of R 3 is determined by three Euler angles : two angles (spherical coordinates θ and φ) are needed to define the normal orientation n and a third angle ψ (between the major axis and e θ on Fig. 7) allows to define the rotation around n of the fracture plane. Three families are considered in this example. The aspect ratio ω i (see Fig. 3) and the tangential permeability are the same for all the families ω = 10 -3 . The other characteristics are the following:

family 1: θ 1 = 75 1 2 ; φ 1 = 0 1 2 ; ψ 1 = 0 1 2 ; η 1 = 0.5 -family 2: θ 2 = 90 1 2 ; φ 2 = 90 1 2 ; ψ 2 = 0 1 2 ; η 2 = 1. -family 3: θ 3 = 90 1 2 ; φ 3 = 45 1 2 ; ψ 3 = 45 1 2 ; η 3 = 0.5 The three family orientations and shapes are represented on Fig. 8. In addition, the three fracture volume fractions are assumed to be equal i.e. one third of the total microfracture volume fraction. The effective permeability of such a fractured medium is computed for different values of the latter (from 0 to 6‰). Moreover, the same computations are carried out in the presence of a family cross-cutting the r.v.e. of normal e 3 (horizontal plane) and of fixed volume fraction 1‰. According to the conclusions of section 3, the normal permeability of the cross-cutting family is zero whereas the tangential permeability is chosen of the same order of magnitude as that of the micro-fractures. In both cases (with or without the family of long fractures), the eigenvalues of K hom are plotted against the total fracture volume fraction on Fig. 9. As regards the case without the family of long fractures, a percolation threshold is detected as expected. Beyond a fracture volume fraction of about 3.1‰, the three eigenvalues of K hom reach the order of magnitude of ω K t and the principal directions remain constant. It is worth precising that the direction related to the highest permeability is characterized by the angles θ = 9 1 2 and φ = 174 1 2 , which is close to the vertical axis e 3 . In the presence of the cross-cutting family, the percolation threshold disappears and three phases can be observed on Fig. 9:

-f < 3‰
The permeability is essentially controlled by the long fractures i.e. the directions related to the two highest permeabilities are close to the direction of the plane (e 1 , e 2 ) but the micro-fractures play nevertheless a role. Indeed, in the orthogonal direction, the permeability starts to increase as soon as the micro-fracture volume fraction increases. This result can be somehow justified by the fact that even few micro-fractures create a connection between the cross-cutting planes because the spacing between these ones is of the same order of magnitude as the length of the micro-fractures. -3‰≤ f ≤ 4‰ This phase is a transition during which the effects of the crossing fractures and those of the microfractures are of the same order of magnitude. The computations also show a progressive change of orientation of the principal directions. In particular, the direction related to the highest permeability initially in the horizontal plane rotates and tends towards the major direction of the case without cross-cutting fractures.

f > 4‰

The permeabilities linearly increase with the fracture volume fraction. Due to the density of the micro-fractures, the effect of the cross-cutting family gets slightly attenuated. Nevertheless, it has still to be considered within the abscisse range of Fig. 9 for the two lowest permeabilites. Indeed, the associated directions are close to the horizontal plane in which the cross-cutting family still plays a role for f < 7‰.

• 2D case : comparisons with computations on explicit discrete fracture networks The same methodology can be used in the 2D case for which we consider for example two families of micro-fractures (modeled as flat cylinders of elliptic section). The angle between the orientations is 45 1 2 and the lengths, the aspect ratio and the densities are identical. The permeability is computed for several values of the total density and compared to numerical simulations on explicit discrete fracture networks by means of the software FRACA developed by IFP. Among many functionalities aiming at analyzing the characteristics of fractured reservoirs, this software allows to randomly generate fractures according to length, orientation and conductivity laws and compute an effective permeability according to the methodology presented in [START_REF] Bourbiaux | A rapid and efficient methodology to convert fractured reservoir images into a dual-porosity model[END_REF]], [START_REF] Sarda | Hydraulic characterization of fractured reservoirs: simulation on discrete fracture models[END_REF], [START_REF] Lange | Hydraulic characterization of faults and fractures using a dual medium discrete fracture network simulator[END_REF] and [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models : Implementation and benefits[END_REF]. In the present case, the length, orientation and conductivity laws are deterministic. The number of generated fractures is computed with respect to each required density while the length remains constant and the fracture centers are randomly chosen on a square (see Fig. 10). In FRACA, the fractures are considered as an infinitesimal space between rectangles, the aspect ratio being therefore defined as the ratio between the aperture and the length. In order to ensure a relevant comparison between the upscaling model (elliptic profile of the fractures) and the numerical simulation (rectangular profile), the length of the rectangular profile, expected to have a straightforward influence on the percolation threshold, is equal to the major diameter of the corresponding ellipse. Besides, the aperture is such that both areas are identical, which means that the aspect ratio of the rectangle profile over the aspect ratio of the corresponding ellipse is equal to π/4. The two principal permeabilites provided by the effective medium theory (the computations are quasiinstantaneous on a 2.8 GHz PC) and the numerical simulations (each cross requires tens of seconds) are plotted on Fig. 11. It appears that both methods give approximately the same value of percolation threshold (about ǫ = 1.85). Nevertheless, the effective medium theory provides a linear evolution of K hom beyond the threshold, which is not expected from the percolation theory [Sahimi, 1995] and the numerical computations here. For higher densities (about twice the threshold value), the theory underestimates the minimal principal permeability computed from the explicit network but both methods give concordant estimates of the maximal permeabilities in the chosen density range (not much than three times the threshold). However, for much higher densities (ǫ ≫ 1), the spacing between fractures (2) becomes very small compared to the length if the latter remains constant i.e. the increase of density is only due to an increase of the fracture number. In that case, the fractures should not be considered as micro-fractures but long ones with respect to a lower scale r.v.e.. Moreover, one can show on the numerical computations of Fig. 11, that the slopes of the permeabilities tend towards those corresponding to a network of crosscutting fractures recalling that K hom is then analytically obtained by the second term of (37).

Conclusion

The computation of the macroscopic permeability of a fractured medium is not an easy task. Indeed, natural networks may show spatial correlations preventing the definition of a r.v.e. on which effective properties could be determined. Nevertheless, in the case of a random fractured medium such that a r.v.e. exists, an estimate of the macroscopic permeability is derived in this paper. Two kinds of fractures have been considered : fractures cross-cutting the r.v.e. and micro-fractures. The determination of the macroscopic permeability of such a medium with several families has been proposed by means of the self-consistent scheme. The application of the latter relies first on the substitution of the fractures by equivalent Darcy media with a tangential permeability deduced from the solution of a Poiseuille flow between two planes and an undefined normal permeability. An important step of this scheme consists in computing the second-order Hill tensor. The algorithm leading to the expression of this tensor has been provided for any anisotropic reference medium and any ellipsoidal inclusion. Implementing then the self-consistent scheme, it has been shown that the relevance of the fictitious normal permeability depends on the fracture type: for fractures cross-cutting the r.v.e., any value can be taken whereas for micro-fractures, the same value as the tangential permeability may be chosen not to artificially prevent connections between fractures. Results have finally been provided in the 2D and 3D isotropic and anisotropic cases. The effective permeability of a 2D anisotropic network has been compared to numerical simulations with a good agreement. The need to take carefully into account the characteristic length of each fracture family has also been put in evidence, being understood that all the phases taking part to a given upscling process have all their characteristic length of the same order of magnitude. If the density is of the order of 1 or less, the characteristic length is the fracture length and the family is considered as a micro-fracture family whereas if the density is large, the characteristic length is the spacing and the family is cross-cutting the r.v.e..

Appendix. Hill tensor in the 3D anisotropic case

This paragraph aims at providing the expression of the Hill tensor P I o related to an ellipsoid I of any shape characterised by the second-order tensor H (16) embedded in an infinite medium of any anisotropic permeability [START_REF] Dormieux | Microporomechanics[END_REF]:

P I o (x) = hess x I G o (x -x ′ ) dΩ x ′ = ∂ 2 ∂x i ∂x j I G o (x -x ′ ) dΩ x ′ e i ⊗ e j ( 66
)
where the integration is performed with respect to the variable x ′ whereas the differentiations concern the variable x. It will be shown in the following that P I o (x) is uniform over the ellipsoid I, which constitutes the transposition of Eshelby's fundamental result [Eshelby, 1957] to the second-order tensor problem class. The determination of P I o (66) requires to solve (65). To this end, it is first convenient to simplify (65) by introducing the orthogonal tensor as well as the new variable of R 3

Q diagonalizing K o t Q • Q = 1 and t Q • K o • Q = D =
y = R -1 • x = D -1/2 • t Q • x =   3 j=1 1 K o j e j ⊗ e j   • t Q • x (69)

Figure 4 :

 4 Figure 4: Matrix ellipsoid shape

Figure 7 :

 7 Figure 7: Angular coordinates of a fracture

Figure 9 :

 9 Figure 9: Eigenvalues of the permeability of a 3D network with and without a family of long fractures crosscutting the r.v.e.

Figure

  Figure 10: 2D discrete fracture network

Figure 11 :

 11 Figure 11: Permeability of a 2D fracture network of two families -angle between the two directions=45 1 2

  K o . Let us first introduce the Green function G o for the reference permeability K o i.e. the function satisfying: div x (K o • grad x G o ) = δ 0 and lim x →∞ G o (x) = 0 (65) where δ 0 denotes the Dirac distribution at point 0 of R 3 . It can be shown that the Hill tensor P I o can write by means of G o [

  e j ⊗ e j and R = Q • D 1/2 (68)

Acknowledgements

The author gratefully acknowledges fruitful discussions with Jean-Marc Daniel, Marie-Christine Cacas and Martin Guiton (IFP).

where the obvious notation D -1/2 = D 1/2 -1 has been used. Using the change of variable (69), the left hand side of the first equation of (65) writes:

In addition, a classical result of the theory of distributions allows the same change of variable on the right hand side:

Finally, the first equation of ( 65) write:

The solution of ( 72) satisfying the remote condition ( 65) is known:

Introducing ( 73) in (66), it comes:

or, with x = R • y and x ′ = R • y ′ :

where the integration is performed over the ellipsoid J of equation:

The classical newtonian potential is recognized in (75):

The potential Φ can simply be expressed by means of the eigenvectors and eigenvalues of the characteristic tensor t H ′ • H ′ of J . Let us then introduce the orthogonal tensor S such that:

Using these notations, it can be shown that the potential (77) is a quadratic form of y for y ∈ J i.e. x ∈ I ( [Kellogg, 1929], [Eshelby, 1957]):

where ∆ I , I A , I B and I C are defined with respect to A, B and C in ( 20)-( 29). Finally, the homogeneity of P I o within I is proved by introducing (79) in (77):