Discrete Fracture Models in Porous Media

Updated 11 December 2025

Discrete Fracture Models are mathematical and computational techniques that simulate flow and transport in fractured porous media by representing fractures as lower-dimensional entities.
They employ various discretization strategies, including finite volume, hybrid-dimensional finite elements, and embedded models, to accurately capture matrix-fracture coupling and complex geometries.
Numerical benchmarks validate DFMs with second-order pressure convergence and extensions that address nonlinear flows, coupled geomechanics, upscaling, and machine learning surrogates for practical real-world applications.

A discrete fracture model (DFM) is a class of mathematical and computational techniques for simulating flow (and, in some cases, transport, deformation, or chemical processes) in fractured porous media, wherein fractures are explicitly represented as geometric entities of codimension one (e.g., surfaces in 3D or lines in 2D) embedded within the larger porous matrix. This hybrid-dimensional paradigm treats flow within the fractures and the matrix separately, with domain-specific governing equations and coupling laws, supporting accurate simulation of high-contrast heterogeneity and complex topologies characteristic of geological fracture networks. DFM approaches, including classical finite volume and finite element formulations, dual virtual element methods, mimetic finite difference, and embedded or projection-based models, form the computational backbone of modern studies on subsurface resource management, contaminant transport, and geomechanics in fractured media.

1. Mathematical Formulation and Limit Models

DFM theory systematically reduces the full-dimensional flow and transport equations to lower-dimensional models as the thickness (aperture) of fractures approaches zero, while scaling their permeability accordingly. Given Ω ⊂ ℝⁿ, a background matrix domain with a collection of (n–1)-dimensional smooth submanifolds Γ representing fractures, one begins with Darcy’s law and mass conservation in both Ω\f (matrix) and an ε-thickened neighborhood Ω_f^ε (fracture zone):

Matrix: uₘ = –Kₘ∇pₘ, –∇·uₘ = qₘ in Ω\Γ
Fracture: u_f = –K_f∇ₜp_f, –∇ₜ·u_f = q_f on Γ

The limit ε → 0, with fracture conductivity K_f scaling as ε^α, yields a rigorous classification of hybrid-dimensional limit models parameterized by α:

For α = –1 (critical scaling), the classical "lower-dimensional Darcy" limit is obtained, where the fracture behaves as a codimension-one manifold supporting a Darcy law for tangential flow, coupled to the matrix by continuity of pressure and jump in normal flux.
For α = 1 ("permeable barrier"), the interface reduces to a Robin-type jump condition relating discontinuity in pressure to flux across the barrier; this is fundamental for modeling flow-blocking features such as mineralized or cemented fractures.
Regimes for α < –1, –1 < α < 1, and α > 1 correspond to highly-conductive, vanishingly-thin, and impermeable limits, respectively, with rigorous H¹/L² convergence guaranteed under ellipticity and smoothness assumptions (Hörl et al., 2023).

These asymptotic regimes underpin all DFM discretizations and inform their coupling strategies for matrix and fracture variables.

2. Core Discretization Strategies

DFMs are implemented via a variety of spatial discretizations, distinguished by the treatment of fracture geometry, mesh conformity, and degree-of-freedom (DoF) allocation:

a) Finite Volume DFMs:

Classical box method ("Box-DFM"): The domain is discretized via a primary (primal) mesh (e.g., Delaunay triangles/tetrahedra), with a secondary (dual) "box" mesh centered at mesh vertices. Matrix flow is approximated with vertex pressures and flux balances over dual volumes. Fractures aligned with mesh edges contribute additional codimension-one flux terms. Matrix–fracture coupling is enforced either by superposing flux contributions along fracture intersections or by modifying transmissibility matrices. The inclusion of low-permeable (barrier) interfaces is achieved with minimal DoF overhead by splitting vertex DoFs at barrier crossings and imposing jump conditions (Xu et al., 29 Apr 2024).

b) Hybrid-Dimensional Finite Elements and Virtual Elements:

Mixed-dimensional FE/VEM approaches embed fractures as Dirac-delta measures in the global permeability tensor, allowing for non-conforming meshes and natural incorporation of curved or spatially varying fractures/barriers. The variational formulation is posed in H¹(Ω) and/or H¹(Γ), with additional interface integrals derived from the δ-support of fracture manifolds (Xu et al., 2021). VEMs provide robust local mass conservation and accommodate complex (polyhedral) cell shapes, advantageous for realistic networks and after mesh coarsening (1711.01818, Fumagalli et al., 2016).

c) Embedded and Projection-Based DFMs:

Embedded (EDFM/pEDFM/LEDFM): Fractures are embedded into an existing coarse or structured background matrix grid, decoupling mesh generation from fracture geometry. Matrix–fracture connections are assembled via geometric overlap and local upscaling of fine-scale subproblems. pEDFM further employs projection operators to robustly represent both high-conductivity and barrier-type fractures, particularly in complex 3D corner-point grids. LEDFM increases accuracy near barriers and tips by solving local fine-scale problems to obtain transmissibilities (HosseiniMehr et al., 2021, Losapio et al., 2022).

d) Mimetic Finite Difference (MFD) and Domain-Decomposition:

MFD approaches on conforming multidimensional meshes enforce strong local conservation, symmetry, and exact matrix–fracture coupling by collocating face variables or by Schur complement elimination of small intersection cells (Hyman et al., 2021, Stefansson et al., 2017).
Stabilized hybridized formulations permit individually meshed fractures and traces, using Lagrange multipliers and natural-norm stabilization to ensure stability and convergence without requiring inf-sup compatibility among discretized spaces (Berrone et al., 26 Jul 2024).

3. Barrier Representation and Pressure Discontinuity

Blockage by mineralized or otherwise low-permeability fractures is central to many DFM applications. DFM extension for barriers involves:

Replacing classic interface conditions of pressure continuity by explicit jump conditions relating pressure difference [p⁺ – p^–] to the normal flux, with transmissibility scaled by barrier thickness and permeability (e.g., T_b = k_b |e|/a in Box-DFM) (Xu et al., 29 Apr 2024).
In hybrid-dimensional formulations, barriers appear in the resistance tensor as Dirac-distributed terms, yielding additional interface integrals in the weak formulation (Xu et al., 2021).
For two-phase and non-linear flow, models allow pressure discontinuity and upwind transmission of phase flux, incorporating an explicit damaged rock layer of tunable thickness and porosity for enhanced numerical stability and physical realism (Droniou et al., 2016).

High-fidelity representations of both conducting and blocking fractures, as well as the interplay between tangential and normal modes of transport, are enabled by these mechanisms.

4. Numerical Benchmarks, Accuracy, and Computational Performance

Comprehensive validation on standard benchmarks and realistic synthetic domains confirms the high-fidelity and efficiency of state-of-the-art DFMs:

Second-order convergence in pressure for single-phase flows (Box-DFM, MFD, FE/VEM) with proper handling of barriers; robust representation of sharp pressure jumps and accurate flux exchange at matrix–fracture interfaces (Xu et al., 29 Apr 2024, Xu et al., 2021, Stefansson et al., 2017).
Sub-optimal but predictable rates on non-matching meshes, with error rates matching theoretical predictions (e.g., h^{1/2-ε} for non-matching stabilized three-field formulations), and robust conditioning when employing proper stabilization weights (Berrone et al., 26 Jul 2024).
Minimal DoF overhead due to barrier representation (e.g., split vertices only at crossings), retention of sparse symmetric positive-definite matrix structure, and applicability of standard preconditioned iterative solvers or direct Cholesky/AMG approaches.
Hierarchical domain-decomposition frameworks and Schur complement techniques efficiently eliminate small or lower-dimensional intersection cells, maintaining stability and accelerating convergence (Stefansson et al., 2017).

Benchmarks include canonical domains with vertical or slanted barriers, regular or complex networks mixing high-permeability fractures and low-permeability barriers, outcrop-inspired networks, and large 3D stochastically-generated testbeds, supporting general conclusions about method scalability, accuracy, and physical fidelity (Hyman et al., 2021, Xu et al., 29 Apr 2024).

5. Extensions: Nonlinear Flow, Coupled Problems, Upscaling, and Surrogates

DFMs now support a range of advanced physical and computational extensions:

Nonlinear and multiphase flow: Models accommodating general two-phase Darcy flows with pressure jumps and accumulation in thin damaged layers at interfaces; rigorous convergence for broad classes of spatial schemes via gradient discretization (Droniou et al., 2016).
Coupled geomechanics and THM: Stratified discretizations enabling strong two-way coupling between flow and fracture deformation, using either fully lower-dimensional DFM formulations with FV/MPSA discretization (Stefansson et al., 2020) or hybrid Embedded Discrete Fracture Models (EDFM) with strong discontinuity enrichment for mechanics (Shovkun et al., 2020).
Upscaling and reduced-order modeling: NLMC and multiscale basis construction techniques achieve order-of-magnitude DoF reduction with <1% error, by constructing nonlocal coarse operators through local saddle-point systems, supporting massive scale simulations (Vasilyeva et al., 2018).
Machine learning surrogates: Pre-trained deep convolutional neural network surrogates predict upscaled hydraulic conductivity tensors directly from local DFM realizations, achieving accuracy comparable to numerical homogenization while providing significant computational speedup (up to 28×) for use in MLMC workflows (Špetlík et al., 9 Jan 2024).

These extensions facilitate efficient simulation of highly heterogeneous, realistic systems with millions of fractures, while capturing complex coupled mechanisms relevant to subsurface resource management, geological storage, and contaminant migration.

6. Practical Recommendations, Limitations, and Future Directions

Best practices and guidelines extracted from contemporary literature include:

Use minimal-DoF extensions (e.g., barrier splitting in Box-DFM) for domains dominated by low-permeability features when tangential flow is negligible; otherwise, employ extended DFMs that support mixed-mode interface conditions (Xu et al., 29 Apr 2024).
Prefer hybrid-dimensional or embedded models for complex geometry, 3D domains, or when non-conforming mesh generation is essential (Hyman et al., 2021, HosseiniMehr et al., 2021).
For scenarios with significant tangential flow within barriers or detailed multiphase processes, select DFM variants incorporating full interface DoFs or explicit modeling of tangential transport.
Stabilized hybridized and dual VEM approaches are advantageous when fracture networks exhibit intricate intersection topology or when computational mesh flexibility is paramount (Berrone et al., 26 Jul 2024, 1711.01818).
Limitations remain in perfectly impermeable interfaces (singular limit for resistance terms), degenerate cases at coarse mesh resolution (e.g., pressure continuity vs. jump at intersections), and unresolved physical processes in damaged layers for multiphase flow.
Future work points toward improved treatments of tangential barrier transport, multiphase and non-Newtonian flows, adaptive hp-refinement, parallel-in-time and domain-decomposition preconditioners, as well as integration of statistical and machine-learning-based upscaling within real-time workflows.

DFMs thus constitute a mature, versatile framework rigorously grounding the simulation of fractured porous media under a wide spectrum of physical contexts, with continuing innovation in numerical analysis, algorithmic implementation, and high-performance computing (Xu et al., 29 Apr 2024, Xu et al., 2021, Hyman et al., 2021, Losapio et al., 2022).