Relationship Between Permeability and Resistivity of Sheared Rock Fractures: The Role of Tortuosity and Flow Path Percolation

The fluid‐flow properties of fractures have received increasing attention regarding the role of geofluids in the genesis of slow and fast earthquakes and recent advances in geoengineering developments. Geophysical observations are promising tools to remotely estimate crustal permeability changes; however, quantitative interpretations are limited by the rock‐physical models' paucity for fractures. This study investigated changes in permeability, resistivity, and their respective relationships at elevated stress by performing numerical simulations of different fracture models with varying fracture size, roughness, and shear displacement. Numerical results and microscopic flow analysis demonstrate that permeability–resistivity relationships are controlled by percolation and are less dependent on fracture geometric characteristics. Our finding suggests that the permeability evolution of fractures can be formulated with resistivity changes independent of both fracture size and microstructure, the trends of which can be predicted using Archie's exponent. The extension to the electro‐mechanical relationship further derives the potential applications of estimating stress changes.

stress and the shear displacement of faults.Given that in situ stresses are never constant throughout geoengineering developments or on the geological time scale, these impacts must be considerable in natural settings even in the absence of major seismicity.
Although direct measurement of permeability changes is a challenging task, geophysical observation using electrical resistivity is a potential approach.In a wide range of fields, electromagnetic surveys have observed resistivity changes that are associated with subsurface stress changes, such as geothermal activity, fluid production, and seismicity (Aizawa et al., 2011;Didana et al., 2017;Johnson et al., 2021;Mazzella & Morrison, 1974;Park, 1991;Peacock et al., 2012Peacock et al., , 2013)).It would be advantageous if changes in permeability caused by variations in the subsurface stress could be linked to electrical resistivity, which could be monitored remotely.Several models have been proposed on the relation between permeability and resistivity of porous media (Bernabé et al., 2011;Cai et al., 2019;Ghanbarian & Male, 2021;Glover et al., 2006;Hilfer, 1992;Katz & Thompson, 1986;Revil & Cathles, 1999).Such models allow the estimation of permeability from resistivity data (Anderson et al., 1985;Becker et al., 1982;Katayama et al., 2020;Slagle & Goldberg, 2011;Yamaya et al., 2022).However, some field observations have detected unusual changes in permeability and resistivity (Park, 1991;Xue et al., 2013) that could not be explained by experimental findings due to the existence of a macroscopic fracture in the fields.More recent observations also evoke the importance of macroscopic fractures as they contribute to the anomalously high permeability in the underthrust, which may be related to the pore pressure supply and the generation of slow earthquakes (Heuler et al., 2017;Hirose et al., 2021).However, the paucity of data using fractured rock mass hinders attempts at updating the rock-physical models that can be used to evaluate changes in permeability.A pioneering study by Stesky (1986) reported that permeability and resistivity in saw-cut fractures are linearly related in a logarithmic coordinate.Conversely, numerical studies based on rough-walled fractures have reported a nonlinear relationship resulting from the heterogeneous flow distributions (Kirkby et al., 2016;Sawayama et al., 2021a), similar to the previously reported percolation for fracture permeability (Berkowitz, 1995;Pyrak-Nolte et al., 1988).However, how the fracture geometric characteristics affect the percolation and its impact on the permeability-resistivity relationship remains unclear.
Thus, this study employed a digital rock physics approach (Andrä et al., 2013;Sawayama et al., 2022) to determine permeability and resistivity simultaneously while assessing the evolution in local connections and percolation of preferential flow paths with the elevated normal stress and shear displacement.We discuss the role of local flow behavior on permeability, resistivity, and their respective relationships, and suggest potential practical applications of our findings.

Methodology
We conducted a series of numerical simulations of fluid flow and electrical conduction within deformed rock fractures.The fracture surface topography is known to be self-affine fractal and can be created numerically (Brown, 1995;Matsuki et al., 2006).Following these studies, synthetic fracture models with controlled roughness and shear displacement were generated based on the fractional Brownian motion (Text S1 in Supporting Information S1).To investigate the effect of the fractal characteristics, the fractal dimension D and surface roughness s (i.e., the root mean square height of the surface topography) were varied as D = 2.0, 2.2, 2.4 and s = 0.1, 0.2, 0.3 mm referring to natural faults (Brown & Scholz, 1985;Power & Durham, 1997;Power et al., 1987).Note that such fractal nature has been confirmed as valid for natural fault roughness from local to global scales (Candela & Brodsky, 2016;Candela et al., 2012;Renard & Candela, 2013).Based on these values, we prepared a synthetic fracture model of varying fracture size L = 24, 48, 96 mm with a 0.1-mm square grid, where a well-known scaling law for fracture roughness was incorporated (Matsuki et al., 2006): ) 3− . (1) Here, the subscript zero represents the reference model (L 0 = 24 mm, s 0 = 0.1, 0.2, 0.3 mm).For each fracture model, we numerically simulated the elastic-fully plastic deformation of the asperities using a half spacebased dry contact model with up to 100 MPa of normal loading (Polonsky & Keer, 1999).The corresponding algorithm is summarized in Text S2 in Supporting Information S1.The boundary was treated as periodic to eliminate the boundary effect (Petrovitch et al., 2014).We confirmed that the constructed aperture, that is, the gap of the deformed rough surface, follows the scaling law (Figure S2 in Supporting Information S1) observed in natural joints and faults (Kim & Sanderson, 2005;Schlische et al., 1996;Schultz et al., 2008).The variations in mean aperture value and the aperture distribution against stress are shown in Figure S3 in Supporting Information S1.Furthermore, additional 13,000 simulations (for a total of 14,530) were performed by varying the random seeds to investigate the stochastic fluctuations of our simulation results (Text S1 in Supporting Information S1).
We then simulated 3D fracture flows using the lattice Boltzmann method, where we adopted a multi-relaxation-time D3Q19 model (Ahrenholz et al., 2008;Jiang et al., 2014;Sawayama et al., 2021a).The provision of a constant body force from the inlet to the outlet boundaries and the periodic boundary along the fracture plane allowed us to simulate the fracture flows in the directions perpendicular to the shear direction (Figure S4 in Supporting Information S1).Then, the fracture permeability was estimated from the macroscopic flow velocity, assuming Darcy flow and negligible matrix permeability (∼10 −20 m 2 ).Subsequently, we analyzed resistivity using the finite element method (Garboczi, 1998;Sawayama et al., 2021b).Similar to the fluid flow, we simulated the local field of electric currents in the direction perpendicular to the shear direction, and calculated corresponding resistivity.In this analysis, fluid and solid conductivities of 5 and 10 −5 S/m, respectively, were modeled (Sawayama et al., 2018).We also evaluated the bulk permeability and resistivity of the fractured rock mass having a constant thickness (n z = 24 mm) to consider the porosity effect.This enables us to correlate the permeability-resistivity relationship within the same elementary volume.Notably, the relationship does not change regardless of the thickness value even though their absolute values change, when matrix permeability and resistivity are presumed to be negligible.Then, we calculated the formation factor F by dividing the bulk resistivity ρ b by the fluid resis tivity ρ w , (F = ρ b /ρ w ) such that the effects of both the temperature and salinity on fluid resistivity can be neglected.

Changes in Permeability and Resistivity as Elevated Normal Stress
We present the evolution of permeability, electrical resistivity, and the associated changes of the local flow paths of various fracture data sets as elevated normal stress (Figure 1).The results demonstrate the dependencies of fracture size (Figures 1a and 1b), shear displacement ratio (Figures 1c and 1d), and surface roughness characteristics (Figures 1e and 1f) on these properties.Permeability and resistivity of joints (i.e., no shear displacement) are less dependent on fracture size, as previous studies have reported (Ishibashi et al., 2015;Sawayama et al., 2021b).In addition, faults (i.e., sheared fractures) exhibit scale-dependent characteristics, where permeability increases and resistivity decreases remarkably with fracture size (Figures 1a and 1b).The stress-dependent changes in permeability and resistivity are significant in fractures with larger shear displacements (Figures 1c and 1d).The roughness increases permeability and decreases resistivity, and these properties are nearly constant under different fractal dimensions (Figures 1e and 1f).The higher roughness of the fracture surfaces (i.e., rms height s) increases the aperture because it enhances the mismatch between two surfaces, and thus leads to higher permeability and lower resistivity (Figure S5 in Supporting Information S1).The roughness effect will change depending on the shear displacement (Figure S5 in Supporting Information S1) and the two orthogonal directions (Figure S6 in Supporting Information S1).The effects of shear displacement and roughness become insignificant under higher stress conditions (Figures 1c-1f).
The panels at the bottom of Figure 1 show representative simulation results for the evolutions of fracture flow within the heterogeneous aperture distribution with and without shear displacement.Note that the flow rates have been vertically summed in every z-direction and the regions with >1% of the maximum is visualized such that the dominant paths in the fractures can be projected onto the x-y plane.At the lowest stresses, dominant flow paths induce channeling that covers most of the open space (Figure 1, images i and iv).As the stress increases, larger portions of the fracture surfaces are in contact; thus, the channeling paths decrease progressively and become narrow (Figure 1, images ii and v).Further normal loading diminishes the dominant flow paths.Accordingly, the flow paths from the inlet to the outlet lose their connection (Figure 1, images iii and vi).This phenomenon is nearly the same in both joints and faults.The stress level at which this transition occurs can be defined as the percolation threshold that considerably controls the transport behaviors of rocks (Berkowitz, 1995;Ghanbarian et al., 2014;Guéguen et al., 1997;Pyrak-Nolte et al., 1988).Above this percolation threshold, permeability, and resistivity demonstrate different trends against changes in stress.

Effect of Percolation and Tortuosity on Resistivity
It is known that resistivity has a power law relationship with porosity (Archie, 1942): where φ is porosity, and a and m are empirical coefficients.Previous studies pointed out that m, the so-called cementation exponent, is dependent on pore connectivity (Glover, 2009;Guéguen & Palciauskas, 1994).Based on the local fluid paths (Figure 1, images i-vi) in each condition, we examine whether the dominant flow paths (i.e., >1% of the maximum flow rate; Figure 1) are fully connected throughout the model using the image processing technique (i.e., bounding box).Consequently, we have classified all results into fully connected and partially disconnected flow regimes of the dominant paths (Figure S7 in Supporting Information S1). Figure 2a plots the formation factor and porosity relative to these flow regimes.Overall, our results in the fully connected flow regime demonstrate a linear trend in the log-log plot (Figure 2a), which is consistent with Archie's exponent m = 1.2.Conversely, the results of the partially disconnected flow regime deviated from the linear line.This nonArchie behavior suggests that the connectivity of the dominant flow path was lost at this state.Note that the threshold value varied slightly (but not significantly) with different fracture geometric characteristics (i.e., the fracture size, shear displacement, and surface roughness characteristics).By introducing such percolation porosity, that is, φ p , Equation 2 can also be expressed as follows: The percolation porosity obtained via curve fitting of all the results was 0.17% (dashed lines in Figure 2), which is greater than granular materials (∼1%-3%, Colombier et al., 2017;Revil et al., 2014) because most of the pore spaces contributed to the flow in our fracture case relative to the porous media.Given such a small porosity in a fractured rock mass (average of 0.8% in our case), the low percolation porosity value is considered reasonable (Hunt & Sahimi, 2017).
Below this percolation threshold, the formation factor was not constant even for similar porosity data, which may reflect the different connectivity (Gueguen & Dienes, 1989), although quantitative analysis of the connectivity is difficult after disconnection of the dominant flow paths.Thus, we evaluated the tortuosity τ 2 of the flow paths because this is comparable to the connectivity c as the following relation (Text S3 in Supporting Information S1).
This relation demonstrates that the tortuosity should be equal to or less than 2.5 while the flow path is connected (c = 1).Figure 2b shows the effect of tortuosity on Archie's relation.Notably, the data sets of τ 2 ≦ 2.5 are consistent with the linear line (m = 1.2), representing the fully connected flow regime (the red line in Figure 2a).Above this threshold value, the tortuosity enhances the formation factor even when the fractures have the same porosity, thereby exhibiting nonArchie behavior (i.e., the resistivity deviates from the linear line).Higher tortuosity corresponds to a longer current path (Text S3 in Supporting Information S1), thereby enhancing the formation factor.These results clarify the role of tortuosity on resistivity.Additionally, such a percolation trend also affects the permeability-resistivity relationship, which is discussed in the following section.

Relationship Between Permeability and Resistivity
Despite of the dependencies of fracture size, shear displacement and surface roughness characteristics on permeability and resistivity (Figure 1), we found that their respective correlations exhibit stable trends invariant of these characteristics (Figure 3).The relationship between permeability k and the formation factor (hereafter referred to as the k-F relationship) can be correlated as follows with two different slopes α.  ∝  − (5) Overall, the slope α remains constant regardless of the different fracture geometric characteristics.Based on the equivalent channel model, changes in bulk permeability and formation factor at elevated stress can be interpreted theoretically as changes in tortuosity and aperture d (Paterson, 1983;Sawayama et al., 2021b; J. B. Walsh & Brace, 1984).
Focusing on the changes in permeability and formation factor from the arbitrary references (i.e.,   and   ) such as monitoring purposes, Equation 5 can be written as follows: Substituting Equation 6and Equation 7into Equation 8 derives the following (the detailed derivation is described in Sawayama et al. (2021b): where the subscript zero indicates an arbitrary reference value.This demonstrates that the slope α represents the sensitivity of the tortuosity against aperture closure and is bounded between values of 1 and 3. To this point, the inflection of the k-F relationship (at k = ∼10 −11.5 ) may represent the transition of their different sensitivities to tortuosity under fully connected and partially disconnected flow regimes, that is, the effect of percolation because of the channeling of the limited flow paths.Our results (Figure 3) for higher permeability (i.e., the fully connected flow regime) are plotted on the straight line of α = 2.7 (close to 3), which indicates that the change in tortuosity is very sensitive to aperture change.As shown in Figure 3, the fully connected flow regime is consistent with Archie's relation.Given the similarity of Equation 2and Equation 7, we can derive the following relation (Text S4 in Supporting Information S1): Substituting   = 1.2 into Equation 10 leads to   = 2.7.Meanwhile, our results show   = 1.0 below the percolation threshold (i.e., the partially disconnected flow regime).This is the lower limit expected based on the equivalent channel model (Equation 9), which indicates that the change in tortuosity is insensitive to aperture closure (Brown, 1989).In both cases, our findings could suggest that the permeability change can be estimated without knowing the absolute value of the permeability but with the observed resistivity data and the slope   , which can be derived from Archie's exponent (Equation 10) in the fully connected flow regime or is unity in the partially disconnected flow regime.Notably, additional 13,000 simulations with varying random seeds validated the rigidity of the proposed relationship at our simulation scale (Figure S8 in Supporting Information S1).

Implications
We have presented how permeability, resistivity, and porosity are affected by fracture closures and their geometrical features.In the field, by assuming a constant Archie's exponent, resistivity data are frequently interpreted as porosity (Becker et al., 1982;Chesley et al., 2021;Naif et al., 2016).We have clarified that the resistivityporosity relationship can change depending on the tortuosity or connectivity, which indicates that we must consider Archie's exponent depending on the local connection of the flow paths.This is particularly important when a macroscopic fracture is present.In contrast, the permeability-resistivity relationship is irrelevant to tortuosity and depends on only percolation.These findings suggest that, compared to porosity, the resistivity observation is less biased in the estimation of fracture permeability.
The presence and reactivation of mesoscale fractures and the associated changes in permeability are important for developing enhanced geothermal systems (EGS).In practice, resistivity observations in EGS projects detect resistivity changes (approximately a few percent) (Didana et al., 2017;Johnson et al., 2021;Peacock et al., 2013) that are associated with the hydraulic stimulation of preexisting fault systems.Under the assumption that α takes a value between 1 and 3 in Equation 5, these resistivity changes correspond to a few to several tens of percent increases in permeability.Note that the variability of α is entirely due to the percolation and is independent of fracture size, shear offset, fractal dimension, and fracture roughness.Below the percolation threshold (i.e., the partially disconnected regime), α = 1 indicates that tortuosity change is nearly constant under normal stress (i.e., an aperture opening/closure).The impermeable faults may correspond to this regime, where the significance of the stagnant flow area may accelerate the mineral precipitation (Sibson, 1987) or short-circuited flow paths (Hawkins et al., 2018;Okoroafor et al., 2022) (Figure 1, images iii and vi).After the percolation, that is, the fully connected regime, a greater α value indicates that tortuosity has decayed significantly and fallen to 2.5, thereby representing connectivity of c = 1 (Equation 4).Above this threshold, permeability hopefully is enhanced by approximately 2.7 orders in the same resistivity change.Although estimating the percolation transition is particularly important for practical use in the field, Kirkby et al. (2016) have argued that the percolation threshold is independent of matrix permeability but dependent on the matrix-to-fluid resistivity ratio (in our case, 5 × 10 5 ) and fracture spacing.Note that our model assumes simple single fractures; thus, the changes in transport properties in intersecting fracture networks (Kirkby & Heinson, 2017;Wu et al., 2021) should be investigated further.
The success of linking the k-F relationship may be further discussed with the mechanical properties as the hydromechanical relationship is theoretically formulated (Pyrak-Nolte & Nolte, 2016;Wang & Cardenas, 2016).Hence, we attempted to incorporate the mechanical property into the resistivity model to further expand the applications.Provided the negligible matrix porosity and constant thickness of the fractured rock mass, the change in porosity of a single fracture is comparable to the aperture change.The stress-induced aperture change (i.e., displacement) can be given by the effective normal stress σ divided by the fracture specific stiffness κ, and thereby Equation 2 can be rewritten as follows: where the subscript zero indicates an arbitrary reference value.Figure 4 plots the relationship between the normalized formation factor and stress-induced aperture closure based on our fault data sets considering the lowest stress condition.Overall, the relationship lies on a straight line (m = 1.2) independent of the fault size and its microstructures at σ < ∼70 MPa (i.e., percolation threshold), demonstrating the validity of Equation 11 in the fully connected flow regime.This proposed formula of the electro-mechanical relationship implies that resistivity data possibly predict either stiffness or stress state.Given that crustal stress can be considered constant (i.e., on relatively short time scales), the observed changes in resistivity with changes in the effective normal stress represent changes in pore pressure.Note that a small fluctuation in the stress (shown in color in Figure 4) is due to the stress-dependent fracture specific stiffness, which varies with the roughness and size of the fault (Worthington, 2007).Importantly, the stress prediction line remains valid across various fracture data sets and additional 13,000 simulations conducted with varying random seeds (Figure S9 in Supporting Information S1).Given the in situ stiffness through a seismic survey (Hunziker et al., 2020;Minato & Ghose, 2016), Equation 11 can be used to predict the changes in pore pressure.In field measurements, a promising method to estimate in situ pore pressure is by applying the empirical fitting in terms of crack porosity reduction to seismic profiles (Kitajima & Saffer, 2012;Tonegawa et al., 2017;Tsuji et al., 2008Tsuji et al., , 2014)).However, the commonly used effective medium theories for elasticity ignore the macroscopic fracture due to the noninteraction assumptions.The proposed formula using formation factor is preferable to estimate the changes in pore pressure associated with fault movements.Although further verification is required for practical application, we expect that this may provide new insights into estimating pore pressure in faulted regions using electrical resistivity monitoring.

Conclusions
This study investigated changes in permeability, resistivity, and their respective relationships at elevated normal stress by performing numerical simulations of various fracture models with varying fracture size, surface roughness, and shear displacement.Based on the deformed rough surfaces simulated using a half space-based dry contact model, we calculated the fracture flow using the lattice Boltzmann method and the electrical resistivity using the finite element method.The numerical results demonstrate that changes in fracture permeability and electrical resistivity at elevated normal stress are affected by fracture size, roughness, and shear displacement.Nevertheless, the relationships between permeability-resistivity exhibit less dependence on these geometric characteristics.The log-log plot of the permeability-resistivity shows inflection at the fracture permeability of approximately 10 −11.5 m 2 , and the microscopic flow analysis further supports that this inflection is related to the percolation of the dominant flow path.The trend of the relationship could be correlated with Archie's cementation exponent based on the equivalent channel model.These findings suggest that the permeability evolution of faults can be formulated with changes in resistivity regardless of the fracture size and its microstructures.In contrast, the resistivity-porosity relationship can vary depending on the tortuosity or connectivity, suggesting that the resistivity observation is less biased when directly estimating fracture permeability compared to porosity.In addition, an electro-mechanical relationship is proposed to estimate the changes in the stress state associated with fault movements.Although these proposed formulas are limited to single fractures and require further verification for practical application, they may provide new insights into electrical resistivity monitoring for estimating changes in permeability and pore pressure in fracture-dominated areas.

Figure 1 .
Figure 1.Plots of permeability and resistivity changes as a function of normal stress (a-f) and representative images derived from the simulation showing the local flow distributions (color) overlayed on the aperture distributions (grayscale) for (i-iii) joints and (iv-vi) faults.The stress levels of images (i-vi) are shown inpanel (a).L, δ, s, and D represent the fracture length, shear displacement, roughness, and fractal dimension, respectively.

Figure 2 .
Figure 2. Correlations between formation factor and porosity with respect to (a) flow path percolation and (b) tortuosity.Note that the formation factor indicates bulk resistivity normalized by fluid resistivity.The solid and dashed lines are approximation lines based on Equation 2 (m = 1.2) and Equation 3 considering a percolation porosity of 0.17%, respectively.

Figure 3 .
Figure 3. Correlation between formation factor and permeability compiled for all fracture data sets.The symbols' color and shape represent different fracture length scales, shear displacement, roughness, and fractal dimension values.The solid and dashed lines denote prediction lines based on Equation 5 for α = 2.7 and α = 1.0, respectively.

Figure 4 .
Figure 4. Correlation between fracture aperture changes and formation factor normalized to the values at σ = 5 MPa compiled from all data sets of faults.The colors and shapes of the symbols represent effective normal stress and shear displacement, respectively.The solid line denotes the prediction lines based on Equation 1.