Firstly, the contact ratio of tight granite fractures was measured to prepare data for test sample preparation. The density of the granite specimen is 2601.67 kg/m³, the porosity is 1.183%, and the matrix permeability is 0.336 mD. Intact granite fractures are created with rough surfaces resembling the natural state using the Brazilian splitting method. Compared to computed tomography (CT) scanning, 3D laser, and reverse mode methods, the pressure-sensitive film offers a more convenient and cost-effective testing solution. The measurement principle relies on the dyeing characteristics of pressure-sensitive film, and its applicability has been validated by Nemoto et al. (2009) and Li et al. (2023).

Figure 2 illustrates the measurement process of the fracture contact ratio. The pressure-sensitive film consists of microcapsules, a color-developing layer (A-film), and a base layer (B-film), and is placed on the fracture surface. Initially, a mechanical device applies a normal load to the fractures, maintaining the target pressure for 2 min. When the pressure on the microcapsules exceeds the rated threshold, they rupture, releasing a color-developing material that reacts with the agent in the color-developing layer, forming red spots of varying brightness. The red image is then converted into grayscale, and the frequency characteristics of the grayscale values are analyzed.

Finally, a grayscale threshold is defined following the method proposed by Choi et al. (2019) and Nemoto et al. (2009). Grayscale values exceeding the threshold are defined as 1 (non-contact pixels), while values below the threshold are defined as 0 (contact pixels). Consequently, the grayscale image is converted into a binary image. As a result, the fracture contact ratio δ (%) under varying normal stress conditions can be calculated using the following formula: δ = ∑ N p x = 0 ∑ N p x × 100 % (1) where N p x = 0 represents the number of contact pixels, N p x represents the total number of pixels in the image.

The normal stress ( σ n ) is set at 5, 10, 15, 20, and 30 MPa. Without losing generality, samples of three different sizes were measured. Each load is measured three times independently to ensure result reliability. Table 1 summarizes the measurement scheme for the contact ratio.

Figure 3a clearly illustrates the spatial features of the interlaced distribution of contact domains and voids within the fractures. The similar spatial structures observed at different normal pressures confirm the repeatability and reliability of the pressure-sensitive film method. Figure 3b shows the relationship between the contact ratio and normal stress, revealing a positive correlation. The average contact ratios at normal stresses of 5, 10, 15, 20, and 30 MPa are 16.7%, 32.8%, 44.5%, 52.3%, and 66.4%, respectively. These results are similar to those reported by Li et al. (2023), Song et al. (2021), and Watanabe et al. (2009).

It is important to emphasize that the distribution of the contact ratio under the same stress conditions is discrete, reflecting the strong anisotropy of natural fractures. The contact ratio data typically follow a normal distribution. The binary image clearly reveals the presence of contact asperities in rough fractures. These randomly distributed contact protrusions occupy the flow channel space, constraining the fluid flow paths and significantly affecting the hydraulic properties of the fracture.

The measurement results in Section 2 provide the data needed to prepare fracture specimens with varying contact ratios. A 304 stainless steel sheet, with a diameter of 6 mm and a thickness of 0.1 mm, is attached to the smooth fracture wall to simulate natural contact asperities. Smooth fractures were created by symmetrically cutting a cylindrical core with a diameter of 50 mm and a length of 100 mm. Note that we focus is solely on the role of contact asperities in seepage; the effects of variable aperture and roughness are excluded for the smooth fracture sample. The distribution of contact asperities on the fracture wall is determined using a random placement method to simulate the heterogeneous nature of natural fractures. The contact ratio (δ) is calculated from the number of discs ( N ) in the fracture using the following formula: δ = N · A c A f × 100 % (2) A c = π r c 2 (3) A f = L · D (4) where A c is the area of a single disc; r c is the radius of the disc, 6 mm; L and D are the length and width of the fracture, respectively; A f is the nominal area of the fracture, 500 mm².

Fracture specimens with varying contact ratios are obtained by adjusting the number of discs (N) according to Equation 2. Based on the granite contact ratio measurements (Figure 3), five contact ratio levels were established: 10%, 20%, 30%, 40%, and 50%. To minimize the uncertainty due to random placement, each contact ratio level was repeated 20 times, resulting in a total of 100 contact fracture specimens. Figures 4a–e show the distribution of contact asperities at different contact ratios. As the contact ratio increases, the flow channel narrows, hindering fluid passage through the fractures.

The process for preparing the fracture specimens is illustrated in Figure 5. First, digital geometries are generated using a random placement method. A laser cutter is then used to cut the stainless steel sheets according to the digital geometry. The cut sheets (50 mm × 100 mm) serve as molds to position the contact asperities, while the cut discs simulate the asperities. The stainless steel mold is carefully aligned with the flat fracture surfaces, and the discs are placed accordingly. Once the target contact ratio is achieved, the mold is removed, leaving the discs to form an uneven arrangement of contact asperities. Finally, the two fracture halves are rejoined, aligned, and secured with heat shrink tubing to prevent the discs from shifting. This process is repeated to prepare fracture specimens for all contact ratios.

Single-phase seepage measurements were conducted on samples S1-S5 using a custom core seepage testing system, as illustrated in Figure 6. The system comprises an injection pump, a core holder, and a data acquisition module. The injection pump is an HXH-100B dual-cylinder model with constant flow and pressure capabilities, providing a maximum flow rate of 29.7 mL/min and an accuracy of 0.3%. In constant pressure mode, the pump applies a specified confining pressure to the core holder, while in constant flow mode, it maintains a consistent flow rate through the sample. Fluid pressure at the inlet is monitored by a PT131 pressure transmitter, which has a range of 0–1 MPa and an accuracy of 0.1%. The core holder accommodates cylindrical rock samples with a diameter of 50 mm and a length ranging from 50 to 200 mm, and is designed to withstand a maximum confining pressure of 30 MPa. The injection pump injects distilled water into horizontally placed, pre-fractured samples, while sensors measure the resulting flow pressure and rate.

Based on the capabilities of the injection pump and core holder, five confining pressure levels and six injection flow rates were set for each contact fracture specimen. Table 2 summarizes the experimental scheme. A total of 100 samples were tested. The main steps of the contact fracture seepage experiment are as follows:

1. Attach stainless steel discs to the fractures using contact asperities generated by a random placement method, and encapsulate the core with heat shrink tubing.

Assuming the fluid is viscous and incompressible, and the flow is steady and laminar, the volumetric flow rate Q f in a smooth fracture is positively correlated with the fracture aperture, as described by the cubic law: Q f = − D d h 3 12 μ w ∇ p (5) where D is the fracture width; μ w is the dynamic viscosity of water; d h is the hydraulic aperture of the fracture; ∇ p is the fluid pressure gradient: − ∇ p = p o u t − p i n L (6)

By measuring the Q f − ∇ p curve of the fracture, the equivalent hydraulic aperture d h obtained from the experiments can be determined using Equation 5: d h = 12 μ w Q f D ∇ p 3 (7)

Then, based on the cubic law, the permeability of a single smooth fracture can be calculated: k f = d h 2 12 (8)

Unfortunately, the fracture interface is irregular, with unevenly distributed asperities, which disrupts the permeability estimation theory based on the cubic law. Efforts have been made to correlate the hydraulic aperture d h with the mechanical aperture d f (Gao M. et al., 2022), roughness (Zheng et al., 2022), contact ratio δ (Zimmerman et al., 1992), or tortuosity (Tsang, 1984; Seybold et al., 2020) in order to correct the estimation bias of the standard cubic law. Specifically, Walsh (1981) proposed a fracture equivalent hydraulic aperture prediction model that does not rely on flow tests: d h d f 3 = 1 − δ 1 + δ (9)

Furthermore, considering the variability in the aperture of the fracture voids, Zimmerman and Bodvarsson (1996) proposed a prediction model based on contact ratio, mechanical aperture d f , and aperture variance s: d h d f 3 = 1 − 1.5 s 2 d f 2 1 − 2 δ (10) where 1 − 2 δ is considered as the correction term. Based on the differences in the correction term coefficient, Yeo (2001) proposed another prediction model: d h d f 3 = 1 − 1.5 s 2 d f 2 1 − 2.4 δ (11)

Since the experimental samples were designed without considering the rough characteristics of natural fractures, Equations 10, 11 can be simplified as follows: d h d f 3 = 1 − β δ (12) where β is defined as the correction factor for the contact ratio. Kirkpatrick (1973) proposed β = 2 based on the effective medium theory. When the contact ratio exceeds 50%, there is a deviation in the prediction of fracture flow.

We recorded the relationship between the volume flow rate and the pressure gradient for fractures with different contact ratios under varying confining pressures, as shown in Figure 7. The gray lines represent 20 samples from the randomly distributed asperity model, where each column corresponds to an increase in confining pressure (left to right), and each row corresponds to an increase in contact ratio (top to bottom). Figure 7 shows that under the experimental conditions, the volume flow rate and the water pressure gradient exhibit a linear positive correlation, indicating that fluid flow within the fractures follows darcy’s law and exhibits laminar flow behavior. The water pressure gradient ranges from 0 to 4 MPa/m, with the majority concentrated between 0 and 2.5 MPa/m. For fractures with a contact ratio above 20%, the water pressure gradient is consistently below 1 MPa/m.

When the contact ratio is below 20%, as the confining pressure increases, slight compressive deformation reduces the fracture channels, causing the water pressure gradient to rise and its slope to steepen. Even with a consistent contact ratio, the varying distribution of contact asperities alters the water pressure gradient, reflecting the resistance encountered by the fluid, as shown by the scattered data points in Figure 7. In scenarios with low contact rates (e.g., 10%, 20%), increased confining pressure results in greater dispersion of the measured ∇ p -Q data, indicating that confining pressure is the primary factor driving the differential flow of fracture elements in the system. However, in high contact rate scenarios (e.g., 50%), the dispersion of ∇ p -Q data is less affected by confining pressure, with the differential flow in the fracture system primarily driven by the distribution of contact asperities.

The average water pressure gradient from 20 realizations was used to assess the effect of the contact ratio on the water pressure gradient, as shown in Figure 8. As the contact ratio increases, the average water pressure gradient initially decreases and then increases. When the confining pressure increases from 5 MPa to 25 MPa, the decrease in the average water pressure gradient during the declining phase becomes more pronounced, with a marked increase in the slope of the curve. In contrast, the increase during the rising phase is less significant. The relationship between the average water pressure gradient and contact ratio is nonlinear in both phases. Specifically, during the phase where ∇ p decreases with δ , the decay gradient of ∇ p increases with both flow rate and confining pressure. In contrast, during the phase where ∇ p increases with δ , the decay gradient of ∇ p increases with the flow rate, but its sensitivity to confining pressure is not significant. Overall, the water pressure gradient’s response to the contact ratio becomes more pronounced.

Based on flow rate and pressure gradient data, Figure 9 illustrates the evolution of the equivalent hydraulic aperture with contact ratio under different confining pressures. As the contact ratio increases, the equivalent hydraulic aperture initially increases and then decreases, which aligns with the pressure gradient trend shown in Figure 8. Notably, within the contact ratio range of 10%–50%, the equivalent hydraulic aperture of the fractures ranges from 30 to 80 μm. This is significantly larger than the test results for straight fractures without contact asperities. A possible reason for this is that the contact asperities maintain voids in the fracture channels through which fluid can flow, reducing the closure effect of the fractures under confining pressure.

Figure 10 presents the evolution of the equivalent permeability of fractures with confining pressure under different contact ratios. As confining pressure increases, the equivalent permeability of the fractures gradually decreases and stabilizes, consistent with previous studies (Chen et al., 2017). Notably, at higher contact ratios (e.g., S4, S5), the influence of confining pressure on fracture permeability reduction becomes negligible. According to Tsang and Witherspoon (1981), and Zimmerman et al. (1992), when δ < 25%, the position of the contact body has a negligible effect on the fracture permeability. However, Figure 10 shows that the permeability data for lower contact ratios are more scattered, suggesting that fracture permeability at low contact ratios is more significantly influenced by the contact position than at higher contact ratios.

To demonstrate the impact of contact ratio on the permeability of fractured pores, Figure 11 shows the trend of the arithmetic mean permeability, calculated from 20 realizations, as it evolves with confining pressure. It is well documented that fractures permeability follows an exponential decay as confining pressure increases. Therefore, the following exponential function was used to model the relationship between permeability and confining pressure (Rutqvist et al., 2002; Lei et al., 2021): k f = a + b e c σ c (13) where a , b , and c are parameters related to the contact ratio.

Figure 11 shows that fracture permeability in scenarios with lower contact ratios (e.g., 10%) is more sensitive to changes in confining pressure compared to higher contact ratios (e.g., 40%–50%). Specifically, as the confining pressure increases from 5 MPa to 25 MPa, the average permeability of fractures at a 10% contact ratio decreases from 4.07 × 10⁻¹⁰ m² to 1.77 × 10⁻¹⁰ m² (a 56.5% reduction), while the average permeability at a 50% contact ratio case decreases from 2.45 × 10⁻¹⁰ m² to 2.08 × 10⁻¹⁰ m² (a 17.7% reduction). Table 3 lists the fitting results of Equation 13, with high R ² values further confirming the exponential relationship between permeability and confining pressure.

Traditional contact-based equivalent hydraulic aperture prediction models have been widely proposed (Kirkpatrick, 1973; Walsh, 1981; Zimmerman and Bodvarsson, 1996; Yeo, 2001). However, their applicability under compression conditions requires further examination. Figure 12a shows the relationship between d h / d f 3 and contact ratio. The gray lines represent theoretical data from the equivalent hydraulic aperture prediction model (Equations 9, 12), with the correction factor β in Equation 12 set to 0.5, 1, 1.5, 2, and 2.4. Additionally, we measured the seepage data for fractures with regular contact asperities, which are also shown in Figure 12a. The results indicate a consistent negative correlation between d h / d f 3 from the prediction model and the contact ratio. This is expected, as the prediction model does not account for the compressive deformation effects induced by confining pressure on the flow channels. Therefore, the primary limitation of such prediction models is their neglect of stress. In contrast to classical prediction models, the experimental data in this study show a strong negative correlation under higher contact ratio conditions ( δ > 20%), with the correction factors β (R ² > 0.7) taking values of 1.82, 1.88, 1.91, 1.94, and 1.96 at 5 MPa, 10 MPa, 15 MPa, 20 MPa, and 25 MPa, respectively.

Interestingly, under lower contact ratio scenarios, d h / d f 3 exhibits a positive correlation with the contact ratio and no longer follows the evolutionary trend of traditional prediction models. This anomalous inversion effect was also observed in fractures with regular contact asperities, although it was less pronounced in these cases. Figure 12b further illustrates how d h / d f 3 evolves with the contact ratio. It shows that, for all fracture specimens, d h / d f 3 initially increases and then decreases as the contact ratio increases. This inversion effect is controlled by confining pressure. For example, compared to the confining pressure of 25 MPa, the phenomenon of d h / d f 3 exhibiting an inversion effect is not significant under a confining pressure of 5 MPa. Additionally, there is a critical contact ratio, δ c r i , in the contact ratio range of 20%–30%. This implies that above the critical contact ratio, the traditional prediction model remains applicable; otherwise, it is no longer valid.

Figure 11 confirms that fracture permeability is sensitive to both contact ratio and confining pressure. Typical thermal-hydraulic-mechanical models in fractured rock masses account only for the stress effect, neglecting the uneven compression process caused by contact ratio, including the inversion effect shown in Figure 12. Consequently, a new fracture permeability model that incorporates both contact ratio and stress is needed. Table 4 presents the exponential relationship between permeability and confining pressure, as described by Equation 13. The coefficients a , b , and c are clearly related to the contact ratio. Therefore, we further examined how these coefficients depend on the contact ratio, with the results shown in Figure 13. Under the dual influence of contact ratio and stress, coefficients a and c exhibit reversal phenomena, whereas coefficient b does not show such reversal. Thus, a staged linear formula was applied to fit the anomaly of the coefficient-contact ratio relationship before and after the critical contact ratio δ c r i , and the results are listed in Table 4. Here, cases below the critical contact ratio are defined as low contact ratio δ l o w , while those above are defined as high contact ratio δ h i g h . It should be noted that the exact value of the critical contact ratio δ c r i is unknown in this study, so the data for δ = 30 % was used for linear fitting in both stages.

The fitting results allow us to develop an empirical prediction model for contact-type fracture permeability, which simultaneously considers contact ratio and confining pressure, as follows: k f = a 1 + a 2 δ + b 1 + b 2 δ e c 1 + c 2 δ σ c (14) where, Table 4 presents the values of a i , b i , and c i ( i = 1 , 2 ) under different applicable conditions. When the effect of confining pressure is not considered, i.e., σ c = 0, the prediction model can be simplified to the generalized contact ratio correction model, which is consistent with Equation 12.

Using the prediction model, Figure 14 compares the measured permeability with the predicted permeability for all test cases, including those with regular contact asperity fractures. The linear regression fit (red line) for all data is used to evaluate the model’s performance. It is observed that the linear fit closely matches the 1:1 line. The root mean square error (RMSE) between the observed log₁₀(k _f ) and the predicted log₁₀(k _f ) is 0.057, indicating that the predictive model performs well. Additionally, nearly all the data points are distributed within the 95% prediction interval. These results suggest that the developed prediction model accurately captures the relationship between permeability, confining pressure, and contact ratio, demonstrating high reliability. However, it can also be observed that some of the discrete data points deviate significantly from the 1:1 line, indicating the complexity of permeability in fractured rock masses. For regular contact asperities, the predicted permeability is generally underestimated. The model performs better for cases with higher contact ratios compared to those with lower contact ratios.