Unified approach to discretization of flow in fractured porous media

In this paper, we introduce a mortar-based approach to discretizing flow in fractured porous media, which we term the mixed-dimensional flux coupling scheme. Our formulation is agnostic to the discretizations used to discretize the fluid flow equations in the porous medium and in the fractures, and as such it represents a unified approach to integrated fractured geometries into any existing discretization framework. In particular, several existing discretization approaches for fractured porous media can be seen as special instances of the approach proposed herein. We provide an abstract stability theory for our approach, which provides explicit guidance into the grids used to discretize the fractures and the porous medium, as dependent on discretization methods chosen for the respective domains. The theoretical results are sustained by numerical examples, wherein we utilize our framework to simulate flow in 2D and 3D fractured media using control volume methods (both two-point and multi-point flux), Lagrangian finite element methods, mixed finite element methods, and virtual element methods. As expected, regardless of the ambient methods chosen, our approach leads to stable and convergent discretizations for the fractured problems considered, within the limits of the discretization schemes.


Introduction
Flow in porous media with thin inclusions is an important process both within subsurface and industrial materials. Our main focus herein is on the subsurface, where the thin inclusions represent fractures, and the fracture space can be either open or filled. We will thus simply refer to fractured porous media in what follows. However, thin inclusions may also be engineered in artificial porous media for the purpose of fluid flow control.
Fluid flow in fractured porous media is a dominating process in several subsurface applications, ranging from geothermal energy production, shale gas recovery and nuclear waste deposits. As such, accurate and reliable numerical representations have been an important topic of research, and Rainer Helmig has been a major contributor to the field for more than three decades. Existing discrete representations of fractured porous media fall in two major categories, depending on whether the fractures conform to the underlying discrete grid representing the porous materials. So-called "unfitted" discretizations, wherein the fractures are allowed to be arbitrary with respect to the grid, have seen significant research and developments in recent years (see e.g. [1,2]). Our focus herein is, in contrast, on discretizations where the discrete grid resolves the fractures, which are conceptually simpler than unfitted discretizations.
Early research into numerical simulation and conforming discretization of fractured porous media was spear-headed by among others Rainer Helmig and his collaborators [3]. This early work was centered around lowest-order finite element discretizations. Later, it was understood that local conservation properties were important for discretization methods for flow in porous media, and conforming discretizations of fractured porous media were developed based on control volume approaches [4,5], mixed finite element methods [6,7], mimetic finite differences [8] and virtual element methods [9]. See also [10] for a comparison study.
A recent development in the mathematical representation of fractured porous media is the modeling and interpretation of fractures as lower-dimensional manifolds [11,7,12]. This concept allows for the introduction of mixed-dimensional partial differential equations (md-PDEs), wherein partial differential equations are defined, in a coupled sense, both in the porous material, lower-dimensional fractures, and yet lower-dimensional intersections. In this abstraction, it can be shown that the mathematical models for fractured porous media, can be cast in a rich functional-analysis framework, ensuring wellposedness, and thus existence and uniqueness, of solutions [13].
In this manuscript, we revisit conforming discretizations of fractured porous media within the context of md-PDEs. We show, by introducing explicit coupling variables in the spirit of mortar methods [14,11,7,15], an abstract framework for constructing a conforming fracture discretization from any discretization of non-fractured porous media. We term this approach the mixed-dimensional flux-coupling (MDFC) method. Viewed from the discretization within each dimension, the coupling between dimensions takes the form of standard boundary value problems, thus any implementation that can account for Dirichlet and Neumann boundary data can be applied to fractured media with minimal adaptations. Our approach thus unifies the various previous developments reviewed above.
We concretize the abstract framework by applying it to well-known discretizations from literature, establishing (in some cases for the first time) that these discretizations are well-posed. To illustrate the versatility of the framework, we provide numerical examples showing how five different discretization methods for non-fractured porous media can be applied as discretization methods for fractured porous media. Of these discretizations, when using mixed finite elements or standard finite elements for the non-fractured media, we recover earlier methods referenced above. In the case of finite volume (both two-point and multi-point flux) and virtual element methods, our approach effectively leads to a discretization scheme not previously discussed in literature. Our numerical examples, which include a 2D case where we use non-matching grids between the dimensions and a relatively complex 3D case, highlight the convergence properties and stability of MDFC even for degenerating parameters.
The remaining manuscript honors the following structure: In section 2, we introduce our novel approach to unifying discretization methods for fractured media. Thereafter, in section 3, we show the stability of the approach theoretically, which emphasizes the conditions required between the (in principle non-matching) grids discretizing the matrix and fractures. Numerical examples and verification are presented before concluding the paper.

Modeling fractured porous media
In this section we introduce our model for fractured media, first by a single fracture, and then extended to general fracture networks.

Domain with a single fracture
Flow in (fractured) porous media can lead to complex and non-linear governing equations. However, at the heart usually lies a second-order partial differential equation, which upon linearization (i.e. within a Newton iteration) thus takes the classical form for a pressure 3 and flux 3 Here we denote by Ω 3 the (3-dimensional) porous medium, and by and its Neumann and Dirichlet boundaries, respectively. We denote by Ω ± 2 Ω 3 the boundary of Ω 3 as seen from the positive (resp. negative) side of Ω 2 , and the outer normal vector is always denoted . The Dirichlet boundary data is set to zero for notational convenience. We emphasize the structure of the governing equations as composed of a conservation law (2.1), and a constitutive (Darcy) law (2.2). In equations (2.1-2.6) we have marked variables by a superscript '3' to emphasize that they belong in 3 dimensions, the necessity of the precision will be clear below. Note that the flux from the (2-dimensional) Neumann boundary is denoted by a superscript '2'. Throughout the manuscript, we will use to denote right-hand sides, which with the chosen sign convention represents fluid extraction. Similarly, we may consider a single fracture as a (2-dimensional) manifold Ω 2 , whereon the governing equations can in the linearized case be expressed as [16] In equations (2.7-2.8), we denote by a double-strike the tensor operating tangentially (parallel) to the manifold and emphasize that the differential operators are 2-D by a subscript. We note that in equation (2.7), two extra terms arise. These represent the outflow from the fracture into the porous medium on the two sides of the fracture (denoted + and -). As above, fracture variables are indicated by a superscript '2' for clarity.
Considering still the case of a single fracture, equations (2.1-2.10) lead to a system of equations where 2 is a variable internal to the system. We thus complete the model with a constitutive law for 2 , which takes the Darcy-like form (see e.g. [7]) We remark that the within-fracture permeability || and the transverse permeability ⊥ may in practice scale with the aperture and its inverse, respectively.
Equations (2.1-2.11) form a closed and well-posed system of equations for a porous medium including a fracture (see e.g. [8]). More generally, we note that we write these equations in a unified way, in that for = {2,3} Equations (2.12-2.17) make sense with the convention that since there is no 4-dimensional domain in the model, the terms 3 = 0 and || 3 = .
From physical considerations, it is customary to consider all boundaries of the fracture as Neumann boundaries with = 0, except where the boundary coincides with an outer boundary of the full domain. However, these restrictions are not necessary from a mathematical or numerical perspective, and we will retain the slightly more general formulation in order to avoid extra notation for distinguishing between internal and external boundaries of fractures.

Extension to general fracture configurations
Equations (2.12-2.17) are written in a way that naturally generalizes also to fracture intersections, both the 1-D line intersections as well as the 0-D point intersections of three fractures [17,6]. We introduce some extra notation to this end. Let each domain (matrix, fracture, or intersection) be indexed by number and dimension, i.e. Ω is domain number ∈ , having dimensionality . We consider a total of subdomains of various dimensionality. This subdivision is illustrated in Figure 2.   [7,1,6] and references therein). These equations have been identified as a second-order system of mixed-dimensional partial differential equations, for which existence and uniqueness theory has been developed under fairly mild assumptions on the geometry [13]. In this work we will only consider planar fractures, but with no restrictions on their intersections or interaction with the boundary.
In order to simplify notation in the following, we consider the dimension associated with each subdomain, = ( ), to be specified, and introduce the compound variables =

Variational formulation
By shifting indexes on the trace term in (2.24), we identify the symmetric and coupling terms as For non-degenerate coefficients, equations (2.24-2.25) are well-posed by standard saddle-point theory [18], and in the remaining manuscript we will only consider this case. Nevertheless, we remark that, following similar arguments as exposed in [6], it can be shown that significant degeneracy of coefficients can be permitted, at the cost of introducing weighted spaces. In particular, it is of interest to also allow for fractures where the tangential permeability is negligible. Equations (2.24-2.25) are well-posed in this sense, since if for a given domain Ω , the permeability can degenerate in the sense of ,|| → 0, as long as ,⊥ remains bounded from below for all ∈̂. However, now the pressure is only in 2 due to the inf-sup condition for ( , ) [6]. This implies that this weakly continuous formulation for fractured porous media is robust both for arbitrarily thin fractures and can also be applied to blocking fractures. We summarize the above discussion as follows: Let an 2 -like norm on ℋ 1 × ℒ 2 be defined as Furthermore, let the set of indexes be refined such that ∈ if ,|| > 0 and ∈ if ,|| = 0. Then we introduce space as Note here that we use a circle above the function space to indicate homogeneous Dirichlet boundary conditions. Then the equations for flowing and blocking fractures can be written as find ( , ℷ) ∈ × ℒ 2 such that The solution of (2.30) is characterized by the following Lemma.

Proof.
For the two cases in the proof for and , respectively, we indicate variables in these domains by similar subscripts. Then formally, equations (2.26) take the form Here, Δ represents the 1 bilinear forms on Ω , ⊥ represents the 2 bilinear forms om Ω , while Σ are the duality pairings in (2.27). The upper-left 3x3 system is coercive due to the conditions of the proof. Furthermore, we obtain the well-posedness of the full system, since it is easy to show that the Σ terms are inf-sup stable between 2 spaces, indeed Since one may simply choose , = . The coercivity of the upper left 3x3 system together with inf-sup for the Σ terms is sufficient for stability of the full system by abstract saddle-point theory [18]. □

Remark 2.2
Lemma 2.1 is not optimal in the sense that it is fairly easy to extract 1 regularities on all domains ∈ , and the restrictions on ⊥ can be somewhat relaxed. However, as we are primarily interested in the numerical implementation in this contribution, we have chosen to keep Lemma 2.1 as simple as possible. Readers interested in the functional analysis for equations of this type are referred to the papers referenced in the introduction.
It is important to note that the main objective of exposing the equations for flow in fractured porous media on the form (2.26-2.27), is that it highlights the specific domain-decomposition like structure of the problem. Indeed, we note that on each sub-domain (be it porous media, fracture, or fracture intersections), we have a fairly standard elliptic partial-differential equation. These are coupled via interface variables, , . This structure is key to design general and flexible discretization approaches, as introduced in the next section.

Discretizations for fractured porous media
Our exposition of the mathematical model for fractured porous media emphasizes two main aspects of the model, namely the second-order elliptic PDE within each domain, and the flux-coupling terms. Numerous discretization methods have been constructed for second-order elliptic differential equations -many of these are bespoke to the particular challenges associated with flow in highly heterogeneous porous media (for an introduction, see the books [19,20,21]). Herein, we will prove that any stable discretization for flow in (fixed-dimensional) porous media can be applied to fractured porous media through the framework introduced in the preceding section.
We subdivide this section in three parts, in order to provide the mixed-dimensional flux coupling (MDFC) discretization framework, its abstract analysis, and a concrete example using finite elements.
To be precise, we consider each domain Ω and its Neumann boundary Γ = Ω ∪ ∈̌Ω as endowed with a numerical discretization (note that Γ includes all boundaries to lower-dimensional manifolds). We will only consider linear discretizations, however the approach should be applicable also to non-linear discretizations (for a recent contribution in this direction from Helmig's group, see [22]). We do not require that a discrete grid be defined, however we let the discrete representation of 2 (Ω ) and 2 (Γ ) be denoted as ℎ (Ω ) and ℎ (Γ ), respectively. For domains ∈ , i.e. where the fractures are permeable with ,|| ≥ 0,|| , the solution operator of the numerical discretization of the heterogeneous elliptic equation on a given domain ∈ can be stated as . This solution operator maps sinks and Neumann data to pressures and pressure traces, as made precise below. Here, we recall that we for notational simplicity only consider homogeneous boundary conditions on the Dirichlet boundaries, and as such suppress the Dirichlet boundary data. For domains ∈ , the solution operator is void, as there is no differential equation on these domains.
We will use the natural requirement that the numerical discretizations provided are consistent approximations in the following sense: Let ∈ , and let [ , ] = ( , ), for ( , ) ∈ ℎ (Ω ) × ℎ (Γ ), then this quadruplet of variables approximates the solution to the elliptic differential equation The precise interpretation of ≈ will depend on the chosen numerical method. We note that standard methods such as finite volume, finite element, mixed-finite element and spectral methods all fall within this framework, where the approximation implied by the ≈ signs of equations (3.1-3-4) can for most numerical methods be characterized by grid regularity, material parameters, grid resolution, etc. By assumption, we consider only stable numerical methods, in the sense of a negative eigenvalue-spectrum for the numerical solution operators , with potentially a single degenerate eigenvalue for subdomains where Ω = Ø, and we will denote the smallest (i.e. most negative) nondegenerate eigenvalue of as − . Furthermore, the system (3.1-3.4) is self-adjoint, so that in many cases the numerical method will be symmetric (see Section 3.3 below for the case of finite elements).

MDFC: A unified discretization of fractured porous media
To provide a discretization for fractured systems, a grid is introduced on the lower-dimensional manifolds Ω on which the boundary flux variables , will be defined. We emphasize that this mortarlike grid can be chosen independently of any grid used by the numerical methods and ̂+ 1 , thus we impose a minimum of restrictions on the grids. Nevertheless, note that this construction ensures that the flux variables on either side of a fracture (or either sides of fracture intersections) are conforming with each other. The precise relationships between the admissible grids , as implied by the numerical methods , will be made clear below. For the sake of symmetry, we also define grids for the Neumann data on Ω .
To formulate discrete methods for fractured porous media, we represent the flux variable as piecewise constant on the mortar grid , thus , ∈ 0 ( ) and ∈ 0 ( ) (higher-order approximations are also possible, but the regularity of the problem does not seem to justify this). We introduce projection operators in order to move between the degrees of freedom of the numerical methods and the mortar grids . We first define the compound operator projecting from the coupling variables on the mortar grids to the subdomain degrees of freedom and conversely from the numerical variables to the coupling variables Now, our MDFC discretization framework for fractured porous media takes the form: For given numerical discretizations : Find , ∈ 0 ( ), for all ∈ and ∈̂ such that subject to the discrete constraints: The dummy variables and have the interpretations of sinks and fluxes due to the interactions with other domains, respectively. In contrast, the variables and are the pressure and pressure traces after projection onto the grids . The variable is the pressure trace projected onto the Neumann boundaries, and is not used with the boundary conditions considered herein (but would be relevant with Robin boundary conditions).
This MDFC scheme has a particularly simple interpretation: For each subdomain ∈ , can be interpreted as a generalized Neumann-Dirichlet map, in the sense that it maps boundary fluxes (which also take the apparent form of sources for neighboring domains of − 1) to Dirichlet data (where conversely, for < , the internal values are considered Dirichlet data for neighboring domains of dimension + 1). As such, equation (3.8) resolves the internal differential equations in each subdomain, equations (3.9) is the projection of variables from the flux grids to the numerical boundary (and source) data, while equation (3.7) simply states that the flux , between a fracture and its surroundings should satisfy a form of Darcy's law, depending on the difference in pressure of the fracture and the pressure at the boundary of the surroundings. Equations (3.7-3.9) are thus a Schurcomplement formulation of the discrete problem.

Abstract analysis
Let the discretization methods corresponding to the solution operators be collected in a linear system, i.e. we state equation (3.8) on the form: (3.10) Similarly, we denote the compound projection operators Π and Π . Furthermore, denote by the discrete divergence operators from equation (3.9), which sums flux variables associated with a fracture while retaining Neumann boundary data i.e.
Finally, let the diagonal mass matrix associated with the inner product (( ,⊥ ) −1 , ) appearing in equation (3.9) be denoted κ −1 . Then we can eliminate the subdomain variables from the discrete system (3.7-3.9) to obtain a Schur-complement system only in terms of the flux variables, i.e.
From the Schur complement form, we immediately obtain the following result:

Lemma 3.1
Let all subdomain discretization methods be negative definite for ∈ (i.e. Ω ≠ Ø for all ∈ ), and furthermore let the assumptions of Lemma 2.1 hold. Then if the projection operators are negative transposes, such that Π T = −Π , the Schur-complement system (3.12) is stable, with no degenerate eigenvalues.
Proof: By the choice of , ∈ 0 ( ), the κ −1 matrix is diagonal, and has positive eigenvalues bounded below by 0,⊥ −1 . Thus, it is sufficient to show that the remaining term has non-negative eigenvalues. But since is negative definite by the assumption of the lemma, then T Π Π = −(Π ) T Π will be non-negative definite. The result follows since the right-hand side operator is bounded by the assumption of the Lemma. □ In order to allow for fractures (and intersections, etc.) which do not have a Dirichlet boundary, the arguments of Lemma 3.1 must be refined. To this end, let ̅ be the subset of which do not have a Dirichlet boundary. For these domains, we have a pure Neumann problem, and equations (3.8) are expected to constrain the solutions up to a constant (pressure). For the analysis, we therefore introduce an auxiliary constant pressure ̅ for each domain ∈ ̅ , and introduce the modified numerical methods ̃∶ [ ℎ (Ω ), ℎ (Γ )] ∖ ℝ → [ ℎ (Ω ), ℎ (Γ )] ∖ ℝ, i.e., the solution corresponding to equations (3.1-3.4) with a compatibility condition (fluxes and sinks must sum to zero), and the additional constraint that the pressure has mean value zero. For ≠ ∖ ̅ , the solution operator is unaltered, ̃= .

Lemma 3.2
Let all subdomain discretization methods ̃ be negative definite for ∈ , and furthermore let the assumptions of Lemma 2.1 hold. Furthermore, let ∖ ̅ contain at least one domain. Then if the projection operators are negative transposes, such that Π T = −Π , the saddle-point system (3.14-3.15) is stable, with no degenerate eigenvalues. This result is obtained by considering (all) such that ∈ ∖ ̅ . Construct a rooted tree(s) from spanning all subdomains (this can always be done for connected domains). Then for leaves (i.e. terminal nodes of the tree) we set , = ̅ , where is the parent of (we use the sign convention that , = − , if is in ̌, and it is sufficient to consider , constant). Proceeding in this manner recursively, let be a node in the tree and let , be determined for all branches extending from . Then set , = ̅ − ∑ , . Proceeding until the root of the tree, we see by construction that (Π ) ̅ , ℷ = ̅, and that ‖ℷ‖ ≤ ‖̅‖, where increases with the depth of the tree(s) . For a finite geometry, is therefore bounded by the geometry of the fracture network, and independent of the discretization methods. The solvability and bounded eigenvalues of (3.14-3.15) then follows from standard theory [18]. □ In practice, it is of course also of interest to obtain values for the discrete solutions , and not only the flux exchanges ℷ. This result is slightly more subtle, in a similar sense as Lemma 2.1. To prepare, we write equation (3.12) in the same form as used in the proof of Lemma 2.1.
Here, the linear operators are the inverses of , and represent the linear discretizations underlying the numerical solution. Hence, equation (3.17) is also structurally similar to the natural implementation of the methodology. It is also important to note that the form (3.17) is agnostic to whether a domain is in ̅ , thus from the perspective of implementation, it will in many cases not be necessary to introduce special treatment of these domains as in Lemma 3.2. We now obtain a similar result as for the continuous case, in the sense that Proof: The proof is identical to Lemma 2.1 in the continuous case. □ We make the following remarks regarding Theorem 3.3 and its implications for MDFC: 1. All standard numerical methods for elliptic partial differential equations will satisfy condition c) in the theorem, thus essentially any numerical method can be applied to fractured porous media through the MDFC approach given in Section 3.1. 2. There are no restrictions on the grids in relation to the numerical methods as long as the fracture permeabilities ,|| do not degenerate. In particular, for grid-based numerical methods , non-matching grids, both coarser and finer, can be used between the external domain and , and furthermore into the internal domain. 3. In practice, conditions c) and d) of the theorem state that for subdomains where ,|| degenerates, the discrete representation of must not be finer than , . This is similar to the typical conditions encountered in traditional mortar methods [15]. 4. In the special case where is chosen as the mixed-finite element method, analysis shows that spatially degenerating ,|| can be allowed, thus circumventing the binary structure of Lemma 2.1 and Theorem 3.3 [6].

Corollary 3.4
A sequence of grids { } , numerical methods { } and projection operators {Π} , where increasing is understood to enumerate finer grids, will be a convergent approximation to equation (2.30), provided the approximations to equations (2.18-2.23) are consistent.
Proof: Since the problem is linear, stability and consistency are sufficient for convergence. □

Worked example: Finite element methods
In order to make the presentation more concrete, we consider the finite element method with continuous linear Lagrange elements in the framework presented above. Thus, for each Ω let be the corresponding grid, with nodal degrees of freedom.
Then for ∈ , the elements of the sub-matrices of are simply given by the inner products of , ∈ 1 ( ) with Neumann data implemented as natural boundary conditions through the duality pairing The Neumann boundary conditions are exactly dual to the evaluation of traces, and thus the operator will be self-adjoint. Standard finite element theory further guarantees that the required bound on the eigenvalues of holds independent of grid spacing with [23] ≤ ( ,|| ) −1 (3.19) Since the solution and its trace live in finite-dimensional subspaces of 2 , the projection operators become defined in the standard way, i.e. for , ∈ 0 ( ) the projection Π , ∈ 1 ( ) satisfies It is therefore clear that Π T = Π . Thus, all the conditions of Theorem 3.3 are satisfied, provided that the grids are no finer than whenever ,|| → 0.
We note that the finite element approximation could also be obtained directly from Section 2 by simply using the finite-dimensional spaces and the bilinear forms defined in equations (2.26-2.27). Thus equations (3.7-3.9) with the numerical methods defined by equations (3.14-3.15) and projection operators defined by equation (3.15) is equivalent to the symmetric and bilinear saddle-point problem: Find ( ,ℎ , ,ℎ ) ∈ 1 ( ) × 0 ( ) such that ( ,ℎ , ,ℎ , , ) + ′ ( ,ℎ , ) + ′ ( , ,ℎ ) = ( , ) for all ( , ) ∈ 1 ( ) × 0 ( ) (3.21) This discretization is consistent within each domain (for shape-regular grids), thus it represents a consistent approximation to equations (2.30) whenever the boundary data is resolved. Note that for matching grids between the mortar space and the finite element spaces, this does not hold, since a "checkerboard"-type oscillation in 2 is projected to zero by Π . Thus while the lowest-order finite element variant of MDFC is stable for matching grids, it requires that the grids are coarser than the grids chosen for resolving the elliptic partial differential equations in order to be a convergent numerical discretization for equations (2.18-2.23).
While the approach as stated above is sufficient, in the sense of obtaining a stable and convergent discretization, we also remark that an improved method would likely be obtained by honoring the structure of from section 2.3, and thus using ,ℎ ∈ 0 ( ) for ∈ . In particular, this would eliminate the projection errors associated with the low-permeable fractures. This highlights the flexibility of the framework to accommodate different discretizations in the different domains, bespoke to the physical processes.

Example calculations
To confirm the theory derived above, we propose two synthetic test cases in which the ambient space is two-and three-dimensional, respectively. Out of the range of numerical methods to which the MDFC applies, we consider five discretization schemes, summarized below.
Two mixed methods are employed, namely the mixed finite element (RT0), and the dual virtual element method (VEM). The mixed finite element, considered in [6], is given by Raviart-Thomas elements of lowest order for the fluxes and piecewise constants for the pressure in all dimensions. On the other hand, VEM [9] employs a single degree of freedom per face for the fluxes without explicitly specifying the basis functions and represents pressures as piecewise constants. Thirdly, employing nodal-based, linear Lagrange elements in all dimensions leads to the primal formulation (P1) as presented in Section 3.3. This is the only method considered in this work which does not respect local mass conservation. Finally, two finite volume methods are considered, the two-point flux approximation (TPFA) and the multi-point flux approximation scheme (MPFA) [24].
In line with the spirit of the theory presented in this work, the coupling between dimensions employs a flux mortar variable, defined as piecewise constants on a separately generated, lower-dimensional grid. All computations are performed using the open-source simulation tool PorePy [25,26].

Two-dimensional fracture system
The first example, obtained from [6], consists of a unit square with five one-dimensional fractures as given in figure 2. Immersed in the top half of the domain are two intersecting, conductive fractures with permeability ⊥ = 10 4 and || = 1. Note that due to the dimensionless scaling, this corresponds to fractures that are equally conductive in the parallel direction (in terms of volume per unit pressure drop) to the full porous unit square domain. Below are two half-immersed blocking fractures ( ⊥ = 1, || = 10 −4 ) and finally, a conductive fracture separates the lower right corner. The boundary conditions are chosen as a unit pressure drop from top to bottom and no-flow conditions on the sides. The matrix permeability is set to 1.
This example is designed to contain all the elements that constitute challenges for numerical methods for fractured porous media: The two intersecting fractures represent both 1D and 0D domains which have no contact with the boundary, thus the numerical methods on these domains will contain a degenerate eigenvalue (i.e. the pressure solutions are only defined up to a constant). Moreover, the low-permeable and horizontal fractures are expected to lead to singularities in the solution in the 2D domain. Finally, in the lower corner there is a domain which intersects both a Dirichlet and a Neumann boundary.

Figure 3:
The contour lines and color scale of the reference solution on the domain given in Figure 2. The different qualitative aspects of the solution between the conductive and blocking fractures can be clearly seen.
In terms of mesh generation, the one-dimensional fracture grids match the trace of the adjacent twodimensional grids. The mortar grid is then constructed at each fracture to have approximately 75% of the number of elements compared to the inner, lower-dimensional mesh.
Qualitatively, all numerical methods produce the same pressure distributions. Aside from artifacts due to the coarseness of the grid, all methods produce solutions which are visually indistinguishable from the figure 3. We turn to a more quantitative measure in order to expose differences between the discretizations. Since the only common property between the methods is the mortar variable, we compute its 2 -error with respect to a fine-scale solution obtained using the RT0 method. In case of convergence, the rate will be limited to first order with respect to the mesh size, since the mortar variable is represented by piecewise constants.
The results of this convergence test are shown in figure 4. For the one-dimensional mortar variables, very similar behavior is observed for the methods RT0, VEM, and MPFA, exhibiting stable and linear convergence. The two remaining methods show lower than first-order convergence on average. For P1, we speculate that this is due to its lack of local mass conservation, since the error is measured in a flux variable. For TPFA, this deviation is likely due to the lack of consistency in the method (i.e. the approximation error to equations (3.1-3.4) does not necessarily go to zero with grid size). We emphasize that all methods are robust and stable from a linear algebra perspective on all grids. The error in the mortar variables defined at the zero-dimensional intersection is analyzed in figure 4 (right). These results are slightly more sporadic since an accumulation of errors can occur from the higher dimensions, and since this essentially represents a point evaluation of the solution. Moreover, the grids used in the computations are not nested and mesh sensitivities of the method with may be the cause of these effects. Nevertheless, the overall trend in all methods is a decrease in error as the mesh becomes finer. The mixed finite element methods exhibit a more monotone decay in comparison to the finite volume methods, likely due to the fact that the reference solution is calculated with the RT0 method on a finer grid.

Stability
It is of interest to verify the claims of Theorem 3.3. In particular, we wish to address whether the discrete representation leads to a linear system which has a lower bound on condition numbers, which is independent of grid resolutions for non-degenerate parameters, and allows for degenerate parameters in the sense of conditions a)-d) in the proof. We have chosen the condition number of the Schur complement system (3.12) as a proxy for the stability of the method, arguing (as in the preceding section) that the condition number of the full system will depend on the particular features of the numerical methods and grids utilized outside of the fractures to an extent where it is difficult to make a fair comparison. In order to emphasize grids and parameters, we simplify the example from Section 4.1 by omitting the fractures which do not touch the boundary, and replacing the no-flow boundary conditions on the sides of the domain by a linear pressure variation. We can then consider Theorem 3.1 purely in terms of the mortar variables , . Furthermore, in order to reduce the parameter space, we will let the remaining three fractures have the same parameters ⊥ and || .
We fix the grid in the 2D domain with a resolution corresponding to the second-coarsest grid (approximately 4.5k triangles) in the convergence test of Section 4.1. Then in addition to the two fracture parameters, we introduce two grid parameters: The relative resolution of the outer grid to the mortar grid, and the relative resolution from the mortar grid to the fracture grid.
Our aim is to see how the lowest eigenvalue of the discrete Schur-complement system (3.12) depends on the fracture parameters and grid parameters. To this end, we have conducted a suite of simulations for all methods, exploring the full 4D parameter space. We observe that the results are completely independent of || and the ratio of the mortar grid to the inner grid. When varying the perpendicular permeability ⊥ , the results depend primarily on whether the mortar grid is finer or coarser than the outer grid, and weakly depends on the ratio. These results are summarized in Figure 5. methods are also stable for coarse mortar grids for large values of ⊥ . This result reflects the fact that for coarse mortar grids, the Neumann-Dirichlet maps stabilize the system, and that numerically there is an inf-sup condition on Π such that T Π Π has a lowest eigenvalue. We note however, that this does not hold for the continuous system given in equation (2.30), since the trace spaces for the pressure are not rich enough to control the mortar space. This explains why stability is lost on fine mortar grids for all methods, and is also reflected for the P1 variant of the method, which has a worse stability constant for high tangential permeability even for matching grids (see also discussion in section 3.3).. Thus, in all cases and for all grids, the MDFC method is stable, with eigenvalue bounded from below by the continuous problem.
Based on these computations, we summarize that for non-degenerate parameters, all discretizations lead to stable systems for the mortar variable, independent of grid resolution between matrix, fluxvariable, and the fractures. For degenerate fracture flow || , all methods remain stable. Finally, for degenerate fracture cross-flow ⊥ , the results are in accordance with Theorem 3.3. Finally, we consider simulations in a 3D problem. The computational domain is taken as the unit cube, and the fracture network for this example is reported in Figure 6 (left). The latter consist of 9 fractures with a structure similar to the Benchmark 1 in [Flemisch2017], extended to 3D. The matrix permeability is the identity tensor. We introduce the scaling factor = 10 −4(3− ) , for each lower dimensional object the normal permeability is given by ⊥ = 10 4 / and the tangential by || = 10 4 . Flow is forced diagonally across the domain by specifying a pressure value of 1 at boundaries characterized by ( , , ) < 0.4, and similarly a pressure of −1 at boundaries with ( , , ) > 0.8. On all other boundaries, no-flow conditions are assigned. For illustration, the numerical solution computed using RT0 is reported in Figure 5 (right).

Three-dimensional Example
To compare the numerical schemes, we investigate numerical convergence of the mortar variables in the same way as in Section 4.1. Three simplex grids are considered, with cell counts of about 3.5k, 4.5k and 10k tetrahedrals, together with a suitable number of triangles, line elements and points. For simplicity, we consider only matching grids in this case. Since P1 is not convergent for matching grids (see Section 3.3. and the discussion in Section 4.2), we exclude this variant of MDFC from our results. Errors in the mortar variables are computed relative to a reference solution obtained with RT0 on a grid with about 190k tetrahedral cells. The resulting error decay is depicted in Figure 7. The simulation confirms the findings in section 4.1: MPFA, RT0 and VEM all exhibit at least first order convergence for all dimensions, while TPFA again suffers from lack of consistency on the ambient grid, thus the low accuracy of the numerical method pollutes the flux variable.

Conclusions
We have developed a new, unified, approach to discretizing fractured porous media, termed Mixed-Dimensional Flux Coupling. The MDFC approach allows for arbitrary numerical discretizations to be used both for the porous media and the fractures. We have supported the development by both theoretical analysis, as well as numerical examples using five different numerical methods.
Several of the limitations included in this work appear to be possible to overcome. In particular, we expect that extension to non-linear discretizations [22] to be straight-forward in practice. Moreover, due to being agnostic of the numerical methods used, our theoretical results are not optimal nor exhaustive, and a more explicit treatment of the precise characteristics of the numerical methods chosen for the various components of the problem is known to provide more nuanced results [6].
In applications, coupled problems are of particular interest. In particular, the fluid flow is often coupled to transport of either mass or energy. Preliminary work in this direction is ongoing, and we expect that the MDFC framework proposed herein will accommodate such coupled problems. We conclude by noting the importance of open-source code availability. The methods developed herein have been implemented in PorePy, and both methods and the scripts used to generate the presented results are available in the public domain at time of publication [26].