Ahmad Pouya
email: ahmad.pouya@enpc.fr

Minh-Ngoc Vu

Fluid flow and effective permeability of an infinite matrix containing disc-shaped cracks

Keywords: porous media, crack, steady-state flow¸ Poiseuille's law, Darcy's law, effective permeability, upscaling

A basic equation governing the steady-state flow around a single crack in an infinite porous body is given. The flow through the crack obeys to the Poiseiulle's law and the matrix has an anisotropic permeability. A semi-analytical solution is established for this equation in the case of elliptical disc-shaped crack. This solution takes a closed-form expression for the case of superconducting circular cracks. The results are compared to those obtained for flattened ellipsoidal inclusions obeying to the Darcy's flow law, which are in some works supposed to represent the cracks. It is shown that the flow solution for an elliptical disc-shaped crack obeying to the Poiseuille's law is different from that obtained as the limiting case of flattened ellipsoidal inclusions. The results are then used to establish dilute Mori-Tanaka and selfconsistent estimates of the effective permeability of porous media containing Poiseuille's type elliptical cracks.

Introduction

The effective permeability of fractured rocks and micro-cracked porous materials is a great interest for many environmental and industrial applications such as hydrology, petroleum engineering, nuclear waste disposal, geothermal expansion…. Numerous studies dealt with permeability of fracture network in three-dimensional media. Koudina et al. [ [START_REF] Koudina | Permeability of three-dimensional fracture networks[END_REF] presented a method to generate a triangular mesh on three-dimensional network of polygonal fracture surfaces, and used the finite volume method to compute the permeability of fractured media. The numerical tool developed in this way was employed to study the effective permeability of fractures network with statistical distribution [ [START_REF] Mourzenko | Macroscopic permeability of three-dimensional fracture networks with power-law size distribution[END_REF]. However, the rock matrix was supposed to be impervious in these works which does not correspond to some real cases. Therefore, other theoretical and numerical investigations have focused on this problem [[3, [4[5[6, [7, [8]. A fundamental problem for modelling the effective permeability of cracked materials is the determination of the steady state flow around a single crack in an infinite homogeneous matrix. In this work, we are interested in the case of cracks and fractures that can be represented geometrically by zero thickness surfaces with non pressure jump between their two faces, which corresponds to an assumption of zero cross section resistance for the fracture. Therefore, we do not consider the cases of non-equilibrium flow [ [START_REF] Simunek | Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone[END_REF] or of impervious fracture [ [START_REF] Manzocchi | Faults and fault properties in hydrocarbon flow models[END_REF]. In many physical studies, the Poiseuille's law is chosen preferentially to represent the flow in the cracks. The paper devotes to determine flow solution around a single crack obeying to the Poiseuille's law that is the key issue for determining the effective permeability.

The flow in cracks, embedded in an infinite matrix, has been determined in some works by using a model of ellipsoidal inclusion, in which the flow is governed by Darcy's law [ [5[11]. It is implicitly supposed in these works that the flow in the crack with the Poiseuille's law is obtained as the limiting case of a flattened ellipsoidal inclusion with thickness tending to zero. In this study, we show that this is not true, .the solution obtained for the flow in and around an elliptical disc-shaped crack obeying to the Poiseuille's law is different from that obtained as a limit of ellipsoidal inclusions obeying to Darcy's law. Only some partial results, for instance, the pressure field in the case of superconducting cracks, are similar in the two problems.

To determine the solution of flow in a matrix containing a single elliptical disc-shaped crack, we use recent theoretical results established for the more general problem of multiple intersecting fracture surfaces [ [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF]. The potential solution obtained for this general problem, as an extension of preceding works for the 2D case [[13 [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF], is based on singular integral equations allowing to reduce the dimension of the problem from 3D to 2D and, hence, particularly interesting for numerical methods. However, it will be seen in the sequel that for the case of a single elliptical disc-shaped crack, it allows also deriving semi-analytical solutions as an alternative way to the use of harmonic functions [START_REF] Green | The distribution of stress in the neighborhood of a flat elliptic crack in an elastic solid[END_REF][START_REF] Selvadurai | On the hydraulic intake shape factor for a circular opening located at an impervious boundary: Influence of inclined stratification[END_REF] or of dual integral equations methods [START_REF] Sneddon | The distribution of stress in the neighborhood of a crack in an elastic solid[END_REF][START_REF] Selvadurai | Advective transport from a penny-shaped crack and an associated uniqueness theorem[END_REF].

The solution obtained for the single crack will be then used to estimate the effective permeability of the cracked porous materials by some classical upscaling methods.

Notations :

In the sequel, light-face (Greek or Latin) letters denote scalars; underlined letters denote vectors, bold-face letters designate second rank tensors or double-index matrices. The scalar product of two vectors a and b is labelled as a.b. For second rank tensors, the tensor transposed from A is denoted A T , the matrix product is labelled as AB and the determinant as |A|. The operation of A on a is labelled as A.a. The convention of summation on repeated indices is not used for Latin indices (i, j, k…) that number surfaces, lines, etc., and are noted indifferently as subscript or upperscript. ∇ represents the gradient and ∆ the Laplace operator for a scalar field and (∇.) the divergence for a vector field.

Governing equations

An infinite porous 3D body Ω containing an elliptical disc-fracture surface Γ is considered (Figure 1). The points of Ω are denoted by x = (x 1 , x 2 , x 3 ) and, when they are located on the crack surface, by z = (z 1 , z 2 , z 3 ). The surface Γ is defined by a function z(s) from 2 → 3 parameterized by s = (s 1 , s 2 ). Fluid flow through the matrix is governed by Darcy's law:

∀x∈Ω-Γ; v(x)= -k.∇p(x) (1)
where, k is the permeability tensor of porous matrix supposed to be uniform; v(x) and p(x) the local velocity and pressure, respectively.

The local mass balance in the matrix reads:

∀x∈Ω-Γ; ∇.v(x)= 0 (2)
Figure 1: An infinite porous medium containing a flat elliptical disc shaped crack

The flow in the fracture surface is characterized by a pressure field p(s) and an infiltration field q(s) in this surface. The total infiltration q(s) through a fracture section is determined by the integral of fluid velocity v over this physical section of the thickness 2e as follows:

∀s∈Γ ( ) ( , ) e e q s v s y dy - = ∫ (3)
The crack-matrix mass exchange at a point on the fracture surface excluded its boundary is established by considering the mass balance in a small volume surrounding this point and reads [ [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF]:

∀x∈Γ-L; ( ) . ( ) . ( ) 0 s v x n s q s + ∇ = (4)
where, . ( ) s q s ∇ is the surface divergence and L the set of points on the crack boundaries. In addition, analysis of mass balance around a point on the crack boundary leads to the following mass balance condition at this point: ∀z(s)∈L;

( ). ( ) 0 q s m s =

In this equation, q is the infiltration vector in the crack and m is the outward unit normal vector of the boundary line that is tangent to the crack plane.

The body is submitted at its infinite boundary to a pressure field p ∞ (x):

lim x →∞ [p(x) -p ∞ (x)] = 0 (6)
p ∞ satisfies ∆p ∞ = 0. For permeability upscaling, a constant pressure gradient is supposed to be prescribed to the matrix at its infinite boundary. Designating by A the far-field pressure gradient, we have:

p ∞ (x) = A.x (7)
The equations to be satisfied by p are the boundary condition [START_REF] Pouya | Permeability of 3D discontinuity networks: New tensors from boundary-conditioned homogenisation[END_REF], mass balances (2), ( 4), [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF] and the constitutive equation [START_REF] Koudina | Permeability of three-dimensional fracture networks[END_REF]. The solution obtained by Pouya [[12] expresses the pressure field in the whole matrix as a function of the infiltration vector in the fracture as follows:

( ) ( )

1 3/2 1 ( ). . ( ( )) 1 ( ) ( ) 4 ( ) . . ( ) q s x z s p x p x d s x z s x z s - ∞ - Γ - = + π   - -   ∫ k k k (8)
In this equation, ds represents the surface element on Γ. This solution satisfies well the farfield condition [START_REF] Pouya | Permeability of 3D discontinuity networks: New tensors from boundary-conditioned homogenisation[END_REF]. Also, noting that the integrant function in the right-hand side integral is equal to ( ). (1 / )

x q s r -∇ with r = [(x-z).k -1 .(x-z)] ½ , the mass balance (2) can be shown for x non situated on the fracture. To establish that ( 4) is well satisfied by ( 8) is little more technical and we refer to [ [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF] for its demonstration.

It should be noted that the only significant assumptions limiting the validity of this solution are that matrix permeability k is assumed to be uniform and that there is no pressure gap between the two sides of the crack surface. Although we consider in the sequel only the case of a plane crack surfaces, the solution (8) remains valid for curved-surface fractures, and, in its general form given in [ [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF], for multiple intersecting fractures. Also, in equation ( 8), the constitutive flow law in the fracture is not constrained. A linear relation between the infiltration vector and the pressure gradient, as an extension of the Poiseuille's law for the viscous fluid flow between two parallel plates, is widely used in the literature to represent the flow in the crack:

( ) ( ). ( ) s q s s p s = - ∇ c ( 9
)
The conductivity tensor c depends on the fluid dynamic viscosity µ, the hydraulic aperture of the fracture e, and also the roughness of fracture surfaces. In most cases, c is supposed isotropic (in the fracture surface) and replaced by a scalar c. In particular, the Poiseuille's law, widely used for flow through rock fractures [ [START_REF] Guéguen | Introduction to Physics of Rocks[END_REF], determines c from a laminar flow between two parallel plates separated by a constant distance e and leads to c = e 3 /(12µ). However, c can be anisotropic in some cases as shown experimentally by Gentier et al.[ [START_REF] Gentier | Anisotropy of flow in a fracture undergoing shear and its relationship to the direction of shearing and injection pressure[END_REF]. Introducing ( 9) into (8) leads to an integral equation involving only the unknown pressure field. A simple analysis shows that the solution for a constant pressure gradient along the disc surface, corresponding to the Eshelby-type ellipsoidal inclusion with Dracy's law, does not satisfy this equation. The comparison between the Posieulle-type crack and the flattened ellipsoidal inclusion with Darcy's law will be discussed further in this paper.

However, in this paper we are interested in the case of superconducting cracks corresponding to the limiting case c→∞. This is an important limiting case, currently addressed in the literature because it allows determining the upper bound of the effective permeability for fixed crack density and geometry.

Superconductive elliptical crack-inclusion

In the case of superconducting fracture, the pressure remains constant along the fracture surface. We consider an elliptical disc shaped crack D and denote by d 1 and d 2 the halfdiameters of the ellipse, by the unit vectors e 1 and e 2 its principal directions and by e 3 the unit normal to the crack plane. A coordinate system is defined with its origin located at the crack centre and its axes parallel to e 1 , e 2 and e 3 . The equation of the ellipse reads x.B.x = 1 for x.e 3 = 0. We note:

λ 1 = 1/d 1 2 , λ 2 = 1/d 2 2 , B = λ 1 e 1 ⊗e 1 + λ 2 e 2 ⊗e 2 (10)
The general result needed for permeability upscalling is the integral of the infiltration over the disc surface. This integral is a linear function of the far-field gradient A and can be written as Ψ Ψ Ψ Ψ.A, where the tensor Ψ Ψ Ψ Ψ is a function of c, k and B. By application of the linear transformation [ [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF], the problem can be changed into that of an ellipsoidal disc in an isotropic matrix. Nevertheless, a numerical method is required to solve the general problem for the case of finite and anisotropic conductivity tensor c. In the following, we consider the particular case of superconducting crack, c→ ∞, that allows us to derive an analytical or semi-analytical solution.

The equation [START_REF] Markov | Permeability of the Fluid-Filled Inclusions in Porous Media[END_REF] written for a single elliptical crack D embedded in an infinite porous matrix with isotropic permeability k yields:

( ). ( ) 1 ( ) . 4 D q x p x A x d k x ξ ξ - = - ξ π ξ - ∫ (11)
In this equation, dξ represents the surface element on D. The general solution can be obtained by superposition method. It requires three elementary solutions for A parallel to e 1 , e 2 and e 3 .

A simple analysis shows that for A parallel to e 3 the solution is a uniform gradient pressure not perturbed by the presence of the crack, i.e., p(x) = A.x in the whole body. To build the complete solution, it is sufficient to determine the solution for A parallel to the x-axis. Hence, we try to solve the problem for p ∞ (x) = A 1 x 1 .

The condition of superconducting crack imposes that the pressure p must be constant in the crack surface, and because of symmetry considerations, this constant is zero. Therefore, designating, for a field q(ξ), by:

3 ( ).( ) ( ) D q x I q d x ξ ξ - = ξ ξ - ∫ (12)
the solution q must satisfy: To find a solution for this problem, we try first the following expression which is inspired by the 2D solution given by Pouya & Ghabezloo [ [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF]:

I (q) = 4πkA 1 x 1 (13) x 2 x 1 d 1 d 2 x ξ θ ρ
0 1 1 ( ) 1 . . q h e ξ = -ξ ξ B ( 14
)
As shown in Appendix A, the expression found for I is a linear function of x 1 for q = q 0 which allows us to build the solution for the equation [START_REF] Liolios | A solution of steady-state fluid flow in multiply fractured isotropic porous media[END_REF]. The general expression found for this solution has recourse to the complete elliptic integrals of first and second kind K(λ) and E(λ) defined by:

( ) / 2 2 2 0 1 s i n π θ= θ λ = -λ θ ∫ d K , ( ) / 2 2 2 0 1 s i n π θ= λ = -λ θ θ ∫ E d (0 ≤ λ≤ 1) (15)
The semi-analytical or iterative numerical methods used to calculate these complete elliptic integrals K(λ) and E(λ) can be found in mathematical handbooks [ [START_REF]The Process of the Arithmetic-Geometric Mean[END_REF].

Moreover, we define:

( ) ( ) /2 2 1 1 1 2 2 2 0 1 2 /2 2 2 2 1 2 2 2 0 1 2 co s , co s s i n s i n , co s s i n d d π θ= π θ=  λ θ φ λ λ = θ  λ θ + λ θ    λ θ  φ λ λ = θ  λ θ + λ θ   ∫ ∫ (16) It can be noticed that, if d 1 ≥ d 2 , denoting by 2 2 1 1   λ = -    d d , we have: 1 1 2 2 1 2 2 2 2 1 1 1 2 2 2 1 2 2 1 ( , ) ( , ) ( ) ( , ) ( , ) ( ) E d d d d K  φ λ λ + φ λ λ = λ    φ λ λ + φ λ λ = λ  (17)
This system of equations allows the determination of φ 1 and φ 2 as function of K(λ) and E(λ). Then, for q = q 0 , it is shown in Appendix A that:

I = -2h 1 π φ 1 (λ 1 ,λ 2 ) x 1 ( 18
)
As a result, ( 13) is verified if we take 1

1 1 2k h A - = φ
, and therefore:

0 1 1 1 2 ( ) 1 . . k q A e ξ = - -ξ ξ φ B ( 19
)
According to this expression, the flux q 0 in the crack has a constant direction parallel to the far-field gradient, as in the case of ellipsoidal Darcy type inclusions. However, q 0 is not the real solution of the flow in the disc. As a matter of fact, it must be noted that the solution of the equation ( 13) is not unique. If the vector ω(ξ) satisfies ∇.ω = 0 in the disc and ω.n = 0 on its boundary where n is outward unit normal vector to the boundary, then we calculate:

( )

1 1 1 ( ) ( ). ( ). . 0 I x d x n ds x d - - - Γ ∂Γ Γ ω = -ω ξ ∇ ξ - ξ = -ξ - ω ξ + ξ - ∇ ω ξ = ∫ ∫ ∫ (20)
Consequently, q = q 0 + ω also satisfies [START_REF] Liolios | A solution of steady-state fluid flow in multiply fractured isotropic porous media[END_REF]. This degeneracy of the solution can be removed if we suppose that the superconducting crack is the limit case of a crack with the conductivity c = σc 0 where c 0 is constant and σ→∞. This implies that c 0 -1 q must be a gradient field (pressure gradient) and satisfy the condition of symmetry of crossed derivatives. This condition allows us to determine a unique q saved that c 0 is symmetric and positive definite. In the case of a superconducting crack considered as the limit case of an isotropic conductivity crack, this condition implies that q has to be a gradient field. It can be checked that q 0 given by ( 19) is not a gradient field. A complementary part ω(ξ) should be added to q 0 to make it a gradient field satisfying [START_REF] Liolios | A solution of steady-state fluid flow in multiply fractured isotropic porous media[END_REF]. Thus, the solution of the infiltration field q reads:

q(ξ) = q 0 (ξ) + ω(ξ) (21)
where, the complementary part ω has to satisfy the following conditions:

∀ ξ∈D; ∇.ω = 0 (22) ∀ ξ∈∂D; ω.n = 0 (23) ∀ ξ∈D; ∂ 1 (q 0 2 +ω 2 ) =∂ 2 (q 0 1 +ω 1 ) (24)
The last condition assures that q is a gradient field. Mathematically, equation ( 22) is satisfied if we take:

ω 1 = ∂ 2 ϕ, ω 2 = -∂ 1 ϕ (25)
with ϕ an arbitrary function. Then, the condition (24) yields:

1 2 1 2 1 . . s kA ∆ ϕ = ∂ -ξ ξ φ B (26)
where, ( ) 25) into [START_REF] Gervais | Modélisation géométrique d'un réseau de fractures dans un massif rocheux stratifié. Application aux carrières marbrières de Comblanchien (Côte d'Or[END_REF] and rearranging the terms, we find:

2 2 1 2 s ∆ ϕ = ∂ + ∂ ϕ . Substituting (
∀ ξ∈∂D; ∇ϕ.t =0 ( 27
)
where t is the unit tangent vector to the crack boundary ∂D. This equation implies that ϕ is constant on the boundary, hence:

∀ ξ∈∂D; ϕ = ϕ 0 (28)
with, ϕ 0 a constant. A solution for the Poisson's equation ( 26) over the elliptical domain exists for every Dirichlet boundary condition and in particular (28) with constant ϕ 0 . We can take ϕ 0 =0 because the addition of a constant value to ϕ does not change the solution ω. The solution of equation ( 26) with ϕ=0 on the boundary can be determined numerically in the general case and analytically in some special cases. Two particular cases in which the solution can be derived analytically are the following:

Case 1: Limit of infinitely long elliptical disc

We suppose the limit case of elliptical disc with infinite extension in the direction 2. Geometrically, this corresponds to an infinite band of width 2d 1 parallel to the direction 2.

The solution, for a far-field gradient A 1 parallel to e 1 corresponds to that of the 2D case with a crack of length 2D in the plane (x 1 , x 2 ) parallel to x 1 . The solution can therefore be compared to that given for this case by Pouya & Ghabezloo [ [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF]. In the 3D problem, d 2 →∞ and λ 2 = 0 and

0 2 2 1 1 1 1 ( ) 2 q k d A e ξ = -
-ξ is a gradient field. Thus, ω = 0. The exact solution is then

2 2 1 1 1 1 ( ) 2 q k d A e ξ = -
-ξ which is exactly the solution found in 2D case by Pouya & Ghabezloo [ [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF].

Case 2: Circular crack with radius R

In this case,

d 1 = d 2 = R. We find 1 2 4R π φ = φ = and 0 2 2 1 1 8 ( ) kA q R r e ξ = - π with r 2 =ξ.ξ.
In order to find the expression of the complementary function ω, we look for a solution of (26) in the following form:

² ϕ(ξ) = χ(ξ) + ∂ 2 ζ(r) (29)
where, ζ(r) satisfies:

2 2 1 8 s kA R r ∆ ζ = - - π (30)
Then, χ has to satisfy ∆χ =0. The general solution of (30) is:

( ) ( ) 2 2 3/ 2 2 2 2 3 2 2 1 1 2 2 2 8 1 1 1 3 6 9 kA R R r r R R r R Ln R r C Lnr C R R r   - - ζ = - - -- - - + +   π + -     (31)
where, C 1 , C 2 are two constants. Thus, ( )

3/2 2 2 1 2 2 1 2 8 1 3 kA R r C r ξ   ∂ ζ = - - -   π   (32)
Substituting this expression in (29) and taking into account the boundary condition ϕ(R)=0, we find that ( )

3/2 2 2 1 2 1 2 1 1 2 2 8 8 1 3 kA kA R r C C r R ξ ξ   ϕ = - - - -   π π   (33)
The constant C 1 is determined by the condition that q has to be finite in the crack or that ϕ is not singular when 0 → r . This imposes

3 1 3 R C =
. Substituting by this expression in ( 25) and ( 21) we find the solution of infiltration in the circular disc under the far-field pressure gradient A 1 as follows:

( ) ( ) ( ) 2 2 2 3 2 2 3 2 2 1 1 1 2 1 2 4 3 1/ 2 3 2 2 2 2 1 1 2 1 2 2 4 2 8 ( ) 3 3 8 2 3 kA R q R r R R r e r r kA R R r R r e r r   ξ ξ -ξ   ξ = - -- - - +     π     ξ ξ ξ ξ     - - - - + -     π     (34)
The two components q 1 (ξ), parallel to the far-field gradient and q 2 (ξ), normal to this gradient, are represented in the Figure 3 for A 1 =1 and k =1.

Figure 3: Closed-form solution of infiltration q in fracture of unit radius in an infinite isotropic matrix under a far-field pressure gradient in x 1 -direction: q 1 component parallel to the far-field gradient, q 2 orthogonal to the far-field gradient

The elementary solution for a far-field pressure gradient A 2 parallel to e 2 is obtained from (34) by index changing 1↔2. As mentioned above, the solution for the general case of a far-field pressure gradient A is obtained by superposition of the three elementary solutions. By introducing the notation:

.

r = ξ ξ , 1 ˆr ξ = ξ (35) T = δ δ δ δ -n ⊗ n ( 36
)
where n is the unit normal vector to the crack plane, the general solution obtained for q has the following expression:

( )

3 3 2 3 2 2 2 2 2 3 2 2 2 8 ˆ( ) 3 2 2 3 k q Rr R R r r R r R R r A r -      ξ = -+ - + -+ - - ξ ⊗ ξ       π       T δ δ δ δ (37)
It is important to note here a fundamental difference between this result and what could be obtained as the limiting case of flow solution in a flattened ellipsoidal inclusion. The comparison between the two models is presented farther in the paper. But it is already possible to note the flow in the elliptical crack given by (37), as seen also in the Fig. 3, is not uniform and parallel to the farfield gradient contrary to the case of ellipsoidal inclusions in which the fluid velocity is found to be uniform and parallel, at least if the matrix and inclusion are isotropic, to the farfield gradient. More precisely, for instance for a far-filed pressure gradient parallel to the x-axis, the solution given by (37) has a non vanishing component in the y-axis direction as shown by the expression (34) (see the figure 3 in the paper). However, for the same problem modelled by a flattened spheroidal inclusion, the fluid velocity in the inclusion would be parallel to the x-axis, so with non component in the y-axis direction. This shows the fundamental difference between the two types of flow and the fact that the solution for the Poiseuille crack can not be obtained as the limiting case of ellipsoidal inclusions with Darcy's law.

The solution (37) could be probably obtained by other mathematical methods for instance based on the use of harmonics functions instead of the potential solution [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF]. But for the difference underlined here-above, it can not be obtained as a limiting case of solutions presented in the literature for ellipsoidal inclusions obeying Darcy's law.

The solution (37) represents an original and fundamental result that opens the way for further theoretical investigations concerning flow in cracked porous materials.

Pressure field in the matrix

The solution of the infiltration q through the crack in (37) allows the determination of the pressure field in the matrix by using equation [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF].

For the discussions carried on farther in this paper, it is useful to decompose q in q 0 and ω parts according to [START_REF]The Process of the Arithmetic-Geometric Mean[END_REF] and to write:

0 3 ( ) ( ) .( ) 1 ( ) . 4 D q x p x A x d k x   ξ + ω ξ ξ -   = - ξ π ξ - ∫ (38)
The equation [START_REF] Gentier | Anisotropy of flow in a fracture undergoing shear and its relationship to the direction of shearing and injection pressure[END_REF] shows that the contribution of ω to the integral vanishes. Besides, we can deduce from [START_REF] Guéguen | Introduction to Physics of Rocks[END_REF], by superposition and by noting that a far-field gradient parallel to e 3 does not induce flow in the fracture, that, for a general far-field gradient A, q 0 is given by: 0 ( ) 2 1 . . .

q k A ξ = - -ξ ξ B Φ Φ Φ Φ (39)
where,

1 1 1 2 1 1 2 2 = ( ) ( ) e e e e - - φ ⊗ + φ ⊗ Φ Φ Φ Φ (40)
Substituting ( 40) into (38) yields:

3 1 . . 1 ( ) . ( ) . 2 D p x x x d A x   -ξ ξ   = - ξ - ξ   π ξ -   ∫ B Φ Φ Φ Φ (41)
The pressure field solution given in this way is an explicit integral function of the geometrical parameters of the elliptical disc contained in the tensor B. It is interesting to notice that this solution is independent of the matrix permeability k. This is due in particular to the assumption of superconducting crack: the pressure solution would depend on k for crack with finite conductivity. Anyway, the equation (41) allows us to determine at least numerically the pressure field around an elliptical crack by direct integration.

For the case of circular disc (B = R -2 T) a closed-form expression for the integral (41) can be derived for a point located on the plane of the crack, i.e. x 3 = 0. We find (see Appendix B):

( )

3 1 . . ( ) sin 2 D x d x R x -ξ ξ π ξ - ξ = ψ -ψ ξ - ∫ B ( 42
)
where, ψ is, as shown in Figure 4, the angular diameter of the disc seen from the position x.

x 2 Due to superconducting crack, p(x) = 0 when x is in the crack. We also have Φ Φ Φ Φ = (4R/π)T and the expression of p for the points exterior to the disc is given by the following expression:

x 1 R x ξ ϕ ρ ψ u v ρ 2 ρ 1
x.e 3 = 0 ;

( )

1 ( ) 1 sin . p x A x   = -ψ -ψ   π   (43)
It can be checked that when x → ∞ then ψ→0 and so p(x) → A.x according to (43). The infinite boundary condition is recovered in this way. Also, when x is on the boundary of the disc, ψ = π/2 and we find the pressure condition in the disc, p(x) = 0. Figure 5 displays the variation of pressure field in the plan containing the superconducting fracture with a far-field pressure gradient A=(1,0,0). It could be noted that this pressure solution could be, in principle, obtained as the limiting case of the pressure field in an infinite matrix containing a superconducting ellipsoidal inclusion. The two problems obey, as a matter of fact, to the same equations and boundary conditions: constant pressure in the crack or inclusion, constant farfield pressure gradient and mass conservation equation in the matrix. The solution could be then obtained as a classical application of potential theories using harmonics functions. However, we have not found in the literature any explicit expression of this solution as given by (43).

Comparison between Darcy-ellipsoidal inclusions and Poiseulle-crack

The crack inclusion has been sometimes considered as the limiting case of a penny-shaped inclusion with aspect ratio (ratio between its thickness and diameter) tending to zero. However, a further investigation shows this similarity holds only for some partial results and in restricted cases and this difference seems to be ignored in the literature.

To show this difference, let a circular penny-shaped inclusion be considered with the revolution axis taken as the third axis of coordinates and with semi-diameters

d 1 = d 2 = d, d 3 << d.
The inclusion has an isotropic permeability k * embedded in the matrix of isotropic permeability k. The fluid velocity v in the inclusion is found to be uniform and its value can be deduced, for instance, from the results given by Shafiro and Kachanov [ [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF]. We consider the limit when ε = d 3 /d→0 and k * →∞. The contribution of the inclusion to the equivalent permeability of the material is made through the integral of the fluid velocity in the inclusion:

* ( ) Q v d Ω = ξ ξ ∫ ( 44
)
This integral is a linear function of the far-field gradient A:

Q = -Ψ Ψ Ψ Ψ.A (45)
Shafiro & Kachanov [ [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF] gave the following general expression of Ψ Ψ Ψ Ψ for an arbitrary value k * of the inclusion and ε:

1 1 * * * * ( ) 1 ( ) 1 1 4 2 k k k k k k n n n n k k - -       Ω -πε - πε     = - + ⊗ + + - ⊗         Ω           Ψ δ - Ψ δ - Ψ δ - Ψ δ - (46)
where, Ω and Ω * represent respectively the matrix reference volume and the inclusion volume.

Considering the limiting case of a superconducting crack by taking the limit ε →0 and k * →∞ Shafiro & Kachanov [ [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF] proposed the following expression for this limit:

3 1 6 ( - ) 3 R k n n ν ⊗ Ψ = δ Ψ = δ Ψ = δ Ψ = δ (47)
where, ν is the crack density. Nevertheless, a careful study shows that the limit is not unique and depends on the value of εk * which is degenerate when ε →0 and k * →∞. Let us suppose that k * =bk/ε where b is a fixed ratio, hence, εk * = bk remains constant when ε →0 and k * →∞.

The limit of Ψ Ψ Ψ Ψ from the Shafiro & Kachanov [ [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF], is obtained to be:

3 16 ( - ) 3 (4 / ) ν ⊗ + π I b R k n n b Ψ = Ψ = Ψ = Ψ = (48)
The limit expression (47) corresponds then to the assumption b →∞. At the boundary, we note v m the velocity in the matrix and v i the velocity in the inclusion. The velocity in the inclusion is constant for ellipsoidal inclusion. We note v i = ve 1 . At the boundary, the mass balance yields v i .n = v m .n, where n is the unit normal vector to the ellipsoid, and then, v m .n = vn 1 where

1 1 2 2 2
(1 )

x n d r ε = - -ε .
The general value of velocity v the inclusion, for small ratio ε=d 3 /d , can be deduced from the expression given by Shafiro & Kachanov [ [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF] for the integral (44) since v is constant in the inclusion. We find, in first order development with respect to ε:

* * 1 4 k k v k k k - = - πε - + (49)
The equivalent infiltration vector in the crack section is deduced from the relation [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF] where the thickness e is given, for the ellipsoidal inclusion, by [START_REF] Mourzenko | Macroscopic permeability of three-dimensional fracture networks with power-law size distribution[END_REF] 2 e d r = ε and v is constant and given by (49). The limit of the equivalent infiltration q = 2ev for ε→ 0 is found to be:

2 2 0 8 li m 8 kc q d r c d k ε→ = - - π + (50)
where, c = 2d 3 k * . We note that for the limit case c→∞, the equation (50) yields a solution different from that given by (37). As a matter of fact, it yields only the expression of q 0 in the equation ( 19) for circular disc (φ 1 =π/4d), which is not the complete solution for infiltration in the Poiseiulle's crack since it misses the complementary part ω.

It can also be checked that the expression (50) does not satisfy the equation [START_REF] Shafiro | Anisotropic effective conductivity of material with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF]. It is then deduced that, in the general case of finite conductivity c, the solution of flow in the crack with Poiseuille's law is not given as the limit of thickness vanishing ellipsoidal inclusions with Darcy's law. This implies also that the pressure fields in the inclusion and the matrix are different for the two problems.

Only in the case of superconducting crack, the pressure fields in the two problems are the same since the pressure is constant in the crack or the inclusion. In this case, the integral of q over the inclusion (determining the contribution to the effective permeability, see the next section) is also the same for the two models since it depends only on q 0 and not on the complementary part ω.

In conclusion, the limit of the flattened inclusion with Darcy's flow does not give in the general case the same result of pressure and infiltration fields that the model of Poiseuilletype crack.

Application to the effective permeability of cracked porous materials

The effective permeability k of the heterogeneous domain is calculated from the average velocity and average pressure gradient in this domain. For a domain Ω containing a distribution of zero-thickness fractures Γ j , the equivalent average velocity is defined by:

V ≡ 1 j j v d q ds Ω Γ     + Ω     ∑ ∫ ∫ ω (51)
It can be shown that if a pressure condition corresponding to a uniform macroscopic gradient is applied on the boundary of Ω, i.e., p(x) = A.x on ∂Ω, the average pressure gradient G is equal to A:

1 G p d A Ω = ∇ Ω = Ω ∫ (52)
We assume that the permeability tensor of the porous matrix is constant and equal to k.

Introducing (1) and ( 52) into (51) and using V = -k .G yield:

k = k + k f (53)
where, k f is the contribution of the cracks to the effective permeability obtained by the following relationship:

1 j j q ds Γ Ω ∑ ∫ = -k f .A (54)
The dilute Mori-Tanaka scheme estimates the left-hand side of (54) as the sum of the integral of q over each individual crack by neglecting the effects of cracks interactions. As seen in section 3, the infiltration in the crack is expressed by q = q 0 + ω. According to [START_REF] Seyedi | A two-scale model for simulating the hydromechanical behavior of faults during CO2 geological storage operation[END_REF], ω can be written as: ω = e 3 ×∇ϕ. Using this relation and (28), we calculate:

0 3 3 ( ) . 0 D D D D d s e d s e n d l t d l ∂ ∂ ω ξ = -∧ ∇ϕ = -∧ ϕ = -ϕ = ∫ ∫ ∫ ∫ (55)
Therefore, only q 0 contributes to the integral over the crack surface. From the results given in Section 3, we deduce:

0 4 ( ) ( ) . 3 k Q q ds q ds S A Γ Γ Γ = ξ = ξ = - ∫ ∫ Φ Φ Φ Φ (56)
where, the tensor Φ Φ Φ Φ defined by (40) and

1/2 1 2 / S d d Γ = π = π B
represents the area of the ellipse Γ. Comparing this relation to (54), the Mori-Tanaka estimate of k f for an isotropic matrix containing a family of elliptical cracks is found as follows:

f 4 3 k S Γ = ν < > k Φ Φ Φ Φ (57)
where, ν is the crack density, i.e: number of cracks (centres) per unit volume and <S Γ Φ Φ Φ Φ > is the average value of S Γ Φ Φ Φ Φ calculated over the cracks population. For an isotropic distribution of circular cracks with radius R, we find:

f 3 3 2 9 k R = ν k δ δ δ δ (58)
and the following expression for the dilute Mori-Tanaka effective permeability:

3 3 2 ˆ1 9 MT k k R   = + ν     (59)
The self-consistent estimate is obtained by replacing, when calculating the contribution of cracks to the effective permeability in equation ( 58), the matrix permeability by the unknown effective permeability. Therefore, equation (59) leads to the following equation for the selfconsistent estimation ˆsc k :

3 3 2 ˆ9 s c s c k k k R = + ν (60) Furthermore, 3 ˆ3 2 1 9 s
c k k R = - ν (61)
The dilute Mori-Tanaka (59) estimation is found to give the first order development for small ν of the self-consistent estimation (61) (Figure 7). By introducing the dimensionless crack

Extension to anisotropic matrix

The linear transformation presented in Pouya [ [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF] can be used to extend the solutions obtained in the previous sections for elliptical crack in an isotropic matrix to the case of anisotropic matrices. Considering a reversible tensor M, we define the transformed variables

x ∼ , v ∼ (x ∼
) and p ∼ (x ∼ ) as follows:

x ∼ =M.x , v ∼ (x ∼ ) = M.v(x) , p ∼ (x ∼ ) = p(x) (62)
Therefore, x ∼ =M.x transforms the domain Ω to a domain Ω ∼ , the cracks Γ to a new geometry Γ ∼

The flow problem is then transformed to a new problem corresponding to matrix permeability.

Let us consider an infinite homogeneous porous medium embedding an elliptical crack with a geometry defined by two principal axes e 1 and e 2 and half-diameters d 1 and d 2 . A pressure gradient A is imposed at the infinite boundary of the domain. We refer to Pouya [[12] for the complete set of transformation relations for different variables in which the tensor of permeability is taken the following expression:

k ∼ = M k M T (63)
The transformation to the problem studied in previous section is applied by taking

M = -1
k . In the transformed problem, we find an isotropic matrix with unit permeability: k ∼ = δ. The equation of the elliptical crack ξ.B.ξ = 1 is equivalent to ξ ɶ .M - The expressions λ ɶ 1 and λ ɶ 2 , and of φ ∼ 1 and φ ∼ 1 that can be deduced from λ ɶ 1 and λ ɶ 2 by the same relations ( 16) are required for writing the solutions of pressure and velocity in the transformed problem. However, if only the pressure field in the matrix in equation ( 41) or the application to effective permeability in equation ( 56) is envisaged, it is easier to express directly Φ Φ Φ Φ as function of B (see Appendix C):

f g B Φ = δ + Φ = δ + Φ = δ + Φ = δ + (66)
where, f and g are scalar functions of B = B:δ δ δ δ and β = |B|.

As a matter of fact, the pressure transformation p ∼ (x ∼ ) = p(x) leads to A ∼ .x ∼ = A.x at infinity; hence, A ∼ = k .A. Moreover, 0 ( ) q ξ ɶ ɶ in the transformed problem is given by (39) such as:

0 ( ) 2 1 . . .Α q k ξ = - -ξ ξ Φ Φ Φ Φ ɶ ɶ ɶ ɶ ɶ ɶ ɶ ɶ B ( 67
)
where 1 k = ɶ and:

= ( , ) ( , ) f B g B β β ɶ ɶ ɶ ɶ ɶ ɶ B Φ δ + Φ δ + Φ δ + Φ δ + (68)
The general transformation rule for infiltration vectors is (see in Pouya [[12]):

0 0 1 1 . . . k k q q n n - = ɶ ( 69
)
Replacing by these values function of initial variables yields: and g given in Appendix C. In this way, q 0 (ξ) is calculated explicitly for arbitrary elliptical shaped crack in a matrix with arbitrary anisotropic permeability k. The pressure field in the matrix is deduced from:

0 ( ) 2 .k. 1 . . . k k q n n A ξ = - -ξ ξ ɶ B Φ Φ Φ Φ
( ) ( )

0 1 3/2 1 ( ). ( ) 1 ( ) . 4 . . D q x p x A x d x x - - ξ ξ - = + ξ π ξ - ξ - ∫ k k k (70)
For extension of the results concerning the effective permeability, it is sufficient to write (57) for the transformed problem:

f 4 3 k S Γ = ν < > k ɶ ɶ ɶ ɶ ɶ Φ Φ Φ Φ (71) Noting that 1 k = ɶ , / ν = ν = ν M k ɶ , / S Γ = π B ɶ ɶ , f f = k k k k ɶ
, and using the expressions obtained here above for ɶ Φ Φ Φ Φ and B ɶ , we find:

f 4 3 S Γ = ν < > k k k ɶ Φ Φ Φ Φ (72)
Substituting ( 72) into (53) allows us to determine the effective permeability of a porous body with anisotropic matrix permeability k containing a distribution of elliptical superconducting cracks by using the Mori-Tanaka or the self-consistent estimate.

In practice, field observation permits us to determine the cracks geometries (size, aperture and orientation) and also their density [ [START_REF] Gervais | Modélisation géométrique d'un réseau de fractures dans un massif rocheux stratifié. Application aux carrières marbrières de Comblanchien (Côte d'Or[END_REF]. An estimation of the crack conductivity is also required that can be obtained by a Poiseuille's model. However, in our work, only superconducting disc-shaped cracks are considered. The crack network is assumed to be randomly distributed. With these data, the models given in Section 6 are suitable to estimating the effective permeability.

Conclusions and perspectives

The equation governing the steady state flow around a crack surface in an infinite porous body, deduced from previous investigations [ [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF], was considered in this paper. An analytical solution was derived for this equation in the case of a superconducting and elliptical disc shaped crack.

Cracks have been very often modelled as flattened ellipsoidal inclusions obeying Darcy's law because for this problem analytical results are available. The question raises then if this model can represent well a crack with zero thickness and obeying to the Poiseuille's law. In the literature, it is often implicitly responded positively to this question: it is assumed that the flow field in and around a crack with Poiseille's law can be obtained as the limiting case of the flow field in and around a flattened ellipsoidal inclusion with Darcy's law when its thickness tends to zero. The results obtained in the present work allowed the comparison. We could show that the equivalence between the two models does not hold in the general case.

Only some partial results in restricted cases, mainly the case of superconducting cracks, can be obtained as a limit case. In general, for cracks with finite conductivity, the pressure gradient through the mean crack surface is constant in ellipsoidal inclusion model but not in Poiseuille's type crack. These results change fundamentally the vision we can have for the flow in a Poiseuille's type crack.

The results obtained for a single crack have been used to estimate the effective permeability of cracked porous materials. The dilute Mori-Tanaka model is a very simple model that does not take into account the interaction between cracks. The expression of the effective permeability is a linear function of the crack density according to this model, for anisotropic matrix permeability and elliptical disc cracks. The crack interaction is taken into account implicitly in the self consistent model. We established the exact expression of the effective permeability according to this model for circular cracks in a matrix with isotropic permeability. This non linear model presents a singularity for a critical value of crack density that is sometimes interpreted as representing the percolation threshold of crack network. This critical value is then compared to analogous values given in the literature, and found to be very close to that obtained by numerical simulation for the percolation threshold of crack networks.

The analytical solutions given in this paper open the way for deeper theoretical investigations.

The effective permeability of micro-cracked materials can take benefit from these results as well as the study of crack interactions within a porous body. Also it is worth reminding that a fundamental assumption underlying the present work was the infinite transversal conductivity of cracks that removes the possibility of pressure jump between to sides of the crack. But in some physical cases, cracks act as an impervious membrane with a pressure discontinuity across the crack surface [ [3, [10, [24]. The approach used by Martin [ [START_REF] Martin | On potential flow past wrinkled discs[END_REF] could be used in this case to extend some of the results of the present paper to these types of cracks.

Analytical infiltration solution through a single crack obtained here was restricted to the case of superconducting cracks. For the case of Poiseiulle-type cracks with finite conductivity, semi-analytical and numerical calculations are required that we have investigated in an ongoing work and we hope present the results in a future paper. Moreover, numerical study allows us to take into consideration explicitly the crack interaction and the crack intersection. The three-dimensional effective permeability model that can be obtained in this way would allow extending the two-dimensional application to CO2 storage [ [START_REF] Seyedi | A two-scale model for simulating the hydromechanical behavior of faults during CO2 geological storage operation[END_REF] to a three-dimensional modelling. where:

2 1 . . ( ) 1 ( ) m m d J d e g e ρ ρ ε ε ε -ξ ξ γ ρ θ = ρ = -ρ + ρ ρ ∫ ∫ B (A.3) with: 2 2 1 2 co s s i n γ = λ θ + λ θ , 2 2 1 2 1 2 ( s i n co s ) H x x = γ -λ λ θ + θ 2 1 1 2 2 co s s i n , x x e g H H λ θ + λ θ γ = = (A.4)
We take the variable α:

sinα = eρ+g (A.5)
When ρ varies between ε and ρ m , α varies between α ε = Arcsin(ρε +g) to π/2. In this interval, cosα ≥ 0, and variable change ρ→α allows calculating:

/2 2 1 ( ) (s i n ) s i n g J g d e g ε π ε α=α   γ - θ = - α + α   α -   ∫ (A.6)
and find:

[ ] 2 (1 ) ( ) ( ) co s ( / 2 ) ε ε ε ε γ - γ θ = θ - α + π -α g J F g e e (A .7)
Where: where the function φ 1 (λ 1 ,λ 2 ) is defined by the equation ( 16) in the main text.

Furthermore, using the same notation with x 1 =x 2 =0, and thus, g=0 and e=γ, we have:

Appendix B: Integral calculation for a point outside the circular disc

As seen in Figure 4 in the text, the integral in the left-hand side of (42) can be re-written in the polar coordinate system as follows:

∫

For the same variable change (A.5), α varies between -π/2 to π/2 when ρ varies between ρ 1 and ρ 2 . By the same method than in the previous case, and using (A.11) with, this time, |g|>1, we calculate: For a given B, we define: The solutions of 1 φ , 2 φ can be written as:

[ ] 2 2 1 1 2 1 1 1 . ( ) ( ) ( , ) K E B λ -λ -λ φ = λ -λ = λ χ β λ (C.3) 2 2 1 2 2 2 2 1 1 1 ( ) ( ) ( , ) K E B   λ -λ -λ φ = - λ + λ =   λ χ β λ λ   (C.4)
Remarking that two tensors Ф and B have the same principal directions e 1 , e 2 , we can write:

B σ τ Φ = δ + Φ = δ + Φ = δ + Φ = δ + (C.5)

3

 3

Figure 2 :

 2 Figure 2: Geometrical parameters defining the collocation point x and the integration point ξ in the elliptical crack surface

  π on the boundary. Hence, the solution of the equation ∆χ =0 is solution for ϕ is:

Figure 4 :

 4 Figure 4: Parameters in polar coordinate system

Figure 5 :

 5 Figure 5: Pressure field around a single superconducting circular disc crack in an infinite matrix prescribed by a constant far-field pressure gradient parallel to x 1 direction.

3 Figure 6 :

 36 Figure 6: Darcy ellipsoidal inclusion Let us consider now the ellipsoidal inclusion having k * tends to infinity and d 3 tends to zero.At the boundary, we note v m the velocity in the matrix and v i the velocity in the inclusion. The velocity in the inclusion is constant for ellipsoidal inclusion. We note v i = ve 1 . At the boundary, the mass balance yields v i .n = v m .n, where n is the unit normal vector to the

Figure 7 :

 7 Figure 7: Dilute Mori-Tanaka and self-consistent estimates of effective permeability

  the expressions the functions f

9 ) 12 )

 912 When changing θ to θ+π, the functions γ, e, H remain unchanged, but g is changed to -g, α ε to α' ε = Arcsin(ρε -g) and the u ε to u' ε = tan (α' ε /2). It can be shown that: in the right-hand side must be considered as the Cauchy principal value and can be shown to take the following values:Replacing by (A.4) for γ, g and e in this expression and after integration we obtain:

  . Then the symmetries with respect to ϕ→-ϕ imply that we can write:

4 ) 2 B

 42 Appendix C: Linear relation between two tensors B and ФFor a general second order tensor B, we note B = B:δ δ δ δ and β =|B| (determinant of B). = λ + λ and β =λ 1 λ 2 . Therefore, λ 1 and λ 2 are the solutions of:

  and K(λ) and E(λ) the complete elliptical integrals of first or second kind.

  1,T BM -1 . ξ ɶ = 1 or

	. . ξ ξ B ɶ ɶ ɶ = 1 with:
	B ɶ	=		k B k							(64)
	The tensor B ∼ defines the geometry of the crack in the transformed problem. In order to calculate the half-diameters d ∼ 1 , d ∼ 2 and the principal directions e ∼ 1 , e ∼ 2 (unit vectors) of the
	transformed ellipse, we need to write:
	B ɶ	= λ ɶ	1	e ɶ	1	1 ⊗ + λ e ɶ ɶ	2	e ɶ	2	⊗	e ɶ	2	(65)

Appendix A: Integral calculation for a point inside the elliptical disc

Replacing by the [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF] in the main text in [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution based on singular integral equations[END_REF], the following equation (A.1) is found in which B is diagonal with eigenvalues λ 1 , λ 2 , and ξ.B.ξ = λ 1 ξ 1 2 + λ 2 ξ 2 2 . We use for ξ the polar coordinates (ρ,θ) in the local coordinate system having x as origin (see Figure 3). We have ξ 1 = x 1 +ρ cos θ, ξ 2 =x 2 +ρ sinθ. Thus, we have:

It is obvious that this integral must be considered as the Cauchy principal value:

where σ, τ are two scalars and

From this relation we deduce that:

So, φ 1 and φ 2 can be determined from the system of equations :

From these equations, we determine σ and τ function of φ 1 and φ 2 and so, finally, function of B and β:

( )

Finally, we define two scalar functions f and g as follows:


J P Stewart

J Douglas

C Di

Y Bozorgnia

N A Abrahamson

D M Boore

K W Campbell

E Delavaud

M Erdik

P J Stafford

Selection of a Global Set of GMPEs for the GEM-PEER Global GMPEs Project

Keywords: ground motion prediction equations, attenuation relations, Global Earthquake Model, seismic hazard assessment, tectonic regionalization

Ground-motion prediction equations (GMPEs) relate a ground-motion parameter (e.g. peak ground acceleration) to a set of explanatory variables describing the source, wave propagation path and site conditions. In the past five decades many hundreds of GMPEs for the prediction of PGA and linear elastic response spectral ordinates have been published. An accompanying paper discusses the pre-selection of GMPEs undertaken within the framework of the Global Earthquake Model (GEM) Global GMPEs project, coordinated by the Pacific Earthquake Engineering Research Center (PEER). Here, we discuss the following step undertaken to reduce the long list of pre-selected models down to a more manageable number for global hazard assessments. The procedure followed, consisting of an examination of the multidimensional (e.g. magnitude, distance and structural period) predicted ground-motion space in various ways and published quantitative tests of the GMPEs against observational data not used for their derivation, is discussed and illustrated for subduction zones.

INTRODUCTION

Ground-motion prediction equations (GMPEs) relate a ground-motion parameter (e.g. peak ground acceleration, PGA) to a set of explanatory variables describing the source, wave propagation path and site conditions. These independent variables invariably include magnitude, source-to-site distance and local site conditions, and often style-of-faulting (mechanism). There have been significant efforts recently to include additional parameters to model ground motions more realistically. In the past five decades many hundreds of GMPEs for the prediction of PGA and linear elastic response spectral ordinates have been published, which are summarized in a series of public reports by the second author (e.g., [START_REF] Douglas | Ground motion prediction equations 1964-2010[END_REF]. Therefore, the seismic hazard analyst is faced with the difficult task of deciding which GMPEs to use for a given project. This decision is a critical step in any hazard assessment because the resulting predicted spectra are strongly dependent on the GMPEs chosen. This paper discusses the selection of GMPEs undertaken within the framework of the Global Earthquake Model (GEM) Global GMPEs project, coordinated by the Pacific Earthquake Engineering Research Center (PEER). The overall GEM-PEER project is described by Di Alessandro et al. (2012). The first step of the GMPE selection process was made by pre-selecting from all the available models the most robust GMPEs. As described by [START_REF] Douglas | Compilation and critical review of GMPEs for the GEM-PEER global GMPEs project[END_REF], this Task 2 of the GEM-PEER project led to the choice of roughly ten GMPEs for each of three major tectonic regimes (shallow crustal earthquakes in active regions; subduction zones; stable continental regions, SCRs). For applications within GEM, the list of the pre-selected GMPEs in Task 2 is too long as the final selected GMPEs will be used by GEM for seismic hazard analysis of the entire world. Therefore, for practical reasons (e.g. calculation times) the selected number of GMPEs for each tectonic environment should be in general less than that used for local and regional hazard analyses. However, it is important that epistemic uncertainty in ground motion prediction is accounted for in the Global GMPEs project by selecting a range of GMPEs that cover the center, body and range of opinion.

It should also be noted that GMPE development is a continuously evolving research area, and new and/or updated GMPEs are developed as more empirical and simulated data become available and our knowledge of ground-motion hazard expands. Thus, the set of GMPEs proposed within Task 3 should not be viewed as a very long-term recommendation and it is subject to change.

In this paper, we describe the procedure adopted in Task 3 of the GEM-PEER project to select a relatively small set of about three recommended GMPEs for each major tectonic regime. The procedure is illustrated for the example of subduction-zone GMPEs.

SELECTION PROCESS AND FACTORS CONSIDERED

The GEM Task 3 core working group (the paper authors) met by web conference on two occasions to discuss the selection of GMPEs for subduction zones. In addition, discussions continued by email within the working group and occasionally with GMPE developers for clarification of certain aspects of their models. The GMPEs were selected from the candidate list produced in Task 2 (provided in [START_REF] Douglas | Compilation and critical review of GMPEs for the GEM-PEER global GMPEs project[END_REF], and the desired number of GMPEs selected was three. All members of the working group were present for at least one of the meetings, and all were included on the relevant correspondence including copies of slides, meeting minutes, and instructions for providing input to the working group chairs (John Douglas and Jonathan Stewart).

In the first meeting, the working group decided on the criteria to be considered in the selection of GMPEs from the list provided in Task 2. There was agreement that relevant criteria for consideration in GMPE selection include:

1. Giving more emphasis to GMPEs derived from international data sets than from local data sets. Exceptions can be made when a GMPE derived from a local data set has been checked internationally and found to perform well. 2. Giving more emphasis to GMPEs that have attributes of their functional form that we consider desirable, including saturation with magnitude, magnitude dependent distance scaling and anelastic attenuation terms. 3. If we have multiple GMPEs that are well constrained by data but exhibit different trends, it is desirable to capture those trends in the selected GMPEs to properly represent epistemic uncertainty.

In the selection process, we decided not to down-weight GMPEs with difficult to implement parameters (e.g., basin depth terms or depth to top of rupture), because those issues can be overcome with appropriate parameter selection protocols (e.g., [START_REF] Kaklamanos | Estimating unknown input parameters when implementing the NGA ground-motion prediction equations in engineering practice[END_REF]. We also decided not to down weight GMPEs that either lack site terms or whose modelling of site response is non-optimal (e.g., lack of nonlinearity) because GMPEs can be evaluated for a reference rock site condition in hazard analysis and site effects subsequently added in a hybrid process [START_REF] Cramer | Site specific seismic hazard analysis that is completely probabilistic[END_REF][START_REF] Goulet | Pitfalls of deterministic application of nonlinear site factors in probabilistic assessment of ground motions[END_REF].

To assist the working group members in making their selections according to the criteria above, two major compilations of information were prepared before the web conferences. First, so-called trellis plots were formulated that show spectral shapes for various magnitude and distance combinations, magnitude-scaling trends for different distance bins, distance-scaling trends for different magnitude bins, site effect terms, hypocentral depth-scaling terms, and standard deviation terms. Second, information from the literature was compiled on GMPE-data comparisons, giving emphasis to those studies that undertake formal analysis of residuals to provide insight into GMPE performance. These two data compilations are described further for the example of subduction zones in Sections 3 and 4 of the paper.

Individual members of the working group provided their recommendations for GMPE selection either orally as part of an open discussion and/or in written correspondence to the group facilitator (C. Di Alessandro). The working group chairs reviewed the input received, including their own, and made recommendations.

COMPARATIVE GMPE SCALING (TRELLIS PLOTS)

All of the GMPEs pre-selected in Task 2 [START_REF] Douglas | Compilation and critical review of GMPEs for the GEM-PEER global GMPEs project[END_REF] were programmed within Matlab. The predicted median ground motions and their aleatory variability from these implementations were checked against the graphs in the original references and against the results from previous implementations in other languages. As a standardized way of comparing the behavior of the GMPEs over the entire magnitude-distance (and other independent variables, e.g. site classification) range of interest to GEM, many trellis charts were drawn for the programmed GMPEs. These charts seek to display the multi-dimensional (magnitude, source-to-site distance, structural period etc.) predicted ground-motion space in various ways to understand the considered ground-motion models better. The aim is to help identify outliers with clearly nonphysical behavior but also to help guide the selection of models to capture the true epistemic uncertainty (e.g. the decay rate appears to be regionally dependent so it is important that this variation is captured).

The first type of these graphs (Figure 1) show predicted response spectra (color/line style coded for each GMPE) where each graph within the trellis has an x-axis of period and a y-axis of pseudospectral acceleration (PSA). The trellis has a super x-and y-axis of magnitude and distance, respectively and each graph within the trellis has its own axes with a common scale. This type of chart allows the experts to see how the spectrum predicted by each GMPE compares to the others over the magnitude-distance range of interest, e.g. are there any models that are consistently high or low or any with a different spectral shape? Because of the requirements imposed by the planned application of the GMPEs within GEM at the physical extremes of magnitude and source-to-site distance, the GMPEs were evaluated from the smallest magnitude considered of importance within the seismotectonic regime of interest (often M w 5) up to close to the largest magnitude that we feel possible in each of the different seismotectonic regimes and to the closest and farthest distance thought important on a global scale. Dotted lines are used for predictions for magnitudes and distances outside the limits of applicability stated by the GMPE developers or the range of data used for their derivation. However, since the goal of the GEM Global GMPEs project is to propose ground-motion models that work over all ranges of interest to GEM even the dotted lines were inspected by the experts. The idea is to thoroughly examine the models even outside their 'comfort zone' [START_REF] Bommer | On the selection of groundmotion prediction equations for seismic hazard analysis[END_REF]. The second type of graph plotted within trellis charts (Figure 2) are plots of predicted PSA against magnitude within a trellis chart with super x-and y-axes of period and source-to-site distance, respectively. This directly shows the magnitude scaling of ground motions. There are theoretical reasons why magnitude scaling is nonlinear and numerous observational studies have provided evidence for it, which is particularly important when magnitude is pushed very high (M w >8). The third set of trellis charts (Figure 3) are similar to the previous type but show scaling with distance for different magnitudes and periods. These plots show the decay rate for the various models, which can vary, for example, because of different anelastic attenuation representing variable crustal structures. Specifically for intraslab subduction earthquakes, the fourth type of chart shows the scaling with focal depth within a trellis chart of super axes of period and magnitude (not shown here). The final set of trellis charts (Figure 4) are similar to the first set (predicted median spectra) but show the predicted inter-event, intra-event and total aleatory variabilities (sigma) as a function of magnitude, distance and period.

Examining the trellis charts for the subduction-zone GMPEs shows that the [START_REF] Kanno | A new attenuation relation for strong ground motion in Japan based on recorded data[END_REF] model is an outlier, particularly at long periods, when evaluated for very large interface earthquakes (Figure 1) because linear magnitude-scaling is assumed (Figure 2). This suggests that this model is not a good candidate since this behavior may lead to erroneous hazard analyses for locations where very large events are possible. In addition, the sigmas associated with this model are higher than other models (Figure 4) suggesting that the functional form is too simple to model the behavior of subduction ground motions. Predictions from the [START_REF] Atkinson | Empirical ground-motion relations for subduction zone earthquakes and their application to Cascadia and other regions[END_REF] for interface events are typically a lower bound on estimates from the other considered GMPEs (Figure 1), except at long distances from very large earthquakes where the flat decay curve leads to high predicted PSAs (Figure 3). Predictions from the other GMPEs are more grouped particularly within the rough center of the distribution of available data from subduction zones (M w 6 to 7 and R from 50 to 150km) (Figure 1). The different GMPEs predict magnitude-dependent distance decay and nonlinear magnitude-scaling for M w >7.5 (Figures 2 and3).

GMPE-DATA COMPARISONS

In most cases, GMPEs are developed from the regression of strong motion data, so model-data comparisons are integral to the process by which they are prepared. Nonetheless, GMPE-data comparisons were considered a critical component of the selection process. GMPEs derived for SCRs are generally based on ground-motion simulations and, therefore, model-data comparisons are even more important for these equations. The value of these comparisons is often derived from the comparison data set being beyond the parameter space considered for the original GMPE. For example, the data may be derived from a different region from that used in the original model development, which can be useful for studies of model applicability to the data region and regional variations of ground motions generally. Another significant example specific to subduction zones is the recent availability of data sets from large-magnitude earthquakes (M w 8.8 Maule Chile and M w 9.0 Tohoku Japan) well beyond the upper bound magnitudes available during GMPE development.

Most GMPE-data comparisons in the literature consist of plots of ground motion intensity measures versus distance along with GMPE median trend curves. Plots of this type have limited applicability for formal analysis of GMPE performance because it can be difficult to judge trends when the data span a very wide range on the y-axis and because event-specific bias (event terms) are seldom taken into account. Accordingly, we generally restrict our literature compilation to studies that include a formal analysis of residuals into the GMPE-data comparisons.

The two general methods of residual analyses performed are: (1) the maximum likelihood approach of [START_REF] Scherbaum | On the use of response spectral reference data for the selection and ranking of ground motion models for seismic hazard analysis in regions of moderate seismicity: The case of rock motion[END_REF] and its extension to normalized intra-and inter-event residuals distributions by [START_REF] Stafford | An evaluation of the applicability of the NGA models to ground motion prediction in the Euro-Mediterranean region[END_REF], which are intended to judge the overall fit of model to data; and (2) analysis of intra-and inter-event residuals specifically targeted to investigations of GMPE scaling with respect to magnitude, distance, and site parameters [START_REF] Scasserra | Database for earthquake strong motion studies in Italy[END_REF]. Results of the latter approach applied to the Maule (Chile) and Tohoku (Japan) data [START_REF] Boroschek | Strong Ground Motion Attributes of the 2010 M w 8.8 Maule Chile Earthquake, Earthquake Spectra[END_REF][START_REF] Stewart | Implications of M w 9.0 Tohoku-oki Japan earthquake for ground motion scaling with source, path, and site parameters[END_REF] were particularly useful in identifying different distance attenuation trends from these large events. In the case of the Maule Chile data, the relatively slow distance attenuation of the [START_REF] Atkinson | Empirical ground-motion relations for subduction zone earthquakes and their application to Cascadia and other regions[END_REF] model provided a good fit to the data; whereas the Tohoku data attenuated relatively fast with distance and was better matched by the model of [START_REF] Zhao | Attenuation relations of strong ground motion in Japan using site classification based on predominant period[END_REF].

CONCLUSIONS

In this brief article we have presented the method used to select the final set of ground-motion models to be proposed by the GEM-PEER Global GMPEs Project, using the subduction-zone GMPEs as an illustrative example. The Task 3 core working croup (the authors of this paper) came up with consensus selections based on the types of criteria and information shown above and following intense discussion for the three main tectonic regimes (shallow crustal seismicity in active regions, stable continental regions and subduction zones). Our reasons for proposing each GMPE, and why others were not selected, were detailed in written documents and presentations, along with the material that led to our decisions. At the Global GMPEs plenary meeting in Istanbul on 17 th to 18 th May 2012, to which all experts of the project were invited, these arguments for each choice were carefully presented and the experts' feedback sought during the second day of the meeting, which focussed on this key step of the project. The wider Task 3 group consists of roughly 30 experts with worldwide experience of ground-motion modelling from dozens of countries. Based on the feedback from these experts final sets of GMPEs for the different regimes were defined and proposed to GEM, for use in their hazard assessments. The interested reader is referred to the final reports of the Global GMPEs Project for details of our recommendations.

Figure 1 .

 1 Figure 1. Trellis chart showing predicted PSAs for pre-selected subduction GMPEs for various interface earthquake scenarios for rock site conditions. Dashed lines indicated where the scenario falls outside the magnitude-distance range of validity of the model.

Figure 2 .

 2 Figure 2. Trellis chart showing magnitude-scaling of predicted PSAs for pre-selected subduction GMPEs for various structural periods and source-to-site distances for rock site conditions. Dashed lines indicated where the scenario falls outside the magnitude-distance range of validity of the model.

Figure 3 .

 3 Figure 3. Trellis chart showing decay of predicted PSAs for pre-selected subduction GMPEs for various structural periods and magnitudes for rock site conditions. Dashed lines indicated where the scenario falls outside the magnitude-distance range of validity of the model.

Figure 4 .

 4 Figure 4. Trellis chart showing inter-(between) and intra-(within) event and total natural log standard deviations of the pre-selected GMPEs for various interface earthquake scenarios.

ACKNOWLEDGMENTS

This study was funded by the GEM Foundation as part of the Pacific Earthquake Engineering Research Center's (PEER's) Global GMPEs project. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the sponsoring agencies.


S F Matar

R Pöttgen

A F Al Alam

N Ouaini

Ab initio investigations of the electronicstructure and chemical bonding of Li 2 ZrN 2

The electronicstructure of the ternary nitride Li 2 ZrN 2 is examined from abinitio with DFT computations for an assessment of the properties of chemicalbonding. The compound is found insulating with 1.8 eV band gap; it becomes metallic and less ionic upon removal of one equivalent of Li. The chemical interaction is found mainly between Zr and N on one hand and Li and N on the other hand. While all pair interactions are bonding, antibonding N-N interactions are found dominant at the top of the valence band of Li 2 ZrN 2 and they become less intense upon removal of Li. From energy differences the partial delithiation leading to Li 2-x ZrN 2 (x=∼1) is favored.

Introduction

In recent decades several classes of non-oxide lithium compounds were investigated for their use as negative electrode materials in lithium-ion batteries [START_REF]Lithium-Ion Batteries: Science and Technologie[END_REF] and [START_REF]High Energy Density Lithium Batteries[END_REF]. Besides silicides [START_REF] Huggins | Lithium alloys electrodes[END_REF] and [START_REF] Pöttgen | [END_REF], sulfides [5] and [6] and phosphides [6], lithium transition metal nitrides [6] and [7] were examined as promising anode material. Along with the many phase analytical and electrochemical studies of such materials, comparatively few investigations focused on the electronicstructure and chemicalbonding, e.g. for Li 15 Si 4 [8] or Li x FeS 2 [9].

In the course of our systematic studies of structure-property relationships and chemicalbonding of binary and ternary lithium-based silicides ( [START_REF] Pöttgen | [END_REF], [10], [11], [12], [13], [14], [15] and [16], and references therein) we recently investigated the electron precise silicides YLiSi (≡Y 3+ Li + Si 4-) [16] and La 2 Li 2 Si 3 (≡2La 3+ 2Li + 2Si1 2-Si2 4-) [15]. Such materials are potential candidates when searching for lithium mobility. Indeed, the calculations for YLiSi [16] indicated the possibility of an at least partial delithiation Li 1-x YSi, while keeping the yttrium-silicon substructure. Extending this work we were interested in the bonding patterns of related lithium nitridometalates Li x [T y N z ]. Out of the huge number of lithium nitridometalates [17], [18], [19], [20], [21] and [22] we choose the structure of Li 2 ZrN 2 [23] and [24]. This nitridometalate was first synthesized by Palisaar and Juza in 1971 from a reaction between lithium nitride Li 3 N and nitrogen excess ZrN 1.22 in an ammonia flow at high temperature [23]. The structure was determined by powder X-ray diffraction without further characterization of its physical properties. Later on, Niewa and Jacobs reinvestigated Li 2 ZrN 2 by neutron powder diffraction [24]. They fully confirmed the original structural model, however, with better precision, especially for the weakly scattering (X-rays) lithium atoms. Li 2 ZrN 2 has been studied with respect to the electrical conductivity, by solid state NMR spectroscopy [25], and coulometric titration [START_REF] Zheltonozhko | Proceedings of the Electrochemical Society 93-9 on Molten Molten Salt Chemistry and Technology[END_REF], indicating that this nitridometalate might be suitable for high-temperature lithium batteries. Keeping these favorable properties in mind we studied the electronicstructure and the properties of chemicalbonding and potential delithiation of Li 2 ZrN 2 within the quantum theoretical framework of the density functional theory (DFT) [START_REF] Hohenberg | [END_REF] and [28].

Computational methodology

Two computational methods within the DFT were used in a complementary manner. The Vienna abinitio simulation package (VASP) code [29] and [30] allows geometry optimization and total energy calculations. For this we use the projector augmented wave (PAW) method [30] and [31], built within the generalized gradient approximation (GGA) scheme following Perdew, Burke and Ernzerhof (PBE) [32]. The conjugate-gradient algorithm [START_REF] Press | Numerical Recipes[END_REF] is used in this computational scheme to relax the atoms. The tetrahedron method with Blöchl corrections [31] as well as a Methfessel-Paxton [START_REF] Methfessel | [END_REF] scheme were applied for both geometry relaxation and total energy calculations. Brillouinzone (BZ) integrals were approximated using the special k-point sampling of Monkhorst and Pack [35]. The optimization of the structural parameters was performed until the forces on the atoms were less than 0.02 eV/Å and all stress components less than 0.003 eV/Å 3 . The calculations are converged at an energy cut-off of 400 eV for the plane-wave basis set with respect to the k-point integration with a starting mesh of 4×4×4 up to 8×8×8 for best convergence and relaxation to zero strains.

Then all-electron calculations, equally based on the DFT with the GGA-PBE functional [32], were carried out for a full description of the electronicstructure and the properties of chemicalbonding. They were performed using the full potential scalar-relativistic augmented spherical wave (ASW) method [36] and [START_REF] Eyert | The Augmented Spherical Wave Method-A Comprehensive Treatment[END_REF]. In the ASW method, the wave function is expanded in atom-centered augmented spherical waves, which are Hankel functions and numerical solutions of Schrödinger's equation, respectively, outside and inside the so-called augmentation spheres. In the minimal ASW basis set, we chose the outermost shells to represent the valence states and the matrix elements were constructed using partial waves up to l max +1=3 for Zr and l max +1=2 for Li and N. Self-consistency was achieved when charge transfers and energy changes between two successive cycles were below 10 -8 and 10 -6 eV, respectively. The Brillouin zone integrations were performed using the linear tetrahedron method within the irreducible wedge [31]. The calculations are carried out assuming spin degenerate configuration. Besides the site projected density of states, we discuss qualitatively the pair interactions based on the overlap population analysis with the crystal orbital overlap population (COOP) [38]. In the plots, positive, negative, and zero COOP magnitudes indicate bonding, antibonding, and non-bonding interactions, respectively. We note that another scheme for describing the chemicalbonding, the ECOV (covalent bond energy) criterion based on both the overlap and the Hamiltonian populations is also accessible within the ASW method [START_REF] Eyert | The Augmented Spherical Wave Method-A Comprehensive Treatment[END_REF]. It provides similar qualitative results to the COOP.

Crystal structure

The trigonal structure (Table 1) of Li 2 ZrN 2 [23] and [24] is sketched in Fig. 1. Along the c-axis it consists of a succession of Zr-N-Li layers, with Li tetrahedrally coordinated to N while Zr is found in an octahedral nitrogen environment. The shortest interatomic distances occur for Li-N (2.098 Å) and Zr-N (2.253 Å), compatible with the sum of the covalent radii of 1.93 Å and 2.15 Å (r Li =1.23 Å, r Zr =1.45 Å and r N =0.70 Å). Thus significant bonding can be expected for Li-N and Zr-N in the ternary nitride. This is explained further in next sections.

Results of geometry optimizations, energies and electron localization

Geometry optimization was carried out for Li 2 ZrN 2 starting with the literature data [24]. Table 1 provides the experimental and calculated crystal structure characteristics. The geometry optimized lattice parameters and atomic positions are found in fairly good agreement with the experimental X-ray [23] and neutron powder determination data [24], especially for the internal coordinates of N and Li. In order to examine the effect of lithium de-intercalation, two additional calculations were done assuming a departure of one out of the two lithium atoms leading to the composition "LiZrN 2 " and a hypothetic lithium-free compound "ZrN 2 ". The results are given in 3rd and 4th columns of Table 1.

The full geometry relaxation keeps the trigonal structure within space group P3¯m1 with a reduction of the c/a ratio due to the departure of Li. This is especially visible for the fully de-lithiated compound resulting into "2□ Li ZrN 2 ", where □ Li designates a lithium vacancy, characterized by a largely reduced c/a ratio of 1.28. The internal atomic positions are also modified but the major changes occur for the cell volume which is only slightly reduced upon the departure of one lithium but largely reduced by ∼10 Å 3 for the departure of two lithium atoms and the reduction of the c/a ratio from 1.62 for "□ Li LiZrN 2 " down to 1.28 for "2□ Li ZrN 2 ". Also it can be noted that the departure of one lithium leads to a shorter Li-N distance and a reduction of the Zr-N separation. This translates a stronger overall bonding in "□ Li LiZrN 2 ". These results would indicate less perturbation of the structure with one lithium equivalent departed and large changes when a second lithium equivalent is removed. The calculated total energies can be used to assess the free energy of the process. For this purpose we further calculate lithium metal with the same potentials and method to obtain a total energy per atom of -1.909 eV. Using the total energy values in Table 1, the Gibbs free energy change in the deintercalation process of Li 2 ZrN 2 is:

For the removal of one lithium equivalent:

∆G 1 =E(Li 2 ZrN 2 )-E(Li)-E(□ Li LiZrN 2 )=-2.47 eV.
For the removal of a second lithium equivalent:

∆G 2 =E(□ Li LiZrN 2 )-E(Li)-E(2□ Li ZrN 2 )=-3.49 eV.
Both ∆G 1 and ∆G 2 values are larger than the ones obtained for the binary silicide Li 15 Si 4 (∆G=-1.136 eV [8]) and the ternary silicide YLiSi (∆G=-1.69 eV [16]). They correspond, with inverse sign, to the average de-intercalation voltage, V. It is interesting to note the 1.4 times larger magnitude for the removal of the second lithium equivalent, which is concomitant with a smaller Li-N distance, leading to a significant stabilization of the mono-delithiated compound. The de-intercalation voltage for the removal of one lithium equivalent, 2.47 V, is intermediate between the highest voltage magnitudes for the lithium intercalation compounds LiTO 2 (T=3d transition metal) [39] with values amounting to ∼3-4 V and Li x Si binaries such as Li 15 Si 4 with ∼0.30 V [8]. Then Li 2 ZrN 2 is found intermediate between oxide intercalation compounds and silicides [8] and [16].

In order to further assess these features, we analyze the charge density issued from the self consistent calculations using the AIM (atoms in molecules theory) approach [40] developed by Bader who devised an intuitive way of splitting molecules into atoms as based purely on the electronic charge density. Typically in chemical systems, the charge density reaches a minimum between atoms and this is a natural region to separate atoms from each other. Such an analysis can be useful when trends between similar compounds are examined; they do not constitute a tool for evaluating absolute ionizations. The results of computed charge changes (∆Q) are such that they lead to neutral compounds: Li 2 ZrN 2 : Li: +1; Zr: +3.28; N: -2.64, leading to an average charge transfer of ±5.28 electrons; "LiZrN 2 ": Li: +1; Zr: +3.33; N: -2.16, leading to an average charge transfer of ±4.34 electrons; and "ZrN 2 ": Zr +3.4; N: -1.7 with a charge transfer of ±3.4 electrons.

These values should be compared with the formal ones assuming a full departure/gain of electrons in the valence states of the atomic constituents: Li (2s 1 ); Zr (5s 2 , 4d 2 ) and N (2s 2 , 2p 3 ). With the three compositions considered, this leads to 6, 5 and 4 electrons transfer. Despite their larger magnitudes versus Bader charges, the values follow the same trend of decrease along the series of studied compositions. Then Li 2 ZrN 2 has a closed shell like configuration with all quantum states occupied as follows for the valence states: 2×Li (2s 0 ); 1×Zr (5s 0 , 4d 0 ) and 2×N(2s 2 , 2p 3+3 ), i.e. with a departure of 6 electrons from Li/Zr and an occupation of all the six 2p states of the two nitrogen atoms. This is opposite to "□ Li LiZrN 2 " and "2□ Li ZrN 2 ", respectively with one and two less electrons, which should present holes in the valence band with consequences on the electronic properties such as the metal/insulating character determined from the density of states (DOS)-vide infra. Thus the passage from Li 2 ZrN 2 to hypothetic "□ Li LiZrN 2 " and then to "2□ Li ZrN 2 " signals the decrease of the ionic character (removal of Li) which is consequently largest for "2□ Li ZrN 2 ". It needs to be mentioned here that "ZrN 2 " itself does not exist; only N excess ZrN is known (ZrN 1.22 ) [23]. This would indicate a destabilization of the compound if more than one lithium equivalent is removed, as the calculations pointed out. Consequently partial delithiation leading to Li 2-x ZrN 2 (x=∼1) is favored.

Finally the electrons transfer from Li to N can be illustrated by the electron localization function ELF [41] obtained from real space calculations: ELF=(1+χσ²) -1 . In this expression the ratio χσ=Dσ/Dσ 0 where Dσ=τσ-s-1/4( ρσ)²/ρσ and Dσ 0 =3/5(6π²) 2/3 ρσ 5/3 correspond, respectively, to a measure of Pauli repulsion (Dσ) of the actual system and to the free electron gas repulsion (Dσ 0 ) and τσ is the kinetic energy density. The ELF function is normalized between 0 (zero localization, blue areas) and 1 (strong localization, red areas) with the value of 1/2 corresponding to a free electron gas behavior (green areas). Fig. 2 shows the ELF contour plots for horizontal planes containing Zr (a), N (b), Li (c) and a vertical plane comprising Li and N. Whereas there is no electron localization around Zr (a) and Li (c, d) and green areas between them, there is a significant localization (orange-red areas) around N (b). This illustrates furthermore the charge density analysis above. The vertical plane along the c-axis (Fig. 2d) shows the layer-like structure of the compound with ∼0 and ∼1 ELF around Li and N, respectively.

Electronicstructure and chemical bonding

Using the calculated lattice parameters in Table 1 we analyze the electronicstructure and the chemicalbonding using all electron calculations with the full potential ASW method. The site projected density of states (PDOS) for Li 2 ZrN 2 , "□ Li LiZrN 2 " and "2□ Li ZrN 2 " are shown in Fig. 3.

From the presence of a band gap of 1.8 eV at the top of the valence band (VB) E V , the zero energy along the x-axis is considered with respect to E V and Li 2 ZrN 2 is a small gap insulator, in agreement with its light-green color [24]. This follows from the closed shell like electronic configuration discussed above along with Bader charges. The VB shows a large contribution from nitrogen due to twice as much N versus Zr in the structure and to the charge transfer from Li and Zr to N. Its s states are low lying at ∼-13 eV and the p states are found from -5 eV up to E V . Those states mix with itinerant Zr d states within the VB; Zr being an early 4d element, its d band is centered within the empty conduction band (CB). Li PDOS resemble the valence states of Zr and N and quantum mixing is expected between them as it is illustrated in next section of chemicalbonding. It is worth noting the dominating N p states at the top of the VB which differentiate the nitrides from the oxides whose VB top is dominated by the T states.

The removal of one equivalent of lithium leading to "□ Li LiZrN 2 " involves one electron less so that the top of the VB is rigidly shifted within the nitrogen states (creation of a hole) as it can be shown in Fig. 3b. The compound is now metallic and energy reference is with respect to the Fermi level (E F ). This feature is also observed for "2□ Li ZrN 2 " in Fig. 3c and a metallic like character is also observed. The decrease of the ionic character, as suggested from the Bader charges, is illustrated by a slightly broadened VB particularly for the s-like blocks, especially in the lithium-free hypothetic compound. The resemblance between the DOS in the three panels shows that lithium plays little role in the change of the structure which is mainly ensured by Zr and N binding. Similar changes in the electronicstructure occur for the various intercalated phases of ZrN versus Li x ZrN and ZrNCl versus Li x ZrNCl. For details we refer to a recent overview [42]. This is further exhibited in Fig. 4 showing the COOP for the three compounds Li 2 ZrN 2 (a) "□ Li LiZrN 2 " (b) and "2□ Li ZrN 2 " (c) accounting for atom multiplicities (Table 1). The whole range of the VB is of bonding character (positive COOP) for the main interaction, Zr-N as well as for Li-N. The Li-Zr bond is negligible as it can be expected from the electropositive nature of both and their large separation (∼2.8 Å). The N-N interaction is of half bonding/half antibonding character so that the large N PDOS at the top of the VB is due to antibonding (repulsive) N-N states. Interestingly, the N-N antibonding character is of less magnitude from Li 2 ZrN 2 (Fig. 4a) to "□ Li LiZrN 2 " (Fig. 4b) and "2□ Li LiZrN 2 " (Fig. 4c) as the Fermi level is moved down within the VB (cf. DOS section). This is also concomitant with the larger binding energy obtained from above energy differences calculated for Li 2 ZrN 2 and "□ Li LiZrN 2 ".

Fig. 1 .

 1 Fig. 1. The trigonal structure of Li 2 ZrN 2 . The edge-sharing LiN 4 tetrahedra and one ZrN 6 octahedron are emphasized.

Fig. 2 .

 2 Fig. 2. Electron localization function ELF plots for Li 2 ZrN 2 over 4 adjacent cells, showing (a) the basal plane at z=0 containing Zr, (b) the plane comprising N, (c) the plane containing Li, and (d) vertical plane with Li and N. Blue, green and red areas designate weak, free electron like and strong localizations respectively (see text). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Fig. 3 .

 3 Fig. 3. Site projected density of states (PDOS) for Li 2 ZrN 2 (a), □ Li LiZrN 2 (b) and 2□ Li ZrN 2 (c) for all atom multiplicities.

Fig. 4 .

 4 Fig. 4. Chemical bonding in Li 2 ZrN 2 (a), "□ Li LiZrN 2 " (b) and "2□ Li LiZrN 2 " (c) for all atom multiplicities.

Acknowledgments :

This work was financially supported by the Deutsche Forschungsgemeinschaft (PAK 177) and the Bundesministerium für Forschung und Technologie (LiVe-Lithium-Verbundstrukturen within the program LIB 2015). Support from the Centre Supérieur de la Recherche (CSR)-USEK-Lebanon is also gratefully acknowledged.