Unified analysis of discontinuous Galerkin approximations of flows in fractured porous media on polygonal and polyhedral grids

1.   Introduction

Many geophysical and engineering applications, including, for example, fluid-structure interaction, crack and wave propagation problems, and flow in fractured porous media, are characterized by a strong complexity of the physical domain, possibly involving thousands of fault/fractures, heterogeneous media, moving geometries/interfaces and complex topographies. Whenever classical Finite-Element-based approaches are employed to discretize the underlying differential model, the process of mesh generation can be the bottleneck of the whole simulation, as classical finite elements only support computational grids composed by tetrahedral/hexahedral/prismatic elements. To overcome this limitation, in the last decade a wide strand of literature focused on the design of numerical methods that support computational meshes composed of general polygonal and polyhedral (polytopic, for short) elements. In the conforming setting, we mention for example the composite finite element method ^[1,2], the mimetic finite difference method ^{[3,4,5,6,7,8]}, the polygonal finite element method ^[9], the extended finite element method ^[10,11,12], the virtual element method ^{[13,14,15,16,17]} and the hybrid high-order method ^{[18,19,20,21]}. In the setting of non-conforming/discontinuos polygonal methods, we mention, for example, hybridizable discontinuous Galerkin methods ^{[22,23,24,25]}, composite discontinuous finite element methods ^[26,27], non-conforming VEM ^[28,29,30], gradient schemes ^[31] and the polytopic discontinuous Galerkin method ^{[32,33,34,35,36,37,38,39,40,41,42,43,44]}.

Within this framework, we focus our attention on the problem of modelling the flow in a fractured porous medium, which is fundamental in many energy or environmental engineering applications, such as tracing oil migration, isolation of radioactive waste, groundwater contamination, etc. Fractures are regions of the porous medium that are typically characterized both by a different porous structure and by a very small width. The first feature implies that fractures have a very strong impact on the flow, since they can possibly act as barriers or as preferential paths for the fluid. The second feature entails the need for a very large number of elements for the discretization of the fracture layer and, consequently, a high computational cost. For this reason, one popular modelling choice consists in a reduction strategy, so that fractures are treated as $(d-1)$-dimensional interfaces between $d$-dimensional porous matrices, $d = 2, 3$. In particular, we refer to the model for single-phase flow developed in ^{[45,46,47,48]}. Here, the flow in the porous medium (bulk) is assumed to be governed by Darcy's law and a suitable reduced version of the law is formulated on the surface modelling the single fracture. The two problems are then coupled through physically consistent conditions to account for the exchange of fluid between them. We remark that this model is able to handle both fractures with low and large permeability. Moreover, its extension to the case of two-phase flows has been addressed in ^[49,50], while the case of a totally immersed fracture has been considered e.g., in ^[51].

Even if the use of this kind of dimensionally reduced models avoids the need for extremely refined grids inside the fracture domains, in realistic cases, the construction of a computational grid aligned with the fractures is still a major issue. For example, a fractured oil reservoir can be cut by several thousands of fractures, which often intersect, create small angles or are nearly coincident ^[52]. In line with the previous discussion, in order to avoid the limitations imposed by standard finite element methods, various numerical methods supporting polytopic elements have been employed in the literature for the approximation of the coupled bulk-fracture problem. For example, in ^[8,52] a mixed approximation based on the use of mimetic finite difference method has been explored; in ^[53,54] a framework for treating flows in discrete fracture networks based on the virtual element method has been introduced, and in ^[55] the hybrid high-order method has been employed. We also mention that an alternative strategy consists in the use of non-conforming discretizations, where fractures are allowed to arbitrarily cut the bulk grid, which can then be chosen fairly regular. In particular, we refer to ^[49,56,57], where an approximation employing extended finite element method (XFEM) has been proposed and to the recent work ^[58], where the use of the cut finite element method (CUTFEM) has been explored.

Recently in ^[59], in the setting of conforming discretizations, we developed a numerical approximation of the coupled bulk-fracture model based on polytopic discountinuous Galerkin (PolyDG) methods. In particular, the intrinsic "discontinuous" nature of DG methods allows very general polytopic elements because of the freedom in representing the underlying (local) polynomial space. Indeed the degrees of freedom are not "attached" to any geometric quantity, (vertexes, edges, etcc), so that mesh elements with edges/faces that may be in arbitrary number and whose measure may be arbitrarily small compared to the diameter of the corresponding element are naturally supported with a solid theoretical background. This approach is then very well suited to tame the geometrical complexity featured by most of applications in the computational geoscience field. Moreover, since the interface conditions between the bulk and fracture problem can be naturally formulated using jump and average operators, the coupling of the two problems can be naturally embedded in the variational formulation.

The goal of this paper is to extend the results obtained in ^[59], where the proposed discretization based on PolyDG was in a primal-primal setting. Indeed, when dealing with the approximation of Darcy's flow there are two possible approaches: primal and mixed. The primal approach considers a single-field formulation with the pressure field of the fluid as only unknown. The mixed approach describes the flow not only through the pressure field, but also through an additional unknown representing Darcy's velocity, which is of primary interest in many engineering applications. The primal setting has the advantage of featuring less degrees of freedom and leads to a symmetric positive definite algebraic system of equations that can be efficiently solved based on employing, for example, multigrid techniques ^[39,44,60]. In this case, Darcy's velocity is afterwards reconstructed taking the gradient of the computed pressure and multiplying it by the permeability tensor. However, this process usually entails a loss of accuracy and does not guarantee mass conservation, see for example ^[61,62]. For this reason, the mixed setting is sometimes preferred. In this case Darcy's velocity is directly computed, so that a higher degree of accuracy is achieved, together with local and global mass conservation. However, the drawback of this approach is the complexity of the resulting scheme, which leads to a generalized saddle point algebraic system of equations. From the above discussion we may infer that, according to the desired approximation properties of the model, one may resort to either a primal or mixed approximation for the problem in the bulk, as well as to a primal or mixed approximation for the problem in the fracture. Our aim is then to design and analyze, in the unified framework of ^[63] based on the flux-formulation, all the possible combinations of primal-primal, mixed-primal, primal-mixed and mixed-mixed formulations for the bulk and fracture problems, respectively. In particular, the primal discretizations are obtained using the symmetric interior penalty discontinuous Galerkin method ^[64,65], whereas the mixed discretizations are based on employing the local DG (LDG) method of ^[66], both in their generalization to polytopic grids ^{[36,37,38,40,41]}. Moreover, the coupling conditions between bulk and fracture are imposed through a suitable definition of the numerical fluxes on the fracture faces. Such an abstract setting allows to analyse theoretically at the same time all the possible formulations. We perform a unified analysis of all the derived combinations of DG discretizations for the bulk-fracture problem. We prove their well-posedness and derive a priori $hp$-version error estimates in a suitable (mesh-dependent) energy norm. Finally, we present preliminary numerical experiments assessing the validity of the theoretical error estimates and comparing the accuracy of the proposed formulations on simplified geometries with manufactured solutions.

The rest of the paper is organized as follows. In Section 2 we introduce the model problem; its weak formulation is discussed in Section 3. The discretization based on employing PolyDG methods is presented, in the unified setting of ^[63], in Section 4. In Section 5, we address the problem of stability and prove that all formulations, namely primal-primal (PP), mixed-primal (MP), primal-mixed (PM) and mixed-mixed (MM) are well-posed. The corresponding error analysis is presented in Section 6. Several numerical tests, focusing, for the sake of brevity, on the primal-primal (PP) and mixed-primal (MP) cases, are presented in Section 7 to confirm the theoretical bounds. Moreover, we assess the capability of the method of handling more complicated geometries, presenting some test cases featuring networks of partially immersed fractures.

2.   Model problem

For simplicity, we consider the case where the porous medium is cut by a single, non immersed fracture. The extension to the case of a network of disjoint fractures can be treated analogously. The case where the fracture is partially or totally immersed in the domain is more complex to analyze, and we refer to ^[51,52] for its discussion. Nevertheless, the capability of our method to deal with networks of partially immersed fractures will be explored via numerical experiments in Section 7.4 in the mixed-primal setting (MP). Finally, the case of a network of interecting fractures will be the object of a future work and we refer to ^[59] for preliminary numerical results (in the primal-primal setting) showing that our method is able to handle also such cases.

The porous matrix is represented by the domain $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, which we assume to be open, bounded, convex and polygonal/polyhedral. Moreover, following the strategy of ^[47], we suppose that the fracture may be described by the $(d-1)$-dimensional $\mathcal{C}^{\infty}$ manifold (with no curvature) $\Gamma \subset \mathbb{R}^{d-1}$, $d = 2, 3$. This approach is justified by the fact that the thickness of the fracture domain is typically some orders of magnitude smaller than the size of the domain. Since we are assuming that $\Gamma$ is not immersed, it separates $\Omega$ into two connected disjoint subdomains, $\Omega \setminus \Gamma = \Omega_1 \cup \Omega_2$ with $\Omega_1 \cap \Omega_2 = \emptyset$. We decompose the boundary of $\Omega$ into two disjoint subsets $\partial \Omega_D$ and $\partial \Omega_N$, i.e., $\partial \Omega = \partial \Omega_D \cup \partial \Omega_N$, with $\partial \Omega_D \cap \partial \Omega_N = \emptyset$, and we define $\partial \Omega_{D, i} = \partial \Omega_D \cap \partial \Omega_i$ and $\partial \Omega_{N, i} = \partial \Omega_N \cap \partial \Omega_i$, for $i = 1, 2$. Finally, we denote by $\textbf{n}_i$, $i = 1, 2$ the unit normal vector to $\Gamma$ pointing outwards from $\Omega_i$ and, for a (regular enough) scalar-valued function $q$ and a (regular enough) vector-valued function $\textbf{v}$, we define the classical jump and average operators across the fracture $\Gamma$ as

where the subscript $i = 1, 2$ denotes the restriction to the subdomain $\Omega_i$. Note that, since we are assuming that the fracture is continuously differentiable, it holds $\textbf{n}_1 = - \textbf{n}_2$. We also denote by $\textbf{n}_{\Gamma}$ the normal unit vector on $\Gamma$ with a fixed orientation from $\Omega_1$ to $\Omega_2$, so that we have $\textbf{n}_{\Gamma} = \textbf{n}_1 = -\textbf{n}_2$. In Figure 1 we report an example of domain cut by a single fracture.

We can now introduce the governing equations for our model. In the bulk, we suppose that the flow is governed by Darcy's law. The motion of an incompressible fluid in each domain $\Omega_i$, $i = 1, 2$, with pressure $p_i$ and velocity $\textbf{u}_i$ may then be described by:

where $f \in L^2(\Omega)$ represents a source term, $g_D \in H^{1/2}(\partial \Omega_D)$ is the Dirichlet boundary datum and $\boldsymbol{\nu} = \boldsymbol{\nu}(x) \in \mathbb{R}^{d \times d}$ is the bulk permeability tensor, which we assume to be symmetric, positive definite, uniformly bounded from below and above and with entries that are bounded, piecewise continuous real-valued functions.

On the manifold $\Gamma$ representing the fracture, we formulate a reduced version of Darcy's law in the tangential direction (we refer to ^[47] for a rigorous derivation of the model). To this aim we assume that the fracture permeability tensor $\boldsymbol{\nu}_{\Gamma}$, has a block-diagonal structure of the form

when written in its normal and tangential components. Here, $\boldsymbol{\nu}_{\Gamma}^{\tau} \in \mathbb{R}^{(d-1) \times (d-1)}$ is a positive definite, uniformly bounded tensor (it reduces to a positive number for $d = 2$). Moreover, we assume that $\boldsymbol{\nu}_{\Gamma}$ satisfies the same regularity assumptions as those satisfied by the bulk permeability $\boldsymbol{\nu}$. Setting $\partial\Gamma = \Gamma \cap \partial \Omega$, $\partial \Gamma = \partial \Gamma_D \cup \partial \Gamma_N$, introducing the fracture thickness $\ell_{\Gamma} > 0$ and denoting by $p_{\Gamma}$ and $\textbf{u}_{\Gamma}$ the fracture pressure and velocity, the governing equations for the fracture flow are

where $f_{\Gamma} \in L^2(\Gamma)$, $g_{\Gamma} \in H^{1/2}(\partial \Gamma)$, $\boldsymbol{\tau}$ is vector in the tangent plane of $\Gamma$ normal to $\partial \Gamma$ and $\nabla_{\tau}$ and $\nabla_{\tau} \cdot$ denote the tangential gradient and divergence operators, respectively.

Finally, following ^[47], we close the model providing the interface conditions to couple problems (2.2) and (2.4) along their interface. Given a positive real number $\xi \ne \frac{1}{2}$ that will be chosen later on, the coupling conditions read as follows

where $\beta_{\Gamma} = \frac{1}{2 \eta_{\Gamma}}$ and $\alpha_{\Gamma} = \frac{2}{\eta_{\Gamma}(2 \xi-1)}$ and $\eta_{\Gamma} = \frac{\ell_{\Gamma}}{\nu_{\Gamma}^n}$, $\nu_{\Gamma}^n$ being the normal component of the fracture permeability tensor, see (2.3).

In conclusion, the coupled bulk-fracture model problem is the following:

3.   Weak formulation

In this section we introduce the weak formulation of our model problem (2.6) and prove its well-posedness. We start with the introduction of the functional setting.

3.1.   Functional setting

We will employ the following notation. For an open, bounded domain $D \subset \mathbb{R}^d$, $d = 2, 3$, we denote by $H^s(D)$ the standard Sobolev space of order $s$, for a real number $s \geq 0$. For $s = 0$, we write $L^2(D)$ in place of $H^0(D)$. The usual norm on $H^s(D)$ is denoted by $|| \cdot||_{s, D}$ and the usual seminorm by $| \cdot|_{s, D}$. We also introduce the standard space $H_{div}(D) = \{\textbf{v} : D \rightarrow \mathbb{R}^d : \, ||\textbf{v}||_{0, D} + || \nabla \cdot \textbf{v}||_{0, D} < \infty\}$. Given a decomposition of the domain into elements $\mathcal{T}_h$, we will denote by $H^s(\mathcal{T}_h)$ the standard broken Sobolev space, equipped with the broken norm $|| \cdot||_{s, \mathcal{T}_h}$. Furthermore, we will denote by $\mathbb{P}_k(D)$ the space of polynomials of total degree less than or equal to $k \geq 1$ on $D$. The symbol $\lesssim$ (and $\gtrsim$) will signify that the inequalities hold up to multiplicative constants which are independent of the discretization parameters, but might depend on the physical parameters.

Next, we introduce the functional spaces for our weak formulation. For the bulk pressure and velocity, we introduce the spaces $M^b = L^2(\Omega)$ and $\textbf{V}^b = \{ \textbf{v} \in H_{div}(\Omega): \, \left[\kern-0.15em\left[ { \textbf{v} } \right]\kern-0.15em\right]|_\Gamma \in L^2(\Gamma), \{\textbf{v}\}|_\Gamma \in [L^2(\Gamma)]^d, \textbf{v} \cdot \textbf{n}|_{\partial \Omega_N} = 0 \}$, and equip the space $\textbf{V}^b$ with the norm $||\textbf{v}||_{\textbf{V}^b}^2 = ||\textbf{v}||_{0, \Omega}^2 + || \nabla \cdot \textbf{v}||_{0, \Omega}^2 + ||\left[\kern-0.15em\left[ { \textbf{v} } \right]\kern-0.15em\right]||_{0, \Gamma}^2 + ||\{\textbf{v}\}||_{0, \Gamma}^2.$

Similarly, for the fracture pressure and velocity we define the spaces $M^\Gamma = L^2(\Gamma)$ and $\textbf{V}^\Gamma = \{\textbf{v}_\Gamma \in H_{div, \tau}(\Gamma): \, \textbf{v}_\Gamma \cdot \boldsymbol{\tau}|_{\partial \Gamma } = 0\}$. The norm on $\textbf{V}^\Gamma$ is given by $||\textbf{v}_\Gamma||_{\textbf{V}^\Gamma}^2 = ||\textbf{v}_\Gamma ||_{0, \Gamma}^2+ ||\nabla_\tau \cdot \textbf{v}_\Gamma||_{0, \Gamma}^2.$ Finally, we define the global spaces for the pressure and the velocity as $M = M^b \times M^\Gamma$ and $\textbf{W} = \textbf{V}^b \times \textbf{V}^\Gamma$, equipped with the canonical norms for product spaces. In order to deal with non-homogeneous boundary conditions, we also introduce the affine spaces $V^b_g = \mathcal{L}_g+ \textbf{V}^b$ and $V^\Gamma_g = \mathcal{L}_{g_\Gamma} + \textbf{V}^\Gamma$, where $\mathcal{L}_g \in H_{div}(\Omega)$ and $\mathcal{L}_{g_\Gamma} \in H_{div, \tau}(\Gamma)$ are liftings of the boundary data $g$ and $g_\Gamma$, respectively. We can then define the global space $\textbf{W}_g = \textbf{V}^b_g \times \textbf{V}^\Gamma_g$.

3.2.   Weak problem

We can now formulate problem (2.6) in weak form as follows: Find $(\textbf{u}, \textbf{u}_\Gamma) \in \textbf{W}_g$ and $(p, p_\Gamma) \in M$ such that

where the bilinear form $A(\cdot, \cdot): \textbf{W}_g \times \textbf{W}_g \rightarrow \mathbb{R}$ is defined as $A((\textbf{u}, \textbf{u}_\Gamma), (\textbf{v}, \textbf{v}_\Gamma)) = a(\textbf{u}, \textbf{v})+a_\Gamma(\textbf{u}_\Gamma, \textbf{v}_\Gamma)$ with

and the bilinear form $B(\cdot, \cdot): \textbf{W}_g \times M \rightarrow \mathbb{R}$ is defined as $B((\textbf{v}, \textbf{v}_\Gamma), (q, q_\Gamma)) = b(\textbf{v}, q) + b_\Gamma(\textbf{v}_\Gamma, q_\Gamma) + d(\textbf{v}, q_\Gamma)$, with

Finally the linear operators $F^\textbf{u}(\cdot): \textbf{W}_g \rightarrow \mathbb{R}$ and $F^p(\cdot): M \rightarrow \mathbb{R}$ are defined as

Next, we prove that formulation (3.1) is well-posed. For the sake of simplicity, we will assume that homogeneous Dirichlet boundary conditions are imposed for both the bulk and fracture problems, i.e., $g_{D, i} = 0$, $i = 1; 2$, and $g_\Gamma = 0$ and that the domain and fracture are smooth enough. The extension to the general non-homogeneous case is straightforward. Note that the existence and uniqueness of the problem can be proven only under the condition that the parameter $\xi > 1/2$.

Theorem 3.1. Suppose that $\xi > 1/2$. Then problem (3.1) admits a unique solution.

Proof. For the proof we follow the technique of ^[47]. First, we define the subspace of $\textbf{W}$, $\widehat{ \textbf{W}} = \{(\textbf{v}, \textbf{v}_\Gamma) \in \textbf{W}: \, B((\textbf{v}, \textbf{v}_\Gamma), (q, q_\Gamma)) = 0 \; \forall (q, q_\Gamma) \in M\}$. To show existence and uniqueness of the solution of (3.1), we only need to show that $A(\cdot, \cdot)$ is $\widehat{\textbf{W}}$-elliptic and that $B(\cdot, \cdot)$ satisfies the inf-sup condition, that is

First, we prove that $A(\cdot, \cdot)$ is $\widehat{\textbf{W}}$-elliptic. Since for elements in $(\textbf{v}, \textbf{v}_\Gamma) \in \widehat{\textbf{W}}$ we have $\nabla \cdot \textbf{v} = 0$ in $L^2(\Omega)$ and $\nabla_\tau \cdot \textbf{v}_\Gamma = \left[\kern-0.15em\left[ { \textbf{v} } \right]\kern-0.15em\right]|_\Gamma$ in $L^2(\Gamma)$, the norm $||(\textbf{v}, \textbf{v}_\Gamma)||_{\textbf{W}}$ is equivalent to $||\textbf{v}||_{0, \Omega}^2 + ||\textbf{v}_\Gamma||_{0, \Gamma}^2 +||\left[\kern-0.15em\left[ { v } \right]\kern-0.15em\right] ||_{0, \Gamma}^2 + || \{\textbf{v}\}||_{0, \Gamma}^2$. Owing to the regularity properties of the permeability tensors $\boldsymbol{\nu}$ and $\boldsymbol{\nu}_\Gamma$, this implies that

Note that the hidden constant also depends on the parameter $\alpha_\Gamma$, and that we need to assume $\alpha_\Gamma > 0$, or, equivalently, $\xi > 1/2$, for the inequality to hold true.

To show that $B$ satisfies the inf-sup condition, given $(q, q_\Gamma) \in M$, we construct, exploiting the adjoint problem, $(\textbf{v}, \textbf{v}_\Gamma)\in \textbf{W}$ such that $B((\textbf{v}, \textbf{v}_\Gamma), (q, q_\Gamma)) = ||(q, q_\Gamma)||^2_M$ and $||(\textbf{v}, \textbf{v}_\Gamma)||_{\textbf{W}} \lesssim ||(q, q_\Gamma)||_M$. Given $(q, q_\Gamma) \in M$, let $(\phi, \phi_\Gamma)$ be the solution of

If we set $\textbf{v} = (\textbf{v}_1, \textbf{v}_2)$ with $\textbf{v}_i = - \nabla \phi|_{\Omega_i}$, $i = 1, 2$, and $\textbf{v}_\Gamma = -\nabla_\tau \phi_\Gamma$, we obtain $\nabla \cdot \textbf{v} = q$, $\nabla_\tau \cdot \textbf{v}_\Gamma = q_\Gamma$ and $\left[\kern-0.15em\left[ { \textbf{v} } \right]\kern-0.15em\right]|_\Gamma = 0$, since $\textbf{v} \in H^1(\Omega)$. This implies that $(\textbf{v}, \textbf{v}_\Gamma) \in \textbf{W}$ and $B((\textbf{v}, \textbf{v}_\Gamma), (q, q_\Gamma)) = ||q||_{0, \Omega}^2 + ||q_\Gamma||_{0, \Gamma}^2 = ||(q, q_\Gamma)||^2_M$. Finally, from elliptic regularity, we have $||(\textbf{v}, \textbf{v}_\Gamma)||_{\textbf{W}}^2 = ||\nabla \phi||^2_{0, \Omega} + ||\nabla_\tau \phi_\Gamma||^2_{0, \Gamma} + ||q||^2_{0, \Omega} + ||q_\Gamma||^2_{0, \Gamma} + ||\{ \nabla \phi \}||_{0, \Gamma}^2 \lesssim ||q||^2_{0, \Omega} + ||q_\Gamma||^2_{0, \Gamma}$, and this concludes the proof.

4.   Numerical dicretization based on PolyDG methods

In this section we present a family of discrete formulations for the coupled bulk-fracture problem (3.1), which are based on discontinuous Galerkin methods on polytopic grids. In particular, since we can choose to discretize the problem in the bulk and the one in the fracture either in their mixed or in their primal form, we derive four formulations that embrace all the possible combinations of primal-primal, mixed-primal, primal-mixed and mixed-mixed discretizations. The mixed discretizations will be based on the local discontinuous Galerkin method (LDG) ^[66,67,68], while the primal discretizations on the Symmetric Interior Penalty method (SIPDG) ^[64,65], all supporting polytopic grids ^{[36,37,40,41]}. The derivation of our discrete formulations will be carried out following the same strategy as in ^[63], so that it will be based on the introduction of the numerical fluxes, which approximate the trace of the solutions on the boundary of each mesh element. In particular, the imposition of the coupling conditions (2.5a)-(2.5b) will be achieved through a proper definition of the numerical fluxes on the faces belonging to the fracture.

First, we introduce the notation related to the discretization of the domains by means of polytopic meshes. For the problem in the bulk, we consider a family of meshes $\mathcal{T}_h$ made of disjoint open polygonal/polyhedral elements. Following ^[36,37,40], we introduce the concept of mesh interface, defined as the intersection of the $(d-1)$-dimensional facets of two neighbouring elements. We need now to distinguish between the case when $d = 3$ and $d = 2$:

● when $d = 3$, each interface consists of a general polygon, which we assume may be decomposed into a set of co-planar triangles. We assume that a sub-triangulation of each interface is provided and we denote the set of all these triangles by $\mathcal{F}_h$. We then use the terminology face to refer to one of the triangular elements in $\mathcal{F}_h$;

● when $d = 2$, each interface simply consists of a line segment, so that the concepts of face and interface are in this case coincident. We still denote by $\mathcal{F}_h$ the set of all faces.

Note that $\mathcal{F}_h$ is always defined as a set of $(d-1)$-dimensional simplices (triangles or line segments). As in ^[36,37,40], no limitation on either the number of faces of each polygon $E \in \mathcal{T}_h$ or on the relative size of the faces compared to the diameter of the element is imposed.

We consider meshes $\mathcal{T}_h$ that are aligned with the fracture $\Gamma$, so that any element $E \in \mathcal{T}_h$ cannot be cut by $\Gamma$ and it belongs exactly to one of the two disjoint subdomains $\Omega_1$ or $\Omega_2$. This implies that each mesh $\mathcal{T}_h$ induces a subdivision of the fracture $\Gamma$ into faces, which we denote by $\Gamma_h$. It follows that we can write

where $\mathcal{F}_h^B$ is the set of faces lying on the boundary of the domain $\partial \Omega$ and $\mathcal{F}_h^I$ is the set of interior faces not belonging to the fracture. In addition, we write $\mathcal{F}_h^B = \mathcal{F}_h^D \cup \mathcal{F}_h^N$, where $\mathcal{F}_h^D$ and $\mathcal{F}_h^N$ are the boundary faces contained in $\partial \Omega_D$ and $\partial \Omega_N$, respectively (we assume the decomposition to be matching with the partition of $\partial \Omega$ into $\partial \Omega_D$ and $\partial \Omega_N$).

The induced discretization of the fracture $\Gamma_h$ consists of the faces of the elements of $\mathcal{T}_h$ that share part of their boundary with the fracture, so that $\Gamma_h$ is made up of line segments when $d = 2$ and of triangles when $d = 3$. Note that, in the $3D$ case, the triangles are not necessarily shape-regular and they may present hanging nodes, due to the fact that the sub-triangulations of each elemental interface is chosen independently from the others. For this reason, we extend the concept of interface also to the $(d-2)$-dimensional facets of elements in $\Gamma_h$, defined again as intersection of boundaries of two neighbouring elements. When $d = 2$, the interfaces reduce to points, while when $d = 3$ they consists of line segments. We denote by $\mathcal{E}_{\Gamma, h}$ the set of all the interfaces (that we will also call edges) of the elements in $\Gamma_h$, and we write, accordingly to the previous notation, $\mathcal{E}_{\Gamma, h} = \mathcal{E}_{\Gamma, h}^I \cup \mathcal{E}_{\Gamma, h}^B$, with $\mathcal{E}_{\Gamma, h}^B = \mathcal{E}_{\Gamma, h}^D \cup \mathcal{E}_{\Gamma, h}^N$.

For each element $E \in \mathcal{T}_h$, we denote by $|E|$ its measure, by $h_E$ its diameter and we set $h = \max_{E \in \mathcal{T}_h} h_E$. Given an element $E \in \mathcal{T}_h$, for any face $F \subset \partial E$, with $F \in \mathcal{F}_h$, we define $\textbf{n}_F$ as the unit normal vector on $F$ that points outward of $E$. We can then define the standard jump and average operators across a face $F \in \mathcal{F}_h \setminus \mathcal{F}_h^B$ for (regular enough) scalar and vector-valued functions similarly to (2.1). We also recall a well-known identity ^[65] for scalar and vector-valued functions $q$ and $\textbf{v}$ that are piecewise smooth on $\mathcal{T}_h$:

where we have used the compact notation $\int_{\mathcal{F}_h} \cdot = \sum_{F \in \mathcal{F}_h} \int_F \cdot$ and jump and average operators on a boundary face $F \in \mathcal{F}_h^B$ are defined as $\left[\kern-0.15em\left[ { q } \right]\kern-0.15em\right] = q \textbf{n}_F$ and $\{\textbf{v}\} = \textbf{v}$.

Analogous definitions may be also set up on the fracture. In particular, given an element $F \in \Gamma_h$, with measure $|F|$ and diameter $h_F$, for any edge $e \subset \partial F$, with $e \in \mathcal{E}_{\Gamma, h}$, we define $\textbf{n}_e$ as the unit normal vector on $e$ pointing outward of F (it reduces to $\pm 1$ when $d = 2$). Finally, standard jump and average operators across every edge $e$ can be defined for (regular enough) scalar and vector-valued functions and an analogous version of formula (4.1) can be stated for piecewise smooth function on the fracture mesh $\Gamma_h$.

We have now all the ingredients to introduce the discrete formulation of model problem (3.1).

4.1.   Discrete formulation

For simplicity in the forthcoming analysis, we will suppose that the permeability tensors $\boldsymbol{\nu}$ and $\boldsymbol{\nu}_{\Gamma}$ are piecewise constant on mesh elements, i.e., $\boldsymbol{\nu}|_E \in [\mathbb{P}_0(E)]^{d \times d}$ for all $E \in \mathcal{T}_h$, and $\boldsymbol{\nu}_{\Gamma}|_F \in [\mathbb{P}_0(F)]^{(d-1) \times (d-1)}$ for all $F \in \Gamma_h$. First, we introduce the finite-dimensional spaces where we will set our discrete problem. We set

Note that, to each element $E \in \mathcal{T}_h$ is associated the polynomial degree $k_E \geq 1$, as well as to each face $F \in \Gamma_h$ is associated the degree $k_F \geq1$. We remark that the polynomial degrees in the bulk and fracture discrete spaces just defined are chosen independently of each other.

We first focus on the problem in the bulk. Multiplying the first and second equations in (2.2) by test functions $\textbf{v} \in \textbf{W}_h^b$ and $q \in Q_h^b$, respectively, and integrating by parts over an element $E \in \mathcal{T}_h$, we obtain

In the spirit of ^[63], we start the derivation of our DG discretization from these equations. Adding over the elements $E\in \mathcal{T}_h$, the general discrete formulation for the problem in the bulk will then be: Find $p_h\in Q_h^b$ and $\textbf{u}_h \in \textbf{W}_h^b$, such that for all $E \in \mathcal{T}_h$ we have

where the numerical fluxes $\hat{p}_E$ and $\hat{\textbf{u}}_E$ are approximations to the exact solutions $\textbf{u}$ and $p$, respectively, on the boundary of $E$. The definition of the numerical fluxes in terms of $p_h$, $\textbf{u}_h$, of the boundary data and of the coupling conditions (2.5a)-(2.5b) will determine the method. Using identity (4.1), we get

where we have introduced $\hat{p} = (\hat{p}_E)_{E \in \mathcal{T}_h}$ and $\hat{\textbf{u}} = (\hat{\textbf{u}}_E)_{E \in \mathcal{T}_h}$. The numerical fluxes $\hat{p}$ and $\hat{\textbf{u}}$ must be interpreted as linear functionals taking values in the spaces $\Pi_{E \in \mathcal{T}_h} L^2(\partial E)$ and $[\Pi_{E \in \mathcal{T}_h} L^2(\partial E)]^d$, respectively. In particular, this means that they are, in general, double-valued on $\mathcal{F}_h^I \cup \Gamma_h$ and single-valued on $\mathcal{F}_h^B$. We also observe for future use that, after integrating by parts and using again identity (4.1), Eq. (4.2) may also be rewritten as

We now reason analogously for the fracture. Multiplying the first and second equations in (2.4) by test functions $\textbf{v}_{\Gamma}$ and $q_{\Gamma}$, respectively, integrating by parts over an element $F \in \Gamma_h$ and summing over all the elements in $\Gamma_h$ we obtain the following problem: Find $p_{\Gamma, h} \in Q_h^{\Gamma}$ and $\textbf{u}_{\Gamma, h} \in \textbf{W}_h^{\Gamma}$ such that for all $F \in \Gamma_h$ we have

Here, we have introduced the numerical fluxes $\hat{p}_{\Gamma, F}$ and $\hat{\textbf{u}}_{\Gamma, F}$. Again, the idea is that they represent approximations on the boundary of the fracture face $F$ of the exact solutions $p_{\Gamma}$ and $\textbf{u}_{\Gamma}$, respectively. Note also that here $\hat{\textbf{u}}$ is the numerical flux approximating the bulk velocity on $\Gamma_h$. Using identity (4.1), we get

We point out that, in all previous equations, the gradient and divergence operators are actually tangent operators. Here, we have dropped the subscript $\tau$ in order to simplify the notation.

In the following, we explore all possible combinations of primal-primal, mixed-primal, primal mixed and mixed-mixed formulations for the bulk and fracture, respectively.

4.2.   Primal-Primal formulation

In order to obtain the primal-primal formulation, we need to eliminate the velocities $\textbf{u}_h$ and $\textbf{u}_{\Gamma, h}$ from equations (4.2)-(4.3) and (4.5)-(4.6). To do so, we need to express $\textbf{u}_h$ solely in terms of $p_h$ (and $p_{\Gamma, h}$), and $\textbf{u}_{\Gamma, h}$ solely in terms of $p_{\Gamma, h}$. As in ^[63] this will be achieved via the definition of proper lifting operators.

We start by focusing on the problem in the bulk. In order to complete the specification of the DG method that we want to use for the approximation, we need to give an expression to the numerical fluxes. We choose the classic symmetric interior penalty method (SIPDG). Moreover, coupling conditions (2.5a)-(2.5b) are imposed through a suitable definition of the numerical fluxes on the fracture faces. Since we want a primal formulation, the definition of $\hat{p}$ and $\hat{\textbf{u}}$ will not contain $\textbf{u}_h$. The numerical fluxes are defined as follows:

Here, we have introduced the discontinuity penalization parameter $\sigma$. In particular, $\sigma$ is a non-negative bounded function, i.e., $\sigma \in L^{\infty}(\mathcal{F}_h^I \cup \mathcal{F}_h^D)$ and its precise definition will be given in Definition 5.2 below. Moreover, we have used the notation $\sigma_F = \sigma|_{F}$, for $F \in \mathcal{F}_h^I \cup \mathcal{F}_h^D$. We remark that, with this choice, the numerical flux $\hat{p}$ is doubled valued on $\Gamma_h$ and single valued on $\mathcal{F}_h^I \cup \mathcal{F}_h^B$.

Using the definition of the numerical fluxes, it follows that

so we rewrite (4.4) as

At this point, we proceed with the elimination of the auxiliary variable $\textbf{u}_h$ from our equations. To this aim, we introduce the lifting operator $\mathscr{L}_b^{ SIP}:[L^1(\mathcal{F}_h^I \cup \mathcal{F}_h^D)]^d \rightarrow \textbf{W}_h^b$ defined by

Similarly, we define the lifting $\mathcal{G}_b(g_D) \in \textbf{W}_h^b$ of the Dirichlet boundary datum $g_D$ as

Thanks to the introduction of the lifting operators, Eq. (4.8) may be rewritten as

Since $\nabla Q_h^b \subseteq \textbf{W}_h^b$, we can write

where $\nabla p_h + \mathscr{L}_b^{ SIP}(\left[\kern-0.15em\left[ { p_h } \right]\kern-0.15em\right]) + \mathcal{G}_b$ can be seen as a discrete approximation of the gradient $\nabla p$.

We can then rewrite Eq. (4.3) as

Using the definition of the discrete gradient (4.11), of the lifting operators (4.9) and (4.10) and of the numerical flux $\hat{\textbf{u}}$ (4.7), we have

Now we move our attention to the problem in the fracture. We define the numerical fluxes $\hat{p}_{\Gamma}$ and $\hat{\textbf{u}}_{\Gamma}$ in order to obtain a symmetric interior penalty approximation as follows:

Again, we have introduced the discontinuity penalization parameter $\sigma_\Gamma \in L^\infty(\mathcal{E}_{\Gamma, h}^I \cup \mathcal{E}_{\Gamma, h}^D)$ and we set $\sigma_e = \sigma_\Gamma|_e$ for $e \in \mathcal{E}_{\Gamma, h}^I \cup \mathcal{E}_{\Gamma, h}^D$. Its precise definition will be given in Definition 5.3 below. Next, as before, we introduce the lifting operator $\mathscr{L}_{\Gamma}^{ SIP}:[L^1(\mathcal{E}_{\Gamma, h}^I \cup \mathcal{E}_{\Gamma, h}^D)]^{d-1} \rightarrow \textbf{W}_h^{\Gamma}$ and the lifting of the boundary datum $\mathcal{G}_{\Gamma}(g_{\Gamma, D}) \in \textbf{W}_h^{\Gamma}$ defined by

Integrating by parts and using (4.1), we can rewrite Eq. (4.5) as

Again, since $\nabla Q_h^{\Gamma} \subseteq \textbf{W}_h^{\Gamma}$ elementwise, we can write

Plugging this last identity and the definition of the numerical fluxes $\hat{\textbf{u}}$ (see (4.7)) and $\hat{\textbf{u}}_{\Gamma}$ (see (4.13)) into Eq. (4.6), we obtain

In conclusion, summing Eqs. (4.12) and (4.16) we obtain the following discrete formulation: Find $(p_h, p_h^{\Gamma}) \in Q_h^b \times Q_h^{\Gamma}$ such that

where $PP$ stands for primal-primal and where $\mathcal{L}_h: Q_h^b \times Q_h^{\Gamma} \rightarrow \mathbb{R}$ is defined as $\mathcal{L}_h^{ PP}(q, q_{\Gamma}) = \mathcal{L}_b^{ P}(q)+ \mathcal{L}_{\Gamma}^{ P}(q_{\Gamma})$ and $\mathcal{A}_h^{ PP}: (Q_h^b \times Q_h^{\Gamma}) \times (Q_h^b \times Q_h^{\Gamma}) \rightarrow \mathbb{R}$ is defined as

with

and

We remark that we have recovered the formulation already obtained in ^[59] (in its not strongly consistent version), the only difference being that the bilinear form for the problem in the fracture is in the shape of SIPDG method, instead of classical conforming finite elements.

4.3.   Mixed-Primal formulation

In this section, we discretize the problem in the bulk in its mixed form. To this aim, we use the local discontinuous Galerkin (LDG) method ^{[66,67,68,69]}. The LDG method is a particular DG method that can be included in the class of mixed finite element methods. However, the variable $\textbf{u}_h$ can be locally solved in terms of $p_h$ and then eliminated from the equations, giving rise to a primal formulation with $p_h$ as only unknown.

In what follows, we first derive the formulation of our method in a mixed setting. After that, we recast it in a primal setting, in order to perform the analysis in the framework of ^[63,69]. However, we remark that the mixed formulation is the one that will actually be implemented for the numerical experiments of Section 7. As far as the problem in the fracture is concerned, we work again in a primal setting, using the SIPDG method for the discretization.

In the bulk, we define the numerical fluxes as

Here, $\textbf{b} \in [L^\infty(\mathcal{F}_h^I)]^d$ is a (possibly null) vector-valued function which is constant on each face. It is chosen such that

with $B \geq 0$ independent if the discretization parameters. Moreover, $\sigma$ is the penalization parameter introduced in (4.7), whose precise definition will be given in (5.2) below. Note that the numerical flux $\hat{p}$ does not depend on $\textbf{u}_h$. This will allow for an element-by-element elimination of the variable $\textbf{u}_h$, generating a primal formulation of the problem. We also point out that the definition of the numerical fluxes on the fracture faces is the same as in the primal SIPDG setting.

With this definition of the numerical fluxes, and after integration by parts as in (4.4), Eq. (4.2) becomes

while Eq. (4.3) turns into

If we discretize the problem in the fracture with the SIPDG method, we obtain the following discrete mixed problem: Find $\big((p_h, \textbf{u}_h), p_h^{\Gamma} \big) \in Q_h^b \times \textbf{W}_h^b \times Q_h^\Gamma$ such that

where

and $\mathcal{A}_{\Gamma}^{ P}(\cdot, \cdot)$ and $\mathcal{L}_\Gamma^{ P}(\cdot)$ are defined as in (4.19) and (4.22), respectively. Also note that we have $\mathcal{I}((p_h, p_{\Gamma, h}), (q, q_{\Gamma})) = \mathcal{I}_1(p_h, q, p_{\Gamma, h})+ \mathcal{I}_2(p_h, p_{\Gamma, h}, q_{\Gamma})$.

We now focus on rewriting the problem in the bulk in a primal form, taking advantage of the local solvability of LDG method. We proceed as in the SIPDG case and introduce an appropriate lifting operator, $\mathscr{L}_b^{{ LDG}}:[L^1(\mathcal{F}_h^I \cup \mathcal{F}_h^D)]^d \rightarrow \textbf{W}_h^b$, defined by

From Eq. (4.24) we obtain

where $\mathcal{G}_b$ is the lifting of the Dirichlet boundary datum defined in (4.10). Eq. (4.25) now becomes

Using again the definition of the lifting $\mathscr{L}_b^{{ LDG}}$ and the identity (4.29), we obtain

Summing Eqs. (4.30) and (4.16) we obtain the following discrete formulation: Find $(p_h, p_h^{\Gamma}) \in Q_h^b \times Q_h^{\Gamma}$ such that

where $MP$ stands for mixed-primal and where $\mathcal{A}_h^{ MP}: (Q_h^b \times Q_h^{\Gamma}) \times (Q_h^b \times Q_h^{\Gamma}) \rightarrow \mathbb{R}$ is defined as

and $\mathcal{L}^{ MP}_h: Q_h^b \times Q_h^{\Gamma} \rightarrow \mathbb{R}$ is defined as

with

Note that the mixed formulation (4.26) is equivalent to the primal formulation (4.31) together with the definition of the lifting operator (4.28) and Eq. (4.29).

4.4.   Primal-Mixed formulation

We now want to approximate the problem in the fracture in mixed form, employing the LDG method and the problem in the bulk using the SIPDG method. We define the numerical fluxes as follows

Here, $\textbf{b}_\Gamma \in [L^\infty(\mathcal{E}_{\Gamma, h}^I)]^{d-1}$ is a vector-valued function that is constant on each edge and it is chosen such that $||\textbf{b}_\Gamma||_{\infty, \mathcal{E}_{\Gamma, h}^I} \leq B_\Gamma$, with $B_\Gamma \leq 0$ independent of the discretization parameters. Eqs. (4.5) and (4.6) now become

where we have also used the definition of the numerical flux $\hat{\textbf{u}}$ on $\Gamma_h$ (see (4.7)) to rewrite $-\left[\kern-0.15em\left[ { \hat{\textbf{u}} } \right]\kern-0.15em\right] \; = \; \alpha_\Gamma(\{p_h\}-p_{\Gamma, h})$. For the bulk we proceed as in the primal-primal section using for the discretization the SIPDG method. We then obtain the following primal-mixed problem: Find $\big(p_h, (p_h^{\Gamma}, \textbf{u}_{\Gamma, h})\big) \in Q_h^b \times Q_h^\Gamma \times \textbf{W}_h^\Gamma$ such that

where

and $\mathcal{A}_b^{ P}(p_h, q)$ and $\mathcal{L}_b^{ P}(q)$ are defined as in (4.18) and (4.21), respectively.

Aiming at rewriting the problem in the fracture in primal form, we introduce the lifting operator, $\mathscr{L}_\Gamma^{{ LDG}}:[L^1(\mathcal{E}_h^I \cup \mathcal{E}_h^D)]^d \rightarrow \textbf{W}_h^\Gamma$, defined by

From Eq. (4.33) we obtain

where $\mathcal{G}_\Gamma$ is the lifting of the Dirichlet boundary datum defined in (4.15). Eq. (4.34) now becomes

We obtain the following primal formulation: Find $(p_h, p_h^{\Gamma}) \in Q_h^b \times Q_h^{\Gamma}$ such that

where $PM$ stands for primal-mixed and where $\mathcal{A}_h^{ PM}: (Q_h^b \times Q_h^{\Gamma}) \times (Q_h^b \times Q_h^{\Gamma}) \rightarrow \mathbb{R}$ is defined as

and $\mathcal{L}^{ PM}_h: Q_h^b \times Q_h^{\Gamma} \rightarrow \mathbb{R}$ is defined as $\mathcal{L}_h^{ PM}(q, q_{\Gamma}) = \mathcal{L}_b^{ P}(q)+ \mathcal{L}_{\Gamma}^{ M}(q_{\Gamma})$ with

4.5.   Mixed-Mixed formulation

Finally, if we approximate both the problem in the bulk and in the fracture with the LDG method, we obtain the following formulation: Find $(p_h, p_{\Gamma, h}) \in Q_h^b \times Q_h^\Gamma$ and $(\textbf{u}_h, \textbf{u}_{\Gamma, h}) \in \textbf{W}_h^b \times \textbf{W}_h^\Gamma$ such that

This formulation, together with the definition of the lifting operators (4.28) and (4.36) and of the discrete gradients (4.29) and (4.37) is equivalent to the following: Find $(p_h, p_{\Gamma, h}) \in Q_h^b \times Q_h^\Gamma$ such that

where $MM$ stands for mixed-mixed and where $\mathcal{A}_h^{ MM}: (Q_h^b \times Q_h^{\Gamma}) \times (Q_h^b \times Q_h^{\Gamma}) \rightarrow \mathbb{R}$ is defined as

and $\mathcal{L}^{ MM}_h: Q_h^b \times Q_h^{\Gamma} \rightarrow \mathbb{R}$ is defined as $\mathcal{L}_h^{ MM}(q, q_{\Gamma}) = \mathcal{L}_b^{ M}(q)+ \mathcal{L}_{\Gamma}^{ M}(q_{\Gamma}).$

Next, we perform a unified analysis of all of the derived DG discretizations for the fully-coupled bulk-fracture problem. We remark that the analysis will be performed considering the mixed LDG discretizations recast in their primal form, following ^[69]. For clarity, in Table 1 we summarize the bilinear forms for all the four approaches.

The bulk, fracture and interface bilinear forms are defined in:

5.   Well-posedness of the discrete formulations

In this section, we address the problem of stability. We prove that the primal-primal (PP) (4.17), mixed-primal (MP) (4.31), primal-mixed (PM) (4.38) and mixed-mixed (MM)(4.41) formulations are well-posed. We remark that all these formulations are not strongly consistent, due to the presence of the lifting operators. This implies that the analysis will be based on Strang's second Lemma, ^[70].

We recall that, for simplicity in the analysis, we are assuming the permeability tensors $\boldsymbol{\nu}$ and $\boldsymbol{\nu}_\Gamma^\tau$ to be piecewise constant. We will employ the following notation $\bar{\boldsymbol{\nu}}_E = |\sqrt{\boldsymbol{\nu}|_E}|^2_2$ and $\bar{\boldsymbol{\nu}}_F^\tau = |\sqrt{\boldsymbol{\nu}_\Gamma^\tau|_F}|^2_2$, where $|\cdot|_2$ denotes the $l_2$-norm.

To consider the boundedness and stability of our primal bilinear forms, we introduce the spaces $Q^b(h) = Q^b_h + \tilde{Q}^b$ and $Q^\Gamma(h) = Q^\Gamma_h + \tilde{Q}^\Gamma$ where $\tilde{Q}^b = \{ q = (q_1, q_2) \in H^1(\Omega_1) \times H^1(\Omega_2)\} \cap H^2(\mathcal{T}_h)$ and $\tilde{Q}^\Gamma = H^1(\Gamma) \cap H^2(\Gamma_h)$. We remark that all the bilinear forms $\mathcal{A}_h^{**}(\cdot, \cdot)$ are also well-defined on the extended space $Q^b(h) \times Q^\Gamma(h)$.

Further, we introduce the following energy norm on the discrete space $Q_h^b \times Q_h^\Gamma$

where

Note that $|||\cdot|||$ is also well defined on the extended space $Q^b(h) \times Q^\Gamma(h)$.

Since our discretization employ general polytopic grids, we need introduce some technical instruments to work in this framework ^{[36,37,38,40,41]}. In particular, we need trace inverse estimates to bound the norm of a polynomial on a polytope's face/edge by the norm on the element itself. To this aim, we give the following

Definition 5.1. A mesh $\mathcal{T}_h$ is said to be polytopic-regular if, for any $E \in \mathcal{T}_h$, there exists a set of non-overlapping (not necessarily shape-regular) $d$-dimensional simplices $\{S_E^i \}_{i = 1}^{n_E}$ contained in $E$, such that $\bar{F} = \partial \bar{E} \cap \bar{S_E^i}$, for any face $F \subseteq \partial E$, and

with the hidden constant independent of the discretization parameters, the number of faces of the element $n_E$, and the face measure.

We remark that this definition does not give any restriction on the number of faces per element, nor on their measure relative to the diameter of the element the face belongs to.

Assumption 5.1. We assume that $\mathcal{T}_h$ and $\Gamma_h$ are polytopic-regular meshes.

With this hypothesis, we can state the following inverse-trace estimate that is valid for polytopic elements ^[38,41].

Lemma 5.2. Let $E$ be a polygon/polyhedron belonging to a mesh satisfying Definition 5.1 and let $v \in \mathbb{P}_{k_E}(E)$. Then, we have

where the hidden constant depends on the dimension $d$, but it is independent of the discretization parameters, of the number of faces of the element and of the relative size of the face compared to the diameter $k_E$ of $E$.

The second fundamental tool to deal with polytopic discretizations, is an appropriate definition of the discontinuity penalization parameter, which allows for the use of elements with arbitrarily small faces. Taking as a reference ^{[36,37,38,40,41]}, we give the following two definitions for the bulk and fracture penalty functions:

Definition 5.2. The discontinuity-penalization parameter $\sigma: \mathcal{F}_h \setminus \Gamma_h \rightarrow \mathbb{R}^+$ for the bulk problem is defined facewise as

with $\sigma_0 > 0$ independent of $k_E$, $|E|$ and $|F|$.

Definition 5.3. The discontinuity-penalization parameter $\sigma_\Gamma: \mathcal{E}_{\Gamma, h} \rightarrow \mathbb{R}^+$ for the fracture problem is defined edgewise as

with $\sigma_{0, \Gamma} > 0$ independent of $k_F$, $|F|$ and $|e|$.

Now we have all the technical tools to work in a polytopic framework. Next, we will state and prove some estimates that will be instrumental for the proof of the well-posedness of our discrete formulations. We start deriving some bounds on the lifting operators, with arguments similar to those of ^[68,69,71]. Note that all the results hold true on the extended spaces $Q^b(h)$ and $Q^\Gamma(h)$.

Lemma 5.3. Let $\mathscr{L}_b^{ SIP}(\cdot)$ be the lifting operator defined in (4.9). Then, for every $q \in Q^b(h)$ it holds

Proof. Denoting by $\boldsymbol{\Pi}_{\textbf{W}_h^b}$ the $L^2$-projection onto $\textbf{W}_h^b$, by definition of the lifting operator $\mathscr{L}_b^{ SIP}$ and Cauchy-Schwarz inequality, we have

Using the triangular inequality, the definition of the penalization coefficient $\sigma_F$ (5.2), the inverse inequality (5.1), the assumptions on the permeability tensor $\boldsymbol{\nu}$ and the continuity property of the $L^2$-projector we have

This proves the desired estimate.

Lemma 5.4. Let $\mathscr{L}_\Gamma^{ SIP}(\cdot)$ be the lifting operator defined in (4.14). Then, for every $q_\Gamma \in Q^\Gamma(h)$ it holds

Proof. Same arguments as in in the proof of Lemma 5.3.

Lemma 5.5. Let $\mathscr{L}_b^{ LDG}(\cdot)$ be the lifting operator defined in (4.28). Then, for every $q \in Q^b(h)$ it holds

Proof. We proceed as in the proof of Lemma 5.3. By definition of the lifting operator $\mathscr{L}_b^{ LDG}$ and Cauchy-Schwarz inequality, we have

From (5.5) we know that $(a) \lesssim \frac{1}{\sigma_0^{1/2}} ||\sigma_F^{1/2} \left[\kern-0.15em\left[ { q } \right]\kern-0.15em\right] ||_{0, \mathcal{F}_h^I \cup \mathcal{F}_h^D}$, while using similar arguments and bound (4.23) on $\textbf{b}$, we can prove that

so that $(b) \lesssim \frac{B}{\sigma_0^{1/2}} ||\sigma_F^{1/2} \left[\kern-0.15em\left[ { q } \right]\kern-0.15em\right] ||_{0, \mathcal{F}_h^I \cup \mathcal{F}_h^D}$. This concludes the proof.

Lemma 5.6. Let $\mathscr{L}_\Gamma^{ LDG}(\cdot)$ be the lifting operator defined in (4.36). Then, For every $q_\Gamma \in Q^\Gamma(h)$ it holds

Proof. Same arguments as in in the proof of Lemma 5.5.

Using these results, we can now prove that the bilinear forms for the bulk problem $\mathcal{A}_b^P(\cdot, \cdot)$ and $\mathcal{A}_b^M(\cdot, \cdot)$ are continuous on $Q^b(h)$ and coercive on $Q_h^b$, as well as the fracture bilinear forms $\mathcal{A}_\Gamma^P(\cdot, \cdot)$ and $\mathcal{A}_\Gamma^M(\cdot, \cdot)$ are continuous on $Q^\Gamma(h)$ and coercive on $Q_h^\Gamma$.

Lemma 5.7. $\mathcal{A}_b^{ P}(\cdot, \cdot)$ is coercive on $Q_h^b \times Q_h^b$ and continuous on $Q^b(h) \times Q^b(h)$, that is

provided that $\sigma_0$ is chosen big enough.

Proof. We start with coercivity. Taking $p = q \in Q_h^b$, we have

From Young inequality we have

so that, using the bound on the lifting (5.4), we obtain

for $\sigma_0$ big enough.

Continuity directly follows from Cauchy Schwarz inequality and the bound on the lifting (5.4).

Lemma 5.8. $\mathcal{A}_\Gamma^{ P}(\cdot, \cdot)$ is coercive on $Q_h^\Gamma \times Q_h^\Gamma$ and continuous on $Q^\Gamma(h) \times Q^\Gamma(h)$, that is

provided that $\sigma_{0, \Gamma}$ is chosen big enough.

Proof. Analogous to the proof of Lemma 5.7.

Lemma 5.9. $\mathcal{A}_b^{ M}(\cdot, \cdot)$ is coercive on $Q_h^b \times Q_h^b$ and continuous on $Q^b(h) \times Q^b(h)$, that is

Proof. We start with coercivity. From Young's inequality and the bound on the lifting (5.10) we have, for every $0 < \varepsilon < 1$,

with $C = \frac{ (1 + B)}{\sigma_0 \varepsilon}(1-\frac{1}{\varepsilon})$, so that $\mathcal{A}_b^{ M}(\cdot, \cdot)$ is coercive for every choice of the parameters $\sigma_0 > 0$ and $B > 0$^*. Continuity is again a direct consequence of Cauchy Schwarz's inequality and the bound on the lifting (5.6).

^* More in detail: we need $1+C > 0$, with $0 < \varepsilon < 1$. We obtain $1 + (1- \frac{1}{\varepsilon}) \frac{(1+B)^2}{\sigma_0} > 0$, that is $\varepsilon > \frac{1}{1 + \frac{\sigma_0}{(1+B)^2}} = \tilde{C}$, being $0 < \tilde{C} < 1$ for every possible choice of $\sigma_0 > 0$ and $B > 0$.

Lemma 5.10. $\mathcal{A}_\Gamma^{ M}(\cdot, \cdot)$ is coercive on $Q_h^\Gamma \times Q_h^\Gamma$ and continuous on $Q^\Gamma(h) \times Q^\Gamma(h)$, that is

Proof. Analogous to the proof of Lemma 5.9.

Employing Lemmas 5.7, 5.9, 5.8 and 5.10, we can easily prove the well-posedness of all of our discrete problems, as stated in the following stability result.

Proposition 5.11. Let the penalization parameters $\sigma$ for the problem in the bulk and in the fracture be defined as in (5.2) and (5.3), respectively. Then, the fully-coupled discrete problems PP (4.22), MP (4.31), PM (4.38) and MM (4.41) are well-posed provided that $\sigma_0$ and $\sigma_{0, \Gamma}$ are chosen big enough for the primal formulations.

Proof. In order to use Lax-Milgram Lemma, we prove that the bilinear forms $\mathcal{A}_h^{PP}(\cdot, \cdot)$, $\mathcal{A}_h^{MP}(\cdot, \cdot)$, $\mathcal{A}_h^{PM}(\cdot, \cdot)$ and $\mathcal{A}_h^{MM}(\cdot, \cdot)$ are continuous on $Q^b(h) \times Q^\Gamma(h)$ and coercive on $Q_h^b \times Q_h^\Gamma$. We have, from Cauchy-Schwarz inequality

so that coercivity and continuity are a direct consequence of the definition of the norm $||| \cdot|||$ and of Lemmas 5.7, 5.9, 5.8 and 5.10. The continuity of the linear operators $\mathcal{L}^{PP}_h(\cdot)$, $\mathcal{L}^{MP}_h(\cdot)$, $\mathcal{L}^{PM}_h(\cdot)$ and $\mathcal{L}^{MM}_h(\cdot)$ on $Q^b(h) \times Q^{\Gamma}(h)$ can be easily proved by using Cauchy-Schwarz's inequality, thanks to the regularity assumptions on the forcing terms $f$ and $f_{\Gamma}$ and on the boundary data $g_D$ and $g_\Gamma$.

6.   Error analysis

In this section we derive error estimates for our discrete problems.

6.1.   Approximation results

The tool at the base of DG-method error analysis are $hp$-interpolation estimates. Here, we summarize the results contained in ^{[36,37,38,40,41]}, where standard estimates on simplices are extended to arbitrary polytopic elements.

First, we give the following definitions.

Definition 6.1. A covering $\mathcal{T}_{\#} = \{T_E\}$ related to the polytopic mesh $\mathcal{T}_h$ is a set of shape-regular $d$-dimensional simplices $T_E$, such that for each $E \in \mathcal{T}_h$, there exists a $T_E \in \mathcal{T}_{\#}$ such that $E \subsetneq T_E$.

Assumption 6.1. ^{[36,37,38,40,41]} There exists a covering $\mathcal{T}_{\#}$ of $\mathcal{T}_h$ (see Definition 6.1) and a positive constant $O_{\Omega}$, independent of the mesh parameters, such that

and $h_{T_E} \lesssim h_E$ for each pair $E \in \mathcal{T}_h$ and $T_E \in \mathcal{T}_{\#}$, with $E \subset T_E$.

Moreover, there exists a covering $\mathcal{F}_{\#}$ of $\Gamma_h$ and a positive constant $O_{\Gamma}$, independent of the mesh parameters, such that

and $h_{T_F} \lesssim h_F$ for each pair $F \in \Gamma_h$ and $T_F \in \mathcal{F}_{\#}$, with $F \subset T_F$.

We can now state the following approximation result:

Lemma 6.2. ^{[36,37,38,40,41]} Let $E\in \mathcal{T}_h$ and $T_E\in \mathcal{T}_{\#}$ denote the corresponding simplex such that $E \subset T_E$ (see Definition 6.1). Suppose that $v \in L^2(\Omega)$ is such that $\mathscr{E} v|_{T_E} \in H^{r_E}(T_E)$, for some $r_E\geq0$. Then, if Assumption 5.1 and 6.1 are satisfied, there exists $\widetilde{\Pi}v$, such that $\widetilde{\Pi}v|_E \in \mathbb{P}_{k_E}(E)$, and the following bound holds

Moreover, if $r_E > 1/2$,

Here, $s_E = \min(k_E+1, r_E)$ and the hidden constants depend on the shape-regularity of $T_E$, but are independent of $v$, $h_E$, $k_E$ and the number of faces per element and $\mathscr{E}$ is the continuous extension operator as defined in ^[72].

Proof. See ^[36] for a detailed proof of (6.1) and ^[38] for the proof of (6.2). Clearly, analogous approximation results can be stated on the fracture spaces, since Assumptions 5.1 and 6.1 are both satisfied.

6.2.   Error estimates

For each subdomain $\Omega_i$, $i = 1, 2$, we denote by $\mathscr{E}_i$ the classical continuous extension operator (cf. ^[72], see also ^[59]) $\mathscr{E}_i: H^s(\Omega_i) \rightarrow H^s(\mathbb{R}^d)$, for $s \in \mathbb{N}_0$. Similarly, we denote by $\mathscr{E}_\Gamma$ the continuous extension operator $\mathscr{E}_\Gamma: H^s(\Gamma) \rightarrow H^s(\mathbb{R}^{d-1})$, for $s \in \mathbb{N}_0$. We then make the following regularity assumptions for the exact solution $(p, p_\Gamma)$ of problem (3.1):

Assumption 6.3. Let $\mathcal{T}_{\#} = \{T_E\}$ and $\mathcal{F}_{\#} = \{T_F\}$ denote the associated coverings of $\Omega$ and $\Gamma$, respectively, of Definition 6.1. We assume that the exact solution $(p, p_{\Gamma})$ is such that:

A1. for every $E \in \mathcal{T}_h$, if $E \subset \Omega_i$, it holds $\mathscr{E}_i p_i|_{T_E} \in H^{r_E}(T_E)$, with $r_E \geq 1+d/2$ and $T_E \in \mathcal{T}_{\#}$ with $E \subset T_E$;

A2. for every $F \in \Gamma_h$, it holds $\mathscr{E}_\Gamma p_\Gamma|_{T_F} \in H^{r_F}(T_F)$, with $r_F \geq 1+(d-1)/2$ and $T_F \in \mathcal{F}_{\#}$ with $F \subset T_F$.

Assumption 6.4. We assume that the normal components of the exact fluxes $\boldsymbol{\nu} \nabla p$ and $\ell_\Gamma \boldsymbol{\nu}_\Gamma^\tau \nabla p_\Gamma$ are continuous across mesh interfaces, that is $\left[\kern-0.15em\left[ { \boldsymbol{\nu} \nabla p } \right]\kern-0.15em\right] = 0$ on $\mathcal{F}_h^I$ and $\left[\kern-0.15em\left[ { \ell_\Gamma \boldsymbol{\nu}_\Gamma^\tau \nabla p_\Gamma } \right]\kern-0.15em\right] = 0$ on $\mathcal{E}_{\Gamma, h}^I$.

From Proposition 5.11 and Strang's second Lemma the following abstract error bound directly follows.

Lemma 6.5. Assuming that the hypotheses of Proposition 5.11 are satisfied, for all the discrete problems PP (4.17), MP (4.31), MM (4.41) and PM (4.41) the following abstract error bound holds

where the residual $\mathcal{R}_h^{**}$ is defined as

with $** \in \{PP, MP, MM, PM\}$.

Note that, irrespective of the numerical method chosen for the discretization (PP, MP, PM or MM), the residual can always be split into two contributions, one deriving from the approximation of the problem in the bulk and one deriving from the approximation of the problem in the fracture, i.e.,

It follows that, to derive a bound for the global residual, we can bound each of the two contributions separately. With this in mind, we state and prove the next two lemmas.

Lemma 6.6. Let $(p, p_\Gamma)$ be the exact solution of problem (3.1) satisfying the regularity Assumptions 6.4 and 6.3. Then, for every $w \in Q^b(h)$ and $w_\Gamma \in Q^\Gamma(h)$, it holds

Proof. First, we prove (6.4). Let $\boldsymbol{\Pi}_{\textbf{W}_h^b}$ be the $L^2$-orthogonal projector onto $\textbf{W}_h^b$, then, integrating by parts elementwise, using the fact that $p$ satisfies (2.2) and recalling that, from Assumption 6.4, $\left[\kern-0.15em\left[ { \boldsymbol{\nu} \nabla p } \right]\kern-0.15em\right]$ vanishes on $\mathcal{F}_h^I$, we obtain the following expression for the residual $\mathcal{R}_b^P$:

Employing the Cauchy-Schwarz's inequality and the definition of the norm $|||\cdot|||$, we then obtain

If we still denote by $\widetilde{\Pi}$ the vector-valued generalization of the projection operator $\widetilde{\Pi}$ defined in Lemma 6.2, we observe that

To bound the term (1), we employ the approximation result stated in Lemma 6.2. We obtain

Exploiting, the boundedness of the permeability tensor $\boldsymbol{\nu}$, the inverse inequality (5.1), the $L^2$-stability of the projector $\boldsymbol{\Pi}_{\textbf{W}_h^b}$ and the approximation results stated in Lemma 6.2, we can bound term (2) as follows:

which concludes the proof of (6.4).

Proceeding as above we obtain the following expression for the residual $\mathcal{R}_\Gamma^P$:

Estimate (6.5) can then be proven with analogous arguments.

Lemma 6.7. Let $(p, p_\Gamma)$ be the exact solution of problem (3.1) satisfying the regularity Assumptions 6.4 and 6.3. Then, for every $w \in Q^b(h)$ and $w_\Gamma \in Q^\Gamma(h)$, it holds

Proof. We focus on the proof of (6.6), estimate (6.7) can be obtained likewise. Recalling that $\boldsymbol{\Pi}_{\textbf{W}_h^b}$ denotes the $L^2$-orthogonal projector onto $\textbf{W}_h^b$, the residual $\mathcal{R}_b^M$ has the following expression:

where we have used the identity $\mathscr{L}_b^{ LDG}(\left[\kern-0.15em\left[ { p } \right]\kern-0.15em\right]) = -\mathcal{G}_b$ and the continuity of $\boldsymbol{\nu} \nabla p$ across internal faces (Assumption 6.4). From Cauchy-Schwarz and triangular inequalities and the bound on the coefficient $\textbf{b}$ (4.23), we have

where we recall that, with a slight abuse of notation, $\widetilde{\Pi}$ still denotes the vector-valued generalization of the projection operator $\widetilde{\Pi}$ defined in Lemma 6.2. The thesis now follows from the boundedness of the permeability tensor $\boldsymbol{\nu}$, the inverse inequality (5.1), the $L^2$-stability of the projector $\boldsymbol{\Pi}_{\textbf{W}_h^b}$ and the approximation results stated in Lemma 6.2.

Theorem 6.8. Let $\mathcal{T}_{\#} = \{T_E\}$ and $\mathcal{F}_{\#} = \{T_F\}$ denote the associated coverings of $\Omega$ and $\Gamma$, respectively, consisting of shape-regular simplexes as in Definition 4, satisfying Assumptions 6.1. Let $(p, p_{\Gamma})$ be the solution of problem (3.1) and $(p_h, p_{\Gamma, h}) \in Q_h^b \times Q_h^{\Gamma}$ be its approximation obtained with the method PP, MP, MM or PM, with the penalization parameters given by (5.2) and (5.3) and $\sigma_0$ and $\sigma_{0, \Gamma}$ sufficiently large for the primal formulations. Moreover, suppose that the exact solution $(p, p_{\Gamma})$ satisfies the regularity Assumptions 6.4 and 6.3. Then, the following error bound holds:

where the $\mathscr{E}p$ is to be interpreted as $\mathscr{E}_1 p_1$ when $E \subset \Omega_1$ or as $\mathscr{E}_2 p_2$ when $E \subset \Omega_2$. Here, $s_E = \min(k_E+1, r_E)$ and $s_F = \min(k_F+1, r_F)$, and the constants are defined according to the chosen approximation method as follows:

Proof. From Lemma 6.5 we know that the error satisfies the following bound

We estimate the two terms on the right-hand side of (6.8) separately. We can rewrite term $I$ as

Again we consider each of the three terms separately. To bound term ($a$), we exploit the two approximation results stated in Lemma 6.2 and obtain

Using analogous interpolation estimates on the fracture we can bound term ($b$) as follows:

Finally, for term ($c$), we have

Exploiting the approximation result (6.2), we obtain

Similarly, we have

Finally, using the interpolation estimates for the fracture, we deduce that

In conclusion, combining all the previous estimates, we can bound the term $I$ on the right-hand side of (6.8) as follows:

Finally, the desired estimates follow from the combination of (6.9), together with the bound on Term $II$ deriving from what observed in (6.3) and Lemmas 6.6 and 6.7.

Finally, from the above result we can derive some error estimates also for the velocities $\textbf{u}$ and $\textbf{u}_\Gamma$.

Theorem 6.9. Let all the hypotheses of Theorem 6.8 hold. Let $(\textbf{u}, \textbf{u}_\Gamma) \in \textbf{W}_g$ and $(p, p_\Gamma) \in M$ be the solution of problem (3.1). Then:

● if $\big((p_h, \textbf{u}_h), p_{\Gamma, h} \big) \in Q_h^b \times \textbf{W}_h^b \times Q_h^\Gamma$ is its approximation obtained with the MP method (4.26), it holds

● if $\big(p_h, (p_{\Gamma, h}, \textbf{u}_{\Gamma, h})\big) \in Q_h^b \times Q_h^\Gamma \times \textbf{W}_h^\Gamma$ is its approximation obtained with the PM method (4.35), it holds

● if $\big((p_h, \textbf{u}_h), (p_{\Gamma, h}, \textbf{u}_{\Gamma, h})\big) \in Q_h^b \times \textbf{W}_h^b \times Q_h^\Gamma \times \textbf{W}_h^\Gamma$ is its approximation obtained with the MM method (4.40), it holds

where the constants $G_E^{M}$, $G_F^P$, $G_E^P$ and $G_F^M$ are defined as in Theorem 6.8.

Proof. Let $\big((p_h, \textbf{u}_h), p_{\Gamma, h}\big)$ and $\big((p_h, p_{\Gamma, h}), (\textbf{u}_h, \textbf{u}_{\Gamma, h})\big)$ be the discrete solutions obtained with the MP method and with the MM method, respectively. Then, using identity (4.29) and the fact that $\mathscr{L}_b^{ LDG}(\left[\kern-0.15em\left[ { p } \right]\kern-0.15em\right]) = -\mathcal{G}_b$, we can rewrite

From the uniform boundedness of $\boldsymbol{\nu}$, the triangular inequality, the bound on the lifting (5.6) and the definition of the $|| \cdot ||_{DG}$ norm it follows that

In particular, this implies that $||\textbf{u}-\textbf{u}_h||_{0, \mathcal{T}_h} \lesssim |||(p, p_{\Gamma})-(p_h, p_{\Gamma, h})|||$. Similarly, one can prove that, if $\big(p_h, (p_{\Gamma, h}, \textbf{u}_{\Gamma, h})\big)$ and $\big((p_h, p_{\Gamma, h}), (\textbf{u}_h, \textbf{u}_{\Gamma, h})\big)$ are the discrete solutions obtained with the PM method and with the MM method, respectively, it holds $||\textbf{u}_\Gamma-\textbf{u}_{\Gamma, h}||_{0, \Gamma_h} \lesssim |||(p, p_{\Gamma})-(p_h, p_{\Gamma, h})|||.$ The thesis is now a direct consequence of Theorem 6.8.

7.   Numerical experiments

In this section we present some two-dimensional numerical experiments with the aim of validating the obtained theoretical convergence results and of comparing the accuracy of the proposed formulations. In particular, we will focus on the paradigmatic primal-primal and mixed-primal settings. This means that, for the approximation of the problem in the bulk, we will employ the SIPDG method in the first setting and the LDG method in the second setting, while, for the problem in the fracture, we will employ the SIPDG method in both settings. All the numerical tests have been implemented in $\texttt{Matlab}^{®}$ . For the generation of polygonal meshes conforming to the fractures we have suitably modified the code $\texttt{PolyMesher}$ ^[73].

In particular, we present three sets of numerical experiments. The first set (Sections 7.1 and 7.2) is obtained assuming that analytical solutions are known and aims at verifying the a-priori error estimates obtained in Theorems 6.8 and 6.9. The second set (Section 7.3) is derived from physical considerations and aims at testing how different values of the fracture permeability may influence the flow in the bulk. Finally, the last set of experiments (Section 7.4) aims at showing how the method is capable of handling more complicated geometries, specifically networks of partially immersed fractures.

7.1.   Example 1: Analytical solution with constant fracture pressure

In this first test case we take $\Omega = (0, 1)^2$ and choose as exact solutions in the bulk and in the fracture $\Gamma = \{(x, y) \in \Omega: \; \; x+y = 1\}$

It is easy to prove that $\textbf{u}$, $p$ and $p_{\Gamma}$ satisfy the coupling conditions (2.5a)-(2.5b) with $\xi = 1$, $\ell_\Gamma = 0.001$ and $\boldsymbol{\nu} = \boldsymbol{\nu}_\Gamma = \textbf{I}$. Note that in this case $f_{\Gamma} = 0$ since the solution in the fracture is constant and $\left[\kern-0.15em\left[ { \textbf{u} } \right]\kern-0.15em\right] = 0$.

Figure 2 shows three successive levels of refinements for the polygonal mesh employed in this set of experiments. In order to test the behaviour of the energy norm of the error, thus validating the $h$-convergence properties of our methods proved in Theorem 6.8, we compute the quantity ${||p-p_h||_{1, \mathcal{T}_h} + ||p_{\Gamma}- p_{\Gamma, h}||_{1, \Gamma_h}}$ (Figure 3(a) and 3(b)). We also want to validate the results of Theorem 6.9 for the velocity, computing $||\textbf{u}-\textbf{u}_h||_{L^2(\Omega)}$ (Figure 3(e)). In addition, we test the behaviour of the $L^2$-norm of the error for the primal unknowns, i.e., $||p-p_h||_{L^2(\Omega)} + ||p_{\Gamma}- p_{\Gamma, h}||_{L^2(\Gamma)}$ (Figure 3(c) and 3(d)). All the plots in Figure 3 show the computed errors as a function of the inverse of the mesh size $h$ (loglog scale). Each plot consists of four lines: Every line shows the behaviour of the computed error for a different polynomial degree in the bulk (we consider uniform degree $k_E = k = 1, 2, 3, 4$ for all $E \in \mathcal{T}_h$). For the fracture problem we always choose uniform quadratic polynomial degree, i.e., $k_F = k_\Gamma = 2$ for all $F \in \Gamma_h$. On the left we report the results obtained with the (MP) approximation scheme and on the right with the (PP) scheme. We observe that, for both methods, the convergence rates are superoptimal with respect to the expected rate of $min(k, k_\Gamma)$ (they coincide with the bulk polynomial degree $k$). This is probably due to the constant nature of the solution in the fracture. Indeed this behaviour will not be observed in the next set of experiments, cf. Section 7.2, where the solution in the fracture is trigonometric. Moreover, Figure 3(c) and 3(d)) show that one order of convergence is gained for the $L^2$-norm. Finally, one can clearly see that the level of accuracy achieved by the (PP) and (MP) schemes is the same.

7.2.   Example 2: Analytical solution with non-constant fracture pressure

Next, we consider the domain $\Omega = (0, 1)^2$ and the fracture $\Gamma = \{(x, y) \in \Omega: \; \; x = 0.5\}$. We reproduce the numerical test already proposed in ^[59] for the primal-primal setting, choosing the exact solutions in the bulk and in the fracture as follows

We impose Dirichlet boundary conditions on the whole $\partial \Omega$ and also on $\partial \Gamma$. Notice that coupling conditions (2.5a)-(2.5b) are satisfied provided that we take $\boldsymbol{\nu} = \textbf{I}$, and $\beta_\Gamma = 2$, that is $\boldsymbol{\nu}^n_\Gamma / \ell_\Gamma = 4$. The source terms are then chosen accordingly as

In the experiments, we set the components of the fracture permeability tensor $\boldsymbol{\nu}_{\Gamma}^{\tau} = 10^2$ and $\nu_{\Gamma}^n = 4 \cdot 10^{-2}$, we set the fracture thickness $\ell_\Gamma = 10^{-2}$ and the closure parameter $\xi = \frac{3}{4}$. As before, in order to test the $h$-convergence properties of the primal-primal and mixed-primal schemes, we compute the quantity ${||p-p_h||_{1, \mathcal{T}_h} + ||p_{\Gamma}- p_{\Gamma, h}||_{1, \Gamma_h}}$. In Figure 4 we show the computed errors as a function of the inverse of the mesh size $h$ (loglog scale), together with the expected convergence rates. As before, we fix the polynomial degree for the fracture problem $k_\Gamma = 2$, and we vary the polynomial degree for the problem in the bulk taking $k = 1, 2, 3, 4$. Figure 4(a) encloses the results obtained with the mixed-primal approximation, while Figure 4(b) with the primal-primal approximation. In each plot, the four lines describe the behaviour of the energy norm of the error for a different polynomial degree in the bulk. Notice that in this case, for both the (PP) and (MP) method, the theoretical convergence rates are clearly obtained, coinciding with $\min(k, k_\Gamma)$ (no superconvergence is observed). In particular, the convergence rate is equal to $1$ when the approximation in the bulk is linear, i.e., when $k = 1$ and it is equal to $2$ in all the other cases, since we always choose a quadratic approximation for the problem in the fracture. Moreover, we remark that also in this case the (PP) and (MP) methods achieve the same level of accuracy.

7.3.   Example 3: Discontinuous fracture permeability

Next, we reproduce some numerical experiments first presented in ^[47]. We examine two test cases with bulk domain $\Omega = (0, 2) \times (0, 1)$ and fracture domain $\Gamma = \{(x, y) \in \mathbb{R}^2: x = 1, \, 0\leq y \leq 1\}$. In the first case, we consider a fracture with constant permeability, while in the second case we consider a fracture with lower permeability in its middle part, thus presenting a discontinuity. In particular:

(a) Case 1: constant permeability: The permeability tensor in the fracture is given by $\boldsymbol{\nu}_{\Gamma}^n = \boldsymbol{\nu}_{\Gamma}^\tau = 100$. The bulk permeability $\boldsymbol{\nu}$ is chosen to be constant and isotropic, i.e., $\boldsymbol{\nu} = \textbf{I}$. We impose Dirichlet boundary conditions on the left and right side of the bulk domain and homogeneous Neumann conditions on the top and bottom sides. On the fracture boundaries we impose Dirichlet boundary conditions.

(b) Case 2: discontinous permeabilty: The fracture $\Gamma$ is subdivided into two areas having different values for the permeability tensor: In the initial and ending part of the fracture $\Gamma_1 = \{ (x, y) \in \Gamma, \; 0\leq y \leq 0.25 \; \mbox{and} \; 0.75 \leq y \leq 1\}$ the permeability tensor $\boldsymbol{\nu}_{\Gamma_1}$ is defined as $\boldsymbol{\nu}_{\Gamma_1}^n = \boldsymbol{\nu}_{\Gamma_1}^\tau = 1$, while in the middle part $\Gamma_2 = \{ (x, y) \in \Gamma, 0.25\leq y \leq 0.75 \}$ the permeability is low and is defined as $\boldsymbol{\nu}_{\Gamma_2}^n = \boldsymbol{\nu}_{\Gamma_2}^\tau = 0.002$. The bulk permeability tensor is chosen again equal to the identity matrix, i.e., $\boldsymbol{\nu} = \textbf{I}$. In the bulk, we impose the same boundary conditions as in the previous test case, while at the fracture extremities we impose homogeneous Neumann conditions.

The two geometrical configurations are shown in Figures 5(a)-5(b), together with the boundary conditions. For both test cases we take the fracture thickness $\ell_\Gamma = 0.01$ and the model parameter $\xi = 2/3$. Moreover, we discretize the problem in the bulk taking as polynomial degree $k = 1$ and the problem in the fracture taking $k_\Gamma = 2$.

The obtained results are shown in Figure 6. For both cases (constant at the top, discontinuous at the bottom) we report the pressure field and Darcy velocity in the bulk (here the grid is very coarse only for visualization purposes) and the value of the pressure along the fracture. In the first case, since we have taken $\boldsymbol{\nu}_{\Gamma}^n = \boldsymbol{\nu}_{\Gamma}^\tau = 100 > 1$, we can observe that the fluid has the tendency to flow along the fracture. In the second case, one can see that the part of the fracture with low (normal) permeability acts as a geological barrier, so that the fluid tends to avoid it and we can observe a jump of the bulk pressure across it. Our results are in agreement with those obtained in ^[47].

7.4.   Example 4: Network of partially immersed fractures

With this last set of numerical experiments we investigate the capability of our discretization method to deal with more complicated geometrical configurations, considering a network of partially immersed fractures. Our reference is, in this case, ^[51], where the mathematical model developed in ^[47] has been extended to fully immersed fractures. In ^[59] we showed that our method in a primal-primal setting is capable of efficiently handling the configuration. Here, we reproduce the same numerical experiments to demonstrate that this holds true also in a mixed-primal setting.

In order to deal with immersed fractures, we need to supplement our model (2.6) with an equation describing the behaviour of the fracture pressure at the immersed tips. Following ^[51], we impose a homogeneous Neumann condition, thus assuming that the mass transfer across the immersed tips can be neglected, i.e., $\boldsymbol{\nu}_\Gamma^\tau \nabla_\tau p_\Gamma \cdot \boldsymbol{\tau} = 0$ on $\partial \Gamma$. At the extremities of the fractures that are non-immersed, i.e., $\partial \Gamma \cap \partial \Omega$, we impose boundary conditions that are consistent with those imposed on $\partial\Omega$ in that point.

We consider the bulk domain $\Omega = [0, 1]^2$ cut by a network made of four partially immersed fractures: $\Gamma_1 = \{(x, y) \in [0, 1]^2: x \geq 0.3, y = 0.2 \}$, $\Gamma_2 = \{(x, y) \in [0, 1]^2: x\leq 0.7, y = 0.4 \}$, $\Gamma_3 = \{(x, y) \in [0, 1]^2: x\geq 0.3, y = 0.6 \}$ and $\Gamma_4 = \{(x, y) \in [0, 1]^2: x\leq 0.7, y = 0.8 \}$. We perform two numerical experiments. In both of them, the fractures $\Gamma_2$ and $\Gamma_4$ are impermeable ($\boldsymbol{\nu}_\Gamma^\tau = \nu_\Gamma^n = 10^{-2}$), while $\Gamma_1$ and $\Gamma_3$ are partially permeable. In the first configuration, we consider for $\Gamma_1$ and $\Gamma_3$ the permeabilities $\nu_\Gamma^n = 10^{-2}$ and $\boldsymbol{\nu}_\Gamma^\tau = 100$, while in the second, we consider $\nu_\Gamma^n = 10^{-2}$ and $\boldsymbol{\nu}_\Gamma^\tau = 1$. Moreover, we vary the imposed boundary conditions as illustrated in Figure 7.

In both the experiments we consider an isotropic bulk permeability tensor i.e., $\boldsymbol{\nu} = \textbf{I}$ and we assume that all the fractures have aperture $\ell_\Gamma = 0.01$. The flow is only generated by boundary conditions, since we take all the forcing terms $f = f_\Gamma = 0$. Finally, we choose as model parameter $\xi = 0.55$.

To obtain our results, we employed cartesian grids featuring approximately the same number of elements as those employed in ^[51] and such that the immersed tips of the fractures coincide with one of the mesh vertices. For the approximation of the problem in the bulk and in the fracture we chose the polynomial degrees $k = k_\Gamma = 2$. In Figure 8, we show the results obtained for the two test cases with a mesh of $26 051$ elements. In particular, we report the pressure field in the bulk with the streamlines of the velocity (left), the value of the bulk pressure along the line $x = 0.65$ (middle) and the pressure field inside the four fractures (right). Our results are in perfect agreement with those obtained in ^[51] and in ^[59], thus showing that, also in a mixed-primal setting, our method is able to efficiently handle this configuration.

8.   Conclusions

In this paper we have proposed a formulation based on discontinuous Galerkin methods on polygonal/polyhedral grids for the simulation of flows in fractured porous media. In particular, we have designed and analysed, in the unified framework of ^[63] based on the flux-formulation, a polyDG approximation for all the possible combinations of primal-primal, mixed-primal, primal-mixed and mixed-mixed formulations for the bulk and fracture problems, respectively. The novelty of our method relies on the imposition of coupling conditions between bulk and fracture through a suitable definition of the numerical fluxes on the fracture faces. We have proved in an unified setting the well-posedness of all the formulations and we have derived a priori $hp$-version error estimates in a suitable (mesh-dependent) energy norm, whose validity has been assessed performing some preliminary numerical experiments, focusing on the paradigmatic primal-primal and mixed-primal methods. In our test cases we have also compared, in a simplified setting, the performance of our approximation schemes. In particular, we have shown that the same level of accuracy may be obtained irrespective of the chosen method. The other the factors that should be taken into account when choosing which one between the (PP), (MP), (PM), and (MM) setting to employ, are summarized in the following.

● The desired accuracy in the approximation of the physical quantities (pressure and velocity) according to the application at hand. The primal approach considers a single-field formulation with the pressure field of the fluid as only unknown, so that velocity may only be afterwards reconstructed taking the gradient of the computed pressure and multiplying it by the permeability tensor. This process usually entails a loss of accuracy and does not guarantee mass conservation, see for example ^[61,62], so it may not be suitable for those Engineering applications (such as oil recovery or groundwater pollution modeling) where the simulation of the phenomenon requires very accurate approximation of the velocities of the involved fluids in order to be effective. In such cases, the mixed approach should be preferred, so that Darcy's velocity can be directly computed and a higher degree of accuracy can be achieved, together with local and global mass conservation.

● The numerical linear algebra resulting from the discretization. The primal setting has the advantage of featuring less degrees of freedom and leads to a symmetric positive definite algebraic system of equations, which can be efficiently solved based on employing, for example, multigrid techniques ^[39,44,60]. The mixed approach leads to a so-called generalized saddle point system of equations. For example, in the MP setting, problem (4.26) translates into an algebraic system of the form

where, referring to Eq. (4.27):

- $M$ is the mass matrix related to the bilinear form $\mathcal{M}_b$;

- $B$ is related to the bilinear form $\mathcal{B}_b$;

- $S$ is related to the bilinear form $\mathcal{S}_b$;

- $C_1$ is related to the terms involving only bulk unknowns contained in the interface bilinear forms $\mathcal{I}_1$ and $\mathcal{I}_2$;

- $C_2$ is related to the terms involving both bulk and fracture unknowns in $\mathcal{I}_1$ and $\mathcal{I}_2$;

- $A_\Gamma$ is related to the fracture primal bilinear form $\mathcal{A}_{\Gamma}^{ P}$;

- $A_\Gamma^R$ is related to the terms involving only fracture unknowns contained in the interface bilinear forms $\mathcal{I}_1$ and $\mathcal{I}_2$, which result in a reaction term for the fracture problem.

In this case the Schur complement of $M$, i.e. the matrix $D+B^TM^{-1}B$, with $D = \begin{bmatrix} S +C_1 & C_2\\ C_2^T & A_\Gamma + A_\Gamma^R \end{bmatrix}$, can be computed explicitly, thanks to the fact that, in the DG setting, $M$ either has a block-diagonal or a diagonal structure, depending on whether a modal or nodal expansion is chosen to span the discrete space. We also point out that the mixed case involves a (three-time) larger number of unknowns, so that efficient memory handling may be needed.

● The conditioning of the resulting systems. The behaviour of the condition number and the possibility of building efficient preconditioners is certainly of fundamental importance in the choice between the primal or mixed approach, especially for three-dimensional problems. We did not address this issue in the present work. However, one may refer to ^[76], where preconditioners for the system stemming from the discretization of problem (2.6) in mixed-mixed form, employing mimetic finite differences, are constructed. In the DG setting, for the $hp$-preconditioning of systems stemming from the discretization of diffusion problems, the reader may refer to ^[71] for standard grids and to ^[39,44,78] for the extension to the polytopic setting.

Finally, we mention that our method may be extended in order to deal with network of intersecting fractures. To this aim, the mathematical model needs to be complemented with some suitable physical conditions at the intersections, prescribing the behaviour of the fluid. One possible choice is to impose pressure continuity and flux conservation as in ^[52,74]. From the DG-discretization point of view, the key point to deal with this case is the generalization of the concepts of jump and average at the intersections. We refer to ^[75] for a rigorous analysis of the method in the primal-primal setting and to ^[59] for a preliminary numerical test case involving a totally immersed network of fractures.

Acknowledgments

Paola F. Antonietti and Chiara Facciolà have been partially supported by SIR Starting grant n. RBSI14VT0S "PolyPDEs: Non-conforming polyhedral finite element methods for the approximation of partial differential equations" funded by MIUR. Paola F. Antonietti and Marco Verani have been partially funded by the Italian research project PRIN n. 201744KLJL. All the authors have been partially supported by INdAM-GNCS.

Conflict of interest

The authors declare no conflict of interest.