Hierarchical Mixed-Dimensional Elliptic Equations

• Hierarchical mixed-dimensional elliptic equations are models coupling PDEs across submanifolds of varying dimensions via hierarchical and interface conditions.
• They employ variational formulations, exterior calculus, and stable discretization strategies to deliver robust numerical results in complex domains.
• Advanced functional a posteriori error estimators provide fully computable bounds for both potentials and fluxes, guiding efficient adaptive mesh refinement.

Hierarchical mixed-dimensional elliptic equations model processes in domains composed of substructures of various spatial dimensions interacting through coupling conditions. These systems emerge naturally in applications such as flow in fractured porous media, composite and network-embedded materials, and domains with boundaries or interfaces of variable codimension. Their theory encompasses the extension of exterior calculus, variational principles, functional spaces, and error estimation frameworks to collections of embedded manifolds (of dimensions 0 to $n$) interconnected in a hierarchical, directed-acyclic structure.

1. Geometric and Analytical Framework

A hierarchical mixed-dimensional geometry consists of an open, contractible domain $Y \subset \mathbb R^n$ partitioned into disjoint connected submanifolds ${\Omega_i^d}$ of dimension $d=0,1,\dots,n$ and interfaces ${\Gamma_j}$ of dimension $d_j = d_i - 1$. These components are arranged in a "forest" structure where every lower-dimensional manifold lies in the boundary of one or more higher-dimensional parent manifolds. The resulting structure is a directed acyclic graph, with recursive boundary maps and smooth coordinate charts linking manifolds.

Mixed-dimensional objects—differential forms, functions, or vector fields—are defined as tuples indexed over the full collection of $d$-dimensional manifolds. Functional spaces are constructed as broken Sobolev products (e.g., $H^1_0(\Omega)=\prod_{i\in I} H^1_0(\Omega_i)$), capable of capturing both continuity within each submanifold and interface conditions between them. The energy and inner product structures aggregate contributions from all manifolds and interfaces, yielding norms and inner products such as

$$(\mathfrak{a},\mathfrak{b})_{\mathfrak{F}} = \sum_{j\in \mathfrak{F}} \int_{\Omega_{s_j}} a_j \wedge \star b_j,$$

where $\mathfrak{a} = [a_j]_{j\in\mathfrak{F}}$ are mixed-dimensional $k$-forms (Boon et al., 2017).

2. Mixed-Dimensional Elliptic Model Formulations

The canonical mixed-dimensional elliptic system (e.g., for Darcy flow) is formulated by seeking pressure fields $p_i$ and fluxes $u_i$ in each subdomain $\Omega_i$, and interface fluxes $\lambda_j$ on each $\Gamma_j$, subject to: $$\begin{aligned} &\nabla_i\!\cdot u_i - \sum_{j\in \check S_i}\lambda_j = f_i & \text{in } \Omega_i, \ &u_i = -K_i \nabla_i p_i & \text{in } \Omega_i, \ &\lambda_j = -\kappa_j(p_- - p_+) & \text{on } \Gamma_j, \ &u_{\hat{\jmath}}\cdot n_{\hat{\jmath}} = \lambda_j & \text{on } \partial_j\Omega_{\hat{\jmath}}, \ &u_i \cdot n_i = 0 \quad (p_i = g_{D,i}) & \text{on } \partial_N\Omega_i\,(\partial_D\Omega_i). \end{aligned}$$ The primal weak form is posed in $H^1_0(\Omega)$ for pressures; the mixed (saddle-point) form is set in $(u,p) \in \mathbf X \times H^1_0(\Omega)$, with $\mathbf X$ a mixed-divergence-constrained space (Varela et al., 9 Dec 2025, Varela et al., 2021).

Exterior calculus for such geometries introduces mixed-dimensional differential operators:

• The semi-discrete exterior derivative $\mathfrak{d}$ combines tangential derivatives on each manifold with jump operators across interfaces.
• The coderivative $\mathfrak{d}^*$ acts as the adjoint under the $L^2$-inner product, aggregating codifferential and trace terms. The mixed-dimensional Hodge–Laplacian is then

$$\Delta_{MD} \mathfrak{a} = \mathfrak{d} \mathfrak{d}^* \mathfrak{a} + \mathfrak{d}^* \mathfrak{d} \mathfrak{a},$$

with associated variational principles and well-posedness established via Poincaré–Friedrichs inequalities adapted to the hierarchical structure (Boon et al., 2017).

3. Discretization Strategies: Matching and Non-Matching Grids

Each subdomain $\Omega_i$ and interface $\Gamma_j$ is discretized using independent conforming simplicial meshes, which may be unaligned ("non-matching") on their coupling boundaries. For non-matching grids, transfer grids are defined as the common refinement

$$\mathcal T_{\Gamma_j,X} = \{ K_\Gamma \cap K_X : K_\Gamma \in \mathcal T_{\Gamma_j}, K_X \in \mathcal T_X \},$$

where $X$ is a higher- or lower-dimensional adjacent mesh. Stable discrete projection operators (Scott–Zhang quasi-interpolants and $L^2$ projections) are constructed to transfer both primal potentials and dual (flux) variables between mismatched meshes. These operators, proven to be $L^2$- or $H^{1/2}$-stable, ensure numerical stability and consistency in saddle-point discretizations (Varela et al., 9 Dec 2025).

Hierarchical discretizations permit both mortar-type and conforming grid designs, and adaptation to problem-specific hierarchical embeddings, such as those encountered in multiscale fractured media (Varela et al., 2021).

4. Functional A Posteriori Error Estimation

Hierarchical mixed-dimensional elliptic equations admit functional-type a posteriori error estimators that provide fully computable, locally efficient, and guaranteed two-sided bounds for both the primal (potential) and dual (flux) variables. The general majorant is given by

$$\mathcal M^\oplus(q,v,f) = \biggl( \sum_{i,K} \eta_{\rm DF,\parallel,K}^2 + \sum_{j,K} \eta_{\rm DF,\perp,K}^2 + \sum_{i,K} \eta_{R,K}^2 \biggr)^{1/2},$$

where

$$\begin{aligned} \eta_{\rm DF,\parallel,K} &= \| K_i^{-1/2}(v_{0,i}+\sum_{j}\mathcal R_j^*\nu_j) + K_i^{1/2} \nabla_i q_i \|_K, \ \eta_{\rm DF,\perp,K} &= \| \kappa_j^{-1/2}\nu_j + \kappa_j^{1/2}(q_- - q_+) \|_K, \ \eta_{R,K} &= \frac{h_K}{\pi\,c_K} \| f_i - \nabla_i\cdot(v_{0,i}+\sum_{j}\mathcal R_j^*\nu_j) + \sum_{j}\nu_j \|_K. \end{aligned}$$

This estimator remains robust for both matching and non-matching grids, providing upper and lower efficiency indices closely tracking the true error even in the presence of complex coupling and local mesh misalignment (Varela et al., 9 Dec 2025, Varela et al., 2021).

Error estimation strategies can condition on the level of conservation enforced (none, subdomain, grid-level, or exact), with sharper bounds if local or exact conservation is achieved. All constants involved—such as local Poincaré constants—are computable from the mesh and geometric data (Varela et al., 2021).

5. Elliptic Theory in Domains With Mixed-Dimensional Boundaries

Mixed-dimensional elliptic operators are well developed for domains with boundaries or interfaces of varying dimension and even irregular (e.g., fractal) structure. Under geometric and measure-theoretic constraints (NTA topology, doubling measures), one constructs degenerate elliptic operators $L u = -\operatorname{div}(A(x)\nabla u)$ weighted by the volume or Hausdorff content of each stratum. Function spaces include weighted Sobolev spaces for the bulk and Besov-type spaces for boundary data. Trace and extension operators are defined, ensuring well-posedness of weak solutions up to complex boundaries (David et al., 2020).

Key results include boundary and interior regularity estimates (Caccioppoli, De Giorgi–Nash–Moser, Harnack inequalities), Dirichlet problem solvability, representation via harmonic measure and elliptic Green's functions, and robust comparison principles for solutions.

6. Applications and Numerical Validation

The theoretical and numerical frameworks are validated on community benchmark problems in fractured porous media, including 3D networks with non-matching fracture and matrix meshes. Numerical experiments demonstrate:

• The majorant and effectivity indices are insensitive (within 1–5%) to grid matching level.
• Localized error indicators correlate with regions of physical singularity (e.g., fracture tips, well neighborhoods).
• Mesh refinement reduces global error estimates monotonically, and spatial plots confirm physical intuition regarding error localization.
• For both matching and non-matching grids, the error estimator reliably guides adaptive refinement (Varela et al., 9 Dec 2025).

Tested discretizations include lowest-order mixed finite elements, mixed virtual elements, multi-point flux, and two-point flux approximations. The functional error framework applies independently of discretization, provided mass conservation properties are incorporated (Varela et al., 2021).

7. Theoretical Properties: Well-Posedness, Regularity, and Open Questions

The mixed-dimensional variational principle achieves strict convexity and coercivity under mild geometric assumptions and positive-definite coefficient data, leading to unique minimizers and well-posed weak formulations (Boon et al., 2017). The Poincaré–Friedrichs inequalities and Hodge decompositions generalize to the setting of hierarchical manifolds, allowing for stability and inf–sup theory analogous to classical PDEs.

Open research directions include the efficient numerical computation of mixed-dimensional Poincaré constants for arbitrary geometries, extension to non-linear or time-dependent PDEs, and generalization to higher (>3) spatial dimensions. Theoretical challenges remain in quantifying regularity and sharpness of error bounds in cases of minimal geometric regularity or highly singular coupling (Varela et al., 2021).

---

Key References:

"Functional Analysis and Exterior Calculus on Mixed-Dimensional Geometries" (Boon et al., 2017); "A posteriori error estimates for hierarchical mixed-dimensional elliptic equations" (Varela et al., 2021); "A posteriori error estimates for mixed-dimensional Darcy flow using non-matching grids" (Varela et al., 9 Dec 2025); "Elliptic theory in domains with boundaries of mixed dimension" (David et al., 2020).