MIMETIC FINITE DIFFERENCE APPROXIMATION OF FLOWS IN FRACTURED POROUS MEDIA

We present a possible framework for the numerical simulation of flow in fractured porous media that couples mimetic finite differences for the porous matrix with a finite volume scheme for the flow in the fractures. The resulting method is theoretically analyzed in the case of a single fracture. Moreover, several numerical experiments show the capability of the method to deal also with complicated networks of fractures. Thanks to the implementation of rather general coupling conditions, it encompasses both “conductive fractures”, i.e., fractures with high permeability and “sealed fractures”, i.e., fractures with low permeability which act as a flow barrier. Mathematics Subject Classification. 65N30, 35Q86, 76S05. Received April 15, 2015. Revised October 19, 2015. Published online May 23, 2016.


Introduction
The simulation of underground flows is of great interest for a large number of applications, ranging from energy production to water resources management: oil fields exploitation, geothermal energy, nuclear waste and carbon dioxide storage, and groundwater contamination. In all the aforementioned applications, at very different length scales, the heterogeneity of the porous medium has a major impact on the flow. Geological applications are often characterized by the presence of layers of different materials, with permeability that can span several orders of magnitude within the domain of interest. Moreover, tectonic stresses, or sometimes human activities, produce fractures at different space scales, ranging from micro-fractures up to large fractures and faults. While the local effect of small fractures on the permeability of the rock can be accounted for by homogenization, large features should be explicitly included in the model. The contribution of fractures to the overall flow may differ.
Keywords and phrases. Mimetic finite differences, flow in fractured porous media.
The aim of this paper is to study whether the Mimetic Finite Difference method can be succesfully used to simulate flows in fractured porous media. Our physical model is based on the (mixed form of the) Darcy's equations for the porous medium (or bulk) flow coupled with the (primal form of the) Darcy's equations for the fracture flow. The resulting system of equations is then closed imposing suitable physically consistent coupling conditions along the bulk/fracture interfaces. We present the weak formulation of new coupled problem and prove its well-posedness. We remark that our approach is different from the one usually employed in the literature, where either the mixed or the primal form of the Darcy's equations are considered for both bulk and fracture flow. Our choice is motivated by the fact that the coupling conditions at the interfaces between bulk and fracture flow involve only the fracture pressure. Therefore, the fracture velocity is not required in practice and a primal formulation can be employed within the fracture network, without loosing any information. Moreover, this brings some simplification in the analysis and allows the direct use of existing codes for discrete fracture networks, like the one described in [44], that employ the primal formulation of the Darcy's problem. From the numerical viewpoint, we propose to employ MFD methods to discretize the bulk equations and a two-point finite volume scheme for the fracture network. To show the effectiveness of the proposed approach we numerically test it on two-dimensional test cases.
The paper is structured as follows. In Section 2 we introduce the governing equations. The weak formulation of the coupled problem and its well-posedness are addressed in Section 3 in the case of a single fracture that cuts the domain into two disjoint subdomains. The extension of the proposed model to the case of a network of partially immersed fractures is discussed in Section 4. The numerical method is presented in Section 5 where it is also shown that the resulting fully coupled discrete problem admits a unique solution. In Section 6 we present a set of two-dimensional numerical experiments showing the robustness and optimal convergence properties featured by the proposed scheme. Finally, in Section 7 we draw some conclusions.

Model problem
This section is devoted to the presentation of the governing equations. Throughout the paper we will adopt the standard notation for Sobolev spaces. More precisely, for D ⊂ R d , d = 2, 3, and for a real number s ≥ 0, H s (D) will denote the standard Sobolev space of order s, endowed with the usual norm · H s (D) and seminorm |·| H s (D) . For a (d − 1)-dimensional manifold Γ ⊂ D, we denote by H s (Γ ) the usual Sobolev surface space, cf. [38]. For s = 0 we will write L 2 (·) instead of H 0 (·). In the following the symbol will signify that the inequality holds up to a multiplicative constant that is independent of the discretization parameter.
To describe an incompressible fluid flow problem in a d-dimensional fractured porous media, d = 2, 3, we need the following ingredients: (i) the governing equations for the porous medium (or bulk) flow; (ii) the governing equations for the fracture flow; (iii) suitable physically consistent coupling conditions along the bulk/fractures interfaces; cf. Figure 1 for a sketch of the mathematical model in a three-dimensional configuration (d = 3). Let Ω ⊂ R d , d = 2, 3, be an open, bounded, convex polygonal/polyhedral domain representing the porous matrix. To keep the presentation as simple as possible, we will assume that there is only one (d − 1)-dimensional manifold Γ ⊂ R d−1 representing a fracture that cuts Ω into two disjoint subregions, say Ω 1 and Ω 2 , and that the measure of Γ is uniformly bounded, i.e., |Γ | 1. We remark that the extension to the case of a finite number of (possibly intersecting) fractures or "partially immersed" fractures can be handled in a similar way but in such a case the mathematical model and the functional setting are much more complex; this will be quickly discussed in Section 4. We assume that the boundary of Ω is decomposed into two non-intersecting subsets, i.e., ∂Ω = Γ D ∪ Γ N , with Γ D ∩ Γ N = ∅, and Γ D = ∅, and set Γ D,i = ∂Ω i ∩ Γ D and Γ N,i = ∂Ω i ∩ Γ N for i = 1, 2. Finally, n Γ denotes the unit normal vector to Γ with a fixed orientation (from Ω 1 to Ω 2 ) and τ Γ denotes the R d×d−1 matrix whose columns form an orthonormal basis for the tangent space at each x ∈ Γ .

Governing equations
We first present the governing equations for the bulk flow. To this aim, let K ≡ K(x) ∈ R d×d be the bulk permeability tensor, which is assumed to satisfy the following (classical) regularity assumptions: (i) K is a symmetric tensor whose entries are bounded, piecewise continuous real-valued functions;

Decoupling
Coupling Bulk F re actur ra F F net tw rk or re r tur (ii) there exists κ , κ > 0 such that Given a function f ∈ L 2 (Ω) representing a source term or a sink and g D ∈ H 1/2 (Γ D ), we consider the Darcy's law to model the motion of a incompressible fluid in each domain Ω i , i = 1, 2, with pressure p i and velocity u i : where f i = f | Ωi , K i = K| Ωi , i = 1, 2, and n denotes the outward unit normal vector to ∂Ω. The second ingredient is represented by the governing equations for the fracture flow. We consider a reduced model consisting in modeling the fracture as a (d − 1)-dimensional manifold immersed in an d-dimensional object. Roughly speaking, the reduced model can be obtained writing the Darcy's equations on the fracture in the normal and tangential components and then integrating the tangential component along the thickness Γ ≡ Γ (x) of the fracture domain, which is typically some orders of magnitude smaller than the size of the domain Ω. We refer to [48] for more details. The fracture flow is then characterized by a permeability tensor K Γ , which is assumed (i) to have a block-diagonal structure of the form where κ τ Γ is a (d − 1) positive definite tensor (it reduces to a positive number for d = 2). By this we mean that the material contained in the fracture (before the model reduction process has been carried out) has a permeability that can be written in the stated form, and which may vary along the fracture but not across the thickness of the fracture. We will see later on that in the reduced model κ n γ and κ τ Γ play very different roles and the effective normal and tangential permeability of the fracture scale differently with respect to the fracture thickness; (ii) to satisfy the same condition stated in (2.1) for the bulk permeability, for x ∈ Γ . We recall that n Γ denotes the unit normal vector oriented from Ω 1 to Ω 2 so that n Γ ≡ n 1 , n 1 being the unit normal vector that point outward from Ω 1 . With the above notation, we define, for a regular enough function v, its jump and average across Γ as , and denoting by p Γ and q Γ the fracture pressure and flux, respectively, the governing equations for the fracture flow read which will be used in the following as a starting point to derive the two-point finite volume scheme.
Finally, we provide the interface conditions to couple problems (2.2) and (2.3) or, equivalently, problems (2.2) and (2.5). Following [48], let ξ ∈ [0, 1], then the coupling conditions are given by Motivated by the fact that the coupling conditions (2.6) involve only the pressure in the fracture p Γ and not the flux q Γ , we will focus on the coupled problem (2.2)-(2.5), where the primal form of the Darcy's equations is considered in the fracture. This is different from the approaches considered, for example, in [48], where the mixed form of the Darcy's equations is solved in the fracture. In the next section we then present the weak formulation of problem (2.2)-(2.5) supplemented with the coupling conditions (2.6) and discuss its well posedeness. This will be instrumental to set up the approximation scheme that will be discussed in Section 5.

Weak formulation and its well-posedness
We introduce the following spaces with the associated norms Note that we are requesting more regularity on the velocity v than mere H(div, Ω 1 ) × H(div, Ω 2 ). This is required to accommodate the Robin-type interface condition given by (2.6a), see [48,50]. It can be shown that the spaces Q, W and V 0,∂ΓD are Hilbert spaces with scalar product inducing the above stated norms. For further use, we introduce the space W 0,ΓN = {v ∈ W : v i · n = 0 on Γ N,i , i = 1, 2}. Next, let the bilinear forms a ξ : where we assume that the index i varies in Z/2Z, i.e., 2 + 1 = 1. With the above notation, the weak formulation of problem (2.2)-(2.5) complemented with the coupling conditions (2.6) reads as follows: find u = (u 1 , u 2 ) ∈ W 0,ΓN , p = (p 1 , p 2 ) ∈ Q, and p Γ ∈ V 0,∂ΓD such that Next, we show that formulation (3.1) is well-posed. We first note that for any ψ ∈ L 2 (Γ ) the solution is unique and may be decomposed uniquely as q Γ = q 0 Γ + q 1 Γ (ψ), where q 0 Γ and q 1 Γ (ψ) are solutions of respectively. Furthermore, we have where the hidden constants depend on Γ κ * . It may be recognized that the mapping from L 2 (Γ ) to H 1 (Γ ) defined as ψ → q 1 Γ (ψ) is linear and continuous. Consequently, we can introduce A ξ (·, ·) : W × W → R defined as follows and rewrite (3.1a)-(3.1b) as: find u = (u 1 , u 2 ) ∈ W 0,ΓN and p = (p 1 , p 2 ) ∈ Q such that where we have set and we have used a shorthand notation for the integrals in the right hand sides.
Proof. Bilinearity of A ξ (·, ·) is straightforward, as it is sufficient to note that the map q 1 Γ (·) is linear by construction. Continuity is a consequence of the fact that the following inequalities hold for any u, v ∈ W, where the hidden constant in the first bound depends on κ * . In order to prove coercivity, we first note that for v ∈ W we have that Finally, we note that (3.8) By collecting the results in (3.6)-(3.8) we obtain We remark that the condition ξ > 1/2 has been found also by [48].

Proposition 3.2. The bilinear form
Proof. The result is rather standard and follows the lines of the one given in [48]. For the sake of clarity we report a sketch of the proof in the case of Dirichlet boundary conditions, namely, satisfying the stability estimate z H 2 (Ω) q L 2 (Ω) . Then, we define the velocities w i as the restrictions on where the last bound follows from the elliptic regularity estimate ∇z L 2 (Ω) q L 2 (Ω) . The proof is concluded by observing that, thanks to the trace inequality, where the hidden constant depends on |Γ | 1/2 |Ω| −1/2 , which is uniformly bounded since |Γ | 1.
Before proving that problems (3.1) and (3.5) are well-posed, we preliminarily show the following equivalence result. Proof. We only show that if (u, p) ∈ W 0,ΓN × Q is a solution of (3.5) then (u, p, The converse is straightforward, so we omit the proof. We first observe that given u ∈ W 0,ΓN , problems (3.1c) and (3.3) are equivalent by construction, i.e., where p Γ solves (3.1c) and q 0 Γ , q 1 Γ are given by (3.3). Next, we show that if (u, p) is a solution of (3.5) then (u, p, p Γ ) is a solution of (3.1) where p Γ is defined as in (3.9). To this aim, we take the residual of (3.1a)-(3.1c) for u ∈ W 0,ΓN solution of (3.5) and , and observe that it is zero since it is identical to the residual of (3.5). Next, we show that (3.1) has a unique solution. To this aim assume, by absurd, is the solution of (3.5) for the same data. Clearly, Taking q = ∇ · (u * − u) and using the definition of the bilinear form B(·, ·) we obtain Notice that the identity p * Γ = p Γ could have also been proved by using (3.1c). Indeed, from (3.1c) and and using that [ Next, using the previous results and subtracting (3.1a) and (3.5a) we obtain Taking v = u − u * , we obtain u = u * in W 0,ΓN . Finally, we are only left to show that p = p * . We employ (3.1a) and write it for (u, p, p Γ ) and (u * , p * , p * Γ ). Subtracting the two equations term by term we obtain From the inf-sup condition proved in Proposition 3.2 we obtain p = p * , and the proof is complete.

The case of immersed fracture networks
We now consider the case of a two-dimensional network of immersed fractures, i.e., networks of fractures whose endpoints may not intersect the domain boundary. We refer to the case where at least one endpoint is on the boundary as partially immersed network. A thorough analysis of this case is still ongoing work, some preliminary results may however be found in the literature. For instance, in [4] the case of a single immersed fracture is analyzed, but eventually using a primal formulation only for the pressure. In [24] the authors consider the case of a partially immersed network, yet employing different (simpler) coupling conditions, and a gradient discretization method is proposed for the numerical solution. In this work we limit ourselves to describing the mathematical model and showing, by numerical experiments, that the corresponding discretization by a mimetic/finite volume scheme gives satisfactory results.
A possible situation is the one depicted in Figure 2. In this section Γ indicates the network, which is composed by a set of M Γ fractures, i.e., Γ = ∪ MΓ k=1 γ k , each γ k being (in the 2D case that we consider here) an open segment in R 2 . We further assume that, for j = k that is, the fractures may join only at their end points. Here i kj indicates the intersection point between fracture γ k and γ j , which may be empty if the corresponding fractures do not intersect each other. A fracture may reach the boundary, we assume that the angle formed by any couple of intersecting fracture as well as the angle between fractures and domain boundary is bounded from below by a positive angle. As a consequence, the number of fractures intersecting at an intersection point is bounded. We further indicate with I = ∪i kj the set of all intersection points, while we set , the latter set collecting the "free" fracture endpoints that are strictly contained in Ω. It is understood that some of those sets may be empty, while their union is the whole ∂γ k . For s = D, N, F we define I s = ∪∂γ s k , and for a given intersection point i ∈ I we indicate with S i the set of fractures γ k intersecting in i, i.e., the fractures γ k such that ∂γ I k ∩ i = ∅. On each fracture γ k we can identify two sides, indicated by γ + k and γ − k , respectively, and the two associated normals n + k and n − k = −n + k . We also associate to each fracture a unique normal by taking n k = n + k . To simplify notation we indicate by n ± Γ the normal vectors to the network Γ , i.e., n ± Γ (x) = n ± k (x) if x ∈ γ k . Analogously for n Γ . For a function f in Ω \ Γ we indicate with f ± its traces on Γ ± = ∪γ ± k , respectively. This allows to extend the jump and average operators on the network Γ , i.e., We also extend the previous definitions of K Γ , κ n Γ and κ τ Γ to the network Γ in a natural way. Finally, in the following we indicate with τ k the unit tangent of γ k at its end points, pointing outwards w.r.t. γ k . Now, we are ready to As for (2.5) it is rewritten by imposing continuity of flux and pressure at the intersection points, namely Here, equations on Γ are in fact on each γ k . Furthermore, we have imposed a zero flux condition on the fracture endpoints I F , which are immersed in the matrix domain. As for the interface conditions (2.6), they may be formally rewritten in the same way, where it is understood that they should be applied on each γ k . As already pointed out, it is beyond the scope of this paper to give more details on this more complex, yet more realistic, model. We only add that is it possible to formally derive a weak form having the same structure as (3.1a)-(3.1c) and to extend to it the numerical discretization techniques presented in the next section for the simplified situation of a single fracture. Indeed, in the section on numerical result we will show an example that considers a network of fractures.

Numerical discretization
In this section we present a Mimetic/Finite Volume discretization of the fully coupled problem (3.1). As a first step, we introduce the mimetic discretization of (3.1a) and (3.1b), under the assumption that p Γ is given (see Sect. 5.1). Then, we discuss the finite volume discretization of (3.1c) under the assumption that the velocity field u is known (see Sect. 5.2). Finally, in Section 5.3 we present the Mimetic/Finite volume discretization of the fully coupled problem. To keep the presentation as simple as possible, in this section we consider a two-dimensional case, i.e., d = 2.

Mimetic discretization of the bulk problem
In this section we present the mimetic discretization of (3.1a) and (3.1b) under the assumption that p Γ is given. We first introduce some useful notation. Let T h be a partition of Ω into non-overlapping (possibly non-convex) polygons E, which are aligned with the fracture Γ . This induces a natural partition of T h into two disjoint sets of polygons T h,1 and T h,2 such that T h = T h,1 ∪ T h,2 . In practice, T h can be simply built as follows: first Ω is meshed with a Cartesian grid, then the elements across Γ are simply cut in such a way that the resulting polygonal elements are conforming with Γ . This procedure induces also a subdivision of Γ which we call Γ h . The set of all edges of the decomposition T h is denoted by E h . In order to deal with the coupling conditions (2.6), and also in view of the discretization of the equation in the fracture, detailed in the next Section 5.2, we number the edgesê i ∈ Γ h for i = 1, . . . , N Γ and for eachê i we create two edges e i 1 and e i 2 , geometrical identical toê i , which will be associated to T h,1 and T h,2 , respectively, so that each subset T h,i is complemented with its own fracture edges. In other words, any original fracture edgeê i created by the mesh generation procedure described above is replaced, for the bulk problem, by two edges, e i 1 and e i 2 . In view of this discussion, the set E h can be decomposed as follows  [27] we require that the decomposition T h satisfies the following shape regularity assumptions:

A1
The number of edges of any element E is uniformly bounded; A2 There exists τ > 0 such that every element E is star-shaped with respect to every point of a ball centered at a point C E ∈ E and with radius τh E ; A3 For any element E and for any edge e ⊂ ∂E it holds |e| h E , where |e| denotes the length of e. Now let us introduce the mimetic spaces. We denote by Q h the discrete space representing the degrees of freedom of the scalar variables. More precisely, we associate the degrees of freedom of the scalar variable to mesh cells so that for For further use we introduce the jump operator across each fracture edge, i.e., We also introduce two projection operators, denoted by the superscript I, from L 1 (Ω) and H(div, Ω) onto Q h and W h , respectively, as follows where, for any edge e ∈ E h , n e is a unit normal vector assigned to e once and for all, see [33]. We also define the mimetic discrete divergence operator DIV h : where G e E = G e n e · n e E ∈ R, being n e E the unit vector normal to e pointing out of E. This definition is consistent with the Gauss divergence theorem. We also recall that it holds, cf. [29] for details.
Next, we define suitable scalar products onto the discrete spaces Q h and W h . On Q h we set which corresponds to the L 2 (Ω) scalar product for piecewise constant functions. The scalar product in W h is defined by assembling elementwise contributions from each element, i.e., where, by following [27,29], the local scalar product [·, ·] E can be defined in such a way that the following two conditions are satisfied: (S1) Stability: for all G ∈ W h and for every element (S2) Local consistency: for every linear function q 1 on E ∈ T h it holds We are now ready to state the mimetic discretization of problem (3.1a) and (3.1b) under the assumption that p Γ is given. To this aim let be the approximation of the fracture pressure inê i that will be computed by a finite volume scheme, as explained in the next section. Then our formulation reads as: find F h ∈ W h 0,ΓN and p h ∈ Q h such that where the index j varies in Z/2Z. The right hand sides in (5.8) are defined as where g e D = 1 |e| e g D ds, and f = f I is the vector of the mean values of f , defined according to (5.3).

Finite volume discretization of the fracture problem
To discretize (2.5) we employ a finite volume formulation. We consider the partition Γ h of Γ induced by the mimetic discretization introduced in the previous section and we set h Γ,i = |ê i | > 0 and set h Γ = max i h Γ,i . We consider a parametric description of Γ h and indicate with s the arc-length coordinate. We assume to have numbered the fracture edges so that s i is the center ofê i , with s 1 < s 2 < . . . < s NΓ , and s i−1/2 and s i+1/2 are the two end points ofê i . Clearly, h Γ,i = s i+1/2 − s i−1/2 . If we consider an edgeê i fully contained in Γ h , i.e., with i ∈ {2, . . . , N Γ − 1}, the integration of the differential equation in (2.5) overê i gives . As explained in the previous section, for eachê i we have two duplicated edges e 1 1 ∈ E Γ h,1 and e 2 2 ∈ E Γ h,2 , and the jump defined as in (5.2) is constant inê i . We assume that also κ Γ is piecewise constant, with values {κ i , i = 1, . . . N Γ }, and we set In particular, we consider the so-called two-point flux approximation, where being T i+1/2 the so-called transmissibility betweenê i andê i+1 . To evaluate it, let us first consider the case whereê i andê i+1 lie on a straight segment. Since the velocity is well defined at the interface, we have Note that, since κ Γ and ∇ τ p Γ may be discontinuous across elements, we are here taking the left and right limits. Assuming that ∇ τ p Γ is continuous on each cell, by performing a Taylor expansion around s i+1/2 inê i andê i+1 , respectively, we obtain where p i = p Γ (s i ). By manipulating the two expressions to eliminate p i+1/2 , neglecting the o(h Γ ) term, and recalling that H i+1/2 is an approximation of u i+1/2 , we derive the following expression for the trasmissibility Note that if p Γ has continuous second derivative on each cell then the formula for the numerical fluxes is in fact second order accurate with respect to h Γ . Indeed, in this case when eliminating p i+1/2 the term of order h Γ in (5.10) cancels out, leaving a remainder of higher order. Ifê i andê i+1 form an angle ζ i+1/2 , we modify (5.12) according to the recipe suggested in [44], and replace (5.11) with In view of the above discussion, equation (5.9) is then approximated by We have now to handle the boundary cells. If a boundary cell is adjacent to ∂Γ N we still use (5.13), but we set to zero the numerical flux at the corresponding cell boundary. For a Dirichlet condition, we use again a Taylor expansion. If, without loss of generality, we assume thatê 1 is adjacent to ∂Γ D , the resulting equation is (5.14) In the case of networks of fractures we have to cope with the intersection of three or more fractures. Each fracture is meshed independently and we have |S i | cells that meet at intersection i, corresponding to the fractures intersecting at that point. For simplicity (and no loss of generality), let us call themê γ k ,i , with γ k ∈ S i , and indicate with α γ k ,i the corresponding coefficient computed according to (5.12). The transmissibility coefficient T kj for the degrees of freedom associated to C γ k ,i and C γj ,i , with γ k , γ j ∈ S i , is now being ζ kj the angle between γ k and γ j at the intersection. Clearly T kj = T jk .

Fully-coupled problem and its algebraic formulation
Let N f , N p be the dimensions of the discrete spaces W h and Q h , respectively, and recall that N Γ is the number of edges belonging to Γ h . The final algebraic system stemming from the mimetic/finite volume discretization has a saddle point structure. Indeed, problem (5.8) is algebraically equivalent to Finally, b f and b p collect the contributions to the right hand side (we have changed the sign in discrete divergence equation to recover the classical matrix structure for saddle point problems). Since the grid used for the fracture problem is conforming to the mimetic grid, it is immediate to recognize that assembling (5.13) and (5.14) and using the given definitions of the numerical flux produces the linear system where T ∈ R NΓ ×NΓ assembles all the transmissibility terms, and b Γ collects all contributions to the right hand side due to the forcing term and Dirichlet boundary data. We can then write ⎡ Now, thank to the hypothesis on Γ D , the matrix A is symmetric and positive definite, while B T has zero null-space because of the inf-sup condition, cf. Proposition 3.2. Furthermore, we have the following result.  ]ê i = 0 ∀i}. Indeed, imposing a zero jump on the velocity is equivalent to applying the discrete divergence operator to a standard Darcy field on Ω without fractures, and for this case we can use standard results on mimetic finite differences. This observation implies that, by setting N f = dim( W h 0,ΓN ) for any 0 = q h ∈ R Np , we may find a G ∈ R N f such that q h = BG, and thus G T B T q h = q h 2 > 0, being · the standard Euclidean norm. Furthermore, for such a G we have that G T C T q h = 0. Consequently, for any q h = [q T h , q T Γ,h ] T ∈ R Np×NΓ different from zero we may consider the following two cases. If q h = 0 then we can select a G ∈ W h 0,ΓN such that We next show a result on the transmissibility matrix T that we will need later on.

Proposition 5.2. The matrix T is symmetric and semipositive definite (with kernel formed by constant vectors)
Proof. It can be verified that, by construction, the matrix T is symmetric and is a Z-matrix with positive diagonal elements.
If ∂Γ D = ∅ we have that T ii = − j T ij for all i, so it is diagonally dominant with the vector [1, 1, . . . , 1] T in the kernel. Otherwise, at the Dirichlet end points equation (5.14) implies that at least one row i satisfies T ii > j |T ij |. The thesis then follows from standard linear algebra results.
The next result shows that the linear system (5.17) admits a unique solution. Proof. We can exploit a well known result for saddle point problems [18]. Indeed, we can reformulate the governing matrix in (5.17) as where T = 0 0 0 T is symmetric and semipositive definite, cf. Proposition 5.2, and C has been defined as in (5.18). Then, K is non-singular if and only if ker( C T ) ∩ ker( T) = {0}. This is automatically satisfied since ker( C T ) = {0} because of Proposition 5.1.

Remark 5.4.
Note that to prove the well-posedness at the continuous level, we assumed that ∂Γ D = ∅, cf. Section 3. From the analysis of the discrete problem, we can conjecture that this condition can be relaxed and that the coupling conditions are sufficient to guarantee the existence of a unique discrete solution. Indeed, the well-posedness for the problem with a single immersed fracture, where no Dirichlet conditions are imposed on the fracture, has been already obtained in [4], even if by using a different formulation.

Numerical results
In this section we describe the numerical results obtained by employing the previously discussed models. The numerical models have been implemented in a software written in C++ language. For the generation of the meshes we have used of the CGAL library [51], while for matrix manipulation and linear system solution we exploited the Eigen library [41].
Throughout this section the parameter ξ appearing in the definition (2.6) of the coupling conditions has been chosen as ξ = 0.75. Moreover, the fracture thickness has been set equal to Γ = 0.01.

Example 1
The first test case is based on the example proposed in [42], with slight modifications. We take Ω = (−1, 1) × (−1, 1), Γ = (−1, 1) × {0} and impose homogeneous Dirichlet boundary conditions on the whole boundary of Ω, i.e., Γ N = ∅ and Γ D = ∂Ω; the boundary conditions for the fracture problem are imposed accordingly. The bulk permeability tensor is assumed to be the identity matrix, i.e., K = I, whereas the fracture permeability tensor is chosen as K Γ = εI, for a positive real number ε that will be specified later on. Notice that, according to the definition given in (2.4), the permeability of the reduced model (2.5) is given by κ Γ = Γ ε. Moreover, with the above choice of K Γ the parameter η Γ appearing in the coupling conditions (2.6) becomes η Γ = Γ /ε. We take the source terms as so that the exact solution is given by We have tested our numerical scheme on a sequence of unstructured triangular (Grid I) and polygonal grids (Grid II) with granularity h ≈ 1/N for N = 2, 4, 8, 16, 32, 64. Any polygonal decomposition has been obtained from a triangular grid by merging randomly two or more triangles leading then to a decomposition containing elements that can have three, four or five edges, cf. Figure 3a for an example. We measure the relative approximation errors for the pressure and the velocity in the bulk domain as: where p I and u I are the interpolants of the exact solution (pressure and velocity) in the mimetic spaces Q h and W h , respectively. In Figure 4 we report the computed relative errors err p and err v as a function of the meshsize (loglog scale) for ε = 0.001, 1, 1000. The results reported in Figure 4 clearly show that a second order convergence rate for the pressure variable is clearly achieved; this superconvergence effect has been already observed in many cases, see [28] and [15] for examples. As concerns velocity we have at least first order convergence in all cases. The mean convergence rates are summarized in Table 1.

Example 2
In this second test case we investigate the robustness of our scheme with respect to the contrasts between the permeability in the bulk and in the fracture. We let Ω = (0, 1) × (0, 1). On Γ D = {0, 1} × (0, 1) we impose a nonhomogeneous Dirichlet condition g D = y whereas on Γ N = [0, 1]×{0, 1} we impose a homogeneous condition, i.e., g N = 0. The right hand side is chosen as f (x, y) = 4 in Ω. The fracture is given by Γ = {(x, y) ∈ Ω : y+2x = 1.4}. We take the medium to be isotropic and set K = I, and let the parameters κ n Γ and κ τ Γ vary. Note that κ n Γ and κ τ Γ completely determine the permeability of the fracture and therefore its behavior when embedded in the porous medium. We have considered the following cases: (i) κ τ Γ = 1, κ n Γ = 1. In this case the pressure is expected to be almost constant across the fracture Γ because the permeability in the normal direction is equal to that of the surrounding medium and the fluid is not affected by the presence of the fracture. This is exactly the behavior observed in Figure 5a; (ii) κ τ Γ = 1, κ n Γ = 0.01. Here the normal permeability is smaller than that of the surrounding medium, therefore we expect a pressure jump across the fracture, cf. Figure 5b where this behavior is clearly observed; (iii) κ τ Γ = 100, κ n Γ = 1. In this case the fracture pressure is expected to be almost linear and the flow is expected to be directed towards the fracture. This because here the fracture is very permeable in the tangential direction and the fluid is "attracted" by the fracture. This is exactly the behavior observed in Figure 5c.
In all the aforementioned cases we compare the result obtained with the mimetic method with the analogous ones computed with the XFEM approach described in [34] where a mixed finite element formulation is combined with a suitable enrichment on the elements crossed by the fracture. As shown in Figure 5 the results obtained with the two methods are in good agreement in all the cases. Meshes with a similar size h have been used in the two approaches: for the mimetic method we use a polygonal grid obtained with the agglomeration of a constrained Delaunay triangulation, while for the XFEM, since we do not need to honor the fracture geometry, we use a structured triangular grid.

Example 3
In this example we aim at assessing the proposed scheme on a more complex network of fractures. We compare the solution of our scheme with a discretization where finite volumes have been employed to discretize both the bulk as well as the fracture equations. We choose a domain Ω = (0, 1) × (0, 1) containing 10 fractures, see Figure 3b. The finite volume scheme adopted employs a simple two-point flux approximation. It is known that  this scheme can be unaccurate unless special grids are employed. In particular, a sufficient condition for the consistency of two-point flux approximation for an isotropic permeability tensor requires that the normal to any interior edge be directed along the segments joining the centroids of the adjacent elements. To satisfy this requirement as far as possible for the finite volume scheme we have generated the grid using the constrained Delaunay triangulation tool of the CGAL library. For more information on the finite volume method, the interested reader may consult [39], or the recent review [35]. We set homogeneous Dirichlet boundary conditions where ε can vary. We have considered the following test cases: (i) no fractures (Fig. 6a); (ii) ε = 1000 (Fig. 6b); and (iii) ε = 0.001 (Fig. 6c). The corresponding discrete pressures computed with the mimetic finite difference scheme are reported in Figure 6(left) and the analogous results obtained with the finite volume method are shown in Figure 6(right). From these results we can conclude that in all the cases the results produced by mimetic finite differences are consistent with those obtained employing finite volumes. Even when no fractures are present in the domain, the solution computed with the finite volume method seems to exhibit higher pressure peaks than those obtained with mimetic finite differences, cf. Figure 6. This might be due to the fact that on grids of similar density the two-point finite volume method considered here is less accurate and tends to produce smaller transmissibility. It is not the scope of this work to analyze this matter in more detail, since we are showing this example just to illustrate the suitability of the proposed scheme also in the case of networks of fractures.
If ε = 1000 the fractures are much more permeable than the surrounding medium. Indeed, pressure is practically continuous across the fractures and the maximum and minimum values of pressure are slightly lower with respect to the non-fractured case, both for mimetic finite differences and finite volumes. Finally, if ε = 0.001, i.e., the fractures are much less permeable than the bulk and we can observe strong pressure jumps across fractures. Once again, finite volume predicts higher peaks with respect to the mimetic finite difference method but the solutions are in good agreement.

Conclusions
In this paper we have presented a possible framework for the numerical simulation of flow in fractured porous media, by coupling mimetic finite differences for the porous matrix with a model of the flow in the fractures based on a pressure formulation and a finite volume scheme. The reason for using different formulations for bulk and fracture flow is that we obtain a rather effective scheme to describe the flow in the rock matrix while accounting for the presence of the fractures. Even if the analysis has been carried out so far for the problem with a single fracture we have shown with numerical experiments the capability of the scheme to deal with fracture networks.
Mimetic finite differences are indeed a natural choice to deal with this type of problems. The intersection of the fractures with an underlying grid for the matrix produces polygonal elements where we can apply the method directly, without the need of complex mesh generation algorithms to produce a standard grid conforming to the fractures. We wish to point out that, inspired from the model presented in [48], we have implemented rather general coupling conditions between the flow in the fractures and in the solid matrix. They may account for both "conductive fractures", i.e., fractures with high permeability, and "sealing fractures", i.e., fractures with low permeability which act as barrier to the flow. This differs from other works that adopt more simplified coupling conditions, which are justified only for very permeable fractures. Further on going developments include carrying out the full analysis of the problem with fracture networks and extending the implementation to threedimensional problems. In this work, having carried out the numerical experiments only in 2D, we have solved the algebraic system with a direct multi-frontal scheme. Moving to the more challenging three-dimensional problems will require also to investigate suitable preconditioners to accelerate a Krylov-subspace iterative solver.