# Micro-slips in an experimental granular shear band replicate the spatiotemporal characteristics of natural earthquakes

David Houdoux, Axelle Amon, David Marsan, Jérôme Weiss, Jérôme Crassous

► To cite this version:

David Houdoux, Axelle Amon, David Marsan, Jérôme Weiss, Jérôme Crassous. Micro-slips in an experimental granular shear band replicate the spatiotemporal characteristics of natural earthquakes. Communications Earth & Environment, 2021, 2 (1), 10.1038/s43247-021-00147-1. hal-03280560

HAL Id: hal-03280560

https://hal.science/hal-03280560

Submitted on 7 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Micro-slips in an experimental granular shear band replicate the spatiotemporal characteristics of natural earthquakes

David Houdoux<fnref n="1" />, Axelle Amon<fnref n="1" />, David Marsan², Jérôme Weiss³ & Jérôme Crassous<fnref n="1" />✉

Memory effects in seismology—such as the occurrence of aftershock sequences—are implicitly assumed to be governed by the time since the main event. However, experiments are yet to identify if memory effects are structural or time-dependent mechanisms. Here, we use laser interferometry to examine the fluctuations of deformation which naturally emerge along an experimental shear fault within a compressed frictional granular medium. We find that deformation occurs as a succession of localized micro-slips distributed along the fault. The associated distributions of released seismic moments, as well as the memory effects in strain fluctuations and the time correlations between successive events, follow exactly the empirical laws of natural earthquakes. We use a methodology initially developed in seismology to reveal at the laboratory scale the underlying causal structure of this behavior and identify the triggering kernel. We propose that strain, not time, controls the memory effects in our fault analog.

<footnote n="1">¹University Rennes, CNRS, IPR (Institut de Physique de Rennes)—UMR 6251, Rennes, France. ²Université Savoie Mont-Blanc, CNRS, IRD, IFSTTAR, ISTERRE, Le Bourget-du-Lac, France. ³IsTerre, CNRS/Université Grenoble Alpes, Grenoble, France. ✉email: jerome.crassous@univ-rennes1.fr</footnote>
Earthquakes are natural phenomena displaying scale-free statistics¹. Empirical power laws are observed for the distribution of their moments (Gutenberg–Richter’s law), rupture lengths and durations, rupture slips², temporal³ and spatial correlations between earthquakes⁴, which also express through a decaying rate of aftershocks (Omori’s law)⁵, characterized by a scale-free (sub)diffusion⁶,⁷. Understanding the origin of those laws as well as reproducing them at the laboratory scale remains nowadays major issue. From a fundamental point of view, these scaling laws are reminiscent of nonequilibrium dynamics and critical phenomena⁸,⁹ and raise the question of the existence of a universality class to which earthquakes would belong. Among mechanical systems, possible candidates for such a universality class are those for which deformation occurs through avalanches. The mechanical response of those systems is characterized by intermittent dynamics alternating slow elastic loading and rapid sliding and relaxation, leading to jerky dynamics and/or stress drops.

Power-law distributions of slip sizes or relaxed energies have been evidenced experimentally in various systems, such as the stationary propagation of a fractured front in a heterogeneous material¹⁰, compression of heterogeneous materials¹¹,¹², or systems involving frictional sliding on or within a granular media¹³,¹⁴. The common basic ingredients underlying the dynamics of those different systems are the existence of the material disorder and the decomposition of the dynamics in elementary events localized both in space and time, coupled together by elasticity. Progressive evolution of avalanche size and duration statistics has been reported for different heterogeneous materials¹² or granular media¹⁵ upon increasing the loading up to a macroscopic yield or failure stress at which scale-free statistics are observed, arguing for a “stress-tuned” critical behavior fundamentally different from a self-organized critical dynamics characterized by steady-state statistics¹⁶. It has been proposed on the basis of a mean-field model of plasticity that those different systems, as well as deformed microcrystals and earthquakes, could belong to the same class of universality¹⁶,¹⁷. This, however, was mainly addressed from an analysis of the size distribution of avalanches as well as their average shape. In the case of earthquakes, the stress-tuned critical hypothesis was argued¹⁶ on the basis of dependence of magnitude distributions on the slip direction on the fault plane (the rake angle), which gives indirect information about the differential stress acting on the fault¹⁸. However, if seismic moment distributions appear to be indeed exponentially tapered at very large scales, the associated upper corner magnitude was found to be independent of the region or the depth interval considered, or of the plate velocity, i.e., to be rather “universal”¹⁹. On the other hand, an important and ubiquitous feature of brittle deformation in the crust is the existence of aftershocks, which occurrences are also governed by scale-free laws. Much fewer attempts have been done to model those spatio-temporal correlations, through the introduction of a memory mechanism such as a slow healing of frictional properties²⁰, or viscoelastic relaxation²¹. Similarly, lab experiments reproducing the clustering of events in time and space remain scarce.

As a matter of fact, systems displaying avalanches can have also fundamental differences that limit the pertinence of a universal picture. As already noted, systems in a stationary regime must be distinguished from those for which the spatio-temporal dynamics are evolving (stress-tuned). In the case of frictional granular media, besides a progressive increase of the maximum avalanche size¹⁵, the spatial distribution of plastic events evolves during loading: initially homogeneously distributed in the bulk of the material, plasticity progressively localizes to form shear bands at macroscopic yield²²⁻²⁴, in which all the shear rate concentrates

afterward while the rest of the system becomes an elastic “solid”. Identifying clustering in time is only possible within those shear bands when the spatio-temporal organization becomes stationary. For stationary systems, the dimensionality of the active zone is expected to play a role, at least on the value of the critical exponents. One must thus distinguish tri-dimensional systems (e.g., plasticity distributed in the bulk of an amorphous material), from those where the plasticity is confined to a quasi-2D zone (a fault in the case of earthquakes, a shear band in the case of amorphous granular media), and finally quasi uni-dimensional active zones (e.g., a propagating crack front). Practically, in experimental works pertaining to the plasticity of amorphous media, it is not always clear whether the plasticity is broadly distributed in the bulk of the system or if it is localized along a shear band. Even when the geometry of the active zone is identified, most experimental set-ups are unable to fully resolve the spatio-temporal organization of the avalanches which are solely identified and studied through indirect measurements of their sizes such as acoustic emissions (e.g., ref. ¹⁴) or stress drops on a loading curve (e.g., ref. ¹⁵).

To address the challenging issue of reproducing an analog of a fault gouge at the lab scale, a straightforward approach consists of imposing the bi-dimensional geometry in stationary conditions by confining a granular material between elastic plates¹⁴,²⁵⁻²⁷. In the vast majority of those experimental studies, quasi-periodical stick-slip events with a typical size are observed, indicating that finite-size effects dominate the dynamics²⁸. The slip events then involve the whole length of the shear band and the dynamics lose their universal features. In addition, those macro-slip events are characterized by a reverse asymmetry of the activity compared to earthquakes, with foreshocks of increasing size as approaching the macro-instability, but an absence of aftershocks¹⁴, likely resulting from the finite-size effects mentioned above. A recent experiment with a quasi-2D shear cell in a stationary regime displayed intermittent dynamics sharing several features of earthquake dynamics, such as the G-R law and a power-law decay of the rate of aftershocks²⁷. While giving promising results in terms of the analogy between the dynamics of stationary sheared granular materials and that of earthquakes, it did not give a direct characterization of the localization and the spatial extension of the detected events. Consequently, the aftershock characterization amounts to a time correlation analysis of discrete events, without quantifying the underlying causal triggering. We can thus ask: is it possible to build a laboratory analog of a fault gouge where well-identified events would share all the properties of earthquakes, and more particularly their spatio-temporal, scale-free clustering properties arising from stress transfers and the resulting cascades of triggering²⁹⁻³².

Here we present experimental results obtained in a 3D granular system in a post-(macro)yield regime displaying a stationary shear band and in which finite-size effects do not dominate the dynamics. Using an interferometric method of measurement of micro-deformations we provide a direct spatially-resolved measurement of the micro-slip events that govern the frictional motion along the shear band. We are able to measure the localization, the spatial extension, and the magnitude of those events, providing the first direct experimental measurement at the laboratory scale of frictional micro-slips along a fault. We show that the statistics of those events display scale-free behavior in close agreement with earthquake phenomenology. Using a methodology developed for earthquake analysis³⁰, we go another step further compared to previous experimental studies by quantifying the causal triggering between events. We show that this underlying triggering process can explain the observed spacetime correlations in the dynamics, much like it does for earthquakes. We argue on this basis that a frictional shear band in a
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Imaging shear band fluctuations. **a** Schematic of the experimental set-up. The material is submitted to a biaxial stress test. The front face of the sample is imaged on a camera. As illumination is done using coherent light, those images display speckles. **b** Normalized deviatoric stress as a function of macroscopic deformation. $\varepsilon_y$ is the yield strain, and the gray zone is the post-yielding zone analyzed in this study. **c** Relative displacement $u_x$ of two blocks separated by a shear band as a function of the direction perpendicular to the shear band. Symbols are experimental data, the plain line is $u_x(z) = \Delta u_x \times [1 + \tanh(-2z/w)]/2$, with $\Delta u_x = 32$ nm and $w = 22d$. **d** Map of the correlation between two successive speckle images. The color of the pixel is related to the value of the correlation. Blue arrows represent the directions of displacement of the blocks with respect to the shear bands. The dashed area is the Region of Interest (ROI) for performing speckle analysis. **e** Schematic of the shear band separating two sliding blocks as composed of discrete shear events of fixed-width $w$ and of size $L \times L$.</figcaption></figure>

granular material represents a formal analog of tectonic faults, with intermittent dynamics probably belonging to the same universality class.

## Strain fluctuations inside a shear band

**Stationary shear band and strain imaging.** We use for this study an experimental set-up that consists of a biaxial cell filled with a granular material composed of an assembly of glass beads confined into a rectangular box (see Fig. 1a). The two lateral faces are deformable latex membranes which allow us to impose a confining stress $\sigma_1$ to the material. This stress is kept constant during the whole experiment. The material is slowly compressed by moving at a fixed velocity the top face with respect to the bottom face, and the axial macroscopic deformation $\varepsilon_M$, and the applied axial stress $\sigma_3$ are measured (see “Methods”).

The strain fluctuations are imaged using an interferometric technique based on diffusing wave spectroscopy (DWS). For this, the material is illuminated with an extended laser beam, and the speckle images are regularly recorded. We note $\delta\varepsilon_M^*$ the macroscopic strain increment between two successive images, and its value is $\delta\varepsilon_M^* = 5.10^{-7}$ if not otherwise specified. The speckle images are then divided into square zones, and for two successive images, the normalized autocorrelation function of the scattered intensities is calculated for each zone. Associating a color to the value of the correlation at a position, we obtain maps of correlation $g_I(\varepsilon_M, \mathbf{r})$, where $\mathbf{r}$ is the position on the observation plane, as shown in Fig. 1d. High correlation $g_I \approx 1$ (white pixels) indicates that beads are uniformly translated without relative motions, whereas low correlation $g_I \ll 1$ (dark pixels) is the signature of bead relative motions. In addition to this interferometric correlation technique, we use a conventional digital image correlation method on the speckle pattern: the displacement of zones of the speckle pattern between different images are measured, giving access to the displacement field. This measure is used to determine the mean relative velocity of blocks when shear bands are formed.

Starting from an initial condition of a material submitted to an isotropic confining pressure, the material is slowly compressed. The beginning of the compression is associated with a plastic flow spatially distributed into the sample and to an increase of the stress difference $\sigma_3 - \sigma_1$ (see Fig. 1b). We analyzed this plastic flow in previous works²² (see also Supplementary Movie 1) and we do not look further to this initial stage here. When the deformation $\varepsilon_M$ of the material exceeds a yield strain $\varepsilon_Y \approx 4.5\%$, the stress difference $\sigma_3 - \sigma_1$ is roughly constant (see Fig. 1b), and

the deformation localizes into the material, with the formation of one or two linear shear bands³³. The shear is not localized close to a moving mechanical boundary as it is the case for a Couette cell (where the shear band appears at rotor), or in a gouge confined between two rigid blocks. Here, the shear band emerges spontaneously in the system. Its orientation is linked to intrinsic properties of the material (it is linked to the deviatoric stress at failure through the Mohr-Coulomb relationship) and not to geometrical constraints. Its width is the result of the self-organization of the system: the flowing material forming the band and the solid material surrounding it is the same and this separation of phase emerges spontaneously in the system.

**Average vs. instantaneous strain.** The map of the correlation of the scattered intensity can be linked to the shear motion of the sliding blocks at each side of a band. This may easily be seen qualitatively: for this we consider a correlation map obtained in the stationary regime (see Fig. 1d), i.e., when the stress difference is in a plateau phase, and $\varepsilon_M > 5\%$. The correlation is close to 1 into the four triangular zones partitioned by decorrelated boundaries. This indicates that the material is split into four rigid blocks separated by deformed zones.

To obtain quantitative information about the shear field inside the band, we assume (this hypothesis will be discussed just below) that the motion of the beads around a point $\mathbf{r}$ is mainly a shear $\gamma_m(\mathbf{r}, \varepsilon_M) = \partial u_x / \partial z$, where $\mathbf{u}$ is the local displacement of a block with respect to another one, and $(x, z)$ are local coordinates associated to a band (see Fig. 1d for axis definition). By making this assumption, we neglect other components of the strain tensor and uncorrelated motion of the beads. If the beads move accordingly to this shear field, the decorrelation can then be related to the local shear as (see Supplementary Methods 2 and Supplementary Figs. 2–4):

$$ \gamma_m(\mathbf{r}, \varepsilon_M) = -\gamma_0 \ln[g_I(\mathbf{r}, \varepsilon_M)] \quad (1) $$

with $\gamma_0 = 2.6 \times 10^{-4}$ a constant given by the optical properties of the material. The time-averaged local deformation is defined as:

$$ \bar{\gamma}_m(\mathbf{r}, \varepsilon_M) = (1/\Delta\varepsilon_M) \int_{\varepsilon_M}^{\varepsilon_M + \Delta\varepsilon_M} \gamma_m(\mathbf{r}, \varepsilon_M') d\varepsilon_M' \quad (2) $$

The hypothesis of local shear can be quantitatively tested. For this, we integrate (2) along a direction perpendicular to the shear band, and we obtain $u_x(z) - u_x(-\infty) = \int_{-\infty}^{z} \bar{\gamma}_m(z', \varepsilon_M) dz'$. Figure 1c shows the displacement field across the shear band,
<figure><img src="image_4.png" /><figcaption>Fig. 2 Strain fluctuations and correlations. a, b Granular experiment. a Spatiotemporal correlation function $C(\delta x = 0, \delta \varepsilon_M)$ as a function of $\delta \varepsilon_M$ see Eq. (3). b Spatiotemporal evolution of the local strain $\gamma_T(x, \varepsilon_M)$ (see text) in a shear band. c, d Californian earthquake catalog. c Spatiotemporal correlation function $C(\delta x = 0, \delta t)$ as a function of $\delta t$ (see Supplementary Methods 1). d Spatiotemporal evolution of the seismic moment $M$ resulting from earthquakes from the Californian catalog projected along the main direction of the fault (see Supplementary Methods 1 and Supplementary Fig. 1 for details).</figcaption></figure>

demonstrating that the strain is concentrated into a narrow zone of width $w = 22d$. This value is in agreement with the reported values of shear band thickness in biaxial experiments³⁴. The spreading of light due to the diffusive transport of photons inside the granular material is expected to occur on a scale of the order of $l^* \sim 3d$, which is small compared to our measured value of $w$. This value of $w = 22d$ is, therefore, representative of the thickness of the actually sheared zone.

For our macroscopic strain increment $\delta\varepsilon_M^* = 5.10^{-7}$, the relative sliding of the two blocks $\Delta u_x = u_x(\infty) - u_x(-\infty)$ if found to be $\Delta u_x = 32$ nm. This value is close to the imposed value $\Delta u_x = 25$ nm that we may estimate from the displacement of the top plate and the orientation of the shear band. This agreement confirms the hypothesis that the decorrelation of the speckle pattern is mainly due to the shear motion of the beads. The small difference between the imposed and the measured value of displacements may have different origins. Shearing a granular media usually produces some dilation or compaction, creating a decorrelation of the scattered light. Also, some non-affine motion of the beads occurring in the sheared zones may produce decorrelation³⁵. Finally, our optical model of perfectly spherical beads ignores possible bead rotations.

**Memory effect in strain fluctuation.** We now analyze the fluctuations of the local shear. Since deformation is located in the shear band, we consider the transverse-averaged shear deformation $\gamma_T(x, \varepsilon_M) = \frac{1}{2w} \int_{-w}^{w} \gamma_m(\mathbf{r}, \varepsilon_M) dz$. Figure 2b shows $\gamma_T(x, \varepsilon_M)$ into the $(x, \varepsilon_M)$ plane. We can clearly see that the shear is heterogeneous both in space (i.e., along the shear band) and intermittent. In order to analyze those fluctuations, we introduce the normalized spatiotemporal correlation function:

$$C(\delta x, \delta \varepsilon_M) = \frac{\langle \gamma_T(x, \varepsilon_M) \times \gamma_T(x + \delta x, \varepsilon_M + \delta \varepsilon_M) \rangle}{\langle \gamma_T(x, \varepsilon_M) \rangle \langle \gamma_T(x + \delta x, \varepsilon_M + \delta \varepsilon_M) \rangle} - 1 \quad (3)$$

where $\langle \cdot \rangle$ is an average both on deformation and position (see “Methods” for details). Figure 2a is a plot of $C(\delta x = 0, \delta \varepsilon_M)$ as a function of the $\delta \varepsilon_M$. There is clearly memory in the deformation. Moreover, this correlation function decays as a power law $C(\delta x = 0, \delta \varepsilon_M) \sim (\delta \varepsilon_M)^{-\theta}$, with $\theta = 0.74$. Binning in space and time the moment along the San Andreas Fault system from the Californian earthquake catalog (see Supplementary Methods 1 for details), very similar spatiotemporal patterns (Fig. 2d) and correlations (Fig. 2c) are obtained. The analogy between our shear band and tectonic faults in terms of scaling laws of seismic moments, temporal clustering, and aftershock triggering, is thoroughly analyzed in what follows.

## Results

**Definition of shear transformation events.** As we saw in Fig. 2a, the local shear $\gamma_T(x, \varepsilon_M)$ exhibits important fluctuations, alternating activity, and quiescent phases. The concept of shear transformation zone has been introduced to deal with the flow of disordered materials: spatial zones reorganize, creating a local shear. We define such zones from the light scattering data. For this, we apply a threshold to the image $\gamma_m(\mathbf{r}, \varepsilon_M)$, and we use a particle detection algorithm to obtain individual shear events (see Supplementary Methods 3 and Supplementary Fig. 5 for details about threshold and detection algorithm). Events are numbered, and each event $i$ is associated with a macroscopic deformation $\varepsilon_{M_i}$ at which the event occurs, a position $\mathbf{r}_i$ (defined as the barycenter of $\gamma_m$), a surface $\Sigma_i$ on the image, and a mean microscopic shear $\gamma_i$. We also define a “seismic moment” as $M_i = \mu u_i L_i^2$ with $\mu$ the shear modulus of the granular assembly, $u_i$ the shear displacement, and $L_i^2$ the shear surface (see Fig. 1e). For a shear zone of width $w$, we have $\Sigma_i = wL_i$ and $u_i = w\gamma_i$, and thus $M_i = \mu\gamma_i\Sigma_i/w$. The shear modulus of the granular assembly may be estimated from mean-field theory of granular elasticity (equation (14) of ref. 36) and is $\mu = 200$ MPa for the pressure of 30 kPa. The shear bandwidth $w = 22d$ is measured from the mean strain (Fig. 1c). The moment magnitude is defined as $m_{w,i} = \frac{2}{3} \log_{10}(M_i) - 6.07$, with $M_i$ expressed in³⁷ Nm. The energy dissipated during an event is $E_i = \tau u_i L_i^2$, where $\tau$ is the shear stress along the shear band and is then directly linked to the moment: $M_i/\mu = E_i/\tau$. The value of $\tau = 33$ kPa may be obtained from the principal stresses $\sigma_1$ and $\sigma_3$ at failure using Mohr–Coulomb construction.

**Scaling laws of events.** We now look at the statistical characterization of the events. We consider for this the sequence of events occurring on one half shear band for macroscopic deformation $6\% \le \varepsilon_M \le 10\%$ (gray zone on Fig. 1b). The total number of counted events is $N_{tot} \approx 1.1 \times 10^5$. The minimum moment is dependent on the threshold and is $M_{min} = 6 \times 10^{-7}$ Nm, whereas the largest events have $M_{max} = 0.05$ Nm. The probability density function of energy $dN/dM$ is plotted on Fig. 3a, and decays as a power of energy $dN/dM \sim M^{-\beta}$, with $\beta \approx 2.1$. Although our moments stand in a range roughly 20 orders of magnitude below that of earthquakes, this behavior is similar to the empirical Gutenberg–Richter’s law. Indeed, the number $N(M)$ of earthquakes with a moment magnitude larger than $m_w$ is $\log_{10}(N(m_w)) = a - bm_w$, leading to $dN/dM \sim M^{-(1+(2/3)b)}$. The value of $b$ for faults is usually $\approx 1.0^1$, leading to a slope $\beta \approx 1.66$. It may be verified (see Supplementary Methods 4 and Supplementary Fig. 6) that the distributions remain the same during the whole post-failure regime, and that they do not depend on the occurrence of large events. In this sense, our model gouge is in a stationary state. The
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Scaling law of events. **a** Probability density $dN/dM$ of events of moment $M$. Dotted line is a power-law $-M^{-2.1}$. **b** Mean deformation $\langle y \rangle$; as a function of the moment. The dotted line is a power-law $-M^{0.17}$. **c** Ratio between the energy dissipated in events of moment $\ge M$ and the total dissipated energy.</figcaption></figure>

mean deformation of an event of moment $M$ is defined as
$\langle y \rangle = \frac{\sum_{M+dM>M_i>M} \gamma_i}{\sum_{M+dM>M_i>M}}$, for a small $dM$. As shown on Fig. 3b, this quantity is relatively constant $\langle y \rangle \sim M^{0.17}$. The broad distribution of the values of the moments $M_i$ (or equivalently of the relaxed energies $E_i$) is then mainly due to a broad distribution of sizes $L_i$, but not of deformation $\gamma_i$. In other words, the stress drop in an event $\Delta\tau_i = \mu\gamma_i$ has always the same order of magnitude. This is consistent with what is generally considered for earthquakes. Indeed, compilations of earthquake data argue for a scaling$^{38}$ $M \sim L^3$, while $M = \mu L^2 u = \mu\gamma L^3 = \Delta\tau L^3$, hence implying a constant stress drop. In our experiment, $\Delta\tau_i/\mu = \gamma_i \approx 10^{-4}$, this ratio is relatively close to the one typically obtained for earthquakes where$^{39}$ $\Delta\tau/\mu \approx 3 \times 10^{-5}$. We should however mention that, in our case, the increase of $\langle y \rangle$ i.e., of $\Delta\tau$, with the seismic moment, though weak, is significant (Fig. 3b). In the case of earthquakes, a potential similar scaling would be indiscernible, owing to the large uncertainty on the estimation of the average slip and the variety of geophysical contexts. We may also consider the relative fraction of the shear stress which is relaxed during an event: $\Delta\tau_i/\tau = \mu\gamma_i/\tau$. For large events, $\langle y \rangle \approx 2.10^{-4}$, whereas $\langle y \rangle \approx 3.10^{-5}$ for small events, hence showing that the events relax typically $\sim 0.1-1$ of the mean stress.

We finally look at the ratio of energy dissipated by events of moment greater than $M: E_{ev}(M) = \frac{\sum_{M_i \ge M} E_i}{\sum_{M_i \ge M} E_i}$, compared to the total dissipated energy $E_{dis}$ (see Supplementary Methods 5 for the details). Figure 3c shows that typically half of the energy is dissipated in events of moment $M > 10^{-4}$ Nm, while all detected events account for a value $R_s = 90\%$ of $E_{dis}$. $R_s$ is an equivalent of the seismic coupling defined in tectonics$^1$. The fact that the value is close to 1 means that only a minor part of the energy is dissipated into a continuous motion, distinct from the detected intermittent shear transformation events. For tectonic faults, the coupling can vary considerably with the geophysical context but is generally strong for interplate continental faults$^{40}$.

**Temporal organization of events.** The statistical laws governing the succession of shear transformations may be analyzed within the framework of the statistical laws of natural earthquakes.

We first look at the rate of events occurring at the same position as a particular event (so-called mainshock). Figure 4a shows the rate of events occurring after (“aftershocks”) of before (“foreshocks”) “mainshocks” of magnitude $m_w \ge -9$ (total number of mainshocks $\approx 13 \times 10^3$). Only aftershocks or foreshock events occurring at the same position ($\pm 15d$) are counted. For large delays, the rate of events is constant, corresponding to a background $(dN/d\epsilon)_{bg}$ rate of events uncorrelated to the mainshocks, while, at small delays, the rate excesses the background rate. The excess rate of aftershocks $(dN/d\epsilon)_{exc} = (dN/d\epsilon) - (dN/d\epsilon)_{bg}$ decays with the macroscopic deformation as $(dN/d\epsilon)_{exc} \sim \delta\epsilon_M^{-1}$ (see Fig. 4b). This behavior is reminiscent of Omori's law$^5$ which states that the rate of seismic events occurring after a mainshock decays with the time $t$ to the mainshock as $n(t) \sim t^{-1}$. Note however that such analysis amounts to blindly characterize time correlations between events. In particular, unlike what is generally done in earthquake analysis, the magnitude of the “mainshock” is not prescribed to be larger than that of its “aftershocks”. Such correlation analysis is a signature, but not a formal quantification, of causal triggering, which is explored in details later. An opposite evolution, reminiscent of an inverse Omori's law, characterizes, on average, an increasing rate of foreshocks before mainshocks, though with a smaller rate (Fig. 4a) that expresses an asymmetry of time clustering.

The productivity law describes the number of excess events in response to an event of magnitude $m_w$. For this, we integrate the total number of aftershocks in excess to the background: $N_{exc} = \int_{\delta\epsilon_M}^{\infty} (dN/d\epsilon)_{exc} d(\delta\epsilon_M)$. Figure 4c shows the evolution of $N_{exc}$ with the mainshock magnitude $m_w$, and we find that $N_{exc} \propto 10^{\alpha.m_w}$, with $\alpha \approx 0.44$. This result may be compared with the productivity law for natural earthquakes, where the number of aftershocks $n_{AS} \propto 10^{\alpha.m_w}$, with $\alpha$ in the range $0.6-1.2^{41}$. In striking contrast, the number of foreshocks appears independent of the mainshock magnitude (Fig. 4c). This is in full agreement with a previous analysis of seismic foreshocks showing that such precursory activity before any event is a mere statistical consequence of cascades of triggering$^{42}$. The triggering of deformation events in our system is thoroughly analyzed below.
<figure><img src="image_6.png" /><figcaption>Fig. 4 Temporal organization of events. a Rate of events $dN/de_M$ occurring after (●) and before (○) an event of magnitude $m_w \ge -9$ as a function of the deformation increment $\delta\varepsilon_M$. The dotted line is the background rate. b Rate of aftershock (●) or foreshock (○) events in excess to the background level as a function of $\delta\varepsilon_M$. The line is a $\delta\varepsilon_M^{-1}$ decay. c Number of events occurring in excess to background activity after (●) or before (○) a main-shock event of magnitude $m_w$. The line is a power law $N_{exc} \sim 10^{0.4m_w}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5 Distribution of inter-occurrence times. a Distribution of inter-occurrence deformation between successive events, for different magnitudes. b Distribution of re-scaled inter-occurrence deformation $\delta\varepsilon_M/\delta\varepsilon_M^*$. The plain line is the gamma distribution $P(x) \propto x^{q-1} \exp(-x/B)$ with $q = 0.6$ and $B = 1.9$.</figcaption></figure>

The distribution of events during the loading may be further characterized by considered the first-return deformation probability $P(\delta\varepsilon_M)$ which is analogous to the first return time probability for earthquakes. For this, we measure the macroscopic deformation $\delta\varepsilon_M$ between successive events, occurring at the same position ($\pm 15d$) along the band. Only events of magnitude greater than $m_w$ are considered. Figure 5a shows the distribution of inter-occurrence deformations for different moment thresholds. The distributions decay with a power-law, followed by an exponential decay. As shown on Fig. 5b, those distributions may be properly collapsed by considering for every magnitude $m_w$, the normalized deformation of $x = \delta\varepsilon_M/\overline{\delta\varepsilon_M}$, where $\overline{\delta\varepsilon_M}$ is the mean deformation between successive events. As indicated in Fig. 5b, the distribution may be well approximated by a Gamma distribution

$$P(x) \propto x^{q-1} \exp(-x/B) \quad (4)$$

with $q = 0.6$ and $B = 1.9$. $P(x)$ decays as a power-law with exponent $\approx 0.4$ up to values of $x \approx 1$, then exponential decays take place. This

behavior is not surprising, as it has been shown to be a mere consequence of triggering dynamics characterized by a GR distribution, the Omori's, and productivity laws$^{43}$. Our results are very similar to observations for tectonic seismicity where $q \approx 0.67^{44}$, or micro-seismicity $q \approx 0.74^{45}$. Such scaling laws are also observed in fracture experiments$^{46}$.

## Discussion

**Shear band viewed as a minimal model of the gouge.** Starting from an initially homogeneous assembly of beads, our system organizes spontaneously to reach a stationary regime where all the deformation is concentrated along shear planes. The analysis of the statistical properties of the strain fluctuations along those planes shows a strong quantitative and qualitative analogy with the statistical characteristic of natural earthquakes along a fault: the shear band may be viewed as a simplified gouge. We discuss here why this scale-free organization of deformation along a gouge is not observed in other laboratory systems.
<table><caption>Table 1 Quantitative comparison between natural gouges and granular shear band.</caption><thead><tr><th>Property</th><th>Earthquakes</th><th>Our experiment</th></tr></thead><tbody><tr><td>Temporal correlation function: C(δx = 0, δt) - δt<sup>-θ</sup></td><td>θ ≈ 0.70</td><td><sup>a</sup>θ ≈ 0.74</td></tr><tr><td>Moment or energy distribution: dN/dM - M<sup>-β</sup></td><td>β ≈ 1.7</td><td>β ≈ 2.1</td></tr><tr><td>Aftershocks rate: dN<sub>as</sub>/dt - t<sup>-P</sup></td><td>p ≈ 1.0</td><td><sup>a</sup>p ≈ 1.0<sup>b</sup>; 1.7<sup>c</sup></td></tr><tr><td>Productivity law: N<sub>as</sub> ~ 10<sup>am<sub>w</sub></sup></td><td>α ≈ 0.8</td><td><sup>a</sup>α ≈ 0.4<sup>b</sup>; 0.24<sup>c</sup></td></tr><tr><td>Recurrence time distribution: P(x) ~ x<sup>q-1</sup> exp(-x/B)</td><td>q ≈ 0.7</td><td>q ≈ 0.6</td></tr><tr><td>Stress drop/shear modulus: Δτ/μ</td><td>≈3.10<sup>-5</sup></td><td>≈10<sup>-4</sup></td></tr><tr><td>Branching ratio</td><td>0.8-1<sup>53</sup></td><td>≤1</td></tr></tbody></table>

Many mechanical systems other than crustal faults exhibit crackling noise when plastically deformed<sup>10,11,15,27,46–48</sup>. In particular, the statistical properties of deformation or mechanical stress fluctuations follow power laws which reveal the absence of any particular scale in the system at least in an extended inertial (scale) range. Those fluctuations arise from individual deformations which interact elastically to create complex scale-free dynamics. In the case of plasticity of amorphous materials before the macroscopic yielding<sup>15,16,48,49</sup>, the criticality of the system is related to the approach of the yield point. In this case, the plastic events are expected to be initially randomly distributed throughout the bulk of the sample and not localized in structures analog to natural faults. The brittle failure of an amorphous material is another configuration where crackling noise is observed<sup>10</sup>. In this case, the plasticity occurs in a damaged zone close to the propagating crack tip, and stationary plastic deformation cannot be defined. In some experiments such as compression of disordered materials<sup>46</sup> or some granular experiments<sup>13,27</sup>, the spatial extension of events and their localization is unknown.

In order to obtain a zone of intense plasticity in a stationary regime, one can shear an artificial gouge, made of a granular material confined between elastic plates of very different elastic modulus (much stiffer or much softer)<sup>14,25,26</sup>. At first glance, this appears to be a reasonable realization of a natural gouge which consists of highly crushed rocks confined between rigid material. However, in those cases, avalanches of various sizes are not observed, but instead, the dynamics are dominated by macro-slips involving all the sliding interface. De Geus et al.<sup>28</sup> have recently proposed a numerical model of frictional contact consisting of an amorphous layer confined between two elastic blocks, in which scale-free dynamics and large macro-slips events implying the whole interface coexists. Such a competition between an avalanche regime and a periodic stick-slip is reminiscent of the Parkfield segment of the San Andreas Fault, confined between a creeping zone and an unloaded segment, where large and pseudoperiodic earthquakes have been observed<sup>50</sup>. In our experiment, the dynamics of the model gouge are dominated by scale-free avalanches and we do not observe such macro-slips. This difference of behavior may arise from the difference in the confinement of the gouge. When the gouge is confined between elastic plates, there is an important contrast of mechanical properties between the gouge (which is an elasto-plastic granular material), and the plates (which are perfectly elastic). We may then expect that the plates transmit integrally the mechanical stress over all the gouge, leading to macroscopic slip events. In our experiments, materials that compose the fault and the surrounding medium are the same: both consist of the same glass beads. Given the applied pressure in our experiment, we do not expect any bead crushing, and this is in agreement with optical observations. So, the mechanical properties of the material are

probably very close inside and outside the shear band, and they are both elasto-plastic. So, the material outside the gouge does not behave as a rigid block transmitting the mechanical stress on all the interfaces. This is probably why we do not observe any large macro-slips but only localized avalanches displaying scale-free dynamics.

In summary, since the stationary shear band emerges from bulk material, we are able to observe a scale-free stationary dynamics occurring in a confined space. The shear band of granular material has the right dimensionality (the 2D shear plane in 3D space) and the right mechanical properties to accurately model a complex gouge at the laboratory scale. As a consequence, we directly observe shear events distributed along the shear plane. The statistical properties of those events are summarized in Table 1, and their size distributions, temporal and spatial organizations, as well as correlations of displacement, are very similar to the ones observed for natural earthquakes. We demonstrate below that the analogy can be pursued one step further through a thorough analysis of triggering.

**Triggering of deformation events.** Correlations in the deformation field and among deformation discrete events are found both in time and in space and obey power-law regimes that highlight the scale invariance of the system. However, correlation is distinct from causality, which in the present context is equivalent to triggering, i.e., how the occurrence of a deformation event mechanically triggers subsequent deformation events. The underlying causal structure can be inferred from the data using methods that have been developed in seismology<sup>30,31,51</sup> or in social science<sup>52</sup>. We find that triggering obeys a scale-free productivity law, so that the number *N* of directly triggered events, per mainshock, depends on the magnitude as *N* ~ 10<sup>0.24m<sub>w</sub></sup> (Fig. 6a), along with an Omori-like kernel, albeit with a relatively steep decay, the density of triggering events decreasing with time *t* after the trigger as *t*<sup>-P</sup> with *P* in the 1.6–1.8 interval, cf. Fig. 6b. Departure from these power laws is observed for the biggest events, that produce relatively more aftershocks in the early times, but are then followed by a clear activity shutdown, both features being likely due to a finite size effect and exhaustion of the stressed, ready-to-fail patches along with the deformation band after such large events. It is customary in the framework of these models to define a so-called branching ratio, which measures the capacity of a perturbation to sustain itself over potentially an infinite time (if the branching ratio is close to 1) or instead to die off quickly (if it is close to 0). This ratio can be estimated as the number of directly triggered aftershocks per “mainshock” averaged over all the events of the catalog. We here find that the branching ratio is very close to 1, implying that the background (i.e., non-triggered) activity consists of a few % at most so that most of the activity is made of events triggered by preceding

<footnote><sup>a</sup>With the substitution $t \to \epsilon_M$.</footnote>

<footnote><sup>b</sup>Count of events.</footnote>

<footnote><sup>c</sup>Triggering kernel.</footnote>
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6** Triggering kernel. **a** Productivity law giving the mean number *N* of triggered events for a trigger of magnitude *m*<sub>w</sub>. The best power-law fit in *N* ~ 10<sup>0.24m<sub>w</sub></sup> obtained when discarding the last point (biggest events) is shown in magenta. **b** Triggering kernels in time, conditioned on the moment of the trigger. We consider the same 8 "classes" of seismic moments as in (a). Power law decays in t<sup>-p</sup>, with p = 1 and p = 2, are shown for visual guidance. **c** Correlation in time, as in Fig. 4b, for two instances of synthetic datasets and for the real dataset.</figcaption></figure>

events, highlighting the dominance of triggering and therefore
clustering in the dynamics of deformation events. This is also
fully consistent with earthquakes, for which the branching ratio
ranges between 0.8 and 1<sup>53</sup>.

The productivity in 10<sup>0.24m<sub>w</sub></sup> found here is distinct from the
10<sup>0.40m<sub>w</sub></sup> scaling observed when stacking the activity past all m<sub>w</sub> ≥
9 events as in Fig. 4c; this is due to the fact that, in the latter case,
the stacking mixes causally-triggered sequences (e.g., if A triggers
B that triggers C, then both B and C will show up in the counting
of subsequent activity, while in Fig. 6 B is counted for A while C is
counted for B). We thus checked that this mixing does indeed re-
create the observations of Fig. 4. To do so, we exploit the fact that
the causal structure can be formulated as a linear model, that is
simply amenable to simulations<sup>30,31,51</sup>. We thus simulate
synthetic datasets of deformation events based on this model
and its basic ingredients: (1) seismic moments are independent,
identically distributed, and follow the Gutenberg-Richter-like
marginal distribution of the experiment; (2) a small proportion
(about 5%) of the events occur randomly in space and time, and
correspond to "spontaneous" events, i.e., events that are not
triggered by previous events; (3) the 95% other events are
triggered events from previous "mainshocks", their distribution
relative to the time and position of the mainshock following the
kernels observed for the experiment dataset (e.g., the temporal
kernel of Fig. 6b). Generating such synthetic datasets, we find that
the correlations (i.e., stacked rates) seen for the real data are
indeed well recovered (they are within the natural fluctuations of
simulation outcomes), demonstrating that these correlations
effectively emerge from more fundamental triggering kernels, cf.
Fig. 6c.

**Structural vs. temporal memory effects.** Strain correlation functions (Fig. 2), rate excess after main-shocks (Fig. 4a), and the causal triggering kernel (Fig. 6) all indicate the existence of scale-free memory effects in our system. In the context of seismology, memory effects are remarkably revealed by the existence of aftershock sequences which are quantified by Omori's law stating that the number of aftershocks decays as the inverse of the time elapsed since the mainshock¹. This law assumes implicitly that the time is the physical variable that governs the memory. The origin of such temporal dependence is however unclear. Several mechanisms such as the temporal dependence of microscopic friction law⁵⁴⁻⁵⁷, sub-critical crack growth, the occurrence of afterslip⁵⁸, or poro-elasticity and the evolution with time of pore fluid pressure⁵⁹ have been proposed as a possible sources of temporal memory effects controlling earthquake occurrences. However, the direct links between any time-dependent microscopic mechanism and memory effects in seismicity are still debated.

In our experiment, we can test whether time is indeed the right
parameter to describe the memory effect that we observe. For this,
we performed experiments at different macroscopic deformation
rates. Figure 7a shows the normalized correlation functions of
the microscopic deformation expressed as a function of the time
increment. If every experiment shows a memory effect, the
magnitude of memory depends on the strain velocity. At a fixed
time delay δt, the correlation function decreases with the velocity.
This reveals that time does not seem to be the right parameter to
describe memory. This may be evidenced by plotting the
correlation functions as a function of the macroscopic strain
increment δε<sub>M</sub> Fig. 7b. In this plot, the curves collapse,
demonstrating that the correlation function decays with the
strain increment and not with the time increment. This
independence on the shear rate may also be evidenced by
considering the similarity of the size distributions of events on
Fig. 7c.

This laboratory observation is evidently not in contradiction
with Omori's law. Indeed, the driving velocity of a given fault is
constant on the temporal scale of human observations. So,
<figure><img src="image_11.png" /><figcaption>**Fig. 7** Effect of deformation rate. Comparison of three experiments performed at different compression velocities: red crosses $\dot{\varepsilon}_M = 3.5 \times 10^{-6} \text{ s}^{-1}$; green-filled triangles $\dot{\varepsilon}_M = 1.1 \times 10^{-5} \text{ s}^{-1}$; open black circles: $\dot{\varepsilon}_M = 3.5 \times 10^{-5} \text{ s}^{-1}$. **a** The normalized correlation function of the deformation (3) as a function of the time increment $\delta t$. **b** Same data plotted as a function of the strain increment $\delta e_M$. **c** Probability density $dN/dM$ of events of moment M. (strain increment $\delta e_M^* = 1.5 \times 10^{-6}$, event threshold (see Supplementary Methods 3) $\gamma_s = 0.1 \gamma_0$).</figcaption></figure>

describing memory effects in terms of time increment or in terms of strain increment is then equivalent for natural faults. The fact that the memory is strain-dependent rather than time-dependent in our system suggests that the memory could be linked to structural/topological rearrangements within the granular medium, inducing a redistribution of local stress and possibly triggering slip events. It is then not surprising that the macroscopic deformation may be the parameter that governs the plasticity around a given position in the material. This also raises the question of the potential role of such geometrical rearrangements in the “time” correlations characterizing natural seismicity. In that case, such mechanisms could combine with truly time-dependent, thermally activated processes such as sub-critical crack growth, to explain memory effects in earthquake occurrences. Interestingly, slip-dependent and time-dependent memory effects combine as well in the classical Rate-and-State friction laws⁶⁰,⁶¹ that remains nowadays a classical framework of interpretation of earthquake physics⁵⁵,⁶².

## Conclusion

By looking at the intermittent strain fluctuations, we showed that a shear band inside an athermal disordered material is an analog of a natural fault: the deformation consists of many micro-slip occurring along a plane, and their collective dynamics is characterized by statistical properties remarkably consistent with the

empirical laws of seismology. This analogy with natural faults is obtained when the fluctuations are observed after macroscopic yielding of the granular medium, when a steady-state regime takes place. The laboratory and natural fluctuations observed are then characteristic of a critical behavior after yielding, which is presumably different from the stress-tuned criticality observed for many systems before yielding. Our statistical analysis of micro-slips also quantifies the causal triggering between events, and reveals that this underlying triggering mechanism is at the root of the space-time correlations in the dynamics, as it has been previously shown for natural earthquakes.

Despite its simplicity, our experimental model recreates the whole set of spatio-temporal characteristics of earthquakes: moment distribution, recurrence time distribution, productivity law, decay of aftershocks rate, and triggering kernel. This simply suggests that those laws may be reproduced by using simple numerical or theoretical systems of frictional particles. Hence, the analysis of such systems should allow understanding the organization of the microscopic stress field close to the shear band. Experimentally, the possibility to control a model fault in the laboratory also opens the road to many studies, such as the effect of mechanical noise on the size distribution of slip events, the influence of the elastic properties of the surrounding material with respect to that of the band, or a study of size effects.
<figure><img src="image_12.png" /><figcaption>**Fig. 8** Geometry of the shear cell. **a** Schematic drawing of the geometry. The sample is the gray box, green and red planes are shear planes, and blue arrows show the displacement directions of the blocks of material with respect to the bottom block. **b** Orientations of the shear. Blue arrows represent the directions of displacement with respect to the shear bands.</figcaption></figure>

## Methods

The experimental setup consists of a biaxial compressive test in-plane strain conditions already described<sup>63</sup>. The granular material is composed of glass beads of diameter $d = 70–110 \ \mu\text{m}$ and initial volume fraction ≈0.60. The granular material is enclosed between a latex membrane (85 × 55 × 25 mm<sup>3</sup>) and a glass plate. A pump produces a partial vacuum inside the membrane, creating confining stress $-\sigma_3 = 30 \ \text{kPa}$.

The sample is placed in the biaxial apparatus where a compression is imposed in the vertical direction (noted 1 on Fig. 8a) while the lateral sides (noted 3 on Fig. 8a) are maintained under a constant stress $\sigma_3$. Two fixed planes ensure that there is no deformation along the direction 2. At the top, a moving plate exerts a compression of the sample at fixed velocity and the bottom plate is fixed. The velocity of the motor is of the order $v \sim \mu\text{s}^{-1}$ (different velocities are used, see main text for precise values), leading to a deformation rate $\epsilon_M = v/L \sim 10^{-5} \ \text{s}^{-1}$ where $L = 85 \ \text{mm}$. The origin of deformation $\epsilon_M = 0$ is taken at the beginning of the contact between top plate and sample, but its precise value is without important because analysis considers only deformation increments. The stress applied at the top of the sample is defined as $-\sigma_1 = -\sigma_3 + F/S$, where F is the force measured by a sensor fixed to the upper plate, and $S = 55 \times 25 \ \text{mm}^2$ is the section of the sample that we consider as constant during all the experiment. When the stress plateau is attained, one or two shear bands form into the sample (green and red planes on Fig. 8a).

Deformations are observed through the front glass plate using DWS, a method already described elsewhere<sup>35,64</sup>. A laser beam at 532 nm is expanded to illuminate the entire sample. The light undergoes multiple scattering inside the granular material and we collect the backscattered rays. The latter interfere and form a speckle pattern. The image of the front side of the sample is recorded by a 2352 × 1728 pixels camera. Speckle images are then subdivided into square zones of size 16 × 16 pixels that we call metapixel. The correlation between two images 1 and 2 is then calculated for each metapixel, and a map of correlation of 147 × 108 metapixels is obtained. The size of the metapixel corresponds to $6.0d \times 6.0d$ on the sample.

Average $\langle \cdot \rangle$ in Eq. (3) is computed in the following way. The total length of the spatiotemporal series is of $7 \times 10^4 \delta\epsilon_M^*$. To smooth out the fluctuations of the strain rate at the level of the studied band, correlation functions are computed by slices of length $3 \times 10^3 \delta\epsilon_M^*$ in deformation interval and on position. Then the resulting functions are averaged on the whole duration of the series.

## Data availability

Optical correlations functions (raw 8-bits images) at https://doi.org/10.5281/zenodo.4525025.

Received: 25 August 2020; Accepted: 12 March 2021;

Published online: 14 May 2021

## References

1. Scholz, C. H. *The Mechanics of Earthquakes and Faulting* (Cambridge University Press, 2019).

2. Mai, P. M. & Beroza, G. C. A spatial random field model to characterize complexity in earthquake slip. *J. Geophys. Res.* **107**, ESE 10-1–ESE 10-21 (2002).

3. Kagan, Y. Y. & Jackson, D. D. Long-term earthquake clustering. *Geophys. J. Int.* **104**, 117–133 (1991).

4. Kagan, Y. Y. Earthquake spatial distribution: the correlation dimension. *Geophys. J. Int.* **168**, 1175–1194 (2007).

5. Omori, F. On the aftershocks of earthquakes. *J. Coll. Sci.* **7**, 111–200 (1894).

6. Tajima, F. & Kanamori, H. Global survey of aftershock area expansion patterns. *Phys. Earth Planet. Inter.* **40**, 77–134 (1985).

7. Marsan, D., Bean, C. J., Steacy, S. & McCloskey, J. Observation of diffusion processes in earthquake populations and implications for the predictability of seismicity systems. *J. Geophys. Res.* **105**, 28081–28094 (2000).

8. Sethna, J. P., Dahmen, K. A. & Myers, C. R. Crackling noise. *Nature* **410**, 242–250 (2001).

9. Salje, E. K. & Dahmen, K. A. Crackling noise in disordered materials. *Annu. Rev. Condens. Matter Phys.* **5**, 233–254 (2014).

10. Barés, J., Dubois, A., Hattali, L., Dalmas, D. & Bonamy, D. Aftershock sequences and seismic-like organization of acoustic events produced by a single propagating crack. *Nat. Commun.* **9**, 1253 (2018).

11. Baró, J. et al. Statistical similarity between the compression of a porous material and earthquakes. *Phys. Rev. Lett.* **110**, 088702 (2013).

12. Vu, C.-C., Amitrano, D., Plé, O. & Weiss, J. Compressive failure as a critical transition: experimental evidence and mapping onto the universality class of depinning. *Phys. Rev. Lett.* **122**, 015502 (2019).

13. Zadeh, A. A., Barés, J., Socolar, J. E. & Behringer, R. P. Seismicity in sheared granular matter. *Phys. Rev. E* **99**, 052902 (2019).

14. Rivière, J., Lv, Z., Johnson, P. & Marone, C. Evolution of b-value during the seismic cycle: Insights from laboratory experiments on simulated faults. *Earth Planet. Sci. Lett.* **482**, 407 – 413 (2018).

15. Denisov, D., Lörincz, K., Uhl, J., Dahmen, K. A. & Schall, P. Universality of slip avalanches in flowing granular matter. *Nat. Commun.* **7**, 10641 (2016).

16. Uhl, J. T. et al. Universal quake statistics: from compressed nanocrystals to earthquakes. *Sci. Rep.* **5**, 16493 (2015).

17. Dahmen, K. A., Ben-Zion, Y. & Uhl, J. T. A simple analytic theory for the statistics of avalanches in sheared granular materials. *Nat. Phys.* **7**, 554 (2011).

18. Schorlemmer, D., Wiemer, S. & Wyss, M. Variations in earthquake-size distribution across different stress regimes. *Nature* **437**, 539–542 (2005).

19. Kagan, Y. Y. Seismic moment distribution revisited: II. Moment conservation principle. *Geophys. J. Int.* **149**, 731–754 (2002).

20. Ben-Zion, Y., Dahmen, K. A. & Uhl, J. T. A unifying phase diagram for the dynamics of sheared solids and granular materials. *Pure Appl. Geophys.* **168**, 2221–2237 (2011).

21. Jagla, E. A., Landes, F. P. & Rosso, A. Viscoelastic effects in avalanche dynamics: a key to earthquake statistics. *Phys. Rev. Lett.* **112**, 174301 (2014).

22. Le Bouil, A., Amon, A., McNamara, S. & Crassous, J. Emergence of cooperativity in plasticity of soft glassy materials. *Phys. Rev. Lett.* **112**, 246001 (2014).

23. Gimbert, F., Amitrano, D. & Weiss, J. Crossover from quasi-static to dense flow regime in compressed frictional granular media. *Europhys. Lett.* **104**, 46001 (2013).

24. Karimi, K., Amitrano, D. & Weiss, J. From plastic flow to brittle fracture: role of microscopic friction in amorphous solids. *Phys. Rev. E* **100**, 012908 (2019).
25. Geller, D. A., Ecke, R. E., Dahmen, K. A. & Backhaus, S. Stick-slip behavior in a continuum-granular experiment. *Phys. Rev. E* **92**, 060201 (2015).

26. Lieou, C. K. C. et al. Simulating stick-slip failure in a sheared granular layer using a physics-based constitutive model. *J. Geophys. Res.* **122**, 295–307 (2017).

27. Lherminier, S. et al. Continuously sheared granular matter reproduces in detail seismicity laws. *Phys. Rev. Lett.* **122**, 218501 (2019).

28. de Geus, T. W. J., Popović, M., Ji, W., Rosso, A. & Wyart, M. How collective asperity detachments nucleate slip at frictional interfaces. *Proc. Natl Acad. Sci. USA* **116**, 23977–23983 (2019).

29. Helmstetter, A., Kagan, Y. Y. & Jackson, D. D. Importance of small earthquakes for stress transfers and earthquake triggering. *J. Geophys. Res.* **110**, B05S08 (2005).

30. Marsan, D. & Lengline, O. Extending earthquakes' reach through cascading. *Science* **319**, 1076–1079 (2008).

31. Marsan, D. & Lengliné, O. A new estimation of the decay of aftershock density with distance to the mainshock. *J. Geophys. Res.* **115**, B09302 (2010).

32. Kagan, Y. *Earthquakes: Models, Statistics, Testable Forecasts* (Wiley-Blackwell, 2014).

33. Nguyen, T. B. & Amon, A. Experimental study of shear band formation: bifurcation and localization. *Europhys. Lett.* **116**, 28007 (2016).

34. Vardoulakis, I. & Sulem, J. *Bifurcation Analysis in Geomechanic* (CRC Press, 1995).

35. Amon, A., Mikhailovskaya, A. & Crassous, J. Spatially resolved measurements of micro-deformations in granular materials using diffusing wave spectroscopy. *Rev. Sci. Instrum.* **88**, 051804 (2017).

36. Makse, H. A., Gland, N., Johnson, D. L. & Schwartz, L. Granular packings: nonlinear elasticity, sound propagation, and collective relaxation dynamics. *Phys. Rev. E* **70**, 061302 (2004).

37. Kanamori, H. & Brodsky, E. E. The physics of earthquakes. *Rep. Progr. Phys.* **67**, 1429–1496 (2004).

38. Wells, D. L. & Coppersmith, K. J. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. *Bull. Seismol. Soc. Am.* **84**, 974–1002 (1994).

39. de Arcangelis, L., Godano, C., Grasso, J. R. & Lippiello, E. Statistical physics approach to earthquake occurrence and forecasting. *Phys. Rep.* **628**, 1–91 (2016).

40. Bird, P. & Kagan, Y. Plate-tectonic analysis of shallow seismicity: apparent boundary width, beta, corner magnitude, coupled lithosphere thickness, and coupling in seven tectonic settings. *Bull. Seismol. Soc. Am.* **94**, 2380–2399 (2004).

41. Hainzl, S. & Marsan, D. Dependence of the Omori–Utsu law parameters on main shock magnitude: Observations and modeling. *J. Geophys. Res.* **113**, B10309 (2008).

42. Helmstetter, A. & Sornette, D. Foreshocks explained by cascades of triggered seismicity. *J. Geophys. Res.* **108**, 2457 (2003).

43. Saichev, A. & Sornette, D. "Universal" distribution of interearthquake times explained. *Phys. Rev. Lett.* **97**, 078501 (2006).

44. Corral, A. Long-term clustering, scaling, and universality in the temporal occurrence of earthquakes. *Phys. Rev. Lett.* **92**, 108501 (2004).

45. Davidsen, J. & Kwiatek, G. Earthquake interevent time distribution for induced micro-, nano-, and picoseismicity. *Phys. Rev. Lett.* **110**, 068501 (2013).

46. Ribeiro, H. V. et al. Analogies between the cracking noise of ethanol-dampened charcoal and earthquakes. *Phys. Rev. Lett.* **115**, 025503 (2015).

47. Abed Zadeh, A., Barés, J., Socolar, J. E. S. & Behringer, R. P. Seismicity in sheared granular matter. *Phys. Rev. E* **99**, 052902 (2019).

48. Miguel, M.-C., Vespignani, A., Zapperi, S., Weiss, J. & Grasso, J.-R. Intermittent dislocation flow in viscoplastic deformation. *Nature* **410**, 667–671 (2001).

49. Lin, J., Gueudré, T., Rosso, A. & Wyart, M. Criticality in the approach to failure in amorphous solids. *Phys. Rev. Lett.* **115**, 168001 (2015).

50. Bakun, W. H. & Lindh, A. G. The Parkfield, California, earthquake prediction experiment. *Science* **229**, 619–624 (1985).

51. Lengliné, O., Frank, W., Marsan, D. & Ampuero, J.-P. Imbricated slip rate processes during slow slip transients imaged by low-frequency earthquakes. *Earth Planet. Sci. Lett.* **476**, 122–131 (2017).

52. Mohler, G. O., Short, M. B., Brantingham, P. J., Schoenberg, F. P. & Tita, G. E. Self-exciting point process modeling of crime. *J. Am. Stat. Assoc.* **106**, 100–108 (2011).

53. Nandan, S., Ouillon, G., Wiemer, S. & Sornette, D. Objective estimation of spatially variable parameters of epidemic type aftershock sequence model: application to California. *J. Geophys. Res.* **122**, 5118–5143 (2017).

54. Dieterich, J. H. Time-dependent friction as a possible mechanism for aftershocks. *J. Geophys. Res.* **77**, 3771–3781 (1972).

55. Scholz, C. H. Earthquakes and friction laws. *Nature* **391**, 37–42 (1998).

56. Gomberg, J., Beeler, N. & Blanpied, M. On rate-state and coulomb failure models. *J. Geophys. Res.* **105**, 7857–7871 (2000).

57. Karner, S. L. & Marone, C. The effect of shear load on frictional healing in simulated fault gouge. *Geophys. Res. Lett.* **25**, 4561–4564 (1998).

58. Perfettini, H. & Avouac, J.-P. Modeling afterslip and aftershocks following the 1992 landers earthquake. *J. Geophys. Res.* **112**, B07409 (2007).

59. Nur, A. & Booker, J. R. Aftershocks caused by pore fluid flow? *Science* **175**, 885–887 (1972).

60. Dietrich, J. H. Modeling of rock friction: 1. Experimental results and constitutive equations. *J. Geophys. Res.* **84**, 2161–2168 (1979).

61. Baumberger, T. & Caroli, C. Solid friction from stick-slip down to pinning and aging. *Adv. Phys.* **55**, 279–348 (2006).

62. Helmstetter, A. & Shaw, B.E. Afterslip and aftershocks in the rate-and-state friction law. *J. Geophys. Res.* **114**, B01308 (2009).

63. Le Bouil, A. et al. A biaxial apparatus for the study of heterogeneous and intermittent strains in granular materials. *Granul. Matter* **16**, 1–8 (2014).

64. Erpelding, M., Amon, A. & Crassous, J. Diffusive wave spectroscopy applied to the spatially resolved deformation of a solid. *Phys. Rev. E* **78**, 046104 (2008).

## Acknowledgements

The authors thank David Amitrano, Anaël Lemaitre, Damien Vandembroucq for scientific discussions, and Patrick Chasle for technical support. We acknowledge the funding from Agence Nationale de la Recherche Grant no. ANR-16-CE30-0022.

## Author contributions

J.C., A.A., and J.W. designed research, D.H., A.A., D.M., J.W., and J.C. performed research, J.C., A.A., and D.M. analyzed data, and J.C., A.A., D.M., and J.W. wrote the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-021-00147-1.

**Correspondence** and requests for materials should be addressed to J.C.

**Peer review information** Primary handling editor: Joseph Aslin

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_13.png" /></figure>

**Open Access** This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021