Numerical investigation of the effects of normal stress and shear displacement on fracture closure and permeability evolution in rough granite fractures

Fracture closure and permeability evolution in rough granite fractures considering normal stress and shear displacement are critical for numerous underground engineering applications, including enhanced geothermal systems, radioactive waste repositories, and underground energy storage. Despite extensive research on stress-dependent fracture permeability, the combined effects of fracture surface roughness, normal stress, and shear displacement on fracture closure and permeability evolution remain insufficiently quantified. In this study, a permeability evolution model incorporating roughness, normal stress, and shear displacement was proposed based on the finite element method. A series of three-dimensional fracture models with four roughness levels (JRC = 10.5, 12.2, 13.6, and 16.7) were reconstructed from high-resolution laser-scanned tensile fractures. Based on the finite element method, fracture closure simulations were performed under normal stresses ranging from 0.1 to 40 MPa and shear displacements ranging from 0 to 10 mm, resulting in a total of 144 simulation cases. The resulting aperture fields were subsequently used as geometric inputs for steady-state flow simulations based on the cubic law. The results show that fracture aperture decreases nonlinearly with increasing normal stress, while increasing shear displacement and fracture roughness enlarge the initial aperture through shear-induced dilatancy. An empirical model relating the initial aperture to JRC and shear displacement was established. Fracture permeability exhibits a strong negative correlation with normal stress and positive correlations with both shear displacement and roughness. Moreover, hydraulic aperture and fracture aperture display a strong linear correlation across all loading conditions, indicating that permeability variations are primarily governed by changes in the initial aperture. These findings provide a quantitative framework for predicting hydro-mechanical behavior in fractured granite reservoirs.

Fluid flow through fractured rock masses is a fundamental issue that governs the performance and safety of numerous underground engineering projects, including enhanced geothermal systems (EGS), geological disposal of high-level radioactive waste, underground energy storage, and hydrocarbon extraction (Zhao and Brown, 1992; Rutqvist and Stephansson, 2003; Selvadurai et al., 2015; 2018; Zhong et al., 2023). As essential units forming complex fracture networks, single rough fractures serve as the basis for understanding fluid flow behavior and associated hydraulic properties in fractured rock masses (Zhong et al., 2021; Jiang et al., 2022). Therefore, accurately characterizing the evolution of contact states during fracture closure is fundamental to understanding the coupled hydro-mechanical behavior of fractured hard rocks.

The evolution of fracture permeability is a coupled hydro-mechanical process governed by the interaction among fracture surface roughness, normal stress, and shear displacement (Rong et al., 2020; Chen et al., 2021b; Li et al., 2024; Cao, 2026). Fracture surface roughness controls the initial contact geometry and the spatial distribution of apertures, thereby determining the available flow pathways (Brown, 1987). Normal stress promotes fracture closure by compressing opposing surface asperities, which reduces the mean aperture and constricts flow channels (Walsh, 1981; Chen et al., 2015). In contrast, shear displacement modifies the matching condition between the upper and lower fracture surfaces, potentially generating dilatancy and additional void space that enhances permeability (Yeo et al., 1998; Fang et al., 2018). These factors, however, do not act independently (Sawayama et al., 2020). The influence of roughness on permeability is strongly affected by the applied normal stress, while shear-induced dilation is progressively suppressed by compressive closure. Despite the practical importance of these interactions, quantitative models that simultaneously account for roughness, normal stress, and shear displacement remain limited, particularly for three-dimensional rough fracture geometries subjected to coupled normal and shear loading.

Extensive experimental investigations have been conducted to characterize fracture permeability under varying stress conditions. Selvadurai (2015) examined the permeability evolution of a fracture in Barre granite under axial stress and observed pronounced permeability hysteresis after complete stress relief. Huo et al. (2024) and Huo and Meguid (2025) systematically investigated the stress-induced permeability evolution in Lac du Bonnet granite using steady-state flow tests combined with high-resolution laser scanning, demonstrating that fracture permeability exhibits a quadratic exponential decline with increasing confining stress and that permeability hysteresis persists after unloading. Fang et al. (2018) conducted shear-flow experiments on fabricated fractures with specified roughness and reported that permeability enhancement associated with dilatancy is largely controlled by fracture surface roughness. Zhong et al. (2023) combined flow tests with numerical simulations to investigate fluid flow through fractures surrounded by permeable matrix, and proposed a threshold for the ratio of matrix to fracture permeability beyond which matrix contributions should not be ignored. Collectively, these studies have significantly advanced the understanding of fracture closure mechanisms and their impact on permeability. However, significant technical challenges remain in the physical experimental investigation of coupled fracture shear-flow behavior. After intact rock is split to form mated fracture surfaces, the application of shear displacement causes relative sliding between the originally interlocked fracture surfaces, leading to substantial changes in the aperture field. When flow tests are further superimposed under such shear conditions, the continuous evolution of the fracture surface contact state during shearing makes it difficult to maintain stable fluid boundary sealing conditions. Consequently, there are few reported laboratory experimental studies on the coupled shear-flow behavior of fractures.

In addition, significant advances in numerical modeling have enabled detailed investigation of fluid flow through rough fractures under complex loading conditions. Chen et al. (2021b) adopted experimentally measured fracture geometries under increasing confining stresses to establish three-dimensional real fracture models, and solved the full Navier-Stokes equations to examine the impacts of stress-induced void and contact alterations on flow characteristics. Their results showed that the streamlines become more channeled with observable tortuosity under high stress, and that the increase in contact area enhances the non-linearity of pressure drop within the fracture. Rong et al. (2020) proposed two geometric parameters to characterize fracture geometry and established empirical models linking the nonlinear coefficient to these parameters based on hydraulic tests conducted on 50 fractured granite samples. Li et al. (2024) established a shear-flow model based on the discrete element method and pore network model (DEM-PNM) to investigate permeability evolution of natural fractures in granite during shearing, demonstrating that fracture roughness plays a dominant role in permeability evolution. Despite these advances, most existing numerical studies have focused on either normal stress or shear displacement in isolation, and systematic investigations that examine their coupled effects on three-dimensional rough fracture models reconstructed from real rock surfaces remain scarce.

To address this gap, this study employs the finite element method (FEM) to systematically investigate the effects of fracture surface roughness, normal stress, and shear displacement on fracture aperture distribution and permeability evolution. Specifically, based on the point cloud data of rough granite fracture surfaces obtained through high-resolution three-dimensional laser scanning, fracture models with different shear displacements (0, 2, 4, 6, 8, and 10 mm) are constructed by offsetting the upper and lower surfaces. A series of fracture closure simulations under normal stresses ranging from 0.1 to 40 MPa are then conducted for fractures with four different roughness levels, followed by steady-state flow simulations. The objective of this study is to establish a quantitative model for fracture closure and permeability evolution under the synergistic control of roughness, shear displacement, and normal stress, thereby providing a parametric basis for coupled hydro-mechanical analysis of fractured rock masses.

Sample preparation and fracture induction

Fresh granite samples were collected from a quarry in Suizhou City, Hubei Province, China. Optical microscopy analysis showed that the rock consisted primarily of approximately 33% quartz, 34% feldspar, 19% amphibole, and 14% accessory minerals. To determine the basic mechanical properties of the intact granite, uniaxial compression and Brazilian splitting tests were conducted in accordance with the testing procedures recommended by the International Society for Rock Mechanics (Ulusay, 2014). The measured mechanical parameters were subsequently used as input data for FEM simulations. The main test results are summarized in Table 1.

To obtain specimens containing rough fracture surfaces, homogeneous and intact granite blocks were processed into cubic specimens with side lengths of 100 mm using a high-precision rock cutting and grinding machine. The end surfaces were carefully polished to ensure adequate flatness and perpendicularity. Artificial fractures were generated using a specialized loading fixture designed according to principles similar to those of the Brazilian splitting test (Meng et al., 2016). Two V-shaped loading attachments were mounted on the upper and lower platens of the testing machine, as illustrated in Figure 1a. Prior to loading, a preload of 1 kN was applied by slowly lowering the upper platen to ensure firm contact between the fixture and the specimen. Subsequently, servo-controlled loading was imposed at a constant displacement rate of 0.01 mm/min until tensile failure occurred along the predetermined splitting plane. This procedure produced a single, clean, and through-going fracture without noticeable secondary cracking (Selvadurai, 2015). The loading configuration effectively confined fracture propagation to the intended plane and preserved good matching conditions between the opposing fracture surfaces. By repeating the Brazilian splitting process on multiple cubic specimens, four groups of artificial fractures with different surface roughness characteristics were successfully prepared. A representative specimen is shown in Figure 1a.

Three-dimensional reconstruction of fracture surfaces

Fracture surface characterization

To accurately characterize the geometric morphology of the fracture surfaces and provide reliable inputs for simulations, a three-dimensional laser scanner (SIMSCAN 30, maximum spatial resolution of 0.02 mm) was employed to acquire point cloud data of the fracture topography. The scanned data were exported in STL format and processed using GOM Inspect software to reconstruct the fracture surface geometry. A spatial resolution of 0.4 mm was adopted for the generation of the digital fracture models (Figure 1a) (Candela and Brodsky, 2016). Figure 2 presents the reconstructed fracture surfaces with four different geometric characteristics.

To quantitatively evaluate fracture surface roughness, the Joint Roughness Coefficient (JRC) proposed by Barton and Choubey (1977) was adopted. Two-dimensional profiles for roughness analysis were extracted parallel to the splitting direction at intervals of 1.0 mm. A total of 11 profiles were obtained from each fracture surface, as illustrated in Figure 1b. The JRC value of each profile was calculated using Equations 1, 2. The overall JRC of a fracture surface was defined as the arithmetic mean of the JRC values of all extracted profiles, which served as the quantitative index for characterizing the global roughness of the fracture surface (Huang et al., 2017; Ye and Ghassemi, 2018; Hu et al., 2022).

The JRC of each surface is quantitatively described as follows: where Z 2 is the root mean square of the first deviation of the profile, and L is the evaluation length. Based on this method, four groups of fracture specimens with different roughness levels were obtained, with measured JRC values of 10.5, 12.2, 13.6, and 16.7, respectively.

Fracture aperture characterization

Characterization of the fracture aperture, defined as the spatial separation between the upper and lower fracture surfaces, is essential for constructing three-dimensional realistic fracture models. The fracture aperture was determined from the height difference between the opposing fracture surfaces. First, the scanned point cloud data were imported into Rhino software. The upper fracture surface was generated by duplicating the lower surface and translating it along the surface normal direction. In this study, the translation distance was set to 0.5 mm to provide an initial state. This approach facilitates mesh generation, as direct use of independently scanned surfaces may lead to misalignment between the upper and lower meshes, complicating the numerical analysis. The resulting initial fracture model possessed identical roughness morphologies on the upper and lower surfaces, which remained parallel and separated in the unloaded state. Subsequently, to construct fracture models corresponding to different shear displacements, the upper surface was subsequently shifted horizontally relative to the lower surface along the shear direction by 0, 2, 4, 6, 8, and 10 mm, respectively. This procedure generated a series of fracture geometries representing different shear displacement conditions. The upper and lower surfaces were represented by the spatial functions Z i, j t o p x, y and Z i, j b o t t o m x, y, and stored as point cloud datasets, as shown in Figure 3.

Because the original scanned points were irregularly distributed within the fracture region, the point cloud data were interpolated onto a regular grid using MATLAB with a grid resolution of 0.2 mm. The local fracture aperture D i, j was calculated from the elevation difference between the upper and lower surfaces at corresponding grid nodes using Equation 3 (Chen et al., 2017; Huo et al., 2024):

The mean fracture aperture was then defined as the arithmetic average of all nodal apertures using Equation 4: where m and n denote the numbers of grid nodes in the x and y directions, respectively. The interpolation onto a uniform grid minimizes errors associated with positional mismatches in the original point cloud data (Raven and Gale, 1985; Carl, 1995).

Note that any negative aperture values caused by numerical interpolation errors or local point cloud noise were set to zero to satisfy the physical non-negativity constraint of the aperture field. After the aperture calculation was completed, the aperture field data were archived and the corresponding mean value was taken as the representative aperture D m for that fracture condition. All processed geometric models are ultimately exported in STL format to ensure compatibility with the COMSOL Multiphysics.

Closure simulations in rough fractures

To investigate the evolution of fracture aperture fields considering normal stress and shear displacement, three-dimensional numerical models were established in COMSOL Multiphysics using the reconstructed fracture surface geometries. The simulations systematically calculated the spatial distribution of fracture apertures for combinations of six normal stress levels (σ n = 0.1, 5, 10, 20, 30, and 40 MPa) and six shear displacement levels (Δ l = 0, 2, 4, 6, 8, and 10 mm), resulting in a total of 144 simulation cases for the four fracture roughness groups. Statistical analyses of the aperture fields obtained under each loading condition were conducted to clarify the respective roles of normal stress and shear displacement in controlling fracture closure behavior. The resulting aperture fields were subsequently used as geometric inputs for the fracture flow simulations, thereby enabling quantitative evaluation of permeability evolution.

A series of FEM simulations was performed with the Solid Mechanics module (Huo et al., 2024). The rock matrix was assumed to behave as a linear elastic material, and the corresponding mechanical parameters are listed in Table 1. Contact pairs were defined between the upper and lower fracture surfaces to prevent nonphysical interpenetration and to realistically simulate the closure process. The computational domain was discretized using four-node tetrahedral elements, which are well suited for complex irregular geometries. To match the resolution of the scanned fracture models, the minimum element size in the vicinity of the fracture surfaces was set to 0.2 mm. In addition, an adaptive mesh refinement strategy was employed to achieve a balance between computational accuracy and solution efficiency (Huo and Meguid, 2025), as shown in Figure 4b.

As shown in Figure 4a, the lower domain was fixed to eliminate rigid motion, while displacement boundary conditions were imposed on its side surfaces. Preliminary simulations showed that direct application of stress loads on the top boundary often caused convergence difficulties in the contact analysis. Therefore, a displacement-controlled loading method was adopted to achieve the target normal stress indirectly. Specifically, incremental vertical displacements were applied to the top boundary, and after each loading increment reached static equilibrium, the average reaction force on the top surface was extracted. The corresponding equivalent normal stress was then calculated by dividing the reaction force by the loaded area.

When the equivalent normal stress reached one of the prescribed target values (σ n = 0.1, 5, 10, 20, 30, or 40 MPa), the normal gap distribution between the upper and lower contact surfaces was extracted and defined as the fracture aperture field for that condition. For different shear displacement conditions, the fracture geometries generated previously with offsets of 0, 2, 4, 6, 8, and 10 mm were directly imported into the model. The same displacement-controlled normal loading procedure and aperture extraction process were then repeated for each shear configuration. Consequently, a complete set of fracture aperture fields was obtained for all combinations of fracture roughness, shear displacement, and normal stress, providing the geometric basis for the subsequent flow simulations.

Flow simulations in rough fractures

Fluid flow through rock fractures has been extensively investigated under steady-state conditions (Happel and John, 1959; Crandall et al., 2010). Under the present flow conditions, inertial effects are negligible relative to viscous forces, and gravitational influences on the local fracture-scale flow field are assumed to be insignificant. Because the simulations were conducted at constant temperature, the injected water was treated as a Newtonian incompressible fluid with constant viscosity. Under these assumptions and within the laminar flow regime, Darcy’s law can be applied. Consequently, the fracture was represented using the parallel-plate approximation, for which the fracture permeability k f (m 2) is related to the hydraulic aperture D h (m) by:

According to Darcy’s law, the flow velocity can be expressed by Equation 6. where K f is the hydraulic conductivity (m/s), h is the hydraulic head (m), and l is the flow path length (m).

Considering mass conservation, the volumetric flow rate through the fracture must equal the injected flow rate Q, which can be expressed by Equation 7. where W is the fracture width. Thus, K f can be expressed as: where h 1 and h 2 are the hydraulic heads at the inlet and outlet, respectively.

The relationship between permeability and hydraulic conductivity is: where ρ is the water density, g is the gravitational acceleration, and μ is the dynamic viscosity of water. In this study, these parameters are taken as 1,000 kg/m 3, 9.8 m/s 2, and 0.001 Pas, respectively.

Combining Equations 8, 9, the fracture permeability can be expressed as:

Based on the above equations, the flow field was solved in COMSOL Multiphysics. After the fracture aperture field had been obtained for each stress state, fluid flow through the corresponding fracture geometry was simulated by solving the Laplace equation on the three-dimensional fracture domain. The hydraulic boundary conditions are illustrated in Figure 4a. Constant hydraulic heads of 10 m and 0 m were prescribed at the inlet and outlet boundaries, respectively, while all remaining fracture-wall boundaries were assigned zero-flux conditions (∇ H · n = 0) (Zhang and Nemcik, 2013; Develi and Babadagli, 2015; Chen et al., 2017; Zou et al., 2017a; 2017b). After solving for the steady-state velocity field, the normal flux across the outlet boundary was integrated to determine the total discharge. The equivalent fracture permeability under each loading condition was then calculated using Equation 10.

Fracture aperture evolution

Relationship between equivalent normal stress and normal displacement

In this study, a displacement-controlled loading method was adopted to achieve the target normal stress indirectly, thereby improving the numerical stability of the contact analysis. Figure 5 shows the relationship between equivalent normal stress and normal displacement for fractures with different JRC and shear displacements.

The results show that the equivalent normal stress increases monotonically with increasing normal displacement for all fracture groups. However, the stress-displacement curves exhibit two distinct loading stages. (1) Nonlinear compaction stage (σ n < 5 MPa): σ n increases slowly with normal displacement. During this stage, surface asperities progressively come into contact, leading to a continuous increase in the real contact area between the opposing fracture surfaces. Because contact is initially concentrated at isolated asperity peaks, the overall normal stiffness of the fracture remains relatively low. (2) Linear elastic closure stage (σ n > 5 MPa): σ n increases rapidly with normal displacement. At this stage, the major contact asperities have already formed stable contact zones, and the fracture exhibits a significantly higher and nearly constant normal stiffness.

The transition between these two stages reflects the progressive evolution of the fracture contact network. Under low normal stress, closure is dominated by asperity engagement and the redistribution of contact areas, whereas under higher normal stress, deformation is primarily controlled by the elastic compression of established contact regions. This behavior is consistent with the characteristic nonlinear closure response commonly observed in rough rock fractures (Luo et al., 2023).

Initial fracture aperture

To quantitatively characterize the initial contact state of the fracture, the initial fracture aperture D 0 was defined as the mean fracture aperture measured under σ n = 0.1 MPa. This low stress level approximately represents the reference state of the fracture prior to significant closure. Figure 6a illustrates the relationship between D 0 and shear displacement for fractures with different JRC.

The results show that, for all fracture groups, D 0 increases monotonically with shear displacement, while the growth rate gradually decreases as shear displacement increases. The relationship can be described by the following empirical equation: where a is a fitting parameter related to the geometric characteristics of the fracture. The coefficients of determination R 2 for all fitting curves exceed 0.90, indicating that the proposed relationship provides an excellent representation of the results.

Figure 6b further presents the relationship between the fitting parameter a and J R C. A strong linear correlation is observed and can be expressed as:

Combining Equations 11, 12, the initial fracture aperture can be expressed as a unified function of JRC and shear displacement:

This empirical formulation directly links D 0 to the JRC and shear displacement, providing a quantitative basis for the subsequent characterization of fluid flow and permeability evolution in rough rock fractures.

Aperture distribution characteristics and statistical analysis

The FEM simulations directly provide the contact stress distribution on the fracture surfaces for each loading condition. Figure 7 presents the contact stress contours for the fracture with JRC = 16.7 and Δ l = 4 mm under different normal stresses. The stress distribution exhibits pronounced spatial heterogeneity, with localized high-stress regions corresponding to asperity contact zones on the opposing fracture surfaces. As normal stress increases from 0.1 MPa to 40 MPa, both the contact area and the peak contact stress increase continuously. This behavior indicates that normal stress exerts a strong control on fracture closure and the evolution of asperity contact conditions (Cardona et al., 2021).

The corresponding spatial distributions of the fracture aperture field under the same conditions are shown in Figure 8. Similar to the contact stress field, the aperture field is highly heterogeneous. Low-aperture regions coincide with asperity contact locations, whereas regions with relatively large apertures are interconnected and form potential preferential flow channels. With increasing normal stress, the overall aperture decreases markedly, and the high-aperture channels become progressively compressed, sometime even completely closed.

To quantitatively evaluate the effects of normal stress, shear displacement, and surface roughness on D m, statistical analyses were conducted for all simulated aperture fields. Taking the case of JRC = 16.7 and Δ l = 4 mm as an example, Figure 9 presents the relative frequency distributions of fracture aperture under different normal stresses. The aperture distributions under all loading conditions can be approximated satisfactorily by normal distribution functions. The mean fracture aperture D m was obtained from the fitted normal distributions. As the normal stress increases, the distribution curves shift systematically toward smaller aperture values, indicating a monotonic reduction in D m.

Figure 10 summarizes the variation of D m with normal stress for all simulation cases. The results show that D m is negatively correlated with the normal stress, and positively correlated with the shear displacement. In addition, fractures with larger JRC values exhibit a more pronounced response to loading, indicating greater sensitivity of rougher fractures to both normal compression and shear-induced geometric changes.

The evolution of the fracture aperture with normal stress for all simulation cases can be uniformly described by: where D 0 is the initial aperture determined by Equation 13, K n is the normal stiffness of the fracture (taken as 17.5 GPa/m in this study). The average value of R 2 exceeds 0.8, indicating that Equation 14 can adequately describe the closure behavior of rough fractures considering normal stress and shear displacement.

Permeability evolution

Fluid flow within rough fractures is jointly controlled by fracture roughness, normal stress, and shear displacement (Sawayama et al., 2020). To illustrate the influence of these factors, Figure 11 presents the flow velocity field for the fracture with JRC = 16.7 and Δ l = 4 mm. As the normal stress increases from 0.1 MPa to 40 MPa, D m decreases from 1.80 mm to 0.88 mm, while the maximum flow velocity decreases from 0.76 m/s to 0.25 m/s.

The reduction in flow velocity is primarily caused by stress-induced fracture closure. As the fracture closes, the contact area between opposing asperities increases and the flow paths become more tortuous. Consequently, the effective flow channels are progressively constricted and distorted, resulting in increased hydraulic resistance and a reduction in the overall flow capacity (Chen et al., 2021b; Egert et al., 2021). Similar heterogeneous pressure distributions have been observed in laboratory experiments, numerical simulations, and field investigations (Yeo et al., 1998; Cappa et al., 2019; Ji et al., 2020).

Evolution of fracture permeability

The variation of equivalent fracture permeability k f with normal stress was calculated for different JRC and shear displacements, as shown in Figure 12.

The main observations can be summarized as follows: For all simulation cases, k f decreases monotonically with increasing normal stress, and the rate of decrease gradually diminishes. This trend is consistent with the nonlinear fracture closure behavior described in Section 3.1.

Under a given normal stress, larger shear displacement results in higher permeability, indicating that the aperture enlargement caused by shear offset significantly enhances fracture transmissivity.

Fractures with higher JRC values exhibit greater sensitivity to shear displacement. As shear displacement increases, the permeability enhancement becomes more pronounced for rougher fractures.

To quantitatively characterize the combined effects of fracture roughness, normal stress, and shear displacement, the permeability curves for all simulation cases were fitted. k f can be expressed as:

The average coefficient of determination for all fitting cases exceeds 0.80, demonstrating that Equation 15 adequately captures the permeability evolution of rough fractures considering normal stress and shear displacement.

Relationship between hydraulic aperture and fracture aperture

The hydraulic aperture D h, calculated from Equation 5, exhibits trends consistent with those of k f, as shown in Figure 13. Specifically, D h decreases with increasing normal stress and increases with both shear displacement and fracture roughness.

Figure 14 presents the relationship between D h and D m. Across all simulation cases, a strong positive linear correlation is observed, with a correlation coefficient exceeding 0.90. This linear relationship remains valid for different combinations of roughness, normal stress, and shear displacement.

These results indicate that, within the parameter range considered in this study, fracture roughness, shear displacement, and contact area primarily affect D h indirectly through their influence on D 0. Meanwhile, D h and D m maintain an approximately linear proportional relationship. This finding has important practical implications, as D h can be estimated directly from its geometric aperture using a linear scaling relationship, thereby substantially simplifying permeability prediction when only geometric measurements are available.

The simulation results demonstrate that fracture closure and permeability evolution in rough granite fractures are governed by the coupled effects of fracture surface roughness, normal stress, and shear displacement. This section discusses the individual contributions and underlying mechanisms of these controlling factors, followed by an evaluation of the engineering implications and the limitations of the present study.

Effect of fracture roughness on permeability evolution

Fracture surface roughness plays a dominant role in controlling permeability evolution considering normal stress and shear displacement. As shown Figure 12, under identical loading conditions, fractures with higher JRC consistently exhibit higher permeability. This behavior is fundamentally associated with the influence of roughness on the initial contact state and aperture distribution. Rougher fractures possess larger asperities that generate greater void spaces between opposing fracture surfaces after shear displacement. These enlarged voids form preferential flow channels, thereby reducing hydraulic resistance and enhancing permeability.

The results further indicate that rougher fractures are more sensitive to variations in both normal stress and shear displacement. As shown in Figure 10, the reduction in aperture with increasing normal stress is more pronounced for fractures with high J R C values, particularly during the nonlinear compaction stage (σ n < 5 MPa). This increased sensitivity arises because rougher fractures have a wider initial aperture distribution, providing greater potential for stress-induced modification of the flow path geometry. Similarly, the permeability enhancement caused by shear displacement becomes increasingly significant with increasing roughness (Figure 12), as larger asperities generate stronger dilatancy effects during shearing. This observation agrees with the findings of Fang et al. (2018) and Li et al. (2024), who reported that fracture roughness is the primary factor governing permeability changes during shear.

However, the influence of roughness is strongly modulated by the applied normal stress. Under high normal stress, progressive closure of large apertures reduces the hydraulic contrast between fractures with different roughness levels, thereby suppressing the sensitivity of permeability to JRC. Similar stress-dependent roughness effects have been reported by Sawayama et al. (2020).

Effect of normal stress on permeability evolution

Normal stress exerts a first-order control on fracture permeability by governing fracture closure and reducing the available flow pathways. For all simulation cases, permeability decreases monotonically with increasing normal stress (Figure 12), consistent with the closure behavior described in Section 3.1.3. As the normal stress increases from 0.1 MPa to 40 MPa, D 0 decreases by approximately 50%–60%, depending on fracture roughness and shear displacement, and the corresponding permeability reduction exhibits a strong nonlinear trend.

The two-stage closure behavior identified in Section 3.1.1 has direct implications for permeability evolution. During the nonlinear compaction stage (σ n < 5 MPa), asperity contacts develop rapidly and the contact area increases substantially. Even modest increases in contact area can significantly disrupt the connectivity of preferential flow channels, leading to a sharp decline in permeability. The contact stress distributions (Figure 7) show that stress concentrations form at discrete asperity contacts, which act as hydraulic barriers that partition the flow field into isolated channels. Once the fracture enters the linear elastic closure stage (σ n > 5 MPa), the primary contact framework has already been established. Additional loading mainly causes elastic compression of existing contact regions, and the rate of permeability reduction becomes progressively smaller. This transition is consistent with the observations of and Huo and Meguid (2025), who reported that the most significant permeability reduction in granite fractures occurs during the initial loading stage.

Effect of shear displacement on permeability evolution

Shear displacement has a significant positive effect on fracture permeability by altering the matching condition between the upper and lower fracture surfaces. As shown in Figure 12, k f increases with increasing shear displacement under a given normal stress and roughness condition. This enhancement is primarily driven by shear-induced dilatancy. As the fracture surfaces slide relative to each other, interlocked asperities become misaligned, creating additional void space and enlarging the effective aperture. The relationship between D 0 and shear displacement quantified in Equation 13 directly reflects this mechanism.

The magnitude of permeability enhancement depends strongly on fracture roughness. For fractures with JRC = 16.7, the permeability increases from Δ l = 0 mm to Δ l = 10 mm is substantially larger than that for fractures with JRC = 10.5 under the same normal stress. Rougher surfaces contain larger asperities that generate greater vertical dilation for a given horizontal offset. The resulting increase in aperture not only raises the mean aperture but also promotes the formation of interconnected high-permeability pathways. This interpretation is consistent with the findings of Li et al. (2024), who showed that the variance of the aperture distribution is an important contributor to permeability enhancement during shearing.

Nevertheless, the permeability enhancement induced by shear displacement is progressively suppressed by increasing normal stress. As shown in Figure 12, the separation between permeability curves corresponding to different shear displacements decreases at higher normal stresses. This indicates that compressive closure counteracts shear-induced dilation by reducing the enlarged apertures. The net permeability therefore reflects a dynamic balance between shear-induced dilation and stress-induced closure.

Implications for engineering

The findings of this study have important implications for underground engineering applications involving fluid flow through fractured rock masses, particularly enhanced geothermal systems (EGS), geological disposal of radioactive waste, and underground energy storage.

For EGS reservoir stimulation, fracture roughness, in-situ stress, and shear displacement should be carefully characterized and monitored because they strongly influence permeability evolution. The empirical relationships developed in this study—namely, the relationship between D 0 and JRC and shear displacement, and the relationship between k f and normal stress and D 0 —provide practical tools for estimating the hydraulic performance of stimulated fractures. In particular, the strong linear correlation between D m and D h enables permeability to be estimated directly from geometric measurements obtained through borehole imaging or core logging. In addition, the pronounced dependence of permeability on shear displacement indicates that hydro-shearing or natural fault reactivation can effectively enhance reservoir permeability, provided that the associated induced-seismicity risk is properly managed (Jia et al., 2022; Ji et al., 2022). These results suggest that integrating fracture roughness characterization with in-situ stress evaluation can improve the prediction of permeability enhancement during reservoir stimulation and support the optimization of injection strategies in EGS projects.

Limitations

Although the numerical model developed in this study successfully captures the closure and permeability evolution of rough fractures considering normal stress and shear displacement, several simplifying assumptions may introduce systematic deviations when extrapolating the results to field-scale conditions. These limitations are acknowledged and discussed below to provide appropriate context for the application of the proposed models.

First, the granite matrix was modeled as a linear elastic material, which neglects localized plastic deformation, asperity yielding, and micro-crushing at contact points under high normal stress. These inelastic processes can significantly influence fracture closure and permeability reduction (Nguyen and Selvadurai, 1998; Selvadurai and Nguyen, 1999; Chen et al., 2021a; Meng et al., 2026). Consequently, the present model tends to underestimate fracture closure and correspondingly overestimate shear-induced dilatancy and permeability enhancement, particularly at elevated stress levels. A rigorous quantification of the exact error magnitude would require direct experimental calibration under equivalent stress paths or the implementation of elastoplastic constitutive models, which we identify as a critical direction for future refinement.

Second, the flow simulations in this study are based on the cubic law and the Darcy flow assumption, which are applicable to laminar flow at low Reynolds numbers. According to Equation 14, the fracture aperture is jointly controlled by fracture roughness, normal stress, and shear displacement. When the aperture is relatively large, the average flow velocity is low, and the corresponding Reynolds number remains within the laminar Darcy flow regime (Re < 10). However, when the aperture is small, localized nonlinear flow may occur (Re > 10), where the flow regime can transition between laminar and turbulent, and inertial effects may become non-negligible. Under such conditions, a nonlinear flow model, such as the Forchheimer model, would be more appropriate (Chen et al., 2015; Rong et al., 2020; Zhong et al., 2023). Future investigations of high-rate injection scenarios should employ the full Navier-Stokes equations to capture potential inertial effects.

Third, this study conducted numerical simulations of fracture closure and permeability evolution based on a single rough fracture and derived an empirical model for permeability evolution. However, for engineering applications such as EGS and underground energy storage, the empirical permeability model cannot directly guide the optimization of field-scale hydraulic fracturing and hydro-shearing parameters. Future work will incorporate rock samples with varying sizes and properties to further refine and validate the empirical model. Subsequently, the permeability evolution law derived from a single fracture will be extended to complex fracture networks to improve the accuracy of permeability predictions.

Despite these limitations, the methodology and empirical relationships developed in this study provide a useful framework for understanding and predicting permeability evolution in rough fractures. The identified limitations also offer clear directions for future research on the hydro-mechanical behavior of fractured rock masses.

In this study, the surface morphology of tensile-induced granite fractures was obtained through high-resolution three-dimensional laser scanning, and numerical models of rough granite fractures with four roughness levels (JRC = 10.5, 12.2, 13.6, and 16.7) were established. A systematic numerical investigation was conducted to examine fracture closure and permeability evolution considering normal stress and shear displacement. Using COMSOL Multiphysics, FEM simulations were performed for six normal stress levels (σ n = 0.1, 5, 10, 20, 30, and 40 MPa) and six shear displacement levels (Δ l = 0, 2, 4, 6, 8, and 10 mm), enabling quantitative evaluation of the controlling effects of fracture roughness, normal stress, and shear displacement on aperture evolution and permeability. The principal conclusions are summarized as follows: A quantitative relationship between the initial fracture aperture, fracture roughness, and shear displacement was established. The initial aperture increases with both JRC and shear displacement. The empirical model provides a practical means of estimating the initial aperture of rough fractures.

Fracture apertures exhibit pronounced spatial heterogeneity due to stress-induced asperity contact and redistribution of void space. The aperture frequency distributions can be satisfactorily described by normal distribution functions. The mean fracture aperture decreases nonlinearly with increasing normal stress, while fracture roughness and shear displacement influence aperture evolution primarily through their effects on the initial aperture.

Fracture roughness, normal stress, and shear displacement affect permeability mainly by modifying the fracture aperture. Permeability decreases with increasing normal stress and increases with both shear displacement and fracture roughness. An empirical permeability model incorporating JRC, normal stress, and shear displacement was developed, providing a parametric framework for predicting permeability evolution in fractured rock masses.

The hydraulic aperture and fracture aperture exhibit a strong linear correlation across all simulation cases, independent of roughness, normal stress, and shear displacement. This indicates that within the investigated parameter range, these factors primarily influence permeability indirectly through changes in the initial aperture. The identified linear scaling relationship offers a convenient approach for estimating fracture permeability directly from geometric aperture measurements.