Laila Amir
email: l.amir@uca.ma

Michel Kern
email: michel.kern@inria.fr

Zoubida Mghazli
email: zoubida.mghazli@uit.ac.ma

Jean E Roberts

Intersecting fractures in porous media: mathematical and numerical analysis

Keywords: Porous media, fractured media, mixed finite element, domain decomposition

This paper studies a model for flow in a fractured porous medium with intersecting fractures. The fractures are treated as lower-dimensional manifolds and then physical transmission conditions express the pressure jump and the continuity of the flux across the fractures. Specific attention is borne to the conditions when several fractures intersect. The resulting system is discretized with mixed finite element, and the well-posedness of both the continuous and the discrete problems are proved. Then a domain decomposition method is formulated, so that the problem is reduced to the set of fractures, and a simple preconditioner is proposed. Numerical results exemplify the performance of the method.

Introduction

Porous media are characterized by their heterogeneity across multiple scales. Fractures are extreme examples of heterogeneity, as they are structures within the porous medium with one dimension (their width) much smaller that the other two-dimensions. However they play an important role in the description of the flow as their permeability may differ greatly from that of the surrounding medium. The permeability in the fractures may be much larger than that of surrounding medium so that the fractures act as hydraulic conductors, providing easy pathways for fluid flow. On the other hand, if the permeability in the fractures is much smaller than that of the surrounding medium, then the fractures act as barriers that prevent flow across them. In this paper we will be mostly interested in the former case, but the model presented is general.

The combination of the two characteristic features of fractures (small width and large influence on the flow) makes the simulation of flow and transport in fractured porous media quite challenging. Accordingly, this topic has been the subject of numerous studies, with several numerical methods being proposed over the past few years.

Originally fractures were included in the flow model by means of so-called multi-continua methods, such as double porosity [START_REF] Douglas | Dual porosity models for flow in naturally fractured reservoirs[END_REF][START_REF] Douglas | Derivation of the double porosity model of single phase flow via homogenization theory[END_REF], or stochastic models [START_REF] De Dreuzy | Efficient algorithms for the determination of the connected fracture network and the solution to the steady-state flow equation in fracture networks[END_REF][START_REF] De Dreuzy | Anomalous diffusion exponents in continuous two-dimensional multifractal media[END_REF]. Such techniques are still in current use, as seen for example in [START_REF] Kumar | Formal upscaling and numerical validation of unsaturated flow models in fractured porous media[END_REF]48,[START_REF] Noetinger | A quasi steady state method for solving transient darcy flow in complex 3d fractured networks accounting for matrix to fracture flow[END_REF][START_REF] Jerbi | A new estimation of equivalent matrix block sizes in fractured media with two-phase flow applications in dual porosity models[END_REF]. These methods made the study of large fracture networks feasible without the need for precise description of the geometry of individual fractures.

The development of advanced discretization techniques and meshing capabilities have enabled numerical methods that will describe the fracture network at the level of individual fractures. The size of the network, broadly characterized by the number of fractures that can be included in the model, is a defining feature of the works. We refer to the recent article [START_REF] Berre | Flow in fractured porous media: A review of conceptual models and discretization approaches[END_REF] for a survey of different possible approaches, and will only quote a subset of the works with which we are most familiar.

One can roughly distinguish between Discrete Fracture Networks (DFNs), where only the fractures are included, and mixed-dimensional models where the interaction between the threedimensional porous medium and the fractures is taken into account. The former class of models rests on the observation that fractures account for most of the flow, so that neglecting the flow in the porous medium may be justified. This simplification makes it possible for DNF models to accommodate a much larger number of fractures, as the resulting model is intrinsically twodimensional. [START_REF] Erhel | Flow simulation in threedimensional discrete fracture networks[END_REF][START_REF] Pichot | A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks[END_REF][START_REF] Pichot | A mortar BDD method for solving flow in stochastic discrete fracture networks[END_REF][START_REF] Formaggia | A reduced model for Darcy's problem in networks of fractures[END_REF][START_REF] Berrone | An optimization approach for large scale simulations of discrete fracture network flows[END_REF][START_REF] Hyman | dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport[END_REF][START_REF] Karimi-Fard | A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features[END_REF][START_REF] Fernando Benedetto | The virtual element method for underground flow simulations in fractured media[END_REF][START_REF] Fumagalli | Conforming, non-conforming and non-matching discretization couplings discrete fracture network simulations[END_REF][START_REF] Hédin | A hybrid high-order method for flow simulations in discrete fracture networks[END_REF] Following [START_REF] Alboin | Domain decomposition for flow in fractured porous media[END_REF][START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF][START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF][START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF], mixed-dimensional models treat the fractures as lower dimensional manifolds within the bulk medium. This usually results in a smaller number of fractures, as the whole three-dimensional volume has to be discretized. However, the simulation of such models remains challenging: one needs to choose between a conforming or a non-conforming mesh, and in the latter case, how to express the transmission conditions for pressure and flux across the fractures. A large number of numerical methods have been used [START_REF] Nordbotten | Modeling, structure and discretization of hierarchical mixed-dimensional partial differential equations[END_REF][START_REF] Nordbotten | Unified approach to discretization of flow in fractured porous media[END_REF][START_REF] Keilegavlen | Implementation of mixed-dimensional models for flow in fractured porous media[END_REF][START_REF] Ahmed | A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models[END_REF][START_REF] Sandve | An efficient multi-point flux approximation method for discrete fracture-matrix simulations[END_REF][START_REF] Frih | Modeling fractures as interfaces with nonmatching grids[END_REF], among which one finds finite volume methods [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF][START_REF] Ahmed | Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model[END_REF][START_REF] Faille | Model reduction and discretization using hybrid finite volumes for flow in porous media containing faults[END_REF], extended finite elements (XFEM) [START_REF] Del Pra | Well posedness of fully coupled fracture/bulk Darcy flow with XFEM[END_REF][START_REF] Flemisch | A review of the XFEM-based approximation of flow in fractured porous media[END_REF], and mixed finite elements [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF][START_REF] Amir | Décomposition de domaine pour un milieu poreux fracturé : un modèle en 3D avec fractures qui s'intersectent[END_REF][START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF][START_REF] Berrone | Flow simulations in porous media with immersed intersecting fractures[END_REF][START_REF] Arrarás | Mixeddimensional geometric multigrid methods for single-phase flow in fractured porous media[END_REF]. Models that specifically take into account the behavior of the flow at fracture intersection include [START_REF] Amir | Décomposition de domaine pour un milieu poreux fracturé : un modèle en 3D avec fractures qui s'intersectent[END_REF][START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF][START_REF] Berrone | Flow simulations in porous media with immersed intersecting fractures[END_REF][START_REF] Schwenck | Dimensionally reduced flow models in fractured porous media: crossings and boundaries[END_REF]. In most of the works cited above, the number of fractures in the domain remains limited, and is in any case much smaller than for DFNs. We note that a fairly general, and quite elegant, formulation for hierarchical problem has been proposed and analyzed in [START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF] (see also [START_REF] Nordbotten | Modeling, structure and discretization of hierarchical mixed-dimensional partial differential equations[END_REF][START_REF] Nordbotten | Unified approach to discretization of flow in fractured porous media[END_REF]). The models have been extended to deal with non-linear flow in the fractures (Forchheimer model), see [START_REF] Frih | Modeling fractures as interfaces: a model for Forchheimer fractures[END_REF][START_REF] Ahmed | A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models[END_REF] and multiphase flow, see [START_REF] Brenner | Hybriddimensional modelling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions[END_REF][START_REF] Gläser | A hybrid-dimensional discrete fracture model for non-isothermal two-phase flow in fractured porous media[END_REF][START_REF] Gläser | A discrete fracture model for two-phase flow in fractured porous media[END_REF][START_REF] Jerbi | A new estimation of equivalent matrix block sizes in fractured media with two-phase flow applications in dual porosity models[END_REF][START_REF] Fumagalli | Numerical modelling of multiphase subsurface flow in the presence of fractures[END_REF][START_REF] Volker Reichenberger | A mixeddimensional finite volume method for two-phase flow in fractured porous media[END_REF][START_REF] Hoteit | Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures[END_REF]. A posteriori error estimates were studied in [START_REF] Mghazli | Guaranteed a posteriori error estimates for a fractured porous medium[END_REF][START_REF] Hecht | A residual a posteriori error estimators for a model for flow in porous media with fractures[END_REF]. Domain decomposition (DD) methods appear naturally in the simulation of fractured media. The fractures define the subdomains between them, and the physical transmission conditions across the fractures give transmission conditions for the DD method. DD was first proposed as a possible formulation in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF][START_REF] Alboin | Domain decomposition for flow in fractured porous media[END_REF], was studied in more details in [START_REF] Flauraud | Domain decomposition for an asymptotic geological fault modeling[END_REF] and more recently in [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF][START_REF] Ahmed | A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models[END_REF]. The original mixed-dimensional problem is reduced to one posed only on the fractures, which is somewhat non-standard from the DD viewpoint, as it has to account both for the subdomain solutions, and also for the flow inside the fractures.

Finally, we note that a set of benchmarks for the simulation of flow in fractured media has been proposed in [START_REF] Berre | Verification benchmarks for single-phase flow in three-dimensional fractured porous media[END_REF]. This paper deals with a model of fractures in a porous medium, and looks at the case where fractures can intersect. The model is concerned with so-called large fractures, where the number of fractures remains small, but their length are significant in comparison with the size of the overall domain. The model is an extension of the one described in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF] (and most of the papers cited above), and was already used in [START_REF] Amir | Décomposition de domaine pour un milieu poreux fracturé : un modèle en 3D avec fractures qui s'intersectent[END_REF], as well as [START_REF] Berrone | Flow simulations in porous media with immersed intersecting fractures[END_REF] or [START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF] to which this work is closely related. While our model is maybe not as general as the one introduced and analyzed in [START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF], we felt it was still worthwhile to give a fairly direct analysis of the scheme for conforming grids, as the machinery presented in [START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF] is quite demanding.

The problem is discretized by means of mixed finite element. The present paper complements [START_REF] Amir | Décomposition de domaine pour un milieu poreux fracturé : un modèle en 3D avec fractures qui s'intersectent[END_REF] by giving the details for the well posedness of both the continuous and the discrete problem. In both cases, checking the hypotheses for the Babuska-Brezzi theory involves carefully constructing the solution to well-chosen auxiliary problems. Then a DD method is proposed along the lines of [START_REF] Ahmed | A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models[END_REF] or [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF] (the basic idea for the formulation goes back to [START_REF] Alboin | Domain decomposition for flow in fractured porous media[END_REF]). The reduction to an interface problem builds a Dirichlet to Neumann map. It was shown in [START_REF] Cowsar | Balancing domain decomposition for mixed finite elements[END_REF] how to build this map for the case of mixed finite element, by using the hybrid formulation. One would then follow the same pattern for the interface problem (this is what is done in [START_REF] Erhel | Flow simulation in threedimensional discrete fracture networks[END_REF] for DFNs), however the Dirichlet to Neumann map is global, and the mixed problem can no longer be reduced to a positive definite problem for the multipliers. For this reason, the interface problem is discretized by means of cell centered finite volumes, as in [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF].

A natural question when one deals with fractured media is whether to use a conforming or non-conforming mesh. In this work, we make the somewhat restrictive choice to use conforming meshes across the fractures. This is what enables a comparatively simple analysis in Sections 3 and 4. The numerical part of the work could be extended to non-matching grids along the lines of [START_REF] Frih | Modeling fractures as interfaces: a model for Forchheimer fractures[END_REF][START_REF] Frih | Modeling fractures as interfaces with nonmatching grids[END_REF][START_REF] Flemisch | A review of the XFEM-based approximation of flow in fractured porous media[END_REF][START_REF] Mghazli | Guaranteed a posteriori error estimates for a fractured porous medium[END_REF][START_REF] Hecht | A residual a posteriori error estimators for a model for flow in porous media with fractures[END_REF], and the analysis in [START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF] handles the non-conforming case by using mortars. The authors of [START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF] note that the choice of mortar variable (pressure trace of flux) influences the analysis. In that context, our work is closer to [START_REF] Frih | Modeling fractures as interfaces: a model for Forchheimer fractures[END_REF] or [START_REF] Hecht | A residual a posteriori error estimators for a model for flow in porous media with fractures[END_REF]. Indeed, the mortar variable is supported by the fractures, which already have a natural mesh. Additionally [START_REF] Hecht | A residual a posteriori error estimators for a model for flow in porous media with fractures[END_REF] shows that the error due the nonconformity can be controlled through a posteriori error estimation.

The outline for the rest of this paper is as follows: the physical problem and the underlying geometric setting are described in Section 2. In Section 3, existence and uniqueness for the continuous problem is proved, while Section 4 does the same for the approximate problem. A DD method is described in Section 5 and numerical examples are presented in Section 6.

Flow in a porous medium: problem with fractures

We consider the flow of an incompressible fluid inside a domain containing a network of fractures that possibly intersect. We denote by Ω a convex domain in R d (here d = 3), and by ∂Ω its boundary.

Fractures are thin structures in the subsurface that can act either as channels with a higher hydraulic conductivity than that of the porous matrix, or on the contrary as impermeable barriers that prevent flow. As noted in Section 1 we adopt the point of view of representing the fractures as lower dimensional structures. Among the references listed there, see in particular [START_REF] Alboin | Domain decomposition for flow in fractured porous media[END_REF][START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF][START_REF] Berrone | Flow simulations in porous media with immersed intersecting fractures[END_REF][START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF][START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF][START_REF] Nordbotten | Modeling, structure and discretization of hierarchical mixed-dimensional partial differential equations[END_REF][START_REF] Nordbotten | Unified approach to discretization of flow in fractured porous media[END_REF][START_REF] Schwenck | Dimensionally reduced flow models in fractured porous media: crossings and boundaries[END_REF].

We assume that Ω is decomposed into a set of subdomains Ω i , i = 1, . . . , N , such that the boundary of each subdomain Ω i is made up of the union of planar polygons, that may either be embedded in Ω or be part of the boundary ∂Ω. We denote by I = {1, 2, . . . , N } the set of (indices of) subdomains, and for each i ∈ I, we let

S i = {j ∈ I\{i} : Ω i ∩ Ω j = ∅}
be the set of its neighbors. Then,

Ω = i∈I Ω i , Ω i ∩ Ω j = ∅.
We define a fracture as a face of a subdomain Ω i that is embedded in Ω. For each fracture there exists exactly two subdomains, say Ω i and Ω j , that contain the fracture in their intersection, and the fracture can be denoted by γij :

γ ij = Ω i ∩ Ω j for i ∈ I, j ∈ S i .
We note that our definition of fracture is in line with that of Boon et al. [START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF], but is different from that of Berrone et al. [START_REF] Berrone | Flow simulations in porous media with immersed intersecting fractures[END_REF] where what we call a fracture is called a face, and a fracture is the union of several faces. As a consequence of defining fractures as delimiting subdomain boundaries, fractures cannot end "in the middle" of a subdomain, that is a fracture boundary is either part of the domain boundary, or is part of the boundary of another fracture. This assumption may appear restrictive, but it could easily be relaxed by introducing additional planes that support "inactive fractures", as is done in [START_REF] Berrone | Flow simulations in porous media with immersed intersecting fractures[END_REF]. One would then have to distinguish between active and inactive fractures, and this only involves a change in the transmission conditions across the fractures.

The set S of all fractures can be identified with S = {(i, j) : i ∈ I, j ∈ S i and i < j}. It will occasionally be useful to count each fracture twice, and this is the purpose of the set S = {(i, j) : i ∈ I, j ∈ S i }.

Fractures may intersect along edges belonging to two or more subdomains. If σ is such a segment, we denote by T σ = {(i, j) ∈ S, σ ⊂ γij } the set of all fractures (identified by the pair of indices) containing σ as part of their boundary. We denote by T the set of all intersections:

T :=    σ/ ∃R ⊂ S, σ = (i,j)∈R (∂γ ij )    Last, we have for all intersections σ ∈ T , γ ij = (γ ij ∩ Ω) σ for (i, j) ∈ T σ ,
and from now on, the fracture will be the open set γ ij . Finally we let γ = ∪ (i,j)∈ S γ ij be the set of all fractures.

For the example depicted on Figure 1, we have three subdomains, so that I = {1, 2, 3}, [START_REF] Ahmed | A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models[END_REF][START_REF] Ahmed | Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model[END_REF], (2, 3)}, while S = {(1, 2), (2, 1), (1, 3), (3, 1), (2, 3), (3, 2)}.

S 1 = {2, 3}, S 2 = {1, 3}, S 3 = {2, 1}, also S = {(1, 2),
There is only one intersection, σ = γ 12 ∩ γ 23 ∩ γ 31 , so that T σ = S. We note that, in order for σ to be an edge of Ω 1 , this subdomain should be considered as a polyhedron with 7 faces (γ 12 and γ 31 are considered two different fractures).

We may now write the flow equations in the subdomains and in the fractures, as well as the transmission conditions across the fractures and in the intersections. We do not repeat the In this work, we rely on the model in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF], which takes into account both conductive and impervious fractures.

For each subdomain Ω i , we find the usual flow equations (mass conservation and Darcy's law):

div( u i ) = f i on Ω i (1)
u i = -K i ∇p i on Ω i (2)
where p i u i , f i and K i are respectively the pressure, the Darcy velocity, the source term, and the permeability tensor.

The conditions on the exterior boundaries of Ω i can be of either Dirichlet or Neumann types:

p i = p di := p d|Γ Di in Γ Di (3)
u i • ν i = 0 in Γ N i , (4)
where

Γ Di = Γ i ∩ Γ D , Γ N i = Γ i ∩ Γ N and p d ∈ H 1/2 00 (Γ D
) is given. We suppose that K i is a diagonal tensor and that there exists 0 < K min and K max < ∞ such that, for all i ∈ I

0 < K min ξ 2 0 ≤ ξ T K i ξ ≤ K max ξ 2 0 , ∀ξ ∈ R d , ξ = 0. (5)
At the interfaces γ i,j , (i, j) ∈ S, instead of imposing the continuity of pressure as in a standard transmission problem, we impose as in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF] the condition

p i = p i,j + d i,j 4K n i,j ξ u i • ν i -(1 -ξ) u j • ν j in γ i,j , (6)
where p i,j is the pressure in the fracture γ i,j , K n i,j is the permeability of the fracture in the normal direction of that one, d i,j is the width of the fracture and ξ is a parameter that can be changed to obtain various alternatives. We suppose that d i,j is constant for all i and j and satisfies

0 < d m ≤ d i,j ≤ d M < +∞ (7)
Remark 1. Equation (6) takes into account a pressure jump across the fracture. It is obtained by integrating across the fracture in the direction normal to the fracture, cf [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF]. Note that, as the width d i,j of the fracture is quite small, if in addition the fracture is highly permeable (in the direction normal to the fracture) i. e. if K n i,j is sufficiently large, then the second term of the right-hand side of equation ( 6) is negligible, and equation (6) may be replaced by an equation expressing the continuity of the pressure across the fracture:

p i = p i,j = p j on γ i,j . (8)
This is the model which was studied in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF], and that we have used for the numerical examples in this article.

For the second transmission condition on the interfaces γ i,j , instead of continuity of flow, one has, as in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] and [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF], a system of equations that models the flow in the fractures. The discontinuity of the flux across a fracture yields a source term in the equations describing the flow in that fracture:

div γ i,j ( u i,j ) =f i,j + ( u i • ν i + u j • ν j ) in γ i,j (9)
u i,j = -d i,j K τ i,j ∇ γ i,j p i,j in γ i,j (10)
p i,j =p dγ i,j on ∂γ i,j ∩ ∂Ω D ( 11
)
u i,j • ν i,j =0 on ∂γ i,j ∩ ∂Ω N , (12)
where • K τ i,j is the tangential component of the permeability tensor of the fracture γ i,j . We assume that K τ satisfies a coercivity hypothesis similar to [START_REF] Alboin | Domain decomposition for flow in fractured porous media[END_REF]. [START_REF] Boffi | Mixed finite element methods and applications[END_REF] for the definition of that space) is the pressure given on its Dirichlet boundary,

• p dγ i,j ∈ H 1/2 (∂γ ij ∩ ∂Ω D ) (see
• and ν i,j is the outward-pointing, unit, normal vector field on ∂γ i,j .

The two operators div γ and ∇ γ are, respectively, the surface divergence and the surface gradient on γ.

Finally, on each intersection σ ∈ T , the continuity of pressure and the continuity of flux are imposed:

p i,j = p σ in σ, ∀(i, j) ∈ T σ ( 13
) (i,j)∈Tσ u i,j • ν i,j = 0 in σ, (14)
where p σ is the (unknown) pressure in σ.

The weak formulation of the problem and its well-posedness

Before giving the weak formulation of our problem, let us first recall some functional spaces. If

O = Ω is an open subset of R d (d = 2, 3) or O = γ is a (d -)-dimensional manifold in R d with = 1, 2, we define H 1 (O) := {v ∈ L 2 (O) / ∇ O v ∈ (L 2 (O)) n }, (n = d or ).
where ∇ O is the gradient with respect to the derivative in Ω or in γ. The space of traces on the boundary of O of functions in H 1 (O) is denoted by H 1/2 (∂O), endowed with the norm

φ 1 2 ,∂O = inf r∈H 1 (O) r| ∂O =φ r 1,O , ∀φ ∈ H 1/2 (∂O).
If δ is a part of ∂O, we denote

H 1/2 (δ) := {φ| δ / φ ∈ H 1/2 (∂O)}. (15)
Let δ D and δ N be disjoint parts of ∂O such that ∂O = δ D ∪ δ N and let E N,0 : H 1/2 (δ N ) -→ L 2 (∂O) be the extension operator defined by (see [START_REF] Gatica | A simple introduction to the mixed finite element method[END_REF])

E N,0 (η) := η on δ N 0 on δ D , ∀η ∈ H 1/2 (δ N ).
We define

H 1/2 00 (δ N ) := {v| δ N / v ∈ H 1 (O)
and v| δ D = 0}. There holds (see [START_REF] Gatica | A simple introduction to the mixed finite element method[END_REF])

H 1/2 00 (δ N ) = {η ∈ H 1/2 (δ N ) / E N,0 (η) ∈ H 1/2 (∂O)}, with the norm η 1 2 ,00,δ N := E N,0 (η) 1 2 ,
∂O . An important role will be played by the space

H(div , O) := { v ∈ L 2 (O) n / div O v ∈ L 2 (O)}, with n = d or n = d -,
where div O is the divergence operator with respect to the derivative on O. The normal traces of functions in H(div , O) belong to the space H -1/2 (∂O) defined as the dual space of H 1/2 (∂O), and we define H

The weak formulation

To give the weak formulation of the problem, we need to define two Hilbert spaces : W, a space for the velocity and M, a space for the pressure. Let

W = i∈I W i × (i,j)∈ S W ij , M γ = i∈I M i × (i,j)∈ S M ij M T = σ∈T H 1/2 00 (σ) (16)
and

M = M γ × M T
where for i ∈ I and (i, j) ∈ S

W i = { v i ∈ H(div , Ω i ) / v i • ν i ∈ L 2 (γ ij ), ∀j ∈ S i and v i • ν i | Γ N i = 0}, W ij = { v ij ∈ H(div , γ ij ) / v ij • ν i | ∂γ ij ∩∂Ω N = 0}, M i = L 2 (Ω i ), M ij = L 2 (γ ij ),
equipped with the norms for u ∈ W and r =

(r γ , µ T ) ∈ M u 2 W = i∈I u i 2 0,Ω i + div u i 2 0,Ω i + (i,j)∈ S u i,j 2 0,γ ij + div γ ij ( u i,j ) 2 0,γ ij + (i,j)∈ S u i • ν i + u j • ν j 2 0,γ ij , r 2 M = i∈I r i 2 0,Ω i + (i,j)∈ S r i,j 2 0,γ ij + σ∈T µ σ 2 00, 1 2 ,σ .
To obtain the mixed formulation we define two bilinear forms α ξ : W × W -→ R and β T :

W × M × M T -→ R by α ξ (u, v) = i∈I Ω i K -1 i u i • v i + (i,j)∈ S γ i,j (d i,j K τ i,j ) -1 u i,j • v i,j + (i,j)∈S γ i,j d i,j 4K n i,j (ξ u i • ν i -(1 -ξ) u j • ν j ) v i • ν i β T (v, r) = i∈I Ω i div v i r i + (i,j)∈ S γ i,j (div γ i, j v i,j -( v i • ν i + v j • ν j )) r i,j - σ∈T µ σ ,   (i,j)∈Tσ v ij • ν γ i,j   00, 1 2 ,σ
, and let the linear forms L f : M -→ R and L d : W -→ R be defined by

L f (r) = i∈I Ω i f i r i + (i,j)∈ S γ i,j f i,j r i,j L d (v) = - i∈I Γ Di p di v i • ν i - (i,j)∈ S p dγ ij v i,j • ν i,j ∂γ i,j ∩∂Ω D .
With these definitions of spaces and forms, the weak formulation of the problem given by equations ( 1)-( 4), [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF], and ( 9)-( 14) can be written as follows:

   Find u ∈ W, p = (p γ , λ T ) ∈ M such that α ξ (u, v) -β T (v, p) = L d (v), ∀v ∈ W, β T (u, r) = L f (r), ∀r ∈ M. (P)

Existence and uniqueness of the solution

We may now state and prove the main result of this section.

Theorem 3.1. Assume that the permeabilities K i on the subdomains Ω i , the permeabilities

K τ i,j
in the fractures γ i,j and the thickness d i,j are bounded above and below by a positive constants K min and K max , and assume also that ξ > 1 2 . Then problem (P) has a unique solution (u, p) ∈ W × M The proof of Theorem 3.1 uses the Ladyzhenskaya-Babuška-Brezzi (or LBB) theory for saddle point problems, and as usual consists of two main steps:

1. Prove the ellipticity of the bilinear form α ξ over the kernel of β T . This is done in Lemma 3.2 below;

2. Prove that β T satisfies the inf-sup condition, and this is done in Lemma 3.3. This will require the construction of an element of W, and each piece will be obtained as the solution of an appropriate auxiliary problem.

Let V be the space defined by

V = {v ∈ W : β T (v, r) = 0, ∀ r ∈ M}. Lemma 3.2. If ξ > 1 2 then the bilinear form α ξ (•, •) is V-elliptic. Proof. For v ∈ V, we have β T (v, r) = 0 for all r ∈ M.
By taking successively µ σ = 0, r i,j = 0, and r i = 0, for all σ and all i and j, we obtain in the distribution sense that

• div v i = 0 in L 2 (Ω i ), for i ∈ I, • div γ i,j v i,j = v i • ν i + v j • ν j in L 2 (γ ij ), for all (i, j) ∈ S, • (i,j)∈ S u i,j ν i,j = 0 in H -1/2 00 (σ), for all σ ∈ T . The norm of v ∈ V becomes v 2 W = i∈I v i 2 0,Ω i + (i,j)∈ S v i,j 2 0,γ ij + 2 v i • ν i + v j • ν j 2 0,γ ij and α ξ (v, v) = i∈I Ω i K -1 i v i • v i + (i,j)∈ S γ ij (d i,j K τ i,j ) -1 v i,j • v i,j + (i,j)∈ S γ i, j d i,j 4K n i,j ξ( v i • ν i ) 2 -2(1 -ξ)( v j • ν j )( v i • ν i ) + ξ( v j • ν j ) 2 .
By ( 5) and ( 7), the two first terms can be bounded from below. The integrand in the last term define a quadratic form whose eigenvalues are 1 and 2ξ -1.

We obtain

α ξ (v, v) ≥ K -1 max ( i∈I v i 2 0,Ω i + d m (i,j)∈ S v i,j 2
0,γ ij ) + d m 4 K -1 max min(1, 2ξ -1)
(i,j)∈S Proof. To show that β T satisfies the inf-sup condition, it suffices to show that there exists a positive constant C, such that given r 0 = (r 0 γ , µ 0 T ) ∈ M, there is an element v 0 ∈ W such that

( v i • ν i 2 0,γ ij + v j • ν j 2 0,γ ij ) ≥ K -1 max   min(1, d m )   i∈I v i 2 0,Ω i + (i,j)∈ S v i,j 2 0,γ ij   + d m 2 min(1, 2ξ -1)   (i,j)∈ S v i • ν i + v j • ν j 2 0,γ ij     ≥ C v 2 W where C = K -1 max min{min(1, d m ),
β T (v 0 , r 0 ) = r 0 2 M v 0 W ≤ C r 0 M .
For i ∈ I, the boundary of

Ω i is such that ∂Ω i = γ ik ∪ γ ji ∪ Γ Di ∪ Γ N i for k = i and j = i,
in case we have an intersection of fractures. We will use the following notations:

γ 1 i := γ ik and γ 2 i := γ ji . Given r 0 i in L 2 (Ω i ), for i ∈ I, let ϕ i the variational solution of the auxiliary problem        -∆ϕ i = r 0 i , in Ω i ∂ϕ i ∂ν ∂Ω i = 0, on γ 1 i ∪ γ 2 i ∪ Γ N i ϕ i = 0, on Γ Di .
By the Lax-Milgram lemma there exists a unique solution ϕ i ∈ H 1 0,Γ Di , and by Poincaré's inequality we have

|ϕ i | 1,Ω i := ∇ϕ i 0,Ω i ≤ C i r i 0,Ω i . (17)
Let

v 0 i = ∇ϕ i , then v 0 i ∈ L 2 (Ω i ) 3 , div v 0 i = r 0 i ∈ L 2 (Ω i ) and so v 0 i ∈ H(div , Ω i ) and moreover we have v 0 i • ν i | γ i = 0 for = 1, 2 and v 0 i • ν i | Γ N i = 0, which gives on the one hand that v 0 i • ν| γ i ∈ L 2 (γ i ) for = 1, 2
, and on the other hand v 0 i • ν i | Γ N i = 0, and so v 0 i is in W i . We now define the remaining components v 0 ij , (i, j) ∈ S so that for all σ ∈ T with (i, j) ∈ T σ div v 0 i,j = r 0 i,j and

(i,j)∈Tσ µ 0 σ , v 0 i,j • ν| σ = µ 0 σ 2 00,1/2,σ ,
where r 0 ij ∈ M ij and µ 0 σ ∈ H 1/2 00 (σ) are given as part of r 0 . This will be done by solving a mixed Dirichlet-Neumann problem, and we follow [36, Sec. 2.4.2], from which we recall the following results. For any σ ∈ T , H 1/2 00 (σ) is a Hilbert space with the inner product

ξ, η 00, 1 2 ,σ := E N,0 (ξ) , E N,0 (η) 1 2 ,σ , ∀ξ, η ∈ H 1/2 00 (σ)
where • , • 1 2 ,σ is the inner product in H 1/2 (σ). We consider the Riesz mapping R 00,σ :

H -1/2 00 (σ) -→ H 1/2 00 (σ) such that , η σ = R 00,σ ( ), η 00,1/2,σ , ∀ ∈ H -1/2 00 (σ), η ∈ H 1/2 00 (σ).
We fix an intersection σ ∈ T and consider for (i, j) ∈ T σ the fracture γ ij . Let ϕ i,j ∈ H 1 0,∂γ i,j \σ be the solution of

γ i,j ∇ γ ϕ i,j • ∇ γ ψ dx = γ i,j r 0 i,j ψ dx- 1 |T σ | R -1 00,σ (µ 0 σ ) ; ψ| σ σ , ∀ ψ ∈ H 1 0,∂γ i,j \σ (γ i,j ), (18)
where |T σ | is the number of elements of T σ . By the Lax-Milgram lemma, this problem has a unique solution and according to the respective continuous dependence result and the duality between H -1/2 00 (σ) and H

1/2 00 (σ) there exists constants C ij > 0 such that

|ϕ i,j | 1,γ ij ≤ C ij r 0 ij 0,γ ij + µ 0 σ 00, 1 2 ,σ . (19)
Now, let v 0 i,j = -∇ γ ϕ i,j . Then v 0 i,j ∈ H(div γ , γ i,j ) since div v 0 i,j = r 0 i,j , and we have

v 0 i,j • ν| σ = - 1 |T σ | R -1 00,σ (µ 0 σ ) and v 0 i,j • ν| ∂γ ij ∩Γ N = 0 (equivalently R 00,σ ( v 0 i,j • ν)| σ = 1 |T σ | µ σ ) so that v 0 i,j ∈ W ij .
We now have all the pieces in place to finish the proof.

β T (v 0 , r 0 ) = i∈I Ω i div v 0 i r 0 i dx + (i,j)∈ S γ ij (div γ v 0 i,j -( v 0 i • ν i + v 0 j • ν j )) r 0 i,j - σ∈T µ 0 σ ,   (i,j)∈Tσ v 0 ij • ν γ i,j   00, 1 2 ,σ = i∈I Ω i |r 0 i | 2 + (i,j)∈ S γ ij |r 0 i,j | 2 + σ∈T < µ 0 σ ; R -1 00,σ (µ 0 σ ) > 00, 1 2 ,σ = i∈I r 0 i 2 0,Ω i + (i,j)∈ S r 0 i,j 2 0,γ ij + σ∈T µ 0 σ 2 00, 1 2 ,σ = r 0 γ 2 Mγ + µ 0 T 2 M T = r 0 2 M . Also, v 0 2 W = i∈I v 0 i 2 0,Ω i + div v 0 i 2 0,Ω i + (i,j)∈ S ( v 0 i • ν i + v 0 j • ν j ) 2 0,γ ij + (i,j)∈ S v 0 i,j 2 0,γ ij + div γ ij ( v 0 i,j ) 2 0,γ ij = i∈I |ϕ i | 2 1,Ω i + r 0 i 2 0,Ω i + (i,j)∈ S |ϕ ij | 2 1,γ ij + r 0 i,j 2 0,γ ij .
By inequalities ( 17) and ( 19) we obtain

v 0 2 W ≤ (1 + 3C 2 0 ) i∈I r 0 i 2 0,Ω i + (i,j)∈ S(2 + C 2 ij ) r 0 i,j 2 0,γ ij +   (i,j)∈ S C 2 ij   σ∈T µ σ 2 00, 1 2 ,σ ≤ C r 0 γ 2 Mγ + µ 0 T 2 M T ≤ C r 0 2 M .
where C = max

   (1 + 3C 2 0 ), (2 + C 2 ij ),   (i,j)∈ S C 2 ij      . We deduce sup r∈M β T (v, r) v W ≥ β T (v 0 , r 0 ) v 0 W = r 0 2 M v 0 W ≥ 1 √ C r 0 M
which is the inf-sup condition.

As mentioned at the beginning of the Section, the proof of Theorem 3.1 is an immediate consequence of Lemmas 3.2 and 3.3.

To conclude this section, we briefly discuss a possible alternative formulation. The continuity on σ∈ T of the normal traces, condition [START_REF] Berrone | An optimization approach for large scale simulations of discrete fracture network flows[END_REF], was imposed in the equations through a Lagrange multiplier λ σ that represents a pressure on σ. As a minimum of regularity (H -1/2 00

) for these normal traces has been considered, the multiplier, which is in the dual space, must be more regular (in H 1/2 00 ). If we now want to impose less regularity on the multipliers λ σ , for example for reasons linked to the approximation as in section 4, we need to impose more regularity on the normal traces. We can suppose that the sum of the normal traces is in L 2 (σ), and introduce this condition in the functional space. We then replace the space W ij by

W ij = v ij ∈ H(div , γ ij ) / v ij • ν ij ∈ L 2 (σ), ∀σ ∈ T, v ij • ν i | ∂γ ij ∩∂Ω N = 0 (20)
and require the pressure to be only in L 2 (σ)

M T := σ∈T L 2 (σ). (21)
In general, we denote with a hat the objects associated with the "new" spaces. In the definition of β T (•, •), the duality bracket becomes an integral:

µ σ ;   (i,j)∈Tσ v ij • ν γ i,j   00, 1 2 ,σ = σ µ σ   (i,j)∈Tσ v ij • ν γ i,j   ds.
We then define the corresponding weak form of the problem

     Find u ∈ W, p = (p γ , λ T ) ∈ M such that α ξ (u, v) -β T (v, p) = L d (v), ∀v ∈ W, β T (u, r) = L f (r) ∀r ∈ M. ( P)
The proof of the existence theorem 3.1 goes through with the L 2 -norm instead of H

Approximation by mixed finite element

In this section we present the discrete setting for approximating problem ( P).

We start with meshes on the fractures, and we consider meshes on the subdomains that are conforming across the fractures, that is, given subdomains Ω i and Ω j with a common fracture γ ij , the trace of the meshes of Ω i and Ω j on γ ij coincide, and they coincide with the mesh on γ ij . Additionally, the meshes of the fractures are assumed to match along the intersections. We denote the mesh on subdomain Ω i by T ih and by E ih the set of faces of T ih . Similarly, the mesh on a fracture γ ij is denoted by T ijh and the set of edges of T ijh is E ijh . Each intersection σ inherits a mesh T σh from the fractures that contain it.

Problem P will be approximated by mixed finite elements. We will use a combination of the lowest order Raviart-Thomas-Nédélec spaces in the subdomains and lowest order Raviart-Thomas spaces in the fractures. For numerical purposes, it will be more convenient to use the formulation given in Proposition 3.4, so that we only work with functions, and not elements of dual spaces. The approximation spaces M h,i , W h,i , M h,i,j and W h,i,j are subspaces respectively of M i , W i , M ij and W ij such that (M h,i , W h,i ) is the lowest order Raviart-Thomas-Nédélec space on Ω i , and (M h,i,j , W h,i,j ) is the lowest order Raviart-Thomas space on γ i,j (see [START_REF] Gatica | A simple introduction to the mixed finite element method[END_REF]Chap. 3] or [START_REF] Boffi | Mixed finite element methods and applications[END_REF]Sec. 2.3] for the precise definitions). The subspace M h,σ of L 2 (σ) is composed of functions that are constant on each cell of the mesh. The approximation spaces W h and M h , are the finite dimensional subspaces of W and M defined by

W h = i∈I W h,i × (i,j)∈ S W h,ij ⊂ W, M h = i∈I M h,i × (i,j)∈ S M h,ij × σ∈T M h,σ ⊂ M.
The global discrete problem is written as:

   Find u h ∈ W h , p h ∈ M h such that α ξ (u h , v h ) -β T (v, p h ) = L d (v h ), ∀ v h ∈ W h β T (u h , r h ) = L f (r h ), ∀ r h ∈ M h . (P h )
In the rest of this section, we prove that the discrete problem (P h ) is well posed, and we prove that the discrete solution converges towards the continuous one. Similarly to Section 3, we first state the results, prove two lemmas for the two technical results needed, and then put everything together.

Theorem 4.1. Under the hypotheses of Proposition 3.4, the discrete problem (P h ) admits a unique solution (u h , p h ) ∈ W h × M h , and the following convergence result holds:

lim h→0 u -u h W + p -p h M = 0.
Remark 2. Assuming more regularity for the continuous solution would give an optimal convergence rate (cf [START_REF] Gatica | A simple introduction to the mixed finite element method[END_REF]Sec. 4.2] or [START_REF] Boffi | Mixed finite element methods and applications[END_REF]Prop. 7.1.2].

To prove Theorem 4.1, it suffices to show the analogs of Lemma 3.2 and 3.3, and to show that the constants in the inequalities can be chosen to be independent of h.

First we show that α

ξ : W h × W h -→ R is uniformly V h -elliptic, where V h = {v h ∈ W h : β T (v h , r h ) = 0, ∀r h ∈ M h }, (22)
and that the ellipticity constant can be chosen independently of h. Generally V h ⊂ V, but if the inclusion holds, then uniform ellipticity of α ξ (•, •) follows. We can prove this inclusion by building an operator P h = (P γh , P T h ) :

M γ × M T -→ M γh × M T h
such that, for all r = (r γ , µ T ) ∈ M there holds (cf [START_REF] Fortin | An analysis of the convergence of mixed finite element methods[END_REF]Lemma 4.2], or [START_REF] Boffi | Mixed finite element methods and applications[END_REF]Prop. 5.4.2]

β T (v h , r -P h r) = 0, ∀v h ∈ W h . ( 23
)
Lemma 4.2. The bilinear form α ξ :

W h × W h -→ R is uniformly V h -elliptic.
Proof. We take for P γh the projection operator P γh = (P ih ) i∈I , (P ijh ) (i,j)∈ S : M γ -→ M γh defined by the following degrees of freedom

for i ∈ I, K (r i -P ih (r i )) dx K = 0, ∀K ∈ T ih ( 24
)
for (i, j) ∈ S,

E (r ij -P ijh (r ij )) dx E = 0, ∀E ∈ T ijh . (25)
Similarly, the projection operator

P T h = (P σh ) σ∈T : M T -→ M T h is defined by for σ ∈ T, e (µ σ -P σh (µ σ )) ds e = 0, ∀e ∈ T σh . (26)
Let r = (r γ , µ T ) be given in M. We have, for all

v h ∈ W h β T (v h , r -P h r) = i∈I Ω i div v i (r i -P ih (r i )) + (i,j)∈ S γ i, j (div γ i, j v i, j -( v i • ν i + v j • ν j )) (r i, j -P ijh (r i, j )) - σ∈T σ (µ σ -P σh (µ σ ))   (i,j)∈Tσ v ij • ν γ i,j   = i∈I K∈T ih K div v i (r i -P ih (r i )) + (i,j)∈ S E∈T ijh E (div γ i, j v i, j -( v i • ν i + v j • ν j )) (r i, j -P ijh (r i, j )) - e∈T σh e (µ σ -P σh (µ σ ))   (i,j)∈Tσ v ij • ν γ i,j   .
Now, for v ∈ W h , we know that the Raviart-Thomas-Nédélec elements are such that div v i is a constant in each K ∈ T ih and v i • ν i is a constant in each E ∈ E ih , in particular on the elements of T ijh for some j; similarly div γ i, j v i, j is a constant in each E ∈ T ijh and v ij • ν γ i,j is a constant in each e ∈ E ijh , in particular on the elements of T T h . Then ( 23) is an immediate consequence of ( 24), ( 25) and ( 26), and we obtain the uniform V-ellipticity of a ξ (•, •).

To prove that the discrete inf-sup condition holds uniformly we use a characterization analogous to the last property.

Lemma 4.3. The discrete inf-sup condition holds on W h × M h with a constant independent of h.

Proof. It is known (see [START_REF] Fortin | An analysis of the convergence of mixed finite element methods[END_REF]Prop. 4.2], or [START_REF] Boffi | Mixed finite element methods and applications[END_REF]Prop. 5.4.2]) that in order to prove the Lemma, it is sufficient to show that there exists an operator Π h such that, for all v ∈ W

β T (v -Π h (v), r) = 0, ∀r ∈ M h (27)
and a constant C > 0 independent of h such that

Π h (v) W ≤ C v W , ∀v ∈ W. (28)
Additional regularity assumptions are required in order to define the interpolation operator associated with W. If s > 2 and we denote

W s := W ∩ i∈I L s (Ω i ) × (i,j)∈ S L s (γ ij )
, then one may define (see [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF][START_REF] Boffi | Mixed finite element methods and applications[END_REF])

Π h := (Π ih ) i∈I , (Π ijh ) (i,j)∈ S : W s -→ W h .
by the following degrees of freedom

for i ∈ I, E ( v i -Π ih ( v i )) • ν E ds E = 0, ∀E ∈ E ih (29)
for (i, j) ∈ S,

e ( v ij -Π ijh ( v ij )) • ν e ds e = 0, ∀e ∈ E ijh . (30)
and the stability condition (28) holds for v ∈ W s . Let v ∈ W s be given. For all r ∈ M h , we have

β T (v -Π h (v), r) = i∈I Ω i div ( v i -Π ih ( v i )) r i + (i,j)∈ S γ i,j div γ i, j ( v i,j -Π ijh ( v ij ))r i, j - (i,j)∈ S γ i, j [( v i -Π ih ( v i )) • ν i + ( v j -Π jh ( v j )) • ν j ] r i, j - σ∈T σ µ σ   (i,j)∈Tσ ( v ij -Π ijh ( v ij )) • ν γ i,j   .
By decomposing the integrals over the meshes, we obtain

β T (v -Π h (v), r) = i∈I K∈T ih K div ( v i -Π ih ( v i ) r i + (i,j)∈ S E∈T ijh E div γ i, j ( v i,j -Π ijh ( v ij ))r i, j - (i,j)∈ S E∈T ijh E [( v i -Π ih ( v i )) • ν i + ( v j -Π jh ( v j )) • ν j ] r i, j - σ∈T e∈T σh e µ σ   (i,j)∈Tσ ( v ij -Π ijh ( v ij )) • ν γ i,j   .
Notice that the last two terms in the expression of β T (v -Π h (v), r) vanish since r ij is constant in E and µ σ is constant over e. On the other hand, since r i is constant in each K ∈ T ih an integration by part and (29) give

K div ( v i -Π ih ( v i )) r i = ∂K ( v i -Π ih ( v i )) • ν ij r i = E⊂∂K E ( v i -Π ih ( v i )) • ν ij r i = 0.
Similarly, since r i, j is constant in each E ∈ T ijh we have by integration by part and ( 30)

E div γ i, j ( v i,j -Π ijh ( v ij ))r i, j = e∈∂E e ( v i,j -Π ijh ( v ij )) • ν ij r i, j = 0.
This shows that [START_REF] Faille | Model reduction and discretization using hybrid finite volumes for flow in porous media containing faults[END_REF] holds and together with (28) establishes the discrete inf-sup condition.

It is now standard to deduce Theorem 4.1 from the two lemmas above, and obtain the convergence result. Note that the operator Π h in Lemma 4.3 is only defined for functions in the smaller space W s but convergence holds with only the hypotheses needed for the existence result in Proposition 3.4 by a density argument (see for instance [START_REF] Gatica | A simple introduction to the mixed finite element method[END_REF]Lemma 4.5]).

Domain decomposition and reduction to an interface problem

In this section we use domain decomposition techniques to eliminate all volume unknowns and to reduce the problem to one set on the fractures.

We will only consider a simplified setup, with only one intersection, so that T = σ. Also, we consider only the model of very permeable fractures with continuity of the pressure, and thus replace equation ( 6) by ( 8). This hypothesis had already been used in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF] for introducing a domain decomposition method in a medium with a single fracture. Formally, it corresponds to the limit K n i,j → ∞. We also assume that all fractures have only Dirichlet boundary conditions A Neumann-Neumann domain decomposition method for mixed finite element has been presented and analyzed in [START_REF] Cowsar | Balancing domain decomposition for mixed finite elements[END_REF] (see also [START_REF] Pencheva | Balancing domain decomposition for mortar mixed finite element methods[END_REF]). In that work it was shown that using the mixed hybrid formulation led to a positive definite Dirichlet to Neumann map (or Schur complement at the discrete level). We follow this approach to eliminate the volume unknowns in the subdomains and reduce the problem to the fractures. In order to obtain a positive definite system, we then solve the resulting problem by a cell-centered finite volume scheme as already hinted at in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF] and detailed in [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF]. This is possible because the pressure traces on the boundary for the problem in the subdomains coincide with the fracture cell pressure unknowns.

To define the Dirichlet to Neumann map, we use linearity to decompose, for each subdomain Ω i , i ∈ I, the solution of problem ( 1), ( 2) with boundary conditions (3), ( 4) and ( 8) as u i = u 0 i + u * i and p i = p 0 i + p * i , with ( u 0 i , p 0 i ) and ( u * i , p * i ) the respective solutions of div

u 0 i = 0 on Ω i u 0 i = -K i ∇p 0 i on Ω i p 0 i = 0 in Γ Di u 0 i • ν i = 0 in Γ N i p 0 i = p i,j in γ i,j . (P 0 i ) div u * i = f i on Ω i u * i = -K i ∇p * i on Ω i p * i = p di in Γ Di u * i • ν i = 0 in Γ N i p * i = 0 in γ i,j . (P * i )
and we let χ i = ( u * i • ν i ) |γ i . We now define the Dirichlet to Neumann map (aka Steklov-Poincaré operator) by

S i : j∈S i M hij → j∈S i M hij λ i := (λ ij ) j∈S i → ( u 0 i • ν ij ) |γ ij j∈S i
where u 0 i is the solution of problem (P 0 i ) with boundary data p ij = λ ij . Then for each fracture γ ij , (i, j) ∈ S, we can rewrite equation [START_REF] Arrarás | Mixeddimensional geometric multigrid methods for single-phase flow in fractured porous media[END_REF] as

S i (λ ij ) + S j (λ ij ) -div γ i,j u i,j = f i,j + χ i + χ j on γ i,j (31)
together with Darcy's law [START_REF] Fernando Benedetto | The virtual element method for underground flow simulations in fractured media[END_REF], the boundary conditions

λ ij = p dγ in ∂γ ∩ ∂Ω λ i,j = λ σ in σ, i ∈ I, j ∈ S i
and condition [START_REF] Berrone | An optimization approach for large scale simulations of discrete fracture network flows[END_REF] over the intersection σ.

We define the global Steklov-Poincaré operator by extending λ i by 0 to γ\γ i , and we do the same for χ i :

S(λ) = I S i (λ) and χ = I χ i .
We note that because the Steklov-Poincaré operators are global, it is not possible to use a hybrid formulation to reduce the problem to the set of edges as was done in [START_REF] Erhel | Flow simulation in threedimensional discrete fracture networks[END_REF][START_REF] Pichot | A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks[END_REF]. Indeed, one can no longer eliminate locally the cell unknowns in Darcy's law. This may be linked to the fact, noticed in [START_REF] Hoteit | Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures[END_REF], that for a discretization of the coupled model by mixed-hybrid finite element, both the face and edge unknowns must be kept for the fractures, in addition to the cell pressure unknowns. In the domain decomposition context, this was also noticed in [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF], where finite volumes are used to discretize the interface problem, and we follow that choice here. For this reason, we return to the primal formulation of the problem over the fractures. Alternatively, one could keep the pressure-velocity formulation with mixed finite element, as is done in [START_REF] Ahmed | A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models[END_REF], at the cost of losing the positive definite character of the problem, and then the resulting system can be solved by GMRES. We write the system over the global interface:

S(λ) -div γ d γ K τ γ ∇ γ λ = f γ + χ in γ λ = p dγ in ∂γ ∩ ∂Ω λ i,j -λ σ = 0 in σ, (i, j) ∈ S (i,j)∈ S -d i,j K τ i,j ∇ γ λ i,j |σ • ν i,j = 0 in σ, (P * γ )
where d γ , K τ γ and f γ are defined in an obvious way. System (P * γ ) is equivalent to the original problem (P), but involves only the unknowns on the fractures. In comparison with domain decomposition problem, system (P * γ ) features the additional term -div γ d γ K τ γ ∇ γ λ that takes into account the flow in the fractures. The flow in the intersection is represented by the last equation in (P * γ ). We discretize the reduced system by a cell-centered finite volume method (see [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF] or [START_REF] Eymard | Finite volume methods[END_REF]). For each fracture γ ij , (i, j) ∈ S, we write the balance equation for each cell K ∈ T ijh

S(λ) |K + L∈N (K) τ KL (λ K -λ L ) + e⊂K e∈E ijh ∩∂γ ij τ K,e λ K -p dγ i,j + e∈T σh τ K,e (λ K -p σ,e ) = f K , (32)
where τ KL and τ K,e are the transmissivities across the edge between cells K and L and at the boundary, and where it is understood that the second or third sums may be empty for interior cells. We emphasize that there is only one unknown per segment on the intersection, expressing the pressure continuity in equation [START_REF] Berrone | Flow simulations in porous media with immersed intersecting fractures[END_REF]. The finite volume unknown at the center of each grid cell coincides with the mixed finite element pressure unknown, which makes the coupling between the two methods particularly convenient.

We also write an equation for the continuity of the flux across the intersection σ:

K⊃e τ K,σ (λ K -λ σ,e ) = 0, ∀e ∈ T σh (33)
where the inner sum is over all the fractures that contain the intersection.

The system formed by equations ( 32) and ( 33) is a square system where the unknowns are one pressure unknown for each cell, and one for each edge in the intersection. The resulting matrix is symmetric and positive definite. Indeed, the Steklov-Poincaré operator is positive semi-definite, and this is also the case for the part coming from the finite volume, as one can see by the usual discrete integration by part technique. Last, the Dirichlet boundary conditions make the whole system definite.

The system can be solved by a conjugate gradient method. At each iteration, evaluation of the operator on a vector λ requires the solution of a Dirichlet problem in each subdomain, given the values on the fractures. The pressure unknowns on the fracture are the values of the pressure traces at the center of the faces of the volume mesh for each subdomain.

A preconditioner is required to accelerate the convergence of the conjugate gradient algorithm. The Neumann-Neumann preconditioner and its two-level variants (see [START_REF] Tallec | Domain decomposition methods in computational mechanics[END_REF][START_REF] Toselli | Domain decomposition methods -algorithms and theory[END_REF][START_REF] Tarek | Domain decomposition methods for the numerical solution of partial differential equations[END_REF]) have shown their efficiency for non-overlapping domain decomposition methods. However in this case, it is likely that a simpler preconditionner wil lbe be more effective: if we compare the differential order of the terms in (P * γ ), we see that the divergence term is second order, while the Steklov-Poincaré operator are first order operators1 . So, unless the value of d K τ γ is very small, we should expect that a preconditioner based on inverting the differential part of the operator should be efficient. This is confirmed experimentally in Section 6.2.

6 Numerical results

Description of the test cases

In this part we present some numerical results to illustrate the properties of the domain decomposition method presented in Section 5. We consider a set of tests that highlight the effect of the fractures on the performance of the domain decomposition algorithm. In all the numerical experiments, we consider a domain Ω represented in Figure 2 that contains in addition to the 3 fractures of the problem model (Figure 1) a single isolated fracture. Three tests cases are considered. In all three cases, we impose a pressure equal to 1 on top of the fractures and a pressure equal to 0 on the bottom of the fractures. The permeability of the bulk is taken equal to 1, while that of the fractures is 100 times as large, to illustrate the effect of permeable fractures.

Horizontal flow

In the first experiment, a pressure gradient is imposed form left to right (p = 1 on the left boundary, p = 0 on the right, and no flow on the other boundaries). We show the influence of the fractures by comparing the results obtained without any fractures, and then with the fractures as in Figure 2. In the first case, the domain decomposition method becomes the usual one. Then boundary conditions on the fractures, are Dirichlet conditions, with p = 1 at the top and p = 0 at the bottom.

The permeability in the fracture multiplied by its thickness is equal to 100, so that the fluid tends to flow rapidly along the fracture.

Figure 3 on the left represents the pressure obtained when no fractures are present (we represented each subdomains separately). The flow is only horizontal, the pressure varies linearly from left to right. If we take into account the fractures, the situation change, as shown on Figure 3

Vertical flow

For the third test we impose a vertical gradient of pressure. Figure 5 represents the geometry and the boundary conditions for this test case: we impose a vertical gradient of the pressure, p = 1 on the top boundary, p = 0 on the bottom boundary and a zero flux condition on the other boundaries.

Figure 6 compares the pressure fields obtained without the fractures (on the left), and with the fracture (on the right). One sees a linear gradient from top to bottom on the left figure,

Performance of the preconditioner

We briefly discuss the performance of the preconditioner introduced at the end of Section 5.

The left part of Figure 7 compares the convergence of the conjugate gradient without and with the preconditioner. As can be seen, the simple preconditioning by the flow in the fracture alone is quite efficient. The right part of Figure 7 shows how the behavior of the preconditioner varies with the permeability of the fracture. As expected from the motivations for the preconditioner given in Section 5, the preconditioner becomes more effective for larger values of the fracture permeability.

Conclusion

In this paper we presented a model for the flow in a fractured porous medium in 3 dimensions that takes into account the intersecting fractures. In this model, we assimilate the fractures to interfaces through which the pressure is continuous but flux is discontinuous. We proved the well-posedness of the model for a general distribution of the fractures, as well as the convergence of a mixed finite element approximation. We used a domain decomposition method to reduce the problem to an abstract equation posed only over the fractures. The use of an iterative method to solve the system on the global interface requires an appropriate preconditioner. Possible directions for extending this work include modeling transport, and handling nonconforming meshes, either using mortar elements [START_REF] Pichot | A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks[END_REF][START_REF] Pichot | A mortar BDD method for solving flow in stochastic discrete fracture networks[END_REF][START_REF] Boon | Robust discretization of flow in fractured porous media[END_REF][START_REF] Pencheva | Balancing domain decomposition for mortar mixed finite element methods[END_REF] or controlling the nonconformity by using a posteriori error estimates as in [START_REF] Hecht | A residual a posteriori error estimators for a model for flow in porous media with fractures[END_REF][START_REF] Mghazli | Guaranteed a posteriori error estimates for a fractured porous medium[END_REF].

31 Figure 1 :

 311 Figure 1: Geometry of domain

-1/ 2

 2 00(Γ N ) as the dual of H 1/2 00 (Γ N ). We denote by < • , • > 00, 1 2 ,Γ N the duality between H 1/2 00 (Γ N ) and H -1/2 00 (Γ N ).

d m 2 min( 1 ,Lemma 3 . 3 .

 2133 2ξ -1)} The bilinear form β T (•, •) satisfies the inf-sup condition.

Proposition 3 . 4 .

 34 Under the same hypotheses as Theorem 3.1 problem ( P) has a unique solution.

1 Figure 2 :

 12 Figure 2: Geometry and boundary conditions for the first case test

Figure 3 :Figure 4 :

 34 Figure 3: Horizontal flow. Pressure field obtained without (left) and with (right) the fractures

u • n = 0 P i = 1 u • n = 0 Figure 5 :Figure 6 :

 01056 Figure 5: Boundary conditions for the third test case

Figure 7 :

 7 Figure 7: Performance of the preconditioner. Left: log of residual vs number of iterations, in red without preconditioning, in blue, with preconditioning. Right: performance for different values of the fracture permeability

The authors are indebted to Y. Achdou for this remark.