The Spectrum of Slip Behaviors of a Granular Fault Gouge Analogue Governed by Rate and State Friction

Currently, it is unknown how seismic and aseismic slip influences the recurrence and magnitude of earthquakes. Modern seismic hazard assessment is therefore based on statistics combined with numerical simulations of fault slip and stress transfer. To improve the underlying statistical models we conduct low velocity shear experiments with glass micro‐beads as fault gouge analogue at confining stresses of 5–20 kPa. As a result, we show that characteristic slip events emerge, ranging from fast and large slip to small scale oscillating creep and stable sliding. In particular, we observe small scale slip events that occur immediately before large scale slip events for a specific set of experiments. Similar to natural faults we find a separation of scales by several orders of magnitude for slow events and fast events. Enhanced creep and transient dilatational events pinpoint that the granular analogue is close to failure. From slide‐hold‐slide tests, we find that the rate‐and‐state properties are in the same range as estimates for natural faults and fault rocks. The fault shows velocity weakening characteristics with a reduction of frictional strength between 0.8% and 1.3% per e‐fold increase in sliding velocity. Furthermore, the slip modes that are observed in the normal shear experiments are in good agreement with analytical solutions. Our findings highlight the influence of micromechanical processes on macroscopic fault behavior. The comprehensive data set associated with this study can act as a benchmark for numerical simulations and improve the understanding of observations of natural faults.


Fault Slip
Active faults are characterized by a wide range of slip behaviors ranging from aseismic creep to seismic stick-slip that may change spatially along the fault and temporally over the seismic cycle (e.g., Harris, 2017;Peng & Gomberg, 2010). The types of slip are defined by their characteristic timescale which ranges from milliseconds to a few years (Obara & Kato, 2016) and by their characteristic magnitude which is usually defined by seismic moment (Gomberg et al., 2016;Ide et al., 2007). Depending on their characteristics in time and seismic wave forms, the slip events are characterized as seismic (very low frequency earthquakes, tremors, normal earthquake) or geodetic (short-term and long-term slow slip events) events. They can occur simultaneously, that is, within one seismic cycle, at the same locality or in different depth ranges of the same main fault (Bürgmann, 2018). The physical origin of this range of slip modes is still not entirely clear, although several approaches for certain phenomena have been proposed (Chen & Spiers, 2016;Ciamarra et al., 2010;Daniels & Hayman, 2008;Dorostkar & Carmeliet, 2018).
A common methodology to model this wide range of slip behaviors is through a continuum based description that reproduces the kinematics and dynamics of fault activity. The rate-and-state framework provides the possibility to characterize fault behavior, or in a general term "fault rheology," by describing the connection of forces in the system (friction μ) and the external influences such as loading rate v L and stiffness k (Brace & Byerlee, 1966;Dieterich, 1978Dieterich, , 1979aDieterich, , 1979bScholz, 1998). In general, the rate-and-state framework is able to describe most observations that lead to fault (in-)stability and has been derived from experimental observations in the laboratory and a few field observations (Marone, 1998, and references therein). Stick-slip experiments using rock and rock analogs suggest that besides intrinsic material properties (e.g., friction coefficient, slip/velocity weakening), extrinsic parameters like stiffness, normalized loading rate and effective normal stress are key controls of frictional stability (e.g., Heslot et al., 1994;Leeman et al., 2016;Mair et al., 2002;Marone, 1998). Recent studies also highlight that several of the fault intrinsic parameters in the rate-and-state equation are also dependent on extrinsic parameters and not constants as previously assumed (Chen & Spiers, 2016;Van den Ende et al., 2018).

Granular Fault Analogs
In this study, we purely focus on the frictional characteristics of an analogue fault zone which is described with the rate-and-state framework (Dieterich, 1979a(Dieterich, , 2007. Our fault zone is composed of a granular fault core with relatively stiff outer boundaries and dominated by granular mechanics. Other processes that influence the slip modes along a fault zone, which are not realized in our setup, are variations in pore-fluid pressure, changes in material because of comminution, or mineral reactions. Not all slip modes are observed for all active zones which strongly suggests that there is a complex interaction between the processes acting on different scales in space and time. However, there is growing evidence that the various slip modes, for example, slow slip events, are ubiquitous and part of a continuous spectrum of transient deformation (Jolivet & Frank, 2020). Therefore, the knowledge of the complex interactions between the different slip modes is highly relevant for estimating the seismicity rates along plate boundaries and therefore for seismic hazard assessment. Other possible areas of application include soil mechanics and mass movements.
Here we report characteristics of slip events in an analogue fault gouge consisting of spherical soda lime glass beads with a diameter of 300-400 μm. In contrast to similar experiments of Frye and Marone (2002), Anthony and Marone (2005), Ferdowsi et al. (2013), Jiang et al. (2016), and Cui et al. (2017) we explore the low pressure (kPa instead of MPa) regime which is rich in slip behaviors and generates regular stick-slip with more complete stress drops similar to seismic cycles along major faults in a highly reproducible and accessible way. The stress drops in our experiments can reach up to 20% of the mean stress, which can be higher than typical inter-plate and intra-plate events (Scholz, 2019). However, higher stress drops can be achieved when additional mechanisms are lowering the energy dissipation along the frictional interface, for example, fluid-rock interaction (Sibson, 1980). Assuming weak subduction zone interfaces stress drops can be as high as the presumed strength of the interface, leading to complete release of elastic stress during large earthquakes. Several studies established the large diversity in slip modes in such experiments. Changes in stiffness and normal stresses lead to first order changes in frictional stress, such as transition from stick-slip to oscillation and stable sliding (Heslot et al., 1994). Nasuno et al. (1997) found localized precursor phenomena in thin sheared glass beads that precede large slip events. Moreover, the use of a ring-shear tester instead of commonly used direct shear apparatuses allows us to apply an in principle infinite amount of displacement and therefore a large number of events, which leads to a solid database for statistical analysis. Results from a similar apparatus by Cain et al. (2001) show that it is suitable to 10.1029/2021GC009825 3 of 28 measure dilation-compaction cycles and show that the conditions in an annular shear cell lead to dilation during the loading phase and compaction during failure which is similar to the results obtained from rock mechanical tests in a rotary shear apparatus for rocks (Beeler & Tullis, 1997). Previous studies have shown that dilatancy may play an important role for the stick slip characteristics of granular media (Lacombe et al., 2000) and for nucleation of slow slip events (Segall et al., 2010).
For the same material we vary the extrinsic parameters normal stress σ N , loading velocity v L , and stiffness k L . In this parameter space, we monitor the occurrence of slip events and creep, as well as the transitions from one slip mode to another. We call these tests the stick slip (STSL) tests. Furthermore, we characterize the analogue fault gouge with commonly used tests to derive the rate-and-state parameters, such as slide-hold-slide tests (SHS). We compare the findings to first order observations from rock friction experiments and assess the suitability of granular analogue fault gouge for its use in combined analogue and numerical modeling.

Methodology
To simulate fault behavior in various settings we use a granular analogue modeling approach. Previous studies examined granular media under natural pressure conditions, whereas we are using conditions realized by analogue models, being 3-4 orders of magnitude lower (Rosenau et al., 2017). This prevents comminution of the glass beads and ensures constant frictional properties over the test duration, which gives well reproducible results. The terminology used to describe and characterize the data is found in Table 1 and Appendix A. The data analysis is done using a suite of Python scripts that pick events and do statistic calculations. All of which are available as the open source software "RST-Stick-Slipy" from the GFZ git repository (Rudolf, 2021) and are also included in the data publication .

Rate-and-State Friction
The relation between shear stress and normal stress for granular media and many other interfaces is determined by a non-linear combination of mean stress, slip velocity, stiffness, and several non-dimensional parameters. This relationship is termed rate-and-state dependent friction that macroscopically leads to alternating cycles of slip, creep, and locking, called stick-slip (Beeler et al., 1994;Dieterich, 1978Dieterich, , 1979aMarone, 1998;Ruina, 1983;Tullis & Weeks, 1986). This effect is used to describe and explain the various slip behaviors that are associated to earthquakes, for example, slip on faults (Marone, 1998), earthquake nucleation (Dieterich, 1992) and slow slip events (Leeman et al., 2016). In our study we use the relationships and testing procedures defined in Beeler et al. (2001), Marone and Saffer (2015), and Dieterich (2007) as well as adapted methodologies of Corbi et al. (2013), Bhattacharya et al. (2015), and Bhattacharya et al. (2017) to estimate the principal parameters for the rate and state equation. A short description of rate-and-state friction and its application to our study is found in Appendix A2.

Experimental Setup
For the experiments we use a ring shear tester of type "RST-01.pc" (ASTM, 2016;Schulze, 1994) which allows to shear a granular sample in an annular shear cell. The machine and methodology has been verified and calibrated using a standard bulk material (CRM-116 limestone powder) and is extensively used for characterizing granular materials in engineering and analogue modeling (e.g., Klinkmüller et al., 2016;Ritter et al., 2016a;Schulze, 1994).
The granular material is confined in a ring shaped shear cell and sheared against a lamellae-casted lid which also imposes the normal load ( Figure 1). The normal load is adjusted using a motorized weight attached to a lever that pulls the lid from below. This ensures a constant normal load on the sample. Two bars attached to force transducers hold the lid in place and measure the shear forces acting on the lid.
The stiffness of the testing apparatus has important implications on the criticality of the system (Leeman et al., 2015). We measure it by fixing the lid to the shear cell and measuring the force increase while moving the shear cell ( Figure A3). The basic stiffness of the apparatus = 246.0 (RST) is mainly influenced by the stiffness of the load cells that measure shear force which acts in series with the aluminum of the lid, tie rods 10.1029/2021GC009825 4 of 28 and crossbar. Adding springs in between the load cells and tie rods lowers this stiffness to the values reported in

Fast event
Abrupt reduction in shear stress along the shear zone coinciding with a counter rotation of the lid. Has a start "Event start" and end "Event end" defined by characteristic points on the shear curve Slow event Similar to a "Slip event" but with intensity and slip velocity a few orders of magnitude lower Recurrence time t r Time between the end (t e ) of a "Slip event" and the start of the next (t p ). The time span is named interevent phase analogous to the interseismic period for earthquakes Event peak p Maximum shear stress before a "Slip event" Event end e Minimum shear stress after a "Slip event" Onset of dynamic event d Critical point where slip velocity during an event is larger than the loading rate End of dynamic event f Critical point where slip velocity during an event is lower than the loading rate Preslip d p Slip that happens during the acceleration of a slip event between "Event peak" and "Onset of dynamic event" Creep Ratio of "Slip displacement" and "Load point displacement" during the interseismic phase. Due to permanent internal deformation at very low rates, there is a deficit between displacement that is imposed on the sample and released slip during a "Slip event" Onset of creep c Position on the shear stress curve where the reloading deviates from the linear trend (defined as "cyclic reloading stiffness k L ") by more than 1% Table 1 Terminology and Definition of Characteristic Points NI-9219, controlled by custom in-house software). This change was due to the end-of-life of the operating system during the course of this study which lead to hardware incompatibilities with the PCI-based approach. Similar to the other measurements, this data is averaged down to a frequency of 1 kHz.  Based on the setup geometry, shear and normal forces are converted into shear and normal stresses according to ASTM (2016) and lid displacement into volumetric change (dilation/compaction). The shear force F s is measured along a circle with a radius of r s that is spanned by the anchor points of the crossbar. This results in a torsional moment which acts within the sample cell: The moment is applied to the sample material by small lamellae that are attached to the lid of the shear cell. Because of the circular nature, displacement at the outer edges is higher than at the inner edges. Therefore, the shear forces are converted to shear stress using the moment of the crossbar M d (Equation 1), the median radius r m and the cross-sectional area of the lid A d : The median radius r m is also the reference position at which the loading rate v L is defined because it divides the cell into two regions of equal volume. While the normal stresses are also continuously measured, we assume a constant normal stress that is equal to the value set at the start of the measurement. The normal stress σ N is calculated by dividing the applied normal force F N with the lid area A d . Due to internal correction factors which are not disclosed by the manufacturer there is always a slight discrepancy between set normal stress and measured normal stress. Partially, this discrepancy stems from the angle of the tie rods and the lid which exerts additional normal stress onto the sample (pers. comm. D. Schulze).

Material Description
As granular material we use 300-400 μm sized fused soda-lime glass micro-beads supplied from Kuhmichel Abrasiv GmbH (Figure 1e). They are characterized by a relatively low dynamic friction coefficient (μ = 0.47, [Ritter et al., 2016a]) and no measurable cohesion (C = 1 ± 12 Pa, [Ritter et al., 2016a]) as well as a strain hardening-weakening behavior associated with dilation-compaction (Klinkmüller et al., 2016;Lohrmann et al., 2003;Ritter et al., 2016a). Glass beads are frequently used as a rock and gouge analogue material and generate stickslip under laboratory conditions (e.g., Mair et al., 2002). Because they are non-cohesive we can approximate the instantaneous frictional resistance of the fault zone μ as the ratio of shear stress τ to the applied normal stress σ N : To calculate dilation we first detrend the lid displacement d lid by a long term sliding average ̄ to remove the effect of transient loss of material. This effect accounts for roughly 1-2 mm of lid displacement over a complete experimental run which is subdivided into nine different velocities. Additionally, we normalize dilation by the average grain diameter of the glass beads d GB to get a better understanding of the spatial scale of dilation. The resulting quantity (Equation 4) is therefore dimensionless.

Experimental Procedure
Before a test is started, the sieved samples are presheared by 10 mm at a loading velocity of 0.5 which ensures a fully developed shear zone without major post failure weakening (derived from Ritter et al. [2016a] and Ritter et al. [2016b]). The recording of data is started after the preshearing phase and therefore the onset of stick slip motion is not recorded for our tests. Table 2 lists the experimental parameters for the various tests performed for this study. For all tests (STSL and SHS) we use four different normal stresses of 5, 10, 15, and 20 kPa. The major difference between the tests is the stiffness k M and loading rate v L .
The STSL tests are conducted with logarithmically spaced velocity v L from 0.02 to 0.0008 . The duration of each run is roughly 5 times the inverse of the respective loading velocity leading to equal displacement (≈5 mm) and a similar amount of slip events, for comparable conditions. Due to manually controlled test time, some time series might have more or less displacement. Each individual test is carried out at constant normal load.
To limit the influence of stick-slip we perform the stressed SHS tests with maximum machine stiffness and at higher shear rates. To also estimate the rate effect on healing rate we vary the loading rate v L from 0.05 to 0.52 . The hold times were increased logarithmically t hold = {1…1,000} s (in accordance with Equation A7) and a constant load point displacement of 5 mm between the hold phases was applied. Additionally, we did one series of unstressed SHS tests with the same set of parameters as for the stressed SHS test but only with a single load point velocity of = 0.16 . For unstressed SHS tests we reduce the normal stress to zero during hold to evaluate whether the change in state θ is purely time dependent (Aging law) or shows a slip dependent behavior (Slip law). Due to the time needed for the machine to unload and reload the normal stress, which was a few seconds per kPa, we only used hold intervals of t hold = {100…1,000} s. All SHS-cycles were repeated 3 times to be able to estimate the variance of the measurements.

Adjustment to Other Setups
To compare our experimental conditions with other setups in terms of stress conditions and stiffness we utilize the normalized stiffness k N and normalized loading rate v N . Both have an influence on the material behavior because of rate-and-state friction and are dependent on the setup. For our setup the normalized stiffness is calculated from the machine stiffness k M , which is modified using springs, the geometrical factor L M converting shear force to shear strain and the applied normal stress σ N .
From this the normalized loading rate v N is derived with the loading rate v L (Equation 6). It can be interpreted as a non-dimensional stressing rate that describes the increase in stress counteracted by friction over time.
= ⋅ For our setup the values for k N are in the range of 10 −2 to 10 2 mm −1 and for v N varies between 10 −5 and 10 −1 s −1 . This is in the same regime as Nasuno et al. (1997) and at least three orders of magnitude higher than the values achieved by Beeler et al. (1994). Further information and the raw data for the measured stiffness can be found in the data publication associated to this publication .

Results
Here we describe the slip modes qualitatively (Figures 2 and 3) and quantitatively using the asymmetry of the event cycles ( Figure A2). Then we determine the constitutive parameters for our setup and analogue fault gouge, which determine the stick-slip characteristics and slip behavior by slide-hold-slide tests. To compare the data across the individual setups we use the normalized loading rate v N (Equation 6) as a key parameter. This parameter contains the joint influence of normal stress and loading rate and makes it possible to plot all test results together without major overlap. In all cases we observe a distinct influence of both parameters on the individual results. Consequently, we describe the findings with respect to variations in stiffness and normal stress separately, where a clear distinction could be made. We use a similar color and marker code in most plots that show results from the tests (e.g., Figure 4). Normal stress, in some cases loading velocity, is indicated by color while the setup stiffness k M is indicated by markers. All errors in plots or in numeric values are given as twice the standard deviation (2σ) of the respective quantity. A rigorous error propagation is done during data analysis using the Python module "uncertainties" (Lebigot, 2021).

Slip Mode
During a test we monitor the change in apparent friction (Equation 3) over load point displacement ( Figure 2a). Slip events are detected as reduction in apparent friction over a short time frame (Figure 2b). Some show a smaller and slower reduction while other are faster. Additionally, the thickness of the material is recorded and converted to dilation (Equation 4, Figure 2c), which also shows an evolution that follows the evolution of friction to some extend.
The slip mode is qualitatively defined by the evolution of stress during an test (Figures 3a-3d). Low stiffness leads to typical sawtooth shaped curves with very sharp acceleration immediately before failure ( Figure 3b). Increasing the stiffness increases the amount of preslip and slows the acceleration before failure. This is expressed as slightly smoother sawtooth curves but the duration of a slip event is still much shorter than the reloading phase.
For Spring C we find oscillations of weakly irregular shape ( Figure 3c). On average the increasing edge of an oscillation is a bit longer but only by a factor of 1.5-2 and not several orders of magnitude as for the softer springs. Figure 3d shows stable slip behavior, which is only found for high slip rates and intermediate stiffness (Spring C). Another slip mode is observed for the highest stiffness which shows stick-slip cycles with a plateau of high stresses before failure. If the sample is at these high stresses we observe small and slower slip events that occur very close to failure ( Figure 3a). Figures 3e-3h shows an overview of all qualitatively determined slip modes in For all experiments we observe a characteristic succession of dilation and compaction. For perfect stick-slip, slip events lead to strong compaction of Δd ≈ 0.07 grain diameters d GB which then slowly dilates during the interevent period. In the first moments of failure for low stiffness experiments we often observe an initially dilating motion in the first few milliseconds. Experiments with stick-slip also show oscillations and characteristic patterns on a variable scale while experiments in the oscillating regime show a mostly stable pattern of ±0.001d GB . For low stiffness experiments this oscillation is on a scale of ±0.003d GB with a period of 2.5 s. Increasing stiffness leads to a reduction of the oscillation period to values of 0.2 s for Spring B and .1s for RST. With increasing amounts of creep we find additional oscillations in the dilation signal which gradually change their period closer to failure ( Figure 2c). For bimodal experiments, we observe a pattern of abrupt dilation events during the interevent period which is similar for consecutive interevent phases ( Figure 2c). The gradual change in oscillation period is also present for these experiments but on a smaller scale than for lower stiffness.
Quantitatively we describe the slip mode through the average asymmetry r a of the stick-slip curves and its distribution. Asymmetry r a is defined as the ratio of slip duration t s versus the reloading time t r ( Figure A2). We find that at generally low normalized loading rates (v N < 10 0 s −1 , low stiffness and low normal stress, Figures 3a, 3b, and A2a) the asymmetry is very high and has a low variability, although the data set is relatively sparse in that region due to very long reloading phases (t rel > 10 2 s, e.g., Figure A2a3). Pure stick-slip, dominant at low setup stiffnesses, has a very high asymmetry because the reloading phase is very long compared to the duration of an event (e.g., Figure A2a). In general both experiments with low stiffness (Spring A and B), show relatively regular and well defined stick-slip events. This is indicated by relatively flat increases in shear stress and abrupt decays with strong to modest acceleration before failure ( Figure 3a).
Experiments with Spring C show three different slip modes depending on normal stress and loading rate. With increasing loading rate we observe an evolution from bimodal over oscillation to stable. This evolution is clearest for low normal stresses and less apparent for high normal stress (Figure 3). At normalized loading rates v N < 10 0 the slip mode is bimodal with oscillating events preceding asymmetric events (Figures 3a-3c and e.g., Figure A2c1). In the interval v N = {10 −0.5 …10 0.5 } we find low asymmetry with low variability which approaches r a = 1, which is the expression of oscillating events becoming more and more symmetrical. The asymmetry decreases until we find oscillating slip modes at normalized loading rates between 1 and 10 s −1 . Oscillations are characterized by symmetrical increases and decreases of shear stress with an almost sinusoidal character (Figure 3b). In terms of asymmetry this leads to an average ratio r a ≈ 1 with a small variance ( Figure A2c1-3). At higher normalized loading rate the system becomes stable and the asymmetry shows a large variance with a mean asymmetry of r a ≈ 1. Slip under these conditions tends to be chaotic and does not show any characteristic features.
At highest stiffness = 246.0 (RST) we see a different evolution of slip mode with complex sets that are influenced by normal stress and loading rate. For lower rates v N < 10 2 the events split into two different distributions one with high asymmetry and one with lower asymmetry (Figure A2d2-4). As shown in Figure 3d, slow and small events alternate with larger events in a characteristic sequence. The duality of slip modes leads to a bimodal distribution of asymmetry with a distinct separation of small/slow and large/fast events. Normal stress has a strong influence on how well defined the separation between these two is. High normal stress leads to a clear separation which is at least one order of magnitude. At higher rates continuous distribution is observed while retaining a mean that is larger than 1, which still indicates defined stick-slip cycles rather than randomness.
In general, we find that above a certain stiffness the slip mode switches from simple stick slip (Spring A + B) to a more complex pattern of slip modes. Furthermore, high loading rates suppress the evolution of well defined stick slip cycles and lead to oscillations and stable slip modes. Normal stress defines the behavior for low loading rates, which is most apparent at higher stiffness (Spring C + RST). At low normal stress, the slip mode is mainly bimodal which is due to higher amounts of creep. Increased normal stress suppresses creeping mechanisms and forces a change from bimodal slip mode, with slow events to pure stick slip. Additionally, this leads to a shift in mode space so that oscillations occur at higher loading rates (Figures 3b-3d vs. Figure 3a).

Stiffness
In addition to the machine stiffness k M the reloading stiffness k L plays an important role during shear. It is a combination of machine stiffness k M and material stiffness. This property is calculated from the linear part of the reloading phase between slip events (Leeman et al., 2015). We find that for our experiments the influence of the glass beads leads to a reduction stiffness by roughly one order of magnitude (Figures 4a and 4b). Minor differences are observed for low normal stresses or when considering non-dynamic events. For Spring B there are a few outliers and a weak increase in k R is observed due to the influence of normal stress for the lowest normal stress σ N = 5 kPa. Spring C shows different results depending on the type of events considered. If only dynamic events with a slip rate above a critical threshold (comp. Section 3.1) are considered ( Figure 4a) the reloading stiffness k L is one order of magnitude smaller than the machine stiffness. Considering all events (Figure 4b) we see that the reloading stiffness is reduced by up to 2 orders of magnitude and there is a large spread of values probably due to the influence of slower events in bimodal slip mode (Section 3.1).

Event Magnitudes
Comparing the frictional stress drops Δμ for all experiments we find three different groups of stress drop highlighted in Figure 4c. These are distinguished by their magnitude and variability of stress drop, as well as the evolution with increasing normalized loading rate.
The first group occurs at low to medium normalized loading rate and at high stress drops (red area in Figure 4c). The stress drop shows an exponential decrease with a similar slope for most experiments in this region. It consists mainly of fast slip events and experiments at low to intermediate stiffness (Spring A, B, and C). A minor outlier is the Spring A-5 kPa experiment which has slightly higher stress drops but the same slope.
The second group consists exclusively of fast events at the highest stiffness (blue area in Figure 4c). They all plot at high normalized loading rate and show the highest stress drops which is roughly one order of magnitude higher than for the previous group at similar normalized loading rate. The evolution of stress drop shows a different slope, but is also decreasing with increasing normalized loading rate.
The third group are slip events with a low stress drop Δμ < 10 −2 and a large variability in stress drop that may span one or two orders of magnitude (green area in Figure 4c). This group is dominated by slow events. In general, the stress drop is decreasing for increasing normalized loading rate, but the slope is not constant.

Slip Velocities
From the above observations we find that there is a characteristic difference between fast and slow events for experiments with high stiffness. For highest stiffness we observe a clear separation of events into a small and fast cluster leading to a bimodal distribution of asymmetry ( Figure A2d1-4). This difference is highlighted using the average slip velocity during an event as an indicator (Figure 4d). The fastest slip events with ≈ 10 are observed for the lowest stiffness (Spring A, = 5.546 ) which is mainly due to the much larger slip during an event. With increasing stiffness the fast slip velocity is decreasing to a level of 10 −0.5 to 10 −1 . When slow events are present, which is not the case for all experiments at intermediate stiffness (Spring B + C), they are generally one to two orders of magnitude slower than the fast slip events ( ≈ { 10 −2 … 10 −3 } ) . The difference between fast and slow events increases toward the highest stiffness experiments, which also has the highest variability for slow slip velocities. Additionally, the peak slip velocity, that is the highest instantaneous slip velocity during a slip event, increases with increasing stiffness. The peak slip velocities are generally higher or in the same range as the mean slip velocity for fast events. Toward higher stiffness, the peak slip velocity is 2 orders of magnitude higher than the average slip velocity during a fast event. Typical peak slip velocities are in the range of = 1 … 10 .
These slow events are characterized by low stress drop, low stress drop rate and a characteristic occurrence late in the cycle at generally high mean stress ( Figure 5). The relative amount of slow events decreases with increasing normal stress. For low normal stress more than 40% of the total events are found to be slow events, whereas for higher normal stresses it is 5%-10%. Additionally, there is a variation in occurrence with loading velocity. At high loading velocity only very few slow events are detected, while at low loading velocity multiple slow events of increasing size can occur before one main event.
We find that in the series with Spring C ( = 136.6 ) at low normal stress and low loading rate these events show a nearly oscillating pattern of multiple cycles that is occasionally perturbed by a fast slip event. This is also highlighted in the asymmetry where these experiments have a bimodal distribution with one mode at high asymmetry (fast events) and one mode at low asymmetry (slow events). For the highest stiffness (RST, = 246.0 ) the slip rates are similar to Spring B for the fast events, but the slow events are slightly slower and show a higher asymmetry. There are fewer slow events and of smaller magnitude, with an average stress drop that is only 2.6% of the corresponding main event.
The occurrence of slow events shows a specific temporal pattern for the highest stiffness. The temporal distribution of slip events show a log-normal distribution skewed toward the end of the cycle and they do not occur in the first half of a cycle. The probability of occurrence increases from 0.7t r onwards with a mean of 0.92t r to 0.94t r and peaks at ≈ 0.95 (Figure 5a). Then the probability drops abruptly to zero and for all experiments almost no precursor has been detected in the last moments of a cycle. Higher normal stresses shift the onset of occurrence closer to failure with a smaller variability but still with no events immediately before failure. The stress level at which the slow events occur is generally very close to the stress level of the main event ( Figure 5). The curves show a log-normal distribution for low normal stresses which changes to a normal distribution (skew ≈ 0) at higher stress level. For higher normal stresses the slow events occur around 0.98τ d , and for σ N = 5 kPa at higher levels of 0.99τ d . Again we see a strong increase in occurrence up to a certain level of stress and an absence of events at stress levels very close to failure strength (>0.99τ d ).

Interevent Times
In experiments where a unimodal distribution of asymmetry is found, it is straightforward to define the interevent time as the time between the individual events. But for bimodal distributions it is more complex. Therefore we use the term "recurrence time t rec " for the time between any events (denoted by i ) and the term "reloading time t rel " for the time between the fast events (denoted by d ). This results in the following definitions: In general the interevent times decrease with increasing normalized loading rate in an exponential fashion. The interevent times t rec and t rel are essentially the same for low stiffness setups (Spring A + B) because there is only one experiment in these series with slow events. For Spring C there are only few fast events which leads to a large error for the reloading times ( Figure 6a3). Furthermore, there is a strong influence of slow event on the variance of recurrence times for the highest stiffness (RST, Figure 6b4). The power law exponent is slightly lower than n = −1 which means that there is a stronger decrease in reloading or recurrence time than what would be expected if there would be a direct correlation. Only the evolution of reloading time, that is for fast events, for the highest stiffness shows a power law exponent of n = −0.98 ± 0.04 which indicates that the occurrence of fast events is directly proportional to the normalized loading rate.

Rate-and-State Parameters
The healing rate b is determined from the change in peak stress after increasingly larger hold times (Equation A6). We find that all stressed SHS-tests show a positive healing rate (Figure 7a). The mean healing rate from all fits is b = 0.0057 ± 0.0005 which indicates time-dependent strengthening of the granular fault over time. There is no apparent correlation of healing rate b with loading velocity v L (Figure 7c). However, we observe a higher healing rate for a low normal stress of σ N = 5 kPa. Statistically it is not significantly different in the 95% confidence band in comparison with the other series due to a relatively high error for all fits at this normal stress ( 5 = 0.007 ± 0.003 , Figure 7a1). But the individual b values all plot outside the 95% interval of the mean fit of all data sets combined (black dotted line in Figure 7d).
The direct effect a is derived from two approaches. The first uses the offset in y-axis intersect of the peak stress change Δμ peak (Equation A8). This effect is clearly visible in Figure 7a where the average peak stress change increases consistently for increasing loading rates while the slope stays constant (results for 0.05 generally plot lower than those for = 0.5 ). From the average increase in peak stress change with increasing loading velocity we compute a direct effect a = 0.0074 ± 0.0031. Using the previous estimate of b = 0.005 7 ± 0.0005 from the change in peak stress after hold we calculate a first (b − a) = −0.0017 ± 0.0032 from this observation only. The second approach exploits the selection of loading velocities with respect to the hold times so that Equation A7 is fulfilled (after Beeler et al., 2001). The average (b − a) = 0.0087 ± 0.0029 fitting all possible combinations from all experiments (Figure 7e). Using b from the previous estimate we arrive at a direct effect a = −0.0030 ± 0.0030.
Another important observation for rate-and-state friction is the change of stress during the hold phase Δμ hold (Figure 7b1-4 and 7f). The data set is very noisy for this observation. We observe a weak correlation of hold stress change over time with increasing normal stress. At low normal stress (σ N = 5 kPa, Figure 7b1) the data set shows a negative slope which becomes smaller at σ N = 10 kPa ( Figure 7b2) and changes to a positive slope for σ N ≥ 15 kPa (Figure 7b3+4). However, the estimated errors for these fits are quite large and while on average the hold stress change is negative Δμ hold = −0.02 ± 0.07 it is not significantly different from zero (Figure 7f).

Additional Observations From STSL Experiments
We observe an increase in peak strength with increasing reloading time for the STSL experiment series. Plotting the reloading time t rel against peak frictional strength at failure τ p a log-linear increase can be observed (Figure 8a). The observed slope ranges from β = 0.0083 to β = 0.0130 and indicates a time-or rate-dependent healing with a similar order of magnitude as the healing rate b. In addition, we observe a decrease in average frictional strength with increasing loading rate v L (Figure 8b). The average slope of all four stiffnesses is negative ranging from (α − β) = −0.0027 to (α − β) = −0.0067 indicating velocity weakening conditions. The extremely low error for the individual fits is the result of the large amount of data points for each experiment (>10 million) because the complete time series is used for fitting. This drastically narrows the confidence band for the slope. Furthermore, comparing the reloading time t rel with the loading velocity v L we find a power-law dependency with exponents B > −1 (Figure 8c). This shows a longer reloading time than extrapolated for the simple increase in loading velocity which is another indicator for time-or rate-dependent healing.

Similarity of Constitutive Parameters
In rate and state friction three key parameters are determined, the direct effect a, the healing effect b, and the state evolution variable ϕ (Dieterich, 1979a;Marone, 1998). From our type of experiments we can not observe the evolution of friction directly because our system is inherently unstable. This is due to the system stiffness k M which is below the critical stiffness k c .
In this study, we combined various indirect methods to estimate the rate-and-state parameters from SHS tests. We can conclude from the STSL tests that the material is definitely velocity weakening at the conditions in the ring shear tester. This is further confirmed by (b − a) = 0.0087 ± 0.0029 using the peak strength Δμ at variable hold time ℎ2 ℎ1 and variable reloading velocity 1 2 (Equation A7 in this manuscript, see also Figure 9 in Beeler et al. [2001]). Additional observations show a decrease of the average frictional strength with increasing loading rates resulting in (α − β) < 0 (Section 3.5.1) in the same order of magnitude as for the estimates of (b − a). We find that a from the y-intersect method yields a = 0.0074 ± 0.0031. However, this value is not significantly different from the related estimate of b from peak stress change of b = 0.0057 ± 0.0005. As a result the calculated (b − a) = −0.0017 ± 0.0032 is indicating velocity neutral to velocity strengthening because it is not significantly different from zero. When using only the estimates from Δμ versus t h we would even arrive at a negative a = −0.0030 ± 0.0030 that is also not significantly different from zero and would represent physical impossibility (Equations A6 and A7). For rate-and-state it is necessary that all parameters a, b, and D c are positive. Rice et al. (2001) argue that it is for instance essential that a > 0 in order to get well-posed results from the quasi-static analysis of rate-and-state equations with the perturbation method. Better results could be achieved by directly Loading velocity v L compared to reloading time t rel . A negative power-law coefficient that is larger than −1 highlights longer reloading times than normal. The legend in a 1 ) applies to all subplots, the confidence band is derived from the covariance of fit (2 = 2 √ 2 ) . inverting the STSL tests or the SHS tests. With this approach one could also get an estimate for the critical slip distance D c . This was not possible because the machine stiffness was not high enough to produce real steady state slip at our conditions which lead to very unstable inversion results.
Assuming that the change β in peak frictional strength μ p with increasing reloading time t rel (Figure 8) is similar to the healing rate b we find that in the STSL experiments the healing rate b ≈ 0.0111 ± 0.0011 and (a − b) ≈ −0.004 266 ± 0.000007 which yields a direct effect a ≈ 0.0068 ± 0.0011. These values are in the same order of magnitude as the other estimates but show a positive direct effect. These observations however include creep and transient slip events during the reloading phase which influences the estimate of b. The estimated weakening (a − b) additionally includes the effect of the stick-slip cycles which distort the calculation of the mean friction. Overall, the values are in the same order of magnitude as the estimates from our SHS-tests and agree well with literature values (Carpenter et al., 2016;Marone, 1998;Marone et al., 1990).
The healing rate b = 0.0057 ± 0.0005 which is equivalent to a frictional strengthening rate β = 0.0122 ± 0.0005 in log 10 -space is at the upper estimate of natural faults and fault rocks (e.g., Alpine Fault or Scheggia Fault in Carpenter et al. [2016] or other data in Marone et al. [1990] and Marone [1998]). This means that the analogue fault material shows a similar amount of time-dependent healing that is observed for natural samples in rock mechanical tests. The underlying physical process is different for analogue materials, although a certain amount of granular mechanical strengthening due to grain rearrangement is probably also present in a dry natural fault. Therefore, the glass beads are found to be usable for small-scale seismotectonic models under analogue modeling conditions. The change of frictional strength over time and average frictional strength are similar to a natural fault and can be used to simulate seismic cycles with dynamic similarity. Due to the higher healing, which is ≈ 30% higher than most rocks, the analogue seismic cycles can be shorter in comparison with natural examples in order to represent a scaled model.
In addition, by qualitatively matching rate-and-state parameters to our SHS-tests we find that the critical slip distance D c is on the order of 10 −1 mm. Direct fits of the SHS-Tests obtained a D c ≈ 200 μm, by assuming an extremely low loading velocity ( = 10 −324 , smallest float represented in NumPy) during hold to obtain valid results during the hold phase. But because the results for the other parameters a and b by direct fitting using non-linear least squares were not stable, these approximations are not statistically sound. Nevertheless, these values are reasonably close to values found in rock mechanical tests which vary from 2 to 100 μm (Dieterich, 2007, and references therein).
There is no statistically significant difference in the estimate of (b − a) from soft and stiff systems, as expected for a material property. The scaling of strength at the onset of slip is consistent with the findings of Beeler et al. (2001) who show the same type of scaling. The scaling coefficient typically attributed to natural rocks or gouge in the seismogenic zone, is in the same range (0.011-0.015 [Beeler et al., 2001]; ≈0.01 [Scholz, 1998]; 0.001 to 0.01 [Dieterich, 2007]). Other analog model studies have used (b − a) values in the same range to model seismotectonic processes with other materials (gel on sand paper: 0.028 [Corbi et al., 2013]; rice: 0.015 [Rosenau & Oncken, 2009]; cacao, ground coffee, and others [Rosenau et al., 2017]). Therefore, we consider our models to be dynamically similar to the natural prototype, to rock deformation experiments in the MPa-range (e.g., Tullis & Weeks, 1986), and to numerical simulations of rate and state friction (e.g., Ferdowsi et al., 2013).

State Evolution During Hold Phases
To test which of the state evolution laws best describe our experimental data we take a semi-quantitative approach which considers certain observations, such as the evolution of stress during a hold phase, or the behavior in unstressed SHS-tests. During a hold phase the state θ of the system changes according to a certain relationship (Appendix A2) which leads to a change in frictional strength of the fault. Beeler et al. (1994) state that purely time-dependent healing, which is given by the Aging law (Equation A2), is independent of stiffness while the Slip law (Equation A3) shows a dependency on stiffness because it requires active fault slip during healing. For our SHS-experiments we did not systematically vary stiffness but a change with normal stress was observed. The experiments at lowest normal stress show a higher than average healing rate (Figure 7d). Additionally, these experiments show a decrease in stress during the hold phase (Figure 7b1). We attribute this to enhanced creep which is promoted by the low normal stress. As a result, the healing is amplified by larger amounts of slip in our experiments. In terms of constitutive laws this would mean that our material is better characterized by the Slip law than by the Aging law. A simple experimental test is to use additional data from unstressed SHS-tests (Marone, 1998). However, preliminary experiments with unloading showed that it is technically not feasible to do unstressed tests with our testing apparatus because the loading and unloading of the samples take too much time in comparison to the hold durations. We observed less healing for these tests but the data set is very noisy so the results are statistically not relevant.
For most experiments we observe an upwards step in stress immediately after the onset of holding. The step is followed by an exponential decay to a lower residual stress which shows a decay rate λ = 2.1 ± 1.2 s −1 . The stress decays to the residual stress after hold within less than 3 s (<1% difference). This indicates that creep due to shear stress quickly dissipates and the sample is not slipping along the shear surface during hold. The thickness of the granular packages first oscillates at a frequency of around 19 Hz and stabilizes to a constant value within the same time as the shear stress. These observations indicate that there is no measurable slip during hold which means that slip during hold is minute and therefore the Aging law would be more appropriate.
In spite of many points that speak for the Aging law, according to Bhattacharya et al. (2017) it is not sufficient to only look at the evolution of stress during hold phases but also to model the experimental time series directly, by numerically inverting for the rate and state variables to find better approximations. Strong stick-slip effects, probably due to insufficiently high machine stiffness, prevented a direct fit of Equation A1 to the data to estimate rate-and-state parameters from classical velocity stepping. Therefore, we cannot exclude the possibility that other state evolution laws apply and additional experiments with different stiffnesses and numerical modeling of the actual data are needed to clarify this finding.

Micromechanical Processes Creating Small Events
Granular material gains shear strength due to force chains oriented in the direction of the maximum stress (Cates et al., 1998). Depending on the number, length and orientation distribution of such chains shear deformation might be stable or unstable. Stick-slip is therefore interpreted as a cyclic setup and breakdown of force chains, the frequency and size of which should be a function of grain size distribution (Mair et al., 2002). Furthermore, granular materials exhibit so called "jammed states," where jamming is induced at high packaging density or by application of shear stress (Bi et al., 2011). We corroborate this view as large slip events are associated with compaction while the interseismic period is characterized by accelerating creep and dilation (Figures 2b and 2c).
The normal stress is one of the critical factors that control the creep threshold of the system. For low normal stresses it is easier for the grains to rearrange during the creep phase. First, this results in higher background slip of grains that exhibit a much lower normal stress along their contacts and can easily slide along each other. Second, the ratio of normal stress to dilatational stress, that pushes the grains apart when sliding over the rough internal shear zone, is smaller. Therefore, the force chains are less effective in strengthening the material at low confining pressures.
The occurrence of small slip events is in accordance with other studies that show transient effects during the transition of the stick phase to dynamic slip (Ferdowsi et al., 2013;Leeman et al., 2015;Nasuno et al., 1998). Because they are much smaller than the main events it is suggested that the events are the expression of internal reorganization in the granular material. During this internal deformation the grains are jammed and the force chains are rearranged into a more stable configuration. Although creep continues, the newly formed granular package is stronger than the previous package and therefore a short period of quiescence without slip events occurs ( Figure 5). This rearrangement can occur several times during the late interevent phase (multiple slow events before major event in Figure 2b). If the internal structure reaches a critical threshold, probably determined by the contact ratio and packing density, a runoff process starts and the system changes from creeping to dynamical slip. Another possibility is that during creep the local properties in the shear zone could change transiently and therefore result in a change of stability. Leeman et al. (2015) show that the experimental stiffness changes due to fabric formation and conclude that these are important factors for slip mode.
Other studies have shown a similar system behavior that is attributed to intermittent criticality (Ben-Zion et al., 2003). In contrast to the self-organized critical system, intermittent criticality implies a cyclic evolution of the fault zone, whereas the former only gives a general statistic fluctuation around the critical state (i.e., failure criterion). If we apply the concept of intermittent criticality, the small precursors are the expression of small scale stress perturbations along the fault zone. Overall the stress field within the granular fault zone homogenizes by increasing rearrangement of force chains, that explains the increasing frequency of events up to a certain point.
Then the system is largely homogenized and is in a critical state, very close to failure, which is comparable with the state of stress in the lithosphere (Sornette et al., 1990). This behavior has also been observed for the temporal and spatial clustering of smaller earthquakes (Hainzl, 2003).
The behavior of dilation during the interevent cycle is even more complex and it is difficult to assign a direct relation to micromechanical processes. The observed increase in wavelength of the small amplitude oscillations could indicate a smoothing of the internal fault surface, leading to a flatter frictional response. This could be attributed to localization of strain during the interevent period. Immediately after a slip event the grains are moving in a broader region leading to many small motions. Over time the deformation localizes into narrower band where stress is transmitted slower and more regularly. The discrete upward and downward steps might be artificial, or the result of sensor noise. Although this machine has been used for several previous measurements of simple frictional properties, such sharp changes in thickness have not been observed. Additional testing by slowly moving the lid without a sample during shear on the same order of magnitude also did not show such a behavior. However, the strong reproducibility over multiple cycles indicates that mechanical explanations can be valid, too. For example, internal reorganization of the granular packaging leads to discrete conformations of packaging with different densities that are characteristic for each state of the system.

Criticality of Analogue Fault
From the determined rate-and-state parameters in the SHS tests we can now derive the critical stiffness to evaluate how close the STSL series experiments are to the bifurcation from stable and unstable slip. Because normal stress is constant in each series we use a stability criteria according to Dieterich (2007): According to Heslot et al. (1994) the critical stiffness normalized by normal load σ N (in their case slider mass M) is a slowly decreasing function depending on loading velocity v L . Accordingly, we correct for loading rate v L with respect to the loading rate of the SHS-tests 0 = 0.52 and the scaling factor α = 10: Utilizing both estimates for ζ = (b − a) = {0.0043, 0.0087} (Section 3.5 and 3.5.1), three different estimates for D c = {50, 100, 200} μm and normalizing by normal stress the normalized critical stiffness k c ranges between 60 and 340 mm −1 . Comparing this with the normalized machine stiffness k N (Figure 3) we find that most experiments show k > k c and therefore should only show stable sliding (Dieterich, 2007). However, due to the material inside the machine the actual stiffness of the system is lower. If the normalized reloading stiffness k R is used instead (Figure 9), the experiments with lowest stiffness now show k ≤ k c and the experiments with higher stiffness Figure 9. Slip modes in the k − v space represented by the normalized reloading stiffness which includes the material's effect. The critical stiffness k c is calculated from the rate-and-state parameters from SHS-experiments at maximum stiffness (Section 3.5) and explains the transition from unstable to stable slip mode for Spring C.
now show k ≈ k c . The use of in situ measured reloading stiffness is important for describing stick-slip effects because during shear stiffness evolves as a function of shear displacement (Leeman et al., 2015).
Consequently, the experiments with Spring C (△ in Figure 9) now fit with the change of unstable to stable slip when transitioning the stability boundary (gray boxes in Figure 9). It is unclear why the stiffest experiments (RST) still shows unstable slip for higher loading rates. One would expect that with a higher stiffness than Spring C one could see an earlier transition to stable sliding. From these observations we also think that it is safe to assume that D c is on the order of 100-200 μm which is roughly the radius of a glass bead r GB = 150…200 μm. Due to the uncertainties in the estimation of k N and the high uncertainty for D c the values for k c and the location of the stability boundary are not well constrained. The uncertainty for D c might result from the possibility that D c is not constant because the thickness of the active shear zone might change during the experiment which is a primary factor for the scaling of D c (Marone & Kilgore, 1993). Furthermore D c shows a dependency on slip velocity which could not be studied with our setup (Hatano, 2009). Nevertheless, the fits of the SHS tests and STSL tests yield results that seem valid and the stability boundary lies within the same order of magnitude.

Modeling of Slow Slip Events
For natural examples the slip mode is usually identified by the relation of seismic moment M 0 and characteristic duration T (e.g., Gomberg et al., 2016;Ide et al., 2007). Regular earthquakes show a much shorter duration (M 0 ∝ T 3 ) in comparison with slow earthquakes (M 0 ∝ T). The scaling relation leads to a characteristic separation between events of equivalent seismic moment M 0 . We observe a similar separation of slow and fast events which is most prominent for the experiments at highest stiffness. The occurence of small events and larger events was previously also described for baking flour in Leeman et al. (2015). Depending on the actual seismic moment the separation in nature is between 2 and 5 orders of magnitude. In the experiments at highest stiffness we observe a separation in average slip rate of 2-3 orders of magnitude. The difference in peak slip rate is up to 5 orders of magnitude. This separation is smaller but still statistically significant for lower machine stiffness where we also find oscillating slip modes. Similarly, the frictional stress drop that can be seen as a proxy for seismic moment in our experiments, shows a separation of 2 orders of magnitude which indicates that the dynamics are different for slow and fast events.
Several studies highlight the relationship of transient slip events promoting seismic activity (A. Kato & Ben-Zion, 2021, and references therein). The results from our study suggest that glass beads as analogue fault gouge shows a large variety of slip modes not only at high mean stresses as previously found by others (Cain et al., 2001;Cui et al., 2017;Dieterich & Kilgore, 1996;Mair et al., 2002). Therefore, the usage of glass beads in small scale analogue models with intermediate to high stiffness at stresses in the kPa range is suitable to model seismic fault behavior. For example, the glass beads can be used as fault gouge in between two elastic blocks in gel-slider type models (similar to Corbi et al. [2011]) or in more complex models. If the stiffness of the model is adjusted correctly (e.g., in the range of Spring C) several slip modes are possible depending on the normal stress on the fault. The normal stress regime on the fault can then be designed by geometric orientation of the fault with respect to the loading direction. More complex fault geometries, possibly forming fault networks, then lead to transient stress on the individual faults and thereby altering slip mode. These setups can be used to study the complex interplay of fault geometry and slip on individual faults, while retaining dynamic and kinematic similarity. Due to transient stress changes individual analogue faults might change their slip mode from pure stick-slip to creep and thereby changing system dynamics and the activity of other faults in the system.
Especially the temporal distribution of slow events between the occurrence of fast events highlights the possible application of glass beads as analogue fault gouge. The increased probability of slow events toward the end of the reloading time and the high stress level is similar to the behavior observed for large fault systems that show increased seismic activity toward the end of the seismic cycle (Bürgmann, 2018;A. Kato & Ben-Zion, 2021). The abrupt decrease in probability just before failure might be an expression of strong locking due to jamming by shear (Bi et al., 2011) and might favor fast slip events by driving the frictional strength above a critical threshold.
In general we find that the "precursory" phase where the fault shows signs of imminent failure starts relatively early, with the onset of creep at about 50% of the reloading time. This is similar to long preparatory phases of large earthquakes (Bouchon et al., 2013) and to recent findings by Igarashi and Kato (2021). In our case the slow events do not always act as precursors because for certain conditions they occur in a repeating pattern at high stresses (e.g., Spring C at low normal stress) and therefore are reminiscent of "similar earthquakes" (Igarashi & Kato, 2021). Although the stress drop magnitude slowly increases over several repeating events, they do not show a clear threshold for which the slow event grows into a large event. Consequently, they only pinpoint that the fault is close to failure but not that the fault will fail with a fast slip event after a certain type of slow event.
The slip behavior of the granular analogue is the result of the interaction of friction with the complex network of force chains that is created in the sheared bulk material (Cates et al., 1998;Daniels & Hayman, 2008). This micromechanical mechanism creates the macroscopic behavior that can be described with rate-and-state friction.
Our findings support the hypothesis that rate-and-state like dynamics are the expression of processes that emerge close to the criticality boundary. Similar kinematic observations can be made for a range of microphysical processes and conditions (e.g., Hecke, 2009;Kabla et al., 2005;Lemaître & Caroli, 2009;Papanikolaou et al., 2013;Scuderi et al., 2015). Denisov et al. (2016) show that force fluctuations in granular matter, of which our experiments but also fault gouges are a part of, show scaling relations that are akin to critical phenomena. Furthermore, the size and stress distribution of slip events within granular materials at the same conditions as in our experiments are found to be universally related with slip events on multiple scales and follow similar scaling relations (Uhl et al., 2015). The underlying physical mechanisms for these observations can be quite different but yield the same results on a macroscopic scale. This description is similar to rheological equations that purely describe the observed relation of motion and stresses, while the actual deformation mechanism is different in various fluids. Therefore, we find that our model has widespread possibility of application within seismotectonics, engineering and hazard assessment for earthquakes and landslides.
The advantage of using a granular analogue is the simplicity with which observations can be made. The analogue modeling approach features lower stresses which simplifies the design and construction of the testing machine. This increases the available parameter space because it is relatively easy to change the system stiffness using springs. For rock mechanical testing apparatuses the change in stiffness is limited to a smaller range that is either accessible through adding rubber blocks or by artificially changing the servo-hydraulic systems to mimic a different stiffness (Beeler et al., 1994;Leeman et al., 2015) The results from this study can be used to improve current numerical models of granular gouge but can also directly be applied in improved seismotectonic scale models (Blank & Morgan, 2019;Rosenau et al., 2017). The glass beads show a slip behavior that naturally emerges from their frictional properties. This can be exploited for larger analogue models to model fault slip in a geometrically complex fault system. The range of available temporal and spatial scales, as well as the self-consistent scaling behavior allow the application in many fields where rate-and-state friction is a dominant process such as landslides, glacial motion, mass movements and lithospheric deformation (Jerolmack & Daniels, 2019). In comparison to numerical simulations the use of an analogue model allows to inherently link the spatial and temporal scales without having to rely on parametrization and grid based methods. The analogue approach allows to model small scale processes, such as earthquakes within a fault zone over many seismic cycles and over a much larger spatial scale within a shorter period of time than numerical simulations of similar complexity.

Conclusions
We have used an annular shear apparatus to characterize the stick-slip behavior of a granular fault zone analogue composed of glass beads. Using slide-hold-slide tests the rate-and-state properties have been qualitatively evaluated and quantified. The healing rate is found to be b = 0.0057 ± 0.0005. The direct effect a is quantified by two approaches and found to be a = 0.0076 for estimates from the change of peak stress with increasing reloading velocity in SHS-tests, while the procedure by Beeler et al. (2001) with a specific ℎ -ratio yields a = −0.003 0 ± 0.0030. While the combinations of our estimates of a and b are subject to relatively large statistical error the resulting estimates for (b − a) are generally ≥ 0 we can safely conclude that the glass beads show velocity weakening characteristics. Also the observations in the STSL tests and the reduction of average frictional strength μ with increasing loading rate v L support this. Due to the evolution of stress during the hold phase we find the Aging law to be slightly more appropriate for our material but the results from SHS-tests only remain inconclusive. The critical slip distance D c is estimated to be in the sub-mm range but can not be quantified with the presented setup. The effect of machine stiffness k M , loading rate v L and normal stress σ N on the slip mode is studied. We find a large variety of slip modes ranging from pure stick-slip, oscillations to bimodal slip modes within the same experiment by only varying certain extrinsic parameters which fits well with the stability boundary derived from the rate and state parameters. Low stiffness, low loading rates and high normal stresses favor pure stick-slip with small amounts of interevent creep. Higher stiffness, especially in combination with low loading rates, leads to a bimodal distribution of fast, large events that are preceded by slow, small events. The slip events reproduce typical characteristics that have been observed in similar experiments in other experimental setups with different boundary conditions and materials allowing to generalize the observations to natural occurrences of earthquakes. In the experiments, rearrangement in the granular package is the major micromechanical process which distributes and dissipates stress during shear. This drives the system closer to criticality leading to the observed precursory strengthening and the short period of quiescence before a large slip event. We conclude that the small transients can strongly affect the statistical characteristics of a single fault zone system and makes the material suitable for the use in larger analogue modeling setups that model seismotectonic deformation with a higher geometrical complexity. The small scale events during the precursory phase are the expression of distributed fluctuations of the system in a critical state. Further examination of these fluctuations and their correlation with the generation of large events may give important constraints on the predictability of slip events (as suggested by Ben-Zion et al. [2003]). Furthermore, the higher complexity with differing slip modes due to the characteristics of the glass beads could provide additional insights into the system behavior and the interaction of faults in analogue models that are closer to the behavior of a natural fault zone. The results from this study shed light on the micromechanical mechanisms from which rate and state-friction emerges. Therefore it can act as a benchmark for numerical models of fault zones, alleviate the design of more complex analogue models and helps interpreting kinematic natural observations of fault slip.

Appendix A: Data Analysis
The experimental data is examined using a combination of classical event detection, statistics and machine learning. To analyze the occurrence and properties of the slip events we employ a peak detection that is based on a minimum stress drop threshold for each experiment. Then we extract certain characteristic points in the stress curve, these are highlighted in Figure 2b.

Appendix A1. Picking and First Order Properties
For first order characterization of the experiments we use a simple peak detection to find slip events in the stress curve. For this the data is split into sets of equal loading rate, normal stress and stiffness. A fixed threshold for stress drop per set facilitates the detection of sudden changes in shear stress. Fine tuning this value enables the detection of large and fast, but also of small and slow events by searching for positive and negative peaks in the data. The result of peak detection is cross checked by manual inspection of the stress curves and the detected peaks (Figures 2a and 2b).
The point X of maximum stress immediately before failure is denoted by X p indicating peak values, the point of minimal stress after a slip event is indicated by X e accordingly (Figure 2c). In most cases X is replaced by the appropriate physical quantity such as shear stress τ or velocity v. A velocity threshold defines the separation of non-dynamic and dynamic slip events. A slip event is considered dynamic when at any point during a decrease of shear stress the slip velocity v s is higher than the threshold v d . In this study we used the maximum loading velocity of = 0.02 as the threshold. This allows the definition of onset of dynamic slip, denoted by X d , maximum slip X m where slip velocity is at its maximum and the end of dynamic slip X f where the slip velocity drops below the critical value. These points now define a full cycle, which we see as an analogue of a seismic cycle.
The full cycle is defined as the period of time between two slip events that have a dynamic phase with velocities above the threshold v d . During the majority of a full cycle the shear stress increases in a linear relation with load point displacement but deviates to a non-linear relation. This point is the onset of creep and is defined as the point where the linear trend extrapolated from the previous points deviates by more than 1%. The slope of the linear trend, calculated by least squares fitting, also defines the cyclic reloading stiffness k L which is a measure for the overall stiffness of the setup including the bulk stiffness of the granular material. This stiffness is also used for further calculations of the criticality using the rate-and-state framework.
Assuming an overall elastic behavior of the granular material when completely locked, we can estimate the amount of creep either as overall proportion or as instantaneous creep. Overall creep is calculated by linearly extrapolating the shear stress increase over the full cycle using k L as a slope and then relating the predicted and observed point of failure. Similarly the instantaneous creep is calculated by a similar method but doing a point wise calculation.

Appendix A2. The RSD-Formulation
Following Dieterich (2007, and references therein), shear stress τ evolves as a function of effective normal stress σ, load point velocity v L and a set of experimentally derived parameters μ 0 , a, b, θ, D c in relation to a reference load point velocity * : This is a heuristic description of the change in shear stress in response to a change in slip velocity. The parameters in Equation A1 are usually derived experimentally using velocity stepping tests where the system sliding at a given reference load point velocity * under stable conditions (̇= 0) is perturbed by setting a new loading velocity v L . This prompts a direct reaction of shear stress that is proportional to the magnitude of the perturbation * and the constant a ("direct effect"). Following this immediate reaction, shear stress adjusts to a new level defined by the evolution of state over time in relation to the new loading velocity normalized by the characteristic slip distance * and a constant b ("evolution effect"). The evolution of state over time ̇ is defined by choosing one of the following evolution laws:̇= The Aging law (Dieterich, 1978) and Slip law (Ruina, 1983) are the most commonly used state equations, while the Kato law (N. Kato & Tullis, 2001) and Nagata law (Nagata et al., 2012) are more recent developments. All of the above laws contain a slip dependent component which is expressed in the term . Consequently, the Slip law (Equation A3) does not show any healing when there is no slip (v L → 0) which can be tested through unstressed SHS tests where the sample at rest is not under stress and thus slip along grain contacts is hindered. The Aging law (Equation A2) shows constant healing at rest due to the 1 in the equation. This purely time-dependent effect makes the Aging law fit better to experimental data (Beeler et al., 1994). However, for large velocity steps (|log 10 * | > 2) the Aging law shows a linear decay that is dependent on the sign and magnitude of the step which is not in accordance with experimental data. Furthermore, the Aging law does not fit well to the state evolution during a hold phase and needs non-constant a, b and D c (Bhattacharya et al., 2017). The Slip law exhibits a better fit for large velocity steps and for the evolution of state during a hold. This resulted in a reformulated version of the Slip law that accounted for time-dependent healing by N. Kato and Tullis (2001)  ). A further improvement to the previous laws has been proposed by Nagata et al. (2012) which incorporates the relation of the "evolution effect" b and a new constant c, as well as the normalized stressing rate ̇ (Equation A5).

Appendix A2.1. Tests and Derived Quantities
We performed SHS tests with stressed hold phases (Marone, 1998) to determine the direct effect a, rate of healing b and the appropriate state law.
Due to the evolution of state during a hold, the frictional resistance μ p of a granular medium increases with the natural logarithm of the hold time t h and gives rise to the healing rate b (Bhattacharya et al., 2017): This increase in μ p is measured with slide-hold-slide tests. During the first slide phase a steady-state value of μ s is established. Then the machine is stopped for a certain time t h , either under stress or unstressed. Then the sample is resheared and the increase of peak strength with respect to the previously established stable sliding resistance is measured as Δμ p = μ p − μ s . Relating the results to ln t h a linear increase can be measured, which is the healing rate b. Furthermore, the loading velocity v L plays an important role in the magnitude of frictional resistance after loading μ p but should not influence the healing rate b (Beeler et al., 2001). Experiments at increasing v L therefore show the same slope b but increasing Δτ p .
This effect can be exploited to determine the direct effect a from slide-hold-slide tests by using a specific spacing between the realized loading rates (Figure 9 in Beeler et al. [2001]). If the ratio of the loading velocities 1 2 is equal to the ratio of hold times ℎ2 ℎ1 (Equation A7) then the increase in Δμ p is proportional to ⋅ ln 1 2 . We estimated the direct effect a from the average increase in Δμ p over all realized hold times using a set of SHS-Tests that fulfill Equation A7. Other approaches to estimate (a − b) were also tested with our data. Corbi et al. (2013) defines (a − b) using the peak friction μ p and sliding velocity v L (Equation A9). For this we determine the frictional resistance μ p at the peak point τ p just before a dynamic failure because at the peak, there is a plateau of shear stress and therefore the current slip rate of the fault equals the loading rate v S = v L .  Figure A2. Asymmetry for all events and all experiments. The legend in a 1 applies for all plots. Each row represents experiments of the same stiffness which is also indicated by the individual markers. Color highlights the different normal stresses which has an additional influence besides the stiffness and normalized loading rate. Figure A3. Stiffness measurements for all springs and the basic ring shear tester (RST). Data was fit in the region of 20 ≤ F s ≤ 180 N which approximately corresponds to a shear stress range of 1,400 ≤ τ ≤ 13,000 Pa. The given error is 2σ as given by the covariance of the least squares fit. The legend in (a) applies to all plots.