Fractures such as joints and faults are prevalent geological features, forming hierarchical networks of discontinuities in crustal rocks (Ouillon et al., 1996; Bonnet et al., 2001; Lei and Wang, 2016). They play a critical role in various subsurface processes such as stress transfer, pressure diffusion, heat transport, and chemical dissolution (Tsang, 1991; Rutqvist and Stephansson, 2003; Viswanathan et al., 2022; Vaezi et al., 2026), which are highly relevant for many geotechnical engineering-related problems, such as civil infrastructure development, geothermal energy exploitation, critical mineral extraction, nuclear waste disposal, underground energy storage, and geohazard risk management. It is therefore essential to advance both our capability to model these ubiquitous discontinuity structures and our understanding of their influence on subsurface processes under complex geological and geotechnical conditions.

To address the need for modelling complex fracture systems, the discrete fracture network (DFN) approach has been developed since the early 1980s, with initial efforts focusing on fluid flow (Long et al., 1982; Long et al., 1985) and later extended to thermo-hydro-mechanical-chemical processes (Lei et al., 2017; Viswanathan et al., 2022). Unlike continuum approaches, the DFN method requires no a priori assumption of a representative elementary volume, making it particularly suited for modelling hierarchical fractured rocks that may lack a characteristic length scale (Ouillon et al., 1996; Bonnet et al., 2001; Lei and Wang, 2016). Over the past decades, the DFN method has evolved to become a cornerstone in modelling fractured geological media. In the following subsection, a concise overview of the DFN method is presented, outlining its fundamental concepts, key assumptions, and general workflow.

The DFN method treats fractures as distinct features by explicitly representing the geometry and behaviour of each individual fracture within a fracture system. Fractures are modelled as discrete, lower-dimensional objects—typically lines in 2D or discs/ellipses in 3D—embedded within the host rock (Dershowitz and Einstein, 1988). Each fracture is characterised by properties such as orientation, size, aperture, roughness, and mechanical/hydraulic attributes, which are often assigned specific values or sampled from probability distributions.

The workflow of DFN modelling typically involves four major steps: (1) statistical characterisation of fracture systems using field data; (2) stochastic or deterministic generation of fracture networks informed by these statistics; (3) mesh generation followed by process-based simulations; and (4) validation and calibration of DFN models against observational data, if available. In step (1), various fracture properties may be sampled from scanline surveys, outcrop mapping, and borehole logs, with their underlying statistical distributions inferred. In step (2), the DFN generation process can be either stochastic, where fracture properties are sampled from the inferred statistical distributions, or deterministic, where observed fractures are explicitly embedded within the model domain. In step (3), mesh generation creates a computational grid, enabling detailed process-based simulations of thermo-hydro-mechanical-chemical processes. In step (4), validation and calibration involve adjusting model parameters and comparing simulation results with observational data, such as hydraulic/tracer tests and microseismic/deformation measurements, to evaluate model accuracy and reduce uncertainties.

A key challenge of applying the DFN approach for practical engineering problems is realistically representing subsurface fracture systems, given the typically limited site characterisation data from outcrops, boreholes, and/or tunnel walls. Substantial efforts have been made to infer fracture statistics from geological mapping and use them to constrain fracture network generation (Andersson and Dverstorp, 1987; Dershowitz and Einstein, 1988; Kulatilake et al., 1993; Priest, 1993) as well as calibrate DFN models based on available field measurement data (Dverstorp and Andersson, 1989; Dverstorp et al., 1992). In addition, the computational cost of 3D simulations and the difficulty of upscaling fine-scale DFN results to field scales remain critical challenges for the broader application of the DFN method.

Over the past decades, extensive studies employing 2D DFNs have been performed to study different physical processes such as geomechanical deformation, seismic attenuation, fluid flow, solute transport, heat transfer, and chemical dissolution as well as their couplings in fractured geological media (Min et al., 2004; Lei et al., 2014; Lei et al., 2015; Lei et al., 2020; Lei et al., 2021; Lei and Gao, 2018; Lei and Gao, 2019; Kang et al., 2019; Hu and Rutqvist, 2020; Sun et al., 2020; Sun et al., 2021; Lei and Sornette, 2021b; Lei and Sornette, 2021a; Wang et al., 2021; Zhao et al., 2021; Zhao et al., 2022; Jiang et al., 2022; Jiang et al., 2024; Steefel and Hu, 2022; Cao et al., 2024). However, important 3D effects associated with fracture network configuration and multiphysical processes cannot be adequately captured in 2D models. Consequently, significant efforts in recent years have been devoted to the development and implementation of 3D DFN models to more accurately represent these complexities and gain deeper insights into the behaviour of fractured media in 3D space. While previous studies have provided comprehensive reviews of DFN methods and applications (Lei et al., 2017; Berre et al., 2019; Viswanathan et al., 2022; Vaezi et al., 2026), this Mini Review highlights recent advances with particular emphasis on 3D DFN modelling as well as further discusses current gaps and future prospects in the field. Detailed discussions of multiphysics coupling methods are beyond the scope of this paper; readers are referred to more comprehensive reviews (Lei et al., 2017; Viswanathan et al., 2022; Vaezi et al., 2026).

Substantial progress has been made in recent years on 3D DFN modelling of multiphysical processes in fractured geological media. For instance, 3D DFN models have been extensively applied to investigate fluid flow and solute transport in 3D fracture networks, demonstrating that the structural configuration and associated geometrical/topological properties (e.g., fracture density, length, and connectivity) play a critical role in driving the emergence of flow channelling and anomalous transport phenomena (Figure 1a) (Hyman et al., 2019; Hyman, 2020; Kang et al., 2020; Yoon et al., 2023; Davy et al., 2024). More advanced 3D models that incorporate aperture variation within individual fractures have further revealed that fracture-scale heterogeneity exerts a significant influence on flow and transport properties at the network scale (Zou and Cvetkovic, 2020; 2021; Hyman et al., 2021; Sweeney et al., 2023; Osuji et al., 2025). Inspired by observed flow channelling phenomena, graph-based reduction models or channel network models have been developed to represent the flow in 3D fracture networks using interconnected 1D conductance channels, significantly reducing computational demand while preserving the essential characteristics of dominant flow structures (Viswanathan et al., 2018; Berrone et al., 2020; Doolaeghe et al., 2020; Dessirier et al., 2023). Geomechanical simulations based on 3D DFNs have also been performed to study the relationship between fracture network properties (e.g., density, length, and connectivity) and the overall mechanical properties of rock masses (e.g., Young’s modulus, Poisson’s ratio, and stress heterogeneity) (Davy et al., 2018; Lavoine et al., 2024). Advanced geomechanics models have been developed to simulate fracture propagation and interaction in 3D rock masses (Paluszny and Zimmerman, 2011; Wang et al., 2018; Paluszny et al., 2020; Sun et al., 2024). Notably, 3D DFN models have been recently applied to investigate earthquake rupture and associated geomechanical responses in 3D fracture networks (Figure 1b) (Pan et al., 2023; Pan et al., 2024; Gabriel et al., 2024; Palgunadi et al., 2024). Coupled hydro-mechanical simulations have also been conducted to investigate the impact of stress on fluid flow and solute transport (Lang et al., 2018; Sweeney and Hyman, 2020; Darcel et al., 2024) as well as to derive bulk hydro-mechanical properties such as equivalent Biot and Skempton coefficients (De Simone et al., 2023b). Additionally, heat transport in fractured rocks has been simulated using 3D DFNs, with the results highlighting the strong influence of flow velocity heterogeneity and matrix diffusion on thermal behaviour (De Simone et al., 2021; De Simone et al., 2023a). Recent advances have also enabled coupled hydro-chemical simulations of reactive transport in 3D fracture networks, shedding light on the interplay between fluid flow, solute transport, and chemical reaction (Hyman et al., 2022; Hyman et al., 2024; Andrews et al., 2023). Furthermore, 3D DFN models have been applied to a range of geoengineering problems, including fluid injection-induced seismicity (Ucar et al., 2017; Ucar et al., 2018), nuclear waste disposal (Pan et al., 2023; Pan et al., 2024; Leone et al., 2025), slope stability (Bonilla-Sierra et al., 2015), and tunnel excavation (Wang and Cai, 2025) in fractured rock masses. Compared to their 2D counterparts, these 3D DFN models reveal significantly richer physical behaviours of fractured media in 3D space and offer deeper insights into 3D effects associated with structural configurations, fracture interactions, stress fluctuations, aperture variations, and flow channelling, which cannot be adequately captured by 2D models.

While substantial efforts have been made to incorporate site characterisation data to inform or parameterise 3D DFN models, the validation and calibration of 3D DFN-based multiphysics simulations using field monitoring data remain limited, with only a few recent attempts addressing this gap. For example, Zhao et al. (2023) conducted a comprehensive case study of the Gotthard Base Tunnel—a 57 km long, up to 2.5 km deep high-speed railway tunnel in the Swiss Alps. During the tunnel construction, significant ground surface displacements up to ∼10 cm occurred, threatening the nearby critical infrastructures like concrete arch dams. To understand the causal mechanism of these decimetre-scale displacements, a 3D DFN-based computational model has been developed (Figure 2a) to realistically represent the heterogeneous ground conditions constrained by extensive laboratory/site characterization datasets and mechanistically compute the coupled processes of fluid flow, ground deformation, and fault slip at the site scale. The simulation results showed an excellent match to the in-situ monitoring data (Figures 2b,c), illuminating that the surface displacements originate from tunnelling-induced water drainage and rock mass consolidation in the deep subsurface. In addition, efforts have been made to integrate 3D DFN models with field experiments from underground research laboratories, e.g., the Mont Terri Rock Laboratory (Zhao et al., 2024), the Grimsel Rock Laboratory (Ringel et al., 2021), the ONKALO Demonstration Area (Hartley et al., 2018), and the Äspö Hard Rock Laboratory (Cvetkovic and Frampton, 2010). These studies demonstrate the validity and broad applicability of the DFN method across diverse geological settings, including both crystalline and sedimentary formations. They also highlight the significant value of integrating 3D DFN models with field observational data, enabling model validation and calibration to assess performance and build confidence, while elucidating the multiphysical mechanisms underlying complex phenomena observed at the site scale.

Over the past years, various 3D DFN generators and simulators have been developed by different research groups (Davy et al., 2013; Hyman et al., 2015; Alghalandis, 2017; Paluszny et al., 2020; Keilegavlen et al., 2021; Zhu et al., 2022; Pan et al., 2023; Im and Avouac, 2024; Palgunadi et al., 2024), each with a distinct emphasis on some specific aspects. Recently, those teams that originally focused on different aspects are now increasingly working towards a common goal—developing 3D DFN models capable of capturing an expanding range of multiphysical processes, such as mechanical deformation, fracture propagation, fluid flow, heat transport, and chemical reactions. A major obstacle to this endeavour is the multifaceted computational challenge posed by 3D multiphysics simulations. This difficulty is rooted in the DFN modelling concept itself, which aims to explicitly represent individual fractures and their associated multiphysical processes in a complex network, given that fractures in rock are inherently multiscale geological features that span a broad range of length scales. This challenge becomes even more pronounced when attempting to couple different physical processes, such as seismic slip, fluid flow, and chemical reaction, which operate across vastly different time scales. Compromised solutions may need to be developed, such as hierarchical modelling of fracture networks, where the effects of small-scale fractures are incorporated into equivalent matrix properties, and sequential coupling, where different physical processes are simulated over separate time steps and then coupled. Apart from the multiscale complexity, this challenge also partly stems from the difficulty of automated mesh generation for stochastically generated 3D DFNs, which involve a wide range of fracture sizes and exhibit complex topologies. Algorithms should be developed to optimise 3D DFN geometries during model construction, aiming to eliminate features such as tiny dead-ends, narrow intersection angles, and closely spaced fractures, which can degrade mesh quality and lead to numerical instability.

Owing to the significant computational challenges associated with 3D DFN multiphysics simulations, various simplifications and assumptions are often introduced in previous work to reduce model complexity and ensure computational feasibility. For example, the matrix is usually neglected in many 3D modelling studies, assuming that flow and transport primarily occur through interconnected fractures (Hyman et al., 2019; Hyman et al., 2021; Kang et al., 2020; Zou and Cvetkovic, 2020; Zou and Cvetkovic, 2021; Sweeney et al., 2023; Yoon et al., 2023; Davy et al., 2024). Local stresses on individual fractures in 3D fracture networks have often been calculated by simply projecting far-field stresses onto local fracture planes (Ucar et al., 2017; Ucar et al., 2018; Sweeney and Hyman, 2020), without accounting for fracture interaction-induced local stress fluctuations. Furthermore, fracture constitutive behaviour has often been modelled using simplified formulations, such as linear elastic or perfectly elasto-plastic models (Davy et al., 2018; Ucar et al., 2018; Sweeney and Hyman, 2020; Darcel et al., 2024), which do not account for the strongly non-linear deformational responses of rough fractures under normal and shear stress loadings (Lei and Barton, 2022). These simplifications can significantly influence the calculation of fracture shear dilations, which play a critical role in the formation of flow channelling and the emergence of anomalous transport behaviour. Future efforts are needed to assess how these simplifications affect the processes we aim to model and capture, especially in the context of multiphysics coupling, while continuously advancing the realism of 3D DFN models through more refined assumptions and enhanced computational capacity and efficiency.

Most existing 3D DFN modelling studies have focused on theoretical or conceptual investigations aimed at uncovering the complex multiphysical processes in fractured geological media. However, a significant gap remains in translating these advances into practical engineering applications. Bridging this gap requires reliably conditioning 3D DFN models using limited site characterization data and rigorously calibrating them against comprehensive field observations. Common calibration methods include parameter sensitivity analysis, inverse modelling, and Bayesian inference (Dverstorp and Andersson, 1989; Somogyvári et al., 2017; Ringel et al., 2021; Jiang et al., 2023), which enable probabilistic updating of model parameters based on observed data. Data assimilation techniques, such as the Ensemble Kalman Filter, may also be applied to integrate real-time monitoring data with DFN simulations (Elahi and Jafarpour, 2018). A major challenge arises from the large number of input parameters for DFN modelling, which often exhibit complex correlations and nonlinear effects on model outputs. Differentiating epistemic uncertainty (due to incomplete knowledge) from aleatoric uncertainty (due to inherent variability) is essential for meaningful uncertainty quantification (Murph et al., 2024). Machine learning techniques can be employed to accelerate calibration and facilitate uncertainty quantification (Srinivasan et al., 2018; Srinivasan et al., 2021), thereby enhancing the practical applicability of DFN models.

Looking ahead, future DFN modelling research may expand along several other avenues. First, significant knowledge gaps remain regarding chemo-mechanical coupling in fracture networks, particularly in understanding how chemical reactions influence the mechanical properties of individual fractures (e.g., stiffness, friction coefficient) and bulk rock masses (e.g., modulus, strength). Addressing this issue will require integrated efforts combining numerical simulations with laboratory experiments. Secondly, while much of the existing DFN research focuses on quasi-static simulations of multiphysical processes, considering inertial effects becomes crucial when addressing problems related to earthquake triggering, seismic slip, and wave propagation. Significant efforts are needed to address the challenge of simultaneously modeling these transient dynamic processes alongside other more quiescent ones like fluid flow and chemical reaction. Thirdly, fracture systems in crustal rocks often exhibit hierarchical structures, such as large-scale fault zones with numerous subsidiary fractures (Faulkner et al., 2010). Investigating how the collective behaviour of these structures governs larger-scale phenomena could provide valuable insights into developing hierarchical modelling strategies, which may be further inspired by the renormalisation group concept (Sornette, 2006), to address the multiscale challenges inherent in DFN modelling. Lastly, modelling fracture propagation in 3D complex fractured media under coupled multiphysical conditions remains a key challenge and warrants further investigation.

Although DFN models have witnessed significant advances in recent years, longstanding criticisms remain and are likely to persist, concerning their complexity, representativeness, and the challenges associated with model parameterisation and calibration, which in turn raise questions about their reliability for practical use. This stems from the inherent dilemma of studying fractured media, which exhibit varying degrees of complexity across multiple scales. The multiscale nature of fractured media arises from two aspects: (i) presence of fractures spanning multiple length scales and (ii) emergent behaviour at the system level resulting from interactions among numerous components. In this context, “small-scale” refers to either features of limited spatial extent or the level of individual components, while “large-scale” denotes either broader spatial domains or the collective behaviour emerging at the system level. At one end, understanding small-scale details, including numerous fractures and masses, requires navigating a high degree of complexity. At the other, practical applications usually require predictions of only a few large-scale properties, such as modulus, strength, permeability, breakthrough, and reaction rate, for which reducing complexity is key. To address this dilemma, a hierarchical modelling strategy should be adopted, integrating models across complexity levels, from high-to intermediate- and low-complexity models. DFN-based multiphysics simulations fall in the category of high-complexity models that resolve detailed physical processes at high spatial and temporal resolutions, with significant computational and data requirements, providing deep insights into the mechanisms. Intermediate-complexity models, such as graph-based reduction approaches, simplify some of these processes to focus on key phenomena, lowering computational costs and data demands while still capturing essential system dynamics, which makes them particularly suitable for uncertainty analyses and Monte Carlo simulations. Low-complexity models, such as statistical or surrogate models, rely on highly simplified representations to capture general trends, offering rapid predictions and enabling robust calibration. Used in a complementary manner, these models offer a pathway to balance complexity and utility: to unravel the intricacies of fractured media while enabling practical, real-world applications. Lastly, it is worth reaffirming the spirit of the adage that “all models are wrong, but some are useful.” DFN models, despite their inherent complexities and uncertainties, can be extremely useful, especially in combination with other models of intermediate/low-complexity, for both gaining mechanistic understanding and informing engineering decisions. Extensive efforts are required to establish such a hierarchical framework to effectively reconciles model predictions with real-world observations, ultimately enhancing the reliability of DFN models for practical applications.