Preconditioning Techniques for the Numerical Solution of Flow in Fractured Porous Media

This work deals with the efficient iterative solution of the system of equations stemming from mimetic finite difference discretization of a hybrid-dimensional mixed Darcy problem modeling flow in fractured porous media. We investigate the spectral properties of a mixed discrete formulation based on mimetic finite differences for flow in the bulk matrix and finite volumes for the fractures, and present an approximation of the factors in a set of approximate block factorization preconditioners that accelerates convergence of iterative solvers applied to the resulting discrete system. Numerical tests on significant three-dimensional cases have assessed the properties of the proposed preconditioners.


Introduction
The simulation of underground flows in fractured porous media is of great interest for a large number of geophysical applications, such as oil production, CO 2 storage, and groundwater contamination and remediation. It is well-known that the presence of fractures and/or faults strongly influence subsurface flows. The major challenges from the numerical viewpoint are represented by (i) geometric complexity, and (ii) strong heterogeneity of materials at different space scales. While micro-fractures can be accounted for by means of homogeniza-tion/upscaling techniques, large fractures and faults play a more complex role, acting as paths or barriers for the flow, and therefore they have to be included in the model explicitly. Fractures are characterized by a small aperture compared to their typical length and the size of the domain, thus a widely employed approach consists in modeling them as (d − 1)-dimensional interfaces immersed in a d-dimensional porous medium (indicated in the following as the bulk). A reduced (d − 1)-dimensional problem is then solved on the surfaces representing the fractures, with physically-consistent coupling conditions accounting for the exchange of fluid between fractures and porous medium.
From the computational viewpoint, this dimensionally-hybrid setting avoids the need for extremely fine grids to resolve the fracture's scales. Assuming that the fractures are filled by a porous medium with its own porosity and permeability, Darcy's law can be used for modeling both d-dimensional bulk and the (d − 1)-dimensional fracture flow problems. The first dimensionally-hybrid model for flow in fractured porous media has been proposed in [4] in the case of a very permeable fracture. Later on, in [51] it has been generalized to fractures featuring low permeability. These models were derived based on the assumption that there is one single fracture cutting the bulk domain in exactly two non-overlapping subdomains; the extension to the case of a fully immersed fracture has been analyzed in [6]. This dimensionally-hybrid model has also been extended to the case of a two-phase flow in [43,47]. One of the main issues concerning discretization of the flow in heterogeneous media is mesh generation. Indeed, the grid has to be conforming with the fractures, but whenever the number of fractures is very large such a constraint can result in a unaffordable computational burden, particularly whenever fractures feature small intersection angles, or they are nearly coincident. Indeed, in such cases, the conformity constraint can lead either to very fine grids, or to low-quality elements (small angles, high aspect ratios). To overcome this difficulty a possible strategy consists in the use of numerical schemes that can support arbitrarily shaped polygonal and polyhedral meshes, and that guarantee good approximation properties also in presence of low-quality mesh elements.
Indeed, in recent years, the exploitation of computational meshes composed of polygonal and polyhedral elements has become very popular in the field of numerical methods for partial differential equations because the flexibility they offer allows for the design of efficient computational grids when the underlying problem is characterized by a strong complexity of the physical domain. Several conforming and non-conforming numerical discretization methods which admit polygonal/polyhedral meshes have been proposed in recent literature; here, we mention, for example, Mimetic Finite Differences (MFD) [12,18,28,30], high-order polyhedral Discontinuous Galerkin (PolyDG) methods, [10,11,13,32], Virtual Element methods (VEM), [15,16], and the Hybrid High-Order (HHO) method, [33,36].
In this paper, continuing the work initiated in [12] we focus on Mimetic Finite Differences, which in the last years have been successfully applied to a wide range of problems, as for example diffusion-type problems [7,17,28,29], electromagnetism [27], plate equations [19], non-linear and control problems [8,9], and to model two-phase flows [48]. We refer to [14,18] for a comprehensive review on MFD schemes. In the context of numerical modeling of flows in fractured porous media, MFD have been successfully used in [3,12,42]. Several other numerical techniques have been proposed in recent years for fractured media flow, reflecting the importance of the subject for various applications. With no claim of completeness, in addition to the already cited literature we mention here some recent works by several research groups [2,[24][25][26]56].
We will consider the formulation proposed in [12], where a mixed MFD approximation for the coupled Darcy's model is analyzed for a fully immersed fracture network. The linear system of equations stemming from the discretization has the form of a generalized saddle-point system, sometimes referred to as a double saddle-point problem [5], reading for suitable matrices that will be defined in Sect. 1, and where u, p and p Γ contain the approximate solution for the bulk velocity, bulk pressure, and fracture pressure, respectively, while the vectors g, h and h Γ contain the terms arising from the forcing and boundary data.
In this paper we analyze the spectral properties of the system of equations stemming from the considered MFD-FV discretization. More precisely, we prove that, as expected, the condition number of the double saddle-point system depends on the contrast of the permeability in the bulk and in the fractures, and, asymptotically, it grows as O(h −3 ), h being the characteristic mesh size. To reach this result we extended the work of [50] and we make use of a conjecture, so far verified only numerically.
We then address the problem of efficiently solving the above linear system of equations by devising suitable preconditioning algorithms to accelerate the convergence of iterative schemes. We have chosen classical approximate block factorization (ABF) preconditioners because they can be readily implemented in existing codes since they make use of quantities directly available. We propose a technique to construct the approximation of the factors, which, despite its simplicity, has proved to be rather effective.
In the context of preconditioners for fractured porous media simulations, we mention the recent work [31], where a set of norm-equivalent preconditioners [49] are presented, though for a different approximation scheme than the one adopted in this work, some of which show some analogies to the ones proposed here. In the context of saddle-point problems arising in geomechanics, we mention also block preconditioners based on the use of approximate inverses, like the ones adopted in [22]. However, in this work we decided to focus on simpler factorized block approximations straightforwardly implementable in an existing code for flow in fractured porous media.
The paper is structured as follows. In Sect. 1, we introduce the hybrid-dimensional mathematical model and its numerical approximation based on the approach of [12]. In Sect. 2 we discuss the spectral properties of the resulting system matrix. In Sect. 3, we then present some techniques for preconditioning the linear system of equations. The performance of the proposed preconditioners is then assessed in Sect. 4 on three-dimensional test cases. Finally, in Sect. 5 we draw some conclusions.

Mathematical Model and Its Numerical Approximation
The mathematical model we consider in this work follows [42] and [12]. We will recall it briefly, for the sake of completeness. Let Ω ⊂ R 3 be an open, convex polyhedral domain representing a porous medium saturated by a fluid. The medium is fractured, and the fractures are modeled as a collection of planar surfaces. More precisely, with Γ we denote the fracture network given by the union of M Γ fractures γ k , for k = 1, . . . , M Γ , where γ k is a 2dimensional planar open domain embedded in R 3 , i.e. Γ = M Γ k=1 γ k . We indicate with i i j = ∂γ i ∩ ∂γ j the intersection between fracture γ i and fracture γ j , possibly being the empty set, and with I the set of all intersections i i j with non-zero 1measure. Finally, Ω Γ = Ω \ Γ denotes the part of the domain occupied by the rock matrix, which in the following we will indicate as bulk.
Flow in the bulk and in the fractures is assumed to be governed by Darcy's law. We are thus assuming that fractures are filled by a porous medium with different porosity and permeability with respect to the surrounding porous matrix. We denote with K a symmetric positive definite permeability tensor, which we assume to be piecewise constant in Ω Γ , and with f a possible source term. We assume also that the fluid and the medium are incompressible and we neglect gravitational effects. Note that, since the fluid viscosity μ is considered constant, we have defined K = μ −1 K, with K being the actual material permeability tensor.
As for the fractures, we use the reduced model originally developed in [51] for a single fracture, extended to fracture networks as explained in [42]. In particular, we assume that on each γ k we can identify a normal vector n k and that we can decompose the permeabilityK k of the material in the fracture into a normal componentK n k and a symmetric semi-definite tensorK τ k so thatK k =K n k n k ⊗ n k +K τ k , withK n k > 0 and v TK τ k v > 0 if v × n k = 0, whileK τ k n k = 0. Following [51], we may define on each fracture the fracture aperture l k > 0 (assumed constant in each γ k ), effective tangential and normal permeability K τ k = l kK τ k , and K n k = l −1 kK n k , as well as the normal effective resistivity η k = K −1 n k = l kK −1 n k . For all functions defined on each γ k we will use normally the subscript Γ to denote their direct product on the whole Γ , Following [12], we employ a mixed formulation in the bulk, where the unknowns are the Darcy velocity u and pressure p, while in the fracture network we use a primal formulation where the only unknown is the pressure, indicated by p Γ : Γ → R.
To describe the equations and the coupling terms we need to define the jump and average operator across the fractures. Let v : Ω Γ → R be a regular function so that on each Note that we can identify a positive and negative side, γ + k γ − k , of each γ k , so that v + and v − are in fact the traces of v on γ + k and γ − k , respectively. We partition ∂Ω Γ into ∂Ω N and ∂Ω D , where we impose conditions on the normal fluid velocity and pressure, respectively. We assume that |∂Ω D | > 0, where |∂Ω D | here indicates the 2-measure of Ω D . As for the fracture network, we can identify three parts of ∂Γ : ∂Γ N and ∂Γ D are the portions of ∂Γ ∩ ∂Ω where we impose the flux or the pressure, respectively, while ∂Γ 0 is the part of the boundary of the fracture network fully immersed in the bulk. Here, we follow the usual practice of imposing zero flux.
We can now write the differential problems in the bulk and in the fracture networks, complemented by the coupling conditions which will be detailed later on: Here, ∇ Γ · and ∇ Γ denote the tangential divergence and gradient operators, while τ Γ is the unitary normal to ∂Γ parallel to the fracture tangent plane.
We now provide the interface conditions that couple the model for flow in the bulk, cf. (1a), with that in the fractures, cf. (1b). For a ξ 0 ∈ (0, 1/4], they are given by The derivation can be found in the cited references. The closure parameter ξ 0 > 0 depends on the assumption made on the variation of pressure across the fracture aperture when deriving the reduced model. 1 For what concerns the conditions at the intersection between fractures, several solutions are possible. For instance, in [54] and in [41] special conditions were studied to account for possible strong variations of permeability between fractures (however limited to the two-dimensional case). Here, for the sake of simplicity, we assume continuity of pressure and flux balance at each intersection, a common choice for discrete fracture network simulations.

Weak Formulation
To define the functional setting of the problem we first note that for all p ∈ [1, ∞], an element of L p (Ω Γ ) may be identified with an element of L p (Ω), since Γ is of zero measure. Furthermore, we state some regularity assumptions on the data. To this purpose we indicate with * and * positive upper and lower bounds of the corresponding quantity. We then require for all ζ ∈ R d with ζ · n Γ = 0 and a.e. on Γ . We introduce the following functional spaces for pressure and Darcy velocity in the bulk, equipped with the norms ||q|| Q = ||q|| L 2 (Ω Γ ) , and ||v|| 2 W = ||div v|| 2 L 2 (Ω Γ ) + ||v|| 2 L 2 (Ω Γ ) + ||{v · n Γ }|| 2 L 2 (Γ ) + || v · n Γ || 2 L 2 (Γ ) . They are Hilbert spaces with scalar products inducing the stated norms.
To account for boundary conditions, we define W 0,∂Ω N = {v ∈ W : v · n = 0 on ∂Ω N }, where n is the unit outward normal to Ω and v · n is intended in the sense of traces of elements of H div (Ω Γ ). For the forcing term and the boundary data we take f ∈ L 2 (Ω Γ ), φ ∈ L 2 (∂Ω D ) and g regular enough such that there exists an element of H 1 (Ω Γ ) whose trace on ∂Ω N coincides with g.
As for the model of the fracture, on each γ k we define Z k = {q : q ∈ L 2 (γ k ), ∇ τ q ∈ L 2 (γ k )}, and We take φ Γ ∈ L 2 (∂Γ N ) and g Γ regular enough to be a trace on Γ D of functions in Z Γ .
We define the following bilinear forms: for u, v ∈ W , q ∈ Q and q Γ , and, to write the weak formulation of the problem in a more compact form, we define Q = Q× Z , Q 0 = Q× Z 0,Γ D and the form B : We introduce also r φ ∈ W and r φ Γ ∈ Z Γ that represent the lifting of the boundary data on ∂Ω N and ∂Γ D respectively. Following the steps illustrated in [12,42], we can write the following weak formulation of our problem where Theorem 1 Problem 1 is well posed.
Proof Well posedness may be proven by exploiting the results in [42]. Since |∂Ω D | > 0, in the cited work it is shown that a ξ is continuous and coercive in W 0,∂Ω N , and B is continuous and satisfies an inf-sup condition. The continuity of the functionals at the right hand side can be assessed by standard techniques, while d is clearly a semi-positive definite form. Well posedness then follows, see also [23].

Numerical Scheme
Let Ω h be a partition of Ω Γ into non-overlapping polyhedra (bulk cells) P, conforming to the fracture Γ . We define the mesh spacing as h = max P∈Ω h h P , where h P is the diameter of cell P. Let F h be the set of facets of the cells in Ω h , i.e. any f ∈ F h is a facet of a P ∈ Ω h . To be able to represent jumps and average values, the facets laying on Γ will be doubled. More precisely, F Γ h = Γ ∩ F h is formed by pairs of geometrically identical facets (f + , f − ) that cover Γ . This also implies that Ω h induces a natural partition of Γ , that we indicate by Γ h , into planar facetsf, called fracture cells, and for anyf ∈ Γ h there is a couple of bulk facets f + (f), f − (f) ∈ F h that coincide geometrically withf. The set F h may then be partitioned as collects the boundary facets, i.e. Dirichlet and Neumann facets, F 0 h collects the internal facets and F Γ h = F Γ h,+ ∪ F Γ h,− is the union of the decoupled fracture facets. We denote by F h (P) ⊂ F h the facets of a P ∈ Ω h and N f P = card(F h (P)). For the cells in Ω h we make the following assumptions [18]. Assumption 1 There exist two positive real numbers N s and ρ s independent of h such that Ω h admits a conforming sub-partition T h made of tetrahedra such that: -each polyhedron P is star-shaped with respect to a point x P ∈ P and each facet f is star-shaped with respect to a point x f ∈ f. -every polyhedron P ∈ Ω h admits a decomposition T h | P of at most N s tetrahedra. Moreover, the sub-partition T h is simple, i.e., it is built in the following way. Firstly, each facet f is subdivided into triangles by connecting each vertex of f with x f . Secondly, each element P is decomposed into tetrahedra by connecting each vertex of P and each point x f , with f ∈ F h (P), to the point x P ; -every tetrahedron T ∈ T h is regular, i.e. the ratio between the radius r T of the inscribed sphere and the diameter h T is bounded from below by ρ s , i.e.
These assumptions impose some limitations on the shape of the admissible elements, which however are not too restrictive. Indeed, the grid Ω h may contain rather generally shaped elements, like non-convex cells. For the forthcoming analysis, we make the following assumption.

Assumption 2
The meshes Ω h and Γ h are aligned with the discontinuities of the piecewise constant permeability tensors K and K Γ , respectively. Moreover, Γ h satisfies a K-orthogonality property in the sense of [1], i.e. there exists a set of control points {x f } f∈F h such that for any f ∈ F h , x f ∈f and, for any pair of neighboring elements facets f, f ∈ F h sharing an edge σ of normal n, it holds that Kn is parallel to the segment joining the control points x f and x f . 2 We refer to the recent paper [20] for details on effective strategies for mesh generation and implementation aspects in the context of Virtual Element approximations of coupled multi-dimensional flow problems.
We denote by |P| and |f| the volume of the polyhedron P and the area of the facet f, respectively. For every facet f ∈ F h we consider a unit normal vector n f , and in particular, for every couple of fracture facets f + , f − ∈ F Γ h , we set n f + = n f − = n Γ , while we denote by n P,f the unit normal vector on a facet f ∈ F h (P) of cell P, and we set α P,f = n f · n P,f , while we indicate with x P and x f the barycentres of P and f, respectively.
To simplify the set-up of the algebraic system stemming from the discretization of (1), we restrict ourselves to the case ∂Ω N = ∅. The more general case may be treated by standard techniques that do not alter the properties of the numerical scheme (as long as |∂Ω D | > 0).
We approximate the discrete pressure fields with piecewise constant values on each cell, and the bulk velocity with constant normal values on each facet. Let h the vectors of discrete pressure in the bulk and in the fracture network, respectively. We will indicate the global space of discrete pressures as With p P , respectively pf, we indicate the approximation of the pressure on bulk cell P and on the fracture network facetf, that is As standard in MFD methods the degrees of freedom for velocity in the bulk approximate with a constant value the normal velocity across each cell facet. Therefore, if u h ∈ W h is the vector of velocity degrees of freedom, we indicate with its component of u h associated to facet f. In the following, whenever convenient, we will use the notation P i ,f i and f i to indicate the ith bulk cell, fracture cell and bulk facet in Q Ω h , Q Γ h and W h , respectively. The derivation of the mimetic discretization is based on the definition of inner products on Q Ω h and W h . More precisely, we have and, for v h and w h in W h , Here, are the jumps and averages of v h across the facetf, ηˆf is the value of η Γ on facetf, while M and E are the matrices that contribute to the mimetic inner product matrix M c = M + E. Matrix M is the classical MFD inner product matrix, which is built by assembling the cell-wise contributions M P , while E accounts for the contribution due to the coupling conditions (2), and is defined implicitly by As for M P its construction follows [18], and is briefly described. Let P ∈ Ω h and f 1 , . . . , f N f P be the facets in F h (P). We define the matrices Z P ∈ R N f P ×d and R P ∈ R N f P ×d as Here, K P is the value of the permeability K on cell P, c f i , and c P are the barycenters of f i and P, respectively. The elemental mimetic inner product matrix is then given by with It can be proven, see [42], that M c is symmetric positive definite and that the inner product in (8) satisfies the requirements of consistency and stability necessary for a proper MFD discretization. Another step in a MFD discretization is the definition of the discrete divergence operator. For the case of fractured media, it has to represent the discretization of both the divergence operator and the flux exchange term appearing in the equation governing flow in the fractures. We thus define D I V h : W h → Q h and an inner product on Q h so that, for any The latter relation, together with (13) and (7), allows us to construct the matrices Details may be found in [12,42].
As for the equations in the fracture network, we have employed the finite volume scheme with the two point flux approximation described in [45]. For each fracture cellf i ∈ Γ h , we identify the set of edges off i as J (i) = {e i : e i edge off i , e i ⊂ ∂Γ N }, and e i j will indicate the jth element of J (i). Clearly, j ranges from 1 to a maximum of 3 if the facet has no edges on the Neumann boundary. For each e i j ⊂ ∂Γ D we indicate withpˆf (e i j ) the approximated fracture pressure in the only cellf(e i j ) such thatf i ∩f(e i j ) = e i j . While, if e i j ⊂ ∂Γ D ,pˆf (e i j ) indicates the average value of the Dirichlet datum on e i , and eventually it contributes to the right hand side.
The approximation of −∇ Γ · (K Γ ∇ Γ p Γ ) by the finite volume scheme may then be written as − Matrix T, of components T i j , is the so called transmissibility matrix, whose computation is explained in details in the given references. It is symmetric and positive semi-definite (positive definite if ∂Γ D = ∅).

The Algebraic System
From now on, to simplify the notation, we omit the subscript h to indicate the vectors of degrees of freedom. Thanks to the definitions in (8), (14) and (15), we may write the linear system governing the discrete problem as where we recall that the vectors u, p and p Γ contain the approximate solution for the bulk velocity, bulk pressure, and fracture pressure, respectively, while the vectors g, h and h Γ contain the terms arising from the forcing and boundary data. (16) is well posed.

Theorem 2 System
Proof Let us first note that we can reduce the linear system to a classic saddle-point algebraic problem. We define the following block matrices and the vectors π = [ p, The system can be rewritten as Following the steps illustrated in [42], we may prove that M c is symmetric positive definite and ker( B T ) = ∅. Since T is positive semi-definite, well posedness follows from standard results on saddle-point systems, see for instance [21,23]. A similar result may be found also in [12]. 3

Spectral Properties of the Governing System of Equations
In this work we will mainly focus on three-dimensional test cases, since their solution is more challenging. That is why the numerical scheme has been introduced for that setting. However in this section, for the sake of generality, we report the results for a generic space dimension d, with d = 2 or d = 3. It is understood that for the two-dimensional case the assumptions made on the mesh have to be reinterpreted appropriately. Moreover, we will keep the notation used for the three-dimensional setting. In the following, with a b (respectively, a b) we indicate the existence of a positive constant C, independent of h, such that a ≤ Cb (a ≥ Cb), while a b means b a b. We also recall that for any family of meshes {Ω h , h > 0} that satisfy Assumption 1 we have for all P ∈ Ω h and f ∈ Γ h , and where N * is a positive integer independent of h. In all the following derivations, we assume that the mesh satisfies Assumption 1, and consequently the inequalities in (19). We first introduce some norms and norm equivalence results. For all v h ∈ W h we define and while with · we indicate the standard Euclidean norm.

Lemma 1
We have the following inequalities: for a C > 0, and, for h sufficiently small, Proof The equivalence relations in (25) is an immediate application of the definition (23) and (24), and of the inequalities in (19).
The proof of (26) and (27) for a sufficiently small h.
In the following we will indicate by A the global matrix, It is well known, see for instance [38,57], that A has N f positive and N B = N Γ + N P negative eigenvalues. We will indicate by 0 < λ M c min ≤ λ M c max the minimum and maximum eigenvalue of M c , by 0 < μ B min ≤ μ B max the minimum and maximum singular value of B T , and by λ T max > 0 the maximum eigenvalue of T, which clearly corresponds also to the maximum eigenvalue of T. Finally, we define where π k = [p k , p Γ ,k ] T ∈ Q h is the pressure component of the eigenvector corresponding to the kth negative eigenvalue of A. Clearly ζ ≤ 1.
We have the following Proof The result is an extension of the one given in [55]. For the sake of brevity, we report only the part that differs from the cited reference. We indicate with In particular, we consider here a λ < 0. In this case, M c − λI is non-singular and we may build the Schur complement of block M c , by which we can write The part of the proof that differs from that in [55] concerns the estimate of λ −,A min , where we exploit the special structure of T. More precisely, from (33) we may derive that The proof is concluded by integrating this result with the other bounds provided in the cited reference.
We can now proceed to specialize to our problem the various terms in Lemma 2.

Lemma 3
The eigenvalues of M c satisfy asymptotically, i.e for h sufficiently small, the following bounds

Lemma 4 If Γ h is K -orthogonal on each fracture, then
Proof If the grid is K -orthogonal T is symmetric positive semidefinite and diagonally dominant, and the estimate is then easily obtained by bounding the Rayleigh quotient and by the fact that the fracture equations are posed on d − 1 dimensional surfaces embedded in R d .
Numerical experiments in Sect. 4.3 have verified the validity of this estimate, also in cases where the grid is not strictly K -orthogonal.
We state now the following conjecture.
This conjecture has been verified indirectly by examining the behavior of the spectrum of A experimentally. It can be justified by the fact the pressure component π = (p, p Γ ) of an eigenvector of A is an approximation of the pressure component of an eigenfunction We can now consider π as the pressure eigenvector for which the ratio in (30) reaches its maximal value. Thanks to (25), we can then infer that ζ = p Γ 2 p Γ 2 + p 2 h. However, no rigorous proof is currently available.
To estimate the singular values of B T we extend the work in [50] to the case of fractured domains. To this aim, we have to make an additional assumption on the mesh. It is a technical assumption which is satisfied if the mesh does not exhibit "pathological situations". We define the total grid Λ h = Ω h ∪ Γ h , which contains both bulk and fracture polygonal cells. So a generic cell c ∈ Λ h may be a bulk cell or a fracture cell. Let us consider the undirected graph G where the elements of Λ h are the graph nodes and the set of graph edges L h ⊆ Λ h × Λ h are defined by:

Assumption 3
The global mesh Λ h satisfies the following assumption: there exists a family of elementary paths Ψ = {ζ i } i=1,...,N γ which are connected subgraphs of G with nodes of maximal degree 2, such that -for every cell c ∈ Λ h there exists one and only one path ζ i ∈ Ψ with c ∈ ζ i , i.e. the family of paths Ψ defines a partition of Λ h ; -the first and last node, c s,i and c e,i , of any ζ i ∈ Ψ have a facet in F ∂Ω h ; -there exists a constant L * independent from h such that where L i is the number of cells in ζ i .
We can now state the following.

Lemma 5 For h sufficiently small, the singular values of B T satisfy the following bounds
Since the proof is rather technical and lengthy, we have postponed it to the Appendix.

Theorem 3
Let h be sufficiently small, and assume that Assumptions 1, 3 and Conjecture 1 hold. Then, the spectrum of A satisfies where k i are positive constants independent from h but depending on the bounds on the permeability in the bulk and in the fractures. Consequently, the condition number is characterized by the following bound, 4 Proof The bounds to characterize the intervals I − and I + are obtained by applying estimates (36), (35), (38) and (37) to Lemma 2. Indeed, by considering only the leading terms for h sufficiently small we have and Since Concerning the interval of positive eigenvalues we immediately have and, finally Defining the following constants, 1 C * , k 3 = C * and k 4 = C * , we recover the asymptotic spectrum estimate (39) and we can conclude that by which we obtain the estimate of the condition number (40). 4 We may note that without Conjecture 1 we will have in (34) that for λ −,A min the term ζ λ T max is of order h d−3 . This will imply a condition number K 2 (A) of the order h −4 , much more restrictive.
Concerning the asymptotic behavior of the condition number in function of the problem parameters, we may note that the most relevant term is proportional to Therefore, as expected, an important role is played by the permeability contrast between bulk and fractures. We remark that the results here presented are asymptotic, i.e. they hold for an h sufficiently small. The dependence on the problem parameters shows that the leading term of O(h −3 ) in the condition number is more relevant when fractures are much more permeable than the bulk. This has been confirmed in the numerical experiments.
The next result is important for the derivation of suitable preconditioners for the system. We define M c D = diag(M c ), the diagonal part of the mimetic inner product matrix. We have the following where λ Consequently, see [58], the Rayleigh quotient satisfies where v T P and v P Γ are the restriction of the vector of velocity degrees of freedom v ∈ W h on the velocity degrees of freedom of P and to the velocity degrees of freedom of P that lies on Γ , respectively. We have that v T P M c P v P = v T P M P v P + v T P Γ E P v P Γ , where M P and E P are the elemental contributions to the matrices M and E, respectively. Coercivity and stability results for M c P , reported in [42], together with the quasi-uniformity assumption (19) and the definition of E, lead to Let us now indicate with m c j j the jth diagonal element of M c P and we set m c j j = m j j + e j j , where m j j and e j j are the elements of M D,P and E D,P , the diagonal part of M P and E P respectively, and we indicate with e j the jth canonical vector with all components zero apart from [e j ] j = 1. With simple algebraic manipulations, from (11) and (12) we find that Since 2 , and noting that I−Z P (Z T P Z P ) −1 Z T P is a projection matrix, we can claim that While, considering the definition of E given in (10), and indicating with F Γ h (P) ⊂ F Γ h the set of facets of P that lay on Γ h andf the fracture cell corresponding to a f ∈ F Γ h (P), we may since ξ 0 > 0. In conclusion, by exploiting (49), (50) and (51) by which the quotient is asymptotically a constant.

Preconditioning Techniques for the Discrete Generalized Saddle-Point System
In this section we present several ABF preconditioners applied to the problem at hand. For an overview of preconditioning techniques the reader may refer to [59], and [21] for the particular case of saddle-point problems. For an illustration of iterative methods for problems arising from the discretization of partial differential equations the reader may refer to [38]. The objective is to verify the effectiveness of a diagonal approximation of the matrix M c . Indeed, the result of Theorem 4 suggests that approximating M c with its diagonal, could be an interesting and simple choice.
To be computationally effective, the block preconditioners we are going to consider require also an approximation for the Schur complement S = − T − BM −1 c B T . For Mixed Finite Elements for the Stokes problem there is a wide strand of literature on the topic of finding good approximations of S, see for instance [21]. On the other hand, in the context of MFD discretizations we are not aware of any effective strategy. We propose here the approximation simply given by Matrix − S is s.p.d. and rather sparse, so it is easily factored with a sparse Cholesky factorization [34]. In a first set of test cases we will consider also the approximation of M c with its lumped form a technique widely used in low-order finite elements. We will show that this choice is completely inadequate for our class of problems. We have chosen two iterative schemes, MINRES [52] and GMRES [53]. The former is in principle the method of choice for indefinite symmetric systems. However it restricts the choice of preconditioners to symmetric positive definite ones. Therefore, the reduced constraints on the possible preconditioners of GMRES makes it particularly interesting, despite its larger memory requirement. For GMRES we have chosen the left-preconditioned formulation.

Block Diagonal Preconditioner
The first ABF preconditioner we have considered is given by where S if computed according to (53) and M c using (52) or, alternatively, (54). The application of the preconditioner is straightforward, and recalled in Algorithm 1. At each iteration, we compute the preconditioned residual z = [z 1 , z 2 ] T given by Pz = r, where r = [r 1 , r 2 ] T is the current residual.
This preconditioner is symmetric positive definite, therefore it may be used for both GMRES and MINRES. The most demanding part of Algorithm 1 is solving the linear system − Sz 2 = r 2 , which is however made effective by the Cholesky factorization (a possible alternative, not considered in this work, is to use multigrid techniques).

Block Triangular Preconditioner
In this case, we set Given the above definitions, the solution of Pz = r can be obtained by employing the following factorization of its inverse More precisely, at each iteration, computing the preconditioned residual is attained by Algorithm 2. Again, the most demanding part of Algorithm 2 is solving the linear system Algorithm 2 Block triangular preconditioner. Computation of z = P −1 r.
The additional computational effort compared with the diagonal block preconditioner is minimal. However, this preconditioner is not s.p.d. and therefore it has been tested only with GMRES.

Block LU Preconditioner
Now we introduce a preconditioner based on an approximate block LU factorization. For the matrix A the following block LU factorization holds: The idea is again to replace M c and S with the proposed approximations M c and S, so that the preconditioner is now Multiplying the factors we get which shows that P has the classical form of an indefinite preconditioner. 5 Computing the preconditioned residual z = [z 1 , z 2 ] T = P −1 r requires the steps illustrated in Algorithm 3. The computational costs are slightly higher than the ones of Algorithm 2, with an additional multiplication by M

Numerical Tests
In this section we present some numerical experiments. The first illustrates some solutions obtained with the proposed model on two test cases, one with a single fracture cutting the domain and the second one with a more complex network. We then use the set-up of the first test case to verify the spectral properties of the matrices with respect to the theoretical estimates. Finally, we test the performance of the proposed preconditioners using the matrices stemming from the discretization of the two test cases.  Fig. 2, we vary their interaction. We have considered 3 cases: a "sealed" fracture K τ = 10 −2 , K n = 1, where the low normal permeability hinders the flow across the fracture; a "neutral case" with K τ = 10 −2 , K n = 10 2 , where the material filling the fracture has the same permeability as the bulk and we expect negligible pressure and velocity jumps; and a highly conductive fracture, K τ = 1, K n = 10 2 . Note the jump of pressure across the fracture obtained in the case of the "sealed" fracture, Fig. 2-left, and the linear pressure profile in Fig. 2-center, where the fracture has a negligible effect on the bulk flow. More interesting is the case of the conductive fracture in Fig. 2-right, where we have a strong interaction between the fracture and bulk flow, and the fracture becomes a preferential path for the fluid.

Setup B: Fracture Network
We here consider a more realistic case, i.e. a network of fractures. The bulk domain is now Ω = (0, 2) × (0, 1) × (0, 1) and the network Γ consists of seven fractures with several intersections shown in Fig. 3. On the left and right boundary sides of the domain we consider Dirichlet conditions, in particular we fix p = 1 on {x = 0} and p = 0 on {x = 2}, while on the top, bottom, front and back boundary sides of the domain we impose homogeneous Neumann boundary conditions. A polyhedral mesh of diameter h = 0.193 is employed and the resulting dimension of the system is 43,834. In this test we consider larger contrasts between the equivalent permeability K τ and K n and the bulk permeability with respect to the previous test. In the case of sealed fractures we set K τ = 10 −3 , K n = 10 −1 . As shown in Fig. 4 the action of the fractures as barriers implies a strong discontinuity of the pressure across the fractures. The case of conductive fractures, with K τ = 10, K n = 10 3 is shown in Fig. 5: note that in this case the distribution of pressure in the bulk follows the main network direction, since the fractures attract the bulk flow.

Spectral Properties of the Matrix
We consider the geometrical configuration of Setup A, illustrated in Sect. 4.1. First, we study experimentally the maximum eigenvalue of the transmissibility matrix T to verify (36). Then, we estimate the condition number of the global matrix of the problem A to verify the theoretical estimate of Theorem 3. The maximum and minimum eigenvalues have been estimated with the eigs function of MATLAB ® .
We wanted to assess that the theoretical bound (36) on the maximal eigenvalue of T holds (at least approximately) also for more general grids. Since d = 3 the bound states that λ T max should be bounded by a quantity invariant with h.
In Table 1 we have reported the estimated values for different grids andK Γ = k f I, with k f = 10 −3 . We may note that for regular hexahedral grids, which induce a rather regular mesh on Γ , the eigenvalue is in fact constant with a value that respects the bound very closely.  However, even for more general meshes, the value of λ T max is scarcely affected by h, with higher values that probably reflects the lower mesh regularity.
We focus now on the numerical validation of the estimate of the condition number of the system matrix A stated in Theorem 3.
We can observe in Table 2 that the h-dependence of the condition number changes with the fracture permeability and with h, as expected. Yet, the theoretical bound in (40), with is reached only for the highest value of k f , when the coefficient in (46) has the highest value of 10. In the other cases, for the values of h here considered, a less restrictive variation with h has been found. This is probably due to the fact that the meshes are not refined enough. Indeed, the theoretical estimate is asymptotic, and with a mesh size not sufficiently small the other terms defining the intervals I − and I + in (39) may have a greater impact on the condition number. The case of high fracture permeability is indeed the most challenging because of the strong effect of the coupling terms. These results also confirm Conjecture 1, without which the condition number would be O(h −4 ).

Testing the Preconditioners
We now present tests to assess the performance of the preconditioners illustrated in Sect. 3. To make a fair comparison all tests have been made by constructing the matrices, building a random vector w, constructing the left hand side as b = Aw and solving Ax = b. The stopping tolerance has been set so that on each test the final relative error in the 2-norm, i.e x − w 2 / w 2 withx being the approximated solution, is of the order of 10 −6 . To assess the performance in terms of CPU times, whenever possible we compare with the time needed for the solution with the multifrontal direct solver UMFPACK [35]. UMFPACK is a very efficient tool for solving sparse linear systems, so it provides a valuable reference. It is memory demanding though, so we could not use it for the largest examples. The tests were performed on an 2.7 Ghz i7 Intel processor with 16 GBytes RAM.
For the initial tests we also considered a simple diagonal preconditioner, which corresponds to a rescaling of the equations. Since our matrix A has a diagonal block of zeroes, for the corresponding rows no scaling has been performed, i.e. the block has been replaced with the identity. This comparison has been set up only to show the performance gained by the ABF preconditioners versus a trivial preconditioner.
In the following we will denote with Diag the simple diagonal preconditioner, with ABFD_D, ABFTr_D and ABFLU_D the approximate block factorization preconditioners based on block diagonal, block triangular and the block LU factorization with M c approximated by its diagonal part, respectively. When we use the lumped inner-product matrix to approximate M c , we use the subscript L: ABFD_L, ABFTr_L and ABFLU_L.
The first set of tests still considers a single fracture that cuts the whole domain as described in Sect. 4.1, whereas the second deals with the network of fractures presented in Sect. 4.2, and is more realistic and challenging. For GMRES the restart level is set to 100 in all simulations (a value that is never reached if we choose an effective preconditioner).
The fundamental goal of the analyses is the study of the effectiveness of the preconditioners for different parameter values and mesh sizes. In all tests, the bulk permeability is assumed to be K = I, the fracture aperture l Γ = 10 −2 .

Setup A: Single Fracture Case
Here we present the case with a single fracture that cuts the whole domain with the geometric configuration shown in Fig. 2.
We consider first how the different preconditioners behave with respect to the mesh size h in the case of neutral fracture, i.e. K τ = 10 −2 and K n = 10 2 . The number of iterations to reach convergence, along with the corresponding computing time, are shown in Table 3. Note that solving the system with a trivial diagonal preconditioner is extremely inefficient, as expected, and will not be considered in the next tests. We have contemplated it here only to show that the linear system stemming from our problem cannot be solved in practice by an iterative method without resorting to a good preconditioner. Indeed, ABFD_D, ABFTr_D, ABFLU_D provide a significant improvement in performance with timings comparable to the direct solver, particularly in the most refined cases.
The first thing we can note is that the approximation with inner-product lumping is highly ineffective. Therefore, it will not be considered in the next tests. The good results obtained RelTime is the ratio between the actual computing time and the one required for a global solve with the direct method implemented in UMFPACK: note that a small RelTime indicates a good performance of the iterative solver. X means that the iterative solver has not reached convergence in the prescribed number of iterations (20000). We have used polyhedral meshes with the diagonal approximation are in line with the finding of Theorem 4. Indeed, it is well known that approximate block triangular and approximate LU preconditioners perform well if the (1, 1)-block matrix and the Schur complement are replaced by spectrally equivalent matrices [38,55]. GMRES with ABFLU_D outperforms all the other techniques, with a number of iteration rather low and computational time that outperforms UMFPACK for the larger matrix. With the ABFD_D we obtain a slightly better performance with MINRES than with GMRES, as expected. Since memory is also a possible bottleneck we report in Table 4 the memory requirement of UMFPACK, MINRES and GMRES for the cases of Table 3, limited to the block diagonal preconditioner (memory requirement of GMRES with block triangular or LU preconditioners is of the same order). We may note that GMRES is more demanding in terms of memory (as expected) and the memory requested by the direct solver is much higher and grows approximately quadratically with the number of degrees of freedom.
We now examine the robustness of the preconditioner with respect to the parameters of the fracture model K τ and K n . We consider the permeability values corresponding to the three cases of Fig. 2: K τ = 10 −2 , K n = 1 (sealing fracture); K τ = 10 −2 , K n = 10 2 ("neutral" fracture); K τ = 1, K n = 10 2 (conductive fracture). Hereafter, we do not consider the two preconditioners based on the lumping of M c , nor the Diag preconditioner, which have shown a poor performance. For this test a polyhedral mesh of size h = 0.208 has been employed RelTime is the ratio of the actual computing time to the one of a global solve with UMFPACK  R1 is the ratio between the time taken for the solution on a given grid and GMRES with ABFLU_D on the same grid. R2 indicates the ratio of the time taken on a given grid and that for the same method on the 300K grid and the corresponding dimension of the linear system is of size 43031. The results, presented in Table 5, show a substantial invariance with respect to parameter changes.
To complete this first analysis we considered also more demanding problems in terms of number of degrees of freedom, namely 318,209, 615,793 and 1,216,061, denoted in the following as 300K, 600K and 1.2M respectively. The value of the permeability in the rock matrix is unitary and in the fracture the effective permeabilities are K τ = 10 −2 and K n = 10 2 . For practical reasons, we are using here grids made of tetrahedra. In this case UMFPACK is not able to solve the problem because of memory constraints. Therefore, we report two different relative times. The first (indicated by R1) is the ratio between the time taken on a given grid and the time taken by GMRES+ABFLU_D on the same grid. The second (indicated by R2) is the ratio between the time taken on a given grid and the one for the same method on grid 300K.
The results are reported in Table 6. We may note that the results in terms of number of iterations are rather invariant with respect to the mesh size, showing the good scalability of all preconditioners. Moreover, the numbers of iterations are smaller than the ones reported in Table 5. This is explained by the fact that with tetrahedral grids the M c matrix is much sparser. Indeed, while the number of elements per row is of the order of 26 for the polyhedral grids we have considered, it drops to 7 for the tetrahedral grids. We have not investigated this behavior further, but it is consistent with similar findings for the 2D case reported in [44], even if for a different discretization scheme. The choice GMRES with ABFLU_D is also in this case the one performing best in terms of both number of iterations and computing time. We can also observe for all preconditioners a power dependence of the computing time on the total number of degrees of freedom with an exponent of about 1.6.

Setup B: Fracture Network Case
We consider now the configuration illustrated in Sect. 4.2. The larger contrasts between the equivalent fracture permeabilities K τ and K n and the bulk permeability with respect to the previous set-up and the presence of a full network makes the effect of fractures more important. Also in this case, illustrated in Table 7, we can confirm the good performance of GMRES with ABFLU_D which has proved to be robust also with respect to variations of the model parameter. The other techniques show a slight degradation of performance compared with the results of the previous, simpler, setup.
Finally, to verify the robustness with respect to parameter contrast, we performed two further tests with high contrast. In particular, we have set a unitary permeability in the bulk, whereas in the fractures permeability is 8 orders of magnitude smaller or larger, leading to effective permeabilities of K τ = 10 −11 , K n = 10 −5 for the impermeable case, K τ = 10 5 , K n = 10 11 for the permeable case. We note that for these values, the estimated condition number of the full matrix is approximately 2 × 10 9 and 3 × 10 11 , respectively. We are then dealing with rather ill-conditioned problems. The results, shown in Table 8, show that all the proposed ABF preconditioners behaves reasonably well even for large contrast problems, in particular GMRES with ABFLU_D is particularly robust and consistently better performing than the direct solver. We point out that in the test cases of this section, GMRES with ABFD_D has restarted, and this explains the higher number of iterations compared to MINRES (which does not apply restart).

Conclusions
In this work we have assessed the spectral properties of the linear system arising from the mimetic finite differences discretization of a hybrid-dimensional Darcy problem in a fractured porous media. We have then implemented and tested a set of ABF preconditioning techniques for the iterative solution of the discrete problem, proposing a strategy to build the approximation of the factors. For what concerns the spectral analysis the technique adopted is an extension of that proposed in [50] to take into account the hybrid dimensional nature of the problem. The main finding is that the condition number scales asymptotically as O(h −3 ), a more restrictive result with respect to mimetic finite differences applied to the standard Darcy equation. The reason lays in the hybrid-dimensional nature of the problem, in particular the presence of the coupling term in the mimetic inner product matrix M c . To reach this result we had to make an assumption to control the contribution to the spectral estimate coming from the term deriving from the discretization of the problem in the fractures. The conjecture is based on the fact that this term operates only on the 1-codimensional manifold representing the fracture network. We note that this situation is different from that of the standard stabilised Stokes problem, where the (2, 2) block in the Stokes matrix is the discretization of an operator acting on the whole domain. Further work will be needed to prove the conjecture, but numerical experiments seem to confirm it.
We have investigated numerically the effect of different grid size and type, and different values of contrast between bulk and fracture permeability. We have found that rather classic block-diagonal and block-triangular and LU ABF preconditioners, where the mimetic inner product matrix M c is replaced by its diagonal and the same technique is adopted to construct the approximation of the Schur complement, work extremely well, showing -a good scalability with respect to the mesh size h, despite the non-favorable scaling of the condition number, -and robustness also with respect to variations of the problem parameters, in particular GMRES with ABFLU_D.
Further analysis may study the effect of heterogeneity and anisotropy in the permeability tensor. Another line of research would be to test our approximation strategy with the preconditioners for double saddle-point problems recently proposed in [5].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

A Proof of Lemma 5
Proof Any singular values B T is the square root of an eigenvalue of which is a full rank matrix, thanks to the discrete inf-sup condition, as proven in [42]. Let us consider each block more closely in order to characterize the elements of the matrix in detail. The block (1, 1) is BB T ∈ R N P ×N P and we have Now we focus on block (1, 2), i.e. BC T ∈ R N P ×N Γ . We have Block (2, 1) is just the transpose of the block (1, 2), while block (2, 2), given by CC T ∈ R N Γ ×N Γ , is diagonal, with elements We start by seeking an upper bound for the eigenvalues λ( B B T ) of B B T . Thanks to Gershgorin's Theorem, we have where R i are the row circles. From (60)-(62), we note that we have two types of row circles.
Indicating by x i c and r i the center and the radius of circle R i , respectively, we have since N * ≥ 2. The estimate of the smallest eigenvalue requires the additional mesh Assumption 3. To do so we first consider the elements of the vector B T π ∈ R N f . For each f ∈ F 0 h we can identify the two adjacent cells P f 1 and P f 2 so that α P f 1 ,f = 1 and α P f 2 ,f = −1. For f ∈ F ∂Ω h we can identify the cell P f of which f is a boundary facet and for f ∈ F Γ h also the correspondingf ∈ Γ h , i.e the fracture cell that coincides geometrically with f. So, for 1 ≤ i ≤ N f and indicating with f i the ith bulk facet, we have If we now consider the Rayleigh quotient of B B T we may rearrange the terms to get We can characterize the smallest eigenvalue as We now order the elements of Λ h path by path. Therefore, π = [p, p Γ ] T may be partitioned in [π 1 , . . . , π N ζ i ] T , where π i = [π c s,i , . . . , π c e,i ] T are the pressure values (either in the bulk or in the fracture) associated to the cells in path ζ i . For each graph edge e of ζ i we indicate with π e 1,i and π e 2,i the elements of π i associated to the cells at end of the edge. Since |f| h d−1 , we have The last term is equivalent to the Rayleigh coefficient of the block diagonal matrix where E i = tridiag(−1, 2, −1) ∈ R L i ×Li . The eigenvalues of E i can be computed explicitly as Therefore, λ min (E i ) = 2 1 − cos 1 L i + 1 , and, from Assumption 3 on the maximum length of the paths, Finally, For sufficiently small h it holds that 1 − cos h L * +h ≤ h 2 2L 2 * , by which we can conclude that