Fluid Flow in Fractured Rocks: From Multiphysics Paradigms to AI-Driven Predictive Modeling

1 Introduction

1.1 Research Background, Significance, and Core Challenges

Fractured rock masses, dissected by discontinuous structural planes such as joints, faults, and bedding planes, exhibit hydro-mechanical behaviors fundamentally distinct from homogeneous porous media. Their pronounced spatial heterogeneity, geometric discontinuity, and permeability anisotropy govern fluid transport and critically influence the stability and performance of geotechnical and energy system [1]. Accurately characterizing and predicting seepage in such media remains a persistent challenge in geoscience and rock engineering.

In practice, seepage processes often dictate the success or failure of major subsurface projects. For example, in deep tunnels and underground caverns, coupled seepage–stress effects can trigger large deformations or catastrophic water inrush [2]. Similarly, in EGS, unconventional hydrocarbon recovery, and CO2 sequestration, effective reservoir management hinges on understanding how THMC processes dynamically alter fracture networks and their permeability evolution [3].

Despite extensive research on individual aspects—ranging from conceptual modeling to numerical simulation—a comprehensive, up-to-date synthesis that integrates theoretical frameworks, experimental insights, numerical advances, and emerging AI-driven paradigms is still lacking. This review aims to bridge these gaps by: (i) systematically categorizing seepage models; (ii) evaluating experimental scale effects and THMC coupling mechanisms; (iii) comparing mainstream numerical methods; and (iv) critically assessing the role of AI, particularly PINNs. The paper is organized as follows: Section 2 presents the fundamental principles and governing equations of THMC coupling. Section 3 reviews modeling paradigms and numerical methods for seepage. Section 4 synthesizes experimental findings on seepage characteristics. Section 5 demonstrates practical relevance through engineering case studies. Section 6 explores the integration of AI and ML in subsurface flow modeling. Section 7 discusses persistent challenges and outlines future research directions. Finally, Section 8 concludes the review with key takeaways and a forward-looking perspective.

The fundamental challenge underlying these engineering practices lies in the fact that fractures—the dominant pathways controlling fluid flow—exhibit highly random and nonlinear geometric configurations, spatial distributions, and mechanical behaviors [4,5,6]. This renders the traditional equivalent continuum theory, which relies on the assumption of a “Representative Elementary Volume (REV) [7] “invalid in many practical scenarios. Crucially, the relationship between fracture permeability and effective stress is not merely linear; instead, it displays strong nonlinearity, hysteresis, and path dependence [8]. This implies that the hydraulic properties of the rock mass evolve along distinct trajectories during loading and unloading cycles—a complex behavior that cannot be captured by static or linear models. Consequently, advanced multi-field coupling theories and dynamic numerical simulation methods are essential for in-depth investigation [9].

1.2 Structure and Research Objectives of This Paper

This review aims to provide a systematic and critical synthesis of the current state of research on fluid flow in fractured rock masses and the underlying THMC multi-field coupling mechanisms. Its primary objective is not merely to catalog existing achievements, but rather to reveal the consensus, controversies, and emerging frontiers in the field by deeply analyzing the intrinsic connections and contradictions among theoretical models, experimental findings, numerical methods, and engineering practices. In doing so, it seeks to offer a robust theoretical reference and methodological guidance for future fundamental research and engineering applications.

To achieve this objective, the paper follows a clear and coherent logical structure: First, it begins with fundamental principles, articulating the physical mechanisms of THMC coupling that govern fractured rock behavior and their mathematical formulations (Section 2). Next, it systematically reviews prevailing theoretical modeling paradigms—such as ECM and DFN—alongside numerical implementation methods—including the Finite Element Method (FEM), Discrete Element Method (DEM), and Lattice Boltzmann Method (LBM)—and provides a comparative analysis of their applicability and limitations (Section 3). Building on this foundation, the discussion turns to experimental research, synthesizing key findings from single-fracture and triaxial seepage experiments to establish a robust empirical basis for theoretical models (Section 4). Subsequently, the paper examines three representative high-stakes engineering case studies—geological disposal of nuclear waste, EGS, and tunneling projects—to demonstrate how THMC coupling theory addresses cutting-edge, high-risk engineering challenges (Section 5). Following this, it explores how AI and ML are emerging as transformative forces that are reshaping research paradigms in the field, analyzing their synergies and complementarities with traditional approaches (Section 6). Finally, the paper confronts current bottlenecks and outlines future directions, offering forward-looking perspectives on critical issues such as multi-scale modeling, parameter uncertainty, and representation of complex fracture networks (Section 7).

2 Fundamental Principles and Governing Equations of THMC Coupling

2.1 Interlocking Mechanisms of THMC Coupling

The complete THMC behavior of fractured rock masses is fundamentally an interlocked system in which multiple physical fields are tightly intertwined and mutually driven [10]. Its defining characteristic lies not in the independent action of individual fields, but in the causal chains and dynamic feedback loops established among them through well-defined physical laws [11]. A perturbation in any one field inevitably propagates and amplifies through coupling mechanisms, ultimately reshaping the state of the entire system [12]. This intrinsic interdependence necessitates that accurate prediction of fractured rock behavior move beyond the simplified paradigm of single-field analysis and instead adopt a unified, fully coupled theoretical framework.

Within this framework, the four physical fields—thermal (T), hydraulic (H), mechanical (M), and chemical (C)—are interconnected through a series of core coupling mechanisms [13]:

• Mechanical-Hydraulic Coupling: This is the most fundamental and direct coupling mechanism. Changes in the mechanical field—such as stress redistribution caused by excavation or fluid injection—directly alter the physical aperture and connectivity of fractures, thereby drastically modulating the permeability and fluid flux of the hydraulic field. Conversely, changes in the hydraulic field—such as increases or decreases in pore pressure—feed back into the mechanical field via the effective stress principle (σ′=σ−αp), altering the stress state of the rock mass and potentially triggering further deformation or even fracturing [14]. This bidirectional interaction forms a classic positive or negative feedback loop, serving as a primary driver of the system’s dynamic evolution.

• Thermo-Mechanical Coupling: Changes in the thermal field induce thermal expansion or contraction of rock minerals, generating additional thermal stresses within the rock mass and consequently causing strain and displacement in the mechanical field. In EGS, the thermal stress induced by injecting cold fluid into hot rock is precisely the key physical mechanism that activates pre-existing fractures, promoting shear slip (hydraulic shearing) and thereby enhancing permeability [15].
• Hydro-Thermal Coupling: In fractured rock masses, heat transfer depends not only on thermal conduction through the rock matrix but, more significantly, on convective heat exchange by fluid flow. Key fluid properties—such as viscosity and density—are strong functions of temperature; consequently, the spatial distribution of the thermal field directly influences fluid flow patterns (hydraulic field), while fluid movement, in turn, redistributes heat and alters the thermal field itself [16]. This mutual modulation forms a complex, nonlinear feedback system.
• Chemo-Hydraulic Coupling: The chemical field primarily manifests through water–rock reactions. Dissolution or precipitation reactions between fluids and minerals microscopically alter the geometry of fracture walls—such as aperture and roughness—and even their mineralogical composition, thereby irreversibly modifying hydraulic conductivity. This process is particularly pronounced in reactive media such as bentonite buffer layers or shale, and represents a critical factor in assessing the long-term (millennia-scale) safety of nuclear waste repositories [17].

These coupling mechanisms do not operate in isolation; rather, they act synergistically and evolve dynamically. For example, in an EGS reservoir, a complete coupling cycle can be described as: hydraulic injection (H) → thermal gradient (T) → thermal stress (M) → fracture shear dilation (M–H) → increased permeability (H) → enhanced heat extraction efficiency (T) [18]. This chain clearly illustrates how an initial perturbation propagates through interlocked physical mechanisms, triggering a cascade of effects that ultimately reshape the system’s state. Fig. 1 further visualizes the core thermo-hydro-mechanical (THM) coupling paths in fractured rocks, explicitly depicting the dynamic feedback loops between hydraulic pressure, thermal gradient, and mechanical stress during fracture dilation and permeability evolution—directly echoing the EGS coupling cycle outlined above. It is precisely this highly nonlinear, path-dependent, and dynamically evolving nature that constitutes the fundamental challenge in studying seepage in fractured rock masses—and underscores the necessity of developing a comprehensive THMC coupling theory [19].

Figure 1: THM couplings in fractured rocks. Adapted with permission from Reference [20]. Copyright © 2025 Elsevier.

2.2 Governing Equations and Constitutive Relationships

The quantitative description of THMC multi-field coupled behavior in fractured rock masses relies on a set of interrelated partial differential equations. These equations respectively express the conservation laws governing the T, H, M, and C fields, and are unified into a single dynamical system through constitutive relationships [21]. This section systematically presents the core governing equations for each field and their key constitutive models, thereby establishing a solid theoretical foundation for subsequent numerical simulations and engineering analyses.

2.2.1 Mechanical Field

The mechanical response of rock masses is governed by the stress–strain constitutive relationship. For rock in an elastic or elastoplastic state, the relationship between the stress tensor σ i j and the total strain tensor ε k l can be expressed as: σij=Dijkl(εkl−εklp)(1) where D i j k l is the fourth-order elasticity tensor, and ε k l p is the plastic strain component (which can be neglected in purely elastic analyses).

Within a multiphysics coupling framework, the principle of effective stress serves as the fundamental link between the mechanical and hydraulic fields. This principle states that the effective stress σ ′ acting on the rock skeleton equals the total stress σ minus the contribution from pore fluid pressure p : σ′=σ−αp(2) where α is the Biot-Willis coefficient, with values ranging between 0 and 1, representing the efficiency with which pore pressure is transmitted to the rock skeleton. When α = 1 , pore pressure fully acts on the solid skeleton; when α < 1 , part of the pressure is borne by the fluid itself or by isolated, unconnected pores. This equation forms the physical foundation for understanding key phenomena such as “stress-induced permeability changes” [22] and “pore-pressure-induced rock deformation” [23].

2.2.2 Hydraulic Field

Fluid flow in fractured rock masses typically follows Darcy’s law [24], which establishes a linear relationship between the fluid flux (Darcy velocity) q and the pressure gradient ∇ p : q=−kμ∇p(3) where k is the permeability tensor of the rock mass (for anisotropic media), and μ is the fluid dynamic viscosity. It is important to note that permeability k is not a material constant; rather, it is a state-dependent variable that strongly depends on the effective stress state and fracture geometry.

However, in high-velocity scenarios such as EGS or during hydraulic stimulation, inertial effects become significant, leading to non-Darcy flow behavior. Under these conditions—typically when the Reynolds number exceeds 1–10—the linear assumption of Darcy’s law breaks down. To account for this, the Forchheimer equation is widely adopted as a physically grounded extension: −∇p=μkq+βρ|q|q(4) where ρ is fluid density and β is the Forchheimer coefficient, which depends on fracture roughness and aperture distribution. This formulation introduces a quadratic velocity term to capture inertial losses and has been validated in experimental studies of single fractures and fracture networks.

In principle, the Navier-Stokes equations provide the most complete description of fluid motion at the fracture scale, resolving local velocity fields, turbulence, and boundary layer effects. However, their direct implementation in large-scale, fully coupled THMC models remains computationally prohibitive due to the need for extremely fine mesh resolution near fracture walls and the strong nonlinearity introduced into the multi-physics system. Consequently, while Navier-Stokes offers theoretical completeness, the Forchheimer model represents a pragmatic compromise that captures essential non-Darcy physics while retaining numerical tractability in field-scale simulations.

Fluid mass conservation is described by the continuity equation: ∂(ϕρf)∂t+∇⋅(ρfq)=Q(5) where ϕ is porosity, ρ f is fluid density, t is time, and Q represents source/sink terms (e.g., injection or extraction). In most engineering applications, if the fluid is considered incompressible ( ρ f ≈ constant ), the equation simplifies to:

2.2.3 Thermal Field

Heat transfer in rock masses occurs primarily through two mechanisms: thermal conduction and thermal convection. Thermal conduction follows Fourier’s law [25], which states that the heat flux q T is proportional to the temperature gradient ∇ T : qT=−kT∇T(7) where k T is the thermal conductivity tensor of the rock mass.

The system’s total energy conservation equation must account simultaneously for conduction, convection, and possible heat sources: ρcp∂T∂t=∇⋅(kT∇T)−ρfcp,fq⋅∇T+QT(8) where ρ and c p are the equivalent density and specific heat capacity of the rock mass, respectively; ρ f and c p , f are the density and specific heat capacity of the fluid; and is the volumetric heat source term (e.g., heat from radioactive decay). The second term on the left-hand side, − ρ f c p , f q ⋅ ∇ T , represents the heat convection term, which explicitly captures the strong coupling between the hydraulic field ( q ) and the thermal field ( ∇ T ). In systems where flow is dominated by fracture networks, convection often becomes the primary mechanism for heat transport—frequently surpassing conduction in efficiency, especially under active fluid circulation such as in EGS.

Moreover, temperature changes influence the mechanical field through thermal expansion effects. This mechanism is incorporated by introducing a thermal strain component ε i j T into the total strain: εijT=αTΔTδij(9) where α T is the coefficient of thermal expansion, Δ T is the temperature change relative to a reference state, and δ i j is the Kronecker delta. This thermal strain component is directly incorporated into the stress–strain constitutive relationship, thereby closing the TM coupling loop.

Thus, the choice of flow model—whether Darcy, Forchheimer, or Navier-Stokes—must be guided by both the characteristic flow regime and the computational feasibility within the context of multi-field coupling simulations.

2.2.4 Chemical Field

The chemical field primarily manifests through water–rock reactions that induce long-term, irreversible alterations to the system’s physical properties. Although the governing equations—such as reactive transport equations—are not elaborated in detail in this review, their physical significance is critical [26]. Mineral dissolution enlarges pore or fracture space, thereby enhancing permeability, whereas precipitation clogs pore throats, reducing permeability and potentially increasing rock strength [27]. This process is especially pronounced in reactive media such as bentonite buffer layers, shale reservoirs, or host rocks surrounding nuclear waste repositories, and represents a central factor in evaluating the system’s long-term (millennia-scale) performance [28]. The chemical field typically couples indirectly with other physical fields by modifying state variables such as porosity ϕ , permeability k , or rock strength parameters.

2.3 Engineering Relevance of Coupling Mechanisms

Although the physical mechanisms of THMC multi-field coupling in fractured rock masses are universal, the dominant processes and key controlling factors vary significantly across different engineering contexts [29]. Accurately identifying and focusing on the core coupling pathways specific to a given application is a prerequisite for developing efficient and reliable numerical models, and is crucial for optimizing engineering performance and managing risks [30].

In EGS, thermo-mechanical-hydraulic (TMH) coupling serves as the core engine driving system performance [31]. The typical process chain can be summarized as: cold fluid injection H → localized cooling of the rock mass generates thermal stresses (T → M) → induces shear slip along pre-existing fractures or propagation of new fractures M → increases fracture aperture and permeability (M → H) → enhances fluid flux and heat extraction efficiency (H → T) [32,33,34,35]. This positive feedback loop constitutes the physical basis for economically viable heat extraction in EGS. However, it also carries the risk of “thermal breakthrough”—where excessively rapid fluid flow causes a sharp decline in reservoir temperature, thereby shortening the system’s operational lifespan [36]. Fig. 2 intuitively illustrates the stress and temperature dynamics on the fault plane of reservoirs intercalated with low permeability zones, covering both no fault offset (top row) and normal fault offset (bottom row) conditions, while quantifying the relative difference in shear stress predictions between MACRIS and uniaxial solutions over the geothermal doublet’s production lifetime. Consequently, the crux of engineering design lies in precisely managing the intensity and spatial extent of TMH coupling by controlling injection parameters such as pressure and flow rate, thus striking an optimal balance between maximizing heat extraction efficiency and ensuring long-term system stability.

Figure 2: Stress and temperature solutions on the fault plane for a reservoir intercalated with low permeability zones in the case of no fault offset (top row) and normal fault offset (bottom row). The relative difference in shear stress solutions between both approaches is presented over the production lifetime of the geothermal doublet. A positive difference indicates MACRIS to predict a larger stress response compared to the uniaxial solution. The reservoir geometry is outlined in grey/black. Adapted with permission from Reference [37]. Copyright © 2024 Springer.

In the geological disposal of high-level radioactive waste (HLW), hydraulic–mechanical–chemical (HMC) coupling becomes the decisive factor in assessing safety over millennial timescales [38]. After the engineered barrier system (waste canisters and bentonite buffer) initially achieves self-sealing and self-supporting functionality through THM coupling, the long-term safety of the repository relies primarily on the natural barrier properties of the host rock. At this stage, potential migration pathways for radionuclides are entirely governed by the fracture network in the surrounding rock. Long-term groundwater flow H through these fractures is not only influenced by M-controlled evolution of fracture apertures but is also deeply coupled with C-driven mineral dissolution and precipitation reactions. Water–rock interactions can either fill fractures (reducing permeability) or cause alteration and widening (increasing permeability), thereby irreversibly altering the rate and pathways of radionuclide transport [39]. Consequently, HMC coupling models must be capable of capturing these slow yet decisive long-term evolutionary processes, as their predictions directly determine the reliability of safety assessments for the repository.

In deep-buried tunnels and underground caverns, hydraulic–mechanical (HM) coupling is often the direct trigger of engineering hazards [40]. Excavation-induced stress redistribution in the surrounding rock M alters the permeability of the fracture network (M → H). Simultaneously, high-head groundwater H infiltrates the excavation zone through fractures, generating seepage forces and reducing effective stress (H → M), which may lead to large deformations of the surrounding rock, water inrush, or mud gushing disasters. In cold regions where artificial ground freezing (AGF) is employed during construction, TM coupling becomes the dominant mechanism: the advancing freezing front T drives moisture migration and ice lens formation, generating frost heave pressure (T → M) that stabilizes the rock mass. However, freeze–thaw cycles T also cause progressive rock deterioration and induce cyclic variations in the stress state of support structures M, posing a threat to their long-term stability.

In summary, although the fundamental physics of THMC coupling in fractured rock is universal, its engineering manifestation is highly context-dependent. As systematically compared in Table 1, the dominant coupling pathways, governing time scales, and critical control parameters differ markedly across EGS, nuclear waste disposal, and deep tunneling—reflecting their divergent primary objectives (heat extraction, long-term containment, and structural stability) and operational horizons (decades, millennia, and construction-phase control, respectively). While these scenarios share common multiphysics foundations, their risk drivers and modeling priorities are shaped by distinct dominant feedback loops. Successful engineering practice therefore hinges not on applying a universal model, but on selectively emphasizing the relevant subset of coupling processes that dictate system behavior in each specific scenario. This targeted modeling strategy ensures both computational efficiency and physical fidelity, ultimately strengthening the bridge between multiphysics theory and real-world subsurface engineering.

3 Modeling Paradigms and Numerical Methods for Seepage

3.1 Comparative Analysis of Mainstream Modeling Paradigms

The core challenge in modeling fluid flow in fractured rock masses lies in effectively representing their inherent discontinuity and strong heterogeneity. To address this challenge, the research community has developed multiple theoretical modeling paradigms, each based on distinct physical abstractions and mathematical assumptions, resulting in unique trade-offs among computational efficiency, representation accuracy, and suitability for specific applications. The fundamental differences among the three mainstream modeling paradigms—ECM, DFN models, and Dual-Medium models—including their computational efficiency—are summarized in Table 2. This section provides a systematic comparative analysis of these three dominant approaches.

3.1.1 ECM

The ECM is a “top-down” macroscopic modeling approach. Its core idea is to treat the intrinsically discontinuous fractured rock mass, at a sufficiently large scale, as a continuous porous medium with anisotropic hydraulic properties [40]. The theoretical foundation of this model is the “REV” hypothesis, which posits that when the scale of investigation is much larger than the characteristic fracture spacing, the hydraulic and mechanical properties of the rock mass become statistically homogeneous and can thus be described by an equivalent permeability tensor K that characterizes the overall seepage behavior [41]. This homogenization logic is intuitively illustrated in Fig. 3, which demonstrates how a complex fracture network is upscaled into a single equivalent continuous flow channel to achieve macroscopic representation of seepage characteristics.

Figure 3: Schematic diagram of equivalent model (Paterson, 1983). I is the unit length of one element, A is the area of cross-section, q is the flow rate, and q′ is the flow rate in equivalent channel. Adapted with permission from Reference [42]. Copyright © 2015 Elsevier.

Advantages:

• High computational efficiency: By simplifying the complex fracture network into a continuous field, the ECM can be seamlessly integrated with well-established numerical methods such as the FEM or Finite Difference Method (FDM). This makes it suitable for large-scale engineering simulations (e.g., at watershed or far-field scales) with relatively low computational cost.
• Conceptual simplicity and strong engineering practicality: Model parameters—such as equivalent permeability—are readily invertible and calibratable using field pumping tests or tracer tests, making the model intuitive and accessible for engineering practitioners.

Limitations:

• Fragility of the REV assumption: In rock masses with sparse, highly heterogeneous, or extremely anisotropic fracture distributions, a statistically REV may not exist. In such cases, the equivalent parameters lose physical meaning, leading to model failure and unreliable predictions.
• Limited ability to represent spatial heterogeneity: ECM inherently adopts an “averaging” approach to the fracture network, making it difficult to capture the controlling influence of localized high-permeability pathways (such as major fractures or faults) on fluid flow trajectories. This often leads to underestimation of extreme-value risks, such as preferential flow or rapid contaminant transport.
• Lack of path dependence: The model typically assumes that permeability is a single-valued function of the current stress state, making it incapable of capturing the hysteresis and irreversible evolution of permeability during stress loading–unloading cycles.

Therefore, the ECM is most suitable for rock masses with high fracture density and relatively uniform fracture distribution, or for preliminary large-scale, low-resolution feasibility studies and regional impact assessments in the early stages of engineering projects [43]. Its success critically depends on the validity of the REV assumption under the specific geological conditions of the site.

3.1.2 DFN Model

The DFN model adopts a “bottom-up” modeling paradigm. Its core principle is to explicitly represent each individual fracture in the rock mass as a discrete geometric entity in three-dimensional space (typically as planar or curved surfaces) [44]. The model is constructed directly from measured or statistically derived fracture parameters, such as location, orientation, size, aperture, and roughness, thereby providing the most physically faithful representation of the discontinuous structure of fractured rock. Within the DFN framework, fluid flow is assumed to occur predominantly within the fracture network, while the rock matrix is typically treated as a low-permeability or impermeable block [45]. This explicit characterization of flow pathways makes DFN an ideal tool for investigating localized seepage paths, solute transport, and preferential flow phenomena.

Advantages:

• High physical fidelity: The DFN model can accurately capture the geometry, spatial orientation, and interconnectivity of individual fractures, thereby revealing the “channeling effect” created by major fractures or faults [46], a feature that ECMs cannot resolve. Its results possess clear physical meaning and can be directly used to analyze key hydrogeological parameters such as flow path tortuosity and the proportion of dead-end fractures.
• Applicability to complex geological conditions: For rock masses with sparse fracturing, extreme anisotropy, or the presence of large-scale structures (e.g., faults), the DFN model does not rely on the REV assumption, thereby offering greater applicability and predictive reliability.

Limitations:

• High computational cost: Explicitly solving fluid flow within each individual fracture and their interactions results in a very high number of degrees of freedom, leading to substantial computational resource demands. This becomes a major bottleneck, particularly when simulating large-scale or high-density fracture networks.
• Stringent data requirements: Model accuracy is highly dependent on the quality and completeness of input fracture geometry and hydraulic parameters. Acquiring a sufficient volume of statistically representative field data—through methods such as borehole imaging, outcrop mapping, or geophysical inversion—is often costly and technically challenging. Uncertainties in these input data propagate directly into the simulation results, affecting their reliability.
• Difficulty in scale upscaling: Extrapolating fracture parameters from small scales (e.g., laboratory or borehole scale) to engineering scales (e.g., kilometer-scale reservoirs) encounters significant scale effects and issues of parameter representativeness, often introducing systematic errors.

In summary, the DFN model, with its unparalleled geometric fidelity, has become the method of choice for near-field, high-resolution analyses—such as seepage around tunnel linings or simulations of the stimulated zone surrounding wells in EGS. However, its high computational cost and strong dependence on high-quality input data limit its applicability to localized, high-precision simulations of critical zones, rather than large-scale, system-wide assessments [47].

3.1.3 Dual-Porosity/Dual-Permeability Model

The Dual-Medium Model (Dual-Porosity/Dual-Permeability Model) is a hybrid, compromise-oriented modeling paradigm designed to bridge the gap between the macroscopic simplification of ECM and the microscopic complexity of DFN models [48]. Its core principle conceptualizes fractured rock as comprising two interconnected yet physically distinct subsystems: a high-permeability, low-storage fracture network, and a low-permeability, high-storage rock matrix [49]. This model acknowledges that rapid fluid transport through fractures is the dominant flow process, while simultaneously accounting for the matrix’s long-term role as a reservoir for fluid storage and solute diffusion.

In the classical Dual-Porosity Model, fluid flow is assumed to occur only within the fracture network at the macroscopic scale, while the matrix is treated as a stagnant “block” that exchanges fluid with fractures solely through diffusion or very slow imbibition. This formulation is suitable for scenarios where the matrix permeability is extremely low—such as in shale or tight sandstone. The more advanced Dual-Permeability Model further allows for slow but non-negligible fluid flow within the matrix itself, thereby providing a more realistic representation of bidirectional mass exchange between the matrix and the fracture network [50].

Advantages:

• Balanced computational efficiency and fidelity: This near-linear scaling enables large-scale reservoir simulations (e.g., geothermal fields, CO2 storage) where rapid scenario testing is essential.
• Multi-scale fracture handling: Aggregates structural and micro-fractures into a single fracture continuum, parameterizing multi-scale effects computationally tractably.
• Extensibility: Enhanced variants (e.g., EDPDP) incorporate stress-dependent apertures and nonlinear flow laws without significant computational overhead.

Limitations:

• Oversimplified fracture geometry: Loses individual fracture connectivity, making it incapable of capturing channeling effects (unlike DFN).
• Transfer function uncertainty: Empirical “shape factors” are highly sensitive to calibration data, introducing significant prediction errors.
• Anisotropy representation: Anisotropic tensors lack direct physical linkage to fracture geometry (unlike DFN), reducing accuracy in strongly anisotropic systems.

3.2 Evaluation of Numerical Modeling Approaches

Theoretical modeling paradigms provide a conceptual framework for fluid flow in fractured rock masses, while their practical engineering implementation heavily relies on efficient and robust numerical methods. Different numerical algorithms each possess distinct advantages and limitations in handling geometric complexity, physical nonlinearity, and multi-physics coupling. A comparative analysis of these mainstream numerical methods is presented in Table 3. This section systematically evaluates four widely used numerical approaches—FEM, DEM, Boundary Element Method (BEM), and LBM—to inform the synergistic selection of models and algorithms.

3.2.1 FEM and FDM

The FEM and the FDM are the most widely used numerical techniques for solving partial differential equations, particularly dominating the field of geotechnical engineering [51]. Both methods discretize the continuous solution domain into a finite number of elements or grid points, transforming the governing equations into large systems of algebraic equations for numerical solution.

The core strength of FEM lies in its ability to handle complex geometric boundaries and strongly coupled multi-physics processes. Through flexible mesh discretization, FEM can accurately conform to irregular engineering structures such as tunnels and caverns. Moreover, it naturally integrates the governing equations of THMC fields within a unified variational framework, enabling fully coupled THMC simulations [52]. Additionally, FEM offers well-established formulations for modeling material nonlinearity, such as elastoplastic constitutive behavior.

However, FEM faces two major challenges: First, when simulating large deformations or fracture propagation, the mesh can suffer severe distortion, leading to numerical instability and necessitating frequent remeshing, which incurs high computational costs [53]. Second, accurately resolving strong discontinuities such as fractures requires extremely fine meshes, significantly increasing the number of degrees of freedom [54]. To address these limitations, hybrid FEM-DEM approaches have been developed—employing FEM in macroscopically continuous regions while embedding DEM in localized fracture zones—thereby enabling a unified simulation of both continuous deformation and discontinuous failure.

Although the FDM is conceptually simpler and computationally efficient on regular grids, it lacks flexibility in handling complex geometric boundaries and struggles to achieve natural coupling of multi-physics processes. Consequently, its application in high-fidelity simulations of fractured rock masses has been increasingly superseded by the FEM.

3.2.2 DEM

The DEM fundamentally abandons the continuum assumption, treating rock masses as assemblies of discrete blocks or particles that interact through contact mechanics models—such as spring-dashpot systems [55]. This intrinsic feature makes DEM naturally suited for simulating discontinuous processes, including large deformations, progressive failure, and fracture evolution.

In fractured rock research, DEM can explicitly reproduce fracture opening, closure, sliding, as well as the initiation and coalescence of new fractures—without requiring predefined fracture paths. This offers a unique perspective for investigating engineering problems involving strong discontinuities, such as hydraulic fracturing, rockbursts, and slope instability. Its micromechanical framework is transparent, enabling a direct link between macroscopic responses and the evolution of microstructural features.

However, DEM also has significant limitations. First, its macroscopic mechanical parameters—such as elastic modulus and strength—are not intrinsic material properties; instead, they are derived through a calibration process that inversely maps microscopic contact parameters to macroscopic behavior. This process is time-consuming and often non-unique [56]. Second, the computational cost of DEM scales exponentially with system size, posing a severe computational bottleneck for large-scale engineering applications. Consequently, DEM is best suited for small-scale, mechanism-focused studies or, alternatively, as a complement to FEM for high-fidelity analysis of localized high-risk zones.

3.2.3 BEM

The BEM is a dimensionality-reduced numerical technique based on boundary integral equations. Its core idea is to discretize only the boundaries of the solution domain—including fracture surfaces—without requiring meshing of the interior region [57]. This feature enables BEM to achieve exceptional computational efficiency when dealing with problems characterized by a high surface-area-to-volume ratio.

For rock masses where flow is predominantly channeled through fracture networks, BEM can accurately represent fracture geometry and their hydraulic interactions with far fewer degrees of freedom, thereby avoiding the redundant computational effort associated with interior meshing in FEM. Moreover, BEM inherently satisfies far-field (infinite-domain) boundary conditions, giving it a distinct advantage in modeling unbounded problems such as far-field seepage.

However, the application of BEM is also constrained by its inherent limitations. First, its fundamental solutions rely on assumptions of homogeneous, linear-elastic media, making it difficult to effectively handle strongly heterogeneous materials—such as those with multiple lithological layers or complex rock matrices [58]. Second, numerical difficulties arise when dealing with low-angle intersecting fractures due to singularities in the integral equations [59]. Finally, for large-scale problems, BEM typically yields a dense system matrix, whose storage and solution costs increase rapidly with problem size, thereby limiting its applicability to very large simulations [60].

3.2.4 LBM

The LBM is a mesoscale numerical approach rooted in statistical mechanics, which solves the Navier-Stokes equations by simulating the collision and streaming of fluid particles on a regular lattice [61]. LBM offers unique advantages in modeling nonlinear fluid flow within complex geometric boundaries.

For fractured rock masses, LBM can directly simulate flow over digitized, rough fracture walls without relying on simplified geometric assumptions, accurately capturing complex flow phenomena such as eddies, backflow, and flow separation induced by surface roughness. Numerous studies have shown that permeability values obtained from LBM simulations are often significantly lower than those predicted by the cubic law based on the idealized parallel-plate assumption, revealing a systematic bias in traditional models under high flow velocities or complex geometries [62]. Furthermore, LBM demonstrates strong potential in handling complex flow scenarios, including multiphase flow and non-Newtonian fluids.

However, despite its high physical fidelity at the pore and fracture scale, the LBM confronts considerable limitations in computational efficiency and scalability for engineering-scale applications. These constraints arise from several inherent features of the method:

Explicit time integration: LBM relies on an explicit time-marching scheme subject to a stringent stability criterion, which typically demands very small time steps in diffusive regimes. Consequently, it requires significantly more time steps than implicit FEM solvers to simulate the same physical duration.

Substantial memory requirements: Each lattice node must store multiple distribution functions, leading to a large memory footprint per node. In contrast, FEM generally stores only primary field variables along with sparse matrix data, resulting in considerably lower memory usage for comparable spatial resolution.

Computational cost: Although LBM exhibits near-linear scaling per time step with respect to the number of lattice nodes, the overall computational expense grows rapidly with both physical time and domain size. For realistic fracture-scale simulations, this can lead to prohibitive runtimes, hindering parametric analyses or uncertainty quantification.

Parallel scaling limitations: While LBM demonstrates strong scaling on GPU-based systems due to its localized update rules, its weak scaling tends to degrade beyond certain problem sizes because of inter-process communication overhead. This makes it less suited for extreme-scale reservoir simulations compared to FEM with advanced solvers or other particle-based methods with efficient domain decomposition.

Hardware dependence: Achieving practical performance in LBM simulations strongly depends on GPU acceleration. CPU-only implementations are markedly slower, rendering them infeasible for large three-dimensional fracture networks beyond laboratory scale.

Moreover, current LBM implementations primarily focus on hydraulic field modeling; efficiently and robustly integrating LBM into a full THMC multiphysics framework remains a key direction for future research. Consequently, while LBM excels in high-fidelity microscale analysis of single fractures or small fracture clusters, its computational demands currently preclude routine use in large-scale, field-level THMC simulations—unless hybrid strategies (e.g., LBM for critical zones + FEM for far-field) are adopted.

Despite its promising potential, LBM still faces challenges in engineering applications. Algorithms for handling complex boundary conditions—such as dynamically moving boundaries—entail substantial computational overhead, and simulation techniques for time-evolving geometries (e.g., shear-induced slip or chemical dissolution) remain under active development. Moreover, current LBM implementations primarily focus on hydraulic field modeling; efficiently and robustly integrating LBM into a full THMC multiphysics framework remains a key direction for future research.

3.3 Model Selection Strategy and Applicability

In the study of fluid flow through fractured rock masses, the choice of theoretical model is not merely a technical issue but a strategic decision that can determine the success or failure of a project. The core principle lies in identifying and focusing on the dominant physical processes specific to the engineering problem at hand, and seeking an optimal balance among computational cost, data availability, and prediction accuracy. Fig. 4 intuitively demonstrates the practical value of simplified models in this balance: it presents a reference case of a homogeneous isotropic fracture network (a typical simplified idealized model), which can efficiently output key engineering-relevant results such as solute breakthrough curves, pressure difference distributions, and solute concentration fields—providing reliable macroscopic trend assessments for preliminary feasibility studies. An effective decision-making framework should integrate three key dimensions: project scale, data availability, and required accuracy, ensuring that the selected model is both computationally efficient and capable of delivering reliable, engineering-relevant results.

Figure 4: A reference case of isotropic fracture network: (a) aperture distribution with k = 1 and uncorrelated aperture distribution and (b) solute breakthrough curve (BTC) with its median travel time (time at 50% of the BTC), (c) pressure difference distribution in the fracture network above the exit-face pressure, and (d) solute concentration in the network at 0.5 years. Adapted with permission from Reference [63]. Copyright © 2025 Springer.

Decision Framework:

• (1)Project Scale and Spatial Extent:

• Large-scale/far-field analyses (e.g., regional groundwater flow, far-field impact zones of nuclear waste repositories): ECM are preferred. They offer high computational efficiency and are suitable for simulations at the kilometer scale, providing macroscopic, homogenized flow field distributions that meet strategic-level assessment requirements.
• Medium-scale/reservoir-scale problems (e.g., geothermal fields, oil and gas reservoirs): The dual-porosity/dual-permeability model is an ideal choice. It strikes a favorable balance between computational efficiency and physical fidelity, effectively capturing mass exchange between the rock matrix and fracture network, and is well-suited for complex systems with multiscale fracture networks.
• Small-scale/near-field analyses (e.g., tunnel surrounding rock, EGS wellbore stimulation zones, near-field of nuclear waste canisters): DFN models are essential. Only DFN can accurately capture critical mechanisms such as localized high-permeability pathways, fracture connectivity, and shear-induced slip, providing microscale insights necessary for high-fidelity design and safety assessment in high-risk zones.

• (2)Data Availability and Quality:

• Rich and high-quality data: When detailed fracture characterization data are available—such as borehole imaging, outcrop statistics, or geophysical inversion results—the DFN model can fully leverage its advantages, delivering high-fidelity predictions of significant engineering value.
• Limited or moderate-quality data: Under conditions of data scarcity or high uncertainty, ECM or dual-porosity models are more robust. They exhibit lower sensitivity to input parameters and can be calibrated through macroscopic field observations—such as pumping test data—thereby reducing the risk of model failure due to insufficient or uncertain data.

• (3)Required Accuracy and Engineering Objectives:

• High-accuracy, mechanism-focused studies: When the objective is to gain deep insight into flow pathways, solute transport mechanisms, or to assess localized hazard risks (e.g., water inrush or mud bursting), computational efficiency must be sacrificed in favor of DFN models to achieve the most physically realistic simulation of the underlying processes.
• Moderate accuracy, system-level evaluation: For optimizing system performance (e.g., EGS heat extraction efficiency or oil/gas recovery rates) or conducting comparative evaluations of multiple design scenarios, dual-porosity models are generally sufficient to provide reliable decision support.
• Low accuracy, feasibility studies: During early project stages or large-scale reconnaissance, ECM are the preferred choice for initial screening and macroscopic trend assessment due to their speed and low computational cost.

The successful application of any theoretical model fundamentally depends on the close alignment between its core assumptions and the actual geological conditions of the target rock mass. For example:

• The ECM relies on the existence of a REV. If applied to sparsely fractured or highly anisotropic rock masses where REV does not exist, the equivalent hydraulic parameters lose physical meaning, leading to predictions that significantly deviate from reality [64].
• The DFN model assumes that fluid flow occurs primarily within fractures, with the rock matrix treated as impermeable. If the matrix permeability of the target rock mass is non-negligible (e.g., in shale), a dual-permeability model must be adopted; otherwise, the storage capacity and diffusive contribution of the matrix will be significantly underestimated [65].
• The dual-porosity model simplifies the exchange process between fractures and matrix through a “transfer function”. In regions with extremely complex fracture network topology or dominant preferential flow pathways, this simplification may lead to misidentification of critical flow channels [66].

Therefore, model selection should not be a rigid, one-size-fits-all application, but rather a dynamic, iterative process of adaptation. Engineers must deeply understand the physical underpinnings and inherent limitations of each model and, in conjunction with thorough site-specific geological investigations, carefully evaluate the validity of model assumptions in the given context. Only in this way can theoretical models be effectively transformed into powerful tools for solving real-world engineering challenges—rather than becoming “black boxes” that produce misleading results.

4 Experimental Foundations of Seepage Characteristics

4.1 Single-Fracture Flow Experiments

Single-fracture flow experiments form the cornerstone for understanding the hydro-mechanical coupling behavior of fractured rock masses. By applying controlled normal stress, shear displacement, and fluid pressure to a single fracture—either artificial or natural—within rock specimens, researchers can isolate the effects of complex fracture networks and directly observe and quantify the constitutive relationships among stress state, geometric morphology, and fluid flow. Extensive experimental data consistently demonstrate that fracture permeability is not a material constant; rather, it is a highly nonlinear, path-dependent state variable jointly governed by effective stress and shear deformation [67]. Fig. 5 illustrates the evolution of hydraulic and maximum apertures with applied displacement during the simulation of a brittle rock. The hydraulic aperture, a key parameter directly correlating with fracture permeability, exhibits a distinct nonlinear response to shear displacement, further validating the path-dependent nature of permeability in fractured rock masses.

Figure 5: Hydraulic and maximum apertures versus applied displacement obtained in the simulation of a ‘brittle’ rock. Hydraulic aperture obtained with for flow in x-direction. Adapted with permission from Reference [68]. Copyright © 2016, Springer.

4.1.1 Nonlinear Control Mechanism of Effective Stress

Experimental results clearly confirm that as normal effective stress increases, the mean aperture of a fracture decreases nonlinearly while the variance of the aperture distribution increases, causing permeability to drop sharply—often by several orders of magnitude [69]. This phenomenon arises from the progressive increase and compaction of contact points between fracture walls. Mathematically, this relationship is commonly described by exponential decay or power-law functions. More importantly, the process exhibits pronounced hysteresis: during unloading, the recovery of permeability is significantly less than the loss observed during loading, indicating that part of the deformation is irreversible. This nonlinearity and hysteresis fundamentally challenge the applicability of traditional linear elastic models, necessitating the incorporation of damage or plasticity mechanisms into theoretical models to accurately capture the evolutionary history of fracture permeability.

4.1.2 Three-Stage Evolution of Shear Displacement and Roughness

When a fracture is subjected to shear displacement under constant normal stress, its permeability evolution follows a characteristic three-stage pattern:

• (1)Initial compression stage: Minor shear displacement induces “shear closure” of the fracture walls, increasing the contact area and causing a slight decrease in permeability.
• (2)Rapid dilation stage: As shear displacement further increases, microscopic asperities on the fracture walls undergo shear failure and ride-over, triggering macroscopic “shear dilation”. This leads to a significant increase in aperture and an exponential rise in permeability.
• (3)Stable residual stage: After dilation peaks, the fracture walls either become smoother or establish new stable contacts, causing the rate of permeability increase to slow down and eventually stabilize.

Fracture roughness is the key parameter governing this process. Experiments show that rougher fracture walls exhibit greater dilation potential and achieve higher peak permeability. This “dilation-induced permeability enhancement” mechanism forms the physical basis for artificially activating reservoir permeability in engineering applications such as EGS. Accurate prediction of this behavior is therefore critical to project success [70].

4.1.3 Non-Darcy Flow Behavior and the Forchheimer Equation

Under high flow velocities or large aperture conditions, fluid flow through fractures deviates from the classical linear Darcy’s law and exhibits non-Darcy behavior dominated by inertial effects. Typical experimentally observed phenomena include a nonlinear relationship between flow velocity and pressure gradient, the emergence of localized vortices and recirculation zones, and a significant increase in kinetic energy losses. To accurately describe such complex flow regimes [71], the Forchheimer equation has been widely adopted: −∇p=μkq+βρq2(10) where the first term represents the viscous resistance (Darcy term), and the second term represents the inertial resistance, where β is the non-Darcy (Forchheimer) coefficient and ρ is the fluid density. Extensive experimental data confirm that the Forchheimer equation effectively fits the pressure–flow rate relationship under high-velocity conditions, with significantly better goodness-of-fit than idealized models such as the cubic law. This finding underscores that in engineering scenarios involving high-pressure injection or high-velocity flow—such as hydraulic fracturing or EGS production—nonlinear flow models must be employed; otherwise, pressure losses and energy dissipation will be severely underestimated, potentially leading to flawed or failed designs. This transition into nonlinear flow regimes is critically dependent on permeability heterogeneity. As shown in Fig. 6, when the permeability ratio between fault and surrounding rock exceeds a threshold (e.g., ≥10), fluid velocity within the fault increases sharply, leading to rapid velocity variations and turbulent-like flow behavior—hallmarks of inertial-dominated non-Darcy flow. Under such conditions, the fault abruptly shifts from a barrier to a high-conductivity conduit, and the entrainment of sediments by high-velocity fluid can trigger catastrophic water inrush and mud outburst during underground excavation. This phenomenon underscores that non-Darcy effects are not merely a laboratory curiosity but a decisive factor in field-scale engineering safety.

Figure 6: Velocity distribution of monitoring line at different permeability ratios. Adapted with permission from Reference [72]. Copyright © 2019, MDPI.

In summary, single-fracture experiments not only provide essential constitutive parameters but also profoundly reveal the core physical mechanisms governing fracture flow: its behavior is a dynamic outcome of the interplay among stress history, geometric morphology, and flow regime. These insights offer indispensable physical constraints and validation benchmarks for developing high-fidelity theoretical models and numerical simulations. The key influencing factors and their governing relationships have been synthesized in Table 4.

4.2 Triaxial Seepage Experiments

While single-fracture experiments reveal fundamental stress–seepage coupling mechanisms, triaxial seepage experiments—by applying an anisotropic stress field (i.e., three unequal principal stresses, σ 1 , σ 2 , and σ 3 )—more realistically replicate the complex in-situ stress conditions of the deep crust. This enables the investigation of permeability evolution in fractured rock masses under multidirectional stress states and its associated engineering implications.

4.2.1 Exponential Decay under Anisotropic Stress Fields

Under triaxial stress conditions, the permeability of fractured rock masses exhibits strong directional dependence and nonlinear decay. Experimental data consistently show that as the maximum principal stress (typically the axial stress, σ 1 ) increases, permeability in that direction decreases exponentially [73]. This behavior arises from fracture closure and compaction under high confining pressure, which shifts the aperture distribution toward smaller values and drastically reduces the number of effective flow pathways. Notably, even at identical mean stress levels, different stress paths—such as isotropic compression versus deviatoric loading—lead to markedly distinct permeability evolution trajectories, further underscoring its path dependence. Moreover, experiments reveal that the permeability of fractured specimens is typically nearly an order of magnitude higher than that of intact rock, highlighting the dominant role of fractures as primary conduits for fluid flow.

4.2.2 The “Tensile Effect” of Lateral Stress

A key and counterintuitive finding is that lateral stress (confining pressure σ 3 ) can, under certain conditions, enhance fracture permeability rather than suppress it. Triaxial experiments on fractured granite show that, when axial stress is held constant, increasing lateral confining pressure can cause specific fractures to open, thereby increasing permeability. This phenomenon is known as the “tensile effect”. Its physical origin lies in the fact that, under an anisotropic stress field, the normal stress component acting on a fracture plane is not solely determined by the axial load but by the vectorial superposition of all three principal stresses [74]. For fractures with certain orientations, an increase in lateral stress can generate a net tensile stress component along the fracture normal direction, partially or fully counteracting the compressive effect of axial loading and leading to fracture dilation. Fig. 7 further validates this mechanism by presenting the variation of equivalent permeability components ( k x x , k y y , k z z ) of a fractured rock layer under 27 polyaxial stress scenarios ( σ ′ x , σ ′ y , σ ′ z = 5, 10, or 15 MPa). It can be observed that for specific stress combinations (e.g., when lateral stress σ ′ y or σ ′ z increases while axial stress σ ′ x remains constant), the permeability in the corresponding direction ( k y y or k z z ) increases rather than decreases—directly reflecting the tensile effect induced by the vectorial superposition of principal stresses. This discovery fundamentally refutes the oversimplified notion that “confining pressure always reduces permeability” and underscores the necessity of incorporating the full three-dimensional stress tensor and its spatial coupling with fracture geometry in modeling.

Figure 7: Variation of the equivalent permeability components (a) k x x , (b) k y y , and (c) k z z  of the fractured rock layer under various polyaxial stress conditions ( σ ′ x , σ ′ y , σ ′ z  = 5, 10 or 15 MPa, i.e., 27 stress scenarios). Adapted with permission from Reference [75]. Copyright © 2017, Springer.

4.2.3 Stress-Induced “Channeling Effect” and Its Engineering Implications

One of the most profound insights from triaxial experiments is the direct observation and validation of the stress-induced “fluid channeling” effect. When the deviatoric stress ( σ 1 − σ 3 ) is sufficiently high, certain dominant fractures undergo significant shear dilation due to shear slip, forming high-conductivity flow channels, while other fractures tend to close. This heterogeneous deformation causes fluid to no longer flow uniformly through the entire specimen but instead become highly concentrated in a few preferential pathways. Micro-scale flow visualization experiments clearly capture the transition of fluid pathways from a diffuse distribution to highly localized channels [76].

This effect carries significant engineering implications:

• In oil/gas and geothermal development, the “channeling effect” is a double-edged sword. On one hand, it can be deliberately harnessed to enhance reservoir permeability; on the other hand, excessive channeling can lead to “thermal breakthrough” or “gas coning”, drastically shortening the operational lifespan of the system.
• In tunneling and underground engineering, stress redistribution induced by excavation can activate fracture networks oriented in specific directions, creating sudden, high-conductivity water inrush pathways that may result in catastrophic failures.
• In nuclear waste disposal, long-term TM perturbations may alter the connectivity of the host rock fracture network, causing unpredictable redirection of radionuclide migration pathways.

Therefore, triaxial seepage experiments not only provide critical constitutive parameters but also reveal the systemic, heterogeneous, and path-selective nature of fluid flow in fractured rock masses. These findings serve as a clear warning: any model that neglects the interaction between stress anisotropy and the spatial geometry of fractures will fail to capture the “channeling effect”—a core physical process that often determines the success or failure of engineering projects.

4.3 Challenges in Scaling from Laboratory to Field Conditions

Extrapolating high-precision experimental data obtained at the laboratory scale to engineering field scales remains a long-standing and unresolved core bottleneck in the study of fluid flow in fractured rock masses. Numerous studies have confirmed that key parameters—such as permeability and reaction rates—measured in the lab often differ by orders of magnitude from those observed in situ or predicted by field-scale simulations. This “scale effect” is not merely a measurement error but stems from fundamental disparities between laboratory and field conditions in terms of system complexity, boundary constraints, and temporal scales [77].

4.3.1 Root Causes of Orders-of-Magnitude Discrepancies

Laboratory experiments—whether single-fracture tests or small-scale triaxial core tests—typically involve sample volumes ranging from a few cubic centimeters to a few liters. While such micro- or meso-scale experiments enable precise control of variables and reveal fundamental physical mechanisms (e.g., stress–permeability relationships, shear dilation), their results inherently describe an “idealized” or “localized” system. In contrast, in-situ rock masses constitute highly heterogeneous, multiscale, and dynamically coupled systems that encompass millions of fractures, span kilometers in spatial extent, and have evolved over millions of years [78]. Laboratories cannot replicate the following critical field-scale characteristics:

• Topological complexity of fracture networks: In the field, fractures are randomly distributed in three-dimensional space, intersecting and forming intricate, interconnected networks. Single-fracture laboratory experiments entirely neglect fracture–fracture interactions, while small-scale core tests are generally unable to capture large-scale structural features—such as dominant flow channels or dead-end clusters—that govern system-wide flow behavior.
• Irreproducibility of geological history: The current state of in-situ rock masses is the product of a long and complex geological history—including tectonic events, fluid–rock interactions, and stress evolution. Laboratory tests can only simulate a specific, instantaneous loading path and cannot replicate this prolonged, multiphysics evolutionary history.
• Significant disparity in boundary conditions: Laboratory tests employ artificially imposed and highly idealized boundaries (e.g., constant stress or constant flow rate), whereas in-situ rock masses exist within complex, evolving in-situ stress fields, temperature fields, and hydraulic gradients, where boundary conditions vary dynamically in both space and time.

4.3.2 Linking Laboratory Observations to Field-Scale Discrepancies

The triaxial seepage experiments reveal a pronounced stress-dependent anisotropy in permeability: under differential stress ( σ 1 > σ 3 ) , the principal permeability direction rotates and aligns preferentially with the direction of minimum compressive stress, while the magnitude of permeability perpendicular to the maximum stress can decrease by up to two orders of magnitude. While these trends are physically consistent with field observations—such as directional flow dominance in fault zones or stress-aligned fracture reopening—they also highlight a critical limitation. The absolute permeability values measured in our centimeter-scale specimens (typically on the order of 10−18–10−16 m2) are consistently lower than those inferred from pump tests or tracer experiments in similar rock formations at the meter-to-decameter scale (often 10−15–10−13 m2).

This discrepancy is not merely quantitative but reflects a fundamental representational gap: our laboratory samples capture only the matrix and microfracture response within a localized volume, whereas field-scale permeability is dominated by sparse, highly conductive fracture corridors that are statistically unlikely to be intersected by a small core. Consequently, while the directional sensitivity of permeability to stress observed in the lab provides valuable insight into the orientation of potential flow pathways at larger scales, the magnitude must be upscaled using network-based approaches that account for the probability of intersecting major discontinuities. In this sense, the experimental data serve as a constitutive foundation for microscale behavior, but cannot alone predict macroscale transmissivity without explicit integration into a multiscale framework.

4.3.3 The Theoretical Gap from Single Fractures to Complex Networks

Directly “linearly superimposing” or “averaging” single-fracture constitutive relationships—such as the cubic law or the Forchheimer equation—to predict the behavior of complex fracture networks is theoretically invalid [79]. This is because:

• (1)Nonlinearity and path dependence are amplified: The stress–permeability relationship of a single fracture is already nonlinear and path-dependent. In a complex network, deformation of one fracture redistributes stress to neighboring fractures, causing the entire system to exhibit even stronger nonlinearity, hysteresis, and unpredictability. Consequently, the effective permeability of the network is not a simple summation of individual fracture permeabilities.
• (2)Emergence of collective phenomena: Complex fracture networks exhibit macroscopic behaviors that do not exist in isolated single fractures—most notably the “fluid channeling” effect. Under stress, certain fractures dilate and become high-conductivity pathways while others close, causing flow to concentrate intensely along a few dominant channels. This results in pronounced anisotropy and strong non-Darcy flow characteristics at the system scale—phenomena that cannot be directly predicted from single-fracture experiments.
• (3)The mathematical challenge of scale bridging: Deriving macroscopic equivalent constitutive relationships for rock masses from the physical laws governing individual fractures (microscale) remains a central challenge in multiscale modeling. Traditional homogenization theories often fail when applied to strongly heterogeneous and discontinuous media. There is an urgent need to develop new mathematical frameworks—such as nested coupling strategies or data-driven cross-scale correlation models—to rigorously link processes across scales.

Therefore, bridging the scale gap between laboratory and field cannot be achieved merely by enlarging sample size or improving measurement precision. It requires a dual-pronged approach combining theoretical innovation and technological integration: on one hand, developing multiscale coupling theories capable of linking microscale mechanisms to macroscale behavior; on the other, leveraging AI and ML to learn from sparse field data and constrain model parameters, thereby enabling an effective transition from “laboratory physics” to “field engineering”. Successfully addressing this challenge will be a pivotal step in advancing fractured rock seepage research from mechanistic understanding toward precise, predictive capability.

5 Case Studies in Engineering Applications

5.1 Bridging Scales: From Laboratory Constitutive Laws to Field-Scale Fracture Network Models

The engineering applications discussed—nuclear waste disposal, EGS, and deep tunneling—all rely fundamentally on accurate predictions of fluid flow through complex fracture networks under coupled THMC conditions. However, the constitutive relationships governing fracture-scale seepage (e.g., nonlinear Forchheimer flow, stress-dependent aperture-permeability laws) are primarily derived from controlled laboratory experiments on single fractures or triaxial rock cores. A critical challenge lies in how to upscale these microscale, often highly nonlinear, relations to inform macroscale models of heterogeneous fracture systems. Two complementary strategies have emerged to address this scale gap.

First, mathematical homogenization theory provides a rigorous framework to derive effective continuum properties from local physics. By solving boundary-value problems over REVs of a fracture network—subject to lab-calibrated constitutive laws—one can obtain an equivalent, pressure- and stress-dependent permeability tensor k e f f ( p , σ ) . This approach is particularly suitable for densely fractured media where an REV exists, enabling efficient THMC simulations at the repository or reservoir scale while preserving key nonlinearities from laboratory data.

Second, when fracture networks are sparse or highly anisotropic (e.g., in crystalline host rock for nuclear disposal or EGS reservoirs), the REV assumption breaks down. In such cases, nested multi-scale coupling offers a more physically faithful alternative: explicitly modeled fractures in a DFN are assigned constitutive parameters directly calibrated from laboratory experiments (e.g., shear-dilation-permeability curves from triaxial tests), while background matrix or minor fractures may use simplified relations. This strategy ensures that critical flow pathways—governing radionuclide migration or heat extraction—are governed by high-fidelity physics, without sacrificing computational feasibility.

5.2 Nuclear Waste Geological Disposal: Ensuring Long-Term Safety through THMC Coupling

The deep geological disposal of HLW is one of the most challenging long-term engineering endeavors humanity has undertaken. Its design objective is to safely isolate radionuclides from the biosphere over geological timescales of tens of thousands to hundreds of thousands of years [80]. Achieving this ambitious goal critically depends on a profound understanding and accurate modeling of coupled THMC processes [81]. A geological repository is not a static “container” but a dynamic system composed of engineered and natural barriers that continuously evolve under the influence of a thermal pulse. Its long-term safety is, in essence, the outcome of the spatiotemporal interplay of THMC coupling mechanisms.

5.2.1 THM Coupling Mechanisms in the Multi-Barrier System

Modern nuclear waste repositories employ a “multi-barrier” design philosophy, centered on three interdependent subsystems: the waste canister, the bentonite buffer layer, and the surrounding host rock. The functionality of each barrier is critically governed and finely regulated by coupled THM processes [82].

• Waste Canister: As the first line of defense, the waste canister’s primary function is to encapsulate the waste and provide initial mechanical protection. However, its more critical role is that of a powerful heat source. The thermal energy released by radioactive decay (the thermal pulse) serves as the primary driver for the entire system’s evolution. Numerical simulations show that the canister surface temperature can initially reach 80–100°C, with thermal effects extending several meters into the surrounding barriers. This thermal pulse is not merely a “burden” to be mitigated; rather, it acts as the “key” that activates the functional responses of the subsequent barriers.
• Bentonite Buffer Layer: Positioned between the waste canister and the host rock, the buffer layer is composed of highly compacted sodium bentonite. Its safety functions are entirely governed by THM coupling [83]: T → H: The thermal pulse drives moisture migration from the high-temperature zone (near the canister) toward the cooler host rock, accelerating the hydration of bentonite. H → M: Water uptake induces hydration swelling of bentonite particles, generating swelling pressures as high as several megapascals. This pressure provides uniform mechanical support against the excavation-damaged walls of the repository, effectively suppressing relaxation deformation of the host rock and ensuring the canister remains stable during seismic events. T → M: Thermal expansion, combined with hydration-induced swelling, further improves the contact between the buffer and the host rock, eliminating voids and forming a low-permeability seal. Simulation results show that under THM coupling, the porosity of bentonite can decrease from an initial value of 0.3 to below 0.1, and its permeability can reach as low as 10−20 m2—making it a near-perfect barrier against radionuclide migration [84,85].
• Host Rock: As the ultimate natural barrier, the host rock’s role is to provide a long-term stable geological environment. However, its effectiveness does not stem from the intact rock matrix itself, but rather from the integrity and evolution of its internal fracture network [86]. Under THM coupling, the host rock undergoes complex physical changes: T → M: Thermal stresses may cause pre-existing fractures to open or close, and in some cases, even initiate new microfractures. H → M: Groundwater flow through fractures alters local effective stress, thereby influencing the mechanical stability of those fractures. M → H: Any changes in fracture geometry—induced by stress redistribution or thermal stresses—directly and nonlinearly modify fracture permeability, thereby affecting the hydraulic connectivity of the entire system.

5.2.2 Fracture Network in the Host Rock: The Ultimate Controller of Radionuclide Migration

Although engineered barriers provide robust protection in the early stages, over timescales of tens of thousands of years, the potential release pathways for radionuclides are almost entirely governed by the fracture network in the host rock. In low-permeability crystalline rocks—such as granite or tuff—the contribution of the rock matrix to fluid transport is negligible; fractures serve as the sole conduits for fluid flow and radionuclide migration. Consequently, the long-term safety assessment of a geological repository fundamentally hinges on predicting the evolution of the fracture network under coupled THMC processes.

This prediction faces immense challenges:

• (1)Path dependence: The evolution of fracture permeability exhibits strong dependence on stress history. TM perturbations during repository operation may permanently alter the hydraulic properties of fractures [87].
• (2)Long-term chemical effects: Over geological timescales, water–rock chemical reactions (the C field) will slowly but continuously alter fracture surfaces. Mineral dissolution may increase aperture, while precipitation can clog flow pathways. This HMC coupling process ultimately governs the long-term rate of radionuclide migration [88].
• (3)Network topological complexity: Radionuclides do not migrate via uniform diffusion; instead, they travel along “preferential pathways” dictated by fracture connectivity. Accurately capturing this “channeling effect” requires models to go beyond the equivalent continuum assumption and explicitly represent the three-dimensional topological structure of the fracture network [89].

Nature itself provides compelling evidence: the Oklo natural nuclear reactors in Gabon underwent self-sustaining fission nearly 2 billion years ago, yet the resulting transuranic elements—such as plutonium—remain largely immobilized at the original sites to this day [90], with no significant migration observed. This “natural experiment” demonstrates that, under favorable geological conditions, fracture systems in rock can indeed retain radionuclides over billion-year timescales. Our engineering objective is to actively create and sustain such a “natural” state of safety through precise THMC-coupled design.

In summary, the success of geological nuclear waste disposal does not lie in constructing an indestructible “iron coffin”, but in designing a “living system” that actively leverages physical laws—particularly THMC coupling—to self-reinforce and self-seal over time. The deliberate utilization of bentonite swelling mechanisms and the accurate prediction of long-term evolution in the host rock’s fracture network constitute the two foundational pillars of system safety. Table 5 clearly illustrates the THM coupling processes within each barrier, serving as a key to understanding their synergistic function. Any uncertainty in modeling ultimately translates into a loss of safety margin; therefore, developing high-fidelity, multiphysics coupled models capable of resolving fracture-scale physical processes remains the central mission of research in this field.

5.3 EGS: Optimization of Heat Extraction and Risk Management

EGS represent a cutting-edge engineering endeavor in which humans actively modify deep subsurface geological environments to harness clean energy. The core objective of EGS is to inject high-pressure fluid into low-permeability hot dry rock to artificially activate or enhance natural fracture networks, thereby creating a sustainable underground heat exchanger [91]. The success of this process hinges entirely on precise control of the coupled TMH processes. EGS is not merely an energy project; it is a complex geophysical experiment operating under extreme conditions—high temperature, high pressure, and intense in-situ stress. Its design must achieve a delicate balance between maximizing heat extraction efficiency and minimizing engineering risks.

5.3.1 The TMH Coupling Chain Driven by Cold Fluid Injection

The core physical process of EGS operation begins with the high-pressure injection of cold fluid (typically water at around 50°C). This initial hydraulic action (H) immediately triggers a tightly coupled sequence of physical responses:

• (1)TM coupling (T → M): Upon contact between the cold fluid and the hot rock mass (e.g., 200°C granite), a sharp thermal gradient develops. Localized cooling causes thermal contraction of the rock, generating significant thermal stresses that concentrate at pre-existing natural fractures or weak planes [92].
• (2)Mechanical-hydraulic coupling (M → H): When thermal stresses exceed the shear strength of fracture surfaces, shear slip—commonly referred to as “hydraulic shearing”—is induced. This process causes macroscopic shear dilation of the fracture walls, significantly increasing aperture and thereby enhancing permeability by several orders of magnitude. This mechanism lies at the heart of EGS reservoir stimulation [93].
• (3)Hydraulic-thermal coupling (H → T): The enhanced permeability dramatically improves fluid circulation within the reservoir, allowing cold fluid to flow more efficiently through the hot rock matrix, absorb heat, and transport it to the production well—thus completing the closed-loop heat extraction cycle.

This positive feedback loop—“injection → cooling → stress → shear dilation → permeability enhancement → heat extraction”—forms the physical foundation for economically viable EGS operation. However, this loop inherently embodies a fundamental contradiction: the co-dependence of permeability enhancement and thermal drawdown. While increased permeability boosts fluid flow rates, it simultaneously accelerates heat extraction and rock cooling, thereby introducing long-term stability risks to the system.

5.3.2 Quantitative Impact of Key Parameters on System Performance

The ultimate economic viability of EGS hinges on their long-term thermal power output and the Levelized Cost of Electricity (LCOE). Comprehensive numerical simulations and field trials have identified several decisive design parameters that dictate system performance. First, the injection–production pressure difference serves as the primary operational lever for flow control [94]. Second, the rock thermal expansion coefficient is pivotal, as it governs the magnitude of thermally induced stresses within the reservoir [95]. Lastly, the initial fracture permeability represents the reservoir’s “inherent quality”, acting as a critical determinant of long-term efficiency.

These quantitative relationships provide a direct basis for optimizing EGS design. Engineers can perform sensitivity analyses to tailor operational parameters—such as well spacing, injection rate, and fluid viscosity—to specific geological conditions, thereby maximizing energy output. For instance, using supercritical carbon dioxide (SCCO2) as the working fluid can further reduce flow resistance and enhance heat extraction efficiency, owing to its high compressibility and low viscosity [96].

5.3.3 “Thermal Breakthrough” and System Lifetime Management

“Thermal breakthrough” is the primary bottleneck limiting the operational lifespan of an EGS. It refers to the irreversible decline in the outlet fluid temperature at the production well caused by excessive localized cooling of the reservoir. Once thermal breakthrough occurs, the system’s thermal efficiency drops sharply, rendering it economically unviable [97]. The root cause lies in the over-activation of the aforementioned TMH coupling chain—permeability increases too rapidly, allowing the cold front to advance swiftly toward the production well.

To delay thermal breakthrough, several engineering strategies are employed: Optimized well pattern design: Increasing the distance between injection and production wells extends the fluid’s residence time within the hot reservoir. Controlled injection strategies: Pulsed or staged injection is used to avoid sustained high pressure that could cause excessive fracture dilation.

Development of the “regenerative EGS” concept [98]: Integrating the system with the power grid allows production to be paused during periods of low electricity demand, enabling the reservoir to naturally reheat and effectively function as a “seasonal thermal storage”, thereby extending the system’s overall lifespan.

5.3.4 Induced Seismicity Risk and Engineering Trade-Offs

The most controversial risk associated with EGS is its potential to induce felt seismicity. When injected fluid triggers shear slip along fractures, and if the slip is sufficiently large or occurs on seismically active faults, it can release accumulated strain energy, generating seismic events. This risk stands in sharp contradiction to thermal energy extraction efficiency: maximizing permeability enhancement often implies maximizing seismic risk.

Therefore, modern EGS projects must implement a strict “traffic light” control system: green light: microseismic activity remains within safe thresholds, allowing injection to continue or increase. yellow light: microseismic activity intensifies, requiring a reduction in injection rate or pressure. red light: felt seismic events or precursors to larger earthquakes are detected, triggering immediate cessation of injection [99].

The essence of a successful EGS project lies in performing a “minimally invasive intervention” on the reservoir through high-fidelity THM coupled models—precisely activating target fractures while avoiding or stabilizing potentially seismogenic faults. This demands models capable not only of predicting permeability evolution but also of assessing the mechanical stability of shear slip and the potential for seismic moment release [100]. Ultimately, the sustainable development of EGS hinges on our ability to translate the predictive power of physical models into real-time, dynamic risk management strategies, striking an acceptable balance between energy gain and societal risk.

5.4 Tunneling and Underground Engineering: Seepage Stability and Multiphysics Coupling Management

The construction of deep tunnels and large underground caverns represents a core engineering endeavor in humanity’s expansion into the Earth’s subsurface. Its success hinges directly on the accurate prediction and proactive management of the complex coupling among the stress field, seepage field, and temperature field in the surrounding rock mass under excavation-induced disturbance. Unlike “active modification” projects such as nuclear waste disposal or EGS, the primary objective of tunneling is “passive defense”—maintaining long-term stability of the surrounding rock and structural serviceability under the dual threats of high in-situ stress and high hydraulic pressure [101]. Achieving this goal requires engineers to deeply understand and manage the multiphysics coupling chain triggered by excavation unloading, as well as the progressive erosion effect of seepage forces on support structures [102].

5.4.1 Coupling Effect of Excavation Unloading and Pore Water Pressure Redistribution

The essence of tunnel excavation is the creation of a “void” within the in-situ stress field, triggering a dramatic redistribution of stress in the surrounding rock. This M immediately engages in strong coupling with the H:

• Stress-path-dependent permeability change: Excavation unloading creates a low-stress zone (relaxation zone) around the tunnel, where fractures open due to stress release, significantly increasing permeability [103]. This forms “highway channels” that facilitate groundwater inflow toward the excavation.
• Hydraulic gradient-driven seepage: The natural hydraulic head difference in the surrounding rock, after excavation, establishes a strong hydraulic gradient directed toward the tunnel. Guided by the newly formed high-permeability pathways, groundwater flows into the tunnel at rates far exceeding pre-excavation levels, resulting in water inrush or dripping seepage.
• Pore pressure dissipation and sharp increase in effective stress: As groundwater drains toward the tunnel, pore water pressure (p) in the surrounding rock drops sharply. According to the principle of effective stress (σ′=σ−αp), this causes a substantial increase in the effective stress (σ′) acting on the rock skeleton. In high in-situ stress regions, this sudden rise in effective stress may exceed the rock mass strength, triggering new fracturing or further propagation of existing fractures, thereby initiating a vicious cycle: “seepage → stress concentration → fracturing → increased permeability → intensified seepage”.

This coupled process constitutes the fundamental physical mechanism behind water and mud inrush hazards in tunneling [104]. Notably, this coupling chain is distinctly manifested in geological environments containing karst collapse columns (KCPs), where the advance of the working face triggers critical spatial-temporal changes in the surrounding rock mass. As illustrated in Fig. 8, the evolution of the damage zone and redistribution of seepage vectors before and after KCP penetration clearly demonstrate how mining-induced fractures connect with pre-existing discrete cracks in KCPs, forming a continuous water-conducting network; concurrently, stress concentration redirects seepage paths toward the working face, further amplifying the synergy between seepage-driven erosion and mechanical rock failure. Successful engineering design must therefore proactively interrupt this intensified coupling process through targeted measures such as advanced geological prediction (to identify KCPs and fracture networks), grouting curtains (to block preferential flow paths), or dewatering and pressure relief (to reduce hydraulic gradients). High-pressure water injection is one commonly used technique, employed to keep both the hydraulic gradient and stress concentration within safe thresholds [105].

Figure 8: Damage zone and seepage vectors evolution when the working face is advanced before and after the penetrated KCP. Adapted with permission from Reference [106]. Copyright © 2022, Springer.

5.4.2 Long-Term Effects of Seepage Forces on Support Structures and Rock Mass Degradation

Seepage is not only an immediate mechanical threat but also a long-term, cumulative “chronic condition”. Its detrimental effects on engineering structures manifest primarily at two levels:

• Direct mechanical action—seepage force: Flowing groundwater exerts a “seepage force” on fracture walls and support structures, acting in the direction of flow. Under high hydraulic head gradients, this seepage force can reach the megapascal scale—sufficient to mobilize loose rock blocks, erode infill materials, and even exert sustained thrust on shotcrete linings or steel arches, leading to deformation, cracking, and potential instability of the support system [107].
• Indirect chemo-mechanical action—rock mass degradation: Long-term seepage accelerates the chemical weathering and mechanical deterioration of the surrounding rock. Groundwater dissolves cementing minerals (e.g., calcite, clay minerals), reducing rock strength and promoting fracture propagation. Simultaneously, the lubricating effect of water lowers the friction coefficient along fracture surfaces, weakening their shear resistance [108]. Experimental tests further confirm this mechanism: as illustrated in Fig. 9, the friction coefficient of basalt fractures shows a significant decreasing trend with the increase of pore water pressure, and this weakening effect is more pronounced in gouge-bearing fractures. Specifically, Fig. 9b (bare basalt fractures) and Fig. 9c (gouge-bearing basalt fractures) clearly demonstrate that pore water not only lubricates the fracture surface directly but also dissolves cementing minerals in the gouge, further reducing the shear resistance of the rock mass. This “water–rock–mechanics” degradation is slow and insidious, often remaining latent for years after construction before manifesting visibly—yet it can fundamentally undermine the long-term structural safety. Therefore, lifecycle management of tunnels must include continuous monitoring and assessment of seepage-induced corrosion effects.

Figure 9: Variations in the friction coefficient with effective normal stress under various (a) Gouge layer thicknesses; (b) Pore pressures (bare basalt fractures); and (c) Pore pressures (gouge-bearing basalt fractures). The μ  values are taken at a shear displacement of 3 mm and applied to a sliding velocity of 6 μm/s. Adapted with permission from Reference [109]. Copyright © 2023, Elsevier.

5.4.3 Unique THM Challenges in Freeze Construction in Cold Regions

In water-rich weak strata or high-hydraulic-pressure zones, AGF is an effective temporary support technique. The method works by artificially refrigerating the ground to form a closed “frozen soil curtain” around the tunnel, leveraging the high strength and low permeability of frozen soil to isolate groundwater [110]. However, this approach transforms the original HM coupling problem into a far more complex THM three-field coupling challenge.

• The double-edged nature of frost heave pressure: Water expands by approximately 9% upon freezing, generating substantial frost heave pressure (up to several megapascals) in confined spaces [111]. This pressure can be harnessed beneficially to compact and densify the surrounding soil, thereby enhancing its mechanical strength. However, if not properly controlled, frost heave pressure can act adversely on support structures, causing lining cracks or deformation.
• Strong nonlinearity of the phase-change process: The advancement of the freezing front involves a highly nonlinear phase transition, accompanied by latent heat release, moisture migration, and ice lens formation. Minor changes in the temperature field (T) can trigger dramatic responses in the H and M fields. For instance, unfrozen water migrates toward the freezing front under thermal gradients, intensifying localized frost heave; meanwhile, ice formation directly alters the soil’s permeability and strength [112].
• Anisotropy and asymmetric deformation: Due to variations in geological conditions, refrigeration pipe layout, or groundwater flow, the temperature field within the frozen curtain is often non-uniform, resulting in pronounced anisotropy in both frozen soil strength and frost heave pressure. This leads to asymmetric stress and deformation distributions around the tunnel, imposing higher demands on support structure design. Numerical simulations show that neglecting this anisotropy can result in severe misestimation of structural loads.
• Cumulative damage from freeze-thaw cycles: During later construction stages or operational periods, thawing of the frozen curtain can induce thaw settlement. Repeated freeze-thaw cycles inflict cumulative damage on both the rock/soil mass and support structures, progressively degrading their strength and stiffness and continuously reducing the safety factor [113].

Therefore, the success of freeze construction critically depends on high-fidelity THM coupled models. Such models must be capable of simulating the entire process—including ice-water phase change, moisture migration, frost heave pressure generation, and stress redistribution—to provide a scientific basis for optimizing refrigeration schemes, designing support structures, and managing thaw-induced settlement risks [110]. This represents not only a technical challenge but also the ultimate test of multiphysics coupling theory’s predictive capability under extreme engineering conditions.

6 The Transformative Role of AI and ML

6.1 AI/ML-Driven New Paradigm

The rapid advancement of AI and ML is injecting a transformative force into the study of fluid flow and multiphysics coupling in fractured rock masses, accelerating a paradigm shift from traditional, physics-based deterministic modeling toward a new data-driven, pattern-recognition, and intelligent prediction framework [114]. At its core, this shift redefines the rock mass not as a set of partial differential equations to be solved, but as a complex data system rich in hidden information. By uncovering latent correlations and nonlinear patterns within data, AI/ML enables efficient and highly accurate prediction of system behavior [115].

The strength of this new paradigm lies in its powerful generalization capability and its ability to handle high-dimensional, heterogeneous data. Faced with the inherent strong heterogeneity, multiscale nature, and data sparsity of fractured rock masses, traditional physics-based models often rely on numerous simplifying assumptions. In contrast, AI/ML can directly learn complex input–output mappings from raw data—such as borehole images, geophysical logs, microseismic monitoring records, and digital rock cores—bypassing the need for explicit physical modeling and offering an “end-to-end” solution [116].

6.1.1 Typical Applications in Rock Engineering

• (1)Intelligent Rock Mass Classification and Prediction of Mechanical Parameters

Traditional rock mass classification systems (e.g., RMR, Q-system) heavily rely on expert judgment, making them subjective and inefficient [117]. AI/ML algorithms—such as Support Vector Machines (SVM), Random Forests (RF), and Artificial Neural Networks (ANN)—have been successfully applied to automate rock mass classification and enable rapid prediction of key mechanical parameters, such as uniaxial compressive strength (UCS) and elastic modulus [118]. For instance, by inputting image features from drill core samples or geophysical logging curves, these models can output rock mass mechanical properties within seconds, achieving accuracy comparable to laboratory testing and providing real-time, low-cost data support for engineering design.

• (2)Intelligent Permeability Prediction and Digital Rock Core Analysis

Permeability, as a core parameter in seepage analysis, is costly to measure directly and often lacks representativeness. AI/ML techniques—particularly three-dimensional convolutional neural networks (3D CNNs)—can now directly extract microstructural pore features from high-resolution digital rock images and predict macroscopic permeability within seconds, achieving speedups of several orders of magnitude compared to traditional methods such as the LBM or direct numerical simulation [119]. This breakthrough enables large-scale, high-precision reservoir property evaluation.

• (3)Intelligent Characterization and 3D Reconstruction of Fracture Networks

Accurate identification and geometric parameter extraction of fractures are fundamental to DFN modeling, yet traditional image segmentation methods often suffer from errors due to noise and complex backgrounds. Intelligent segmentation algorithms based on RF and deep learning can now identify fracture boundaries with high precision from complex 2D/3D images—such as CT scans or outcrop photographs—effectively overcoming the subjectivity of manual interpretation and the tendency of conventional algorithms to overestimate fracture aperture [120]. Moreover, AI models can integrate multisource data (e.g., seismic attributes, borehole imaging) to construct reservoir-scale 3D fracture intensity models, enabling a leap from “point” observations to “volumetric” characterization.

• (4)Hydraulic Fracturing and Production Forecasting

In oil/gas and geothermal engineering, predicting fracture geometry and production performance after hydraulic fracturing remains a central challenge. Purely physics-based models are limited by high computational costs and difficulty in capturing complex nonlinear behaviors. AI/ML models, by learning from vast historical fracturing datasets—including geological parameters, pumping schedules, microseismic events, and production data—can now predict fracture propagation patterns and ultimate productivity. This capability provides powerful real-time decision support for optimizing injection strategies in the field [121]. To further improve the generalization ability and computational efficiency of AI/ML models, parameter sensitivity analysis is critical for identifying key input factors and simplifying model structures. Fig. 10 presents the sensitivity analysis of two typical regression models (KRR and GPR), which quantifies the impact of core parameters (e.g., Young’s modulus, injection rate) on the prediction accuracy of fracture geometry. The results show that injection rate and closure stress are the most sensitive parameters for both models, providing a theoretical basis for feature selection in model optimization.

Figure 10: The sensitivity analysis of the KRR and GPR. Adapted with permission from Reference [122]. Copyright © 2025, Elsevier.

6.1.2 Maturity Assessment of AI/ML Methods in Practice

It is critical to distinguish between AI/ML approaches that have achieved engineering validation and those that remain emerging or experimental. Based on field deployment, reproducibility, and industry adoption:

• Widely validated methods include RF, SVM, and shallow ANNs. These have been successfully applied across multiple projects for rock mass classification, permeability estimation from logs, and production forecasting in conventional reservoirs [123]. Their robustness, interpretability (to some extent), and low data requirements make them suitable for routine engineering use.
• Emerging but promising methods include 3D Convolutional CNNs for digital rock physics, graph-based models for fracture network reconstruction, and physics-informed surrogates [124]. While they demonstrate high accuracy in research settings, their reliance on large, high-quality 3D datasets and limited field validation restrict current use to pilot studies or academic prototypes.
• Frontier research directions—such as Graph Neural Networks (GNNs) for dynamic DFN evolution or transformer-based models for spatiotemporal forecasting—remain largely conceptual [125], with few real-world demonstrations in deep subsurface environments due to data scarcity and computational complexity.

This tiered perspective ensures that practitioners can select appropriate tools based on project risk tolerance, data availability, and regulatory requirements—especially in safety-critical contexts like nuclear waste disposal or deep tunneling.

In summary, the AI/ML-driven paradigm is not intended to replace physics-based models but rather to serve as a powerful “accelerator” and “amplifier” for them. By automating and intelligently processing data and uncovering hidden patterns, AI/ML dramatically enhances the efficiency and accuracy of our understanding of complex geological systems. The applications listed in Table 6 offer innovative technical pathways to bridge the longstanding “data-model” gap that has challenged traditional approaches.

Nevertheless, it should be acknowledged that high-quality labeled data remain extremely scarce in deep underground engineering environments—often limited to sparse borehole measurements or indirect geophysical proxies. While few-shot learning and transfer learning offer promising avenues to mitigate data scarcity, their applicability hinges critically on the physical similarity between source and target domains. For instance, transferring a model trained on shallow sedimentary rock masses to deep crystalline fractured reservoirs may yield predictions with high confidence but poor physical fidelity, especially when key THMC coupling mechanisms differ significantly. To address this risk, future efforts should embed fundamental conservation laws as hard constraints within the transfer process and combine AI with active learning strategies to prioritize acquisition of informative field data under critical loading or injection conditions. Such safeguards are essential to ensure robustness and reliability when deploying AI in safety-critical subsurface applications.

6.2 Frontiers in Physics-Informed and Hybrid Modeling

Although purely data-driven AI/ML models demonstrate remarkable capabilities in prediction speed and pattern recognition, their “black-box” nature and lack of physical consistency remain fundamental barriers to widespread adoption in high-risk engineering applications. To bridge the gap between data-driven and physics-driven paradigms, two cutting-edge approaches have emerged: PINNs and Hybrid Modeling. These represent a profound shift in AI—from empirical fitting toward physics-embedded learning [126].

6.2.1 PINNs: Embedding Physical Laws into Neural Networks

The core innovation of PINNs lies in treating physical laws not as external constraints, but as intrinsic, inviolable rules directly embedded into the neural network’s loss function. Specifically, the loss function of a PINN consists of two components: (1) a data-fitting term (e.g., mean squared error between observations and network predictions), and (2) a physics residual term (i.e., the residual obtained when the network’s output is substituted into the governing partial differential equations, or (PDEs) partial differential equations). By minimizing the total loss, PINNs simultaneously fit the available data and enforce the solution to satisfy fundamental physical conservation laws—such as mass, momentum, and energy conservation—as well as constitutive relationships.

The core advantages of PINNs lie in:

• Guaranteed physical consistency: Model predictions inherently satisfy physical laws, eliminating non-physical solutions that may arise from purely data-driven approaches, thereby significantly enhancing result credibility and extrapolation capability.
• Mesh-free solution: PINNs do not require discretization of the solution domain into a mesh, thereby fundamentally avoiding the mesh distortion and remeshing challenges that plague traditional numerical methods (e.g., FEM) when handling complex geometries—such as rough fractures or moving boundaries—making them particularly well-suited for strongly discontinuous media like fractured rock masses.
• Data efficiency: By leveraging physical laws as strong prior knowledge, PINNs can produce physically plausible solutions even when observational data are sparse or noisy, significantly reducing reliance on large volumes of high-quality training data.
• High-dimensional and inverse problem solving: PINNs demonstrate unique advantages in handling high-dimensional parameter spaces and solving inverse problems—such as inferring permeability or stress fields from observational data—as their framework naturally supports the simultaneous optimization of state variables and unknown parameters.

6.2.2 Challenges in Enforcing Discontinuities at Fracture Intersections

Despite their theoretical suitability for fractured media, the practical implementation of PINNs faces significant challenges when handling discontinuous flow fields at fracture intersections, where hydraulic head gradients and fluxes can change abruptly due to mismatched transmissivity, aperture heterogeneity, or contact resistance. Standard PINN formulations assume sufficient smoothness of the solution to compute PDE residuals via automatic differentiation—a requirement violated at geometric singularities such as triple junctions or crossing fractures.

To embed such discontinuous interface conditions into the loss function, a common strategy is to explicitly enforce mass conservation and head compatibility at intersection points as soft constraints. Consider two fractures T 1 and T 2 intersecting at point x i , with hydraulic heads h 1 and h 2 governed by 1D flow equations along each fracture. Assuming head continuity (in the absence of skin effects) and flux balance, the interface conditions read: h1(xi)=h2(xi),T1∇T1h1(xi)+T2∇T2h2(xi)=0(11) where T k is the transmissivity of fracture k , and Δ T k denotes the tangential gradient along the fracture direction.

These conditions are incorporated into the total loss as penalty terms: Linterface=λh|h1(xi)−h2(xi)|2+λq|T1∇T1h1(xi)+T2∇T2h2(xi)|2(12) where λ h , λ q are tunable weights. While effective, this approach introduces sensitivity to weight selection and requires dense sampling near intersections—highlighting a key implementation difficulty in real-world DFNs with thousands of junctions. Emerging alternatives, such as graph-based PINNs or domain decomposition, may offer more scalable solutions for complex networks.

These capabilities have been concretely demonstrated in a recent study by Liu et al. [127], who developed an Enriched Physics-Informed Neural Network (E-PINN) specifically for simulating two-phase flow in heterogeneous and fractured porous media—a problem long considered challenging for standard PINNs due to sharp property contrasts and interface discontinuities. To rigorously enforce flux continuity across fracture-matrix interfaces, the authors integrated the Embedded Discrete Fracture Model (EDFM) and replaced automatic differentiation with a finite volume method (FVM) to construct the physics-informed loss, ensuring accurate treatment of conservation laws at material boundaries. Their architecture further incorporates adjacency-aware anchoring, adaptive activation functions, and hard enforcement of initial/boundary conditions, significantly improving training stability and solution accuracy. Validated on 2D and 3D immiscible displacement problems in complex fractured reservoirs, this work provides a robust, mesh-free framework that directly addresses the geometric and physical complexities inherent in fractured rock masses—even under sparse data conditions.

To enhance the training stability and efficiency of PINNs, cutting-edge research has developed several advanced techniques, such as Exponential Moving Average (EMA) for weight stabilization, Residual-based Adaptive Refinement with Distribution (RAR-D) for intelligent sampling, and Bayesian optimization to automatically balance the weights between the physics residual and data-fitting terms [128]. In the context of fracture evolution in rock mechanics, strategies like transfer learning have also been introduced to avoid retraining the network from scratch at each loading step, thereby significantly improving computational efficiency.

6.2.3 Hybrid Models: Synergistic Co-Evolution of Data-Driven and Physics-Driven Approaches

If PINNs embed physical laws within the neural network, hybrid models achieve strategic synergy between data-driven and physics-driven approaches at the system architecture level. The core idea is to use a physics-based model to generate a reliable, interpretable baseline prediction, and then employ a data-driven model (e.g., a neural network) to learn the residual—or the mapping—between the physics-based prediction and real-world observations. This compensates for biases introduced by simplifying assumptions or parameter uncertainties inherent in the physical model.

Typical hybrid frameworks include:

• Low-fidelity physics-based model + high-fidelity data-driven corrector: For example, a computationally efficient ECM is used for preliminary simulation, and a neural network is then trained on the ECM’s inputs and outputs together with field monitoring data to learn the nonlinear patterns of the ECM’s prediction errors. The final prediction is given by “ECM output + neural network correction”, significantly improving accuracy while preserving physical interpretability [129].
• Physics-guided AI surrogates: Outputs from physics-based models (e.g., stress fields, temperature fields) are used as input features for AI models, guiding them to learn more complex, implicit relationships that the physics-based model cannot directly capture—such as the correlation between microseismic events and permeability evolution.
• Multiscale hybrid simulation: Physics-based models are employed at the macroscale, while data-driven models (e.g., CNNs for predicting local permeability) are used at micro- or mesoscales. Cross-scale coupling enables an optimal balance between computational efficiency and local accuracy.

The greatest value of hybrid models lies in their unification of robustness and interpretability. They inherit the strengths of physics-based models in mechanistic understanding and reliable long-term extrapolation, while simultaneously harnessing the power of AI models to handle high-dimensional, nonlinear, and heterogeneous data. This synergistic “1 + 1 > 2” effect makes hybrid modeling the most promising paradigm for tackling the extreme complexity of multiphysics coupling in fractured rock masses.

In summary, PINNs and hybrid models represent the most advanced form of AI for Science in geotechnical engineering. They do not seek to overturn traditional methods but rather to revolutionize and augment their capabilities. In the future, as algorithms continue to mature and computational power grows, physics-informed and hybrid modeling will inevitably transition from cutting-edge research to the core of engineering practice, laying a solid foundation for the next generation of intelligent rock mechanics models—ones that are interpretable, trustworthy, and capable of reliable extrapolation [130].

6.3 Comparative Analysis of AI/ML and Traditional Methods

The rise of AI and ML is not aimed at overturning or replacing traditional physics-based numerical simulation methods, but rather at providing a powerful complementary and enhancing tool. In the study of fluid flow and multiphysics coupling in fractured rock masses, AI/ML and traditional methods (e.g., FEM, DEM, DFN) each possess distinct strengths and limitations. Their true value lies not in competition, but in strategic integration—leveraging their complementary advantages to jointly tackle complex engineering challenges that neither approach can resolve alone [131]. Table 7 presents a systematic comparison between the two across multiple key dimensions.

As shown in the table above, AI/ML and traditional methods are fundamentally complementary rather than mutually exclusive. The “speed” and “data-driven” strengths of AI/ML precisely address the “computational slowness” and “parameter dependency” limitations of traditional methods, while the “physical consistency” and “interpretability” of traditional approaches provide essential constraints and validation benchmarks for the “black-box” nature of AI/ML.

Future research and engineering practice will increasingly adopt the following synergistic strategies:

• (1)AI as an “accelerator” for traditional models: Train AI models as surrogate models for computationally expensive physics-based simulations. For example, a neural network can replace time-consuming DFN or THMC coupled simulations to enable second-level parameter sensitivity analysis or real-time optimization, while the physics-based model is reserved for validation and calibration in critical scenarios.
• (2)Traditional models provide “physical constraints” for AI: Embed physical laws into AI training (e.g., via PINNs) or use outputs from physics-based models as input features for AI, ensuring the physical plausibility of AI predictions and enhancing their generalization capability in data-sparse regions.
• (3)Hybrid modeling frameworks: Construct an integrated architecture combining a “physics-based model + AI corrector”. The physics-based model provides a baseline prediction, while the AI component learns the residual between this baseline and real-world observations, thereby significantly improving predictive accuracy while preserving physical interpretability.
• (4)AI-driven data assimilation and parameter inversion: Leverage AI’s powerful pattern recognition capabilities to efficiently invert key subsurface parameters—such as permeability and stress fields—from sparse, noisy field monitoring data (e.g., microseismic events, temperature, pressure), thereby providing more accurate inputs for traditional physics-based models.

In summary, the integration of AI/ML with traditional numerical methods represents an evolution in fractured rock mass research—from reliance on a “single tool” to the deployment of an “intelligent toolchain”. Table 8 clearly illustrates the respective strengths and weaknesses of AI/ML and traditional approaches, offering an intuitive basis for understanding their complementary relationship. Future breakthroughs will no longer hinge on the extreme advancement of any single method, but on how skillfully we weave this synergistic “physics + data” network to achieve an unprecedented balance among computational efficiency, predictive accuracy, physical robustness, and engineering practicality. This is not merely a choice of technical pathway, but a profound transformation in research paradigm [132].

6.4 Embedding AI into Traditional Simulation Workflows: Practical Integration Pathways for Engineering Applications

While the complementary nature of AI/ML and traditional numerical methods has been widely acknowledged, a critical gap remains between theoretical synergy and practical implementation in real-world engineering projects. To bridge this gap, it is essential to move beyond conceptual comparisons and articulate concrete, actionable strategies for mbedding AI modules into established simulation workflows—such as those based on FEM, DFN, or dual-porosity models—without disrupting existing engineering practices.

Three pragmatic integration pathways have emerged as particularly viable for field deployment:

• (1)AI as a Plug-in Surrogate for Bottleneck Components

Identify computationally intensive subroutines within a traditional workflow (e.g., fracture network flow solver, THMC coupling loop, or parameter sensitivity analysis) and replace them with pre-trained AI surrogates. For instance, a 3D CNN trained on LBM simulations of digital rock cores can instantly predict permeability tensors [133], which are then fed into a reservoir-scale DFN model. This “swap-and-run” approach requires minimal modification to legacy code and delivers orders-of-magnitude speedup while preserving the overall workflow structure.

• (2)AI-Augmented Physics-Based Simulators via Residual Learning

In this hybrid architecture, the physics-based simulator runs as usual to generate a baseline prediction (e.g., pore pressure evolution during tunnel excavation). A lightweight neural network—trained on historical monitoring data (e.g., piezometer readings, extensometer displacements)—learns the residual error between simulation and reality [134]. The final output is the sum of the physics-based prediction and the AI correction. This strategy compensates for model simplifications (e.g., homogenized fracture properties) without sacrificing interpretability, making it suitable for safety-critical applications like dam foundation seepage analysis.

• (3)AI-Driven Real-Time Data Assimilation Loops

Integrate streaming field data (e.g., from fiber-optic DTS/DAS, microseismic arrays, or IoT sensors) into the simulation cycle via AI-enabled inversion. For example, a BNN can continuously update fracture aperture distributions in a DFN model using real-time pressure transients [135], effectively creating a “digital twin” of the subsurface system. This closed-loop framework transforms static simulations into dynamic, adaptive decision-support tools—critical for operations such as EGS reservoir management or CO2 injection monitoring.

Key Enablers and Implementation Considerations:

• Modular software design: Adoption hinges on simulation platforms (e.g., OpenGeoSys, TOUGH–FLAC3D, COMSOL) supporting external AI calls via APIs or co-simulation interfaces.
• Uncertainty quantification: Techniques like Monte Carlo dropout or ensemble methods should be integrated to propagate AI prediction uncertainty into final engineering decisions.
• Data curation pipelines: Field engineers need standardized protocols to convert raw sensor logs into AI-ready features (e.g., time-series embeddings, spatial statistics).

These integration pathways demonstrate that AI does not require a “rip-and-replace” overhaul of existing workflows. Instead, it can be incrementally embedded as intelligent, task-specific modules—enhancing accuracy, speed, and adaptability while respecting the institutional knowledge and validation standards embedded in traditional geomechanical practice.

7 Core Challenges and Future Research Directions

The profound scale-dependent discrepancies between laboratory measurements and field-scale behavior constitute a primary motivation for the development of advanced multiscale modeling frameworks. As shown, key hydraulic and mechanical parameters (e.g., permeability, reaction rates) measured on centimeter-scale rock cores often deviate by orders of magnitude from their in-situ counterparts due to the irreducible complexity of natural fracture networks [136], geological history, and boundary conditions at engineering scales. This “scale gap” renders single-scale models inherently unreliable for field prediction [137]. Consequently, the multiscale coupled modeling strategies discussed in this section are not merely theoretical refinements; they represent a necessary paradigm shift to bridge the chasm between controlled laboratory physics and real-world fractured rock systems.

7.1 Multiscale Coupled Modeling

One of the most fundamental challenges in the study of fluid flow and multiphysics coupling in fractured rock masses is bridging the vast gap between microscopic pore/fracture scales and macroscopic engineering/geological scales [138]. The macroscopic behavior of rock masses is not a simple summation of its components but an emergent property arising from the collective evolution of millions of microscopic physical processes—such as non-Darcy flow, fluid–rock reactions, and fracture wall contact mechanics—under complex geometric constraints. While current mainstream models have matured significantly at individual scales, the core bottleneck limiting predictive accuracy and generalizability remains the lack of a theoretical framework that seamlessly links physical mechanisms across scales and enables a consistent transfer from “microscale mechanisms” to “macroscale responses” [139].

7.1.1 Decoupling between Microscale Mechanisms and Macroscale Behavior

At the microscale (micrometer to millimeter scale), fluid flow through rough fractures has long deviated from the classical linear Darcy’s law, exhibiting non-Darcy characteristics dominated by inertial effects. Flow patterns include vortices, recirculation zones, and high-velocity jets—complex phenomena that require nonlinear models such as the Forchheimer equation for accurate description [140]. Simultaneously, water–rock chemical reactions (dissolution/precipitation) at mineral surfaces proceed slowly at molecular diffusion rates, yet they irreversibly alter fracture geometry, ultimately governing the long-term permeability and mechanical strength of the rock mass at the macroscale. However, in engineering-scale (meter to kilometer) numerical models, these microscale processes are typically either highly parameterized or entirely neglected. For instance, the permeability tensor in ECMs functions as a “black-box” parameter, disconnected from its underlying microscale physical origins. Similarly, while the “shape factor” in dual-porosity models attempts to characterize matrix–fracture exchange, its physical definition remains ambiguous—particularly for sparse fracture networks [141].

This scale disconnect leads to systematic biases in model predictions. Reaction rates or permeability parameters measured in the laboratory often differ by orders of magnitude when extrapolated to field scales. The root cause lies in the inability of macroscopic models to capture the nonlinear cumulative effects of microscale mechanisms and the amplification of spatial heterogeneity across scales.

7.1.2 Building Scale Bridges: Homogenization Techniques and Nested Coupling Strategies

To bridge this gap, future research must move beyond simplistic “scaling-up” approaches and develop cross-scale modeling methods capable of faithfully transferring physical mechanisms across scales. Two cutting-edge directions hold particular promise:

• (1)Advanced Homogenization Techniques: Traditional homogenization methods (e.g., volume averaging) often fail when applied to strongly heterogeneous and discontinuous media. Next-generation homogenization theories must accommodate nonlinear constitutive relationships and dynamic boundary conditions. For instance, the multiscale FEM, based on Representative Volume Elements (RVEs), embeds a microscale RVE model at each macroscopic integration point to compute local nonlinear responses—such as stress-induced permeability changes—in real time. The homogenized results are then fed back into the macroscale model [142]. This approach enables microscopic physical laws—such as Forchheimer flow or chemical reaction kinetics—to be explicitly embedded into macroscopic governing equations in the form of constitutive relationships, thereby achieving faithful transfer of physical mechanisms across scales.
• (2)Nested Coupling Strategies: This approach acknowledges that different physical processes dominate at different scales and employs scale-appropriate models in different regions, dynamically exchanging data to achieve coupling. For example, in the near-wellbore stimulation zone of an EGS, a high-fidelity DFN model or LBM can be used to simulate microscale fracture dilation and non-Darcy flow; in the far-field reservoir, a computationally efficient dual-porosity model suffices. The two domains are bidirectionally coupled through shared boundary conditions—such as pressure, flow rate, and temperature. This “local high-fidelity + global efficiency” architecture captures critical physical details where they matter most while maintaining manageable overall computational cost, making it an ideal solution for large-scale, multiscale engineering problems.

Moreover, data-driven cross-scale correlation models are emerging as a powerful new tool. By leveraging ML—particularly DNNs—these models can directly learn the complex, nonlinear mapping between microscale features (e.g., porosity distribution, fracture roughness) and macroscale effective properties (e.g., equivalent permeability, elastic modulus) from large datasets of microscale simulations or high-resolution digital rock images [143]. Unlike traditional approaches, this method does not rely on explicit physical derivation; instead, it efficiently constructs an “empirical bridge” across scales, making it especially valuable for complex systems where underlying physical mechanisms remain poorly understood.

In summary, the future of multiscale coupled modeling does not lie in pursuing a single “universal” model, but in developing a flexible, open, and integrative multiscale collaborative framework [144]. Such a framework should be capable of intelligently selecting and coupling sub-models at appropriate scales based on the specific problem at hand, ensuring that microscale physical mechanisms are faithfully reflected in macroscale predictions. Only through this approach can we achieve a true cognitive leap—from “seeing a single tree” to “understanding the entire forest”—and lay a robust scientific foundation for the safety and efficiency of deep Earth engineering.

7.2 Model and Parameter Uncertainty

In the study of multiphysics coupling in fractured rock masses, the reliability of model predictions is fundamentally constrained by the accuracy of input parameters. However, obtaining high-fidelity THMC coupling parameters—especially under the extreme conditions typical of deep subsurface environments, such as high temperature, high pressure, and strong chemical corrosivity—represents one of the most severe technical bottlenecks in the field today [145]. This parameter uncertainty is not merely a matter of measurement error; rather, it arises from the combined effects of physical complexity, limitations of experimental techniques, and scale-dependent phenomena. Consequently, it creates a persistent and often irreducible gap between numerical simulations and field observations.

7.2.1 Inherent Difficulties in Parameter Acquisition under High Temperature and High Pressure

Deep rock masses—such as those in EGS reservoirs or host rocks surrounding nuclear waste repositories—typically exist under extreme conditions of several hundred degrees Celsius and tens of megapascals of pressure. Under such harsh thermo-hydro-chemical environments, the constitutive relationships between rock and fluid exhibit strong nonlinearity and path dependence. For example:

• Stress sensitivity of permeability: Under high in-situ stress, fracture closure and shear dilation exhibit highly nonlinear behavior that is strongly dependent on stress history. Conventional laboratory triaxial tests struggle to replicate the complex in-situ stress paths, causing laboratory-derived stress-permeability relationships to fail when extrapolated to field condition [146].
• Temperature dependence of thermal properties: Parameters such as rock thermal conductivity, thermal expansion coefficient, and fluid viscosity and density vary significantly with temperature. In high-temperature zones above 200°C, even minor errors in these parameters can be exponentially amplified through the THM coupling chain, leading to substantial deviations in predicted thermal recovery efficiency or temperature field distribution from actual field behavior [147].
• Scale effects in chemical reaction rates: Mineral dissolution/precipitation rates measured in the laboratory are typically obtained under idealized conditions. At the field scale, however, these rates are strongly controlled by fluid transport, fracture geometry, and mineral spatial distribution, often resulting in discrepancies of several orders of magnitude compared to laboratory values [148].

More critically, many key parameters—such as fracture surface friction coefficient, Biot-Willis coefficient α , and TM coupling coefficients—are inherently emergent properties whose values depend on the system’s overall state rather than being intrinsic material constants [149]. This makes direct measurement at small scales and subsequent extrapolation to engineering scales extremely challenging.

7.2.2 Urgent Need for High-Fidelity In Situ Experiments and Extreme-Environment Sensing Technologies

To overcome this bottleneck, the research focus must shift from “relying on laboratory calibration” to “developing in situ validation capabilities”. Currently, there is an urgent demand for the following two types of technologies:

• (1)High-Fidelity In Situ Coupled Experimentation Techniques:

Current laboratory equipment is mostly designed for single physical fields or mild conditions, making it inadequate for replicating the true multiphysics environments found in EGS or nuclear waste repositories. There is an urgent need to develop next-generation, field-scale physical simulation platforms capable of simultaneously applying high temperature, high pressure, reactive fluids, and mechanical perturbations under controlled conditions to directly observe the dynamic response of fracture networks under coupled THMC processes [150]. For instance, a large-scale true-triaxial seepage–thermal–mechanical coupled testing system could simulate the combined effects of kilometer-scale in-situ stress fields and thermal pulses, providing the most field-representative validation benchmark for numerical models.

• (2)Intelligent Sensing and Monitoring Technologies for Extreme Environments:

Direct field measurement is the ultimate means to reduce parameter uncertainty, yet current sensors commonly suffer from short lifespans, low accuracy, and frequent failure under deep subsurface high-temperature and high-pressure conditions. Future research must deeply integrate with materials science and microelectronics to develop a new generation of intelligent sensors capable of withstanding extreme environments:

• Distributed fiber-optic sensing: Based on Brillouin or Raman scattering principles, this technology enables continuous, real-time monitoring of temperature, strain, and pore pressure along boreholes or tunnels, with spatial resolution down to the meter scale. It is passive, immune to electromagnetic interference, and highly durable in harsh environments.
• Nanofunctional material-based sensors: Examples include temperature-sensitive quantum dots or piezoelectric nanoparticles that can be injected into fracture networks. These nanoscale probes enable in situ monitoring by remotely detecting their response signals to infer local physical fields—effectively acting as “microscopic sensors” embedded directly within the rock matrix.
• Self-powered and wireless transmission technologies: Address the critical challenges of power supply and data retrieval for deep subsurface sensors, enabling long-term, real-time, and unattended monitoring.
• Furthermore, data assimilation techniques will serve as the critical link between field monitoring and numerical modeling. By continuously integrating sparse and noisy field observations—such as microseismic events, temperature, and pressure—into numerical models in real time, data assimilation dynamically updates model parameters and state variables. This significantly reduces predictive uncertainty and enables a paradigm shift from “open-loop simulation” to “closed-loop forecasting” [151].

In summary, model and parameter uncertainty remains the central obstacle preventing fractured rock multiphysics coupling research from advancing from qualitative understanding to quantitative prediction [152]. Overcoming this challenge cannot rely solely on algorithmic improvements; it demands a dual-driven approach combining experimental innovation and sensing technology breakthroughs. A unified, four-in-one parameter calibration and validation framework—integrating laboratory experiments, in situ testing, field monitoring, and data assimilation—must be established. Only through such an integrated system can we provide truly reliable scientific foundations for the safe design and risk management of deep Earth engineering projects [153].

7.3 Characterization of Complex Fracture Networks

The core challenge in accurately modeling natural fractured rock masses lies in realistically and efficiently characterizing their inherent high degree of randomness, complex three-dimensional topological structure, and multi-fracture intersection characteristics. These geological features are not mere geometric details—they are decisive factors that directly control fluid flow pathways, solute transport efficiency, and mechanical stability. Current mainstream modeling approaches universally encounter theoretical bottlenecks and computational difficulties when dealing with such complexity, creating an urgent need for next-generation intelligent algorithms to overcome the dual constraints of representation accuracy and computational efficiency.

7.3.1 Modeling Challenges Posed by Stochasticity, 3D Topology, and Fracture Intersections

Natural fracture networks form complex systems comprising randomly interweaving fractures of diverse scales, orientations, and densities in three-dimensional space. The stochastic nature of their 3D topology is illustrated in Fig. 11 through a comparison between measured fracture traces and simulated ones; the simulated traces are generated by intersecting a sampling window with a 3D stochastic discontinuity model. This correspondence between field observations and numerical simulations not only validates the rationality of the stochastic modeling approach—the primary modeling challenges arise in three key aspects:

Figure 11: Comparison between measured traces and numerical model. (a) In-situ measured discontinuities traces. (b) Simulated traces formed by intersection between sampling window and stochastic model. Adapted with permission from Reference [154]. Copyright © 2015, Springer.

• (1)Geometric stochasticity and statistical characterization: Fracture parameters—such as location, size, aperture, and roughness—follow specific probability distributions (e.g., power-law or log-normal distributions). Although DFN models can generate stochastic networks via Monte Carlo methods, the calibration of their statistical parameters (e.g., fracture density, mean aperture) heavily relies on sparse and costly field data—such as borehole imaging or outcrop mapping. Uncertainty in these data propagates directly into the model, resulting in large variances in simulation outcomes and limiting their applicability for deterministic safety assessments in high-risk engineering projects [155].
• (2)Three-dimensional topological connectivity: The connectivity of fractures in 3D space is the key determinant of the “channeling effect”. An isolated fracture contributes negligibly to global flow, whereas a connected pathway formed by intersecting fractures can carry over 90% of the total fluid flux. Traditional DFN models often fail during mesh generation when handling 3D fracture intersections—either due to geometric conflicts or the creation of severely distorted elements—compromising both computational stability and accuracy.
• (3)Physical modeling of fracture intersections: At fracture intersection points, fluid flow and mechanical behavior undergo dramatic changes. The classical cubic law breaks down, and flow exhibits strong three-dimensional vortices and significant pressure losses. Simultaneously, stress concentration at these junctions greatly increases the likelihood of new fracturing or shear slip. Current models typically simplify intersections as “nodes” and approximate their hydraulic behavior using empirical “intersection rules” or “transfer coefficients”, which lack a solid physical basis and constitute a major source of predictive uncertainty in numerical simulations.

These challenges give rise to a fundamental contradiction: although DFN models can, in principle, represent fracture networks with the highest fidelity, their high computational cost and demanding data requirements make them impractical for large-scale engineering applications. Conversely, while ECM or dual-porosity models offer computational efficiency, their excessive geometric simplifications prevent them from capturing critical localized flow pathways [156].

7.3.2 ML: Solving Inverse Modeling and Parameter Constraint Challenges

ML, particularly DNNs, offers a revolutionary approach to addressing the aforementioned characterization challenges. Its core value lies in the ability to automatically learn the complex, nonlinear mapping between fracture network geometric features and macroscopic system responses—such as pressure, temperature, and microseismic signals—from sparse and noisy observational data, thereby enabling efficient inverse modeling and parameter constraint.

• (1)Inverting fracture networks from observational data: Traditional inversion methods—such as trial-and-error or optimization algorithms—suffer from exponentially increasing computational costs and a tendency to converge to local optima when dealing with high-dimensional parameter spaces (involving tens of thousands of fracture parameters). In contrast, DNNs can establish an end-to-end mapping from “observational data → fracture network parameters” through supervised learning [157]. For instance, Fig. 12 illustrates the complex dynamic evolution of fracture networks captured by high-resolution CT scanning. Such image sequences contain rich information regarding pore structure and geometric topology, presenting a significant challenge for traditional manual interpretation or simple segmentation methods to accurately capture subtle, nonlinear changes. This highlights the specific advantage of ML: AI models can directly utilize these complex CT image sequences as input data, automatically extracting key features to rapidly and accurately identify and quantify the evolution of fracture network patterns. By taking the spatial distribution, magnitude, and timing of microseismic events as input, a trained DNN can predict the orientation, density, and aperture distribution of subsurface fractures—achieving speeds several orders of magnitude faster than conventional inversion techniques.

Figure 12: Pixel density distribution of a CT-scan of fractured concrete. Adapted with permission from Reference [158]. Copyright © 2023, Elsevier.

• (2)Data-driven parameter calibration and uncertainty quantification: DNNs can serve as “surrogate models” to replace computationally expensive physics-based simulations for massive parameter sweeps. Using methods such as BNNs, they can not only predict optimal parameter sets but also quantify the uncertainty of predictions, providing risk-informed decision intervals for engineering applications. For example, in an EGS project, historical injection pressure and production temperature data can be used to train a DNN to invert the reservoir’s permeability and stress fields; the resulting model can then directly inform the optimization of injection strategies for subsequent operational phases [159].
• (3)Intelligent generation and augmentation: When field data are extremely scarce, GANs or VAEs can generate large ensembles of geologically plausible “synthetic fracture networks” based on limited statistical priors. These synthetic datasets can be used to augment training data or enable Monte Carlo–based risk analysis, effectively alleviating data scarcity bottlenecks.
• (4)Cross-scale correlation modeling: ML excels at integrating multi-source, multi-scale data. For example, by using both macro-scale geophysical attributes (e.g., seismic velocity, resistivity) and micro-scale digital rock image features as joint inputs, a DNN can be trained to predict a reservoir-scale Fracture Intensity Index (FSI), enabling an intelligent leap from “point” (borehole) to “area” (reservoir) characterization.

In summary, the accurate characterization of complex fracture networks has evolved from a purely geometric modeling problem into a systems-level endeavor that integrates data, models, and intelligence. ML is not intended to replace physics-based models; rather, it equips them with an “intelligent eye” and “accelerating wings”, enabling automatic extraction of critical information from massive datasets and efficient constraint of model parameters—ultimately allowing us to “see through” and “preview” the subsurface world. Looking ahead, with the advancement of cutting-edge technologies such as Physics-Informed Machine Learning (PIML) and GNNs, we are poised to develop a new generation of intelligent characterization models that faithfully capture geological complexity while maintaining physical consistency and computational efficiency—paving a transformative path toward safer and more effective deep Earth engineering.

7.4 Theoretical Model Innovation

The theoretical foundation of current multiphysics coupling research in fractured rock masses is primarily built upon the framework of classical Continuum Mechanics (CM). However, this theoretical system faces fundamental limitations when applied to naturally discontinuous, multi-fractured media like rock [160]. The continuum assumption requires the material to be defect-free and infinitely differentiable at all scales—a premise fundamentally incompatible with the pervasive presence of discontinuities (joints, faults, bedding planes) in fractured rock and the processes of fracture initiation, propagation, and coalescence under stress. To overcome this theoretical bottleneck, developing new modeling frameworks that inherently accommodate discontinuities has become the most cutting-edge and potentially transformative research direction in the field. Among these, Peridynamics (PD) and Continuum-Discontinuum Coupling Methods represent two of the most promising pathways.

7.4.1 Beyond CM: Inherent Limitations of Classical Theory

Traditional numerical methods—such as the FEM—must either predefine crack paths or rely on complex fracture criteria (e.g., maximum tensile stress criterion, energy release rate criterion) when simulating rock fracturing. This approach is inherently descriptive rather than predictive. When crack paths are unknown or involve complex branching and coalescence, such models often fail. Moreover, FEM suffers from numerical instability in strongly discontinuous problems involving large deformations or fragmentation, due to severe mesh distortion. This necessitates frequent and costly remeshing, which not only increases computational expense but also introduces ambiguity in the physical interpretation of results. These shortcomings stem from the mathematical foundation of FEM—PDEs—which require field variables (e.g., displacement, stress) to be spatially continuous and differentiable, whereas real rock fracture is an inherently discontinuous and nonlocal physical process.

7.4.2 PD: A Revolutionary Breakthrough in Nonlocal Theory

PD is a novel solid mechanics theory based on nonlocal integral equations. Its core idea is to abandon the concept of spatial derivatives and instead describe a material’s mechanical behavior through long-range interaction forces between material points and their neighbors. In the PD framework, a body is discretized into a set of finite-volume “material points”, each connected to other points within a finite distance—known as the “horizon”—via pairwise “bonds”. Material damage is represented by bond breakage: when the relative displacement between two connected points exceeds a critical threshold, the bond breaks. This mechanism naturally captures crack initiation, propagation, and branching without requiring predefined crack paths or additional fracture criteria.

The core advantages of PD lie in:

• Native handling of discontinuities: Crack initiation and propagation emerge naturally from the model solution rather than through artificial intervention. This gives PD unparalleled advantages in simulating complex fracture networks and fragmentation processes.
• Mesh-free and derivative-free formulation: Since PD does not rely on spatial derivatives, it is immune to mesh distortion, making it especially well-suited for simulating large deformations and extreme failure scenarios.
• Nonlocality: The “horizon” in PD inherently incorporates nonlocal material interactions, enabling a more physically realistic representation of phenomena such as stress wave propagation and dynamic fracture.
• To date, PD theory has evolved from the original bond-based formulation to the more general state-based framework and has been successfully extended to multiphysics coupling domains such as TM and hydro-mechanical processes. It has demonstrated significant potential in simulating engineering problems like hydraulic fracturing, rockbursts, and slope instability, and is increasingly regarded as a “game-changer” for the future of theoretical modeling in rock mechanics [161].

Despite its conceptual elegance, the practical deployment of PD in realistic subsurface environments demands explicit integration with THMC coupling physics. In geothermal reservoir stimulation, for instance, thermal contraction induces tensile stresses that concentrate at pre-existing fracture tips—a classical stress singularity under continuum theory. Similarly, during CO2-brine-rock interactions, localized mineral dissolution creates sharp porosity/permeability fronts, while pore pressure diffusion across fracture-matrix interfaces generates discontinuous hydraulic gradients. These multiscale, multiphysics discontinuities challenge the differentiability assumptions of PDE-based models.

PD’s nonlocal integral formulation inherently circumvents this limitation: by replacing spatial derivatives with finite integrals over a horizon, field variables need not be continuous or differentiable. This provides a mathematically consistent framework to resolve stress concentrations at crack tips and abrupt changes in temperature, concentration, or pore pressure without enrichment or regularization.

However, extending PD to fully coupled THMC systems introduces significant numerical challenges. First, formulating nonlocal constitutive laws that consistently couple thermal expansion, fluid pressure, chemical potential, and bond failure remains an open problem. Second, the explicit time integration typically used in PD becomes prohibitively restrictive when stiff diffusive (e.g., heat conduction) or reactive (e.g., dissolution kinetics) processes are included, often requiring time steps orders of magnitude smaller than those in mechanical-only simulations. Third, the dense interaction matrix in 3D nonlocal models leads to O(N2) memory and computational complexity, which is exacerbated when multiple fields (T, H, M, C) each require their own horizon-based operators. Addressing these issues—through implicit PD schemes, adaptive horizons, or hybrid local–nonlocal coupling—is essential for PD to transition from a fracture mechanics tool to a predictive THMC simulator.

7.4.3 Continuum-Discontinuum Coupling Methods: Integrating Macroscopic Efficiency with Microscopic Fidelity

Whereas PD seeks to reconstruct the theoretical foundation from the ground up, continuum-discontinuum coupling methods represent a pragmatic engineering strategy that aims to combine the computational efficiency of continuum models with the fracture-representation fidelity of discrete approaches. The core idea is to apply different modeling scales within the same computational domain: use an efficient FEM in regions far from fractures (“continuum zones”), and switch to discrete models—such as the DEM or extended Finite Element Method (XFEM)—in areas where fracturing is anticipated or has already occurred (“discontinuum zones”). Seamless transfer of data and energy across the coupling interface is achieved through carefully designed algorithms, ensuring physical consistency and numerical stability.

Typical coupling frameworks include:

• FEM-DEM coupling: FEM is used at the macroscale to simulate the overall deformation of the rock mass, while DEM particles are embedded in localized high-stress zones or along predefined weak planes to model fracture opening, sliding, and fragmentation. This approach captures both the global stress field and the detailed local fracture processes, and has been successfully applied in tunnel excavation and hydraulic fracturing simulations.
• FEM-XFEM coupling: The XFEM introduces enrichment functions into the FEM framework, enabling cracks to propagate arbitrarily within elements without requiring remeshing. The FEM-XFEM coupled model can efficiently simulate the evolution of complex crack networks while preserving the original FEM mesh structure.

The advantages of hybrid methods lie in:

• High computational efficiency: High-cost discrete models are employed only in critical regions, significantly reducing the overall computational burden.
• Clear physical mechanisms: Classical CM retains its interpretability in continuous regions, while fracture behavior is accurately captured in discontinuous zones.
• Strong engineering applicability: These methods can be readily integrated into existing commercial software frameworks, facilitating adoption and practical implementation by engineers.
• In summary, the ultimate goal of theoretical model innovation is to establish a universal framework capable of consistently describing the entire spectrum of rock behavior—from continuous deformation to discontinuous fracture. Peridynamics, with its revolutionary nonlocal formulation, offers a fundamentally new mathematical foundation toward this vision. Meanwhile, continuum-discontinuum coupling methods provide practical, engineering-oriented solutions for today’s complex geomechanical challenges. Looking ahead, as computational power grows and theoretical frameworks mature, we may witness the dawn of a “post-continuum mechanics” era—one in which rock “fracturing” is no longer a modeling deficiency, but rather the central and most authentic physical attribute of the material itself [162].

7.5 Efficient Numerical Algorithms

As the physical fidelity and geometric complexity of multiphysics coupling models for fractured rock masses continue to increase, the computational challenges in their numerical solution have grown exponentially. Modern high-fidelity models—such as fully coupled THMC-DFN simulations or peridynamics—often involve millions to hundreds of millions of degrees of freedom, with governing equations characterized by strong nonlinearity, tight multiphysics coupling, and multiscale features. Traditional serial solvers and general-purpose commercial software are no longer adequate for such large-scale, highly nonlinear problems. Consequently, the development of efficient, robust, and scalable next-generation numerical algorithms has become the critical enabler for translating theoretical advances into practical engineering applications. Achieving this goal hinges on two pillars: leveraging advanced parallel computing architectures and strategically integrating AI/ML.

7.5.1 Parallel Computing and Solver Development for Large-Scale, Strongly Nonlinear Problems

Modern high-performance computing (HPC) provides the hardware foundation to address computational bottlenecks, but fully harnessing its potential requires the development of advanced algorithms specifically designed to exploit its capabilities:

• (1)Large-Scale Parallelization and Domain Decomposition Techniques:

For computationally intensive models such as DFN or FEM-DEM, large-scale parallel computing strategies are essential. By employing domain decomposition methods, the entire computational domain is partitioned into multiple subdomains, each assigned to a different computing node for concurrent solution, significantly reducing simulation time [163]. For instance, when simulating a kilometer-scale EGS reservoir, the domain can be divided into hundreds of subregions, enabling near-linear speedup on supercomputing clusters.

• (2)Efficient Solvers for Nonlinear Systems:

THMC-coupled equations fundamentally form a strongly nonlinear system of algebraic equations. The traditional Newton-Raphson method often struggles with convergence in highly nonlinear or ill-conditioned systems due to sensitivity to initial guesses or singularity of the Jacobian matrix. To address this, more robust solver strategies are required:

• Preconditioning techniques: Design efficient preconditioners—such as multigrid methods or incomplete LU (ILU) factorization—to improve the condition number of the Jacobian matrix and accelerate iterative convergence.
• Globally convergent algorithms: Incorporate trust-region methods or homotopy continuation strategies to enhance global convergence, especially when the initial guess is far from the true solution.
• Matrix-free methods: For extremely large-scale problems, avoid explicitly assembling and storing the full Jacobian matrix; instead, approximate its action via automatic differentiation or finite differences, drastically reducing memory consumption.

• (3)Efficient Algorithms for Multiphysics Coupling:

Fully implicit solution of strongly coupled problems (e.g., TMH) leads to a dramatic increase in the size and complexity of the global system matrix. To mitigate this, partitioned solvers or loosely coupled algorithms are essential. For example, a staggered iterative strategy—such as “implicit solution of the mechanical field + explicit update of the hydraulic field” [164]—can significantly reduce the computational cost per time step while maintaining numerical stability.

• (4)Adaptive Meshing and Dynamic Load Balancing:

When simulating localized processes such as fracture propagation or shear slip, Adaptive Mesh Refinement (AMR) dynamically refines the mesh in critical regions—such as fracture tips—while maintaining a coarser discretization in the far field. This strategy optimizes computational resources by preserving accuracy where it matters most without unnecessary overhead elsewhere. Coupled with dynamic load balancing algorithms, AMR ensures uniform computational workload across all processing nodes, preventing performance bottlenecks caused by localized mesh refinement.

7.5.2 Strategic Integration of AI/ML for Acceleration and Optimization

AI/ML is not merely an independent predictive tool; its deeper value lies in its role as a “computational accelerator” and “intelligent optimizer”, deeply embedded within traditional numerical solution workflows to achieve revolutionary gains in computational efficiency:

• (1)AI Surrogate Modeling:

Train DNNs or GPR models as high-speed surrogates for computationally expensive physics-based simulations. For instance, a DNN can replace time-consuming LBM or DFN simulations, reducing a single simulation from hours to seconds—enabling rapid parameter sensitivity analysis, uncertainty quantification, or real-time optimization control [165]. In EGS design, engineers can use such surrogate models to evaluate thousands of design scenarios within minutes, a task that would take months with conventional methods.

• (2)AI-Driven Solver Acceleration:

• Intelligent preconditioners: Use GNNs to learn the structural features of the system matrix and automatically generate optimal preconditioners, significantly accelerating the convergence of linear solvers [166].
• Convergence path prediction: Train ML models to forecast the convergence behavior of nonlinear iterations, enabling dynamic adjustment of step sizes or switching between solution strategies to avoid unproductive iterations.
• Adaptive time-step control: Employ RL-based intelligent controllers that dynamically adjust the time step based on real-time system states, maximizing computational efficiency while maintaining solution accuracy.

• (3)AI-Optimized Model Parameters and Design Variables:

Integrate AI/ML with optimization algorithms—such as genetic algorithms or Bayesian optimization—to establish an “intelligent optimization loop”. For example, in tunnel support design, a surrogate model can be trained with support parameters (e.g., bolt spacing, shotcrete thickness) as inputs and key performance indicators (e.g., rock mass displacement, seepage rate) as outputs. Bayesian optimization then automatically explores the design space to identify the optimal support configuration, drastically reducing the design cycle from weeks or months to hours or days.

• (4)Multifidelity Modeling and Transfer Learning:

Construct a hybrid architecture combining a low-fidelity physics-based model (e.g., ECM) with a high-fidelity AI corrector. The AI component learns the error mapping between the low-fidelity and high-fidelity models (e.g., DFN), enabling near-high-fidelity predictions at a fraction of the computational cost. Transfer learning further enhances efficiency by allowing a model trained in one context—such as a laboratory-scale experiment—to be rapidly adapted to a different context, such as a field-scale engineering application, significantly reducing the need for redundant training and data collection.

In summary, the future of efficient numerical algorithms lies in the deep integration of computational science, applied mathematics, and AI. Neither hardware upgrades alone nor incremental algorithmic improvements are sufficient to meet the ever-growing computational demands. Only by treating AI/ML as a “meta-algorithm”—strategically embedding it into the core components of numerical simulation, from mesh generation and equation solving to parameter optimization—can we build a truly next-generation “intelligent computational engine”. This engine will not only overcome current bottlenecks in large-scale, strongly nonlinear problems but also unlock transformative possibilities for multiphysics coupling research in fractured rock masses, propelling the field from being merely “computable” toward a new era of “real-time optimization” and “autonomous decision-making”.

In summary, research on multiphysics coupling in fractured rock masses is undergoing a pivotal paradigm shift—spanning from the fundamental challenge of characterizing complex fracture networks, through groundbreaking theoretical innovations, to the engineering implementation of highly efficient numerical algorithms. Table 9 systematically outlines these core challenges and future directions, providing a clear roadmap for the technological evolution of the field.

8 Conclusion

8.1 Summary of Key Research Findings

This review synthesizes recent advances in understanding fluid flow and THMC coupling in fractured rock masses, highlighting a unifying principle: hydraulic behavior is inherently nonlinear, stress-sensitive, path-dependent, and dynamically evolving—defying simplistic or static representations.

A central finding is the bidirectional “closure-dilation” mechanism governed by effective stress. Permeability decays exponentially as normal stress closes fractures, yet can surge during shear-induced dilation. Critically, this process exhibits strong hysteresis: permeability trajectories differ during loading versus unloading, meaning the current hydraulic state encodes the full stress history. Consequently, models ignoring path dependence inevitably misrepresent real rock behavior.

This mechanism underpins major engineering phenomena: in EGS, cold-water injection induces thermal stresses that trigger shear slip and fracture activation; in tunneling, excavation unloading can cause shear dilation and water/mud inrush; in nuclear waste disposal, thermal pulses alter fracture connectivity and radionuclide transport pathways. Thus, predicting and controlling dilation is pivotal to safety and performance.

Modeling approaches must be selected based on scale, objective, and data availability:

• ECM offer computational efficiency for large-scale far-field analysis but require a REV;
• DFN models provide high-fidelity near-field simulation of DFNs, ideal for critical zones like EGS wellbores or tunnel walls, albeit at high computational cost;
• Dual-porosity/permeability models balance efficiency and physics for reservoir-scale simulations with multiscale fractures.

Complementary numerical methods (FEM, DEM, LBM, BEM) should be chosen according to dominant physical processes (e.g., large deformation, non-Darcy flow). Crucially, model validity depends on alignment between its assumptions and actual geological conditions—no matter how elegant, a geologically detached model lacks predictive power.

In essence, the field has evolved from empirical descriptions toward an integrated, multiscale, multiphysics framework. Future progress hinges on two challenges: (1) faithfully capturing the dynamic stress–deformation–seepage coupling chain, and (2) reliably upscaling microscale mechanisms to inform macroscale predictions.

8.2 Implications for Engineering Practice

The practical value of this synthesis lies in actionable, scenario-specific modeling strategies grounded in fundamental physics and comprehensive site characterization—the cornerstone for reducing uncertainty.

• Nuclear Waste Disposal: Adopt a far-field–near-field nested coupling framework. Use ECM or dual-porosity models at kilometer scales for regional flow and thermal plume prediction; deploy DFN or PINNs at meter scales around canisters to resolve bentonite hydration, THM coupling, and fracture evolution under thermal pulses.
• EGS: Implement a dynamic optimization approach. During stimulation, use DFN or FEM-DEM to simulate shear-induced dilation from cold injection; during production, shift to efficient dual-porosity or AI surrogate models for long-term forecasting, enhanced by real-time data assimilation (e.g., microseismic inversion via ML).
• Tunneling: Focus on seepage–stress coupling for hazard prevention. Apply HM-coupled FEM during excavation to manage pore pressure and grouting; use THMC models in operation to assess long-term chemical degradation, freeze-thaw damage, and support stability—with uncertainty quantified via probabilistic or ML-based methods.

Across all scenarios, integrated site characterization is indispensable:

• (1)Fuse multiscale data (microscopic to regional) to build “point-to-volume” geological models;
• (2)Conduct multiphysics co-observation (e.g., DTS/DAS for concurrent temperature/strain);
• (3)Employ Bayesian inversion, EnKF, or ML surrogates for real-time parameter calibration;
• (4)Establish a closed-loop “characterization–modeling–validation” cycle that evolves throughout the project lifecycle.

Only through tight integration of advanced modeling and thorough characterization can we move from managing uncertainty toward deterministic prediction—ensuring the safety and sustainability of deep Earth engineering.

8.3 Future Outlook

We summarize core findings and forward-looking priorities in Table 10. This structured overview underscores the interdependence among fundamental mechanisms, modeling strategies, engineering implementation, and emerging interdisciplinary frontiers—providing a roadmap for next-generation research.

To translate the promise of AI into reliable engineering practice, we propose a three-phase research roadmap aligned with data availability, model fidelity, and deployment readiness:

Short-term (1–3 years): Focus on hybridizing validated ML models (e.g., RF, shallow ANNs) with domain knowledge—such as incorporating rock mass rating rules or Darcy’s law as feature constraints. Prioritize standardization of small-scale labeled datasets from borehole logs, core imaging, and microseismic monitoring, and deploy these models in low-risk decision support tasks (e.g., preliminary rock classification or injection zone screening).

Medium-term (3–7 years): Develop physics-informed deep learning frameworks that embed conservation laws (mass, momentum, energy) as hard or soft constraints. Integrate multi-source data (geophysical, geochemical, mechanical) via multimodal architectures, and validate surrogate models against controlled field experiments (e.g., in-situ stimulation tests at underground research laboratories). Active learning strategies should guide targeted data acquisition under critical THMC conditions.

Long-term (7–15+ years): Realize autonomous, self-calibrating digital twins of deep subsurface systems capable of real-time forecasting and risk-aware control. This requires breakthroughs in few-shot domain adaptation across geological settings, explainable AI for regulatory trust, and robust uncertainty quantification under extreme data sparsity. Close collaboration among geoscientists, AI researchers, and regulators will be essential to establish certification protocols for safety-critical applications such as nuclear waste repositories or deep geothermal reservoirs.

The field is entering a transition phase where breakthroughs will stem not from isolated advances, but from deep interdisciplinary integration:

• AI & computational science will shift from auxiliary tools to core simulation engines, with PINNs, GNNs, and operator learning enabling real-time, physics-consistent solutions;
• Materials science will enable extreme-environment sensing (e.g., quantum dots) and refine chemo-mechanical constitutive laws;
• Rock mechanics will transcend continuum limits via peridynamics and hybrid continuum–discontinuum frameworks.

Three technical directions will dominate the next decade:

• (1)Physics-informed AI to embed conservation laws into neural networks, ensuring extrapolatable predictions;
• (2)Hybrid modeling, combining interpretable physics-based baselines with AI correctors for high-accuracy, high-trust forecasts;
• (3)Intelligent multiscale simulation, dynamically linking pore-, fracture-, and reservoir-scale domains via adaptive homogenization and ML-based cross-scale mapping.

Accelerating progress demands global collaboration: strengthening initiatives like DECOVALEX, building open data/model repositories, and establishing unified Verification & Validation (V&V) standards for key applications. Only through shared benchmarks and transparent validation can models gain the reliability needed for engineering decision-making.

We stand at an inflection point—where intelligent, physics-grounded digital twins of fractured rock systems are within reach. Achieving this vision is both a scientific imperative and an engineering necessity for the sustainable development of subsurface resources.