Block Preconditioners for Mixed-dimensional Discretization of Flow in Fractured Porous Media

In this paper, we are interested in an efficient numerical method for the mixed-dimensional approach to modeling single-phase flow in fractured porous media. The model introduces fractures and their intersections as lower-dimensional structures, and the mortar variable is used for flow coupling between the matrix and fractures. We consider a stable mixed finite element discretization of the problem, which results in a parameter-dependent linear system. For this, we develop block preconditioners based on the well-posedness of the discretization choice. The preconditioned iterative method demonstrates robustness with regards to discretization and physical parameters. The analytical results are verified on several examples of fracture network configurations, and notable results in reduction of number of iterations and computational time are obtained.

age and production. It has become clear that the dominating role of fractures in the flow process in the porous medium calls for reexamination of existing mathematical models, numerical methods and implementations in these cases.
Considering modeling and analysis, a popular and effective development is reduced fracture models [9,18,22] that represent fractures and fracture intersections as lower-dimensional manifolds embedded in a porous medium domain. The immediate advantages of such modeling are in more accurate representation of flow patterns, especially in case of highly conductive fractures, and easier handling of discontinuities over the interfaces. This has also allowed for implementation of various discretization methods, from finite volume methods [22,30] to (mixed) finite element methods [18] and other methods [17,19]. These methods mostly differ in two aspects: whether the fractures conform to the discrete grid of the porous medium [9] or are placed arbitrarily within the grid [12,16,31], or whether pressure or flux continuity is preserved. Comparison studies of different discretization methods and their properties can be found in [5,15,28].
Although there is a wide spectrum of discretization methods, little has been done to develop robust and efficient solvers. This aspect of implementation can be very important since applications of fractured porous media usually include large-scale simulations of subsurface reservoirs and the resulting discretized linear systems of equations can become ill-conditioned and quite difficult to solve. The linear system represents a discrete version of the partial differential equation (PDE) operator that has unbounded spectrum. Thus, its condition number tends to infinity when the mesh size is approaching zero. Moreover, the variability of the physical parameters, such as the permeabilities and aperture, can additionally influence the scale of the condition number of the system. Instead of using direct methods, we consider Krylov subspace iterative methods to solve such large scale problems. Since the convergence rate of the Krylov subspace methods depends on the condition number of the system, suitable preconditioning techniques are usually required to achieve a good performance. A recent study on a geometric multigrid method [4] for the fracture problem shows how standard iterative methods can be extended and perform well on mixed-dimensional discretizations, but still there are limitations that need to be overcome for general fractured porous media simulations.
In this paper, we aim to provide a general approach to preconditioning the mixed-dimensional flow problems based on suitable mixed finite element method discretization developed in [9]. Beside introducing the mixeddimensional geometry, the main aspects of the discretization are flux coupling between subdomains using a mortar variable and inf-sup stability of the associated saddle-point problem. Moreover, this framework has been shown to be well incorporated within functional analysis as a concept of mixed-dimensional partial differential equations [8], allowing even further applications in poroelasticity and transport problems.
We propose a set of block preconditioners for Krylov subspace methods for solving the linear system of equations arising from the chosen discretization. Following the theory in [26] and [25], we derive uniform block precondition-ers based on the well-posedness of an alternative but equivalent formulation. Proper weighted norm is chosen so that the well-posedness constants are robust with respect to the physical and discretization parameters but depend on the shape regularity of the meshes. Both block diagonal and triangular preconditioners are developed based on the framework [25,26]. Those block preconditioners are not only theoretically robust and effective but can also be implemented straightforwardly by taking advantage of the block structure of the problem.
The rest of the paper is organized as follows. In Section 2 we first introduce the mixed-dimensional geometry and the governing equations of the single-phase flow in fractured porous media followed by the variational formulation and the stable mixed finite element discretization of the problem. The framework of the block preconditioners is briefly recalled in Section 3 and its application to mixed-dimensional discretization of flow in fractured porous media is proposed and analyzed in Section 4. We verify the theoretical results by testing several numerical examples in Section 5 and finalize the paper with concluding remarks in Section 6.

Preliminaries
In this section, we set up the problem of flow in fractured porous media following [9]. Let Ω n be a domain of the porous medium of dimension n " 2, 3 that can be decomposed by fractures into Ω n i , i P I n . The fractures and their intersections are represented as lower d-dimensional manifolds Ω d i , i P I d , 0 ď d ă n, and inherit the similar decomposition structure as the porous medium Ω n (see Figure 1). Here, we use I d as a local index set in dimension 0 ď d ď n. Furthermore, we define Γ d ij for j P J d i Ď I d as interfaces between Ω d`1 i and adjacent Ω d j . Union over the subscript set I d represents all d-dimensional subdomains, that is Finally, the fractured porous medium domain Ω with interface Γ is defined as Remark 2.1 . Even though the theoretical results in [8,9] allow for a more complex geometrical structure, for the sake of simplicity we restrict the model to domains of rectangular type. That is, we approximate fractures as lines on a plane for n " 2 or flat surfaces in a box for n " 3. However, we allow for any configuration of fractures or fracture intersections within, for example, very acute angles of fracture intersections, multiple intersecting fractures or T-type intersections.
Now that we have set up the dimensional decomposition framework for the fractured porous medium, we introduce the governing laws in the subdomains and fractures. First, notation and properties of the physical parameters are introduced. For the sake of simplicity, we slightly abuse the notation by omitting subdomain subscripts and dimension superscripts in the following definitions. We only keep the indices in certain cases when clarification is necessary.
Assume that the boundary of Ω can be partitioned to BΩ " BΩ D Y BΩ N such that BΩ D X BΩ N " H and BΩ D is of positive measure. We adopt the notation in each dimension 0 ă d ď n , that is The material permeability K and normal permeability K ν tensors are considered to be bounded both above and below, symmetric and positive definite, and we denote ν as an outward unit normal on Ω. Furthermore, let γ d ij be the distance from Γ d ij to Ω d i , which for d " n´1 represents the fracture aperture. The physical parameters K and γ may vary spatially. However, to simplify the analysis, we assume that they are constant on each subdomain in each dimension.
In each Ω d , we introduce the governing Darcy's law and mass conservation, find fluid velocity u d and pressure p d that satisfy where we introduce an additional mortar variable λ d , defined as to account for the mass transfer across each interface Γ d ij , and a jump operator ¨ : L 2 pΓ d q Ñ L 2 pΩ d q as Since there is no notion of interface Γ n or flow in a point Ω 0 , we extend the definion of λ n and u 0 by setting them equal to zero. An additional interface law on Γ d ij is introduced to describe the normal flow due to the difference in pressure from Ω d i to Ω d`1 , Finally, proper boundary conditions are needed. For example, (2.10) Remark 2.2 . In the previous equations, we have used u d as integrated flux and p d as averaged pressure in each Ω d , 0 ď d ď n. Therefore, the scaling with the cross-sectional area ε of order Opγ n´d q due to the model reduction has been accounted for within the permeability parameters K and K ν .

Variational formulation
Now we consider the variational form of the problem (2.5)- (2.10). For any open bounded set ω P R n , let L 2 pωq and H s pωq denote the L 2 space and the standard Sobolev spaces on functions defined on ω, respectively. Also, denote H 1 2 pBωq as the space of L 2 -traces on the boundary Bω of functions in H 1 pωq. Let p¨,¨q ω be the L 2 -inner product and }¨} L 2 pωq the induced L 2 -norm. We define where V d representing the flux function space on Ω d , Q d the pressure space on Ω d , and Λ d the function space of normal flux across interface Γ d . Furthermore, let V d 0 be a subspace of V d containing functions v 0 such that v 0¨ν " 0 on Γ d´1 . In addition, define the extension operator R d : To summarize the formulation, we compose function spaces over dimensions 12) and associate composite L 2 -inner products and induced composite L 2 -norms The system (2.5)-(2.10) in the weak formulation reads: Find pu 0 , λ, pq P V 0ˆΛˆQ that satisfies with g P H 1 2 pBΩ D q and f P L 2 pΩq. As before, functions u 0 0 , v 0 0 , λ n and µ n are set to zero.
We end this section by observing the saddle point structure of the system (2.13). First, let W " V 0ˆΛ be the function space of all flux variables, including mortar variable, and define the mixed-dimensional divergence operator D¨: W Ñ Q as D¨w " D¨ru 0 , λs " ∇¨u 0` λ , w P W . (2.14) Define the following two bilinear forms Then the saddle point form of system (2.13) reads: Find pw, pq P WˆQ such that It has been shown in [9] that the bilinear forms ap¨,¨q and bp¨,¨q are continuous with respect to the following norms for r " rv 0 , µs P W and q P Q, In addition, ap¨,¨q is shown to be coercive on the kernel of bp¨,¨q in [9] as well. Finally, the following theorem states that bp¨,¨q satisfies the inf-sup condition. Theorem 2.1 [9]. Let the bilinear form bp¨,¨q be defined as in (2.15b). Then there exists a constant β ą 0 independent of the physical parameters K, K ν and γ such that Following the classical Brezzi theory [6,11], we conclude that the saddle point system (2.16) is well-posed, i.e., there exists a unique solution of (2.16).

Discretization
We continue this section with discretizing the problem (2.16) by the mixed finite element approximation. Let T d Ω and T d Γ denote a d-dimensional shaperegular triangulation of Ω d and Γ d , and h " max 0ďdďn h d the characteristic mesh where RT 0 stands for lowest-order Raviart-Thomas(-Nédélec) spaces [27,29] and P 0 for the space of piecewise constants. Furthermore, define Then we can define the discrete extension operator R d h : Analogous to the continuous case, we define the discrete composite spaces and the linear operators With W h " V 0hˆΛh , the finite element approximation of the system (2.13) is formulated as follows: Find pw h , p h q P W hˆQh such that, Due to our choice of the finite element spaces, the continuity of ap¨,¨q and bp¨,¨q and the coercivity of ap¨,¨q on the kernel of bp¨,¨q are preserved naturally. To show the well-posedness of the discrete saddle point system (2.22), we need the inf-sup condition to hold on the discrete spaces as well. This has been shown in [9] and is stated in the following theorem. Theorem 2.2 [9]. There exists a constant β ą 0 independent of the discretization parameter h and the physical parameters K, K ν and γ such that Therefore, the finite element method (2.22) is well-posed by the Brezzi theory [6,11].
We finalize this section with the block formulation of the discrete saddle point system (2.22 The well-posedness of the system (2.22) ensures that A is an isomorphism from W hˆQh to its dual W 1 hˆQ 1 h and, therefore, (2.24) has a unique solution pw h , p h q P W hˆQh .

Block preconditioners
In this section, we briefly present the general preconditioning theory for designing block preconditioners of Krylov subspace iterative methods [25,26], which introduces necessary tools for the analysis in the following section.
The block preconditioning framework [25,26] is based on the well-posedness theory. Therefore, we first introduce the setup of the problem. Let X be a real separable Hilbert space and p¨,¨q X represent the inner product on X that induces the norm }¨} X . Furthermore, denote X 1 as a dual space to X and x¨,¨y as the duality pairing between them. Let Lp¨,¨q be a bilinear form on X that satisfies the continuity condition and the inf-sup condition, for α, β ą 0. We aim to construct a robust preconditioner for the linear system where A : X Ñ X 1 is induced by the bilinear form Lp¨,¨q such that xAx, yy " Lpx, yq. The properties of the bilinear form ensure that A is a bounded and symmetric linear operator and the system (3.2) is well-posed. Our goal is to develop block preconditioners for solving (3.2).

Norm-equivalent Preconditioner
Consider a symmetric positive definite operator M : X 1 Ñ X which induces an inner product px, yq M´1 :" xM´1x, yy on X and corresponding norm }x} 2 M´1 :" px, xq M´1 . Naturally, MA : X Ñ X is symmetric with respect to p¨,¨q M´1 and we can use M as a preconditioner for the MINRES algorithm whose convergence rate is stated in the following theorem.
Theorem 3.1 [20]. Let x m be the m-th iteration of the MINRES method preconditioned with M and x be the exact solution, it follows that κpMAq`1 and κpMAq denotes the condition number of MA. As shown in [26], if (3.1) holds and M satisfies, then A and M are called norm-equivalent and κpMAq ď c2α c1β . Thus, if the well-posedness constants α and β and the norm-equivalence constants c 1 and c 2 are all independent of the physical and discretization parameters, then M provides a robust preconditioner.
One natural choice of the norm-equivalent preconditioner is the Riesz operator B : X 1 Ñ X corresponding to the inner product p¨,¨q X pBf , xq X " xf , xy, @f P X 1 , x P X.
It is easy to see that if we choose M " B, then (3.3) holds with constants c 1 " c 2 " 1 and, therefore, the preconditioned system has a bounded condition number If the constants α and β are independent of the discretization and physical parameters, we obtain a robust preconditioner.

Field-of-values-equivalent Preconditioner
In this section, we recall the class of field-of-values-equivalent (FOV-equivalent) preconditioners which allow more general preconditioners than the norm-equivalent ones. Consider a general operator M L : X 1 Ñ X which can be used as a preconditioner for the GMRES method. The following theorem, developed in [13,14], characterizes the convergence rate of the GMRES method.
Theorem 3.2 [13,14]. Let x m be the m-th iteration of the GMRES method preconditioner with M L and x be the exact solution, it follows that where, for any x P X, M L is referred to as an FOV-equivalent preconditioner if the constants Σ and Υ are independent of the physical and discretization parameters. Usually M L provides a uniform left preconditioner for GMRES.
In a similar manner, we can introduce a right preconditioner for GMRES, M U : X 1 Ñ X and consider the preconditioned system By introducing an inner product on X 1 , defined as px 1 , y 1 q M :" xx 1 , My 1 y, we say M U and A are FOV-equivalent if, for any where the constants Σ and Υ are independent of the physical and discretization parameters. Therefore, M U can be used as a uniform right preconditioner for GMRES. In many cases [1,2,25], the FOV-equivalent preconditioners can be derived based on the Riesz operator and the FOV-equivalence can be shown based on the well-posedness conditions (3.1).

Robust Preconditioners for Mixed-dimensional Model
In this section, we design block preconditioners based on the general framework mentioned in the previous section. Consider the finite element approximation (2.22). In this case, define X " W hˆQh associated with the following norm Then, the operator A : X Ñ X 1 in (2.24) is induced by the bilinear form and satisfies the well-posedness conditions (3.1) due to Theorem 2.2, the continuity of the bilinear forms ap¨,¨q and bp¨,¨q, and the coercivity of ap¨,¨q on the kernal of bp¨,¨q. Moreover, the constants α and β are independent of parameters h, K, K ν and γ. The Riesz operator corresponding to the norm }¨} X in (4.1) is where A and B are defined as in (2.24) and I p is the identity operator on Q, i.e., xI p q h , q h y " }q h } 2 Q . The main challenge in implementation of this preconditioner is to solve for the upper block A`B T B that corresponds to I`grad div problem. One way of resolving this is to use auxiliary space theory (see for example [21,24]). However, in our case, additional theory resulting from the mixed-dimensional exterior calculus in [8] is needed, which is the topic of our ongoing work [7]. However, in this paper, we consider an alternative formulation of the problem (2.22) and show the well-posedness with respect to a different weighted norm, which allows for a simpler robust preconditioner.

An Alternative Formulation
In order to introduce the alternative formulation, we need to define a discrete gradient operator D h : Q h Ñ W h such that, for r h " rv 0h , µ h s, where, for w h " ru 0h , λ h s and r h " rv 0h , µ h s, , . -.
Here T d P T d Ω is either a tetrahedron for d " 3, a triangle for d " 2 or a line segment for d " 1. Furthermore, f d P BT corresponds to a face of the element T d , ν f d denotes the unit outer normal of face f d , and φ f d P RT 0 pT d q is the basis function on face f d . Using the discrete gradient operator, an alternative formulation of the system (2.22) is given as follows: Find pw h , p h q P W hˆQh such that, The well-posedness of the alternative formulation (4.5) with respect to the norm (4.1) follows directly from the well-posedness of the original formulation (2.22) because the two formulations are equivalent. However, in order to derive a block preconditioner different from (4.3), we shall consider the same coefficient operator A (2.24) with a different weak interpretation and the wellposedness in a different setting.
The alternative weighted norm we consider for the alternative formulation (4.5) is defined as where }r h } 2 a :" apr h , r h q and }r h } 2 a D :" a D pr h , r h q. In order to show (4.5) (or the operator form (2.24)) is well-posed with respect to this alternative norm (4.6), we need the following two lemmas. The first lemma shows that the forms ap¨,¨q and a D p¨,¨q are spectrally equivalent.
Lemma 4.1 . There exist constants c 1 , c 2 ą 0, depending only on the shape regularity of the mesh T Ω , such that the following inequalities hold, Proof. Recall that }r h } 2 a " aprv 0h , µ h s, rv 0h , µ h sq " pK´1pv 0h`Rh µ h q, pv 0h`Rh µ h qq Ω`p γK´1 ν µ h , µ h q Γ , Obviously, (4.7) holds if pK´1pv 0h`Rh µ h q, pv 0h`Rh µ h qq Ω and pK´1pv 0hR h µ h q, pv 0h`Rh µ h qq D,Ω are spectrally equivalent. Note that Therefore, we can immediately observe that it is enough to show that pK´1pv 0hR h µ h q, v 0h`Rh µ h q T d and pK´1pv 0h`Rh µ h q, v 0h`Rh µ h q D,T d are spectrally equivalent on each element T d , 0 ă d ď n. In addition, by using the scaling argument [10, Section 4.5.2], we only need to show they are spectrally equivalent on a reference elementT d , i.e., We show the proof for d " n " 3. For other cases the proof follows similarly. For d " n " 3, the reference elementT d is a tetrahedron with vertices p0, 0, 0q, p1, 0, 0q, p0, 1, 0q and p0, 0, 1q in the Cartesian coordinates. The local matrix AT d , representing pK´1pv 0h`Rh µ h q, v 0h`Rh µ h qT d , takes the following form By the definition, pK´1pv 0h`Rh µ h q, v 0h`Rh µ h q D,T d is represented by the diagonal of AT d , which we denote as D AT d " K´1 120 diagp18, 16, 16, 16q. To show (4.8) onT d , it is enough to notice that, under our assumption that K is constant on each T d , the generalized eigenvalue problem AT d y " χD AT d y gives all eigenvalues χ ą 0 independent of physical and discretization parameters. Therefore, (4.8) holds withc 1 " ? χ min andc 2 " ? χ max , where χ min and χ max denote the smallest and largest eigenvalue, respectively. The spectral equivalent result (4.7) follows directly by the scaling argument [10, Section 4.5.2] and summing over all T d P T d Ω , 0 ď d ď n. The constants c 1 and c 2 depend on the shape regularity of the mesh due to the scaling argument but do not depend on the physical and discretization parameters.
[ \ Based on the spectral equivalence Lemma (4.1), we have the following infsup condition regarding the discrete gradient D h .

Lemma 4.2 . Let the discrete gradient operator D h be defined as in (4.4).
Then there exists a constant β ‹ ą 0 independent of the discretization and physical parameters such that Proof. Using Lemma 4.1, we have for any Now the result follows taking infimum over all q h P Q h and β ‹ " c´1 2 .
[ \ Based on Lemma 4.1 and 4.2, by Babuska-Brezzi theory [6,11], we can conclude that the alternative formulation (4.5) is well-posed with respect to the norm (4.6) as stated in the following theorem.
It satisfies the continuity condition and the inf-sup condition with respect to |||pr h , q h q|||, i.e., for any pw h , p h q P W hˆQh and pr h , q h q P W hˆQh , with constants α and β dependent on the shape regularity of the mesh but independent of discretization and physical parameters.

Block diagonal preconditioners
The well-posedness Theorem 4.3 provides alternative block preconditioners for solving the linear system (2.24) effectively. To this end, we introduce a linear operators D A : W h Ñ W 1 h which is defined as xD A w h , r h y " a D pw h , r h q for w h , r h P W h . The reason we use the notation D A here is that, by the definitions of ap¨,¨q and a D p¨,¨q, the matrix representation of linear operator D A is exactly the diagonal of the matrix representation of linear operator A. Then, by the definition of the discrete gradient operator D h (4.4), we have D A D h " B T and, therefore, for q h P Q h . Based on the above operator form of the }¨} a D norm, the Riesz operator corresponding to the |||¨||| norm (4.6) is (4.10) As discussed in Section 3.1, B D is a norm-equivalent preconditioner for solving the system (2.24) and we have the following theorem regarding the condition number of B D A. In practice, applying the preconditioner B D implies inverting the diagonal block exactly, which can be expensive and sometimes infeasible. Thus, we consider the following preconditioner where the diagonal blocks H w and H p are symmetric positive definite and spectrally equivalent to diagonal blocks in A and BD´1 A B T , respectively, i.e.
where the constants c 1,w , c 1,p , c 2,w , and c 2,p are independent of discretization and physical parameters. In practice, H w can be defined by a diagonal scaling, i.e., D´1 A and H p can be defined by standard multigrid methods. In general, the choice of H w and H p are not very restrictive, provided it handles possible heterogeneity in physical parameters K, K ν , and γ. M D is a norm-equivalent preconditioner as well. Following [26], we can directly estimate the condition number of M D A in the following theorem.
Theorem 4.5 . Let M D be as in (4.11) and let (4.12) hold. Then it follows that κpM D Aq ď αc2 βc1 , where c 2 " maxtc 2,w , c 2,p u and c 1 " mintc 1,w , c 1,p u. Remark 4.2 . Again, κpM D Aq is bounded independently of h and parameters K, K ν and γ, but remains dependent on the shape regularity of the mesh.

Block triangular preconditioners
In this subsection, we consider the block triangular preconditioners based on the FOV-equivalent preconditioners we discussed in Section 3.2. Here, we analyze the robustness of block triangular preconditioners and show the corresponding FOV-equivalence, which leads to uniform convergence rate of the GMRES method.
The block lower triangular preconditioners take the following form On the other hand, the block upper triangular preconditioners are given as (4.14) Basically, M L and M U are inexact versions of B L and B U when the diagonal blocks are replaced by spectrally equivalent approximations (4.12). Next theorem shows that B L and A are FOV-equivalent.
Theorem 4.6 . There exist constants ξ 1 , ξ 2 ą 0, independent of discretization and physical parameters, such that for every x " pw h , p h q P W hˆQh , x ‰ 0, Proof. By the definition of the linear operators A and D A , we naturally . On the other hand, again using the Cauchy-Schwarz inequality and equivalence of the norms }¨} D´1 A and }¨} A´1 we get for each y " pr h , q h q P W hˆQh , y " 0 with ξ 2 " max 2, 2c´1 1 ( , which concludes the proof. [ \ The next theorem states that if the conditions (4.12) hold then M L and A are FOV-equivalent.
Using the same conditions to show the upper bound, we obtain This gives the upper bound with ξ 2 " maxt2c´1 1,w , c´1 1,p c´1 1 p1`c´1 1,w qu, which concludes the proof.
[ \ Remark 4.3 . Due to Lemma 4.1, the constants ξ 1 and ξ 2 are independent of h and parameters K, K ν and γ, but remain dependent on the shape regularity of the mesh. This means that the convergence rate of the preconditioned GMRES method with preconditioner B L or M L depends only on the shape regularity of the mesh.
Similarly, we can derive the FOV-equivalence of B U and M U with A. Since the proofs are similar to the two previous theorems, we omit them and only state the results here. Theorem 4.8 . There exist constants ξ 1 , ξ 2 ą 0 independent of discretization and physical parameters such that for any x 1 " B´1 U x with x " pw h , p h q P W hˆQh , x ‰ 0, Theorem 4.9 . If the conditions (4.12) hold and }I´H w A} A ď ρ for 0 ď ρ ă 1, then there exist constants ξ 1 , ξ 2 ą 0 independent of discretization and physical parameters such that for any Remark 4.4 . Similarly, the constants ξ 1 and ξ 2 here are independent of h and parameters K, K ν and γ, but remain dependent on the shape regularity of the mesh. This means that the convergence rate of the preconditioned GMRES method with preconditioner B U or M U depends only on the shape regularity of the mesh.

Numerical results
In this section, we propose several test cases to verify the theory on the robustness of the preconditioners derived above. Both two and three dimensional examples emphasize common challenges in fracture flow simulations such as large aspect ratios of rock and fractures, complex fracture network structures and high heterogeneity in the permeability fields.
In each example below, a set of mixed-dimensional simplicial grids is generated on rock and fracture subdomains, where the coupling between the rock and fracture is employed by a separate mortar grid. Since our main objective is to show the robustness of our preconditioners for standard Krylov iterative methods, for the sake of simplicity, we take the mortar grid to be matching with the adjacent subdomain grids. However, the theory in Section 4 shows no restrictions to relative grid resolution between the rock, fracture and mortar grids. Furthermore, in [28] the discrete system remains well-posed with varying coarsening/refinement ratio for non-degenerate (normal) permeability values, which is one of our assumptions. Therefore, we expect that our block preconditioners give similar performance for general grids between the rock, fracture, and coupling part.
To solve the system (2.24), we use a Flexible Generalized Minimal Residual (FGMRES) method as an outer iterative solver, with the tolerance for the relative residual set to 10´6. The block preconditioners designed in Section 4 are used to accelerate the convergence rate of FGMRES. Each preconditioner B D , B L and B U requires inversion of the diagonal blocks corresponding to flux and pressure degrees of freedom, while the spectrally equivalent versions M D , M L and M U approximate the inverses with appropriate iterative methods. For that, we implement both exact and inexact inner solvers. Solving each diagonal blocks exactly means we use the GMRES method with a relative residual tolerance set to 10´1 0 , while in the inexact case it is set to 10´3. Inner GMRES is preconditioned with unsmoothed aggregation Algebraic Multigrid method (AMG) in a W-cycle.
For obtaining the mixed-dimensional geometry and discretization, we use the PorePy library [23], an open-source simulation tool for fractured and deformable porous media written in Python. Our preconditioners are implemented in HAZMATH library [3], a finite element solver library written in C, also where all solving computations are performed. The numerical tests were performed on a workstation with an 8-core 3GHz Intel Xeon "Sandy Bridge" CPU and 256 GB of RAM.

Example: two-dimensional Geiger network
In the first example, we consider the test case presented in the benchmark study [15]. The domain Ω " p0, 1q 2 , depicted in Figure 2, has unitary permeability K " I for the rock matrix and it is divided into 10 sub-domains by a set of fractures with aperture γ. In our case, we set the tangential and normal permeability of the fractures to be constant throughout the whole network, and vary the value from blocking to conducting the flow. The tangential fracture permeability is denoted as K f to avoid confusion with the rock permeability. At the boundary, we impose zero flux condition on the top and bottom, unitary pressure on the right, and flux equal to´1 on the left. The boundary conditions are applied to both the rock matrix and the fracture network. The numerical solution to this problem is also illustrated in Figure 2. Our goal is to investigate the robustness of the block preconditioners with respect to discretization parameter h and physical parameters γ, K f and K ν . To this end, we generate a series of tests in which we vary the magnitude of one of the parameters, while setting others to a fixed value. This also tests the heterogeneity ratios between the porous medium and the fractures, since we keep spatial and physical parameters of the porous medium unitary. We compute and compare number of iterations of the outer solver for both exact and inexact implementations of the proposed preconditioners. This way we clearly see if the stability of the proposed preconditioners depends on one or a combination of given parameters.
The results of these robustness tests on are summarized in Tables 1 -3. We start with setting K f " K ν " I that, together with rock permeability K, gives a global homogeneous unitary permeability field. We also fix the aperture to γ " 10´2. Refining the initial coarse grid by a factor of 2 recursively, Table  1 demonstrates the robustness of all block preconditioners with respect to the mesh size h. Additionally, the different implementations of the preconditioners result in similar behavior of the solver. We notice that the block triangular 1{4  20  13  12  19  10  10  1{8  19  13  11  19  10  10  1{16  19  13  11  19  10  10  1{32  19  13  11  19  10  10  1{64  19  13  11  19 10 10 Table 1: Number of iterations of outer FGMRES solver with exact and inexact block preconditioners for the case study in Example 5.1. Varying mesh size h while aperture is set to γ " 1{100 and all the permeabilities are set to K " K f " Kν " I.  preconditioners B L and B U show a slightly better performance compared the block diagonal B D as expected. The same behavior can be observed for inexact preconditioners M L and M U in comparison to M D . This is expected since the block triangular preconditioners better approximate the inverse of the stiffness matrix in (2.24). It is noteworthy to mention that the action of the block triangular preconditioners is more expensive computationally than the action of the block diagonal preconditioners. Similar performance can also be observed in Table 2, where we scale down the fracture width on a fixed grid of mesh size h " 1{16. Lastly, in Table 3 we test the influence of the heterogeneity in the permeability fields. We keep the mesh size to be h " 1{16 and fracture aperture to be γ " 10´2, while introducing both conducting and blocking fracture network in the porous medium. Again, the robustness is evident in terms of the number of outer FGMRES iterations with both exact and inexact block preconditioners. The block triangular preconditioners, B L , B U , M L , and M U , provide somewhat lower values comparing to their block diagonal counterpart.

Example: two-dimensional complex network
This example is chosen to demonstrate the robustness of the block preconditioners on a more realistic fracture network. Such a complex fracture configuration often occurs in geological rock simulations and the geometrical and physical properties of the fracture network can significantly influence the stability of the solving method. This is especially seen in mpartitioning the frac-imposed on top and bottom, with pressure 1013250 Pa on the left and 0 Pa on the right boundary. For the comparison with the previous example, we refine the mesh size h with respect to the width of the domain L " 600. However, due to the complex fracture structure, it is possible to end up with smaller and badly shaped elements in the rock matrix grid around the tips and intersections of the fractures. For example, see Figure 4 on the left. The coarser the mesh is, the more irregular the elements are, especially when partisioning in between many tightly packed fractures. Therefore, we expect that the solver requires more iterations to converge on coarser meshes. This is evident in the table on the right in Figure 4. We see the reduction of number of iterations when refining the mesh in all the cases, with the lowest number required by the block upper triangular M U . We also notice that the solver manages to provide the correct solution on all given meshes in an acceptable number of iterations. The results are slightly worse than the previous example, but keep in mind that the complex geometry is still an important factor in the mesh structure and, therefore, influences the convergence rate since the shape regularity of the mesh deteriorates. For complex fracture networks, it is beneficial to invest in constructing a more regular mesh of the fractured porous medium and then applying the proposed block preconditioners in the iterative solvers.

Example: three-dimensional Geiger network
This last example considers the simulations of a 3D problem taken from another benchmark study [5], a three-dimensional analogue to the test case in Subsection 5.1. The geometry is extended to the unit cube and the fracture network now consists of nine intersecting planes (see Figure 5). As before, we take the rock matrix permeability K to be the identity tensor, while we vary the tangential K f and the normal K ν permeability, as well as the fracture aperture γ. For a fair comparison with the two-dimensional case, we perform similar robustness tests of the preconditioners to study the effect of mesh refinement, as well as permeability and aperture changes. However, we stick to only inexact preconditioners M D , M L and M U since they are less computationally expensive and perform comparably well, which makes them good choices in practice. The results are presented in Tables 4-6. We can see that the simulations confirm the findings of Section 4: all block preconditioners show robustness with respect to the discretization and physical parameters. The block diagonal preconditioner requires a slightly higher number of iterations to converge compared to block triangular ones, as we saw in the previous example.  1{4  26  18  15  1{8  26  17  15  1{16  24  16  14  1{32  24  16  13  1{64  24  16  12   Table 4: Number of outer iterations of FGMRES solver with exact and inexact block preconditioners for the case study in Example 5.2. Varying mesh size h while aperture is set to γ " 1{100 and all permeabilities are set to K " K f " Kν " I.
In 3D simulations it is also important to study the overall computational complexity of the solving method. For that, we analyze in Figure 6 the required CPU time of the FGMRES solver preconditioned with each block preconditioner M D , M L and M U . All preconditioners show a optimal OpN dof q complexity, where N dof is the number of degrees of freedom of the discretized system. Notice that even though the block triangular pair of preconditioners  1  24  16  14  1{10  24  16  13  1{100  24  16  14  1{1000  26  16  14  1{10000 26 17 14 Table 5: Number of outer iterations of FGMRES solver with exact and inexact block preconditioners for the case study in Example 5.3. Varying aperture γ while mesh size is set to h " 1{16 and all permeabilities are set to K " K f " Kν " I.
require solving a denser system, it is still time-wise less expensive due to a lower number of iterations needed to converge.

Conclusions
We have presented block preconditioners for linear systems arising in mixeddimensional modeling of single-phase flow in fractured porous media. Our approach is based on the stability theory of the mixed finite element discretization of the model which we extended to provide an efficient way to solve large systems with standard Krylov subspace iterative methods. We have thoroughly analyzed the robustness of the derived preconditioners with regard to discretization and physical parameters by proving norm and field-of-value equivalence to the original system. Our theory has also been supported by several numerical examples of 2D and 3D flow simulations. It is noteworthy to mention that even though our analysis depends on a more regular mesh, the numerical results show that the preconditioners still perform well since the mixed-dimensional discretization approach handles fractures independently of the rock matrix and, therefore, generates simpler meshes in most fracture network cases. The large aspect ratios that parametrize the model then become the main stability problem, which we have successfully overcome with the proposed block preconditioners. This is impor- We conclude by recalling that the alternative approach to block preconditioners mentioned in the beginning of Section 3 is a non-trivial extension to this work and a part of an ongoing research.