Hydraulic Transport Through Calcite Bearing Faults With Customized Roughness: Effects of Normal and Shear Loading

Understanding fluid flow in rough fractures is of high importance to large scale geologic processes and to most anthropogenic geo‐energy activities. Here we conducted fluid transport experiments on Carrara marble fractures with a novel customized surface topography. Transmissivity measurements were conducted under mechanical loading conditions representative of deep geothermal reservoirs (normal stresses from 20 to 70 MPa and shear stresses from 0 to 30 MPa). A numerical procedure simulating normal contact and fluid flow through fractures with complex geometries was validated toward experiments. Using it, we isolated the effects of roughness parameters on fracture fluid flow. Under normal loading, we find that (i) the transmissivity decreases with normal loading and is strongly dependent on fault surface geometry and (ii) the standard deviation of heights (hRMS) and macroscopic wavelength of the surface asperities control fracture transmissivity. Transmissivity evolution is nonmonotonic, with more than 4 orders of magnitude difference for small variations of macroscopic wavelength and hRMS roughness. Reversible elastic shear loading has little effect on transmissivity; it can increase or decrease depending on contact geometry and overall stress state on the fault. Irreversible shear displacement (up to 1 mm offset) slightly decreases transmissivity and its variation with irreversible shear displacements can be predicted numerically and geometrically at low normal stress only. Finally, irreversible changes in surface roughness (plasticity and wear) due to shear displacement result in a permanent decrease of transmissivity when decreasing differential stress. Generally, reduction of a carbonate fault's effective stress increases its transmissivity while inducing small shear displacements does not.


Introduction
Fluids are pervasive in Earth's crust. They interact with rocks, modifying their physical properties and deformation mechanisms. In turn, host rocks control the way fluids migrate in the crust either due to natural forcing or to anthropogenic activities (Sibson, 1994(Sibson, , 1996. Rock masses in the brittle crust are pervasively fractured and have permeabilities ranging from around 10 −16 to 10 −17 m 2 (Faulkner et al., 2010;Townend & Zoback, 2000). These permeability values are more than 2 to 3 orders of magnitude larger than those of the intact rock matrix (10 −21 to 10 −19 m 2 ) at depths ranging from~2 to 15 km. Most of the fluid flow needs therefore to be controlled by single fractures or fracture networks. Thus, it is of outmost importance to understand how fractures and faults transport fluids in the subsurface. This is particularly valid for the safe and clean development of underground anthropogenic geo-energy activities such as geothermal energy exploitation (Breede et al., 2013;Violay et al., 2015Violay et al., , 2017. Indeed, a popular strategy for enhancing fluid transport in Enhanced Geothermal Systems (EGS) is fracture hydro-shearing, by fluid injection. It consists on reactivating preexisting faults to increase the deep crystalline reservoir's permeability (Breede et al., 2013;Cladouhos et al., 2010;Cornet, 2016). Nevertheless, the enhancement of fluid flow following stimulations in such reservoirs remains poorly predicted. A too low subsurface production flow rate results in economic losses while too high flow rates can lead to fluid leak off and reactivation of faults located far from the injection wells. The latter was probably the case of the St. Gallen geothermal project (Zbinden et al., 2019) and possibly of several other injection induced seismicity cases (Ellsworth, 2013;Goebel & Brodsky, 2018;Grigoli et al., 2017;Kim et al., 2018;Lengliné et al., 2014;Yeck et al., 2016). The poorly estimated flow rates partly arise due to the difficulties in detecting the fracture networks in the underground and partly due to the difficulties of estimating fluid flow through rough fractures with complex surface topographies, submitted to large stresses. To estimate the fluid transport capacity of a rough crack, it is useful to consider its The application of normal stress has been shown to increase the real contact area between the two fracture walls and to reduce the geometrical aperture and hydraulic transmissivity (Brown 1987b;Brown et al., 1998;Kang et al., 2016;Pyrak-Nolte & Morris, 2000;Renshaw, 1995;Rutter & Mecklenburgh, 2017, 2018Walsh, 1981;Watanabe et al., 2008Watanabe et al., , 2009Witherspoon et al., 1980). After passing a stress value (percolation threshold), the fracture reaches a configuration where further increases in normal stress result in small changes of the hydraulic aperture (thus of transmissivity). Then, fracture transmissivity remains constant due to the formation of preferential channels in between the highly stressed asperities (Brown et al., 1998;Kang et al., 2016;Pyrak-Nolte et al., 1988;Watanabe et al., 2008Watanabe et al., , 2009). The influence of reversible shear loads (in the elastic domain) has been rarely studied experimentally. In some few observations, it is seen that reversible shear loading (elastic loading, with no displacement) can cause a slight decrease in fracture transmissivity (in relatively smooth fractures of hard rock; Faoro et al., 2009;Rutter & Mecklenburgh, 2017, 2018. Most of the efforts have been put to determine the effect of irreversible shear displacement on fracture transmissivity, usually considering large displacements (more than 1-20 millimeters) at low stresses (usually lower than 20 MPa) and/or on rock fractures generated by tensile or shear fracturing as well as on artificial rock proxies (Carey et al., 2015;Chen et al., 2000;Esaki et al., 1999;Ishibashi et al., 2012;Lee & Cho, 2002;Olsson & Brown, 1993;Pyrak-Nolte et al., 1988;Watanabe et al., 2008Watanabe et al., , 2009Wenning et al., 2019;Yeo et al., 1998). From such studies, the usual knowledge with respect to the influence of shear displacement on transmissivity is that it strongly increases hydraulic transport on the fault. In contrast, recent studies (Rutter & Mecklenburgh, 2017, 2018 have shown that for displacements inferior to 1 mm, on real rock samples with smooth surfaces, at high stresses (up to 100 MPa normal stress), the transmissivity rather decreases or remains fairly constant with increasing shear displacement. These types of studies seem more relevant to fault reactivation due to anthropogenic activities (in EGS stimulation for example) particularly because the reactivation of reservoir faults needs to target small displacements to avoid large magnitude seismicity.
In this work, we developed an experimental technique to customize the roughness of hard-rock fracture surfaces by imposing different macroscopic wavelengths in suborthogonal and subparallel directions with respect to the sense of fluid flow. Then, the fluid flow through the wavy-rough fractures was experimentally measured both under normal loading only (up to 40 MPa) and under reversible shear loading (up to shear and normal stresses close to 30 and 70 MPa, respectively). The fractures loaded in shear were then submitted to irreversible shear displacement (up to 1 mm total offset). A numerical procedure that first simulates the normal contact between wavy-rough surfaces and then fluid flow through them was developed and verified with the experimental results. It is noteworthy that the numerical procedure consists on a combination of open-source models. Through the use of the calibrated numerical procedure, we isolated the effects of roughness parameters on fracture transmissivity under normal load. The numerical procedure was also used to isolate the influence of reversible elastic shear loads on the experimentally measured transmissivity. Finally, we evaluated how the small transmissivity changes during irreversible shear displacement can be predicted by a change in geometry of the fracture surface for different applied normal stresses.

Starting Samples With Customized Roughness
The samples were 75 mm long cylinders (36 mm diameter) of Carrara marble. Carrara marble has a well characterized mineralogy (~99 vol% calcite), low porosity (~1%), fine grain (average grain size <0.5 mm), and high homogeneity and isotropy (Chen, 1995;Delle Piane et al., 2015;Pieri et al., 2001). Its mechanical properties (Young's modulus~30 GPa and uniaxial compressive strength~160 GPa; Edmond & Paterson, 1972;Paterson & Wong, 2005) make it a standard in rock mechanics and an ideal material for laboratory testing. For normal loading experiments, two semicylinder's vertical flat faces were ground prior to sample coring, to obtain a perfect semicircular geometry. For shear loading and reactivation experiments, the cylinders were cored first and then saw cut at 30°toward the cylinder's long axis to create an oriented fracture. The fracture's faces were ground flat to ensure perfect contact. Finally, an injection/extraction borehole of 2 mm diameter was drilled in each half sample from the horizontal flat surface (in contact with the top/bottom anvils) to inject/extract fluid directly into/from the fracture. In saw cut configuration, the resulting fracture was of elliptical shape with long axis 2a e ¼ 72 mm and short axis 2b e ¼ 36 mm. The elliptical contacting fracture surfaces had a nominal area A~2036 mm 2 . Prior to loading, the distance between boreholes centers was of 60 mm.
For all experiments, a customized fracture roughness was imposed to each flat surface of the half samples using a vertical-axis milling machine. The machine is composed of three main elements: (i) a table where the half samples were locked and leveled to a horizontal position. (ii) A rotary milling cutter mechanically linked to a rotating spindle whose spin is controlled by a motor. The rotary cutter can be lowered to enter in slight contact with the half sample. (iii) An automatically advancing arm mechanically fixed to the rotary cutting tool. As the rotary cutting tool advances, it periodically removes rock material over the tool blade's edge, making arc shaped grooves on the sample's surface. The grooves' wavelengths are smaller for faster advancement speeds and larger for slower ones, resulting in customized roughness depending on the chosen advancement speed.
The different experimental geometries tested here are detailed in Table 1 and sketched in Figure 1b. In the sample's names, the first subscript denotes large or small macroscopic wavelength (L or S, corresponding to wavelengths of 1.7 and 0.9 mm, respectively) and the second subscript denotes the sense of fluid flow with respect to the macroscopic grooves (P for subparallel and O for suborthogonal). Finally, M f denotes the sample with no imposed macroscopic wavelength where the roughness was manually imposed through #80 grit. In all experiments, fluid flow occurred following the y axis (fracture's long axis). In experiments with shear loading and displacement, shear occurred along the y axis (vertical in the rest of this articles' figures).

Roughness Measurement and Data Processing
The measurement of surface roughness was performed using a 3-D optical profilometer (Contour GTI-3D Optical Microscope by Bruker Nano surfaces Division). The tool uses green light interferometry to determine the surface topography of the sample with an accuracy down to~100 nm. The green light pulse has an area of 1 mm 2 . A motorized base allows sample movement in the x and y directions (minor and major axis of the fracture, respectively). The tool allows automatic scanning of large surfaces by performing several measurements with a given overlap (here of 20%), which are later stitched together to reconstruct a larger surface topography.
Under this configuration, two overlapping areas of sample M LO were analyzed. The first area had a surface of 1 cm × 1 cm and the second area had a 3 cm × 3 cm area. The measurement results showed that the surfaces' statistical properties (radially averaged 2-D power spectral density [PSD] of heights) are transitionally invariant. Hence, it is assumed that taking only a portion of 1 cm 2 of the sample's surface instead of taking the whole area gives statistically the same result. Thus, for time purposes, only an area of 1 cm 2 was analyzed on the profilometer for all other samples.
The following corrections were then applied to the measured data (x, y, and z topography maps): (i) tilt removal. The intrinsic tilt due to leveling error at measurement was removed. (ii) Interpolation of missing points. Missing values are a specific consequence of rough surface optical measurements because the reflected light path can be cut when large slopes are encountered (Jacobs et al., 2017). A 2-D nearest neighbor interpolation technique (Pingel et al., 2013) was used to interpolate the missing data points.
(iii) Correction for sampling artifacts. The sampling theorem states that the minimum wave number to be considered in spectral analysis should be smaller than the Nyquist frequency the total number of points in the sampled domain (N ¼ 5,044) and L m ¼ 10 mm is the length of the measured domain. Thus, the cutoff wavevector (e.g., the maximum wavevector analyzed) should be q cut ¼ 2.5e5 m −1 . A low-pass Gaussian filter was applied to remove all wavevectors higher than q cut in the data. Finally, to evaluate the properties of rough surfaces ( Figure 6), a radially averaged 2-D PSD analysis with radially symmetric Welch windows was performed to avoid artifacts (Jacobs et al., 2017;Kanafi, 2019).
It is important to notice that in the rest of the manuscript the roughness parameters are evaluated with the available data. For example, the Hurst exponent is evaluated on windows smaller than 1 order of magnitude in wavevector, thus leading to an intrinsic error related to the availability of the data (more data could be obtained through higher/lower resolution measurements to complement the data set). Similar difficulties arise on the estimation of the roll-off wavevector for the experimental samples. Notwithstanding the estimated intrinsic errors, the same technique was applied for all the measured samples. Thus, the comparative analysis presented remains robust even though the absolute values of these parameters might not be as accurate as desired.

Experimental Setups and Flow Through Experiments
The experimental set up was an oil-medium triaxial Hoek-cell (Figure 1a) of the Laboratory of Experimental Rock Mechanics (LEMR) at EPFL, Switzerland. The cell can hold 70 MPa (±50 kPa resolution) in confinement pressure (σ 3 ¼ σ 2 ). For flow through experiments, the top and bottom anvils were specifically designed to allow controlled fluid pressures and volumes independently at the top and bottom ends of the samples ( Figure 1a). The pressure/volume controllers have a capacity of 200 cm 3 (±1 mm 3 resolution) in volume and 30 MPa (±10 kPa resolution) in pressure.
One experiment was performed to evaluate the matrix permeability of Carrara marble and have a point of comparison for the fracture fluid flow experiments. Due to the low permeability of the rock matrix, an oscillatory fluid flow method was used under the same experimental setup. Details of the oscillatory fluid flow method can be found in Bernabé et al. (2006) and Acosta and Violay (2020). Matrix permeabilities ranged from~5.99 e-19 m 2 at σ 3 ′ ¼ 8 MPa to 4.92e-20 m 2 at σ 3 ′ ¼ 20 MPa (Appendix Figure A1) and the exponential decay seems consistent with previous literature studies of the permeability of Carrara marble (Chen, 1995;Delle Piane et al., 2015;Zhang et al., 1994) even if, in this study, measurements at low effective stress were not performed for time purposes.

Experiments Under Normal Loading
The half cylinders were clamped together and let to saturate in a vacuum chamber for a minimum of 1 week. Following this, samples were confined to σ 3 ¼ 5 MPa and fluid pressure (p f ) to 1 MPa during a minimum of 120 min for additional (pressurized) saturation. Once fluid and confinement pressures and volumes reached an equilibrium, σ 3 was increased to the target pressure of 43 MPa and p f was changed stepwise to study the effect of effective normal stress (σ ′ 28,30,32,34,36,38,40 MPa; Figure 2a). At each step, a differential pressure Δp f ¼ 0.3 MPa was imposed between the top and bottom ends of the sample and the steady state fluid flow rate (Q in m 3 . s −1 ) was measured (Figures 2a and 2b, insert in panel b shows one example of the steady state flow rate). Because the flow rate on the fracture was more than 3 orders of magnitude larger than in an intact marble cylinder (Appendix Figure A1), it is reasonable to assume that all the flow occurred through the fracture.
For each σ ′ N step, the fluid flow through the rectangular fractures was quantified by the product of the permeability (k in m 2 ) and the effective thickness (t in m) (Rutter & Mecklenburgh, 2017, 2018. kt is called the fracture's hydraulic transmissivity (kt in m 3 ) which can be estimated directly from Darcy's law as follows: with μ the dynamic viscosity of the fluid, w the fracture's width, and L its length.

Experiments Under Shear Loading
The procedure to saturate the samples, place them in the cell, and take them to isostatic loading (σ′ 1 ¼ σ′ 3 ) was the same as for normal loading experiments. For the shear experiments, confining pressure σ 3 was either 15 or 35 MPa and p f was 5 or 15 MPa, respectively (so that the effective confinement σ′ 3 ¼ 10 and 20 MPa). Then, the axial displacement was increased by steps of 0.1 mm at a displacement rate of 10 −6 mm.s −1 . Such a low displacement rate was used to allow fluid pressures equilibrium on the fault during shear loading (i.e., fault drainage). Under saw cut configuration, both shear and normal stresses on the fault increased with increase of axial displacement and were calculated as follows: with θ the angle between the saw cut and the vertical.
The final axial displacement in our experiments was of~1.1 or 1.2 mm. As fault reactivation (e.g., departure from elasticity and initiation of sliding) occurred often slightly after~0.1 mm displacement, the final shear offset was of~1 mm in most experiments. At every displacement step, the piston's position was held constant. The differential stress naturally relaxed as the piston stopped. When the differential stress reached a steady state, a differential fluid pressure of 0.3 MPa was imposed between the injection and extraction boreholes to measure sample's transmissivity (Figures 2c and 2d).
The steady state flow rate Q (Inserts in Figures 2b and 2d) was determined and the hydraulic transmissivity was estimated from the flow lines in a perfect elliptical surface using the dipole image method of Mecklenburgh (2017, 2018) such that the transmissivity (product of permeability and equivalent hydraulic aperture) writes (Rutter & Mecklenburgh, 2017, 2018: with a e the half distance between the injector and extractor boreholes, r 0 the borehole diameter, dP dx the spatial pressure gradient between the boreholes, and B a constant close to unity. It is noteworthy that here, no corrections for the changes in elliptical surface geometry were made (see Tembe et al., 2010) because the total displacement on the saw/cut was <1.2 mm, resulting in a change in nominal contact area lower than 8% which would result in less than 0.5% change in transmissivity.

Experimental Results
In all figures, a schematic legend is presented to show the customized sample's roughness. We note that the samples prepared for flow-through experiments under normal loading are not the same as those for experiments under shear loading; thus, slight differences between the sample's roughness can be found for a same nomenclature. Figure 3 shows the results of flow-through experiments obtained in terms of transmissivity as function of effective confining pressure. The dashed lines represent a first cycle where σ 3 ′ was decreased from 40 to 28 MPa (by increasing fluid pressure). The full lines represent a second cycle where the effective confinement was increased from 28 to 40 MPa. Note that the measured transmissivities were always higher when σ 3 ′ was decreased (first loading) than those of the increasing cycle (second loading). To avoid issues related to this hysteretic behavior of the fractures, the transmissivities that will be used hereafter are those of the increasing effective confinement cycle (second loading) because they are supposed to be representative of the fracture's transmissivity under elastic behavior (Iwai., 1976;Rutter & Mecklenburgh, 2017, 2018Witherspoon et al., 1980). In these experiments, we did not perform more (3 or 4) effective pressure increase and decrease phases to assess the reversibility of the transmissivity. Instead, we assume that the second loading is representative of an elastic behavior of the fractures with loading-unloading phases, as has been shown in previous studies (Rutter & Mecklenburgh, 2017, 2018.

Fluid Flow Through Single Fractures: Normal Loading
The experimental sample M LO (e.g., large wavelength suborthogonal to fluid flow) showed transmissivities ranging from 3.05e-18 m 3 at σ 3 ′ ¼ 28 MPa down to 0.76 e-18 m 3 at σ 3 ′ ¼ 40 MPa. Then, M SO (e.g., small wavelength suborthogonal to fluid flow) showed transmissivities half an order of magnitude smaller than those of M LO (ranging from 0.49e-18 m 3 down to 0.17 e-18 m 3 at σ 3 ′ ¼ 40 MPa). The sample with no imposed macroscopic wavelength M f had transmissivities ranging from 0.23e-18 m 3 down to 0.15 e-18 m 3 . The transmissivities were close to those of M LO and M SO , but the decay with increasing confinement was smaller in this experiment with respect to the samples with macroscopic wavelength. The sample M LP (e.g., large wavelength subparallel to fluid flow) had transmissivities ranging from 0.48e-18 m 3 at σ 3 ′ ¼ 28 MPa down to 0.20 e-18 m 3 at σ 3 ′ ¼ 40 MPa. The transmissivities were almost half an order of magnitude lower than those of M LO and close to those of M SO . Figure 4 (left axis) shows the shear stress versus axial displacement curves of all the conducted experiments. In all cases, shear stress first increased elastically (i.e., in a reversible manner) in response to increases in axial piston displacement (notice that due to the fault orientation and loading configuration, the normal stress also increased [in a reversible manner] on the fault during elastic loading). During this stage, the faults were fully locked (Acosta et al., 2019;Byerlee & Summers, 1975;Ohnaka, 2013;Scholz, 2019). Then, once the shear strength of the faults was reached, reactivation (i.e., initiation of motion) occurred and shear stress versus displacement curves showed a roll over (at displacements d ro~0 .11 − 0.20 mm) until reaching a steady state where shear stress stayed close to constant with increasing axial displacement. During this stage, the faults were unlocked and slipped at a near-to-constant rate. It is noticeable that because axial displacement was increased stepwise to measure transmissivity, the fault showed an increase of shear strength at the start of every new displacement step. This restrengthening behavior usually represented less than 10% stress change with respect to fault's shear strength. The peak values of stress and its relaxation are due to the time dependence of the fault's real contact area (Dieterich, 1979;Dieterich & Kilgore, 1994). It is also noteworthy that all samples with an imposed macroscopic wavelength showed a near constant increase in shear stress with displacement after reactivation occurred (Figures 4a, 4b, and 4d) at friction values (τ=σ′ N ) close to 0.5 (Table 1). The sample without macroscopic wavelength first showed a (very slow) stress drop at reactivation and then a near-constant increase in shear stress with increasing axial displacement ( Figure 4c).  (Byerlee, 1978). It is noteworthy that the "steady state" friction did not show any significant trend with respect to the sample's roughness configuration. A compilation of the values is given in Table 1.

Transmissivity Results
The transmissivities measured during shear loading experiments are shown in Figure 4 (right axes). In experiments conducted at σ′ 3 ¼ 10 MPa (lighter colors), the initial transmissivities (kt 0 ; e.g., with no applied deviatoric stress and at zero axial displacement) ranged from~2.33.10 −16 to 2.42.10 −18 m 2 with maxima and .) During irreversible shear displacement, transmissivity usually slightly decreased overall with local rises (to kt max ss ) and drops (to k t min ss ) of lesser magnitude than the decrease during elastic loading. Finally, after unloading, transmissivity (kt unl ) slightly increased in most cases (except for sample M LP deformed at σ′ 3 ¼ 20 MPa) but was far from being recovered to kt 0 . Transmissivity after unloading (e.g., under no applied differential stress) was usually close to the value found at the onset of reactivation. For experiments conducted at σ′ 3 ¼ 20 MPa, kt 0 was 3 to 42 times larger than at σ′ 3 ¼ 10 MPa (except for M LP where kt 0 was surprisingly 2 orders of magnitude higher at σ′ 3 ¼ 20 MPa). At the onset of reactivation, kt ss σ′ 3 ¼ 20 MPa ð Þwas 13 to 392 times lower than kt ss at σ′ 3 ¼ 10 MPa). Finally, at unloading, kt unl was 11 to 105 times lower at larger confining pressure.

Microstructures
The surface topography maps of intact samples are presented in   Figure 6a shows the PSD curves for the initial experimental surfaces (shown in Figures 5a-5d) as surface topography maps). The 2-D PSD curves presented two sections: a first "flatter" part where the PSD was close to constant with increasing wave number until the roll off wave number q r . In this first part, the macroscopic wavelengths constitute a peak in the PSD curves. With increasing wave number, a second part presenting a power law dependence on wave number was observed. The slope in a log-log plot is −2(H + 1) with H the Hurst exponent (Candela et al., 2012;Jacobs et al., 2017). The Hurst exponent characterizes the power law decay of PSD with increasing wavelength. In that sense, H usually characterizes the fractal dimension of a surface (Candela et al., 2012;Jacobs et al., 2017, and references therein). Finally, the area under the PSD curves represents the root-mean-square height (h RMS ), which is the standard deviation of the heights distribution (Candela et al., 2012;Jacobs et al., 2017, and references therein).

Initial Surface Roughness
Prior deformation, for q < q r , the samples M LO had PSD amplitudes (C(q < q r )~2.10 −19 m 4 ), and the samples M SO and M LP had smaller PSD amplitude prior to roll off (C(q < q r )~7-8.10 −20 m 4 ). Finally, the sample with no macroscopic wavelength M f had the smallest PSD amplitudes prior to roll-off (C(q < q r )~1.10 −20 m 4 ). The roll of wave numbers were the smallest for M LO , M LP , and M SO (q r~6 ,900 rad.m −1 ). q r were larger for the sample with no macroscopic wavelength M f (q r~1 0,000 rad.m −1 ) ( Table 1) Table 2.

Postdeformation Surface Roughness
Fault surface roughness visibly changed in sheared samples. The topography maps after deformation are shown in Figure 5 (at low confining pressure in panels e-h and at high confining pressure in panels i-l). Note that a different sample was used for each confining pressure experiment to initiate loading under similar conditions. Overall, the postmortem samples showed evidence of striation (grooves in the sense of shear), as well as changes in the height distributions ( Figure 5, one column can be observed for comparison). Evidence was found of gouge formation during shearing with pervasive presence of microscopic particles in low height zones. Larger amounts of gouge were generated in experiments at higher confining pressure. It is noticeable that at both effective confinements, the samples M LP (e.g., with grooves subparallel to the shear sense), showed very large changes in the surface characteristics (Figures 5d, 5h, and 5l). There, the initial macroscopic wavelengths were unrecognizable from topography measurements while in all the other experiments, the initial surface topography could be partly recognized.

Journal of Geophysical Research: Solid Earth
To study the statistical properties of these surfaces, the PSD's were computed again on postmortem topography maps and analyzed in the same manner as those of intact surfaces (Figures 6b and 6c). The PSD curves showed a change in shape. Indeed, the previously "flat" part of the PSD's showed an overall slope after deformation, adding additional complexity to the estimation of the roll-off wavevector for postmortem samples. The results are compiled in Table 3 and summarized in Figures 6d-6f. The root-mean-square of heights h RMS decreased after shearing in all cases (except for the sample M LP deformed at σ ′ 3 30 MPa) (Figure 6d). No tendency was observed regarding h RMS with respect to the confining pressure at which the samples were deformed. q r was here estimated where the slope of the PSD curves changed. With that (rough) estimation of q r , an overall decrease was observed for all samples after shearing with the measurements converging on all samples toward values of~5,000-8,000 rad.m −1 (Figure 6e). Finally, the Hurst exponent also showed an overall decrease after shearing for all samples with slightly lower values of H in experiments at higher effective confinement (Figure 6f).

Generation of Artificial Surfaces
Artificial wavy-rough surfaces were independently generated through use of the algorithm by Kanafi (2018). The algorithm uses the roughness parameters (measured on the experimental samples with the profilometer) from the PSD of surface heights (h RMS , H, and q r ) to generate an artificial randomly rough surface with the corresponding properties. In addition to the randomly generated rough surface, the experimental samples had a customized macroscopic wavelength (see section 4), which represents a singularity at a given wavevector in the radially averaged PSD's (Jacobs et al., 2017). The macroscopic wavelength and the corresponding amplitude were evaluated through the profilometer measurements (Table 2 and Figure 5). The final surfaces are the resultant of a random roughness created from the artificial surface generator on top of a sinusoidal macroscopic wavelength estimated from the experimental samples (Figure 7a). The total roughness was adjusted so that the h RMS of the sum of the two surfaces is equal to the true h RMS measured on the experimental sample.

Surface Contact Under Normal Stress
To simulate contact of two opposing surfaces resulting from lithostatic pressure in Geo-energy reservoirs (represented by σ 3 ′ in the experiments), a half space-based, dry contact model from Tribonet was used (Akchurin et al., 2015;Lubrecht & Ioannides, 1991;Tribonet n.d.). The model uses the artificially generated surfaces discretized to either 2,048 × 682 (for model calibration) or 768 × 256 nodes (for parametric analysis). a Parallel or orthogonal to fluid flow. b Determined by counting the number of wavelengths in a given area. c Determined from the difference between the six prepared half samples. d Determined from the roughness profiles. e Determined from the integral of radially averaged 2-D PSD and from the standard deviation of heights. f Determined from the differences between integral of radially averaged PSD and from the standard deviation of heights and the six half samples. g Determined from the flat part of radially averaged PSD. h Determined from the difference between the six prepared half samples. i Determined from the self-similar part of radially averaged PSD. j Determined from the difference between the six prepared half samples and from the maximum and minimum slope taken on a moving window of half an order of magnitude.

Journal of Geophysical Research: Solid Earth
Solid material properties were assigned to the contact bodies (which can differ but are here taken equal) described by a saturating elastic stress-strain relationship (e.g., the deformation is purely elastic until a stress threshold is reached, then stress remains constant with increasing strain). The parameters used here are the material's Young's modulus E (here 30.2 GPa) and the Poisson's ratio ν (here 0.3) for elasticity (measured from a Marble deformation experiment shown in Annex 1), and the yield stress σ y (here 0.2 GPa; Violay et al., 2014), which describes the limit of plasticity (or saturation threshold). The simulated load applied between the half spaces corresponds to the macroscopic normal stress (here σ 3 ′) over the nominal contact area A n ¼ L · w (see Table 1). The contact problem is solved under plane-strain boundary conditions and takes into account the mechanical interactions between microcontacts by the use of a double-continuum convolution integral (Polonsky & Keer, 1999). It calculates how the deflection at each mesh node affects the surrounding nodes. The calculation iterates until convergence of the deflection at all nodes (Akchurin et al., 2015;Lubrecht & Ioannides, 1991;Polonsky & Keer, 1999;and Tribonet n.d.). This calculation allows for fairly realistic simulation of real contact area evolution under normal load. In that sense, the roughness of each half surface is allowed to change under the application of normal load. As outputs, two-dimensional real contact area (A r ) maps and geometrical aperture (e m (x,y)) maps were recovered at each studied effective confining pressure (Figure 7b). Note that, in the aperture maps, a zero-aperture value is not allowed in our procedure. Thus, the contacting zones were replaced with apertures more than 10 orders of magnitude smaller than the mean aperture to avoid numerical issues for fluid flow calculation.

Fluid Flow Calculation
Finally, once that the contact area and geometrical aperture maps were extracted under different normal loads, the flow through the rough fractures was resolved (Figure 7c). A finite volume formulation (Brush & Thomson, 2003;Crandall et al., 2017) was used to solve the Reynolds lubrication equation (Reynolds, 1886). To apply the Reynolds lubrication equation (here simplified to the local cubic law; Zimmerman & Bodvarsson, 1996) in the fracture, the main assumption is that the variations in aperture occur gradually in space over the fracture plane (smoothness hypothesis). This hypothesis seems reasonable because (i) h RMS λ ≪1, thus the vertical variations of roughness with respect to macroscopic wavelength are small, and (ii) L λ > 10, thus the aperture due to macroscopic wavelength is small compared to the total fracture length. The Reynolds boundary layer approximation can therefore be expressed as follows (Brown, 1987b;Jaeger et al., 2009;Watanabe et al., 2008Watanabe et al., , 2009Zimmerman & Bodvarsson, 1996): where ρ and μ are the fluid density and dynamic viscosity, respectively, e m (x, y) is the local mechanical (or geometrical) aperture in the vertical direction (Brush & Thomson, 2003), S is the domain's surface, and b n is the outward unit normal vector to the local element. Details on the discretization and resolution of the mass conservation equation above can be found in Brush and Thomson (2003) and Crandall et al. (2017). While few studies have managed to simulate normal stresses in rough faults followed by the fluid flow through them (Kang et al., 2016), to our knowledge, this is the first study to use a combination of open-source numerical models for (i) generating wavy, rough surfaces, (ii) simulating the contact under effect of normal stress, and (iii) the study of fluid flow through the fractures. This procedure was used due to the relatively low computational needs if compared to a full monolithic coupled computation of deformation and fluid flow (see Shvarts andYastrebov., 2018a, 2018b;and Shvarts, 2019, for details on such approach).

Validation of the Numerical Simulations Toward Experimental Data
The hydraulic transport properties of rough fractures submitted to normal stress are highly dependent on the surface geometry and roughness parameters (section 3; and Chen et al., 2000Chen et al., , 2009Iwai, 1976;Patir & Cheng, 1978;Pyrak-Nolte et al., 1988;Rutter & Mecklenburgh, 2017, 2018Watanabe et al., 2008, Walsh, 1981Witherspoon et al., 1980). We now study the flow through wavy rough fractures submitted to normal loading only in this subsection. As described in section 5.1, first artificial surfaces with roughness parameters similar to those measured experimentally were generated. Then, the contact between the surfaces at given loads was simulated, and finally, the flow through the fractures was computed. The numerical results are presented in Figure 8a corresponding to the samples tested experimentally. For all samples, two types of numerical simulations were conducted. One where the large wavelengths were in peak-to-peak contact (e.g., nonmated surfaces, Figure 8 triangles). In those cases, the resolved numerical transmissivities were more than 2 orders of magnitude larger than those measured experimentally (Figure 3). Another set of numerical simulations was conducted where the large wavelengths were in peak-to-valley contact (e.g., fully mated surfaces, Figure 8a circles). In this case, the resolved numerical transmissivities ranged from 0.50e-18 m 2 at σ 3 ′ ¼ 28 MPa down to 0.20e-18 m 2 at σ 3 ′ ¼ 40 MPa, for sample M LP as an example. The fully mated transmissivity results are in strong compatibility with experimental results. This highlights the strong influence of flow channeling on fracture hydraulic transport capacity which cannot be neglected in the analysis (Shvarts, 2019;Shvarts & Yastrebov, 2018a;Watanabe et al., 2009). Note that in the experiments, the surface mating was difficult to control when setting up the half fractures in the triaxial cell. Nevertheless, the strong compatibility with mated numerical simulations highlights that the experimental surfaces were close to fully mated during experiments.
We observe that for all samples, the numerical transmissivity has a near to exponential decay with increasing normal stress (the plot is log normal) as kt σ ′ N À Á ¼ a:e −bσ ′ N , with a and b two fitting parameters. Note that here, we tested normal stresses that are representative of geo-energy reservoirs, therefore the lower range of stress (<28 MPa) was not explored. The sample's initial surface geometry seems to condition the para- To summarize, samples with grooves perpendicular to fluid flow are more sensitive to normal stress than samples with no roughness which are in turn more sensitive than samples with grooves subparallel to fluid flow. The transmissivity decay with increasing normal stress has been seen with (i) near exponential decays (Iwai, 1976;Pyrak-Nolte et al., 1988;Witherspoon et al., 1980), (ii) logarithmic decay (Walsh, 1981), or (iii) through more complex decays depending on heterogeneous topographies (Watanabe et al., 2008(Watanabe et al., , 2009) and loading paths (Iwai, 1976;Rutter & Mecklenburgh, 2017, 2018Witherspoon et al., 1980).
From the numerical procedure developed in this work, one can expect a near exponential decay of transmissivity at the working experimental normal stresses for the modeled carbonate rock. In turn, the experimental results show a small scatter toward the exponential decay predicted by the model, which can be due to (i) the loading path (e.g., hysteresis; Iwai, 1976;Rutter & Mecklenburgh, 2017, 2018Witherspoon et al., 1980), (ii) imperfections in the experimental contacts with respect to the numerical model (e.g., not perfectly mated experimental surfaces), or (iii) small nonlinearities in fluid flow in real sample surfaces (Zimmerman & Bodvarsson, 1996), which are not considered in the Reynolds lubrication approximation for the In both panels, triangles represent the numerical results from samples with an unmated configuration and circles represent a mated configuration. In panel (b), the black lines represent a slope of 1 and a deviation with a factor 2.5 to that line. Note that more than 90% of the mated data points were contained within that factor. The numerical model developed gives good agreement with the experimental hydraulic transport properties measured under normal loading.
simulations. Overall, the numerical results are remarkably consistent with the experimental data, as shown by the small deviation of data from the 1:1 slope in Figure 8b. In the numerical simulations, more than 90% of the points are contained within a factor 2.5 from the experimental data.

Influence of Roughness Parameters on Fracture Transmissivity Under Normal Loading
The goal of the following section is to isolate the roughness parameters that have stronger control on hydraulic transport properties through a parametric analysis of the numerical procedure. First, the parameters are isolated in absence of macroscopic grooves and then the effect of the large-scale wavelength with grooves perpendicular to the fluid flow is studied.

Flat Rough Surfaces (Absence of Macroscopic Grooves)
First, the roughness properties are analyzed for rough surfaces without imposed wavelength to avoid bias from the macroscopic grooves on all other parameters (H, q r , and h RMS ). Artificial surfaces with L ¼ 75 mm (length) and W ¼ 25 mm (width) with 768 × 256 nodes were created and tested through the same numerical procedure described in section 5.1. The steady parameters (the parameters that do not change when only one of the other parameters is varied) are H ¼ 0.6, q r ¼ 8,000 rad.m −1 , and h RMS ¼ 8 μm. Hurst exponents are varied from 0.5 to 0.9 (values that can be found in natural and experimental faults; Insert shows a sketch of flat rough surfaces and wavy rough surfaces and the flow direction. Note that in panels (b) and (e); as well as in (c) and (f ), the parameters chosen for the analysis can differ. This has been done to explore a narrower parameter range in presence of macroscopic grooves.

10.1029/2020JB019767
Journal of Geophysical Research: Solid Earth Brown & Scholz, 1985;Candela et al., 2009), roll of wave vectors are varied from 1,000 to 10,000 rad.m −1 and h RMS are varied from 2 to 8 μm. Results are shown in Figure 9.
1. Changes from 0.5 to 0.9 in Hurst exponent alone (Figure 9a) do not show large changes in transmissivity (30% to 60% lower for the lowest tested Hurst exponent) of the simulated fractures. 2. From Figure 9b, q r has a strong influence on transmissivity only for q r < 5,000 rad.m −1 . Under our experimental conditions, q r is often comprised between 5,000 and 10,000 rad.m −1 ( Table 3). In that range, q r variation has little influence over the transmissivity response. 3. From Figure 9c, we observe that h RMS has the largest influence over the transmissivity response of the flat-rough fractures.
Increase of 2 to 8 μm results in over 3 orders of magnitude difference transmissivity. The lower h RMS height, the flatter the surfaces, hence, under normal load, the better they match to each other and the smaller their mechanical and hydraulic apertures.

Wavy Rough Surfaces (in Presence of Macroscopic Grooves)
To study the effect of the macroscopic wavelength on fracture's transmissivity, artificial surfaces were generated by overlaying a surface with a macroscopic wavelength and a flat rough surface such that the overall h RMS equals the targeted value. The steady parameters are H ¼ 0.6, q r ¼ 7,000 rad.m −1 , and h RMS ¼ 4 μm and wavelength λ ¼ 1.7 mm. These reference values are taken from experimental observations. Hurst exponents are varied from 0.5 to 0.9, roll of wave vectors are varied from 4,000 to 10,000 rad.m −1 , h RMS are varied from 3 to 5 μm and wavelengths λ from 0.9 to 2.5 μm. Results are shown in Figures 9d-9g. For q r and h RMS , the parameter range explored here is narrower than for flat rough surfaces presented above. Figures 9d and 9e show again that H and q r alone have little influence on the transmissivity of wavy rough surfaces (changes of 30% to 60% in transmissivity for a change of 0.4 in H and 3 to 5 times increase in transmissivity for a variation of 6,000 in q r ). On the other hand, h RMS has strong control over the fracture's transmissivity ( Figure 9f). Indeed, at σ 3 ′ ¼ 35 MPa, the transmissivity is more than 3 orders of magnitude larger from h RMS~3 to 5 μm. We observe that for a change in λ from 0.9 to 2.5 cm, at σ 3 ′ ¼ 35 MPa, the increase in transmissivity is more than 3 orders of magnitude (Figure 9g) highlighting the strong control of the wavelength on the hydraulic transport capacity of the rock fracture.
While the physical reasons explaining the small dependence of transmissivity with H and q r remain unknown, it is noteworthy that the realistic model proposed in this study predicts a strong dependence of transmissivity with h rms and λ. To explore the combined effects of wavelength and standard deviation of heights on fracture's transmissivity, we performed 36 additional simulations (with σ′ N ¼ 28 MPa and for grooves suborthogonal to fluid flow as examples). The results are shown in Figure 10. The evolution of transmissivity is strongly nonmonotonic in this parameter space. A combination of small λ and low h RMS naturally results in low transmissivities. Nevertheless, low λ (500-1,000 μm) and intermediate h RMS (10 μm) can result in regions of lower transmissivities than its surroundings. At intermediate λ (~1,500 μm), an increment of h RMS results in less increase in transmissivity than at higher λ (~2,100 μm) for example. We conclude from this analysis that, for complex topographies, the transmissivity of wavy rough fractures from averaged height values (h RMS only for example) needs to be evaluated with care (Brown, 1987b;Hakami, 1989;Patir & Cheng, 1978;Piggott & Elsworth, 1992;Renshaw, 1995;Zimmerman & Bodvarsson, 1996).

Effect of Shear on Transmissivity: Extrapolation of the Numerical Procedure
As observed in the previous section, the numerical procedure developed here predicts well the transmissivity of rough fractures submitted to normal stress. It allows therefore to explore the influence of each surface roughness parameter on the hydraulic transport properties of the fractures. It is noticeable that the contact model does not include a shear stress component analysis nor a wear analysis. From our knowledge, integrating an elasto-plastic shear component and wear components to the analysis is not a straightforward task and is out of the scope of this study (Aghababaei et Molinari et al., 2018). In section 6.1, we use the developed numerical procedure to simplify the shear loading problem and study the influence of reversible shear loading and irreversible shear displacement on fracture transmissivity.

Influence of Reversible Elastic Shear Loading on Transmissivity (Prior the Onset of Sliding)
Before the shear strength was reached, the faults were fully locked (Acosta et al., 2019;Byerlee & Summers, 1975;Ohnaka, 2013;Scholz et al., 1972), during that phase, the faults were loaded elastically (in a reversible manner, gray shades in Annex Figure 4) in response to increases in axial piston displacement. Elastic loading strongly decreased transmissivity under all tested conditions. The change in transmissivity is due to the increase in both normal and shear stress on the fracture prior to sliding. (Note that this type of loading is more common in natural faults than a shear loading at constant normal stress.) In order to isolate the effect of reversible shear load, we now compare the results from the validated numerical model (for fractures submitted to normal stress only) to the experimental results obtained under both shear and normal stress. To do this, the model results of fracture transmissivity function of normal stress were fitted with an exponential decay as kt σ′ N ð Þ ¼ a:e −bσ′N (full colored lines in Figure 11). The values of a and b are given

Journal of Geophysical Research: Solid Earth
in Figure 11 (top of the panels). Then, we compare the model values to experimental results (measured or interpolated) of fracture transmissivity submitted to both normal and shear stress during reversible (elastic) loading (filled circles in Figure 11; darker colors for higher effective confinement). Note that assuming that the model predicts well the experimental transmissivity under normal stress, all deviations from the (normal loading) model result from the effect of shear loads in this configuration. We observe that reversible elastic shear load has very little effect on transmissivity. It either increases or decreases fault transmissivity with respect to the case where only normal load was applied. When the sample had no macroscopic roughness (M f ; Figure 11c), the transmissivity mostly decreased, in particular when high (normal and shear) stress was applied (experiment at σ′ 3 ¼ 20 MPa). Similar observations were observed in smooth hard-rock surfaces at high stresses (Rutter & Mecklenburgh, 2017, 2018. In the case of sample M LO (Figure 11a), transmissivity rose of almost 2 orders of magnitude under all stress conditions. For sample M SO (Figure 11b), transmissivity was not affected by shear stress at low confinement (first two circles from left to right, in lighter color) but increased at high confining pressure (third and fourth points from left to right). Finally, for sample M LP (Figure 11d), the transmissivity usually decreased except for one point that seems to be an outlier in the transmissivity evolution (Figure 4d; displacement < 0.1 mm).
An increase in transmissivity due to the application of reversible shear load alone should be due to a decrease of the real contact area with increasing shear stress because the hydraulic aperture (h H ) is a function of the real contact area (A r ) as (Walsh, 1981): where h M is the geometrical aperture of the fracture surface. The evolution of the real contact area with increasing shear in a frictional contact is today a controversial issue in the contact mechanics community. On the one hand, several studies have shown a decrease of the real contact area with increasing shear stress prior to (and at the onset of) sliding in mortar rock replicas (Grasselli & Egger, 2003;Park & Song, 2013), on hard polymers (Bayart et al., 2016;Ben-David et al., 2010;Svetlizky & Fineberg, 2014), and on soft polymers (Sahli et al., 2018). On the other hand, experiments in hard, coated polymers (Bay & Wanheim, 1976) and in polystyrene-on-glass contacts (Weber et al., 2019) have shown that the real contact area in turn increases with shear stress and initial displacement due to a mechanical degradation of the asperities (plastic deformation). By isolating the effect of normal stress and shear stress separately on transmissivity, our results show that the application of reversible shear load is not the only factor affecting fracture transmissivity but rather the combination of geometry and stress even if its effect is small. Indeed, both samples with macroscopic grooves suborthogonal to shear sense showed the largest increase of transmissivity. In experiments similar as those conducted here but in smooth, hard-rock samples (Rutter & Mecklenburgh, 2017, 2018, the increase of shear stress at constant normal stress led to a very slight decrease in fault transmissivity at the onset of reactivation. There, the authors inferred that the slight transmissivity decrease was due to asperity collapse during shear loading. In our experiments, it seems therefore reasonable that during elastic loading only (e.g., while no permanent changes in surface topography are expected), the transmissivity varies very little with increasing, reversible shear load. It is possible that part of this small variation comes from a hysteretic behavior of the fracture under reversible loading (e.g., Figure 3).
From this analysis, we conclude that the stress state alone does not fully determine the hydraulic transport capacity of a rock fracture. It seems that there is an interplay between fracture geometry and the state of stress (shear and normal stress applied on the fracture) that determines fracture transmissivity.
To examine the exact contribution of shear stress applied on the fracture, more complex elastic-plastic models are required. First attempts going in this direction are the models of Yastrebov et al. (2017) and those of Shvarts andYastrebov (2018a, 2018b) where the problem is approached through spectral boundary element methods and full monolithic coupling of the mechanical deformation and hydraulic transport in the fractures. These methods should allow the extension to shear stress and plasticity (Frérot et al., 2019) with the geometries used in this study (grooves suborthogonal to fluid flow and shear sense as well as parallel ones).

10.1029/2020JB019767
Journal of Geophysical Research: Solid Earth

The Effect of Irreversible Shear Displacement on Transmissivity
Once that irreversible shear displacement started and the fault reached a "steady state" sliding condition (see Annex Figure 4), with near constant shear and normal stresses, the fault's transmissivity varied only slightly, usually of less than 1 order of magnitude (Figure 4). Our experimental results of transmissivity change during shear sliding are plotted as full lines (with empty circles as markers, left y axes) in Figure 12. There, transmissivity is normalized by the mean fracture transmissivity during shear sliding.
In order to simulate the effect of shear displacement on fracture transmissivity, we simplified the problem to the shift between the two (initially mated) opposing artificial surfaces of a given displacement (Figure 7d). For this model, two artificial rectangular surfaces (32 mm × 64 mm; equivalent to the ellipsis area to simplify the problem) were generated with shifts of 0.1 mm in the y direction and put into contact using the procedure described above. For each sample, a total of 10 surface pairs was computed to evaluate the evolution of transmissivity with increasing displacement up to 1 mm total displacement. In this case, the applied load on the shifted fracture surfaces corresponds to the mean experimentally measured normal stress during shear sliding (σ′ Nss in Table 1). The numerical results from the model are plotted as thick transparent areas (with filled triangles as markers, right y axes) in Figure 12. Again, the values are normalized by the mean fracture transmissivity during shear sliding. On the one hand, in experiments conducted at low confining pressure (σ′ 3 ¼ 10 MPa; Figures 12a-12d), the general tendency of transmissivity change is fairly well captured by the simplified model even though the absolute values can be off by more than an order of magnitude. On the other hand, at high confining pressures ( σ ′ 3 ¼ 20 MPa; Figures 12e-12h), the simplified model is far from capturing the tendency of transmissivity change with ongoing shear slip. We interpret these differences at low and high confining pressure as due to strong changes in surface roughness (due to plastic deformation and wear) with ongoing slip. Indeed, the model presented here is extremely simplified because it does neither consider the evolution of an applied shear stress on the fracture's macroscopic wavelengths nor the production of wear particles (Aghababaei et al., 2016;Milanese et al., 2019;Molinari et al., 2018). The plasticity (Frérot et al., 2019) and wear (Archard, 1959;Molinari et al., 2018;Tesei et al., 2017;Yamashita et al., 2015) processes should be enhanced by higher applied normal and shear stresses, thus generating larger changes in surface topography in our experiments at high confining pressures. In turn, the generation of a third body during frictional sliding is expected to highly contribute to (i) the frictional sliding process (Yamashita et al., 2015;Tesei et al., 2017; a review is given in Scholz et al., 1972) and (ii) the fracture transmissivity (Faoro et al., 2009;Tanikawa et al., 2010;Rutter & Mecklenburgh, 2017, 2018. We note that, a more refined model should consider not only the deformations due to shear stress (Park & Song, 2013) but mostly to wear processes and the formation of a third body (Aghababaei et al., 2016;Milanese et al., 2019;Molinari et al., 2018).
After shearing, the axial stress on all samples was decreased (i.e., both normal and shear stress were decreased to isostatic stress conditions) and transmissivities measured again under isostatic stress conditions (σ ′ 1 ¼ σ ′ 3 ). The transmissivities after unloading (kt unl ) slightly increased with respect to the values measured at the end of the shearing stage (last points in Figure 4). Nevertheless, the transmissivity was far from being recovered to the isostatic initial transmissivity (kt 0 ). The difference between kt 0 and kt unl was usually larger than an order of magnitude, meaning that after shearing, the faults hydraulic transport capacity was much lower even if the stress state was the same. This result can be attributed to the production of wear products (thin gouge) which can obstruct fluid flow in the fault, increasing the contact area and reducing the hydraulic aperture at given stress conditions (Faoro et al., 2009;Rutter & Mecklenburgh, 2017, 2018Tanikawa et al., 2010). In our experiments, this is further supported by the evolution of surface roughness maps after shearing (Figures 5 and 6). Similar observations have been made by Molecular Dynamics simulations of frictional sliding (Spijker et al., 2013) where roughness decreased exponentially during sliding in hard metals. In natural faults, this effect is difficult to quantify. On the one hand, shear reactivation can decrease fault transmissivity due to the formation of wear products at small displacements. At large displacements, the wear products can become large enough that the transmissivity is in turn increased depending on fault roughness, sliding conditions (stress, slip velocity, power density, etc …), and rock properties (Milanese et al., 2019;Molinari et al., 2018;Tesei et al., 2017;Yamashita et al., 2015). On the other hand, shear reactivation can increase transmissivity on "healed" faults (e.g., if the fault's core is composed of glassy products and or

10.1029/2020JB019767
Journal of Geophysical Research: Solid Earth consolidated fault gouge), due to porosity unclogging (Elkhoury et al., 2011) or macroscopic dilatancy (Cox, 2010;Crawford et al., 2008;Faulkner et al., 2010;Jeanne et al., 2018;Tesei et al., 2017;Zoback & Gorelick, 2012). It is possible that in faults with large scale roughness and heterogeneity, large displacements (larger than the characteristic wavelength) lead to macroscopic dilation and increased fracture aperture (Chen et al., 2000;Ciardo & Lecampion, 2019;Watanabe et al., 2008Watanabe et al., , 2009. It is noteworthy that natural faults that have seen several slip episodes present grooves in the shear sense, which become more marked with increasing slip (Candela et al., 2009;Sagy et al., 2007). Nevertheless, the periodicity of such grooves might not show a clean periodicity as the one imposed on our initial surfaces (Tesei et al., 2017). Further experimental work dealing with the evolution of fault roughness and its relation to the mechanics of frictional sliding is needed to understand how the production of wear products interacts with mechanical and hydraulic aperture changes.

Discussion and Implications for Geo-Energy Reservoirs
In our experiments, the fault's transmissivity first strongly decreased due to the application of normal stress, the sensitivity of fault's transmissivity to normal stress depended on fracture geometry (Figures 3, 4, 9, and 10). When shear stress was applied, its influence was very small and was conditioned by the interaction between stress and fault geometries. Then, once faults were submitted to shear sliding (up to 1 mm), the transmissivity changes were small (often less than 1 order of magnitude) and its evolution could be predicted by the fault's geometry in experiments at low normal stress only. During the stimulation of deep geothermal reservoirs, for example, one stimulation strategy to increase the production flow rate is to reactivate faults in shear, expecting that the shear displacement permanently increases transmissivity (Breede et al., 2013;Cladouhos et al., 2010;Cornet, 2016). The strategy consists in injecting pressurized fluids into the fractured reservoir in order to decrease the effective normal stress on the fault, leading to reactivation. From the results in our study, we can observe that the reduction in normal stress can indeed increase the fault's transmissivity and that the magnitude of increase will depend on the fault surface geometry. If faults are flat but rough, the transmissivity increase will be smaller than if faults present some kind of macroscopic wavelength and heterogeneous topography. (Watanabe et al. [2008(Watanabe et al. [ , 2009 present somewhat similar observations in granitic rocks.) In addition, decreasing the fault's effective normal stress will increase the shear-to-normal stress ratio (friction); thus, transmissivity should change under the effect of shear stress, depending on the fault's geometry. This remains valid only if the overpressure in the reservoir is somehow maintained during the stimulation and production phases, reducing the reservoir's effective stress.
The usual idea regarding fault reactivation influence on transmissivity is that slip on rough faults generates a permanent increase in hydraulic transport on the fault (Carey et al., 2015;Esaki et al., 1999;Guo et al., 2013;Ishibashi et al., 2012;Lee & Cho, 2002;Olsson & Brown, 1993;Pyrak-Nolte et al., 1988;Wenning et al., 2019;Yeo et al., 1998;Zambrano et al., 2018). In most of those studies, the shear displacements were in the range of 3 to 20 mm and/or on faults with large roughness (Chen et al., 2000;Wenning et al., 2019). Our results show that when shear slip occurred in the fault (<1 mm), transmissivity remained close to constant with very slight changes, in strong contrast to the usual understanding of shear reactivation described above. Very few other experimental studies have reported decrease in hydraulic transport properties with shear displacement (Faoro et al., 2009;Rutter & Mecklenburgh, 2017, 2018Tanikawa et al., 2010). In deep geothermal reservoirs, stimulations that generate large shear displacements are usually associated with the occurrence of large magnitude seismic events. In the light of our results, inducing small shear displacements will not significantly increase the reservoir's fluid transport capacity. However, if the reservoir's fractures and faults are "clogged" (due to the presence of frictional wear products or by mineral precipitation), it is highly possible that shear reactivation, followed by an unclogging treatment (chemical or hydraulic for example), does increase the fractures transmissivities (Elkhoury et al., 2011). This should be studied in future work.
We can conclude that changing the state of stress will improve fluid production in a larger scale than generating significant offsets in the fault (which can be technologically difficult). This can be done for example through stress preconditioning (Fryer et al., 2018(Fryer et al., , 2020. A novelty from our study is the ability to "predict" the shape of transmissivity changes with fault reactivation at low normal stress. Indeed, as shown in section 6.1, it seems that knowledge of fault surface geometry (even though this is a difficult measure to obtain) on faults submitted to small normal stresses can help estimate the shape of transmissivity change with shear displacement. Further work is nevertheless needed to estimate at which stress levels this prediction stops being accurate and to predict the change at high normal stress. Indeed, by preconditioning the reservoir (reducing the effective stress on it; Fryer et al., 2018Fryer et al., , 2020, the predictions of transmissivity should become more accurate. We can speculate that the use of more complex models that include shear stress and wear processes can strongly improve the prediction of transmissivity with shear displacement (Aghababaei et al., 2016;Frérot et al., 2019;Milanese et al., 2019;Molinari et al., 2018;Shvarts & Yastrebov, 2018a, 2018bShvarts, 2019;Yastrebov et al., 2017). The results obtained for a single rough-fault could then be input into discrete fracture models (DFMs) to evaluate the hydraulic transport in fracture networks (McClure & Horne, 2013).

Conclusions
In this study, we developed an experimental technique to customize the surface roughness of real-rock samples, which can be representative of engineered geothermal reservoirs (and more generally of underground carbonate reservoirs; Delle Piane et al., 2015;DiPippo, 2012). The resulting surfaces had a fully controlled geometry which was precisely measured and analyzed using the radially averaged PSD of heights ( Figures 5 and 6). Then, single fractures of Carrara marble with customized roughness were experimentally loaded under deep reservoir conditions (both under normal and shear loading) to study fluid flow through the wavy-rough fractures (Figures 3 and 4).
A numerical procedure was developed to simulate wavy-rough contacting surfaces (section 5 and Figure 7). It was validated with the experimental data, and it allowed isolating the influence of different roughness parameters on fluid flow in fractures under normal loading. Under normal loading, the macroscopic geometry has a strong influence on the fracture's transmissivity decay (here fitted to an exponential decay) (Figures 3 and 8). Surfaces with grooves suborthogonal to fluid flow are more sensitive to the application of normal stress than samples with grooves subparallel to fluid flow and with no grooves in that order (Figure 8). A parametric analysis of the numerical procedure on samples with grooves suborthogonal to fluid flow showed that changes in the Hurst exponent and in the roll-off wavevector alone have little influence on transmissivity (Figure 9). On the other hand, the standard deviation of heights and the macroscopic wavelength have a strong influence on hydraulic transport capacity (Figure 9). The evolution of transmissivity has a strong nonmonotonic dependence on these two parameters ( Figure 10). Similar experimental fractures were loaded in shear and the numerical procedure was then used to isolate the effect of reversible shear loads on the fracture transmissivity (during elastic deformation; Figure 11). The influence of reversible shear load on fault's transmissivity is almost negligible and its magnitude depends on the fault's geometry; thus, it is not straightforward to estimate how reversible shear loading affects fracture transmissivity. Finally, the effect of irreversible shear displacement on transmissivity was studied with support on the numerical procedure. It was found that, increasing shear displacement (until 1 mm) generally decreases fracture transmissivity with slight variations at each displacement step (of 0.1 mm up to 1 mm; Figure 4). In that case, the transmissivity could be roughly predicted by changing the model geometry at low normal stress (Figures 12a-12d) probably because wear was not too prominent. On the other hand, the geometrical model was completely off when normal stress was high during shear reactivation (Figures 12e-12h), probably due to prominent wear and surface topography changes ( Figures 5 and 6).
From this study, we observed that the main controls on fluid flow through rough fractures are the surface geometry and the stress applied on the fracture (normal to the fracture plane). Regarding surface geometry, we observed that the imbrication of the main wavelengths, the h RMS , and magnitude of the macroscopic wavelengths had a much larger influence on fluid flow than the Hurst exponent and the roll-off wave vector of the PSD of heights. We observed that the application of normal stress significantly reduced fluid transport capacity compared to the application of shear stress. Finally, small irreversible shear displacements (<1 mm) did not increase the fluid transport capacity of the rough fractures.
Further work is needed to develop a numerical procedure that allows simulating fractures under both normal and shear stresses. The target model should also include examination of shear induced plastic deformation and wear processes. This would allow coupling the evolution of fault roughness with the hydraulic transport properties as shear reactivation occurs.
In the light of our results, during EGS stimulations, small shear displacements will not significantly increase the reservoir's permeability (unless porosity unclogging occurs). In turn, reducing the effective stress on the reservoir's fractures and faults will generate large permeability increases and improve the ability to predict its evolution. Figure A1 shows the permeabiity of intact Carrara marble as function of effective stress and compares it to existing data from the litterature Figure B1 shows the stress-strain data from an experiment of triaxial deformation of intact Carrara marble.   . Typical stress-displacement curve and the associated transmissivity. The shaded areas represent what is considered as "reversible loading" in the text. The black dots correspond to the transmissivity measurements reported in Figure 11. An interpolated data point is used when shear reactivation occurs at displacements inferior to 0.2 mm, not allowing two transmissivity measurements during reversible loading. In that case, the transmissivity is estimated through a nearest neighbor interpolation technique.

Data Availability Statement
All the data related to this manuscript are uploaded to Zenodo.org (DOI: 10.5281/zenodo.3935857).