Susceptibility of Microseismic Triggering to Small Sinusoidal Stress Perturbations at the Laboratory Scale

Small transient stress perturbations are prone to trigger (micro)seismicity. In the Earth's crust, these stress perturbations can be caused by various sources such as the passage of seismic waves, forcing by tides, or hydrological seasonal loads. A better understanding of the dynamic of earthquake triggering by stress perturbations is essential to improve our understanding of earthquake physics and our consideration of seismic hazard. Here, we study an experimental sandstone‐gouge‐filled fault system undergoing combined far field loading and periodic stress perturbations (of variable amplitude and frequency) at crustal pressure conditions. Microseismicity—in the form of acoustic emissions (AEs)—strains, and stresses, are continuously recorded in order to study the response of microseismicity as a function of loading rate, amplitude, and frequency of a periodic stress perturbation. The observed AE distributions do not follow the predictions of either a Coulomb failure model, taking into account both constant loading and oscillation‐induced strain rates, or a rate and state model. A susceptibility of the system's AE response to the amplitude of the confinement pressure perturbation is estimated, which highlights a linear relation between confinement pressure amplitude and the AE response amplitude, observations which agree with recent higher frequency experimental results on dynamic triggering. The magnitude‐frequency distribution of AEs is also computed. The Gutenberg‐Richter b‐value oscillates with stress oscillations. Our experiments may help complement our understanding of the influence of low inertia stress phenomena on the distribution of seismicity, such as observations of dynamic triggering and seismicity modulation by tides or hydrological loading.

existence of periodic variations in seismicity rates in various tectonic settings, such as annual variations of seismicity in continental collision zones such as the Himalayas (Ader & Avouac, 2013;Bettinelli et al., 2008;Bollinger et al., 2007), in faults located near subduction zones as in Taiwan (Hsu et al., 2021) and Japan (Heki, 2003), in non-deforming intraplate regions such as the New-Madrid Seismic Zone (Craig et al., 2017), and even in the case of deep-focus earthquakes (Zhan & Shearer, 2015). Observations of tidal variations of seismicity have also been made in various seismotectonic settings, including along shallow thrust faults (Cochran et al., 2004), as well as in seismicity associated with geothermal (Wang et al., 2022) and submarine volcanic activity (Tolstoy et al., 2002). Tidal forces have also been linked to the triggering of tectonic tremors (e.g., Chen et al., 2018;Ide & Tanaka, 2014;Rubinstein et al., 2008) and slow slips (Hawthorne & Rubin, 2010). Moreover, the dynamic triggering of earthquakes by transient oscillatory stress changes caused by seismic waves has largely been accepted (e.g., Anderson et al., 1994;Brodsky & van der Elst, 2014;Hill et al., 1993).
All these observations are difficult to interpret due to differences between studies, with significant changes in tectonic contexts, fault geometry, and oscillation frequency and amplitude. Theoretical and numerical studies have attempted to unify all these observations, invoking nucleation times relative to the period of the stress oscillations as one of the factors differentiating a stress controlled regime from a stressing rate controlled regime (Ader & Avouac, 2013;Beeler & Lockner, 2003;Dublanchet, 2022;Heimisson & Avouac, 2020;Perfettini et al., 2001), or suggesting magnitude-dependent and oscillation-geometry-dependent modulation (Pétrélis et al., 2021). Experimental studies focused on the links between stress oscillations and seismicity have investigated the increased synchronization of the temporal distribution of AEs with periodic stress oscillations prior to macroscopic failure (Chanard et al., 2019;Noël, Pimienta, & Violay, 2019) and either the triggering or the induced clock-advance of stick-slips (Bartlow et al., 2012;Chelidze et al., 2010;P. A. Johnson & Jia, 2005; P. A. Johnson et al., 2008;Lockner & Beeler, 1999;Savage & Marone, 2007. To the best of our knowledge, no study has investigated the influence of stress oscillation parameters on AEs distributions in granular shear experiments, such that observations of seismicity modulation in different natural tectonic contexts lack sufficient experimental references. Furthering our understanding of the role of small periodic stress changes on earthquake triggering could have implications for earthquake hazard assessment. This has been highlighted by observations of magnitude-frequency distribution (b-value) variations due to tides (Ide et al., 2016) as well as by observations of gradually increased correlation of seismicity with tides before large megathrust earthquakes, which could signify that seismic response of faults close to failure grow more susceptible to stress variations (Tanaka, 2010(Tanaka, , 2012. These latter observations have, however, been called into question for their statistical significance (Wang & Shearer, 2015) although experimental results support the validity of the theory on which they are based (Chanard et al., 2019;Noël, Pimienta, & Violay, 2019). These elements act as further motivation to experimentally investigate the question.
In this study, we perform shearing experiments on a granular fault gouge analog at different strain rates, in order to mimic different tectonic loading rates. Confinement pressure oscillations are applied and AEs within the gouge are monitored. The magnitude and phase of these AEs is then computed and analyzed. We start by presenting our experimental set-up and the range of parameters explored within our study. We then present experimental results in the form of stress and friction evolution with slip, and temporal distribution variations and magnitude properties of AEs. The modulation of the temporal distribution of the AEs is then quantified, and both the influence of the loading and oscillation parameters, and the variations of magnitude-frequency distributions with the stress oscillations are discussed. Finally, we consider an upscaled interpretation of our results and address their implications for seismology.

Material and Sample Preparation
Experiments are performed on matching aluminum cylindrical guide blocks ("saw-cuts") featuring corrugated planar surfaces oriented 30° from the axial direction ( Figure 1a and Figure S1 in Supporting Information S1). In order to mimic naturally occurring gouge-filled faults, we place 20 g of gouge-the equivalent of a roughly 5 mm thick layer-between the saw-cuts, which we enclose with aluminum tape. The gouge is generated by crushing pieces of Fontainebleau sandstone-a well-studied sandstone known for its purity and its homogeneous composition of quartz grains of characteristic size 250 μm (Bourbie & Zinszner, 1985)-down to a fine polydisperse gouge, with the largest grain size matching the characteristic grain size ( Figure S2 in Supporting Information S1). The experiment is designed with an artificial fault with gouge to produce enough AEs for statistical significance of the AEs catalog analysis, a criterion not met by experiments using bare rocks as fault analogs. The use of gouge to investigate microseismicity experimentally is not a novel approach (e.g., Bolton et al., 2019Bolton et al., , 2022P. A. Johnson et al., 2013;Mair et al., 2007).

Sample Instrumentation
The sample is equipped with both acoustic sensors and strain gauges to monitor AEs and deformation during experiments. Double component strain gauges (FCB-2-11, 120 Ω by Tokyo Sokki Kenkyujo Co., Ltd, Japan) are glued on each saw-cut to measure strain in both axial and radial directions. The saw-cuts are then inserted into a neoprene jacket. Eight acoustic transducers, consisting of a piezoelectric (PZT) crystal (PI CERAMIC PI255, 5 mm diameter, 2 mm thickness, 1 MHz central frequency) sensitive to P waves (normal to the sample-sensor surface) are encased within brass holders. These are glued directly onto the lower guide block surface, in a wide formation depicted in Figure 1b and Figure S3 in Supporting Information S1. More information on the sensors can be found in Brantut et al. (2011). The array is designed to achieve homogeneous coverage, with sensors as close as possible to the gouge layer to detect the smallest possible AEs. Finally, to ensure proper sealing from the pressurized oil present in the confining chamber, a bi-component epoxy resin is applied on the neoprene jacket around the holes punctured within the jacket for strain-gauges and acoustic transducers inserts.

Apparatus
The sample is inserted into a triaxial apparatus at the École Normale Supérieure in Paris. One important feature of this experimental system is the presence of an auto-compensated confinement chamber, that is, a chamber which counterbalances the confining pressure exerted on the piston. More details can be found in Schubnel et al. (2005). The auto-compensation process depends on both the main confinement chamber and the auto-compensation chamber being at the same pressure. Therefore, changes in the confining pressure need to be gradual in order for the pressures in the main confinement chamber and the auto-compensation chamber to equilibrate.
Piston displacement, applied axial stress (σ 1 ), confining pressure (P c ), as well as differential stress (Δσ = σ 1 − P c ), radial and axial strains measured by strain gauges are recorded with a sampling rate of 10 Hz with a digital encoder (HBM MGC Plus). Resolution on acquired stresses is ∼10 −3 MPa and 10 −6 on strains.
From these measurements, the shear stress τ and normal stress σ N applied in the oblique direction of the gouge layer can respectively be derived as follows: (1 + cos (2 )) + where Δσ is the differential stress, α is the angle between the axial direction and the normal to the gouge plane (here α = 60°) and P c is the confining pressure.
The piston displacement is measured with a Linear Variable Differential Transformer (LVDT) placed atop the piston with a ∼0.1 μm resolution. Throughout the experiments, the piston and the aluminum saw-cuts undergo elastic shortening. The measure given by the LVDT as the piston advances therefore does not correspond to the actual slip occurring within the gouge. Slip δ is determined by correcting for the elastic contributions of both the machine and the saw-cuts, as follows: where d p is the raw displacement measured on the top of the piston by the LVDT, k (m/N) is the machine compliance, F is the applied load (equal to F = Δσπr 2 , with r the sample radius), and L and E al are the aluminum sample length and Young's modulus respectively.

Experimental Procedure
Once sealed inside the confinement chamber, the confining pressure exerted on the sample is raised and kept overnight at 90 MPa in order for the gouge to undergo compaction and comminution. During this stage, several thousands of acoustic emissions (AEs) are detected, showcasing the effective compaction of the gouge via grain crushing and grain contact rearrangement. The confining pressure is then lowered to 50 MPa-corresponding to depths of approximately 2 km under typical lithostatic pressure conditions-for the rest of the duration of the experiment.
The sample is then loaded axially, the piston displacement being hydraulically servo-controlled via a high-pressure syringe pump. Displacement rate is controlled as a proxy by oil injection rate into the piston upper pressure chamber. A high injection rate segment of 50 cc/hr (corresponding to a piston velocity V p of roughly 1 μm/s) is first applied for a displacement of 1.5 mm, to elastically load the system, compact the gouge further and approach a frictional steady-state. We then perform constant injection rate segments, first lowering the velocity three times from around 1 μm/s down to a minimum of approximately 10 nm/s, then raising the velocity by repeating the use of the previous velocity segments in reverse order. A displacement of 0.5 mm is allowed at each velocity segment, and the segment duration t segment is therefore determined by the loading velocity. For each segment, the slip velocity V is also computed (see Table 1). The slowest segment, which lasts approximately a full day, is not repeated twice but is lengthened to ensure the measured displacement during the segment is twice that of faster segments. These velocity segments are undertaken to characterize the frictional stability behavior of the gouge at these strain rates, as well as to investigate the influence of loading rate on the temporal distribution of AEs.
In total, six experiments are conducted following this axial loading procedure. A first experiment is conducted by simply following the loading procedure and is used as the control experiment. During the remaining five experiments, oscillations are introduced in the confining pressure, as detailed in Table 2. The investigated oscillation amplitudes correspond to 5%, 1%, and 0.2% of the imposed confining pressure depending on the experiment.
Confining pressure oscillations are performed via a second servo-controlled hydraulic syringe pump. The oscillations are quasi-sinusoidal (by setting separately the amplitude, rate, acceleration, and period of the pressure oscillation). Three important experimental limitations need to be pointed out: The difference in average slip velocity for segments of identical injection rates is due to the elastic response of the saw-cuts and apparatus, which accommodates part of the displacement of the piston, as well as due to the gouge sample maturing throughout the experiment. b Segment 7a is the last segment of the control experiment. c Segment 7b is the last segment of all oscillation experiments. 1. The confining pressure oscillations generate shear stress oscillations, albeit smaller, due to elastic couplings and to the geometry of the experimental setup. 2. It is impossible with our system to perform confining pressure oscillations, while also maintaining the confining pressure average constant. In consequence, the experiments with oscillations are performed under "constant pressurized volume" conditions, which lead, because of oil leakage, to a small pressure drop at the beginning of our experiment, followed by slow variations in the average confining pressure over the course of the experiments. These slow pressure changes are not significant at the timescale of the imposed oscillations. 3. The period of oscillation is limited by a lower bound due to the presence of an auto-compensation chamber, which precludes the investigation of very short periods. An upper bound to the oscillation period also exists with our experimental procedure, as we wish to observe multiple oscillation periods within each velocity segment. The oscillation period must therefore be shorter than the duration of the shortest velocity segment.
Finally, the apparatus used to impose the stress oscillations does not produce oscillations of exactly the specified period and can deviate from the imposed value by up to 2%. To improve the determination of the oscillation period, a Fourier Transform is applied to the confining pressure measurements. A Gaussian fit is then applied to the power spectrum in a frequency range surrounding the imposed period to determine the true oscillation period with improved precision (Gasior & Gonzalez, 2004). When the confining pressure measurements display discontinuities due to either erroneous values that require removing ( Figure S4 in Supporting Information S1) or software restarts, the overall period is determined by taking an average of the thus derived periods weighted by the length of each segment. Calculations of the period over these segments give period estimates that are in good agreement with one another.

Rate and State Parameters Calculation
As velocity changes are applied during our experimental procedure, the constitutive rate and state parameters of our sample can be retrieved. Indeed, according to the rate and state friction model (Dieterich, 1981), the friction coefficient is defined as follows: is dependent on both the rate of displacement and a state variable θ, such that during a velocity change from V o to V the friction coefficient, is defined as follows: where μ o is the initial friction coefficient measured at V o , A, and B are constitutive parameters and D c is the macroscopic critical slip distance over which the friction coefficient reaches its new steady-state value. The state variable θ is commonly defined with two empirical formulations (Ampuero & Rubin, 2008;Ruina, 1983), known as the slip law:̇= and the aging law:̇= 1 − The rate and state parameters A, B, and D c are determined for the control experiment using the RSFIT3000 Matlab GUI (Skarbek & Savage, 2019).

AEs Detection and Processing
AEs are detected in triggering mode, that is, by applying a voltage-threshold trigger logic, requiring at least 3 of the 8 acoustic sensors to reach a specified voltage (typically 100 mV, once amplified at 45 dB -x150) within a given time window (typically 50 µs) for waveforms to be recorded. Waveforms of 8,192 data points are recorded at 14 bit, with a 10 MHz sampling rate using an eight channel digital oscilloscope. Recording, storage, management, and basic data processing are performed using a licensed software (Insite, Applied Seismology Consulting Ltd.). The threshold trigger criteria are set through trial and error to be just above the background noise level of the acoustic monitoring. The threshold trigger criteria are also set so that the AE productivity rate remains below the maximum triggering rate capability (∼50 AE/s) of the hardware throughout the experiment, at least after the first velocity segment.
Waveforms are first filtered with a two-pass low-pass filter of cut-off frequency 2 MHz to remove noise and better detect AEs. An STA/LTA auto-picking procedure (Trnkoczy, 2009) is manually tuned to each sensor in each experiment after suppressing low signal-to-noise records to precisely determine the arrival time of the P waves resulting from the AEs. In our catalogs, timing of AEs correspond to the P-picks arrival times as the travel times within the sample (less than 10 µs) are negligible compared to the imposed stress oscillations periods. Using the average Root Mean Square (RMS) of the AE waveforms, relative magnitudes M r are calculated following Rivière et al. (2018) as follows: with n the number of working AE sensors. Absolute magnitudes can also be determined by fitting, after careful sensor calibration (Marty, 2020), the average displacement spectra for each AE using the Ω −2 law (Madariaga, 1976).

Coherent Averaging of the AE Temporal Distribution
To determine whether discrete AEs are correlated with the applied stress oscillations over the many periods present in each velocity segment, a coherent averaging distribution of AEs is calculated by considering where each AE lies relative to the phase of the stress oscillation: where ϕ is the phase at which the AE occurs within the oscillation, t is the time at which the AE occurs, t 0 is the reference time of the oscillations defined as the time of a stress oscillation maximum, and T is the stress oscillation period. This average distribution is a useful quantification of the average temporal distribution of AEs.

Magnitude Frequency Distribution
The magnitude distribution of AEs in laboratory experiments has repeatedly (e.g., Lockner, 1993;Marty, 2020;Meredith et al., 1990;Mogi, 1962;Scholz, 1968) been shown to follow the same distribution as earthquakes, that is, the Gutenberg-Richter (GR) law: where N is the number of events of magnitude greater than M that occur, a and b being constants. The b-value of the GR-law describes the ratio of occurrence of large to small magnitude emissions and corresponds to the slope of the magnitude distribution above a certain cut-off magnitude M c , called completeness magnitude, under which catalogs are viewed as incomplete. The b-value is calculated using both a least squaress method and a maximum likelihood estimate (Aki, 1965), the latter being considered more accurate. Indeed, larger magnitude ranges are given less weight, thus accounting in some measure for the uncertainty linked to the small sample size of larger magnitude emissions.

Mechanical Behavior
As axial loading is applied to the sample, the normal and shear stresses applied on the fault increase ( Figure 2). The sample accommodates a small amount of stress through elastic deformation before slip initiates. After the initiation of slip, the stresses continue increasing and tend toward an asymptotic steady-state behavior, corresponding to a friction coefficient of 0.7 (Figure 3a). This value is comparable to those generally observed in the literature for quartz-based gouges (e.g., Bedford et al., 2022;Leeman et al., 2016;Marone et al., 1990;Marone & Scholz, 1989;Samuelson et al., 2009).

Rate and State Parameters
Both the slip law and the aging law of the state variable are used for the calculation of the rate and state parameters, with very similar values being obtained with both laws. The fit of the friction coefficient evolution and rate and state parameters presented in Figures 3b and 3c correspond to the slip law. A − B decreases slightly during the experiment from an initial value of roughly 4e−3 down to a value of roughly 2.5e−3, and D c increases slightly from roughly 60 to 80 µm. Such changes in frictional properties with shear displacements are common in gouge experiments (e.g. Beeler et al., 1996;Dieterich, 1981).
Due to the increase in friction coefficient with slip, a detrending procedure is applied to the friction coefficient in order for the rate and state fitting to only consider the effect of the velocity change. This is done by fitting low order polynomials to the friction data and running RSFit3000 calculation on the fit residuals.
Due to our experimental geometry, the sliding surfaces are not orthogonal to the loading axes. As such, loading velocity changes induce both shear stress and normal stress changes, the latter inducing additional frictional changes due to the Linker-Dietrich effect (Linker & Dieterich, 1992). The rate and state parameters determined here therefore capture the effect of both the velocity change and the normal stress change. Still, the absence of macroscopic stick-slip in our experiments is the confirmation that the main conclusion of this calculation, that the gouge is velocity strengthening at the imposed velocities-that is, that A − B > 0-holds ground (Marone, 1998).

Figure 2.
Stresses evolution during the control experiment. As the piston must first overcome the static friction threshold to advance, the differential stress first increases without any recorded piston displacement. Slip is therefore overcorrected before the onset of slip, with negative slip values appearing at the onset of the experiment. The absence of stick-slip is convenient for the study of the temporal distribution of AEs. Indeed, the clustering of the AEs catalog into simple foreshock-mainshock-aftershock series is much reduced, and allows the presence of many smaller acoustic events, which lead to larger catalogs more suitable for statistical analysis. Whilst the gouge itself also undergoes strain, quantifying it precisely is not possible within our experiment. The gouge is therefore assumed to not undergo major changes in thickness after an initial compaction regime visible in the first few velocity segments of our experiment (Figure 3a), with both the friction coefficient and the AE rate reaching a plateau.

AE Rate and Mechanical Properties
All stress modulation experiments demonstrate the same evolution in AE rate and friction as the control experiment ( Figure S6 in Supporting Information S1 and Figure 3), and exhibit similar numbers of AEs to the control experiment (Table 2). An initial compaction regime characterized by high AE rates and low friction coefficients progressively transitions into a stable regime with consistently low emission rates and high friction coefficients in the later velocity segments.
The imposed stress oscillations induce fluctuations in shear and normal stress (Figure 4 and Figure S5 in Supporting Information S1), AE-rate ( Figure 5), friction coefficient ( Figure S6 in Supporting Information S1), and slip velocity on the fault plane ( Figure S9 in Supporting Information S1). The amplitude of the friction coefficient oscillations remains approximately constant at all velocities within an experiment ( Figure S6 in Supporting Information S1). The amplitude of the shear stress and normal stress oscillations, however, are highly dependent on the slip velocity and on the imposed stress oscillation parameters ( Figure S5 in Supporting Information S1). Most notably, the shear stress oscillations are of the highest amplitude for the highest average slip velocities.

AE Temporal Distribution Relative to the Stress Oscillation
As illustrated in Figure 5, AEs are not distributed uniformly over the stress oscillation but occur preferentially around specific phases, except for the lowest stress oscillation amplitude. This non-uniformity of the AE distributions is herein referred to as the modulation of the AE distribution. The influence of the oscillation amplitude on the distribution is greatly noted, the influence of the average slip velocity is also apparent. The higher the stress oscillation amplitude and the higher the average slip velocity, the more AEs are modulated. The influence of the period seems less significant.

AE Magnitudes
Absolute magnitudes range roughly between −9 and −7 ( Figure 6) and do not vary significantly in distribution across experiments. Comparing absolute and relative magnitude estimates for each AE shows a good agreement between magnitudes sets, both scaling linearly at the first order (Figure 6a), thus verifying the validity and reliability of the absolute magnitudes. This is important as absolute magnitude could not be computed for all AEs due to numerical difficulties when fitting the Ω −2 law reliably on waveforms with poor signal-to-noise ratio. In the following, magnitude statistics are therefore computed using the RMS method (described in Section 2.6), and assigning to relative magnitudes the corresponding absolute magnitude obtained by calibration and the Ω −2 law fitting.
The b-values calculated in our experiments are relatively high (Figure 6b), with little variation being noted from one experiment to the other.
The maximum likelihood estimations consistently yield smaller values between 2.8 and 3.3 whilst the least squares fits range between 3 and 4.5. This larger spread of b-values produced by the least squares method is expected due to its reliance on large emissions, of which the amount is highly variable due to the stochastic nature of AEs. Though rarely found associated with natural background seismicity of which b-values typically fall around 1, high b-values have been observed in seismic swarms (Adhikari et al., 2021) as well as in the presence of fluids (Bachmann et al., 2012;Murru et al., 1999). Such high b-values estimates could also in part be due to the relatively small dynamic range between the completeness magnitude and the maximum magnitude, which can cause overestimation of b-values, as discussed by Geffers et al. (2022)  4. Discussion

Coulomb Failure Model
When considering the effect of periodic stress on seismicity modulation and triggering, Coulomb failure stress is often invoked (e.g., C. W. Johnson et al., 2017;Scholz et al., 2019). Assuming a stable mechanical regime and according to Beeler and Lockner (2003) and Knopoff (1964), at oscillation periods larger than the nucleation time of microseismicity, a nucleation model based on Mohr-Coulomb failure should explain the distribution of the nucleation of microseismicity. In this theoretical framework, an AE occurs whenever the nucleation threshold is locally reached within the stress field. At oscillation periods smaller than the nucleation time, however, microseismicity should show reduced correlation with the oscillations.
Let us here consider the case of the Coulomb failure model, where strain is dissipated solely by nucleation of AEs. Due to the stochastic nature of AE nucleation, to the local unloading of stress by AEs, and to local geometric heterogeneities within the gouge, the stress field in the gouge is not homogeneous. We here assume a randomly . Stresses evolution during oscillation experiments. P c oscillations impose oscillation on the differential, normal, and shear stresses. Holes in the time series correspond to periods of either erroneous or missing data ( Figure S4 in Supporting Information S1).
distributed set of starting stress distributed throughout the sample. The constant background loading rate imposed on our sample is considered equivalent to a constant uniform stressing rate, and harmonic stress oscillations in turn induce harmonic oscillations in the stressing rate.
Within this framework, the theoretical AE distribution can be easily derived (see Figure 7), and is controlled by a single dimensionless parameter :̃= where ΔCFS is the amplitude of the Coulomb stress oscillation, ̇0 is the background Coulomb failure stressing rate, and T is the oscillation period. Assuming a constant stiffness of the apparatus and sample, the background stress rate is proportional to the background loading rate V. Thus, for simplicity, we introduce the following parameter c: such that ∝̃ . The parameter c is hereafter called the Coulomb stiffness because of its unit in MPa/m.
As we apply confining stress oscillations and axial loading, depending on the ratio of these loadings, our sample can experience both positive and negative Coulomb stressing rates. For positive Coulomb stressing rates, the probability density of the AE distribution is proportional to the instantaneous Coulomb stressing rate. If the Coulomb stressing rate is only positive, the average theoretical distribution of AEs is sinusoidal (Figures 7b-7e). For negative Coulomb stressing rates, also known as stress reversal, the phases of the oscillations during which the stress reversal occurs exhibit quiescence, that is, an absence of AEs, as the nucleation stress threshold is not reached (Figure 7f). In the case of stress reversal, the probability density of the AE distribution is equal to zero until the Coulomb stress reaches its previous maximum value.
The response amplitude s of the AE distribution directly correlates with the Coulomb stiffness, such that first order theoretical properties can be derived: 1. For a given value of the Coulomb stiffness, the response amplitude should be constant. 2. The response amplitude increases with the Coulomb stiffness. 3. Quiescence should be observed whenever stress reversal is present.
Due to the complex geometry of our experimental system and to the oscillations of both normal and shear stresses on the fault, the notion of stress reversal relates to oscillations of the Coulomb failure stress, defined as follows: where μ′ is the static friction coefficient. Assuming μ′ = 0.7, a value close to that reached by the effective friction coefficient by the end of each of our experiments and comparable to that reported by the literature (e.g., Scholz et al., 2019), the Coulomb stress oscillates and displays stress reversal during almost every velocity segment ( Figure S8 in Supporting Information S1). This is in line with tidal stress oscillations, which very commonly induce stress reversal (Beeler & Lockner, 2003;Heaton, 1982). An important assumption made here is that the shear and normal stresses at a macroscopic scale apply at a microscopic scale. Indeed, the stress field within the gouge is heterogeneous at the grain scale and highly dependent on the geometry of gouge grains and grain boundaries. However, we are unable to measure stress at this scale and must therefore rely on macroscopic stress, which gives us an average value of the stress throughout the sample.
The amplitude of the Coulomb stress oscillations is proportional to the amplitude of the imposed confining pressure oscillations, such that ΔCFS = 0.7ΔP c as observed in Figure 5. Therefore, the Coulomb stiffness can be expressed as follows:  Beeler & Lockner, 2003;Heki, 2003 (Beeler et al., 2000) as follows: Thus, considering a static friction coefficient of 0.7, the Coulomb stiffness can be generalized as follows: with Δ̃ the imposed stress oscillation amplitude (either ΔP c or ΔP f ).
Sinusoidal least squares fits of fixed periods of the AE phase histograms are calculated, the amplitudes of which are used as proxies for the response amplitude of the AE distribution. Interestingly, our experimental results ( Figure 8) do not match the theoretical properties of the Coulomb Failure model. The response amplitude is not directly correlated to the Coulomb stiffness, the influence of the stress oscillation amplitude being much greater than the influence of the oscillation period for instance, with increases in oscillation amplitude resulting in large increases in response amplitude. Moreover, the influence of velocity is the opposite of the model prediction, as increasing velocity leads to increases in response amplitude. Finally, no consistent period of quiescence is observed in our experiments, despite most of our experimental velocity segments containing stress reversal.
Moreover, the velocity segments where stress reversal occurs the most correspond to the lowest velocity segments, which also happen to be the velocity segments where the response amplitude is at its lowest within each experiment. This implies that at lower velocities, the AE rate becomes oscillation-independent. This could be due to the complex geometry of the stress field within the gouge. Stress oscillations may induce local stress heterogeneities, which may promote nucleation during macroscopic stress reversals. At larger periods and slower velocities, greater viscoelastic relaxation of stress through creep may also take place. Such creep could in turn lead to increased background AE-rates (according to the "preslip" model of earthquake nucleation; Beroza & Ellsworth, 1996) that would also work against any shear stress build-up. This could also be due to the lower phase clustering of AEs at low velocities, which could allow less shear stress build-up throughout the oscillation (Figure 9). These observations could also signify that our experiments are mostly conducted in the high-frequency regime described earlier.
This hypothesis is also supported by the phase of AEs. In most experimental segments, the preferential occurrence phase of AE correlates positively with the phase of greatest Coulomb stress. This is not in line with the Mohr-Coulomb failure model, which predicts that AEs should nucleate preferentially at phases of greatest stressing rate. The only segments where the correlation with the greatest Coulomb stress occurs are the highest velocity segments in the experiments with either the largest oscillation amplitude or the largest oscillation period. This is in line with the theory of two nucleation-time-dependent regimes, as initially described by Dieterich (1987). Both experimental (Beeler & Lockner, 2003) and numerical (Dublanchet, 2022) observations reproduce these behaviors. Likewise, our experiments are done around the transition between the two nucleation regimes.

Susceptibility of the AE Distribution Response Amplitude to Stress Oscillation Amplitude
For moderate stress oscillation amplitudes, numerical Burridge-Knopoff simulations (Pétrélis et al., 2021) observe a linear relation between the response amplitude of the AE distribution and the oscillation amplitude such that: Calculating the susceptibility of the response amplitude to the confining pressure amplitude (Δ̃) , that is, the response amplitude divided by the stress oscillation amplitude: a first-order relationship relative to the Coulomb stiffness is made apparent (Figure 10), confirming as a first-order approximation the linear scaling of the response amplitude with the oscillation amplitude.
Furthermore, data from studies involving competition between fluid-induced pore pressure oscillations and axial loading with different experimental conditions and setups to this study's (Chanard et al., 2019;Noël, Passelègue, et al., 2019;Noël, Pimienta, & Violay, 2019) follow the same trend when the susceptibility of their response amplitudes to the pore pressure oscillation amplitude is plotted against Coulomb stiffness. This result is obtained despite other experiments being conducted with fluid pressure oscillations and using rock saw cuts as slip interfaces rather than gouge. The response amplitude for these experiments is derived following the same procedure as for this study's experiments.

Rate and State Friction Failure Model
The Coulomb Failure model fails to explain the distribution of AEs, let us now focus on the rate and state failure model. In this high-frequency regime, rate and state friction predict the seismicity rate to be independent of period and loading rate (Dieterich, 1994(Dieterich, , 2007Heimisson & Avouac, 2020). The rate and state model is most often times discussed in terms of seismicity rate R. Seismicity rate variations are similar to the response amplitude s used to quantify our experimental results, and as such, we consider that changes to the seismicity rate translate proportionally to changes in s.
In the rate and state friction failure model, the response in seismicity rate on how the oscillation period compares to the characteristic timescale t a of response of the seismicity rate to stress perturbations with = ∕̇0 where ̇0 is the background Coulomb loading rate. As our experiments change the loading rate of our system, t a changes accordingly and varies between roughly 20 s during the fastest velocity segments and 5,000 s during the slowest velocity segments.
If T < t a then the seismicity rate response to small normal changes compared to the normal stress is: with r the background seismicity rate when no stress oscillations are present and M is a scaling factor that does not change the amplitude of the AE distribution response. This approximation is however, only valid for t ≫ t a (Heimisson & Avouac, 2020;Helmstetter & Shaw, 2009). As the mean AE rate is relatively constant in each of our velocity segments and as our velocity segments last for a duration far greater than t a , we consider the conditions under which this expression of the seismicity rate is derived to be met when the oscillation period is shorter than t a . In this case, following Heimisson and Avouac (2020), the AE response to the Coulomb stress oscillations can be expressed as follows: Considering our highest velocity segments where t a > T, the susceptibility to oscillation amplitude of the AE distribution response is linear, as established in Section 4.2. This linearity is compatible with Equation 20 if ΔCFS ≪ Aσ N . This condition, however, does not seem to be met, as we estimate Aσ N ∼ 0.2 MPa, with ΔCFS ranging from around 0.05 to 2 MPa. As such and considering that ΔCFS spans 2 orders of magnitude, we would expect an exponential response of the AE distribution to the oscillation amplitude, which we do not observe. This model also does not explain the dependence on loading velocity or oscillation period of susceptibility to the oscillation amplitude of the AE distribution.
The rate and state friction failure model does, however, describe that at T > t a the AE response correlates with the stressing rate and that at T < t a the AE response correlates with the stress. This is in line with Figure 5 where a correlation between the stressing rate and the AE distribution response is only observable at the highest loading rates and either highest period or oscillation loading rate where t a is at its smallest.
The case of velocity strengthening stress and velocity oscillations was investigated by Ader et al. (2012) considering small Coulomb stress oscillation amplitude relative to (a − b)σ N and small velocity oscillation relative to the background loading rate. Neither the Coulomb stress oscillations nor the velocity oscillations present in our experiments correspond to the conditions considered by this study. Furthermore, as the velocity oscillation amplitude is larger than the background velocity, the slip velocity becomes negative during our stress oscillations. To the best of our knowledge, this is fundamentally incompatible with the rate and state framework and prohibits any solid investigation of the link between shear stress oscillations and velocity oscillation.

b-Value Modulation
Stress oscillations have been shown to have an influence on the magnitude-frequency distribution of natural seismicity (Ide et al., 2016;Scholz, 2015;Tan et al., 2019). An analysis of the variations of b-value throughout oscillations was thus undertaken using our experimental catalogs. All AEs except those contained in the first 1 mm of displacement of each experiment were here considered in order to have large enough catalogs for maximum likelihood b-value estimates to be statistically relevant.
AE sub-catalogs were defined for each experiment by selecting events within a given phase window. This phase window was then moved across the phase space to investigate the variations in b-value. The choice of the size of the phase windows induces a bias in the estimation of the b-values and the associated uncertainties. Larger phase windows contain a larger number of AEs and therefore are associated with smaller uncertainties, but variations of the b-value are smoothed out, whereas smaller phase windows contain fewer AEs and are therefore associated to larger uncertainties but do not suffer from as much smoothing of the variations of the b-value. A compromise thus had to be made when selecting a phase window size to highlight variations of b-values that are significant when compared with the associated uncertainties. The width of the phase window was here set to both π and π/2 to better evaluate the significance of the variations in b-value ( Figure 11).
Oscillations in b-value are observable regardless of the imposed phase window width for experiments 1 (ΔP c = 0.5 MPa, T = 100 s) and 3 (ΔP c = 2.5 MPa, T = 100 s; see Figure S10 in Supporting Information S1). The amplitude of oscillation of the b-value increases with the stress oscillation amplitude, as observed in natural settings (e.g., Scholz, 2015;Tan et al., 2019). However, the different b-value oscillations are not in phase.
According to numerical studies conducted on Burridge-Knopoff and Olami-Feder-Christensen models, b-value oscillations are expected when the modulation frequency is small enough. These b-value oscillation are in phase opposition with the normal stress oscillations (Pétrélis et al., 2021). However, our experimental results do not conform to these numerical results and would require further investigation. The effect of the oscillation period on the b-value modulation is less clear. Indeed, for experiments conducted at ΔP c = 0.5 MPa, the error bars attached to the b-value oscillations are comparable in amplitude to the observable oscillations. Therefore, it is not possible to determine the dependence of the b-value oscillation amplitude to the oscillation period.

Relevance to Natural Oscillatory Stress Phenomena
Our experimental observations show an increase in microseismicity modulation primarily with increased stress oscillation amplitude as well as with both oscillation period and background loading rate. The increase in microseismicity modulation with stressing amplitude is in line with the natural case, where seismicity has robustly been shown to be triggered by large stress variations caused by seismic wave propagation (Brodsky & van der Elst, 2014), and more modestly modulated at longer periods by smaller amplitude of stress variations (e.g., Bollinger et al., 2007;Métivier et al., 2009). Moreover, the increase in modulation with loading period could explain why despite tidal and seasonal or multiannual stressing being of comparable amplitudes, tidal modulation is less commonly observed than seasonal and multiannual modulation (Beeler & Lockner, 2003). Interestingly, our experimental results also show the role of background loading rates in modulating seismicity, which seems consistent with larger-scale observations. Indeed, natural faults appear to express different levels of sensitivity to oscillating stresses, for example, at seasonal periods, with the best constrained modulations being witnessed in tectonic contexts with large strain rates such as the Himalayas (Bollinger et al., 2007) or Alaska (C. W. Johnson et al., 2020). In fact, similarly to experimental results, natural observations of triggered or modulated seismicity by various transient and oscillatory phenomena can be represented by the Coulomb stiffness as a function of the oscillation period, over more than 10 orders of magnitude ( Figure 12). Within this framework, natural observations and experimental results of this study roughly fall in a range of constant Δ̃∕ . As such, the susceptibility of earthquake nucleation to stress oscillations in natural setting seems better probed by experiments applying Figure 11. Gutenberg-Richter b-value calculated for acoustic emissions (AEs) within a moving phase window of π (red) or π/2 (blue) throughout Experiment n°3 (ΔP c = 2.5 MPa, T = 100 s). The b-value is determined with a maximum likelihood estimate.
small amplitude of stress oscillations. However, comparing experiments conducted by different studies is difficult due to different sample geometries, experimental procedures, and sample material. For instance, experiments conducted with intact rocks were conducted at higher Δ̃∕ ratios than experiments using either faulted or gouge samples and higher ∕ ratios than are usually observed in the natural setting but still present important results that could be applicable to systems that undergo small tectonic loading or particularly large stress oscillations.
The scaling of these observations is to be considered with caution. Indeed, the largest tectonic loading rates correspond to highly seismogenic zones (Ide, 2013) which deliver larger seismic catalogs and thereby make them more studied. Additionally, the seasonal modulation of seismicity in the New Madrid Seismic Zone (Craig et al., 2017) remains a troubling observation, as regional Global Navigation Satellite System (GNSS) observations do not indicate any large-scale deformation (Craig & Calais, 2014), although the seismicity could be explained by the presence of localized stress on and around specific faults not captured by the sparse GNSS network.
Experiments also indicate that an increase in the amplitude of GR b-value oscillation is correlated with an increase in stress oscillation amplitude, but experimental limitations preclude the quantification of a phase link between stress and b-value. Yet, this link is coherent with the few natural observations of b-value modulation by stress oscillations (Scholz, 2015;Tan et al., 2019). Given that larger magnitude events have larger nucleation times (Ohnaka, 2000), it would be expected that the response regime of earthquakes to stress oscillations would be magnitude-dependent. This could explain the b-value variations with stress oscillations and is in line with observations in the natural case where differential stress variations are linked to the b-value.
Some limitations of our experiments include assuming a fault geometry which is akin to that of thrust faults, these being the faults on which the modulation of earthquakes by periodic loading seems to be the most readily observed (Cochran et al., 2004;Ide et al., 2016;C. W. Johnson et al., 2017). It would be of interest to also consider the case of strike-slip geometries in the laboratory with an appropriate experimental device, to explore their response to different stressing conditions. Another experimental limitation of this work is the use of oscillation periods corresponding to the range between dynamic triggering and tidal characteristic Coulomb stiffness. The large oscillation periods of most natural oscillatory stress phenomena make experimental work with those periods difficult. Figure 12. Coulomb stiffness as a function of oscillation period for experimental data and select natural seismicity observations from earthquake catalogs. In natural contexts with similar background stressing rates, stress oscillation of similar amplitudes, but different periods elicit different responses from seismogenic zones, for example, directly or delayed triggered seismicity, modulation of seismicity rate, or no observed correlation. Velocities are taken from the references shown in the figure. Periods and stress variations are given in Table S1 in Supporting Information S1.

Conclusions
Confining pressure oscillations experiments-of amplitudes ranging from 0.1 to 2.5 MPa and periods ranging from 20 to 500 s-are conducted on gouge layers with a stress geometry associated with reverse faulting and undergoing slow axial stressing-inducing slip ranging from 1 μm/s to 10 nm/s-to determine the modulation of microseismicity by stress oscillations in tectonic contexts. The resulting response amplitude of the microseismicity distribution is linearly susceptible to oscillation amplitude and increases with velocity. Our results mostly do not conform to the low-frequency Coulomb-failure regime, indicating that the considered oscillation periods are smaller than the nucleation times of microseismicity (here recorded in the form of AEs). The high-frequency rate and state regime also fails in explaining the linear susceptibility of the AE distribution response to stress oscillation amplitude despite the range of parameters probed. However, the susceptibility to stress amplitude does scale with the Coulomb stiffness. Moreover, we show that experimental results tend to follow the same Coulomb stiffness-period scaling law as natural seismicity linked to oscillatory phenomena. We experimentally confirm the absence of quiescence during stress reversal, probably linked to complex stress heterogeneities and viscoelastic relaxations taking place in our granular fault gouge medium. We also confirm the possibility of small oscillatory perturbations modulating the b-value, which is relevant to observations of tide-based modulation of seismicity b-value (Ide et al., 2016). Further investigation of the influence of stress oscillations on b-value could be useful for earthquake hazard assessment, as well as more experimental work to assess the response of AEs in both the high-frequency nucleation-driven regime and low-frequency threshold-driven regime. Finally, one important perspective of this work is to address the specific case of the response and susceptibility of a system undergoing stick-slip motion, for which the microseismicity might be clustered into foreshock-mainshock-aftershock sequences, in order to investigate further how the susceptibility evolves during the seismic cycle at the laboratory scale.

Data Availability Statement
The acoustic emission waveforms, mechanical data and derived data used to construct the figures in the study are available at Zenodo.org via https://doi.org/10.5281/zenodo.7786438 with a Creative Commons Attribution 4.0 International Licence.