Origin of the Co‐Seismic Variations of Elastic Properties in the Crust: Insight From the Laboratory

Abstract Seismological observations highlighted that earthquakes are often followed by changes in elastic properties around the fault zone. Here, we studied the origin of these variations using stick‐slip experiments on saw‐cut granite samples presenting different degrees of bulk damage (i.e., microcracks). Stick‐slip events were induced under triaxial compression configuration with continuous active ultrasonic measurements at confining pressures representative of upper crustal conditions (15–120 MPa). Both the P‐wave velocity (VP) and amplitude (AP) showed drops, concurrently with stress drops, and had a non‐monotonic dependence toward the fault's stress state. Our experimental results suggest that co‐seismic changes in VP were mostly controlled by the elastic re‐opening of microcracks in the bulk, rather than by co‐seismic damage or the formation of fault gouge. Co‐seismic changes in AP were controlled by a combination of elastic re‐opening of microcracks in the bulk and inelastic processes (i.e., co‐seismic damage and gouge formation and dilation).

response of fault zones to loading in terms of seismic properties is expected to vary spatially and temporally, and to be a function of both fault structure and travel paths of the seismic waves (Nishizawa, 1982).
To get insights on co-seismic seismic properties variations throughout the seismic cycle, several experimental studies focused at monitoring the evolution of elastic properties through laboratory friction experiments on artificial faults (Kendall & Tabor, 1971). Yoshioka and Iwasa (2006) already used transmission waves to monitor a brass fault contact evolution under normal and shear stress, finding a clear increase in wave amplitude with the increase of normal and shear stresses and amplitude variations linked with precursory slip due to change of the fault's contact area. Following studies performed with gouge interfaces (Kaproth & Marone, 2013;Scuderi et al., 2016;Tinti et al., 2016) showed both co-seismic and precursory changes in P-wave velocity associated with laboratory earthquakes, attributed to the gouge layer dilation and its change of porosity. Scuderi et al. (2016) explored the complete spectrum of failure modes, from slow to fast earthquakes, showing that not only co-seismic changes but also precursory variations of P-wave velocity occur for each mode of failure. Fukuyama et al. (2018) studied amplitude variation during high-velocity friction experiments. Moreover, Shreedharan et al. (2021) showed clear precursory P-wave amplitude variations occurring with the instability nucleation phase and precursory P-wave velocity variations distorted by the presence of the surrounding bulk material. These observations suggest that the elastic properties of the bulk material surrounding the fault may play a role in seismic velocity drops associated to natural earthquakes, as well as its recovery in the months following the rupture. Indeed, seismic waves velocities are sensitive to a change in the degree of damage (i.e., presence of microcracks) of the medium they travel through (Blake et al., 2013;Brantut, 2015;Griffiths et al., 2018;Guéguen & Palciauskas, 1994;Kuttruff, 2012;Nasseri et al., 2009;Nishizawa, 1982).
Our study aims at understanding how much of the change in seismic properties observed during earthquakes is controlled by co-seismic damage occurring on-(i.e., gouge production) and off-(i.e., formation of microcracks in the fault wall due to seismic rupture) fault, and how much is instead affected by the presence of the initial degree of damage characterizing the bulk material and its response to stress changes. To this end, we conducted stick-slip experiments (Brace & Byerlee, 1966) under a wide range of confining pressures on granite saw-cut cylindrical samples presenting two different degrees of initial bulk damage, to mimic different fault damage zone properties.

Materials
The tested material is La Peyratte granite, a crustal rock presenting a modal composition of 38.5% plagioclase, 28.5% quartz, 20% K-Feldspar, 13% biotite with an average grain size of 800 µm.
Right-circular cylindrical samples were prepared with 38 mm diameter and 78 mm height. Some were thermally treated before the experiments by slowly heating them (5°C/min, to avoid thermal gradients inside the sample (Wang et al., 2013)) up to different target temperatures (650°C) and let cool down to ambient temperature inside the oven overnight, to avoid thermal shock. Target temperatures were chosen above the α-β quartz transition (572°C), allowing intense intra-granular cracking, randomly oriented in the bulk, producing isotropically damaged media (Glover et al., 1995;Pimienta et al., 2019;Wang et al., 2013), with reduced fracture toughness (Kang et al., 2020;Nasseri et al., 2007). To characterize the different samples, density and porosity were measured, obtaining densities of 2.63 g/cm 3 and 2.58 g/cm 3 and porosities of 0.4% and 6.6%, respectively, for non-treated and thermally treated granite at 650°C.
Samples were saw-cut with an orientation of 30° to the vertical axis, creating an artificial fault plane. The fault roughness was imposed by hand using #240 grit sandpaper, generating a smooth fault, optimally oriented for reactivation, avoiding the propagation of new secondary fractures in the surrounding medium (Renard et al., 2020). The lack of secondary fracture formation under this configuration has been verified in previous experimental work (e.g., Acosta et al.

Testing Procedure
Tests were run in an oil medium high-pressure triaxial apparatus, FIRST (installed at LEMR, EPFL). The samples were first submitted to a target confining pressure ( C P ) (15,30,45,60,90,and 120 MPa), with a subsequent increase of axial load. Axial load was applied by controlling the oil flow rate (0.25 or 0.50 ml/min in few cases), pushing the piston, generating a displacement rate of ∼6 · 10 −6 mm/s. For the different samples (non-treated and treated), experiments were conducted starting from the lowest C P and, once the stick-slips series was performed, C P was increased to the following target C P and a new stick-slips series performed, up to the highest target C P . Two displacement transducers were placed beside the sample, measuring locally the sample' shortening and/or the fault slip. Mechanical data were recorded at a frequency of 100 Hz for the whole duration of the tests.

Acoustic Measurements
Active acoustic measurements were recorded during deformation, using acoustic sensors (PZT crystal) placed inside the top and bottom anvils of the triaxial apparatus, with a recording frequency of 100 Hz. The acquisition system setup and the picking procedure were modified and adapted from Acosta and Violay (2020) (refer to the supporting information for details).
Seismic waveforms were used to measure the evolution of P-wave velocity and amplitude along the experiment. Once detected the P-wave arrival time ( P t ) the P-wave velocity ( P ) V was computed as with corrected L the length of the sample, systematically corrected by the elastic shortening and slip occurring. The P-wave amplitude ( P A ) was computed as the difference in amplitude between the first maximum value and minimum value of the P-wave (Figure 1b, inset). Seismic measurements were performed in the vertical direction, parallel to the sample axis (ray path showing the largest variations in wave velocity due to the mechanical anisotropy occurring during differential loading) (refer to the Supporting Information, Figure S1). PAGLIALUNGA ET AL. , ,pulsing direction (in red) and strain gauge location. (b) P-wave arrival time detection; the top panel displays the wave energy evolution with time, the bottom panel displays the detected P-wave arrival time (in red), and P-wave first arrival amplitude (inset, in blue). (c) Seismic waves evolution during a stick-slips series performed for a treated sample at a C P of 15 MPa. Red markers indicate the arrival time detected by the automatic picking. Shown waves are sampled (1:5). In white, the shear stress evolution during the test.

Experimental Results
Stick-slip experiments conducted under different C P were used to investigate seismic properties evolution throughout the seismic cycle. For each of them, the shear stress increased first linearly and, once reached the fault strength, dropped to a residual value ( Figure 2). As expected, the higher the applied C P , the higher the fault strength, stress drop, and resulting slip were observed. Concerning the seismic properties, an increase of P V and P A was observed during the hydrostatic loading up to the target C P . Moreover, both P V and P A responded to the applied differential stress accordingly, increasing during loading and decreasing during unloading. For both the non-treated and treated sample, the increase in P V during differential loading ( P I loading V ) was larger for low C P , and smaller for high C P ( Figure S2). In particular for the non-treated sample, was ∼4.5 × 4 10 V, ∼4.55 × 4 10 V, and ∼1.2 × 4 10 V, respectively at a C P of 15, 45, and 120 MPa. As stress drops occurred, associated to seismic fault slip, a drop in P V as well as in P A was observed.
These co-seismic drops in velocity ( P ΔV ) and amplitude ( P ΔA ) were computed for each stick-slip, and compared with their respective stress conditions ( Figure 3). P ΔV did not show a linear dependence on stress conditions applied to the fault (i.e., normal stress, confining pressure, and shear stress). In the case of non-treated sample, for low C P (15-30 MPa), hence for low  Δ (∼1-3 MPa), P ΔV were larger (∼4-9 m/s) than for events recorded at higher C P (60 MPa) and medium  Δ (∼4-9 MPa), which were ∼2-6 m/s.

Discussion
In our experiments, a non-monotonic P ΔV evolution with shear stress drops was observed (Figure 3), suggesting that distinct physical processes coexist at the origin of velocity changes during stick-slip instabilities, due to the combination of initial bulk damage and loading conditions. In particular, these drops in velocity during stick-slip events could be related to (i) horizontal microcracks re-opening in the bulk, after initial closure during increasing differential stress (Passelègue et al., 2018), due to differential stress reduction or (ii) co-seismic damage induced around the fault during dynamic rupture propagation and fault motion Okubo et al., 2019).
To test these hypotheses, we estimated the maximum possible contribution of microcracks re-opening due to co-seismic stress drop on the associated P ΔV . Such effect is expected to be similar during both loading and unloading of the bulk (if no adhesion is considered on the microcrack, for example, stress-induced microcrack opening/closure is a reversible process). Under these assumptions, P I loading V for each C P can be used to estimate the contribution of microcracks opening following co-seismic stress drops and associated strain release, not considering possible co-seismic damage occurring off-fault. We predicted P ΔV only due to the re-opening of microcracks occurring in the bulk as follows: where  Δ ax is the drop in axial strain measured concurrently with stress drop,  I loading ax the increase in axial strain during differential loading (strain gauge located in the bulk material, far enough from the fault, expected to capture elastic deformation of the bulk). P Δ predicted V for all the events at each Pc for both treated and non-treated samples, showed the same evolution with loading conditions of the ones experimentally observed   P Δ measured V . In fact, by plotting them together ( Figure 4a), a linear dependence between the two is noted, with a slope very close to 1:1, indicating that P Δ measured V are well explained by the co-seismic re-opening of microcracks in the bulk, resulting from the release of strain. This suggests that in our experimental configuration, no significant co-seismic damage was generated during rupture propagation, or it was negligible with respect to the observed velocity variations. Once again, the non-monotonic trend observed as a function of applied stress (Figure 3) is explained by the interplay between C P and P I loading V . For low C P the induced stress drops are very small (∼1-3 MPa/∼3-4 MPa, respectively, for non-treated and treated sample) but the related P I loading V (seen here as the maximum potential velocity drop caused by microcracks opening, at a specific C P ) is very large (∼200/∼390 m/s), generating PAGLIALUNGA ET AL.
Remarkably, while co-seismic P ΔV proved, in our experiments, to be mostly related to the re-opening of microcracks in the bulk and elastic relaxation, we could observe some gouge production on the post mortem samples fault surfaces (in particular on the one of the treated sample), which is expected to have an influence on the seismic properties measured across the sample Shreedharan et al., 2021;Tinti et al., 2016). In particular, the presence of a gouge layer is expected to affect P A , considered as a simplified way to account for attenuation (Lockner et al., 1977) (i.e., the higher the amplitude of the wave, the lower the attenuation and vice versa). Compared to P V , which is mainly affected by elastic processes such as microcracks closure and re-opening, P A is also influenced by the fault's specific stiffness and by the inelastic, dissipative deformation processes occurring on and off-fault (i.e., frictional sliding of microcracks in the bulk and/or gouge particle shearing).
A prediction similar to the one described above (Equation 2) was tempted to test if P A measured in these experiments was also mainly influenced by the bulk properties and stress conditions. Equation 2 was modified and P V was replaced by P A as follows: with P Δ predicted A the predicted amplitude drops and P I loading A the overall increase of P A during differential loading. For the non-treated sample, the prediction works well, with values falling very close to the prediction line of slope 1:1 (Figure 4b). However, for the treated sample, this is true only for the lower C P (15 MPa). For higher C P (45 and 120 MPa), the predicted drops do not mimic the measured ones, the latter being notably larger (up to 400% larger). This might be explained by the change in the fault's contacts and/or by non-elastic processes occurring either in the bulk (i.e., friction caused by shear along microcracks) or on the fault surfaces (i.e., gouge production, shearing, and dilation). Since the expected stress responsible for microcracks shearing is larger than the one expected to activate shearing along the artificial fault, we assume that the non-elastic processes observed are caused by the fault's response to stick-slips. This was verified by analyzing the evolution of P A with cumulative slip ( Figure S3), since (a) gouge production is expected to increase linearly with cumulative fault slip (Archard, 1953)   under high applied stresses. Given that (a) the thermally treated sample is expected to have a lower fracture toughness than the non-treated one (Nasseri et al., 2007), (b) we observed a decrease in P A for consecutive stick-slips only under medium to high C P (i.e., normal load acting on the fault), and (c) that we could observe a large amount of gouge on the post-mortem sample's fault, we ascribed P A behavior to be a function of the gouge production (Frérot et al., 2018) and subsequent gouge particles shearing during fault's slip under these conditions. This looks coherent with the evolution of the fault specific stiffness (Pyrak-Nolte et al., 1990) for the different stress conditions ( Figure S5), which in the case of the treated sample reaches a sort of saturation for the highest C P (120 MPa)(i.e., the gouge once filled all the voids available in the interface and compacted, will not deform any further for higher C P , not influencing F k ).

Implications and Conclusions
Summarizing our interpretation of the results, co-seismic P ΔV seem to be controlled by the combination of bulk properties and applied stress (i.e., re-opening of the microcracks present in the bulk concurrently with stress drop). This does not imply that other phenomena occurring during stick-slips, such as gouge layer dilation, could not contribute to P ΔV itself, but only that their influence, compared to one of the pre-existing microcracks in the bulk, resulted negligible in our experiments. In addition, while co-seismic P ΔA also looked to be controlled by the combination of bulk properties and applied stress when the presence of gouge was not dominant, they were probably controlled by dissipative processes occurring on-fault when the conditions (treated sample and higher applied stress) allowed an important production of gouge, hence a necessary shearing of gouge particles.
Conversely to previous experimental studies (Kaproth & Marone, 2013;Scuderi et al., 2016;Shreedharan et al., 2021), no significant and clear precursory variation of seismic velocity and amplitude was observed. This might be due to several reasons; among the others, the localized nature of the nucleation phenomenon, known to be the cause of observed pre-seismic slip. Depending on the nucleation patch size, either a low or a high-stress perturbation will be induced in its vicinity. A nucleation patch length significantly smaller than fault length is expected under this configuration (Harbord et al., 2017). Assuming this, the stress release during the nucleation of instability is expected to affect only a small fraction of the whole sample, without inducing any strong premonitory change in P V or P A . Another reason could be related to our resolution of the seismic measurements, which may be not high enough to capture precursory changes which remain lost within the error linked to the present measurements.
However, our results could help to better understand in which conditions precursory variations of seismic properties can actually be detected and used to monitor fault's state of stress. It is clear that the wallrock's elastic properties have a huge control on the seismic properties measured across the system. It was recently shown that this distortion is crucial for observations of the aforementioned precursory phase (Shreedharan et al., 2021). For this reason, the luckiest combination to observe this variation would be to encounter a fault composed by a wide gouge layer and with a large nucleation patch. The shearing of gouge particles within the fault zone will strongly affect the seismic amplitude, which is the parameter the most sensitive to inelastic processes.
Moreover, even if a direct comparison remains risky given the differences in the applied conditions and the large uncertainties in the estimation, there are some analogies between the relative variations in P V recorded in this study and the ones observed after real earthquakes. It emerges that the overall range of values measured in our experiments (0.03%-0.35% and 0.16%-1.23% respectively for non-treated and treated sample), performed under stress conditions representative of the upper crust, is comparable to the ranges of values measured after real earthquakes (Brenguier et al. [2008] estimate variations of ∼0.02%-0.07%, Chen et al. [2010] find ∼0.04%-0.08%, Nimiya et al. [2017] find ∼0.4%-0.8%, Qiu et al. [2020] find ∼0.15%-0.25%, Taylor and Hillers [2020] find ∼0.15%). The similarity between our observations and the ones referring to natural earthquakes, suggests that the monitored seismic properties could be controlled by the same factors (i.e., combination of propagation of seismic waves through fault core, damage zone, wallrock). In fact, measurements performed across artificial gouge faults, monitoring P V evolution of the only gouge layer with PAGLIALUNGA ET AL.
10.1029/2021GL093619 7 of 10 no contribution of the surrounding medium, showed much higher relative P V variations ∼1%-4% Tinti et al., 2016) (for a graphical representation refer to the Supporting Information, Figure S8).
Finally, given the impossibility to measure natural seismic variations of the only fault core, monitoring the evolution of seismic velocity along faults surrounded by large damage zones, could be of interest for observing co-seismic changes during shallow earthquakes, since the combination of large and highly damage zones and low-stress conditions lead to an extremely high sensitivity in velocity changes due to stress perturbations, especially at low depths. Moreover, since many earthquakes are preceded by a nucleation stage (Latour et al., 2013;Ohnaka, 2003;Ruiz et al., 2014;Socquet et al., 2017;Tape et al., 2018), which is expected to release part of the stress along the fault, the amplitude evolution may provide, under the aforementioned conditions, some indications about stress evolution along the fault and the proximity to failure. This kind of observations could, yet, be limited by the current spatial resolution of seismological observations and by the knowledge of the damage zones in seismogenic faults.

Data Availability Statement
The raw data can be found at the following address: http://doi.org/10.5281/zenodo.4892328.