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

Abstract

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.

Introduction

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^6,7. 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^8,9 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^11,12, or systems involving frictional sliding on or within a granular media^13,14. 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^16,17. 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^22,23,24, 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^14,25,26,27. 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^29,30,31,32.

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 space-time correlations in the dynamics, much like it does for earthquakes. We argue on this basis that a frictional shear band in a 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 σ₁ 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 ε_M, and the applied axial stress σ₃ are measured (see “Methods”).

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. ε_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 Δ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 × L.

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.1{0}^{-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(ε_M, r), where r is the position on the observation plane, as shown in Fig. 1d. High correlation g_I ≈ 1 (white pixels) indicates that beads are uniformly translated without relative motions, whereas low correlation g_I < < 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 σ₃ − σ₁ (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 ε_M of the material exceeds a yield strain ε_Y ≈ 4.5%, the stress difference σ₃ − σ₁ 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 ε_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 r is mainly a shear γ_m(r, ε_M) = ∂u_x/∂z, where 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}({\bf{r}},{\varepsilon }_{M})=-{\gamma }_{0}{\mathrm{ln}}\,[\,{g}_{I}({\bf{r}},{\varepsilon }_{M})]$$

(1)

with γ₀ = 2.6 × 10⁻⁴ a constant given by the optical properties of the material. The time-averaged local deformation is defined as:

$${\bar{\gamma }}_{m}({\bf{r}},{\varepsilon }_{M})=(1/{{\Delta }}{\varepsilon }_{M})\mathop{\int}\nolimits_{{\varepsilon }_{M}}^{{\varepsilon }_{M}+{{\Delta }}{\varepsilon }_{M}}{\gamma }_{m}({\bf{r}},{{\varepsilon }_{M}}^{\prime} )d{{\varepsilon }_{M}}^{\prime}$$

(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 )=\mathop{\int}\nolimits_{-\infty }^{z}{\bar{\gamma }}_{m}(z^{\prime} ,{\varepsilon }_{M})dz^{\prime}$. Figure 1c shows the displacement field across the shear band, 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^* ~ 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.1{0}^{-7}$, the relative sliding of the two blocks Δu_x = u_x(∞) − u_x(−∞) if found to be Δu_x = 32 nm. This value is close to the imposed value Δ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}\mathop{\int}\nolimits_{-w}^{w}{\gamma }_{m}({\bf{r}},{\varepsilon }_{M})dz$. Figure 2b shows γ_T(x, ε_M) into the (x, ε_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\,$$

(3)

Fig. 2: Strain fluctuations and correlations.

a, b Granular experiment. a Spatiotemporal correlation function C(δx = 0, δε_M) as a function of δε_M see Eq. (3). b Spatiotemporal evolution of the local strain γ_T(x, ε_M) (see text) in a shear band. c, d Californian earthquake catalog. c Spatiotemporal correlation function C(δx = 0, δt) as a function of δ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).

where 〈⋅〉 is an average both on deformation and position (see “Methods” for details). Figure 2a is a plot of C(δx = 0, δε_M) as a function of the δε_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 θ = 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 γ_T(x, ε_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 γ_m(r, ε_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 r_i (defined as the barycenter of γ_m), a surface Σ_i on the image, and a mean microscopic shear γ_i. We also define a “seismic moment” as ${M}_{i}=\mu {u}_{i}{L}_{i}^{2}$ with μ 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 Σ_i = wL_i and u_i = wγ_i, and thus ${M}_{i}=\mu {\gamma }_{i}{{{\Sigma }}}_{i}^{2}/w$. The shear modulus of the granular assembly may be estimated from mean-field theory of granular elasticity (equation (14) of ref. ³⁶) and is μ = 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}{{\mathrm{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 τ is the shear stress along the shear band and is then directly linked to the moment: M_i/μ = E_i/τ. The value of τ = 33 kPa may be obtained from the principal stresses σ₁ and σ₃ 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% ≤ ε_M ≤ 10% (gray zone on Fig. 1b). The total number of counted events is N_tot ≈ 1.1 × 10⁵. The minimum moment is dependent on the threshold and is M_min = 6 × 10⁻⁷ 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 ~ M^−β, with β ≃ 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 ${{\mathrm{log}}\,}_{10}(N({m}_{w}))=a-b{m}_{w}$, leading to dN/dM ~ M^{−(1+(2/3)b)}. The value of b for faults is usually ≃1.0¹, leading to a slope β ≃ 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 mean deformation of an event of moment M is defined as $\langle\gamma \rangle={\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 〈γ〉 ~ 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 γ_i. In other words, the stress drop in an event Δτ_i = μγ_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³⁸ M ~ L³, while M = μL²u = μγL³ = ΔτL³, hence implying a constant stress drop. In our experiment, Δτ_i/μ = γ_i ≈ 10⁻⁴, this ratio is relatively close to the one typically obtained for earthquakes where³⁹ Δτ/μ ≈ 3 × 10⁻⁵. We should however mention that, in our case, the increase of 〈γ〉 i.e., of Δτ, 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: Δτ_i/τ = μγ_i/τ. For large events, 〈γ〉 ≈ 2.10⁻⁴, whereas 〈γ〉 ≈ 3.10⁻⁵ for small events, hence showing that the events relax typically ~0.1–1 of the mean stress.

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 〈γ〉; 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 ≥M and the total dissipated energy.

We finally look at the ratio of energy dissipated by events of moment greater than M: ${E}_{ev}(M)={\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⁻⁴ 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¹. 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⁴⁰.

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 ≥ −9 (total number of mainshocks ≈13 × 10³). Only aftershocks or foreshock events occurring at the same position (±15d) are counted. For large delays, the rate of events is constant, corresponding to a background (dN/dε)_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ε)_exc = (dN/dε) − (dN/dε)_bg decays with the macroscopic deformation as ${(dN/d\varepsilon )}_{exc} \sim \delta {\varepsilon }_{M}^{-1}$ (see Fig. 4b). This behavior is reminiscent of Omori’s law⁵ which states that the rate of seismic events occurring after a mainshock decays with the time t to the mainshock as n(t) ~ t⁻¹. 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.

Fig. 4: Temporal organization of events.

a Rate of events dN/dε_M occurring after (•) and before (∘) an event of magnitude m_w ≥ −9 as a function of the deformation increment δε_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 δε_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 1{0}^{0.4{m}_{w}}$.

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}=\mathop{\int}\nolimits_{\delta {\varepsilon }_{M}^{* }}^{\infty }{(dN/d\varepsilon )}_{exc}d(\delta {\varepsilon }_{M})$. Figure 4c shows the evolution of N_exc with the mainshock magnitude m_w, and we find that ${N}_{exc}\propto 1{0}^{\alpha .{m}_{w}}$, with α ≃ 0.44. This result may be compared with the productivity law for natural earthquakes, where the number of aftershocks ${n}_{AS}\propto 1{0}^{\alpha .{m}_{w}}$, with α in the range 0.6–1.2⁴¹. 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⁴². The triggering of deformation events in our system is thoroughly analyzed below.

The distribution of events during the loading may be further characterized by considered the first-return deformation probability P(δε_M) which is analogous to the first return time probability for earthquakes. For this, we measure the macroscopic deformation δε_M between successive events, occurring at the same position (±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)$$

(4)

with q = 0.6 and B = 1.9. P(x) decays as a power-law with exponent ≈0.4 up to values of x ≈ 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⁴³. Our results are very similar to observations for tectonic seismicity where q ≃ 0.67⁴⁴, or micro-seismicity q ≃ 0.74⁴⁵. Such scaling laws are also observed in fracture experiments⁴⁶.

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}/\overline{\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.

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.

Many mechanical systems other than crustal faults exhibit crackling noise when plastically deformed^{10,11,15,27,46,47,48}. 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^15,16,48,49, 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¹⁰. 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⁴⁶ or some granular experiments^13,27, 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)^14,25,26. 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.²⁸ 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 pseudo-periodic earthquakes have been observed⁵⁰. 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.

Table 1 Quantitative comparison between natural gouges and granular shear band.

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^30,31,51 or in social science⁵². 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 \sim 1{0}^{0.24{m}_{w}}$ (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^−p 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 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⁵³.

Fig. 6: Triggering kernel.

a Productivity law giving the mean number N of triggered events for a trigger of magnitude m_w. The best power-law fit in $N \sim 1{0}^{0.24{m}_{w}}$ 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^−p, 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.

The productivity in $1{0}^{0.24{m}_{w}}$ found here is distinct from the $1{0}^{0.40{m}_{w}}$ scaling observed when stacking the activity past all m_w ≥ 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^30,31,51. 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^54,55,56,57, 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 δε_M 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.

Fig. 7: Effect of deformation rate.

Comparison of three experiments performed at different compression velocities: red crosses ${\dot{\varepsilon }}_{M}=3.5\,{\times} 1{0}^{-6}$ s⁻¹; green-filled triangles ${\dot{\varepsilon }}_{M}=1.1\,\times 1{0}^{-5}$ s⁻¹; open black circles: ${\dot{\varepsilon }}_{M}=3.5\,\times 1{0}^{-5}$ s⁻¹. a The normalized correlation function of the deformation (3) as a function of the time increment δt. b Same data plotted as a function of the strain increment δε_M. c Probability density dN/dM of events of moment M. (strain increment $\delta {\varepsilon }_{M}^{* }=1.5\,\times 1{0}^{-6}$, event threshold (see Supplementary Methods 3) γ_s = 0.1 γ₀).

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, 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^60,61 that remains nowadays a classical framework of interpretation of earthquake physics^55,62.

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.

Methods

The experimental setup consists of a biaxial compressive test in-plane strain conditions already described⁶³. The granular material is composed of glass beads of diameter d = 70–110 μm and initial volume fraction ≈0.60. The granular material is enclosed between a latex membrane (85 × 55 × 25 mm³) and a glass plate. A pump produces a partial vacuum inside the membrane, creating confining stress −σ₃ = 30 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 σ₃. 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 ~ μm s⁻¹ (different velocities are used, see main text for precise values), leading to a deformation rate $\dot{{\epsilon }_{M}}=v/L \sim 1{0}^{-5}$ s⁻¹ where L = 85 mm. The origin of deformation ϵ_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 −σ₁ = −σ₃ + F/S, where F is the force measured by a sensor fixed to the upper plate, and S = 55 × 25 mm² 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).

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.

Deformations are observed through the front glass plate using DWS, a method already described elsewhere^35,64. 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 × 6.0d on the sample.

Average 〈⋅〉 in Eq. (3) is computed in the following way. The total length of the spatiotemporal series is of $7\times 1{0}^{4}\delta {\varepsilon }_{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 1{0}^{3}\delta {\varepsilon }_{M}^{* }$ in deformation interval and on position. Then the resulting functions are averaged on the whole duration of the series.