A posteriori error estimates for hierarchical mixed-dimensional elliptic equations

Mixed-dimensional elliptic equations exhibiting a hierarchical structure are commonly used to model problems with high aspect ratio inclusions, such as flow in fractured porous media. We derive general abstract estimates based on the theory of functional a posteriori error estimates, for which guaranteed upper bounds for the primal and dual variables and two-sided bounds for the primal-dual pair are obtained. We improve on the abstract results obtained with the functional approach by proposing four different ways of estimating the residual errors based on the extent the approximate solution has conservation properties, i.e.: (1) no conservation, (2) subdomain conservation, (3) grid-level conservation, and (4) exact conservation. This treatment results in sharper and fully computable estimates when mass is conserved either at the grid level or exactly, with a comparable structure to those obtained from grid-based a posteriori techniques. We demonstrate the practical effectiveness of our theoretical results through numerical experiments using four different discretization methods for synthetic problems and applications based on benchmarks of flow in fractured porous media.

while the two remaining correspond to 3d benchmark problems [14].

Introduction
Mixed-dimensional partial differential equations (mD-PDEs) arise when partial differential equations interact on domains of different topological dimensions [1].Prototypical examples include models of thin inclusions in elastic materials [2,3,4], blood flow in human vasculature [5,6,7], root water uptake systems [8], and flow in fractured porous media [9,10,11].The latter example has an appealing mathematical structure, in that the model equations allow for a hierarchical representation where each subdomain (matrix, fractures, fracture intersections, and intersection points) only has direct interaction with subdomains of topological dimension one higher or one lower [12].Such hierarchical mD-PDEs are the topic of the current paper.
mD-PDEs are intrinsically linked to the underlying geometric representation, which, in a certain sense, generalizes the usual notion of the domain.One can then define sets of suitable functions (and function spaces) on this geometry, and these sets are then naturally interpreted as mixed-dimensional (mD) functions.Exploiting this concept, one can generalize the standard differential operators to mappings between mD functions and thus obtain an mD calculus.The fact that this mD calculus inherits standard properties of calculus, particularly partial integration (relative to suitable inner products), a de Rham complex structure, and a Poincaré-Friedrichs inequality, was recently established using the language of exterior calculus on differential forms [15].
The inherent geometric generality of hierarchical mD-PDEs also demand the same level of abstraction of a posteriori error estimation techniques.This requirement makes error estimates of the functional type particularly well-suited for the task [16,17,18,19,20,21].The most attractive feature of this approach is that error estimates are derived using purely functional methods [20].The bounds are therefore agnostic to the way approximated solutions are obtained in the energy space, and the only undetermined constants arise from Poincaré-type inequalities [22].
However, unlike other types of error estimates [23,24,25,26,27], this generality makes standard functional estimates of limited applicability to hierarchical elliptic mD-PDEs due to the following reasons: (1) for general fracture networks, the mixed-dimensional Poincaré constant is not easily computable, and (2) since Poincaré constants are proportional to the diameter of the physical domain, residual estimators cannot exhibit superconvergent properties.
To circumvent the aforementioned issues, we exploit the fact that Poincarétype inequalities imply weighted norms [28,29], and use spatially-dependent weights to control the residual norms.We show both theoretically and numerically that this treatment leads to sharper estimates when approximations to the exact solution satisfy mass conservation in a given partition of the domain.
In view of the preceding discussion, our aim is therefore to obtain a posteriori error estimates for the approximate solution to the mD scalar elliptic equation [12,30,15], where the mD Laplace equation for geometries such as those illustrated in Figure 1 is described in detail in Section 3.
We remark that while a broad range of a posteriori error techniques are available for mono-dimensional problems, existing error bounds for mD models are far more scarce.Moreover, the ones available, are restricted to specific cases (e.g., in the context of mortar methods [31,32,33,34] and fractured porous media [35,36,37]) with far less geometric generality than what we present here.Thus, for practical problems, a posteriori error bounds for mD geometries have until now essentially not been available.
The rest of the paper is structured as follows: Section 2 is devoted to introducing the model problem, functional spaces, and variational formulations for the case of a single 1d fracture embedded in a 2d matrix.The section is concluded by providing a first upper bound for the primal variable.In Section 3, we generalize the results from Section 2 to the case of fracture networks and introduce the necessary tools to perform the a posteriori analysis in an mD setting.After reviewing necessary tools from functional analysis in Section 4, in Section 5, we provide our main results starting from a generic abstract estimate and then considering specific cases depending upon the degree of accuracy at which residual terms are approximated.In Section 6, we introduce the approximated problem using mixed-finite element methods and thus make the estimates concrete.Sections 7 and 8 deal, respectively, with numerical validations and practical applications of the derived bounds.Finally, in Section 9, we present our concluding remarks.

Upper bounds for a single fracture
In this section, we introduce the model problem together with functional spaces and the variational formulations for the case of a single 1d line embedded in a 2d matrix, as illustrated in Figure 2. Furthermore, a first upper bound for the primal variable is derived following the classical functional approach.We remark that the case of a single fracture embedded in a matrix has been analyzed before.For example, [35] and [37] proposed error estimators based on the residual approach, whereas [36] obtained guaranteed a posteriori error estimates using the approach of Vohralík [26].
Figure 2: A horizontal 1d fracture embedded in a 2d matrix.Left: Subdomains and interfaces.Right: Boundary conditions.For the fracture, the purple square denotes a no-flux boundary condition, whereas the green square a Dirichlet boundary condition.Note that ∂ 1 Ω 2 , Γ 1 , Ω 1 , Γ 2 , ∂ 2 Ω 2 , all coincide spatially.For illustrative purposes, however, they are placed in different locations.

The model problem for a single fracture
Before writing the set of equations describing general fracture networks, let us first introduce the governing equations of a simpler configuration; that is, a unit square domain Y ⊂ R 2 decomposed as a 1d fracture Ω 1 embedded in a 2d matrix Ω 2 as shown in the left Figure 2. Interfaces Γ 1 and Γ 2 , at each side of Ω 1 , establish the link between Ω 2 and Ω 1 .The model presented below is well-established for these problems, and we point the reader to the references for further justification of this system [38,12,30].
The strong form of the governing equations in Ω 2 reads Here, (1a) is the mass conservation equation, u 2 is the matrix velocity, and f 2 an external source.The fluid velocity is given by the standard Darcy's law (1b), where K 2 is the matrix permeability; a bounded, symmetric, and positive-definite 2 × 2 tensor, and p 2 is the fluid pressure.Equations (1c) and (1d) require that at each side of the internal boundary of Ω 2 , the normal component of u 2 to match the interface (mortar) fluxes λ 1 and λ 2 .To fix the direction of the normal vector on internal boundaries, we require n 2 pointing from the higher-to the lower-dimensional subdomain.No flux conditions are prescribed in (1e), where u 2 • n 2 represents the outer normal flux across ∂ N Ω 2 .Finally, Dirichlet boundary conditions are imposed in (1f), where g D,2 is a prescribed function on the Dirichlet boundary.
In the fracture Ω 1 , the equations are given by In (2a), are the divergence and gradient operators acting in the tangent space of Ω 1 , u 1 is the tangential fracture velocity, the term in parentheses represents the jump in normal fluxes from the adjacent interfaces Γ 1 and Γ 2 onto Ω 1 , and f 1 is an external source.
The tangential velocity u 1 is again expressed via Darcy's law (2b), where in a slight abuse of notation, we use K 1 to refer to the tangential component of the fracture permeability, which is again assumed to be positive and bounded from above.Finally, (2c) and (2d) are the Neumann and Dirichlet boundary conditions, respectively.Again, we use g D,1 to denote a prescribed function on the Dirichlet part of the fracture boundary.
To close the system of equations, we must specify a constitutive relationship for the interface fluxes.Here, we use a Darcy-type law [38], where mortar fluxes are linearly related to pressure jumps with κ 1 and κ 2 representing the effective normal permeability on Γ 1 and Γ 2 , respectively.We restrict our analysis to the case where κ 1 and κ 2 are nondegenerate.Thus, following [12], we further require the existence of two constants γ 1 and γ 2 such that 0 < γ 1 ≤ κ −1 j ≤ γ 2 < ∞ for j ∈ {1, 2}.

Functional spaces and variational formulations
Let us now present the primal weak formulation of the single fracture model from the previous section.To this aim, consider first the energy space with vanishing traces on Dirichlet boundaries and the product spaces Furthermore, let •, • Ωi and •, • Γj denote respectively the L 2 -inner products on Ω i and Γ j , and • Ωi and • Γj the relevant L 2 -norms.Finally, we denote by g = [g 1 , g 2 ] ∈ H 1 (Ω) two functions extending the boundary data into the domains, and thus satisfying tr ∂D Ωi g i = g D,i .We now state the primal weak problem as: Definition 1 (Primal weak formulation for a single fracture).
Refer to Appendix A.1 for the derivation of the primal weak form from the strong form in Section 2.1.We see directly from equation ( 6) that the primal weak form has a minimization structure subject to the stated conditions on K i and κ j , and well-posedness follows by standard arguments.
A dual weak form for the model problem, with explicit representation of the subdomain velocities and mortar fluxes, can also be constructed.We first define the space H(div; Ω i , ∂ X Ω) as the space of L 2 -vector functions on Ω i with weak divergence in L 2 (Ω i ) and zero trace on the part of the boundary indicated by ∂ X Ω.Then, we denote the product spaces of H(div)-functions that are zero on Neumann, and on Neumann and internal boundaries as: Furthermore, we define the L 2 -product spaces on the domains: With these spaces in hand, we consider the standard linear extension operators from internal boundaries onto domains denoted The precise choice of the extension operator R j is not important; however, the natural choice based on the solution of an auxiliary elliptic equation is reasonable [12].We naturally extend the definition of The above constructions allow us to represent subdomain fluxes as where u 0 ∈ V 0 and λ ∈ L 2 (Γ).This motivates the construction of a compound H(div)-type spaces, as This construction will become key when we generalize to more complex geometries in the next section.
Remark 1 (On the regularity of H(div; Ω, Γ)).It is worth remarking that the restriction of space H(div; Ω, Γ) to the domain Ω 2 has slightly enhanced regularity relative to the standard space H(div; Ω 2 ), as this latter space has normal traces which do not lie in L 2 (Γ 1 ) nor L 2 (Γ 2 ).
Definition 2 (Dual weak formulation for a single fracture.).
Refer to Appendix A.2 for the derivation.
Remark 2 (Well-posedness).The variational formulation from Definition 2 can be classified as a saddle point structure, for which well-posedness results have been established for fracture networks, see e.g.Theorem 2.5 from [12].

A first a posteriori error estimate for the primal variable
Having the functional spaces and weak formulations formally introduced, in this section, we provide a first upper bound for an approximation to the primal variable q = [q 1 , q 2 ] ∈ H 1 0 (Ω) + g for the case of a single fracture in the energy norm Theorem 1 (A first upper bound for the primal variable).Let p ∈ H 1 0 (Ω) + g be the solution to the primal weak form (6) with ∂ D Ω 1 non-empty.Then for any q ∈ H 1 0 (Ω) + g, it holds that with where C Ω1 and C Ω2 are the permeability-weighted Poincaré-Friedrichs constants for Ω 1 and Ω 2 : Proof.Refer to Appendix B for the proof.
Remark 3 (Nature of the estimators).The upper bound (15) is a guaranteed upper bound for the deviation between the primal solution p ∈ H 1 0 (Ω) + g and an arbitrary approximation q ∈ H 1 0 (Ω) + g in the energy space.There are three types of contributions to the upper bound: (1) diffusive flux estimators (16a) and (16b) measuring the difference between the approximate fluxes v 0 + Rν ∈ V and fluxes obtained from H 1 0 (Ω)-potentials q, (2) domain coupling estimators (16c) and (16d) measuring how close the approximate normal fluxes ν ∈ L 2 (Γ) are to the jump in H 1 0 (Ω)-potentials q, and (3) residual estimators (16e) and (16f) measuring the difference between the exact source term and the divergence of the approximate flux plus the jump in adjacent approximate normal fluxes.An important detail is that the approximate cross-domain fluxes ν 1 and ν 2 enter into the residual estimators of both the higher-and lower-dimensional subdomain.
Remark 4 (Sharpness of the estimates).The estimates above are in principle sharp, as can be shown by standard arguments [20].However, in practice, we will often have access to additional information about the approximate solution (most commonly if it is derived with a local conservation property).This allows for improvements in the residual estimators (16f) and (16e), as we will show in Section 5.2.
It is clear that even for this fairly simple configuration, the variational formulations (and the analysis in general) can be quite cumbersome.The situation escalates in complexity when intersecting fractures (see Figure 3) are part of the geometric configuration, in particular as the proof of Theorem 1 relies on all subdomains having some non-vanishing Dirichlet boundary.Indeed, when floating subdomains (e.g., fully embedded fractures or isolated rock domains) are present in the fracture network, the standard procedure used in Theorem 1 Left: The domain Y is decomposed into two 2d matrices (Ω 9 and Ω 10 ), four 1d fractures (Ω 5 , Ω 6 , Ω 7 , and Ω 8 ), one 0d fracture intersection point (Ω 4 ), and three 0d fracture end-points (Ω 1 , Ω 2 , Ω 3 ).Note that we allow fractures and other lower-dimensional subdomains to form parts of the boundary of the domain (e.g., Ω 5 with its endpoints Ω 1 and Ω 2 ).Center: Interfaces between subdomains.Right: Subdomain boundaries.Internal boundaries are depicted in red, whereas fracture's boundaries touching the ambient boundary are depicted in green.can no longer be applied directly.Thus, in the remainder of the paper, we deal with these challenges in a more general framework.

Extension to fracture networks
In this section, we extend the single fracture model to account for several subdomains as part of a general fracture network.Our vocabulary is motivated by the physical case of n = 3, where the surrounding rock is composed of simply connected 3d subdomains, fractures are simply connected planar 2d subdomains, the intersection between such fractures are 1d lines, and the intersection between fracture intersections are 0d points (see Figure 3 for an example with n = 2).
We start with the classical description and then introduce the mD notation.The rest of the section is devoted to introducing key tools that are necessary to perform the analysis in an mD setting.

Mixed-dimensional geometric representation
The derivation of a posteriori estimates for generic fracture networks greatly benefits from an mD decomposition of the domain of interest, and we therefore follow the approach of [12].We start by considering an n-dimensional contractible domain Y ⊂ R n , n ∈ {2, 3}, decomposed into m planar, open and non-intersecting subdomains Ω i of different dimensionality i=1 Ω i (see left Figure 3).The partitioning is constrained such that any d-dimensional subdomain (for d < n) is always either the intersection of the closure of two or more subdomains of dimension d + 1, or a cut in a domain of dimension d + 1.This hierarchical structure excludes e.g., a 1d line or a 0d point embedded directly in a 3d domain.
We adopt a structure where neighboring subdomains one dimension apart are connected via interfaces, denoted by Γ j for j ∈ {1, . . ., M }.To be precise, let Γ j be the interface between subdomains indexed by  and  of dimension d and d + 1, respectively.Then Γ j = Ω  (see center Figure 3), and furthermore, we denote the adjacent boundary of Ω  by Γ j = ∂ j Ω .We emphasize that while the internal boundary ∂ j Ω  is defined to spatially coincide with the interface Γ j , which in turn coincides with the lower-dimensional subdomain Ω , their distinction is crucial to define variables properly.
To keep track of the connections from subdomains to interfaces, we introduce the sets Ŝi and Ši , containing the indices of the higher-dimensional (respectively lower-dimensional) neighboring interfaces of Ω i , as illustrated in the right panel of Figure 3.These sets are dual to  and  defined in the previous paragraph, thus for all j ∈ Ŝi , it holds that  = i, while for all j ∈ Ši , it holds that  = i.
We will be interested in defining functions on the above stated partition of the domain and the interfaces.This motivates us to define the disjoint unions A complete mixed-dimensional partitioning, including both subdomain and interfaces, is given by Ω ⊔ Γ.
In order to speak of boundary conditions, we introduce the decomposition of the boundary of Ω.Let ∂Ω be partitioned into its Neumann, Dirichlet, and internal parts.That is, we define Finally, to ensure the existence of a unique solution, we require ∂ D Ω = ∅.

The model problem for a fracture network
Let us now present the model problem valid for m subdomains of dimensionality 0 to n, and M interfaces of dimensionality 0 to n − 1.Our model summarizes the derivations given in recent literature [12,30,39].For all domains Ω i , we consider a scalar pressure p i together with a flux u i in the tangent space of the domain.On all interfaces Γ j , we consider a scalar coupling flux λ j , oriented as positive for flow from the higher dimensional domain Ω .We will, in this section, assume sufficient regularity that the strong form makes sense, and return to the weak formulation in later sections.The governing equations from the previous section then generalize as In (19a), the summation captures the contribution of fluxes from the adjacent interfaces to Ω i , and can be seen as a generalization of the second term in (2a).Note that for d i = n, the set Ŝi = ∅, and thus the jump operator, evaluates to zero in the highest-dimensional domains.Conversely, in (19a), the differential term ∇ i •u i is void whenever d i = 0, as there is no tangent space to a point in all subdomains, and indeed, we will not consider the u i defined on these domains, which justifies why equation (19b) are not applied to 0d domains.
We are now ready to recast the model problem in mD notation, building on the product space structures introduced in Section 2.2.Let us start by defining the mD pressure as the ordered collection of subdomain pressures p := [p i ] ∈ CΩ, i.e., scalar functions on Ω.We now decompose the fluxes as in (11), so that such that u 0,i satisfies u 0,i • n i = 0 for all j ∈ Ši , and where the reconstruction operator is generalized as R j : CΓ j → CΩ  satisfying: This allows us to define the mD flux as the internal (tangential) domain fluxes and (normal) interface fluxes u := [u 0,i , λ j ] ∈ C 0 T Ω × CΓ, i.e., the pairing of sections of the tangent bundle T Ω together with scalar functions on Γ.By the subscript C 0 T Ω, we indicate that both u i • n i = 0 on all ∂ j Ω i , where j ∈ Ši , and also We now define a generalized divergence operator D • (•) : C 0 T Ω × CΓ → CΩ which acts on the mD flux in accordance with the left-hand side of (19a): where q = [q i ] ∈ CΩ is a scalar function for each domain Ω i , defined by: Similarly, we define an mD gradient operator D (•) : CΩ → CT Ω×CΓ acting on the mD pressure in accordance with the right-hand sides of equations (19b) and (19c): where v = [v 0,i , ν j ] ∈ CT Ω × CΓ has the same structure as the mD flux (but without the boundary conditions), such that for all i ∈ {1, . . ., m} and j ∈ {1, . . ., M }, it holds that Recalling that the full flux v i is recovered from equation ( 20), we note that the second term above is simply the gradient on each subdomain.We will, in Section 3.3, further justify the terminology "divergence" and "gradient" due to the fact that these operators satisfy an integration-by-parts property with respect to the suitable inner products, and are thus adjoints (subject to appropriate boundary conditions).
Material parameters are collected into the mD permeability K : CT Ω×CΓ → CT Ω × CΓ, defined such that for then from the model given in equation ( 19), we recognize the desired relationships The second term, corresponding to Darcy's law, can be rewritten in terms of the decomposition u = [u 0,i , λ j ] from equation ( 20) as: The presence of the extra terms arising from the decomposition is analogous to that in (19).We note that the restriction u ∈ C 0 T Ω × CΓ, implicitly places constraints (depending on the material constants K and via the definition of D ) on the admissible pressures p.This space of admissible pressures can be understood as the domain of the restricted operator In view of the mD variables and operators defined above, and subject to the right-hand side data f = [f i ] ∈ CΩ and the boundary data g D = [g D,i ] ∈ C∂ D Ω, a straightforward substitution of definitions shows that problem (19) is equivalent to the concisely stated mD elliptic problem defined for u ∈ C 0 T Ω × CΓ and p ∈ CΩ.

Remark 5 (Internal Neumann boundaries).
For simplicity of exposition, the domain Y is taken as contractible, and Ω i is considered a partitioning of Y .However, the reader will appreciate that these assumptions can be relaxed.Most importantly, from the perspective of applications (as discussed in Section 2.1), some internal interfaces may be modeled as impermeable, i.e. λ j = 0. We refer to the remaining (permeable) interfaces as Ξ ⊂ {0, . . ., M }.The impermeable interfaces are then excluded from the problem, and considered as internal Neumann interfaces.To be precise, we define a reduced disjoint union of interface domains Γ = j∈Ξ Γ j .
The internal Neumann boundaries may partition the domain into disconnected parts.We refer to a subdomain as "Dirichlet-connected", denoted i ∈ ξ if either (1) there exists some j ∈ Ŝi such that  ∈ ξ, or (3) there exists some j ∈ Ši such that  ∈ ξ.This allows us to construct a reduced disjoint union of subdomains All the derivations in the continuation are equally valid for these reduced product domains.
Remark 6 (Extensions to the model equations).The results of this paper can with minor modifications be extended to non-zero Neumann boundary conditions, and with some additional effort to the class of non-planar geometries considered in [15].However, as this generality is typically not needed for applications, we restrict the presentation as indicated above.

Variational formulations in mixed-dimensional notation
Before writing the variational formulations in mD notation, let us first define the relevant mD inner products and norms.Consider the following inner-products and their respective induced norms With these inner products, the previously defined mD divergence satisfy the following integration-by-parts formula [12,15] whenever v ∈ CT Ω × CΓ and q ∈ CΩ.
In the above the restriction to the boundary is denoted T X (•) (for X = D, N ), which depending on context acts as the boundary values of pressure variables, From the product structure in the definition of the C and L 2 spaces, the continuous spaces inherit their density from the individual subdomains to the product spaces on Ω and Γ.We can thus follow standard procedures to obtain weak extensions of the mD differential operators, the boundary restriction (trace) operators, and the corresponding function spaces [40,41,42].We elaborate this below.
Due to the density of , the mD divergence from Section 3.2 is a densely defined unbounded linear operator on the latter space Let us now (temporarily) use the notation (T, dom(T )) to emphasize that an operator T has domain of definition dom(T ), and we denote the adjoint operator with respect to the L 2 inner product by an asterisk.
We recall that the Neumann boundary is incorporated into the definition of the continuous flux spaces C 0 T Ω × CΓ, thus the last term in the integrationby-parts formula (34), is zero.Hence, we can define a weak mD gradient and the corresponding space of weakly mD differentiable functions with zero trace on the Dirichlet boundary H 1 0 by considering the adjoint: Clearly, C 0 Ω ⊆ H 1 0 (Ω), and thus it is appropriate to consider (D , H 1 0 (Ω)) as a weak gradient.Moreover, the domain of definition simply corresponds to the standard H 1 0 (Ω i ) on each domain, where the subscript zero indicates zero trace on all Dirichlet boundaries.Thus H 1 0 (Ω) = m i=1 H 1 0 (Ω i ), which generalizes (5).Considering the integration-by-parts formula again, the weak mD divergence and the corresponding space of flux functions with divergence in L 2 and zero trace on the Neumann boundary H(div; Ω, Γ) can be defined as Again C 0 T Ω × CΓ ⊆ H(div; Ω, Γ), and it is appropriate to consider (D•, H(div; Ω, Γ)) as a weak divergence.This domain of definition of the weak divergence has the interpretation of H 0 (div; Ω i ) on all subdomains Ω i (where the subscript zero indicates zero trace on all boundaries except for Dirichlet boundaries), and L 2 (Γ j ) spaces on all interfaces Γ j .Thus H(div; , which generalizes (12).Due to the above identification of H 1 (Ω) and H(div; Ω, Γ) in terms of product spaces of standard function spaces on subdomains, we extend the definition of the boundary restriction operators T X (•) to trace operators on the weak spaces by requiring that they coincide with the standard trace operators on subdomains.
In the continuation, we will always consider the weak mD gradient and divergence, and denote these simply by D and D•, respectively.Similarly, we will always consider the boundary restrictions as trace operators.The above definitions of weak mD gradient and divergence operators, and their adjoint property on the above weak spaces, has the following statements of the primal and dual weak formulations of equations ( 29) as a direct consequence: Definition 4 (Mixed-dimensional dual weak form).
The above weak forms of the mixed-dimensional elliptic problem are wellposed for bounded coefficients [15], in the sense that there exist positive constants K 0 and K ∞ such that: The solutions of the primal and dual weak formulations are equivalent, and define true solutions p ∈ H 1 0 (Ω) + g and u ∈ H(div; Ω, Γ) against which the approximate solutions will be measured in later sections.

Functional analysis tools
In this section, we summarize the main functional analysis tools we will exploit for the a posteriori analysis.

Poincaré-type inequalities
We recall the following weighted Poincaré inequalities: Lemma 1 (Permeability-weighted Poincaré-Friedrichs inequalities).There exist Here, we denote by qΩi and qK the mean value of q over the subdomain Ω i and an arbitrary d i -simplex K ⊂ Ω i , respectively.
We refer to C Ω,Γ as the mixed-dimensional permeability-weighted Poincaré-Friedrichs constant (whose existence was shown in [15]), C Ωi is the standard subdomain permeability-weighted Poincaré-Friedrichs constant, and C K is a local permeability-weighted Poincaré-Friedrichs constant.
It is important to mention that concrete values of C Ωi are available only for a limited set of geometries, see e.g., [43,44,45].An upper bound exists for convex domains, and thus for a simplex K ⊂ Ω i we have [46,47] where c K is the lower bound on the permeability within K: The importance of this is understood if K is an element of a simplicial partition of Ω i , in which case C K scales with the mesh size h K = diam(K).This allows for super-convergent properties of residual estimators for some locally mass-conservative approximations [26,48,49].We analyze these cases with further details in Section 5.2 and Remark 15.

Conforming flux spaces
It is often possible to verify that an approximate solution v ∈ H(div; Ω, Γ) satisfies some stronger conservation property, that is to say, that there is some space This allows for the construction of stronger a posteriori estimates, and as such, we formalize this concept as a generalization of H(div; Ω, Γ) to "U -conforming flux spaces": Definition 5 (Conforming mD flux space).Let H(div; Ω, Γ; U ) ⊂ H(div; Ω, Γ) be a U -conforming flux space, in the sense of To exploit the conforming flux spaces, we must construct certain projected H 1 (Ω) spaces.Consider therefore U as some subspace of L 2 (Ω) and define U ⊥ to be its orthogonal complement: Moreover, let π U ⊥ be the L 2 -projection onto U ⊥ , such that for any r ∈ L 2 (Ω), π U ⊥ r ∈ U ⊥ satisfies the orthogonality property: Consider now the projected , and let the norm of W be defined as a weighted L 2 -norm with nonnegative weights µ ∈ L ∞ (Ω) which are defined within the class C W with unit Poincaré constants: Indeed, such classes exist in the literature of Poincaré inequalities for weighted norms, see e.g., [29,28].Note that a trivial member of C W is the inverse of the permeability-weighted mD Poincaré-Friedrichs constant µ(x) = C −1 Ω,Γ .As we will see in Sections 5.1 and 5.2, the concrete choice of the space U and the corresponding weights µ will directly impact the strength of the estimates.
Remark 7 (On the space H(div; Ω, Γ; U )).The conforming mD flux spaces allow us to obtain sharper estimates in Section 5.However, it is important to note that the standard case U = L 2 (Ω) is included in our definition, for which the orthogonal complement is void, and the projection π W = I; thus W = H 1 0 (Ω).This and other cases are elaborated in more detail in Sections 5.2.1 to 5.2.4.

Bilinear forms and energy norms
For the a posterior i analysis, we will need the next two mD bilinear forms and their induced energy norms which are related via We also define the full norm for a mixed-dimensional pair of primal and dual variables as Note that the last norm will depend on the eventual choice of µ −1 , which we emphasize must be from the class µ ∈ C W , as defined in the preceding section.

A posteriori error estimates
This section is devoted to obtaining the error bounds for our model problem.
First, we provide general abstract estimates, and later we focus on the evaluation of the different bounds.

General abstract estimates
Let us now present the general abstract bounds.We formalize the main results presented in Section 3 and extend the ones presented in Theorem 1 in the following theorem.
Theorem 2 (General abstract a posteriori error bounds).Let the error majorant be defined as where valid for all q ∈ H 1 0 (Ω) + g and v ∈ H(div; Ω, Γ; U ).Then, the following a posteriori error estimates hold.
where M ⊕ p is the upper bound of the error for the primal variable.(2) Let u ∈ H(div; Ω, Γ) be the solution to (38) and v ∈ H(div; Ω, Γ; U ) be arbitrary.Then where M ⊕ u is the upper bound of the error for the dual variable.(3) Let p ∈ H 1 0 (Ω)+g be the solution to (37) and u ∈ H(div; Ω, Γ) be the solution to (38), and let (q, v) ∈ (H 1 0 (Ω) + g) × H(div; Ω, Γ; U ) be arbitrary.Then, where M ⊖ p,u and M ⊕ p,u are the lower and upper bounds of the error for the primal-dual variable.
Proof.Due to the construction of mixed-dimensional product spaces and the adjoint property of the weak differential operators, the proof from the monodimensional case can (to a large extent) be applied directly [19].A notable deviation from the standard proofs is the use of conforming flux spaces, and the inclusion of the Poincaré-constants in the weights C W .The full proof is included for completeness in Appendix C.
Remark 8 (Non-conforming approximations).Referring again to the general setting of mD calculus, it has been shown that the differential operators form part of a cochain complex, and that an mD Helmholtz decomposition exists [15].Thus, by realizing the above constructions as Hilbert complexes, the above error bounds can be extended also to non-conforming approximations following, e.g., Theorem 4.7 of [21].However, as a main objective of our work is to obtain bounds based on conforming properties of the approximations, we will not pursue non-conforming approximations in this work.

Evaluation of the majorant
The aim of this section is to provide concrete forms of the majorant M(q, v, f, µ) from Theorem 2 depending upon the choices of the weights µ.For this purpose, consider once again the definition of the majorant The estimation of the first term η DF (q, v) is independent of the weights µ.Indeed, by applying (50), it is straightforward to see that The terms η DF ,K and η DF ⊥ ,K measure the diffusive flux error in the tangential and normal directions associated with the subdomain element K ∈ T Ωi and the mortar element K ∈ T Γj , respectively.
To complete the evaluation of the majorant, we are left with the estimation of η R (v, f, µ), which depends on the choices of µ.Recall that this term measures the mismatch in satisfying the conservation equation in each subdomain.To be precise, there are four main types of conforming fluxes; Standard L 2 -conforming, subdomain conservation, grid level (local) conservation, and point-wise.The quality of the residual balance can be verified explicitly before applying the a posteriori estimates, and thus is not considered an assumption in the theory.Below, we make precise the aforementioned cases.

No mass-conservation
Assume nothing is known about the approximation of the residual terms beyond the L 2 structure.We indicate this case by the abbreviation "NC", and set U NC = L 2 , and v ∈ H(div; Ω, Γ; U NC ) = H(div; Ω, Γ).Then U ⊥ NC = 0, which implies that π W = I, and W = H 1 0 (Ω).Then, a priori, we only know the global (mixed-dimensional) Poincaré (40a), i.e., we have no better weight than setting µ(x) = C −1 Ω,Γ for x ∈ Ω.Using (49) and the mD Poincaré inequality (40a), one obtains the following bound, which is the weakest bound available within the class of bounds considered in this paper: Here, η R,K;NC denotes the local residual error for non-conservative approximations.The majorant when mass conservation cannot be guaranteed at any level is then given by, and it follows from the above that this is an upper bound, M ≤ M NC .

Subdomain mass-conservation
Due to the structure of the equations, where interface fluxes are stated explicitly, many approximations will have mass conservation satisfied in a subdomain level, which is in a sense a compatibility condition on the floating domains Ω i .We indicate this case by the abbreviation "SC".In particular, the divergence r = Thus, by definition U ⊥ SC is the space of constants over the floating subdomains Ω i , and the space W is the space of H1 (Ω i ) functions, with zero mean if This case represents an improvement relative to the previous one, in the sense that we can now employ the subdomain Poincaré constants instead of the mD constant.Let us make this point precise in the following lemma.H1 (Ω i ), where Then, µ(x) = C −1 Ωi for x ∈ Ω i belongs to the class C W , where C Ωi is the permeability-weighted Poincaré-Friedrichs constants defined in Lemma 1.
Proof.Using the Poincaré inequality (40c) and the fact that the sum of broken norms is weaker than the full norm, the following result holds sup In view of Lemma 2, η R can be bounded as where η R,K;SC are the local residual estimators for subdomain mass-conservative approximations.The majorant for this case is given by This estimate is sharper than that in the preceding section, since C Ωi ≤ C Ω,Γ , thus whenever the assumptions of this section are satisfied, it holds that M ≤ M SC ≤ M NC .Note that ( 65) is identical in structure to the residual estimators (16e) and (16f) obtained in Theorem 1.However, they are fundamentally different in the sense that (65) do not require all subdomains to have a non-empty Dirichlet part but rather mass to be conserved in each subdomain Ω i .

Local mass-conservation
By choice of numerical method, it is often easy to verify that mass is conserved on an element basis in a subdomain partition.We indicate this case by the abbreviation "LC".As in the preceding section, this implies that the divergence r where T Ωi denotes a finite partition of Ω i (typically a simplicial grid).In this case, U LC contain functions having zero mean on each element K ∈ T Ωi , and from (66) we see that U ⊥ LC = m i=1 P 0 (T Ωi ).We will consider the slightly weaker case, where (66) is only required to hold for all "non-Dirichlet boundary" elements, that is for all elements where ∂K ∩ ∂ D Ω = ∅.This is sufficient for the results from Lemma 2 to be extendable to the grid level by considering the space W Ω = m i=1 K∈TΩ i H1 (K), where H1 (K) is defined in (63).Lemma 2 now applies without modification, and weights µ(x) ≥ C −1 K for x ∈ K are therefore in C W .Moreover, thanks to convexity of simplicial grid elements, the local permeability-weighted Poincaré-Friedrichs constants are now fully computable.This allows us to bound η R,Ω as follows: where η R,K;LC are the local residual estimators for locally mass-conservative approximations.Using the above results, the majorant for locally mass-conservative approximations reads The local residual estimates η R,Ω;LC correspond to the ones previously obtained by [26,49] for mono-dimensional problems subject to a flux equilibration step.Since C K ≤ C Ωi , then, as before, whenever the assumptions of this section are satisfied, it holds that M ≤ M LC ≤ M SC ≤ M NC .
Remark 9 (Fully computable residual estimators).Unlike estimators obtained with residual methods (containing unknown constants [50,51]) or a purely functional approach such as in Sections 5.2.1 and 5.2.2 (containing constants that are generally difficult to determine [20]), estimators such as (67) contain only known local constants depending on the mesh size and material parameters.This justifies the claim that these estimators are fully computable.

Exact mass-conservation
Methods with local mass conservation, as discussed in the previous section, when applied to problems where the RHS data f is zero or piecewise constant, can then often be verified to have an exact (pointwise) conservation property.We indicate this case by the abbreviation "EC", for which f = D • v, so that U EC = 0 and U ⊥ EC = L 2 (Ω).Now, π W = 0 and W = 0. Thus, any finite weights µ are admissible, yet the choice is immaterial since the residual term µ −1 (f − D • v) Ω always evaluates to zero.Consequently, only diffusive-type errors are present in the a posteriori estimation, and the majorant takes the form This case can also be seen as the limiting case of local mass conservation for a family of grid partitions where h K → 0.

Summary of majorants and subdomain errors
With the obtained majorants, we can define the corresponding upper bounds for the errors of the primal, dual, and primal-dual variables.
Definition 6.Let α = NC, SC, LC, EC, corresponding to the flux conformity spaces U α discussed in the preceding sections.Then, in view of the results from Theorem 2 and the majorants (61), ( 65), (68), and (69), the upper bounds for the error in the primal, dual, and primal-dual pair, for arbitrary approximations q ∈ H 1 0 (Ω) + g and v ∈ H(div; Ω, Γ; U α ), are while the lower bound for the error in the primal-dual pair is It is our interest not only to measure local errors, but also to distinguish between subdomain and interface errors.This motivates the definition of the following errors estimators.Definition 7 (Subdomain and interface error indicators).Let α = EC, LC, SC, NC.Then, we will denote by ε Ωi;α and ε Γj the subdomain and interface error indicators, defined by We emphasize that while the majorants provide guaranteed bounds, the subdomain and interface error indicators can only be expected to correlate with the error.
Figure 4: Left: Matching coupling between the grids T Ω  , T Γj , and T Ω .Right: Degrees of freedom involved in the coupling between a 2d higher-dimensional cell, a 1d mortar-cell, and a 1d lower-dimensional cell.Locally, tangential fluxes are approximated using RTN 0 (K), whereas mortar fluxes and pressures using P 0 (K).

Concrete bounds for locally mass-conservative approximations
In this section, we will make the evaluation of the bounds concrete by providing explicit approximations to (38) using the lowest-order mixed-finite element method (MFEM).

Grid partitions
Ultimately, a posteriori estimates are primarily applied to approximations that are defined on computational grids.We therefore, in this section, summarize the relevant notation for grids and the mapping operators between subdomains and interfaces.Let us start by defining the partitions of the domains of interest.To this aim, denote by T Ωi , T Γj , and T ∂iΩj the partitions of Ω i , Γ j , and ∂ j Ω i , respectively.Moreover, let T Ω = ∪ m i=1 T Ωi , T Γ = ∪ M j=1 T Γj , and T ∂I Ω = ∪ m i=1 ∪ j∈ Ši ∂ j Ω i represent the union of all subdomain, mortar, and internal boundary grids.
Here, we only consider simplicial partitions.In particular, we require all elements K ⊂ Ω i to be strictly non-overlapping simplices of dimension d K = d i .We use h K to denote the diameter of K, and define h Ωi = max hK T Ωi , h Γj = max hK T Γj , and h ∂j Ωi = max hK T ∂j Ωi .
We will not at this point place any conditions on the grid partitions, although several aspects of this will be advantageous from the perspective of computation.

Finite element spaces and the approximated problem
Let us introduce the finite element spaces necessary to write the approximated problem.We start by defining a local space for the approximated pressures, mortar fluxes, and tangential fluxes.They are given, respectively by where P 0 and RTN 0 denote the spaces of constants and lowest-order Raviart-Thomas(-Nédélec) spaces of vector functions [52,53].See also Figure 4 for the degrees of freedom involved in the generic coupling between a (higherdimensional) triangle, a mortar line segment, and a (lower-dimensional) line segment.
The composite space for the approximated mD pressure Q h ⊂ L 2 (Ω) and the approximated mD flux X h ⊂ H(div; Ω, Γ) are defined respectively by and (72) While not strictly necessary from a theoretical perspective, in the discrete setting, it is often useful to choose a finite-dimensional reconstruction operator based on the discrete spaces, and we allow for this through the notation R h,j : Λ h,j → H(div; Ω i ), which in practice is often further restricted to R h,j : Λ h,j → V h, .Such discrete reconstruction operators are natural for matching grids, and can also be constructed in the more general case of nonmatching grids, see e.g., [12,54,55].Here Π h : Λ h,j → Λh,j is the L 2 projection from the mortar grid on Γ j to the boundary simplicial partition of Ω .
We have now all the elements necessary to write the finite-dimensional approximation to the dual mixed problem (38).

Definition 8 (Approximated mD dual mixed formulation). Find
Due to the presence of the discrete reconstruction operator, this approximation is conforming whenever Λ h,j = Λh,j , i.e., for matching grids.For non-matching grids, the approximation is still convergent, subject to normal conditions on the mortar grids [12].
Remark 10 (Conservation properties).Whenever equation (73b) is satisfied exactly, then equation (66) holds, and we have local mass conservation for matching grids.Thus, the fluxes lie in the smaller space X h ∩H(div; Ω, Γ; Q ⊥ h,i ), and the results from section 5.2.3 apply.Furthermore, if f i ∈ Q h,i and R h,j : Λ h,j → V h, , then the projection of the source term, and hence the residual error, onto Q ⊥ h,i vanishes.Thus, the local conservation is verified to be pointwise, the fluxes lie in X h ∩ H(div; Ω, Γ; 0) and the results from Section 5.2.4 apply.
Remark 11 (Well-posedness and a priori estimates).The stability and a priori approximation properties of the finite-dimensional system given in (73) has been previously established [12].

Pressure reconstruction
Recall that Theorem 2 requires any approximation to the mD flux to be in H(div; Ω, Γ), whereas approximations to the mD pressure must lie in H 1 0 (Ω)+ g.By the condition that u h ∈ X h ⊂ H(div; Ω, Γ), the solution of equations ( 73) by definition satisfy the first condition.On the other hand, the approximated mD pressure p h is only in L 2 (Ω).We therefore need to enhance the regularity of the approximated pressure and thus obtain a reconstructed pressure.
Definition 9 (Reconstructed pressure).We will call reconstructed pressure ph to any function constructed from the mD pair Remark 12 (On potential reconstruction).Several techniques for obtaining ph are available in the literature.Arguably, the simplest option is to perform an average of the P 0 (K) pressures on local patches and from there construct local affine P 1 (K) functions [56].Other techniques aim at solving first a local Neumann problem to obtain a P 2 (K) post-processed pressure, and then apply interpolation techniques to get energy-conforming potentials [27,57,26,58,59].Any of these choices are compatible with the bounds derived herein.
Remark 14 (Other locally mass-conservative methods).In addition to the MFEM scheme of the lowest-order (RT0-P0), other flux-based numerical methods such as the Mixed Virtual Element Method (MVEM) [60,61] and Cell Centered Finite Volume Methods (CCFVM), including the Two-Point Flux Approximation (TPFA) and the Multi-Point Flux Approximation (MPFA) [62,63], can be analyzed with our framework provided that the fluxes are interpolated in X h and the pressures reconstructed as indicated above.For methods without an explicit flux representation, an additional flux reconstruction step may be needed.
Remark 15 (Superconvergence of the residual estimators).Due to Remark 10, the residual estimators η R,K,LC are superconvergent for lowest-order locally mass-conservative approximations.This property is guaranteed since: (1) local Poincaré constants decay as O(h K ) for simplicial elements and (2) the norm of the residual [64]; leading to an overall rate of O(h 2 K ).We validate the a posteriori bounds and assess their efficiency on a 1d/2d problem (Section 7.2) and a 2d/3d problem (Section 7.3), both with manufactured solutions.The geometric configuration for both problems is shown in Figure 5.Let us denote the fracture as Ω 1 , the matrix as Ω 2 , the left interface as Γ 1 , and the right interface as Γ 2 .Further, assume the existence of an exact, smooth pressure p 2 (x) in Ω 2 .Refer to Table 7 and Table 8 from the Appendix D for the analytical expressions of all variables of interest.

Efficiency indices
Efficiency indices are used to assess the performance of the approximations when exact solutions are available.They are defined as the ratio between the estimated and the exact errors.Here, we consider the following efficiency indices.
Definition 10 (Efficiency indices).Let α = NC, SC, LC, EC and let p ∈ H 1 0 (Ω)+ g and u ∈ H(div; Ω, Γ) be the solutions to (37) and (38), respectively.Then, in view of Theorem 2, the efficiency indices for the primal, dual, and primal-dual pair, for arbitrary approximations q ∈ H 1 0 (Ω) + g and v ∈ H(div; Ω, Γ; U α ), are Remark 16.Optimal efficiency indices (equal to 1) are obtained when the approximations match the exact solutions.Moreover, in general the efficiency indices satisfy the bounds: For the final term, we note that since η R;α ≤ M α , then for α = NC, SC the total efficiency index satisfies I p,u;α ≤ 3, while for local conservation I p,u;LC ≤ 2 + O(h 2 ) and finally for exact conservation I p,u;EC ≤ 2. To this aim, we consider four levels of successively refined combinations of mesh sizes, characterized by

Two-dimensional validation
The global Poincaré constant is obtained numerically by solving the associated eigenvalue problem (see e.g., [66]), giving a value of C Ω,Γ ≈ 0.2251.
Further inspection shows that efficiency indices lie within the expected bounds discussed in Remark 16.In particular, efficiency indices for the primal variable using local weights are very accurate, and only a ∼ 7% deviation with respect to the actual error (for the finest grid) is observed in the case of RT0-P0, MVEM-P0, and MPFA.For TPFA, the efficiency index is worse, as a consequence of the flux approximation being worse.Efficiency indices for the dual variable are in general larger than the ones obtained for the primal variable; this is to be expected for mixed-dual approximations with the relatively simple pressure reconstruction, where the approximated fluxes have relatively good accuracy as compared to the reconstructed pressures.Finally, efficiency indices for the primal-dual variable are less than 2 for all methods in consideration.
Considering now the local error indicators, shown in Table 2, we note that diffusive errors decrease linearly for the matrix, fracture, and interfaces.Likewise, residual errors for the matrix and fracture decrease linearly when the global Poincaré-Friedrichs constant is used.When the local Poincaré-Freidrich constants are used, the residual estimators for the matrix and the fracture decrease quadratically, which goes in agreement with the super-convergent properties discussed in Remark 15.

Three-dimensional validation
For our next numerical validation, we employ the 2d/3d configuration from the right Figure 5.We repeat the same analysis from the previous section, and investigate four refinement levels.The mixed-dimensional Poincaré constant for this configuration corresponds to a value of C Ω,Γ ≈ 0.1838.The results are shown in Table 3 and Table 4.As in the previous validation, we can see that the majorants capture the local and global convergence tendency of all numerical methods.Again, RT0-P0, MVEM-P0, and MPFA give quite similar results, whereas TPFA showcase larger errors.As expected, efficiency indices again lie within the stated bounds from Remark 16.

Numerical applications
In this section, we apply our estimators to numerical approximations of challenging problems solving the equations of incompressible flow in fractured porous media.Importantly, since source terms are zero in both applications, by applying matching grids the residual errors are zero, and we are in the setting of having an exact conservation property from the numerical approximation.From Remark 16, we then know that the efficiency index for the primal-dual error will be less than 2; even if the exact solution and error are both unknown.

Two-dimensional application
In this numerical experiment, we consider the benchmark case 3b from [13].
As shown in the left panel of Figure 1, the domain consists of ten (partially intersecting) fractures embedded in a unit square matrix.The exact fracture The results for M ⊕ u;EC are omitted since they are equal to M ⊕ p;EC .coordinates can be found in Appendix C of [13].Fractures 4 and 5 represent blocking fractures (K = 10 −4 and κ = 1) whereas the others represent conductive fractures (K = 10 4 and κ = 10 8 ).The matrix permeability is set to one.A linear pressure drop is imposed from left (p = 4) to right (p = 1), whereas no flux is prescribed at the top and bottom of the domain.
The benchmark establishes three refinement levels; coarse, intermediate, and fine, with approximately 1500, 4200, and 16000 two-dimensional cells.The structure of the local contributions to the majorant (confer e.g.equation ( 59)) are shown in Figure 6, based on the approximate solution obtained by the MPFA discretization.
In Table 5, we show the errors bounds for the three refinement levels.To avoid numbering domains and interfaces, we refer to the matrix error as ε Ω 2 ,EC , and group the fracture and interface errors by conductive and blocking.For example, ε Ω 1 ,C,EC refers to the sum of the errors of 1d conductive fractures.
An important observation is that the persistent reduction of the majorant M ⊕ p,u;EC , together with the known upper and lower bounds on the efficiency indexes established in Remark 16, provides a post factum verification of the convergence of all the numerical methods considered.
The error estimates suggest that the contribution to the overall error bounds are concentrated, primarily, on highly conductive interfaces (see the column corresponding to ε Γ 1 ,C ).On a more qualitative note, Figure 6 suggests that subdomain diffusive errors are concentrated at the fracture tips and fracture intersections, which is where singularities may typically be encountered [12].

Three-dimensional application
Our last numerical application is based on a modified version of the threedimensional benchmark 2.1 from [14].The domain consists of nine intersecting fractures embedded in a unit cube, as shown in the middle panel of Figure 1.This results in an intricate network with 106 subdomains and 270 interfaces of different dimensionality.
The original benchmark imposes an inlet flux (purple lower corner u = −1) and an outlet pressure (pink upper corner p = 1), and for the rest of the external boundaries null flux.Since we have only detailed our results for zero Neumann boundary conditions, we have replaced the inlet flux by a constant pressure condition (p = 1) and modified the value of the outlet pressure (p = 0).The benchmark assigns heterogeneous permeability to the matrix subdomain, whereas the fractures are assumed to be highly conductive.For the complete description of the benchmark, we refer to [14], and for an impression on how the contributions to the majorant are distributed, see Figure 7.Here we show the error estimates for the whole fracture network obtained with RT0-P0, where it becomes evident that the subdomain diffusive errors are concentrated at the inlet and outlet boundaries; refinement efforts should therefore focus on these regions.
As in Section 8.1, we collect the local errors of subdomains and interfaces of equal dimensionality.The results are summarized in Table 6.As in the previous cases, we have local and global convergence for all four numerical methods.Again, RT0-P0, MVEM-P0, and MPFA show very similar results, while TPFA show larger errors.
As in the 2d case discussed above, the persistent reduction of the majorant M ⊕ p,u;EC , again serves as a verification of the convergence of all four numerical methods.

Conclusion
In this paper, we obtained a posteriori error estimates for mixed-dimensional elliptic equations.Depending upon the level of accuracy at which residual balances can be approximated, we have derived four concrete versions of the majorant; i.e.: for no mass-conservative, subdomain mass-conservative, locally massconservative, and point-wise mass-conservative approximations.Furthermore, we have demonstrated both theoretically and numerically that sharper bounds can be obtained (for locally mass-conservative methods) using local Poincaré constants instead of the global ones.
Our bounds have been thoroughly tested with numerical approximations obtained with four locally mass-conservative methods of the lowest-order, namely: RT0-P0, MVEM-P0, MPFA, and TPFA.We performed a detailed efficiency analysis comparing the use of global and local Poincaré-Friedrichs constants in two and three dimensions.In both validations, our upper bounds reflected the optimal convergence rates of the numerical methods.In addition, we applied our bounds to two-and three-dimensional community benchmark problems exhibiting challenging fracture networks.Again, in both cases, the bounds reflected the limitations and the convergence rates of the methods satisfactory.
To the best of our knowledge, the bounds obtained here are the first of their kind to provide a practical tool to measure the error in numerical approximations to the equations modeling the incompressible, single-phase flow in generic fractured porous media.

A Derivation of variational formulations
Here, we present the derivations for the primal and dual variational formulations for the case of a single fracture immersed in a matrix.
A.1 Derivation of the primal weak form for a single fracture Substitute (1b) into (1a), multiply each term by q 2 ∈ H 1 0 (Ω 2 ), and integrate over Ω 2 .Similarly, substitute (1b), (3a), and (3b) into (2a), multiply each term by q 1 ∈ H 1 0 (Ω 1 ) and integrate over Ω 1 .Add the resulting equations to obtain Using integration by parts, the first term of (77) can be expressed as Here, we use the internal boundary conditions (1c) and (1d) and the definition of the mortar fluxes (3a) and (3b).Analogously, integration by parts allows us to write the second term of (77) as Note that the boundary terms vanish due to the choice of boundary conditions.Finally, we note that the third and fourth terms from (77) can be equivalently written as The proof is completed by substituting (78), (79), and (80) into (77) and grouping common terms.
A.2 Derivation of the dual weak form for a single fracture Let us start with (13a).Multiply respectively (1b) and (2b) by v 0,2 ∈ V 0,2 and v 0,1 ∈ V 0,1 , integrate over the subdomains Ω 2 and Ω 1 , use integration by parts By noticing that the first three terms of (85) add up to the right-hand side of ( 6), and adding the identity ν j , (p 1 − q 1 ) − tr ∂j Ω2 (p 2 − q 2 ) Γj = 0, valid for any v 0 ∈ V 0 and ν ∈ L 2 (Γ) to (85), we obtain − ν j + κ j q 1 − tr ∂j Ω2 q 2 , (p 1 − q 1 ) − tr ∂j Ω2 (p 2 − q 2 ) Γj Recognizing that since K 2 is symmetric positive definite, it can be expressed as , where K is also symmetric positive definite, and therefore self-adjoint.The square-root of the material coefficients can therefore be moved to the second argument of the three first inner products in (86).After applying the Cauchy-Schwarz inequality to each inner product of (86), one gets Applying the permeability-weighted Poincaré-Friedrichs inequality (40b) to the terms p 1 − q 1 Ω1 and p 2 − q 2 Ω2 , the proof of the theorem is completed.

C Proof of Theorem 2
Here, we present the proof of our main theorem, which deals with the general abstract estimates in a mixed-dimensional setting.
Proof.(1) The proof for the bounds for the mD primal variable follows the one presented in Appendix B, modulo its generalization to the mD setting and the use of weighted norms on the residual terms.Start by computing the difference between any q ∈ H 1 0 (Ω) + g and p ∈ H 1 0 (Ω) + g using ( 49), to get Here, we used ( 49), (55), and added the fact that D• and D are adjoints.By exploiting the orthogonality property (46) and then introducing the weights to the second and third terms, (87) can be equivalently written as: Finally, applying the Cauchy-Schwarz inequality to the first and second terms of (88), and then the norm definitions (49), (50), and (47), we arrive at the desired bound: (2) The proof for the bounds for the dual variable is given next.We remark that an alternative proof based on a generalized abstract estimate (see [26], Theorem 6.1) can be used to obtain equivalent upper bounds after its generalization to the mD setting.
We start by adding the square of the primal and dual error to obtain: Here, we used the norm definitions ( 49) and ( 50) together with the mD constitutive relationship (29a).
Using partial integration, mass conservation (29b), and the orthogonality property (46), the second term of (90) can be equivalently written as Using the Cauchy-Schwarz inequality twice and the definition of the weighted norms ( 47), (91) can be estimated as Substituting ( 92) into (90) and applying the Cauchy-Schwarz inequality to the first term, we arrive at, 2 Ω , from which we conclude that (56) indeed holds.
(3) To prove the upper bound for the primal-dual pair, we choose an arbitrary pair (q, v) ∈ (H 1 0 (Ω) + g) × H(div; Ω, Γ; U ), and measure its difference with the exact solution (p, u) ∈ (H 1 0 (Ω) + g) × H(div; Ω, Γ) in the norm (52), to get where we use the bounds (55) and (56).For the proof of the lower bound, we start from the definition of the majorant, to get This completes the proof for the two-sided bounds and the abstract theorem.

D Exact solutions to numerical validations
Herein, we provide the exact expressions for the pressure, velocities, mortar fluxes, and source terms for the numerical validations presented in Section 7.
We will conveniently define the following quantities for notational compactness:

Figure 1 :
Figure 1: Example geometries falling within the context of hierarchical mixeddimensional geometries studied herein.Left figure corresponds to a 2d benchmark problem[13] while the two remaining correspond to 3d benchmark problems[14].

1 Ω 2 Figure 5 :
Figure 5: Geometric setups used for the numerical validations.Left: A 1d fracture embedded in a 2d matrix and the exact pressure solution.Right: A 2d fracture embedded in a 3d matrix.

For our first
validation, we consider the 1d/2d case as shown in the left Figure 5.This validation has two purposes: (1) compare the majorants and efficiency indices obtained using global (no mass-conservation) and local (local mass-conservation) Poincaré-Friedrichs constants, and (2) show the different errors associated with subdomains and interfaces.

Figure 6 :
Figure 6: Two-dimensional benchmark problem and the errors associated with the matrix and fractures for the coarse (left), intermediate (center), and fine (right) grid resolutions.Fractures 4 and 5 are blocking, whereas the others are conductive.The local bounds were obtained using MPFA.The results suggest that subdomain diffusive errors are concentrated around fracture tips and fracture intersections.

Figure 7 :
Figure 7: Subdomain diffusive error contributions to the majorant for the fine grid resolution obtained with RT0-P0.

Funding:
Jhabriel Varela was funded by VISTA -a basic research program in collaboration between The Norwegian Academy of Science and Letters, and Equinor.Additionally, this work was supported in part through the Norwegian Research Council grant 250223.The authors would like to thank W. M. Boon for helpful discussions on this topic.
The results for M ⊕ u;NC and M ⊕ u;LC are omitted since they are equal to M ⊕ p;NC and M ⊕ p;LC .

Table 2 :
Two-dimensional validation: Subdomain and interface errors.

Table 4 :
Three-dimensional validation: Subdomain and interface errors.

Table 5 :
Error estimates for the two-dimensional application.

Table 6 :
Error estimates for the three-dimensional application.
The results for M ⊕ u;EC are omitted since they are equal to M ⊕ p;EC .