Robust Discretization of Flow in Fractured Porous Media

1. Introduction.

• The design goals are mass conservation, grid flexibility, robustness for arbitrarily small apertures, and provable stability and convergence.Conservative flow matters for transport sensitivity; flexible grids can accommodate fracture networks and other features such as heterogeneities and wells.Source: §1 · p. 1
• No reviewed method is identified as meeting all four design goals; XFEM handles independent fracture meshing but cut elements can become arbitrarily small, requiring special stability constructions.The proposed alternative allows non-matching grids, though mesh construction may be more involved in complicated cases.Source: §1 · p. 2
• The paper uses mixed finite elements because flux is an explicit variable, enabling local mass conservation and rigorous analysis.Source: §1 · p. 2
• The method uses fluxes as mortar variables in a mixed finite element framework, coupling fractures and the surrounding domain with a stronger notion of mass conservation.The paper notes that earlier mortar MFE methods using pressure in the mortar space do not allow strong flux continuity.Source: §1 · p. 2
• The hierarchical formulation unifies treatment of the permeable domain, fractures, intersection lines, and intersection points, and is designed for arbitrary dimensions.The method is analyzed as robust with respect to aspect ratio, although the analysis excludes degenerate normal permeability.Source: §1 · p. 2
• The paper presents analysis and numerical experiments in two and three dimensions; a full numerical comparison with alternative discretizations is reported separately in benchmark study [14].Source: §1 · p. 3

2.1. Geometric Representation.

• The domain is partitioned by dimension: fractures are codimension-one subdomains, while intersections are represented as lower-dimensional features.Intersections of d-dimensional subdomains are treated as (d−1)-dimensional domains and excluded from the d-dimensional domains.Source: §2.1 · p. 3, Figure 1 · p. 4
• The framework covers ambient dimensions n=2 or 3, but this work excludes subdomains of codimension two.Domains of any dimension may contact the model boundary; codimension-two cases can be accommodated by adding intermediate-dimensional subdomains.Source: §2.1 · pp. 3–4

2.2. Governing equations.

• Mass conservation and Darcy’s law govern flow in the surrounding porous material and in all lower-dimensional subdomains.Applying Darcy’s law in fractures represents either thin open fractures with Poiseuille flow or fractures filled with material.Source: §2.2 · pp. 4–5
• Spatially varying geometric scaling parameters ϵ and γ are allowed to decrease toward zero, representing closing fractures.The model assumes these parameters are bounded and known a priori.Source: §2.2 · p. 5
• Lower-dimensional conservation combines tangential flux divergence with jumps in normal flux entering across interfaces.The scaled conservation law is ∇·(ϵu^d) + Jˆϵλ^dK = ϵ^2f^d for 1 ≤ d ≤ n−1.Source: §2.2 · pp. 6–7
• Fractures may meet the domain boundary; immersed tips have no-flux conditions, and the analysis requires nondegenerate normal permeability.Where a fracture ends on a Neumann boundary, no flux is imposed; at a Dirichlet boundary, pressure is set to the cross-sectional average of g.Source: §2.2 · pp. 7–8

2.3. Weak Formulation.

• Normal fluxes across interfaces form a mortar space in the mixed weak formulation, coupling tangential fluxes and pressure across domains.The formulation introduces an extension operator to impose the interface flux condition and expresses the resulting system as a saddle-point problem.Source: §2.3 · pp. 8–10
• When fracture aperture and scaling vanish uniformly, the model reduces to domain decomposition with a flux mortar enforcing weak pressure continuity.Source: §2.3 · p. 10
• The weak formulation couples domain and interface variables through bilinear forms, with forcing and prescribed boundary data on the right-hand side.Source: §2.3 · p. 10

2.4. Well-posedness.

• The extension operator has a stated bound controlling its domain-field and scaled auxiliary-pressure terms by interface data.Source: §2.4 · p. 11
• The extension-based inner product makes the space X a Hilbert space, providing completeness for the well-posedness analysis.Source: §2.4 · pp. 11–12
• The analysis notes that avoiding the stated restriction on ν requires weighted Sobolev spaces, as in.Source: §2.4 · p. 12
• The analysis establishes continuity and ellipticity of the bilinear forms and an inf-sup condition for the pressure coupling.Source: §2.4 · pp. 12–13
• The constants in the preceding lemmas are independent of γ and ϵ; aperture dependence is reflected entirely in the norm definitions.Source: §2.4 · p. 14
• Theorem 2.5 establishes existence and uniqueness for the weak problem in the energy norm.The proof applies saddle-point theory after establishing completeness, continuity and ellipticity, and an inf-sup condition.Source: §2.4 · p. 14

3.1. Discrete Spaces.

• Discrete spaces are dimensionally composite and use stable mixed finite-element pairs, with the mortar mesh chosen sufficiently coarse relative to higher-dimensional traces.The coarse-mortar condition helps eliminate a nonzero projection kernel that could cause mortar-flux oscillations, particularly because γ is typically small.Source: §3.1 · pp. 15–16
• Normal fluxes serve as mortar variables, strengthening mass conservation in the coupling between fractures and the surrounding domain.The finite-element problem solves for the surrounding-domain flux, mortar flux, and pressure.Source: §3.1 · p. 16
• The discrete extension operator enforces the projected mortar flux on the fracture boundary and can be supported near that boundary to yield a sparse pattern.For the analysis, a specific extension operator is constructed, while practical choices may instead accommodate the essential boundary conditions for the combined flux.Source: §3.1 · p. 16, §3.2 · p. 18

3.2. Stability and Convergence.

• The discrete problem has a unique solution and a stability estimate; its pressure inf-sup constant is independent of γ, ϵ, and h.These results hold when the mesh and function spaces satisfy the restrictions in Subsection 3.1.Source: §3.2 · pp. 18, 20
• The stability estimate implies existence and uniqueness of the solution because the discrete system is square.The estimate follows after choosing δ1 sufficiently small.Source: §3.2 · p. 21
• The combined flux error is bounded by projection errors and converges at order k + 1 on quasi-uniform grids under the stated regularity assumptions.The estimate scales as h^(k+1) times the solution norms for u, its divergence, λ, and p.Source: §3.2 · p. 24
• The numerical studies introduce a two-dimensional fracture with a tip that gradually decreases to zero and a three-dimensional problem illustrating full-dimensional decomposition.The two-dimensional example is designed to indicate that condition (2.15) may not be necessary.Source: §4 · p. 24

4.1. Two-Dimensional Problem.

• The two-dimensional test includes an intersection, immersed and boundary-crossing fractures, two models of fracture endings, and independently meshed fracture sides.A virtual zero-aperture extension eliminates tangential flow and reduces the method to a mortar scheme with normal flux as the mortar variable.Source: §4.1 · p. 25, Figure 3 · p. 25

4.1.1. Qualitative Results.

• The two-dimensional solutions are oscillation-free and physically consistent despite fracture endings, an intersection, and non-matching grids.High permeability yields nearly continuous pressure, while low permeability creates pressure discontinuities; abrupt endings produce less regular solutions than gradual pinch-outs.Source: §4.1.1 · p. 26
• Violating the aperture-gradient condition produces minor oscillations near fracture tips.The test supports the condition as necessary for the method as well as for the analysis.Source: §4.1.1 · p. 26

4.1.2. Convergence.

• Average convergence is linear in all variables in the two-dimensional test, while point-evaluation rates near the intersection are grid-dependent and erratic.The flux comparison excludes small balls around fracture tips, where flux behavior is irregular.Source: §4.1.2 · pp. 26–27

4.2. Three-Dimensional Problem.

• In three dimensions, the computed variables show at least first-order convergence in each dimension, and flow is predominantly top-to-bottom through conductive fractures.The test domain is decomposed into subdomains, fracture planes, intersection lines, and an intersection point.Source: §4.2 · pp. 27–28, Table 3 · p. 28

4. Numerical Results.

• The proposed mixed finite element scheme uses flux mortars to conserve mass while handling non-matching grids.Its hierarchical domain subdivision supports the discretization.Source: §5 · p. 28

• Hierarchical domain subdivision and dimensionally composite function spaces underpin the method's stability and a priori error analysis.Source: §5 · p. 29
• The analysis shows robustness and convergence for varying, arbitrarily small apertures.Source: §5 · p. 29
• Numerical results suggest the analysis's constraint on degenerate normal permeability may not be needed in practice; the three-dimensional solution is qualitatively consistent with conducting-fracture flow.Source: §5 · p. 29, Figure 5 · p. 29