On the Relationship Between Normal Stiffness and Permeability of Rock Fractures

The hydraulic permeability and normal stiffness of rock fractures are fundamentally related, which can serve as a probe for assessing hydro‐mechanical properties in a variety of geo‐engineering contexts. We present experimentally validated numerical simulations of fracture closure and fluid flow processes in realistic aperture structures with mean in the range of 0.01–2 mm. A relatively simple relationship between permeability and normal stiffness is derived in the effective medium regime, using the elastic modulus of the rock matrix, the mean mechanical aperture, and the relative standard deviation (RSD) of the local mechanical aperture. The established relationship between stiffness and permeability appears independent of scale or orientation.

LI ET AL. 10.1029/2021GL095593 2 of 9 were experimental data compiled by Pyrak-Nolte and Morris (2000), generally dominated by low-conducting fractures. Hence, the currently available relationships between fracture stiffness and permeability (Pyrak-Nolte & Nolte, 2016;Wang & Cardenas, 2016) are derived primarily for tight fractures with relatively low mean apertures. Specifically, mean apertures considered by Pyrak-Nolte and Nolte (2016) are below 10 μm, whereas the mean aperture for the Gaussian model considered by Wang and Cardenas (2016) is 5 μm.
In spite of significant advances in understanding stiffness and permeability of fractures following the works of Petrovitch et al. (2013), Pyrak-Nolte and Nolte (2016), and Wang and Cardenas (2016), there are several important knowledge gaps that need to be addressed. Most fractures for instance in the Scandinavian bedrock considered in engineering applications from tunneling to geological disposal, are of transmissivities in the range 10 −9 to 10 −5 m 2 /s (Ludvigson et al., 2002;Zou & Cvetkovic, 2021), which based on the cubic law corresponds to mean apertures of 0.01 and 0.23 mm respectively. Rock fractures are generated by shear displacement and modified by geochemical processes over time. While aperture genesis has yet to be fully understood, there is a general interest to explore aperture structures beyond those generated by classical statistical models. Specifically, methods have been proposed where two compressed rough-walled natural fracture surfaces can generate a wide spectrum of correlations induced by various mechanical (e.g., slip) and chemical (e.g., erosion) processes (Brown & Scholz 1985, 1986, Pyrak-Nolte & Nolte, 2016. The resulting aperture fields subject to a normal stress are characterized by their unique physical, mechanical, and geometrical properties (Swan, 1983).
The stress concentration on asperities that are in contact in a typical fracture, could result in substantial damage on asperity tips as evidenced by hysteresis in cyclic loading tests and post-experiment surface measurements (Bandis et al., 1983;Brown & Scholz, 1986;. In previous studies of the stiffness-permeability relationship (Petrovitch et al., 2013;Wang & Cardenas, 2016), pure elastic contact models were employed; neglecting plastic deformation will typically underestimate the fracture closure and therefore overestimate the local apertures, which may consequently lead to a biased estimation of the stiffness-permeability relationship. Thus, there is a need to explore effects of combined elastic-plastic deformation on the stiffness and permeability of fractures.
The aim of this study is to address above mentioned knowledge gaps and to establish a stiffness-permeability relationship for fractures with realistic surfaces and with mean apertures typical for the effective medium regime, that is, 0.01-2 mm. An experimentally validated elastic-plastic contact model is employed for calculating the closure of three types of natural tensile fractures with a length scale between 0.08-0.6 m. The permeability of resulting aperture structures and generated Gaussian aperture fields along two orthogonal directions subject to a normal stress up to 100 MPa is estimated and validated by stress-flow tests. Finally, a relationship between permeability and fracture stiffness using measurable matrix/fracture properties, is presented and its physical basis are discussed.

Methodology
Three medium-grained granite samples (MG-1, 2, and 3), three fine-grained sandstone samples (FS-1, 2, and 3) and three medium-grained sandstone samples (MS-1, 2, and 3) with a dimension of 80 × 80 × 80 mm, a cylindrical granite sample (MG-4) with a diameter of 50 mm and a height of 100 mm, and a granite fracture surface sample (MG-5) with a dimension of 800 × 800 mm were prepared. All granite samples are of the same origin. The mechanical properties of tested rock samples were measured using a MTS 815 system (see Table S1 in Supporting Information S1), which were used in fracture closure calculations. The cubic samples and the cylindrical sample were axially split to create fresh tensile fractures. Surfaces of these fractures were optically scanned and digitized with a resolution of 0.1 μm, while the MG-5 was laser-scanned with a resolution of 25 μm (see Figure S1 in Supporting Information S1).
Schematic views of the normal loading and stress-flow experiments are shown in Figure S2 in Supporting Information S1. The experimental setup and procedure for the normal loading test follow the same protocols as described in . Samples MG-3, FS-3, and MS-3 were selected for normal loading tests after being artificially dislocated by 4 mm. To avoid splitting of the entire sample, the maximum normal stress for FS-3 and MS-3 was limited to 5 MPa while the maximum stress for MG-3 was 10 MPa.

of 9
The stress-flow tests employed a conventional core-holder and the fractures were dislocated by an increasing offset of 4, 6 and 8 mm along the axial direction. Under each offset, the confining pressure was increased at a rate of 0.5 MPa/min and temporally stopped at 2, 4, 6, 8, and 10 MPa to allow the implementation of fluid flow tests. A series of flow rates ranging from 1.0 to 150 ml/min was imposed on the inlet for each test using a syringe pump (500D, Teledyne, USA). The differential hydraulic pressure between the inlet and outlet was simultaneously measured by a differential pressure transmitter (EJX-110A, Yokogawa, Japan) with a resolution of 10 Pa.
The fracture closure subject to normal stress is calculated following the framework of the variational principle for contact mechanics (e.g., Li et al., 2015;Tian & Bhushan, 1996;. By assuming that the rock asperity deforms as an elastic-perfectly plastic material and the plastic deformation is confined within a local small area, an elastic-plastic contact model (EPM) was developed based on an elastic contact model (EM). In EPM, the dissipation of energy caused by plastic deformation is added to the total complementary potential energy (see Text S1 in Supporting Information S1 and Tian & Bhushan, 1996 for details); a parameter referred to as hardness of the rock matrix (H) is a threshold value for the local contact stress when plastic deformation occurs (i.e., when σ n > H). The value of H for a specific rock is typically linked to its elastic modulus E or unconfined compressive strength UCS (Kling et al., 2018). Here, H = 0.03E for MG, H = 0.025E for FS, and H = 0.02E for MS were adopted through calibration with experimental data, due to difficulty in physical measurement of the matrix hardness with surface roughness.
EM and EPM are both applied to the tested fracture samples to verify their applicability and for the purpose of validation. The remaining fractures were numerically subject to maximum normal loading of 100 MPa to create a series of aperture fields. Each loading process was solved in a stepwise scheme by applying an incremental normal stress of Δσ n = 0.1 MPa. The normal stiffness k n and the permeability k were calculated by (Pyrak-Nolte & Morris, 2000), where Δδ is the incremental normal displacement, and e h is the hydraulic aperture back-calculated from the measured or calculated differential pressure-flow rate data using the cubic law.
A pair of natural fracture surfaces can only provide a certain pattern of aperture structures. To study the effect of aperture structure on the relationship between stiffness and permeability, a series of unmated synthetic fractures were created by dislocating the originally mated natural fracture surfaces with different offset (displacement or dislocation) lengths l o (0.2, 0.3, 0.5, 1, 2, 4, 6 and 10 mm). These synthetic fractures are more realistic representation compared to randomly generated apertures, since they consist of rough surfaces with real surface roughness. In parallel, a series of surface samples with a side length ranging from 100 to 600 mm were extracted from the large-sized sample (MG-5) and two representative offsets l o (1 and 10 mm) were applied on this sample series. In the computations with the dislocated fractures, an incremental normal loading of 100 MPa was applied and the relative standard deviation (RSD) of the resultant aperture fields at each Δσ n was extracted. Gaussian aperture structures (Wang & Cardenas, 2016) for a 100 × 100 mm sample were also generated to explore the relationship more broadly between permeability and stiffness.
The hydraulic pressure difference between the inlet and outlet of each fracture was computed by solving the Reynolds equation, applying a constant flow rate to the inlet and non-permeable conditions to the lateral boundaries (Koyama et al., 2008;. Permeability k of each fracture along two orthogonal (x-and y-) directions was calculated.

Relating Fracture Stiffness and Permeability
To demonstrate the importance of plastic deformation, we considered simulations with EM and EPM to compare the closure behavior observed in normal loading experiments on three rock fractures (MG-3, FS-3, and MS-3) as shown in Figure 1a. All three cases demonstrate that the EPM well accommodates the experimental results, while the EM yields notably smaller deformation than that of the experiment. The contact area ratio C a increases almost linearly with the increasing normal stress σ n , and the EM produces smaller C a than the case with EPM. The EM treats the contacted asperities as elastic mediums without yielding, 4 of 9 and therefore the local stress acting on each contacted asperity increases monotonically with the increasing normal displacement δ. This results in unrealistically large local stresses and underestimated fracture closures (see Figure S3 in Supporting Information S1). The post-experiment measurements on fracture surfaces reveal significant asperity damages that are well accommodated by EPM (see Figure S4 in Supporting Information S1).  (Zimmerman & Bodvarsson, 1996;Zimmerman et al., 2004).
Using the validated EPM, we then calculated the closure of fractures of different scales subject to a series of offsets l o . All calculated cases indicate that δ increases with the increasing σ n at a declining rate, and asymptotically approaches some limiting value (see Figure S5 in Supporting Information S1). Comparisons among different cases reveal that the magnitude of compression-induced fracture closure depends on the deformability of asperities, and the available volume of internal voids and their distributions that are determined by the scale and the surface morphology of two fracture walls as well as their correlations.
In Figure 2, we show exemplified aperture fields representing void space distributions and flow fields along x-and y-directions for sample MG-1 with two different l o and under four different σ n . The degree of mismatch and the available void space increase as l o increases, resulting in expanded apertures ( Figure 2) and reduced fracture stiffness. The dispersed small contact spots when l o = 1 mm change to a few large, concentrated contact areas for l o = 10 mm, analogous to those observed in fractures during a shear process (Li et al., 2019). As σ n increases, both the number of contacting spots and their area increase, gradually dividing the fluid flow into a series of channels with enhanced tortuosity that is the main source of additional pressure losses reducing the ratio of e h /e m . The shape of contact spots is string-like with the long-axis oriented along the y-direction. Fluid flow in this case is subject to less resistance along the y-direction than that along the x-direction, resulting in anisotropic permeability. In general, the permeability decreases with the increasing normal stress for all cases (see Figure S6 in Supporting Information S1).
The elastic modulus E of the rock matrix or asperities is a controlling factor for normal deformation of fractures, and the RSD (i.e., σ e /e m , where σ e and e m are the standard deviation and the mean of mechanical aperture, respectively) of aperture structure is a direct measure for fluid flow (Wang & Cardenas, 2016). It is also noted that the permeability k scales with 2 h E e , suggesting that the e h ∼ e m interrelation reflects the fluid flow and aperture structure correlation. Following this line, we analyzed the calculated results and find the following relationship between the relative permeability k r and normal stiffness k n as ( where, m E e is the mean value of all calculated mechanical apertures that in this study is approximately 0.5 mm, and α is a dimensionless scaling factor that is found to be 0.003. In Equation 2, E/k n can be interpreted as a representative asperity positively correlated with e h , and its physical meaning is illustrated by Figure S7 in Supporting Information S1. In addition, the RSD is placed in the numerator because it is positively correlated with e h , and the m E e is placed in the denominator for normalization. 10.1029/2021GL095593 6 of 9 All the experimental and calculated data based on either dislocated natural fractures or generated Gaussian apertures generally follow this relationship which is applicable in the effective medium regime (Figure 3). The relationship (Equation 2) is dimensionally consistent and is based on measurable, physically valid parameters. In Figure 3, a critical point is identified for E k e n m   RSD/ 1 (denoted by the green symbol), quantifying the transition between the percolation and effective medium regimes as discussed by Petrovitch et al. (2013). In view of Equation 2, the slope α = 0.003 is obtained as the y-coordinate for the critical point ( Figure 3).

Discussion
One way to test the physical basis of the proposed relationship (Equation 2) is to explore alternative means for inferring the slope α ≈ 0.003 as inferred from Figure 3. Since the mean value of RSD for the investigated fractures is around 1 (see Data Set S1), we assume at first-order that RSD ≈ 1 in Equation 2 / / , we obtain α ≈ 0.003, which indicates that this factor has a physical basis and can be deduced from data generated by different approaches. A mechanical interpretation of the expression α = k n /(E/e h ) = 0.003 is as follows. E/ e h is the stiffness of the rock matrix defined as the stress change required to deform the rock matrix by e h . Equation 2 with α = 0.003 implies that the stiffness of a fracture k n with hydraulic aperture e h is roughly 300 times smaller compared to the stiffness of the rock matrix (E/e h ), which is applicable to fractures with a typical surface roughness. We calculated the transmissivity Q/Δh for all experimental and numerical cases where Q is the flow rate and Δh is the hydraulic head difference. In Figure 4, the correlation between Q/Δh and k n , along with experimental data from different sources (Pyrak-Nolte & Morris, 2000;Thörn et al., 2015) is shown. We identified a general descending trend of Q/Δh with the increasing k n over a range of three orders of magnitude. Each fracture exhibits a unique inclination representing the sensitivity of Q/Δh to k n . The sensitivity increases as the matedness of fractures decreases and becomes particularly remarkable when k n is greater than 10 3 GPa/m, indicating the transition from a percolation to an effective regime (Petrovitch et al., 2014). Data from other sources register a similar trend and range but variable patterns of inclination due to their independent combinations of mechanical and geometrical properties.
It has been shown that neglecting plastic deformation will generally underestimate fracture closure and therefore overestimate local apertures, which may be particularly relevant in the percolation regime where the change of permeability is most sensitive to change in stiffness. Our method for incorporating plastic deformation uses hardness H as the threshold, a parameter that is still difficult to infer independently as a material property. Exploring alternative and possibly more rigorous means for incorporating plastic deformation in the stiffness-permeability relationship of fractures, is one interesting direction for future research. Another interesting direction is to study how real aperture structures generated for instances by shear deformation, affect the stiffness-permeability relationship, compared to the dislocation and pure statistical aperture models used in the present work.

Conclusion
Using high-resolution simulations supported by experimental measurements, we established a new and relatively simple relationship (Equation 2) between fracture permeability k (or 2 h E e ) and fracture stiffness k n . The proposed relationship between stiffness and permeability can be used for first-order estimates of coupled hydro-mechanical properties in natural rock fractures relevant for many geo-engineering applications. Relationship (Equation 2) involves three physical parameters: modulus E that characterizes the deformability of the matrix, m E e as a representative mechanical aperture, and RSD to represent the relative spatial variation in local apertures. It represents linear scaling between two dimensionless parameter groups of physical