The Stress‐Memory Effect of Fracture Stiffness During Cyclic Loading in Low‐Permeability Sandstone

The hydraulic performance and mechanical stability of open fractures are crucial for several subsurface applications including fractured geothermal reservoirs or nuclear waste repositories. Their hydraulic and mechanical properties (fluid flow and fracture stiffness) are both strongly dependent on the fracture geometry. Any change in effective stress impacts aperture and thus the ability of fractures to promote flow. Here, we carried out flow experiments with shear displaced tensile fractures in pre‐loaded, low‐permeability sandstones with two different cyclic loading scenarios with up to 60 MPa hydrostatic confining pressure. During “constant cyclic loading” (CCL) experiments, the fracture was repeatedly loaded to the same peak stress (up to 60 MPa). During “progressive cyclic loading” (PCL) experiments, the confining pressure was progressively increased in each cycle (up to 15, 30, 45, and 60 MPa). The matrix and fracture deformation was monitored using axial and circumferential LVDT extensometers to obtain the fracture stiffness. The fracture geometry before and after the experiment was compared by calculating the aperture distribution from 3D surface scans. Initial loading with confining pressure of the fracture leads to a linear fracture specific stiffness evolution. For any subsequent stress cycles fracture stiffness shifts to a nonlinear behavior. The transition is shown to be related to a stress memory effect, similar to the “Kaiser Effect” for acoustic emissions. Progressive loading of fractures possibly leads to less permeability reduction compared to continuous cyclic loading.

The memory of rocks is described as the capacity of rocks to retain "imprints" from their stress history (Lavrov, 2005). Rocks contain crucial information about the stress history during non-elastic deformation. Loading a rock to a large stress level generates damage or microcracks within the rock. This level of stress can be identified by reloading the rock above the previous stress level (Lockner, 1993), while monitoring acoustic emissions. Acoustic emissions will be present when exceeding the "ancient" stress level. This is known as the "Kaiser Effect" (Kaiser, 1953). It is important to understand if other mechanisms than the generation of microscopic fractures in an intact rock (Holcomb, 1993) lead to a stress memory effect (Lavrov, 2003). In particular, we try to assess here if existing fractures could contain information about their stress history since they behave non-elastically over a larger range of stress owing to the damage of their asperities (Bandis et al., 1983). Bandis et al. (1983) evaluated the mechanical behavior of fractures during cyclic loading. Since then, a large number of data have been published that study in detail the evolution of fracture closure and permeability during cyclic normal loading (e.g., Z. Chen et al., 2000;Cook, 1992;Hofmann et al., 2016;Kluge et al., 2020;Milsch et al., 2016;Watanabe et al., 2009). The loading path in most of these experiments comprises multiple loading cycles to the same peak load. Such experiments allow for a qualitative description of the changes in physical properties by a repetition of the same loading path. This is important in understanding the performance of a reservoir during different injection and production scenarios . Cyclic loading experiments contain more information when varying the stress path of each loading cycle. The pore fluid pressure oscillations technique can be used to analyze the frictional stability of the fault and to analyze the transition from stable to unstable slip by progressively increasing the magnitude of pore fluid pressure from one cycle to another (Noël et al., 2019). This technique can also be implemented by loading the sample with an increasing external confining pressure with the same differential stress through the different cycles, instead of increasing pore fluid pressure at a constant confining stress.
The hydro-mechanical properties of fractures depend on fracture contact area, fracture size, fracture roughness, and loading stress history (Wang & Cardenas, 2016). These physical factors also control the stiffness of a fracture (Pyrak-Nolte & Morris, 2000). The normal stiffness describes how much a fracture closes when being subjected to an increasing normal load with respect to the normal plane. Generally, normal fracture stiffness evolves exponentially with load (Bandis et al., 1983). Two different models are usually applied to characterize fracture stiffness from laboratory experiments: (a) The fracture stiffness characteristic, , which is a parameter that is based on the semi-log closure model for a single loading cycle (Zangerl et al., 2008). It can be used to quantify changes of stiffness in a series of repetitive stress cycles (e.g., Bandis et al., 1983;Crawford et al., 2016;Kluge et al., 2020). In particular, it is useful to describe strain-hardening effects of fractures. However, this parameter is strongly based on a specific model, the semi-log model and can be poorly assessed if the model is not correctly describing the behavior of fractured rock. (b) The specific fracture stiffness,  E , is defined as the ratio of the increment of stress to the increment of displacement caused by a change of the void space in the fracture (Pyrak-Nolte & Morris, 2000). It was shown numerically and experimentally, that it depends on the elastic properties, the fracture geometry, and stress history (e.g., Cook, 1992;Marache et al., 2008;Petrovitch et al., 2013;Pyrak-Nolte & Morris, 2000;Wang & Cardenas, 2016). This property enables to monitor dynamic changes of fracture stiffness for complex stress paths.
Similar to the fracture stiffness, fracture permeability also depends on these parameters. Cyclic loading experiments on hydro-mechanical responses of fractured rocks have a long history and have provided a large number of data (e.g., Z. Chen et al., 2000;Hofmann et al., 2016;Milsch et al., 2016;Watanabe et al., 2009). These studies have focused on a repetition of periodic loading cycles with the same peak stress. Studies with aperture or permeability measurements during non-periodic cyclic loading tests where the cycle's maximum stress is increased from cycle to cycle gained little attention in the past (Lavrov, 2005). Bandis et al. (1983) and Pyrak-Nolte and Morris (2000) showed that larger apertures and therefore more permeable fractures were more compliant than fractures with initially lower aperture and permeability. The question is, whether the mechanical (stiffness) and hydraulic properties (permeability) of fractured rocks are dependent on stress history and if stress cycling alters this relationship (Pyrak-Nolte & Morris, 2000;Pyrak-Nolte & Nolte, 2016). The fracture closure and stiffness are also expected to depend on fracture surface roughness (Akarapu et al., 2011;Persson, 2007). This can only be shown by means of the specific fracture stiffness,  E .
The deformation of asperities at the fracture surface may lead to changes in the fracture topography and consequently the aperture. Previous studies showed difficulties to quantify changes of fracture topography (e.g., Bandis et al., 1983;Vogler et al., 2016;Xia et al., 2003;Yoshioka, 1994;Zou et al., 2020). Further, it is not clear if potential topography changes affect the self-affine scaling properties, such as the power spectral density (Schmittbuhl, Schmitt, & Scholz, 1995;Schmittbuhl, Vilotte, & Roux, 1995), of fracture surfaces. A combination of measuring the fracture stiffness and permeability during cyclic loading experiments with a progressively increasing stress magnitude together with measurements of the fracture roughness evolution might enable to better understand the dependency of these properties and their relation to the stress history.
In this paper, we present results from a set of laboratory experiments on fractured rock samples with a single displaced tensile fracture, being cyclically loaded using two different loading scenarios: constant cyclic loading (CCL) and progressive cyclic loading (PCL). We will first review the experimental results and analyze the fracture stiffness and fracture permeability evolution. These will be discussed in respect of a possible memory effect of fracture stiffness, similar to that in intact rock during plastic deformation. We then elaborate how this impacts the relationship of the hydraulic-mechanical properties of a fracture under cyclic loading conditions and investigate possible fracture deformation mechanisms.

Testing Equipment
The flow-through experiments were carried out in a conventional MTS tri-axial compression cell. The stiff, servo-controlled loading frame (MTS 815, Material Testing Systems Corporation) holds a loading capacity of up to 4,600 kN (load cell calibrated to 2,000 kN, calibration error <1%) and a servo-controlled maximum hydrostatic confining pressure of 140 MPa applied via an oil-filled pressure vessel coupled to an external pressure intensifier. The pore fluid pressure was applied via four Quizix fluid pressure pumps (Model C6000-10K-HC-AT) with a maximum fluid pressure of 69.7 MPa. The differential fluid pressure, which is the difference between in-and outflow pressure, was measured using a differential pressure transducer (GP:50, Model 215; range: 1 MPa; line pressure max. 69.7 MPa; precision:  E 1%). The circumferential strain was measured using an LVDT extensometer chain and the axial strain was measured with two axial LVDT extensometers ( Figure 1). All experiments have been performed at a temperature of 30 C. Temperature was controlled via heat stripes on the outer side of the loading vessel. During confining pressure increase and decrease the temperature of the confining oil increases and decreases, respectively, due to the Joule-Thomson effect. Data were recorded at a frequency of 1 Hz. A detailed description of the machine can be found in Pei et al. (2016). we calculated the sample permeability, s E k , using Darcy's law (Darcy, 1856): Once permeability was measured, the constant inflow rate was changed to a constant pressure mode. By applying a constant fluid pressure of 0.2 MPa at both sides of the sample using only one pump, we measured the pore volume changes during loading. For pre-conditioning, a hydrostatic stress of 65 MPa was applied in a total of seven pressure cycles with a loading rate of 5 MPa/min ( Figure 2). The effective pressure resulting from the confining and pore pressure was calculated following Terzaghi's effective pressure law (Terzaghi, 1925), assuming a linear pressure gradient In six of the cycles, the pore volume change was measured, while the sample permeability (Equation 1) was measured at several hydrostatic pressure levels in a seventh cycle. This ensured the stress-strain curves to where ( ) E P k is the Fourier power spectrum, E k is the wave number, E C is a pre-factor, and E H is the roughness exponent. Since higher frequencies are over-represented in log-log-space, the mean spectra were re-sampled to 20 points (frequency) averaging the data in between for the linear fit. This way, all frequencies were evenly represented providing a better fit. The method was verified on a synthetic fault generated for two dimensional roughness exponents.
To obtain the aperture distribution of the top and bottom surface of one rock sample, the point cloud topography data of both surfaces were correlated. To calculate the aperture distribution of two independently scanned surfaces, both surfaces need to be matched. This was done by matching the best fitting principal planes of the bottom and top surface and applying a grid search algorithm. The surface data of both surfaces were interpolated on a grid with a point distance in x-y-direction of 0.05 mm. The top and bottom surfaces have the same orientation and shared the same grid. The fracture surfaces were displaced to an offset of 0.5 mm, as in the experiment. At each point across the x-y grid, the aperture (vertical distance) between the top and bottom surface was calculated to obtain the aperture distribution.

Procedures for the Fracture Flow Experiments
After generating and analyzing the fracture and its geometry the samples were prepared to perform the fracture flow experiments under the respective loading scenario (constant cyclic loading (CCL) and progessive cyclic loading (PCL)). First, the two samples halves were placed together at a manual offset of 0.5 mm using perforated steel spacers at the opposite side of each sample half ( Figure 1a). Any rotation of the two fracture planes can be ruled out due to the parallelism of the end cap and the spacer. The rigid shear offset of 0.5 mm was chosen based on three criteria: (a) larger than the grain size, (b) a comparable initial mean aperture for all three samples and (c) not too large to be able to measure the sample permeability (Equation 1) with our testing equipment (between 1   18 2 10 m E and 1   12 2 10 m E ). A brass stripe was used to cover the resulting holes caused by the spacers and the displaced fracture to avoid the heat-shrink tube to be punctuated ( Figure 1d).
After the sample was installed in the triaxial cell, the confining pressure was increased to 2 MPa with a loading rate of 0.5 MPa/min. The fractured sample was then saturated for 24 h under vacuum conditions with a constant pore fluid pressure of 0.2 MPa. After the saturation, a constant inflow rate of 2.5-10 ml/min was applied from one side of the sample, while the outflow pressure was kept constant at 0.2 MPa, resulting in an effective pressure according to Equation 2.
The sample permeability was measured over the entire duration of the experiment. We assumed that measured flow rate is the sum of the individual flow rates through the matrix and through the fracture Here, h E a is the hydraulic fracture aperture, s E k is the measured sample permeability, and E r is the sample radius. From the hydraulic aperture we then calculated the fracture permeability, f E k , using the cubic law (Witherspoon et al., 1980): Permeability errors were marginal and can hardly be quantified. The main error source was the frictional pressure losses within the capillary tubes connecting the fluid pumps and the sample, potentially leading to slight pressure changes. A second error source was the accuracy of the differential pressure transducer with an absolute error of 1% of the pressure range (1 MPa), which corresponds to only 1 kPa.
A circumferential LVDT extensometer chain was attached to the center of the sample to measure the bulk strain of matrix and fracture, c E  . To remove the non-elastic strain from the total strain signal, we subtracted the fitted elastic strain of the intact rock sample, , The fracture specific stiffness,  E , was defined as the ratio of the increment of stress to the increment of displacement caused by the deformation of the void space in the fracture. The fracture stiffness,  E , was calculated from the change in fracture closure, E a a n a n , per increment effective pressure increase, ( ) e E p n averaged over an interval of 720 s (6 MPa): p n p n a n a n (11) The interval of averaging only affects the smoothness of the signal to be able to better illustrate the stiffness evolution with changes in stress. Shorter intervals lead to noisy data while longer intervals lead to smoother data.

Constant Cyclic Loading (CCL) Experiment
The CCL experiment at hydrostatic conditions was performed according to Figure 2a, using sample SBT6-BE-04-03. The confining pressure was increased from 2 to 60 MPa at a constant loading rate of 0.5 MPa/ min. At 60 MPa, the confining pressure was held for 20 min before unloading was started at the same rate as during loading. This process was repeated six times to obtain the fracture closure curves for a constant loading procedure with the same peak stress. After two cycles, the system was held at a constant low-level load for about 12 hr.

Progressive Cyclic Loading (PCL) Experiments
The PCL experiments at hydrostatic conditions were performed according to Figure 2b, using samples SBT6-BE-04-09 and SBT6-BE-04-10. In this procedure, we distinguished between pressure cycles and pressure sub-cycles. The sub-cycles describe the stepwise increase of hydrostatic confining pressure from 2 to 15, 30, 45, and 60 MPa. One sub-cycle was therefore the increase from 2 to the respective stress level (15, 30, 45, and 60 MPa) and the decrease or unloading from the respective stress level to 2 MPa. The confining pressure was held for 20 min before unloading was started at the same rate as during loading. That way, the effective pressure was increased by an additional 15 MPa over the previous pressure level to identify potential changes in the fracture closure or opening behavior when exceeding the previous stress level. This progressive loading was repeated two times, that is, two complete cycles, with a hold phase of about 10 hr between the two. The loading rate for the confining pressure was 0.5 MPa/min for loading and unloading.
At the end of the experiments, the samples were removed from the cell and dried at 50 C for at least 24 hr. The fracture surfaces were scanned and analyzed as described in Section 2.3.3.

Sample Volume and Permeability of the Intact Rock
Based on the bulk volume change (Equation 4), the plastic and elastic sample deformation was monitored during six loading cycles up to 65 MPa (Figures 3a-3c). Most of the irreversible plastic deformation was found in the first loading cycle (Figures 3a-3c). After six pressure cycles, no more significant bulk strain changes suggested fully elastic sample deformation. The total volume loss was about 0.27 3 cm E measured by the volumetric strain and 0.39 3 cm E measured by the pore volume change. This corresponded to about 0.1% of the initial porosity of 5%-6%. The seventh loading cycle is not shown here, because the fluid pressure in the sample was higher during the permeability measurements.
In the seventh loading cycle, the intact rock permeability (Equation 1) was measured at four to five pressure levels during loading and unloading ( Figure 3d). The sample permeability, s E k at 2 MPa was reduced from about 1 to 3   15 2 10 m E before pre-conditioning down to 5   18 2 10 m E -1   17 2 10 m E after pre-conditioning. This corresponds to a difference of more than two orders of magnitude. The largest incremental change in permeability was found during the first 15 MPa of confining pressure, while the change becomes smaller after 40 MPa. At effective pressures larger than 40 MPa the permeability was about 1   19 2 10 m E , which was the lower limit of measurable permeability with our machine. The permeability of all measured samples showed a similar behavior. During loading, the permeability was overall higher compared to unloading and the change in permeability was always reversible after pre-conditioning.

Tensile Fracture Generation
The tensile strength, E TS , calculated from Equation 5, for the three samples SBT6-BE-04-03, 09, and 10 was 5.4, 4.0, and 5.5 MPa, respectively. The mode I fracture toughness, Ic E K , was calculated after Guo et al. (1993) and corresponded to 0.78, 0.63, and 0.82 MPa⋅m 1/2 , respectively. The sample halves were taken apart, loose fragments were carefully removed and possible breakouts at the corners were filled with an Araldite-sand mixture. This was done to avoid inward bulging and rupture of the heat-shrink tube during loading. However, it did not lead to additional contact-area between the opposing fracture surfaces and had thus no impact on the strength or stiffness of the fracture.

Fracture Permeability and Aperture During CCL
Comparing the initial sample permeability after pre-conditioning (Figure 3b), to the initial sample permeability containing a fracture at 2 MPa effective pressure (Equation 1), we found an increase from 5.4   18 2 10 m E to about 6.0   13 2 10 m E . In the following, we will refer permeability to the fracture permeability calculated following Equation 8. During the first loading cycle from 2 to 60 MPa of the confining pressure, the permeability reduced from 3.5   10 2 10 m E to 2.9   12 2 10 m E (Figure 4a). While the permeability changes were only minor during the first 20 MPa, the largest permeability decrease was observed above a confining pressure of 20 MPa. During the hold phase, the permeability reduced slightly but reached an almost constant level within 20 min (Figure 4b). The permeability recovery when reducing the effective pressure was slower than the permeability loss during the effective pressure increase. When reaching 2 MPa, the permeability loss was more than one order of magnitude with a permeability of about 2   11 2 10 m E at the end of the first cycle. During the second loading cycle, the incremental permeability decrease was larger at lower stress,  while only minor permeability losses were observed at an effective pressure above 30 MPa. The minimum permeability at 60 MPa in the second cycle was 1.9   12 2 10 m E . In the following cycles, the trend of the second loading cycle was maintained and a reversible fracture permeability was observed. There was almost no permeability reduction during the overnight hold phases (Figure 4b, inset). Smaller and short-term peaks in permeability were related to changes in the flow rate during loading and unloading.
The maximum elastic matrix deformation measured with the axial extensometers was no more than 0.14 mm at 60 MPa ( Figure 5a). Subtracting this from the total radial deformation of about 0.53 mm after the first loading cycle, the actual fracture closure was roughly 75% of the total measured deformation resulting from six cycles. While the matrix compaction was fully elastic within the six loading cycles, the residual fracture closure was 0.28 mm after the first loading cycle at 2 MPa. The following cycles showed further fracture closure by up to 0.37 mm. The incremental fracture closure with increasing stress was rather linear at effective pressures larger than 20 MPa. During unloading, the fracture remained closed until about 10-15 MPa (Figure 5b). This was similar to the permeability data, which shows a nearly constant permeability until about 10-15 MPa.

Fracture Stiffness Evolution During CCL
The fracture stiffness was calculated following Equation 11 at intervals of 720 s. This was done to reduce the noise related to the sensors only with impacts on the signal quality of the fracture stiffness. The results were separated to first calculate the fracture stiffness during the successive loading cycles (Figure 6a). In the first loading cycle, the fracture stiffness increased linearly from about 80 MPa/mm to about 750 MPa/ mm at 60 MPa. The 2 nd E -6 th E loading cycles showed a nonlinear but reversible fracture stiffness evolution. Only minor increases with progressive loading cycles were found, smaller than noise measurements. This reversible fracture stiffness from the second loading cycle was higher compared to the first loading cycle, with values between 450 MPa/mm at 2 MPa and up to 2,300 MPa/mm at 60 MPa. The peak stiffness was reached at about 40-50 MPa with a slight decrease toward the final stress level. During unloading, the fracture displacement (Figure 5a) was almost constant, that is, the fracture remained closed. Therefore, fracture stiffness was undetermined until a stress of about 10-15 MPa.

Permeability and Aperture During PCL
In the following, we refer to loading sub-cycles as the pressure change from 2 to 15 MPa, 2 to 30 MPa, 2 to 45 MPa, and 2 to 60 MPa, as well as a full loading cycle, which is a complete set of the four sub-cycles (Figure 2). We performed a total of two complete loading cycles for samples SBT6-BE-04-09 and 10.
The fracture permeability was first measured at 2 MPa confining pressure with values of about 3.7-3.8   10 2 10 m E (Figures 4b and 4c). This was close to the highest permeability that can be measured in our testing device. The flow rate was set to a maximum of 10 ml/min to avoid turbulent flow conditions. Permeability reduced after every loading cycle, with a larger reduction with increasing pressure in the respective sub-cycle. The permeability reductions after each sub-cycle during the first cycle to 15, 30, 45, and 60 MPa correspond to a factor of 1.2, 1.3, 1.6, and 2.2 compared to the initial permeability. During re-loading, the permeability follows the same permeability trend as during unloading up to 20 MPa. At stresses higher than 20 MPa, the permeability was lower during loading compared to unloading. After the first loading cycle, a permeability loss by a factor of 5.2 was measured for both samples. This corresponds to about 3.5-3.7   10 2 10 m E at the beginning and 0.5-0.7   10 2 10 m E at the end of the first cycle. During the hold phases, a continuous reduction of permeability was observed at stresses above 30 MPa, with an increasing reduction at increasing stress levels (Figures 4d and 4f). In the second loading cycle, the permeability was almost reversible, following the trend during unloading up to about 45 MPa, a further reduction was found. The permeability reduction during the hold phases in the second cycle was minor compared to the initial loading cycle.
The elastic matrix deformation was roughly 30% of the bulk deformation measurement. Therefore, about 70% of the total deformation was related to fracture closure (Figures 5c and 5e). When re-loading the fracture, the aperture follows the unloading path until reaching the previous stress level. After that, the fracture closure was larger, meaning that the slope of the fracture closure versus the effective pressure curve was shallower. This trend was similar in all sub-cycles up to 60 MPa. The total closure in sample SBT6-BE-04-09 and 10 at the end of the second loading cycle was 0.37 and 0.35 mm. During unloading, the fracture remained closed until about 10 MPa (Figures 5d and 5f). This was similar to the permeability data, which shows a nearly constant permeability until about 10 MPa.

Fracture Stiffness Evolution During PCL
The fracture stiffness results during the respective sub-cycles are summarized in Table 2. Both PCL experiments (SBT6-BE-04-09 and 10) showed the same trend and magnitudes of fracture stiffness. The effective stress was increased from 2 to 15 MPa during the first sub-cycle. This resulted in a linear increase in fracture stiffness with increasing effective pressure from around 80 to about 190 MPa/mm (Figures 6b and 6c). During the 2 nd E sub-cycle, the pressure was increased from 2 to 30 MPa. The fracture stiffness initially followed the nonlinear fracture stiffness curve of the 2 nd E to the 6 th E of the CCL experiment (SBT6-BE-04-03), starting from around 500 MPa/mm. Fracture stiffness decreased shortly before reaching the effective stress of 15 MPa. The curve returned back to the linear fracture stiffness curve of the first loading cycle of the CCL experiment. This resulted in a fracture stiffness of around 350-400 MPa/mm at 30 MPa for both samples. During the third sub-cycle, the stress was increased from 2 to 45 MPa. The fracture stiffness followed the nonlinear trend of the 2 nd E to 6 th E loading cycle of the CCL experiment, as well as the second sub-cycle recorded before that. Both samples kept following this trend exceeding a pressure of 15 MPa, but fracture stiffness started to decrease shortly before reaching the previous stress level of 30 MPa. Again, the fracture stiffness returned to the linear trend of the first loading cycle of the CCL experiment, reaching an end value of about  sub-cycle of the PCL experiment. The fracture stiffness deviated from the nonlinear fracture stiffness shortly before reaching 45 MPa, but the reduction was smaller compared to the clear drop of the second sub-cycle before 15 MPa. The fracture stiffness almost reached the end value of about 1,000 MPa/mm of the first loading cycle of the CCL experiment at 60 MPa. At stresses larger than 50 MPa we found less changes in fracture stiffness with increasing stress. The smaller reduction in stiffness when reaching the previous stress level of 45 MPa indicated less fracture closure with increasing stress. The four sub-cycles were repeated in a second cycle of the PCL experiments, increasing pressure from 2 to 15, 30, 45, and 60 MPa, respectively. During all subsequent loading sub-cycles, the fracture stiffness followed a similar nonlinear trend as during the 2 nd E -6 th E loading cycle of the CCL experiment.
Similar to the CCL experiment, the fracture displacement (Figures 5b and 5c) was almost constant during unloading until about 10 MPa, resulting in an undetermined fracture stiffness.

Fracture Roughness Exponent
The roughness exponent was determined using the power spectral density approach as explained in Section 2.3.3. Assuming that the tensile fracturing in sandstones was exclusively intergranular, the higher frequencies reflect the surface roughness of the grains (Boffa et al., 1998). Those frequencies led to deviations in the power spectrum. They were suppressed by a frequency cutoff for length scales that are two times the grain size, that is, two times 0.2 mm. The resulting cutoff frequency additionally marks a break in the slope and a deviation from a linear trend in log-space. Figure 7 shows the power spectrum for the surfaces before and after the flow-through experiments in x-y-directions. The slope of the fitted trend eventually leads to the roughness exponent (Equation 6).
The power spectra indicate a similar roughness exponent for all surfaces independent of the direction (Table 3). We found a mean of 0.58 ( E 0.02) in x-direction and 0.57 ( E 0.02) in y-direction (shear direction) before the experiment. These values are in agreement with roughness exponents of around 0.5-0.6 for sandstones that are commonly slightly lower than the 0.8 found for most other rock types (Boffa et al., 1998). After the cyclic loading experiment, the surface topography was obtained a second time in the same orientation as before. The post-experimental mean values for  the roughness exponent in the x-y-direction were 0.57 ( E 0.01) and 0.58 ( E 0.03), respectively. This indicated no distinct change in the scaling properties of the surface roughness due to the cyclic loading.

Fracture Aperture Distribution and Contact-Area Ratio
From the surface topography from each fracture surface, we calculated the aperture distribution by matching the top and bottom surfaces as explained in Section 2.3.3 (Figures 8a, 8c and 8e; left). The aperture, E a , was taken as the distance between each point across a x-y-grid (point distance 0.05 mm). This was done to calculate the initial aperture distribution of every sample. We considered this aperture distribution as the initial aperture, ini E a , at zero stress. With the assumption of two interpenetrating surfaces under normal load, that is, geometrically overlapping regions are assumed to be in contact without deformation (the "overlap" model; Pei et al., 2005), we calculated the resulting evolution of contact-area ratio, c E R . It is commonly defined as the ratio of the surface area in contact, c E A and the total surface area, t E A : Considering only one contact point at zero stress would lead to an overestimation of the mean aperture when small fragments protrude from the fracture surface. Similar to Wang and Cardenas (2016), we defined a threshold to shift the normal distribution to the left and reduce the aperture. The two fracture surfaces were brought into contact at an initial contact-area ratio, c E R , of 0.1%. Furthermore, we consider the contact-area (zero aperture) as a discontinuity (delta function) in the aperture distribution (Pyrak-Nolte & Morris, 2000) and excluded these values when calculating the mean aperture.  (Table 4). There was a slight trend in the data indicating that the higher the initial mean aperture, the higher the fracture closure at the end of the second complete loading cycle. The maximum contact-area ratio, , c max E R , at the maximum effective stress at the end of the second loading cycle, was between 16% and 20% for all samples (Table 4). The contact-area, that is, zero apertures were marked in red in Figures 8a, 8c and 8e (center). Here, we observed that the layering of the sample (perpendicular to the shear-direction), was visible in the contact-area distribution leading to "contact bands" along the bedding. This possibly resulted in smaller necks for fluid flow at large stress. The contact points were predominantly distributed along the edges of the sample indicating a generally concave shape of the fracture surface. This might result from a combination of the tensile fracture generation during diametrical loading conditions deforming the fracture surface and the finite size effect of the samples.
After the experiment, the mean aperture at no confining pressure was between 0.31 and 0.34 mm for all samples, which corresponds to a total permanent aperture reduction of about 0.09-0.12 mm (Figures 9a, 9c and 9e; right). A kink at the peak of the normal distribution in the post-experimental aperture distribution indicates changes in the fracture topography (Figures 9b, 9d and 9f). Due to the anisotropic deformation, caused by one-sided "contact bands" perpendicular to shearing, we observe a cutoff of the mean aperture. This process was possibly more pronounced in the CCL sample compared to the PCL samples. It was not possible to relate the statistical aperture distribution to differences in measured permeability. Permeability not only relies on the aperture distribution but also on the spatial correlations in the aperture (self-affine property). For the same aperture distribution, with or without spatial correlations, the fracture will have a different permeability.
The mean aperture with increasing shear offset before and after the experiment was compared and is shown in Figure 9a. Here we found a mechanical imprint after cyclic loading, meaning that the mean aperture at zero offset was higher and reduced toward the given shear offset applied during the experiment (0.5 mm). The lowest mean aperture was found at a shear offset of 0.35 mm. At a shear offset above 1 mm, the mean apertures after the experiment are similar to those before the experiment. Although no change in roughness above grain scale was found, a change in fracture topography was visual when comparing the mean aperture with increasing shear displacement.

Fracture Stiffness Evolution and the Stress Memory Effect
In this section, we discuss the evolution of fracture stiffness with increasing stress. We address the debate about the linearity of fracture stiffness with increasing effective pressure, as well as the hysteresis effect during loading and unloading. We then relate our findings to a possible "memory effect" of fracture stiffness during the progressive cyclic loading (PCL) experiments.
The loading and unloading path of fracture closure during loading shows a hysteresis effect ( Figure 5). This is well known (Bandis et al., 1983;Brown & Scholz, 1986;Cook, 1992;Pyrak-Nolte, 1987;Skurtveit et al., 2020;Thörn et al., 2015;Yoshioka, 1994;Zou et al., 2020). During cyclic loading, hysteresis decreases and consequently the displacement between cycles decreases (Bandis et al., 1983;Brown & Scholz, 1986;Pyrak-Nolte, 1987). We see the same behavior in our CCL experiments. During unloading, the fracture remains closed at higher stress and opens only at a stress below 10 MPa, even though applying the same pressure rate of 0.5 MPa/min. This leads to larger and permanent fracture closure magnitudes, especially during and after the first loading cycle.
The fracture stiffness magnitudes were similar for all three experiments ( Figure 6). We assume that the data of all three experiments are repeatable and that the experimental workflow and boundary conditions led to consistent data. However, small deviations in fracture geometry between the samples can cause large deviations in fracture stiffness (Pyrak-Nolte & Morris, 2000). The measured fracture normal closure is largely dependent on the position along the sample. Furthermore, the measured values are dependent on local variations caused by local aperture and contact-area variations (Cook, 1992;Marache et al., 2008). Our calculated stiffness values were in the range of 80-3,000 MPa/mm for a sample scale of 100 mm, similar to experiments at the same effective pressure ranges with sandstone reported by Y. Chen et al. (2017) and Skurtveit et al. (2020).
From the CCL experiments, we found a stress-path-dependent fracture stiffness with a linear trend during the initial loading phase and a nonlinear but reversible trend for all subsequent loading cycles (Figure 6a). Previous studies reported contrasting results: while most authors describe a linear relationship of stiffness and stress with different slopes for different stress magnitudes (Akarapu et al., 2011;Bandis et al., 1983;Cook, 1992;Persson, 2007;Pyrak-Nolte, 1996;Wang & Cardenas, 2016;Zou et al., 2020), some reported a partly nonlinear increase of fracture stiffness (Cook, 1992;Raven & Gale, 1985;Pyrak-Nolte, 1987;Pyrak-Nolte & Morris, 2000). From our data, we see that linearity or nonlinearity of fracture stiffness is not trivial, but depends on the stress history of a fracture. During initial loading, the stiffness trend with increasing stress is linear up to at least 60 MPa. The linear behavior can be caused by multiple rheologies, such as elastic, plastic, and elasto-plastic (Greenwood & Williamson, 1966;Kling et al., 2018;Pei et al., 2005;Persson, 2006;Zou et al., 2020). At higher stresses and depending on the rock type, roughness and host rock properties, a change in slope at a certain stress level during initial loading might be possible (Wang & Cardenas, 2016). This linear behavior is not reversible when re-loading the sample within the same range of stress. During these subsequent loading cycles, the system becomes nonlinear. The non-reversible behavior clearly indicates plastic effects. The nonlinear fracture stiffness trend is characterized by an initially steep increase at effective pressures up to 10 MPa and being reduced before reaching the previous maximum stress level for sandstones. Several repeated loading cycles lead to a slight increase in the nonlinear stiffness trend (Figure 6a). Therefore, an additional contribution due to the plastic component during repeated loading cannot be excluded. Since the maximum stress in both experimental scenarios (CCL and PCL) was 60 MPa, we see a slight reduction of fracture stiffness toward that maximum. Further increase in stress beyond 60 MPa might lead to a return of stiffness toward the linear trend of the first loading cycle. It is not clear at what stress the stiffness might lead a limiting value. The initial fracture stiffness during unloading, that is, the constant fracture closure, indicates a permanent aperture reduction with only a minor recovery at stress levels below 10 MPa (hysteresis effect).
When exceeding the previous stress level during re-loading, however, the nonlinear fracture stiffness trend returns to a linear trend (Figures 6c and 6e). This behavior was observed in our PCL experiments. During unloading the fracture remains closed. When exceeding the previous stress, the change from nonlinear to linear fracture stiffness behavior can be repeated. We conclude, this effect is similar to the "Kaiser Effect" (e.g., Kaiser, 1953) and reveals a stress-memory effect of fracture stiffness during progressive cyclic loading (PCL). Figure 9b summarizes the fracture stiffness evolution of the CCL and PCL experiments for stresses of up to 60 MPa. It is not clear whether fracture stiffness approaches a limiting value independent of the number of cycles at higher stresses. The turning point from nonlinear to linear was visible already before the previous peak load is reached. This is similar to the classic "Kaiser Effect" in uniaxially loading tests with intact rock while monitoring acoustic emissions (Lavrov, 2005). Additionally, the fracture stiffness reaches its initial path with some delay. Lavrov (2005) argued that the stress-memory effect may decay in the course of time, that is, when the time interval between successive loading cycles is increased. Whether the stress-memory effect of fracture stiffness decays over time cannot be seen from our data. This is because the time frame of the experiments was too short. We could show that the stress-memory effect is measurable using saturated samples. Experiments by Lavrov (2003) showed that a change in moisture is critical using acoustic emission when trying to detect the "Kaiser Effect."

Relationship of Mechanical and Hydraulic Properties
The fracture stiffness describes the amount of fracture closure with increasing normal stress and therefore directly affects the hydraulic properties of fractures.
Pyrak-Nolte and Morris (2000) described that fluid flow and fracture-specific stiffness are implicitly related since both depend on the size and spatial distribution of aperture and contact-area, or more generally, the fracture geometry. Additionally, stiffness is not only dependent on stress magnitude, since all of the fractures they tested appeared to behave differently, such that any interrelationship among the fracture properties was obscured. They related this to the formation of new contact area as a direct function of the aperture distribution affecting the fracture normal closure. Albeit no relation to the stress magnitude was found, it was concluded that stiffness is dependent on the stress path. Our data support this assumption, while we also observed similar stiffness values at similar stress states. This is possibly due to the accurate sample selection from one block and the resulting compatibility of the three experiments. Contrary, the variety of trends shown by Pyrak-Nolte and Morris (2000) can be caused by a larger variety in fracture geometries of natural fractures. Attempts to normalize the relation of fracture stiffness and permeability were made by Pyrak-Nolte and Nolte (2016) based on numerical simulations. Unfortunately, the required scale-dependent fracture stiffness and permeability cannot be derived from the bulk measurements we obtained in the laboratory.
The initial large decrease in permeability can be explained by a change of flow regime from sheet-like to channelized flow. In this "percolation regime," the contact-area remains unchanged (Cook, 1992;Pyrak-Nolte, 1987;Zimmerman, 2008) and a residual fracture permeability is observed . It implies that permeability becomes increasingly independent of stress (Pyrak-Nolte, 1987;Petrovitch et al., 2013). Progressive cyclic loading leads to a higher residual fracture permeability compared to constant cyclic loading. This highlights the fact that permeability is also a stress-path-dependent property for this type of rock.
Both, CCL and PCL lead to hysteresis effects in permeability, especially between the first and second loading cycles. We could show that after the first complete loading cycle up to 60 MPa, permeability was permanently reduced, that is, the fracture did not fully re-open. In all following loading cycles, permeability also showed hysteresis effects. Permeability was always reduced to about the same value at the lowest applied normal load. Such a behavior was shown already (e.g., Z. Chen et al., 2000;Hofmann et al., 2016;Milsch et al., 2016;Pyrak-Nolte & Morris, 2000;Watanabe et al., 2009). Applying a progressive loading procedure shows a surprising behavior. When reloading the fracture, permeability reduces less than during initial loading and permeability starts to decrease more when the previous peak load is exceeded, as expected. When comparing the total permeability reduction after the second cycle, progressive loading leads to less permeability reduction compared to constant loading conditions. The main difference in permeability evolution was found at the effective pressure of 20 MPa. We believe that permeability is affected by the allocation of crushed asperities. During constant cyclic loading (CCL), the crushed material remains relatively in place, blocking potential fluid pathways and therefore continuously reduces permeability. During progressive cyclic loading (PCL), the stress is released after reaching a certain stress level before increasing the stress further. That way, fines are potentially flushed out of the sample or allocated to larger fracture void spaces with low flow velocities when reducing the load to a minimum. Since the aperture is self-affine with an roughness exponent of about 0.6 (it is even stronger for H = 0.8), open spaces are extending on larger scales than for a random surface with no self-affine properties (negative H), and accordingly channeling is stronger. This possibly leads to a higher permeability compared to the continuous loading process. This can only be shown by analyzing the effluent or analyzing the fracture morphology before and after the experiments and therefore requires further studies. The resolution of the surface scans was insufficient to detect their possible changes.
The stiffening effect by progressive loading can be explained by three causes: (a) The overall longer duration when a fracture is loaded in progressive cycles. Since more time passes during the loading and unloading procedure, the fracture has longer time for compaction. Asperity deformation might be reduced and fracture consolidation is more effective. (b) Particle transportation causing partial blockage of fluid pathways.
Unfortunately, we could not analyze the effluent on any changes in fluid composition or fines migration. (c) The observed stiffening effect could be sample-dependent. All samples showed the same initial fracture permeability and similar fracture stiffness evolution suggesting a good experimental comparability and reproducibility. In fact, it cannot be ruled out that the permeability deviations during CCL and PCL scenarios can be caused by any variations of the sample by either variation in fracture geometry, asperity strength, clay content, and so on.
During the hold phases, time-dependent permeability reductions are present, but with a limited impact on the overall permeability. Reductions were observed in the first loading cycle up to 30 MPa in the PCL experiments. In all subsequent cycles, permeability remained almost constant. In contrast, the CCL experiments showed a continuous decrease in permeability in all hold phases at 60 MPa (Figures 4b, 4d and 4f). Time-dependent permeability changes during the hold phases contribute by about 14%-15% to the total permeability changes caused by the pressure changes in the PCL experiments. In the CCL experiments creep contributes by about 0.5% although the permeability change was overall larger. We believe that the time frame of 20 min is too short to define a clear exponential or power-law behavior which is typical for mechanical creep for fracture closure (e.g., Im et al., 2018).
The roughness, aperture distribution, and contact-area control the plastic and elastic deformation of asperities, as well as the fluid flow in fractures (e.g., Cook, 1992;Pyrak-Nolte & Morris, 2000;Zou et al., 2020). We did not find a change in roughness by cyclic loading experiments, similar to Yoshioka (1994). The roughness exponent was found to be similar before and after the experiment at frequencies above the length of the grain size. Most changes in surface topography are related to changes at the grain surface scale and/or to an elastic rearrangement of grains near contact-areas. However, calculating the mean aperture at various shear offsets using the post-experimental fracture surface scans (Figure 9b) and the aperture histograms ( Figure 8) revealed a mechanical imprint. It is likely that fractures stiffness is not controlled by the bulk properties of rock, but by the grain properties. The self-affine roughness exponent is not affected by a complex mechanical loading history, proving its universality. The amount of aperture reduction did not exceed 0.12 mm comparing the mean aperture before and after the experiment. While plastic deformation of asperities seems to dominate the initial loading cycle with permanent fracture closure and permeability reduction, the following loading cycles lead to reversible fracture closure. The deformation process is dependent on the distribution of asperities (Zou et al., 2020). Hence, fracture stiffness depends on the shear offset and the resulting evolution of the fracture aperture. Larger or smaller offsets than the 0.5 mm from our experiments possibly lead to a different evolution of stiffness. Especially when considering a percolation threshold that leads to large permeability reductions. The contact-area and is dominated by the elastic properties of rock (Cook, 1992;Pyrak-Nolte & Morris, 2000) after the first loading cycle when most plastic deformation was done.
To be able to apply the experimental results to the field scale, the following aspects need to be considered. Creep should be considered as an important factor for larger systems due to possible temperature changes, rock-fluid interactions, or changes in stress distribution. Permeability can only be recovered by increasing the pore pressure in the fracture or an increased shear offset.
What we cannot address in this study is the mechanics behind each released pressure step and the dependency of a variety of geometries on the stress-paths dependent permeability. These are some potential aspects that should be considered in future studies.

Limitations of the Experimental Data
In the following, we review the assumptions made for our calculations and the limitations that emerged.
The fracture stiffness was calculated from the corrected fracture closure (Equation 11) by using the circumferential extensometer data (Equation 9) according to Bandis et al. (1983). The correction by the elastic strain impacts the slope and magnitude of the fracture closure and must be treated with caution. Regardless, this correction is crucial to be able to determine the strain caused by the deformation of a fracture only. Due to the cylindrical geometry with a diametrical fracture, the measured strain is a length phenomenon depending on the size of the fracture and the matrix surrounding it. Strain is therefore not homogeneous within the sample. The measured change in aperture (Equation 10) also depends on the local position of the extensometer along the sample and is controlled by local closure magnitudes. The roughness measurements (Section 3.5.1) showed that plastic asperity deformation takes place on the grain scale, although we consider bulk measurements of fracture closure and stiffness. Consequently, the aperture measurement must be considered as an indirect measurement.
Similar to the strain distribution, the hydrostatic confining pressure applied to the sample is not distributed equally throughout the sample. Depending on the fracture topography and sample geometry, the stress acting across sample and fracture varies. In the following we estimate the stress acting across the contact-area. From the aperture distribution ( Figure 8) and the measured fracture closure, Δ E a , we calculated the evolution of the contact-area, c E R , as described in Section 3.5.2 (Figures 10a and 10b). Dividing the applied pressure by the computed contact-area, we obtained the stress acting on the fracture contacts at the respective applied effective pressure level (Figures 10c and 10d). During initial loading, the contact stress during the CCL experiment reached its peak of almost 700 MPa at about 4 MPa applied effective pressure and a decrease to about 450 MPa at 60 MPa (Figure 10c). The 450 MPa contact stress exceeds the uniaxial compressive strength of about 57 MPa for the Flechtingen sandstone as measured by Hassanzadegan et al. (2012). In all subsequent loading cycles, the contact stress was reversible without a peak. Most asperity damage is therefore done during initial loading at low applied effective pressures. Applying the same procedure to the PCL data, the contact stress increased until reaching the previous stress level (Figure 10d). The contact stress approached a limiting value of about 350-400 MPa at 60 MPa applied pressure. This suggests, that there is an universal contact stress level, in our case between 350 and 450 MPa, which controls the fracture stiffness at larger applied stresses. What this value represents, for example, the uniaxial compressive strength of a quartz grain, is not clear.
Although this shows the complex stress distribution within the sample, we assumed a homogeneous strain and stress distribution. We also consider fracture stiffness as a bulk property, albeit it is not controlled by the bulk properties of rock. Therefore, the uniaxial compressive strength of the bulk rock is not controlling the asperity strength as described by Milsch et al. (2016). These calculations and considerations reveal that we need to define clearly, what is actually measured and what is assumed in such laboratory experiments.
We further assume that the mechanical processes during cyclic loading (hysteresis) were not affected by the data acquisition rate (1 Hz). This is due to the low loading rate of 0.5 MPa/min and the measured displacement magnitudes within the mm-range. Studying rock samples with high-frequency loading rates and observations on a smaller spatial scale (nm-range) might require an adaption of the acquisition rate to properly determine all dynamic processes (Bella et al., 2021;Scalerandi et al., 2020).

Conclusions
In this study, we were able to demonstrate a novel experimental procedure to depict the fracture stiffness evolution during two different loading scenarios: constant cyclic loading (CCL) and progressive cyclic loading (PCL). Due to the high resolution of the deformation and pressure data, we were able to reveal a stress-memory effect of fracture stiffness during cyclic hydrostatic loading. Measuring the evolution of the hydraulic properties suggested that the permeability is dependent on the stress history. Overall, we suggest the following conclusion to be made from our experimental results: (a) Initial loading of a fracture leads to a linear stiffness evolution. The linear trend is non-reversible when re-loading the fracture within the same stress range. The second and all subsequent cycles show a nonlinear and almost reversible behavior. The responsible micro-mechanical deformation modes (elastic, plastic, and elasto-plastic) in each phase remain to be evaluated. (b) When exceeding the previous stress level, the stiffness evolution turns from a nonlinear to a linear behavior. This suggests a stress-memory effect in fractures similar to the "Kaiser Effect" in intact rocks. (c) The permeability of a fracture is stress-path dependent. Progressive cyclic loading (PCL) potentially leads to a stiffening of the fracture at stress levels below the previous maximum stress. Therefore, the reduction caused by effective stress changes in fractured rocks could potentially be mitigated by a cyclic, step-wise pressure function. (d) The stiffening effect might also hold for larger-scale reservoirs where a reduction in productivity can be related to a decrease in pore pressure after stimulation and during production. We therefore suggest to verify cyclic or step-wise pressure reductions in field tests.
(d) The fracture surface roughness above grain scale remains unchanged applying stress of up to 60 MPa. This supports the universality of the self-affine roughness exponent, since it is not affected by a complex mechanical loading history. Topography changes were indicated by a change in aperture distribution and a mechanical imprint, which reduces the self-propping effect at the given displacement.