Time Dependent Mechanical Crack Closure as a Potential Rapid Source of Post‐Seismic Wave Speed Recovery: Insights From Experiments in Carrara Marble

Seismological observations indicate strong variations in wave velocities around faults both co‐seismically during earthquakes, and post‐seismically. Recovery is commonly associated with a reduction in crack damage. Here, we explore the recovery associated with time‐dependent mechanical closure of cracks. We report results from laboratory experiments conducted on dry cores of Carrara marble at room temperature. We deformed cylindrical samples in the semi‐brittle regime to induce crack damage before subjecting them to hydrostatic and triaxial stress conditions for extended periods of time while recording dilatancy and wave speeds repeatedly. We report wave speed increases of up to 40% of the damage‐induced wave speed drop in samples subject to hydrostatic loading. Moreover, we report the occurrence of significant wave speed increases contemporaneously with time‐dependent creep in triaxially loaded samples. Wave speed recovery during creep is only observed below a threshold creep strain rate, a result we interpret as a transition from brittle to plastic creep with decreasing strain rate. We interpret the wave speed increase in terms of reduced crack density and increased contact area within the crack array, and show that around 40% of the total crack surface has to be closed to justify the observed wave speed recoveries. We propose that mechanical crack closure is driven by the viscous relaxation of the bulk rock under the influence of locked‐in stresses at low confining pressure, and of the external stresses at higher confining pressure. Our study shows that mechanical crack closure is a significant source of time‐dependent wave speed recovery.


Introduction
The slip behavior of crustal faults is largely controlled by the elastic properties of the surrounding country rock. However, these properties are not constant over time and may vary considerably over the seismic cycle. For instance, it has been observed that local elastic wave speeds around faults (a direct proxy of rock elastic properties) drop significantly during earthquakes (see Table 1 and references therein), which reveals the existence of an increase in compliance of the fault host rock due to earthquake ruptures. Since cracks greatly impact rocks elastic properties (and therefore the propagation of elastic waves) (e.g., Birch, 1960;Birch, 1961;Guéguen & Sarout, 2011;Lockner et al., 1977;O'Connell & Budiansky, 1974;Sayers & Kachanov, 1995), this result is commonly associated with the nucleation and/or reactivation of crack damage in the fault damage zone due to dynamic changes of stress or strain (e.g., Brenguier et al., 2008;Olivier et al., 2015;Rubinstein & Beroza, 2004;Taira et al., 2015).
Similarly, a growing number of seismological studies have revealed that, in the post-seismic phase of the earthquake cycle, the co-seismic drop in elastic wave speeds is partly or totally recovered over time (see Table 1 and references therein). Naturally, this observation has been interpreted as a diminution in crack damage intensity over time, a phenomenon often referred to as fault or damage "healing".
In the following, we wish to make a distinction between the overall "healing" phenomenon and the underlying driving processes, and we will therefore restrict our use of the term healing to chemical processes (e.g., diffusion or reactions) that lead to permanent removal of discontinuities from a medium. Mechanical processes involving crack aperture reduction but no removal of the discontinuity will be simply referred to as closure processes. To further clarify our terminology, the general phenomenon of post-deformation increase in wave speed or elastic moduli will be called recovery.
Abstract Seismological observations indicate strong variations in wave velocities around faults both co-seismically during earthquakes, and post-seismically. Recovery is commonly associated with a reduction in crack damage. Here, we explore the recovery associated with time-dependent mechanical closure of cracks. We report results from laboratory experiments conducted on dry cores of Carrara marble at room temperature. We deformed cylindrical samples in the semi-brittle regime to induce crack damage before subjecting them to hydrostatic and triaxial stress conditions for extended periods of time while recording dilatancy and wave speeds repeatedly. We report wave speed increases of up to 40% of the damage-induced wave speed drop in samples subject to hydrostatic loading. Moreover, we report the occurrence of significant wave speed increases contemporaneously with time-dependent creep in triaxially loaded samples. Wave speed recovery during creep is only observed below a threshold creep strain rate, a result we interpret as a transition from brittle to plastic creep with decreasing strain rate. We interpret the wave speed increase in terms of reduced crack density and increased contact area within the crack array, and show that around 40% of the total crack surface has to be closed to justify the observed wave speed recoveries. We propose that mechanical crack closure is driven by the viscous relaxation of the bulk rock under the influence of locked-in stresses at low confining pressure, and of the external stresses at higher confining pressure. Our study shows that mechanical crack closure is a significant source of time-dependent wave speed recovery.

MEYER ET AL.
While the phenomena listed above are all chemical in nature, mechanical processes have also been reported to dramatically and rapidly impact wave speeds and crack anisotropy in rocks. Microcrack opening, and therefore effective elastic properties of rocks, is known to be very sensitive to applied stresses: Reversible crack opening and closure has been evidenced during cyclic loading in granite (e.g., Holcomb, 1981;Passelègue et al., 2018). Therefore, stress variations have a direct impact on elastic properties, which could be source of crack closure as observed in seismological records. However, even under constant stress conditions, cracks close in a time-dependent manner. In dry granite, Scholz and Kranz (1974) reported gradual recovery of dilatancy after uniaxial compression cycles, which they interpreted as time-dependent reverse sliding on "wing-cracks". Similarly, Brantut (2015) showed time-dependent recovery of elastic wave speeds in deformed limestone, also consistent with crack closure due to reverse sliding on internal defects. Elastic wave speed recovery is also commonly observed during stress-relaxation tests, for instance in deformed carbonates (Kaproth & Marone, 2014;Schubnel, et al., 2005). Despite its significance, time-dependent mechanical crack closure remains poorly understood, and the underlying microphysical processes are enigmatic. Here, we propose to palliate this lack by means of experiments on dry cores of Carrara marble in a conventional oil medium triaxial apparatus. We deformed cylindrical samples in the semi-brittle regime to induce diffuse crack damage and subsequently exposed the same samples to constant stress conditions for extended periods of time while repeatedly recording elastic wave speeds. Additionally, we conducted a detailed examination of the microstructures produced during the experiments using both optical and electron microscopy. We report the occurrence of considerable wave speed recoveries in samples left under hydrostatic stress, at room temperature and in dry conditions. Similarly, we report sizable wave speed recovery in samples left under constant triaxial stress conditions and subjected to time-dependent creep. We detail the existence of two distinct creep regimes, a first regime where strain rate is high and wave speeds decrease and a second regime where strain rate is lower and wave speeds recover. We interpret this observation as marking a transition from brittle to plastic creep with decreasing strain rate. Furthermore, we use a model of islands of contact to interpret wave speed recovery and find that a reasonable increase in contact area is sufficient to explain the amplitude of the observed recovery. We hypothesize that viscous relaxation of the grains driven either by internal "locked-in" stresses or by the externally applied stresses controls the closure of cracks. Finally, we conclude that crack mechanical recovery on its own can be a rapid and potent damage recovery mechanism, potentially representing a large part of the wave speed increase around faults. In conjunction with chemical healing, it could lead to the disappearance of cracks from the damage zone over a few days only.

Material and Apparatus
Carrara Marble is a calcite aggregate from Northern Italy, monominerallic in nature and with a very low starting porosity (0.1%, Rutter, 1972). Its fabric is essentially isotropic, with a relatively uniform grain size of around 200 μm. We selected this rock type because its rheology has been extensively studied previously, so that its deformation properties and brittle-ductile transition at room temperature and confinements up to 300 MPa are already known (Fredrich et al., 1989). We cored 40 mm in diameter cylindrical samples from a single quarried block. Samples were then cut to 100 mm in length and the end surfaces ground flat and parallel using a surface grinder. Each prepared sample was then fitted with a pair of radial (horizontal) and a pair of axial (vertical) strain gauges, glued onto the rock surface and used to record axial and radial strain from which volumetric strain in the sample is computed. Samples were finally oven dried at 60°C for at least 48 h prior to testing.
All experiments were performed using an oil medium, conventional triaxial deformation apparatus located in the Rock and Ice Physics Laboratory at University College London (Eccles et al., 2005). Experimental samples are placed in a steel pressure vessel filled with silicon oil. Confining pressure is applied and maintained by an external pump which can reach a maximum pressure of 400 MPa with an accuracy of 0.4 MPa. The load is applied via a 150 tonnes servo-controlled hydraulic actuator. The load is transmitted to the sample with a self-compensated piston fitted with two hemispherical seats to ensure an even vertical loading of the sample. The load is recorded by an external load cell and the total shortening (sample and loading column) is recorded by two external linear variable differential transformers (LVDT). Shortening from the LVDT is corrected from the loading column deformation and divided by the initial length of the sample to access vertical strain.
Samples were isolated from the oil pressurizing medium by an engineered Viton rubber jacket fitted with ports to house 14 P-wave and 2 S r -wave piezoelectric transducers. The geometry of the transducer array allows for the recording of the P wave speed along four directions with respect to the vertical compression axis (28°, 39°, 58°, and 90°) and a single S r wave speed along the horizontal (radial) direction (see, for instance, Brantut et al., 2014, for jacket design details). Transducer signals were transmitted to the outside of the pressure vessel through ceramic lead-throughs and connected to 40 dB pre-amplifiers. The array of amplifiers was in turn connected to a pulser and 50 MHz digital oscilloscopes. Active wave speed surveys were performed by sending a 200 V and 1 MHz pulse to a selected source transducer while recording received waveforms at all the other transducers. Six recorded signals were stacked to increase the signal-to-noise ratio. A full wave speed survey consisted of pulsing each transducer as a source sequentially. The change 10.1029/2020JB021301 4 of 29 in average wave speed between each pair of transducers was computed using waveform cross correlation between a set of manually picked arrival time in a "master" survey, and those of the successive surveys (e.g., Brantut et al., 2014).
Several selected samples were retrieved after the experiments and thin sections were produced in planes both parallel and perpendicular to the loading axis. The sections were produced in two batches at distinct moments in time and by two different manufacturers, limiting our capacity to compare some of them. The thin sections were investigated using both, optical and scanning electron microscopy using the backscattered electron (BSE) mode. BSE imaging and focused ion beam (FIB) cross-sectioning was carried out in a Thermo Fisher Scientific (former FEI) Helios NanoLab G3 machine located in the electron microscopy facility at Utrecht University.

Wave Speed Measurement Accuracy
In our wave speed acquisition process, there exist three sources of error. First, errors on the absolute wave speed values originate from the picking process. Arrival time picks can be made with an accuracy of around 3 μs, so that we estimate the absolute error on the wave speed data to be about ±200 m/s. This error is systematic, in other words, it is constant and offsets all wave speed data in a given data set and does not impact the relative changes in wave speeds.
The cross correlation method employed here is very precise and the error arising from its use, the second source of error in our data, depends on the sampling frequency of the waveforms. Here, we resample all data on a vector with a frequency twice larger than that of the raw data before conducting the cross correlation, i.e., 100 MHz. This results in a resolution on arrival times of the order of 0.01 μs which translates to a resolution of about ±1 m/s on variations in wave speeds.
A third source of error is the change in sample dimensions due to deformation; any given relative shortening or extension of the distance between sensors would translate into the same decrease or increase in travel time, even if the wave speeds remain unchanged. In this study, we do not correct wave speeds for the wave travel path distortion induced by sample deformation; such systematic corrections are complicated by the tendency of samples to become barreled. However, we estimate that these corrections remain small compared to the observed changes, as discussed below. Since we focus on the radial wave speeds, the wave speed data presented here may be overestimated by an amount equal to the sample radial strain. As shown later in this work, the maximum recorded axial strain was about 2.5%, and the maximum radial strain was around 3% during the deformation stage of the experiments. The observed relative change in wave speed is far greater than these recorded relative changes, so the correction here would be negligible. Similarly, the maximum radial strains recorded during constant stress and incremental stress recovery experiments were 3% and 2.5% respectively. These deformations are positive, i.e., the wave speed travel path are becoming longer leading to an underestimation of our recorded wave speed recovery of about 3%.

Recovery Experiment Protocol
Within the overall experimental program, several sets of experiments were conducted with varying protocols: Hydrostatic recovery experiments, constant differential stress recovery experiments and incremental differential stress recovery experiments. These are all listed (see Table 2). All experiments were divided into two stages; first, a deformation stage and then an ensuing recovery stage. The deformation stage was identical for each experiment, with the difference being the experimental protocol applied to the recovery stage only. Between deformation and recovery, the samples were never removed from the pressure vessel. Moreover, both strain and wave speed data were recorded throughout both stages of all experiments.

Deformation Stage
In order to induce microcrack damage in the samples, we deformed them at room temperature, under a confining pressure P c = 40 MPa and at a constant strain rate    5 10 E  s −1 until approximately 2.4% axial strain was accumulated in the sample (ranging from 2.3% to 2.5% depending on the sample). The maximum differential stress Q max reached during deformation was recorded for all samples. After the target strain was 10.1029/2020JB021301 5 of 29 reached, the differential stress was fully unloaded using the same constant strain rate, before proceeding to the recovery stage. The inelastic irrecoverable strain after unloading ranged from 1.7% to 2.0%.

Hydrostatic Recovery
For our hydrostatic recovery experiments, we completely removed the differential load and then immediately (less than 10 s) increased or decreased the confining pressure to a pre-determined target value (Table 2) at a rapid rate (≈0.5 MPa/s). Samples were then left to recover under constant hydrostatic pressure for extended periods of time while strain and wave speeds were monitored.

Constant Differential Stress Recovery
For our constant differential stress recovery experiments, following complete removal of the differential load, we immediately (less than 10 s) reloaded the samples to a pre-determined target recovery differential stress, Q, at a constant strain rate of    5 10 E  s −1 and while keeping the confinement constant at P c = 40 MPa. Q was subsequently held constant for the remainder of the recovery period.

Incremental Differential Stress Recovery
Incremental differential stress recovery experiments were similar to the constant differential stress recovery experiments, we loaded and unloaded the sample at a constant strain rate    5 10 E  s −1 , but here, Q was increased step-wise. Each differential stress step was maintained constant for 24 h before proceeding to the next step. At the end of each stress step, the differential load was decreased to 10% of the sample Q max value, before being immediately (less than 10 s) re-increased to its next increment value. Samples were not fully unloaded in order to ensure that the load did not fall below the inherent machine friction so that there would be no hysteresis in the loading cycle.

Complementary Hydrostatic Wave-Speed Experiments
In order to record the pressure dependency of elastic wave speeds in Carrara marble we conducted additional pressure survey experiments on both undeformed and previously deformed samples. In a first experiment, we pressurized an intact sample step-wise up to P c = 150 MPa in 20 MPa pressure steps at a rate of approximately 0.5 MPa/s. At each step, the pressure was held constant for 10 min (to allow for confinement The notation is as follows: Pc, confining pressure; Qmax, maximum differential stress reached during the deformation stage of the experiment; Q/Qmax -recovery differential stress in percents of the maximum differential stress reached during the deformation stage of the experiment. The creep shortenings presented here correspond to the shortening undergone by the sample during the recovery stage of the experiment at a given differential stress (in the case of the incremental stress experiment). Type refers to the genre of recovery performed on the sample (see Sections 2.3.2, 2.3.3, and 2.3.4 for details). stabilization), and wave speeds recorded every minute during that time. After having reached 150 MPa, the pressure was then decreased step-wise in an identically reversed way down to P c = 10 MPa.
In a second experiment, we subjected a sample to the deformation stage (i.e., vertical deformation at P c = 40 MPa up to 2.5% axial strain, followed by differential unloading), and then increased the confining pressure to 150 MPa in 20 MPa steps at a rate of approximately 0.5 MPa/s. The confinement was then again lowered step-wise with 20 MPa pressure steps. Similarly to the experiment above, each confining pressure step was held for 10 min and wave speeds recorded every minute.
Finally, we selected the sample that had been left to recover at P c = 150 MPa (CMh12), to undergo an additional experimental step after recovery. In this case, the confining pressure was decreased step-wise from the recovery confining pressure down to atmospheric pressure by 10 MPa steps at a rate of approximately 0.5 MPa/s. Once again, the sample was held under constant confining pressure for 10 min at each step and the wave speeds recorded every minute.

Pressure Surveys
In intact Carrara Marble, wave speeds increase noticeably with increasing confinement (Figure 1a). The radial P wave speed increases by 20% of its initial value from around 5,000 m/s at P c = 10 MPa to 6,000 m/s at P c = 150 MPa. Similarly, V Sr increases by 20% from around 3,000 to 3,600 m/s over the same pressure range. At P c = 150 MPa, wave velocities are consistent (within experimental errors) with Voight-Reuss-Hill averages in pure calcite (using bulk and shear moduli of 76.1 and 32.8 GPa, respectively (Chen et al., 2001)), which suggests that the overwhelming majority of voids are closed.
In the deformed sample (Figure 1b), the radial wave speeds also increase with increasing confining pressure, but start from a much lower value of V P = 3,500 m/s at P c = 10 MPa and increase by 31%, up to around 5,100 m/s at P c = 150 MPa. A similar trend is observed for the shear wave speed which increases by 14% from 2,800 m/s at P c = 40 MPa to 3,200 m/s at P c = 150 MPa. In the deformed sample that underwent 9 days of recovery at 150 MPa confining pressure, the radial P wave speed increases from V P = 3,600 m/s at P c = 10 MPa up to around 5,400 m/s at P c = 150 MPa which represent an increase of 50%. The radial S wave speed increases by 40% from V Sr = 2,400 m/s at P c = 10 MPa up to around 3,350 m/s at P c = 150 MPa.
In addition to wave speed measurements, we use the method developed by Walsh (1965) to infer the crack porosity in our samples. This method relies on the fact that, at sufficiently high confinement, all cracks in the rock have closed and the overall volumetric strain response to a change in confining pressure follows the linear elastic compressibility of the bulk rock. By graphically extrapolating this bulk response down to zero confinement, it is therefore possible to estimate the total crack volume in the rock at room pressure. Such estimates provide lower bounds of the total porosity, since equant voids cannot easily be closed even at the highest tested pressure in our experiments. The authors demonstrated that this method is very accurate and could determine crack porosity in granite with an uncertainty of around 0.1%. We applied this method to our pressure survey data by visually fitting a line to the straight part of the curves in Figure 1 and find an initial crack porosity in intact Carrara marble of 0.04% ( Figure 1d). In the deformed sample, the crack porosity at room pressure is found to be 0.33% ( Figure 1e) and in the sample that has undergone 9 days of recovery at P c = 150 MPa, crack porosity it is found to be 0.36% (Figure 1f).

Deformation
Typical data for the deformation stage of the experiments are presented in Figure 2. During deformation, stress initially increases with axial strain quasi-linearly ( Figure 2a). Eventually, at around 0.3% axial strain, the stress rolls over with a large amount of strain hardening until it reaches a maximum (Q max ) of 200 MPa at the end of the loading phase. During unloading, the stress drops rapidly with decreasing strain until it reaches zero at an axial strain of 1.8%. Volumetric strain calculated from strain gauges data shows a slight initial increase of 0.05% (denoting compaction) over the first 0.5% of axial shortening before plummeting to −3.5% (denoting a strong dilation) by 2.5% axial strain ( Figure 2b). During the unloading phase, volumetric strain increases up to a final value of −2.8%. Simultaneously, the wave speeds along all directions of propagation initially increase by a few tens of m/s upon initial deformation, reaching a maximum at around 0.3% axial strain. They subsequently decrease significantly with further deformation ( Figure 2c). The amplitude of the P wave speed reduction with increasing axial strain is greatest for the radial propagation direction (90°) with a decrease of 39% of the maximum wave speed, and least for the sub-axial propagation direction (28°) which decreases by 25%. The S r wave speed decreases in the same fashion by 42% of its maximum value. Moreover, upon unloading, all the wave speeds in the sample increase, albeit in different ways. At 58° and 90°, the rate of increase in wave speeds is greater than the rate of decrease in strain; whereas at 28° and 39° the rate of increase is lower than the rate of decrease in strain. It is interesting to note here that the wave speed variability between all samples prior to deformation was below the absolute accuracy of the arrival time picking method. The deformation stage of the experiment did not lead to macroscopic sign of strain localization and but rather slight barreling ( Figure 3). s −1 . In A. is differential stress against axial strain, in B. is volumetric strain against axial strain and in C. is elastic wave speeds against axial strain for different directions of recording. The angles in brackets are expressed in degrees with respect to the direction of compression. The arrow in A. is a visual cue of the Q max picking method.

Hydrostatic Recovery
During the hold time of the hydrostatic recovery experiments, both P and S r wave speeds increase along all propagation paths ( Figure 4a).
Samples held under hydrostatic stress conditions also show some variations in volumetric strain over time (Figure 4b). At confining pressures below 40 MPa, samples contract by 0.05% of their post deformation volume, at P c = 40 MPa, volumetric strain after 5 days is 0.08% and at confining pressures above 40 MPa, sample contract by around 0.15% of their post deformation volume over the same time.
The rate of wave speed increase is not constant, but typically follows a logarithmic relationship with time ( Figure 4c). In terms of order of magnitude, at P c = 40 MPa, the P and S wave speeds increase by around 100-150 m/s along all propagation angles during the first 24 h of recovery. The amount of recovery increases with increasing confining pressure above 40 MPa, up to around +300 m/s in the first 24 h at P c = 150 MPa. In the sample held at P c = 3 MPa, the increase in P wave speed is significant and comparable to that observed at P c = 20 and 40 MPa, of around +100 m/s in 24 h.  ) computed with the radial P wave speed (V P (90°)) for all samples subjected to hydrostatic recovery (D.). In B, C, and D the numbers adjacent to the curves represent the recovery confining pressure (in MPa).
To better quantify the extent of the wave speed recovery, we normalize the wave speed increases during the holding time by the amplitude of the wave speed drop induced during the deformation stage of the experiment. Thus, we define the recovery factor R as where V is the evolving wave speed, V 0 is the wave speed at the beginning of the recovery period (i.e., after the hydrostatic confining pressure has been set to its recovery value) and V ref is the measured reference wave speed in the intact sample at the confining pressure corresponding to the recovery pressure (Shown in Figure 1a). With this formulation, the recovery is expressed as the recovered wave speed over the total wave speed reduction induced by the deformation episode. For example, a recovery factor R = 50% would indicate that the wave speed in the sample has recovered 50% of the deformation-induced wave speed reduction. In this study, we will calculate R with the P wave speed along the radial direction of recording as it is the most sensitive to the deformation (see Figure 2). The problem of anisotropy and its evolution during recovery will be tackled in a later part of the article.
R increases in a log-linear fashion with time at all tested confinements (Figure 4d). At P c = 3 MPa, the sample recovers 5.5% of the damage-induced wave speed reduction after five days, whereas at P c = 150 MPa, R is equal to 35% after five days.

Constant Differential Stress Recovery
For convenience, in the following, we will express the magnitude of the differential stress acting on any sample during a recovery hold time as the percentage of the sample maximum stress sustained during the deformation stage, Q max . The value of Q max is slightly different for different samples, due to natural sample variability (see Table 2 for all Q max values in MPa and Figure 2a for a visual representation of the picking process).
During the recovery stage of the constant differential stress experiments, all samples show time-dependent creep, with the rate being sensitive to the imposed differential stress Q (Figure 5a). At Q/Q max = 30% and 60%, the axial strain of the samples is 0.01% after 24 h, at Q/Q max = 95%, the axial strain is 0.53% after 24 h; and at Q/Q max = 97.5%, the sample shortens by 0.95% over 24 h. It is important to note here that the vertical, step-like kinks in the axial strain data originate from noise on the actuator and have no link with the time-dependent creep.
Simultaneously, the samples left at Q/Q max = 30% and 60% both contract by 0.025% of their post-deformation volume over 5 days ( Figure 5b). Unfortunately, the strain gauges on the samples left at Q/ Q max = 95% and 97.5% were defective and we therefore do not have volumetric strain data for these samples.
During the recovery stage, wave speeds increase logarithmically with time in all samples held at constant stress Q/Q max up to 97.5% (Figure 5c), with a rate of about 100 m/s over 24 h, which is similar to that in the sample tested under purely hydrostatic condition at P c = 40 MPa. At the two highest differential stresses Q/ Q max = 95% and Q/Q max = 97.5%, over the initial 15 min, the radial P wave speed remains stable or slightly decreases (by around 10 m/s) prior to the recovery, whereas at Q/Q max = 30% and 60%, the wave speed increases monotonically.
The evolution of the recovery factor R during the constant differential stress experiments is similar to that of the wave speeds (Figure 5d). At Q/Q max = 30% and 60%, R increases in a monotonic fashion and reaches 11% and 8% after five days respectively. At Q/Q max = 95% and 97.5%, R initially decreases by 0.1% and 0.4% over the first 15 min of recovery before increasing and finally reaching 8% after 3 days at both differential stresses.

Incremental Differential Stress Recovery
During the incremental differential stress recovery test, the sample showed time dependent axial creep at all differential stresses (Figure 6a), and the axial shortening over time is stress dependent. Over a 24 h period, the sample shortened by 0.07% at Q/Q max = 75%, and by 0.53% at Q/Q max = 99%. At Q/Q max = 102%, creep was unstable and led to sample-scale brittle failure.
Volumetric strain data were somewhat noisy ( Figure 6b). Nevertheless, some trends can be observed. At differential stresses below Q/Q max = 92%, the sample initially showed dilation, and after a few hours volumetric strain stabilized and became slightly compactant. At Q/Q max = 97% and 99%, only dilation was observed.
Contemporaneously, the radial P wave speed in the sample increased at all stress steps, except at Q/ Q max = 102% where it decreases dramatically until failure (Figures 6c and 6d). While the wave speed increase at Q/Q max = 75% was monotonic and of around 130 m/s (R ≈ 6%) in 24 h, at higher differential stress, the wave speed initially remained stable or slightly decreased over the first 15 min of hold time. This initial decrease tended to be more pronounced and of longer duration with increasing differential stress.

Microstructure
In addition to the thin sections prepared from samples after recovery tests, thin sections were prepared from one sample that had undergone the first stage of the experiment only (deformation stage), and from an intact core of Carrara marble. In selected samples, two sections were cut: one in a plane perpendicular to the Intact Carrara marble is characterized by a homogeneous distribution of calcite grains of sizes ranging between 70-220 μm (Figure 7). The grains have well defined boundaries and are generally equiaxed. Twins are present and can be observed in some grains. Intragranular microcracking is mostly absent in the intact microstructure.
In the deformed microstructures, cracking and twinning are ubiquitous. In the radial thin section, the crack array consists almost entirely of intragranular cracks with no particular orientation (Figure 8a). The intragranular cracks seem to mostly originate from the grain boundaries where their aperture is greatest (up to few tens of μm), and generally span the entire length of the grains. While some grains are pervasively cracked, other present little to no sign of cracking. Moreover, the density of twins is noticeably higher than in the intact section and some grains show dense collections of twins. The microstructure in the axial section is similar to that of the radial section with a dense array of wide intragranular cracks and intense twinning ( Figure 8b). Nevertheless, one important difference is that the cracks in the axial section are preferentially oriented parallel to the direction of compression. After 5 days of recovery at P c = 20 MPa (Figure 9a), the microstructure is qualitatively similar to that of the deformed sample (Figure 8), with the presence of a dense array of wide intragranular cracks and pervasive twinning.
After 9 days of recovery at P c = 150 MPa (Figure 9b), intragranular cracks are still present, but apparently with a much lower density than in samples held at lower pressure. Careful optical examination reveals numerous intragranular cracks with extremely narrow apertures. Twinning is still widespread in the thin section and there appears to be no difference between the number of twins in this microstructure and that of the deformed sample ( Figure 8).
After 3 days of recovery at P c = 40 MPa and with a differential stress Q = 200 MPa (Q/Q max = 97.5%), wide microcracks are widespread throughout the microstructure (Figure 9c). The extent and aspect of the crack array is comparable to that of the sample left recover under hydrostatic conditions at P c = 20 MPa. It is important to note here that, the thin sections presented in Figures 9c and 9d are thicker than those shown above. This is due to the manufacturer of the sections being different. For this reason, the microstructure appears darker, cloudier, and overall more cracked because the three-dimensional complexity of the grain boundaries is visible within the thickness of these sections. As a result, twins are indiscernible, rendering a quantitative assessment of their number impossible.
After less than 1 h at P c = 150 MPa (wave speed survey experiment, see Figure 1b and 1e), the microstructure appears extensively cracked (Figure 9d). The slightly larger thickness of the section makes direct observations of crack opening difficult; by comparison with the thin section in Figure 9c, which was of same thickness, we conclude that microcrack opening tends to be reduced in the sample that underwent pressurization to 150 MPa.
To resolve for the finer characteristics of the crack network in the microstructures, we used BSE imaging in a scanning electron microscope ( Figure 10). Due to the tendency of cracks to propagate axially (i.e., parallel to the direction of maximum compression) during triaxial deformation, we focused our attention on radial sections which intersect most of the cracks perpendicularly, consequently giving a better image of crack aperture throughout the microstructure. The BSE image of intact Carrara marble reveals the existence of some open porosity at the boundaries between grains (Figure 10a). Additionally, the intact BSE microstructure further demonstrates the initial scarcity of cracks in Carrara marble.
The deformed BSE microstructure illustrates the extent of cracking during the deformation of Carrara marble (Figure 10b). Intragranular cracks appear relatively straight with little rugosity on their surfaces. Grain boundaries show increased width compared to that in the intact thin section. Moreover, the BSE image exposes the existence of transgranular shear zones, invisible in the optical microstructures, which contain smaller fractured grains (denoted CA in Figure 10b.). This might indicate some level of shear between grains accompanied by grain comminution during the deformation of the sample.
After 5 days of recovery at P c = 20 MPa, the sample shows an extensive array of cracks but with some that have visibly closed (denoted CC in Figure 10c). The closed cracks seem to have a reduced length and aperture. Where cracks have closed, a faint trace can be observed on the surface of the thin section (see cracks indicated CC in Figure 10c). The presence of these traces at the tips of certain partially closed cracks might be an indication that crack closure occurred from the tip inwards (see for instance the crack indicated by CC in left hand side of Figure 10c). Crack closure is more pronounced in certain grains, with some showing pervasive cracking while others are almost back to their undeformed state. Crushed material can still be observed at grain boundaries.
After 9 days of recovery at P c = 150 MPa, the sample shows a strikingly reduced crack array compared to the deformed sample (Figure 10d). The cracks in this sample have almost all closed with only the larger ones showing noticeable aperture. Traces can be observed in grains throughout the microstructure (see CC in Figure 10d). The grain boundaries are seemingly back to their intact state and the crushed grains seem to have become welded to their larger counterparts (see WCA in Figure 10d). ). Note: The sections shown in C. and D. have been made at a different time than those shown in A. and B. and they are generally much thicker, hence they appear darker and more damaged (see text). The arrows and acronyms indicate the following features: C, crack; GB, grain boundary; and T, twin. Bottom left is a visual cue of the orientation of the thin sections with respect to the direction of compression.
After 3 days of recovery at P c = 40 MPa and with a differential stress Q = 200 MPa (Q/Q max = 97.5%), the microstructure remains extensively cracked (Figure 10e). Few of these cracks appear to have closed at their tips and wide transgranular shear zones with reduced grain sizes can be observed. The overall aspect of the fabric is similar to that of the sample that has undergone recovery at P c = 20 MPa.
After less than 1 h of recovery at P c = 150 MPa, the microstructure shows a vastly reduced crack array from the deformed state (Figure 10f). Reduced crack apertures and traces are present in all grains. Grain boundaries are seemingly back to their intact state and there is no transgranular shear zone to be seen. The overall aspect of the fabric is very similar to that of the sample that has undergone recovery at P c = 150 MPa. Additionally, we excavated a cross-section for BSE imaging using focused ion beam-scanning electron microscopy (FIB-SEM) based milling in the sample that experienced recovery at P c = 150 MPa for 9 days. We targeted the tip of a closed crack in order to determine the nature of the faint traces observed in the recovered samples (Figure 11a). It appears that the traces are cracks with extremely reduced aperture, of the order of few tens of nm. Here, we can observe that cracks have some nanoscopic scale rugosity that was not previously visible. The crack shows shear jog-like features where the aperture is greatest, and, in between, we can observe all but closed contacts between the crack surfaces. The crack aperture decreases toward its tip until it totally vanishes in the grain.
Finally, it is important to note that none of the recovered samples show any sign of pressure solution and precipitation (necking at contact asperities, e.g., Beeler & Hickman, 2015), diffusion healing (array of fluid inclusion at crack tips, e.g., Renard et al., 2009) or sealing (veins, e.g., Lee & Morse, 1999).

Wave Speed Anisotropy
The wave speed data reflect the evolution of the elastic moduli of the material during deformation and recovery. Consequently, it is possible to access the elastic anisotropy and observe its evolution. This is of particular interest as the inherently anisotropic triaxial stress conditions imposed during some of our experiments induce a complex evolution of the sample anisotropy during the recovery hold times. To quantify anisotropy, we compute the Thomsen parameters (Thomsen, 1986(Thomsen, , 1995. The three Thomsen parameters (γ, δ, ϵ) describe the wave speed anisotropy (i.e., elastic anisotropy) in a transversely isotropic (TI) medium. Given that we had access to a more complete set of P wave speeds in our experimental set-up, we focus here on ϵ Thom which is linked to the fractional difference between the P wave speed in the radial (horizontal) and axial (vertical) directions. This adimensional parameter represents the strength of the anisotropy in the sample, and the further it diverges from zero and the stronger the anisotropy of the medium is. This approach involves an assumption of weak anisotropy so that ϵ Thom ≪ 1. ϵ Thom is given in terms of the elastic stiffness tensor C ij (in the Voigt notation) of the rock as: We use the dependency of the phase velocity on the incidence angle with respect to the axis of symmetry of the anisotropy of a TI medium θ (Daley & Hron, 1977) to directly access the elastic moduli of the medium such that:   (2 )] , and d is the density of the rock matrix which is assumed to be constant. In our experiments, we most likely record group velocity rather than phase velocity (e.g., Nadri et al., 2012). Here, however, we have not accounted this difference and use our wave speed data directly, because the difference is expected to be small if ϵ Thom ≪ 1. Based on our four measurements of P wave speed and one measurement of the radial S wave speed, we solve Equations 3 and 4 for the C ij s by least-square minimization. We can then use Equation 2 to compute ϵ Thom for all our experiments.
During the deformation stage of our experiments, P wave speed anisotropy initially increases slightly up to 0.3% axial strain before showing a marked decrease down to −0.25 at the end of the deformation stage ( Figure 12a). This deviation of ϵ Thom indicates a strong deformation-induced increase in P wave speed anisotropy. The small initial offset in ϵ Thom seen here is likely to be an artifact of the absolute error in picking arrival times rather than any inherent anisotropy in the microstructure (Figure 10a). During hydrostatic recovery, ϵ Thom increases from its negative post deformation value toward less negative values at all tested confinements, reflecting a reduction in anisotropy (Figure 12b). At P c = 3 and 20 MPa, the anisotropy reduction after five days is around 0.009 whereas at P c = 40 MPa and higher, it lies between 0.015 and 0.02. Interestingly, in hydrostatic conditions and at confinements greater than P c = 40 MPa, there appears to be no correlation between confinement and the amplitude of anisotropy reduction.
When subjected to recovery under constant triaxial stress conditions for extended periods of time, the samples show a clear reduction in anisotropy (Figure 12c). Additionally, the anisotropy shows an initial plateau, or delay, preceding the recovery. The differential stress amplitude does not appear to have a major impact on the amplitude of the anisotropy reduction, which always lies between 0.01 and 0.015 after five days, for all stresses tested. It is interesting to note that the anisotropy reduction is only smaller than for the hydrostatic cases (Figure 12b). In the case of the incremental differential stress experiment, the data are somewhat similar to those of the constant differential stress experiments. At all tested stresses, with the exception of Q/Q max = 75%, the anisotropy reduction shows an initial period of around 2 h where it remains constant. At Q/Q max = 75%, the anisotropy reduction is monotonic and log-linear with time. The amplitude of the anisotropy reduction after 24 h is similar at all Q (between 0.006 and 0.01).

Crack Density
The low starting porosity of our samples, the observations of the post-mortem microstructures and the calculated anisotropy data all suggest that the observed changes in wave speeds are due to changes in the microcrack density within the samples. To access and quantify these changes, we have used our wave velocity data to estimate crack density evolution during deformation and subsequent recovery. We use the model of Sayers and Kachanov (1995) later expanded on by Schubnel and Guéguen (2003) and Schubnel, Benson, et al. (2006) in order to describe the elastic anisotropy that is added to a medium by an array of non-interacting, open, penny-shaped cracks. This approach describes the state of cracking in a rock using a dimensionless crack density ρ and is well suited for a low porosity rock such as Carrara marble (Schubnel, Walker, et al., 2006). For an isotropic crack distribution, ρ is defined as: where N is the number of cracks, l is the mean crack radius and V is the rock volume being considered. In the more complex case of transverse-isotropy, the crack density can be represented by two densities ρ a and ρ r which represent the axial (parallel to the direction of compression) and radial (perpendicular to the direction of compression) crack populations, respectively, and for which N and l will also likely differ. The crack densities are computed using the wave speed data following the procedure described in Brantut, Schubnel, and Guéguen (2011).
During deformation, axial crack density shows an initial plateau followed by a slight drop around 0.5% axial strain and a dramatic increase from 0.08 to 0.65 over the remainder of the deformation (Figure 13a). The radial crack density plateaus at 0 up to around 1% axial strain before increasing to 0.06 over the remainder of the deformation. This result corroborates the fact that, under these pressure and temperature conditions, Carrara marble is in the semi-brittle regime (Fredrich et al., 1989) and most of the damage induced during the deformation stage of our experiment is comprised of sub-axial microcracks. The variations of the axial crack density being tenfold that of the radial density, the former will be the focus of our analysis.
During the hydrostatic recovery hold times, the wave speed increase is interpreted as a log-linear decrease in axial crack density (Figure 13b). While at P c = 3 MPa, the axial crack density decreases by 0.33 over five days, it decreases by only 0.037 only over the same period at P c = 150 MPa. Therefore, it appears that the crack density reduction is faster and of greater amplitude at lower confinement.
The axial crack density decrease under all tested constant differential stress conditions (Figure 13c). The general evolution mirrors that of the wave speeds under the same conditions with an initial plateau at Q/ Q max = 30% and 60% or a slight increase at Q/Q max = 95% and 97.5%, preceding a decrease with time. Similarly, axial crack density decreased under all tested differential stress conditions during the incremental stress experiment. At Q/Q max = 75%, the decrease is monotonic and almost log-linear, whereas at higher differential stress, axial crack density shows an initial stress-sensitive increase.

Island Factor
Until now, we have interpreted wave speed recovery as a reduction in crack density. In the formulation for the crack density where microcracks are assumed to be of a unique size and circular shape, such a reduction implies a decrease either in number of cracks N or in crack average radius l (Equation 6). In practice, our microstructural observations indicate that microcracks have more complex shapes than assumed in the model, and that their apertures are variable.
We observed in the microstructures that cracks are inherently rough and undergo some degree of aperture reduction during the recovery process ( Figure 11b). Therefore, it is reasonable to postulate that, during the recovery process, islands of contacts are formed between opposing faces of the cracks. Such contacts would greatly impact the stiffness of the crack network without changing either N or l (e.g., Kachanov & Sevostianov, 2012;Pyrak-Nolte et al., 1990); in that sense, the crack densities estimated assuming open circular cracks represent "effective" crack densities and, hence, will not correspond to the actual microcrack geometry in the samples.
To verify if the formation of contact islands could explain the amplitude of the recovery observed in our experiments, we use the results from Trofimov et al. (2017) to compute a so-called "crack island factor" which The island factor M links the effective compliance added to a medium by a population of cracks of density ρ 0 containing islands of contact to its equivalent population of open penny shaped cracks of density ρ* such that The model of Trofimov et al. (2017) is based on the assumption that the crack population is comprised of non-interacting penny shaped cracks containing a single contact, embedded in an homogeneous linear elastic matrix. In the case of a contact island of radius a located in a penny shaped crack of radius l, the M factor is expressed in terms of λ = c/l, where c = l − a that relates to the island size, and a term β that relates to the island position in the crack, such that and By assuming the island of contact sits in the center of the crack, the β term vanishes. By further assuming that the recorded crack reduction from its maximum value during deformation (ρ 0 ) is due to the formation of contacts only (no cracks are removed from the medium, i.e., N and l are constant), λ can be computed. From our estimates of λ, we can access the ratio of contact area over total crack area in our samples The contact area ratios computed with the axial crack density data are gathered in Figure 14. During hydrostatic recovery, the initial contact area increases with increasing confinement (Figure 14a and 14b). At P c = 3 and 20 MPa, the initial contact area is zero due to the decrease in P c between deformation and recovery which induced an increase in crack aperture. At P c = 40 MPa and higher, the initial contact area ratio is non-zero and increases with recovery confinement. This is due to the removal of the differential stress and increase in pressure between deformation and recovery, which elastically closes part of the crack population. At all tested confinements, the contact area increases as the logarithm of time and the rate of increase depends positively on confining pressure.
The contact area also increases over time in the presence of a constant differential stress, regardless of its magnitude (Figure 14c). However, the magnitude of the differential stress does correlate with the initial value of the area ratio, with samples subjected to higher stresses showing lower initial contact areas. At Q/ Q max = 95% and 97.5%, the initial contact area is zero possibly due to the differential stress initially maintaining cracks open. Additionally, we observe an initial delay in the contact area growth in the form of a plateau.
During the incremental differential stress experiment, contact growth behaves in a manner analogous to that of the constant stress experiments (Figure 14d). Contact area increases with time at all stresses, however, at higher stresses the amplitude of the growth is reduced. At Q/Q max = 75%, the area increase with time is almost log-linear, whereas at higher differential stress the increase only commences after some time delay during which the area remains essentially constant.

On the Permanency of Time-dependent Crack Mechanical Recovery
Throughout the recovered microstructures, we observed that mechanical crack recovery is marked by a sizable crack aperture reduction of all cracks from the post deformation state, especially at high confinement (Figures 9 and 10). Nevertheless, while the difference between the sample left to recover at P c = 20 MPa and that left at P c = 150 MPa was striking, the difference between the sample left at P c = 150 MPa for 9 days and the sample left at the same confinement for less than one hour was barely discernible. Similarly, the crack porosities calculated from the pressure survey data showed that porosity was higher in a sample that had undergone recovery (0.36%), than in a sample that had simply been subjected to the deformation stage (0.33%, Figure 1).
Additionally, we postulated in the previous section that such an aperture reduction should naturally, given the inherent roughness of the cracks, lead to the formation of contact islands within the crack population.
With these observations, one can speculate whether the contact area created during mechanical recovery survives decompression to leave a permanent signal in the wave speeds and microstructures. In this section, we address that question by comparing data from two depressurization experiments and unloading data from the incremental differential stress experiment.
As mentioned in Section 2.4, a sample was decompressed step-wise after deformation followed by 9 days of recovery at P c = 150 MPa down to room pressure in 10 MPa steps. In this case, the radial P wave speed decreased in the sample from its post recovery value of 5,400 m/s down to 3,600 m/s at P c = 10 MPa (Figure 15). Another sample was deformed and the confinement immediately (less than 10 s) increased to 150 MPa from where it was similarly decreased step wise (i.e., no extended recovery period). In this sam- ple, wave speed decreased from 5,100 m/s at P c = 150 MPa to 3,500 m/s at P c = 10 MPa. We observe a stronger pressure sensitivity of the radial P wave speed after mechanical recovery compared to that after deformation. Such an increase in pressure sensitivity often marks the presence of cracks with higher aspect ratio (i.e., longer, thinner cracks; see for instance Jaeger et al. (2007, Chap. 8)) and generally hint at a greater number of cracks being reopened. Moreover, at 10 MPa confinement, the wave speed difference between the post-recovery and the post-deformation sample is only 100 m/s. Extrapolation of the wave speeds to atmospheric pressure would suggest that there would be little or no difference when the samples were fully depressurized. This quantitative observation is directly corroborated by the microstructures which showed essentially identical fabrics between the two samples (Figures 10d and 10f). Therefore, it appears that mechanical recovery does not leave an unequivocal permanent imprint in post-mortem samples and that the differences observed in Figure 10 were due to the pressure increase only rather than the time spent at high pressure.
During the incremental differential stress recovery experiment, between consecutive stress steps, the sample was partially unloaded to roughly 10% of Q max between consecutive stress steps, before subsequent reloading to the next step (Figures 16a and 16b). During this unloading phase, the radial P wave speed in the sample varies in a non-monotonic fashion. For instance, between the stress steps Q/Q max = 75% and Q/Q max = 85%, over more than half of the unloading process, from Q = 160 MPa down to Q = 60 MPa, the P wave speed drops by around 50 m/s from its post-recovery value of 3,400 m/s, and then increases by 50 m/s over the remainder of the unloading process (Figures 16c and 16d). Concomitantly, volumetric strain decreases by 0.2% (dilation) between Q = 160 MPa and Q = 60 MPa, and then increases by 0.06% (compaction) over the remainder of the unloading process (Figures 16e and 16f). Upon reloading the sample, the P wave speed continues increasing and reaches a maximum value of 3,500 m/s at a differential stress Q = 70 MPa. With further loading, the P wave speed decreases again and reaches 3,300 m/s when Q equates to the target differential stress of 85% of Q max . During the reloading of the sample, volumetric strain increases by a further 0.035% until Q = 100 MPa and then drops by 0.09% until Q reaches the target value of 85% of Q max . Similarly, while unloading the sample, P wave speed anisotropy increased by 0.03 from its initial value at the beginning of the recovery period at Q/Q max = 75% (Figures 16g and 16h). Upon reloading the sample, anisotropy decreased and reached its level prior to the recovery period at Q/ Q max = 75% (i.e., ϵ Thom − ϵ Thom0 = 0) when the target differential load was reached.
This particular behavior of the radial P wave speed and volumetric strain during the stress-stepping process is observed between all stress steps (Figure 17a and 17b). The amplitude of both the P wave speed and volumetric strain reduction during the unloading of the sample is comparable regardless of the differential stress during the preceding recovery hold time and generally occurs over the first 0.15% axial strain of the unloading process. Interestingly, such behavior is not observed during the unloading-loading phase immediately after the deformation stage of the experiment, prior to any recovery. In fact, after deformation, the P wave speed in the sample increases monotonically by almost 250 m/s during the unloading process and the volumetric strain by almost 0.4%. It is worthwhile noting that the radial wave speeds recorded here is comparatively stress insensitive. In fact, studies in granite have shown that crack wave speed anisotropy almost entirely vanishes upon removing the differential stress (Passelègue et al., 2018). This difference might be rooted in the semi-brittle nature of the crack networks created here, but the question remains open.
This unusual P wave speed drop upon unloading the sample after recovery indicates that part of the contact area created during the recovery hold time is lost during the stress perturbation. Moreover, this drop appears to be insensitive to the differential stress acting on the sample during the preceding recovery period. It is reasonable to assume that the recovery hold time is the main control of this phenomenon, but more experiments would be necessary to prove or disprove this hypothesis. In conclusion, a large part of the mechanical crack recovery is reversible by decompression. When samples are subjected to differential stress conditions, cracks can close mechanically, generating contact points that increase the effective stiffness of the rock, but small perturbations of the stress state partly destroy those created contacts. Moreover, the imprint of mechanical recovery in the microstructures is practically indiscernible. Indeed, the effect it had on crack porosity in our samples was so faint that it was lost within the natural sample variability (Figure 1). For this reason, we believe it would not be possible to discern the potential occurrence of mechanical recovery from in-situ field samples alone. Figure 16. Differential stress (A.), radial P wave speed (C.), volumetric strain (E.) and ϵ Thom evolution (G.) against time during the Q/Q max = 75% and Q/Q max = 85% stress steps of experiment CMh20. Zoomed view of the stress (B.), radial P wave speed (D.), volumetric strain (F.) and ϵ Thom evolution during the unloading-loading process between Q/Q max = 75% and Q/Q max = 85% (see Section 2.3.4 for more details on the differential stress stepping process). From left to right, the dashed lines in B, D, F, and H represent the starting point of the unloading process, the starting point of the loading process and the completion of the loading process. Numbers adjacent to the curves represent the differential stress acting on the samples during the recovery hold time in % of the sample Q max .

Brittle to Plastic Creep Transition
We saw that rather than completely hindering wave speed recovery, the application of a differential stress to samples during the holding time instead introduces a time delay to the recovery in the form of an initial plateau or even a small decrease at the highest differential stresses (Figures 5c and 6c). When we plot the changes in radial P wave speed against the creep strain rate derived from the Linear Variable Differential Transformer (LVDT) data (Figure 18), two creep regimes emerge: A fast creep regime where wave speeds decrease, and a slower creep regime where they increase. The transition between the two regimes occurs at some threshold strain rate that ranges from 5 × 10 −7 to 2 × 10 −6 s −1 .
We interpret this change in wave speed evolution at constant confining pressure as a transition from microcracking-dominated fast brittle creep to plasticity-dominated slower creep with decreasing strain rate. At higher strain rates, deformation is accommodated through diffuse brittle microcracking and this results in a decrease in wave speeds and an overall dilatant behavior of the sample. Conversely, at lower strain rates, deformation is accommodated through intragranular plasticity resulting in microcrack closure and an increase in wave speeds. It is likely that at any given strain rate, both brittle and plastic microprocesses are active but in different proportions resulting in a competition between brittle microcracking and plastic crack closure. Moreover, at the threshold strain rate, the effect of both microprocesses on wave speed is perfectly balanced. Such a phenomenon was previously observed in wet shales (Geng et al., 2018) where it was imputed to chemical healing but here we show that it can occur even when only mechanical processes are active.

Mechanical Recovery Processes
Mechanical recovery produces a log-linear increase in wave speeds through an increase in contact area within cracks. This observation is phenomenologically consistent with contact "aging" observed along rock interfaces and commonly cited in rate-and-state interface friction laws (e.g., Dieterich, 1972;Dieterich, 1978;Scholz & Engelder, 1976). Here, we propose that the contact formation and growth is itself driven by the viscous deformation of the bulk rock. In our experiments, in the presence of very low externally Figure 17. Changes in radial P wave speed (A.) and changes in volumetric strain (B.) against axial strain during the unloading process between holding periods of experiment CMh20 (A.). Numbers adjacent to the curves represent the differential stress acting on the samples during the recovery hold time in % of the sample Q max . The black curve marked post-def is the data obtained during the unloading phase between the deformation stage and the first recovery hold time. Figure 18. Changes in radial P wave speed against creep strain rate during the holding time of the constant stress experiments (A.) and during the incremental stress experiment (B.). Numbers adjacent to the curves represent the differential stress acting on the samples during the recovery hold time in % of the sample Q max . Arrows indicate the transition between decreasing and increasing wave speeds. applied stresses (experiments at P c = 3 and 20 MPa), we recorded sizable wave speed recovery, suggesting that one potential mechanism for the closure of cracks and the formation of contacts is the relaxation of internal, grain-scale residual stresses accumulated during deformation, for instance due to internal friction (e.g., Lawn, 1998) or dislocation glide (e.g., Hansen et al., 2019). Such relaxation could notably occur by the means of grain boundary sliding or time-dependent crack backsliding, both of which are known to be operative at low pressures and temperatures (e.g., Brantut, 2015;Scholz & Kranz, 1974). Conceptually, this approach is compatible with slow dynamics models where a superposition of different relaxation microprocesses results in an overall log-linear response at the macroscopic scale (Ostrovsky et al., 2019;Shokouhi et al., 2017;Snieder et al., 2017) At higher confinement, we observed a strong pressure dependency of the mechanical recovery process, with the recovery being greater at higher confining pressure. This result suggest that mechanical recovery can also be driven by externally applied stresses. In this case, we propose that cracks close by the means of dislocation motion around crack tips, twinning and other low temperature plasticity microprocesses (Fredrich et al., 1989). Unfortunately, our current data do not allow us to discriminate these processes at the microscale.
The effect of plastic processes such as dislocation glide or twinning on wave speeds is mostly unknown, however, we can explore whether viscous deformation of the rock medium could lead to crack closure under our laboratory conditions. The geometrical evolution of a spheroidal cavity embedded in a viscous matrix under remotely applied stress has been calculated by Budiansky et al. (1982). Although a solution for a non-linear viscous matrix was derived by these authors, its application is limited to initially spherical pores and not presented in a closed form. Closure of non-spherical voids in power-law materials generally depends on initial void orientation (Lee & Mear, 1992), and numerical solutions that include void shape tracking remain challenging to obtain in a tractable form (Lee & Mear, 1999). Closed-form expressions only exist for the linearly viscous case, which we use here as first approximation with the aim of gaining some basic insights into the closure process.
The configuration analyzed by Budiansky et al. (1982) is that of a linear-viscous matrix of viscosity η containing a spheroidal cavity of semiaxes a and b under remotely applied triaxial stress T and S. T is assumed to be in the radial direction parallel to b and S to be in that of a. The change in shape of the spheroidal is given by (Budiansky et al., 1982): and where ξ = a/b is the aspect ratio of the spheroid, V = 4πab 2 /3 its volume, and For thin cracks (ξ ≪ 1), these expressions simplify and asymptotically become and It arises from these equations that, for compressive stress states (S, T < 0), the viscous flow of the matrix implies that both the crack aspect ratio and the crack volume decrease linearly with time and that cracks fully close within a finite time. Under hydrostatic conditions such as those in our experiments (S = T = P c ), a crack of aspect ratio ξ 0 , will have a closure time given by Consequently, linear viscous relaxation of the rock matrix can lead to the total closure of cracks. Moreover, since the closure time of a crack depends on its initial aspect ratio, the general time-dependent closure behavior of a crack population will be a function of the initial aspect ratio distribution within the population. This could explain the slower rate of crack density reduction observed at higher hydrostatic pressure in Figure 13b. Since the elastic closure pressure P closure of a crack in a medium of Young modulus E is P closure ∼ πEξ 0 (Jaeger et al., 2007, Eq. 8.309), it arises that at higher confinement, cracks with initially low aspect ratio will have closed elastically, which leaves only high aspect ratio cracks in the material, which tend to close over longer timescales. In other words, the remaining cracks in the rock at P c = 150 MPa are slower to close and therefore the crack density decrease is also slower.
To illustrate the potential viscosity involved in the crack closure process, let us now estimate the closure time of a crack of initial aspect ratio 10 −4 using Equation 19. This value of aspect ratio is commonly found in rocks and is therefore representative of our case (Kranz, 1983). For such a crack to close in 1 day at a confining pressure P c = 40 MPa, the surrounding rock matrix would need to have a viscosity of around 10 16 Pa⋅s. It is difficult to assess the plausibility of this estimate from independent data; creep flow laws for calcite have been established (e.g., Renner et al., 2002) but are only valid for steady-state, high temperature deformation. At room temperature and low pressure, calcite behaves as a power-law solid with strain hardening, and estimates of rheological parameters in that case have been obtained by Nicolas et al. (2017). Using estimates obtained from creep data on Tavel limestone (Table 3 of Nicolas et al. (2017)), the apparent viscosity of calcite at room temperature ranges from 10 19 Pa⋅s to 10 14 Pa⋅s at applied stresses from 40 to 150 MPa. Despite the considerable uncertainty of such estimates and the simplification associated with assuming linear viscosity, we find that viscous crack closure is compatible with our experimental results.

Implication for Fault Processes
Our study has demonstrated that mechanical recovery is a fast and potent wave speed recovery source, even at room temperature. At room temperature and under very low confinement, wave speed recovery can reach 6% of the damage induced wave speed drop, and at P c = 150 MPa, it can reach almost 40% of that drop in five days. Therefore, it is possible for the impact of mechanical recovery on wave speeds to exceed that of chemical healing. This is possible not only for the observed wave speed recovery around faults at the field scale, but also for laboratory-scale experiments involving both chemical healing and mechanical recovery (e.g., Aben et al., 2017).
In nature, boundary conditions during the post-seismic phase might lead to stress relaxation, which will produce elastic crack closure effects in addition to the static time-dependent mechanical recovery observed here. The net effect on wave speeds will depend on the regime and mechanisms of crack closure, and is hard to predict in general. Elastic unloading of differential stress tends to increase wave speed (Figure 17), whereas unloading of confining pressure tends to decrease them ( Figure 15). Time-dependent backsliding and creep processes systematically tend to produce wave speed recovery. Existing experimental work under stress relaxation conditions show systematic wave speed increase, so the combined unloading and creep processes seem to lead to net crack closure (Kaproth & Marone, 2014;Schubnel et al., 2005) Our data showed that mechanical recovery occurred even at low confinement and ambient temperature. It is therefore reasonable to assume that this process is active everywhere in Earth's crust. For this reason, data reported by the studies gathered in Table 2 are all likely to be, at least partially, impacted by mechanical recovery, regardless of the sampling depth of the method used. Unfortunately, a direct comparison of our data to field data is not possible and it is not yet possible to quantify the exact contribution of mechanical recovery to seismic observations yet.
Furthermore, the strong pressure dependency of the phenomenon suggests that it would be faster and more important at depth. Moreover, our microstructural analysis of samples that have undergone mechanical crack closure revealed that it involves confining pressure-sensitive crack aperture reduction. Past studies have shown that chemical healing of cracks is highly sensitive to crack geometry, having a faster kinetics in narrow cracks (Brantley et al., 1990;Hickman & Evans, 1987). Therefore, we expect mechanical crack closure to have a positive feedback on chemical healing of crack damage, accelerating even further damage healing at depth. These results tend to imply that co-seismic damage at depths where hydrothermal conditions are found could heal very rapidly, in a matter of hours or days, and much faster than the overall timescale of the inter-seismic phase.
Additionally, such dramatic changes in elastic properties need to be taken into account when modeling post-seismic deformation and fault expression at the surface. Even more so as our experiments suggest that creeping rocks can undergo substantial stiffening, making our data set particularly relevant to creeping faults and post-seismic creep modeling.
Nevertheless, several important questions concerning time dependent mechanical recovery of rocks remain unanswered by our study. Wave speed recovery has been recorded at laboratory timescale (several hours) in different rock types and in settings where chemical healing was unlikely to occur. For example (Schubnel et al., 2005), reported wave speed recovery of about 150 m/s in dry damaged Solnhofen limestone (Kaproth & Marone, 2014), recorded wave speed recovery of more than 1,000 m/s in brine saturated halite gouge, and (Brantut, 2015) recorded wave speed recoveries of about 75 m/s in decane-saturated Purbeck limestone. These observations show that mechanical recovery is not limited to Carrara marble. While we suggest that the micromechanics of crack recovery would remain essentially the same in a different type of rock, the intrinsic properties of another rock type such as viscosity or young modulus, will likely impact both the time-scale and amplitude of the mechanical wave speed recovery. The closure time of a cavity in a viscous matrix is directly proportional to the matrix intrinsic apparent viscosity (Equation 19), and is significantly impacted by power-law creep exponents and parameters (e.g., Lee & Mear, 1999), therefore using a rock with a different viscosity than Carrara marble would produce quantitatively different results.
Similarly, an increase in temperature is expected to change the properties of rocks, thus altering the amplitude and timescale of mechanical crack closure; particularly so when the rheology of the rock considered is primarily controlled by plasticity. Crystal plastic processes are thermally activated, so that elevated temperatures are often associated with a reduction in rock viscosity (e.g., Kohlstedt & Hansen, 2015, and references therein). In the framework of our time dependent viscous closure of cracks (Equation 19), a reduced viscosity implies a reduced closure time and consequently a faster mechanical crack closure. One important parameter that was not investigated here is pore fluid pressure, which is also expected to impact mechanical crack closure. As mentioned above, while it is known that the presence of water will activate chemical healing of cracks which is itself accelerated by mechanical closure, our modeling of the micromechanics ignores the potential effect of pore pressure. Further work with rocks in the presence of pressurized water is needed to fully address these questions.

Conclusion
We investigated the mechanical crack closure microprocesses (as opposed to chemical healing) by the means of hydrostatic and triaxial recovery experiments on deformed dry Carrara marble at room temperature during which we held samples under constant stress conditions for extended periods of time. In all our experiments, elastic wave velocities first increase, and then drop during deformation due to microcrack growth, and recover significantly, by up to several tens of percent, during hold periods of the order of a few days. Microstructural analysis showed that mechanical recovery involves a reduction of crack aperture and the generation of contact points along the microcrack faces. Mechanical recovery is sensitive to applied stresses, and is quantitatively stronger at elevated confining pressure. Time-dependent crack closure and wave speed recovery can be explained by viscous flow of the rock. Perhaps less expected is the observation that significant recovery also occurs after stress removal. We suggest that crack closure at decreasing stresses is due to relaxation of internal stresses accrued during deformation, for instance due to internal friction.
Mechanical recovery is driven by applied stresses, such as lithostatic pressure, and is expected to be facilitated at elevated temperature due to the thermal dependency of rock viscosity. However, as reported by Brantut (2015); Scholz and Kranz (1974), internal stress relaxation is already sufficient to produce significant recovery (up to 10% of the wave speed drop) in only a few days. Therefore, we expect that mechanical recovery can be a prevalent mechanism in the time-dependent changes in fault elastic properties following earthquakes across the whole seismogenic layer. By promoting crack aperture reduction, mechanical closure should also facilitate chemical healing or sealing processes, making earthquake-induced damage short lived.
In the future, the impact of temperature and pore pressure on mechanical recovery will be assessed to better constrain the phenomenon. Additionally, a similar study in a different rock type will be conducted to further strengthen our conclusions.