Continuously sheared granular matter reproduces in detail seismicity laws

TL;DR

This study demonstrates that a simple shear experiment on granular matter can replicate key seismicity laws, including Gutenberg-Richter, Omori, and inter-event time distributions, revealing insights into earthquake dynamics.

Contribution

The paper introduces a shear experiment with granular disks that quantitatively reproduces seismicity laws and features, highlighting the role of force networks in earthquake-like behavior.

Findings

Energy distribution follows Gutenberg-Richter law

Foreshocks and aftershocks follow Omori laws

Inter-event times match real earthquake distributions

Abstract

We introduce a shear experiment that quantitatively reproduces the main laws of seismicity. By continuously and slowly shearing a compressed monolayer of disks in a ring-like geometry, our system delivers events of frictional failures with energies following a Gutenberg-Richter law. Moreover foreshocks and aftershocks are described by Omori laws and inter-event times also follow exactly the same distribution as real earthquakes, showing the existence of memory of past events. Other features of real earthquakes qualitatively reproduced in our system are both the existence of a quiescence preceding mainshocks, as well as magnitude correlations linked to large quakes. The key ingredient of the dynamics is the nature of the force network, governing the distribution of frictional thresholds.

Figures

Full text

For more than a century, fracture and stick-slip frictional sliding have tried to explain the behavior of earthquakes. Brittle fracture induced by shear Reid (1911) was the most accepted model until the sixties. However, a more precise analysis of the radiated waves Benioff (1964), the low amount of stress released by an earthquake in relation to the available one, the high energies needed to shear over a fractured surface, and over all, the lack of healing required to generate a second earthquake at the same location and close in time to the first one, set stick-slip sliding mechanisms as a more plausible explanation of earthquakes Brace and Byerlee (1966). Despite these facts, the subcritical fracture of heterogeneous materials shows naturally a jerky behavior that seems closer to earthquake statistics than frictional sliding, which commonly displays a quasi-periodic stick-slip dynamics. Indeed, several fracture experiments Scholz (1968); Måløy et al. (2006); Baró et al. (2013); Stojanova et al. (2014); Barés et al. (2018) and numerical models Kun et al. (2013, 2014) have reported statistics of events following power-law distributions of sizes that have been compared to the the Gutenberg-Richter law Gutenberg and Richter (1956). The existence of aftershocks that follow the Omori law Omori (1894) are also common in fracture experiments Måløy et al. (2006); Baró et al. (2013); Stojanova et al. (2014); Barés et al. (2018).

Concerning stick-slip frictional sliding, different laboratory experiments have analyzed the sliding dynamics between two solid blocks. From a physical perspective, studies on acrylic blocks have focused on the complex evolution of the frictional strength during the slipping process, describing the behavior as a dynamic fracture problem Rubinstein et al. (2004); Ben-David et al. (2010). Recent friction experiments on rocks have reported results on supershear ruptures Passelègue et al. (2013) and precursory activity prior to stick-slip instabilities Passelègue et al. (2016). Precursory activity to stick-slip instabilities has been also reported in experiments shearing a layer of granular material Rivière et al. (2018). Other relevant results on similar experimental systems include remote triggering Johnson et al. (2008), and the controlled slowing down of the dynamics Scuderi et al. (2016). However, one common limitation of many of those laboratory experiments is the fact that they show a main dynamics consisting in a quasi-periodic stick-slip behavior with a narrow distribution of sizes, which do not correspond with the complex dynamics of real earthquakes described by the laws of seismicity. Other experimental systems have also sheared a granular layer, aiming at mimicking the intermittent behavior of a tectonic fault Daniels and Hayman (2008); Walker et al. (2014). Nevertheless it has been difficult to obtain a distribution of events that resembles the Gutenberg-Richter law Gutenberg and Richter (1956) due to insufficient statistical sets of data.

Here we introduce a shear experiment capable of quantitatively reproducing the main statistical laws describing seismicity (Gutenberg- Richter law Gutenberg and Richter (1956), Omori law Omori (1894), distribution of inter-event times Corral (2004)), as well as sharing many other qualitatively similarities with earthquake dynamics. As far as we know, it is the first time that such a quantitative agreement concerning simultaneously three main laws of seismicity is reported in a shear experiment (see more about the need of quantitative analogies at the Supp. Mat.). Its circular geometry allows the system to run continuously, capturing the considerable statistics required to analyze the dynamics of very large events, which rarely take place.

Experimental System. We study a 2D cylindrical pile confined in between two concentric fixed acrylic cylinders, and bounded by two rough circular rings (Fig. 1a, b, Movie S1), with a dead load placed over the top ring and compressing the granular pile. The top ring is free to move vertically but not to rotate, while the bottom one is slowly rotated with a period of $18.3\bar{3}\leavevmode\nobreak\ \text{hours}$, quasi-statically shearing the granular pile with a linear velocity of $48.84\leavevmode\nobreak\ \text{mm/hour}$ (approximately 12,600 times faster than the San Andreas fault, with an average slip rate of $33.9\leavevmode\nobreak\ \text{mm/year}$ Sieh and Jahns (1984)). Thanks to a lever and a force sensor, we measure the torque $\Gamma(t)$ applied by the granular pile on the top ring. Six piezoelectric pinducers are inserted regularly in the top ring and simultaneously record acoustic emissions (AE). Both measures are done at a rate of 100,000 samples per second. The system is left to evolve for typical times of 24 hours.

During the shear, $\Gamma(t)$ shows an irregular stick-slip like behavior compatible with earthquakes dynamics, consisting in a continuous loading interrupted by intermittent drops with a large distribution of sizes (Fig. 1c). The detection of the torque drops (Fig. 1d) is performed by applying a threshold to mechanical energy variations $\Delta(\Gamma^{2}(t))$ (see Appendix 1 for details). Acoustic events are linked to local releases of energy taking place at the two-dimensional interfaces between grains Dubourg et al. (2017). The analysis of acoustic recordings (Fig. 1e) is based on a Discrete Wavelet Transform (DWT) Chui (1992) resulting in a time-frequency energy distribution map (Fig. 1f), that is then processed to detect peaks, corresponding to the energy of the events (details in the Supp. Mat.). The high number of events detected with both methods, respectively around 2,000 torque drops and more than 1.8 million acoustic emissions for a 24 hours experiment, allows to compute precise statistical characteristics of the system’s behavior and to compare it with the dynamics of earthquakes.

Reproducing main seismicity laws. Gutenberg-Richter: On Fig. 2a we present the probability distribution of both acoustic $E_{\text{ac}}$ and mechanical $E_{\text{m}}$ energies of detected events, on logarithmic intervals. The alignment of both types of energy on the abscissa axis has been obtained from Fig. S1, where synchronous corresponding mechanical and acoustic events are represented. For one given mechanical energy, the associated acoustic one is obtained as the median energy of corresponding acoustic emissions.

The two probability density functions behave like power laws $P(E)\sim E^{-\beta}$ with exponents $\beta_{\text{m}}=1.71\pm 0.04$ and $\beta_{\text{ac}}=1.71\pm 0.01$ for mechanical and acoustics energies respectively, obtained by a maximum likelihood method Clauset et al. (2009); Alstott et al. (2014). The AE energies spread over six decades while the mechanical energies cover only three decades, showing the better sensitivity of the acoustic detection. This power-law behavior is to be compared with the Gutenberg-Richter law Hanks and Kanamori (1979); Choy and Boatwright (1995) which states that the PDF of radiated energies of globally measured real earthquakes follows a power-law with an exponent $\beta=5/3=1.67$.

Inter-event times: By defining a threshold in energy $E_{0}$ (or in magnitude $M=2/3\leavevmode\nobreak\ \log E-2.9$, which is more used in earthquakes studies Choy and Boatwright (1995)) we can analyze the inter-event time between two consecutive events $\tau_{E\geqslant E_{0}}$ (or $\tau_{M\geqslant M_{0}}$). In nature, the rate of seismicity $R_{M\geqslant M_{0}}$, defined as the number of earthquakes larger that a given magnitude $M_{0}$ per unit time varies depending on the region (eg. $R_{M\geqslant 2}\sim 15,000$ earthquakes/year for California Felzer (2008)). However, the distribution of $\theta=\tau_{M\geqslant M_{0}}/\tau^{*}_{M\geqslant M_{0}}$ is a universal function following $f(\theta)\sim\theta^{-0.3}\exp(-\theta/1.5)$ for all seismic zones Corral (2004), where $\tau^{*}_{M\geqslant M_{0}}$ is a characteristic time defined as the inverse of the rate of seismicity $\tau^{*}_{M\geqslant M_{0}}=1/R_{M\geqslant M_{0}}$. Although the rate of $all$ our AE events corresponds to $R=17.35$ events/s (more than 36,000 times higher than Californian earthquakes), which gives a characteristic time $\tau^{*}=57.64$ ms, the $\theta$ distribution of our AE events follows quantitatively the same universal function (Fig. 2b). In the case where all the events are considered, $\theta=10^{-1}$ corresponds to a $\tau=5.764$ ms, which corresponds to the duration of the events (see Fig. 1e). This sets a limit to the inter-event time distributions and we can notice a deviation of the distribution from the power-law behavior for values of $\theta<10^{-1}$ (Fig. 2b).

The universal function tells us that the system “remember” the events of energy $E\geqslant E_{0}$ during a time corresponding to their characteristic time $1/R_{E\geqslant E_{0}}$, where the distribution is a power law. However, this memory, as in the case of earthquakes, is quite weak an carries no predictive capabilities Touati et al. (2009). For longer inter-event times the exponential tail of the distributions indicates that the events are independent.

For large threshold values ($E_{0}$), the distributions deviate from the universal law both for small and large $\theta$ values (Fig. 2b). The increase of short inter-event times is a direct consequence of the increase of the activity associated to aftershocks and foreshocks in the dynamics, which is also a feature of real earthquake data Corral (2004). The increase of long inter-event times is linked to insufficient statistics. In order to verify that, we have analyzed the distribution of inter-event times $\tau_{\ell}$ for large events only ($E_{0}=10^{8}$) (Fig. S2a). It presents two regimes: a clustering of events for short inter-event times ($\tau_{\ell}<\tau^{*}$), and an exponential decay $D(\tau_{\ell})\sim\exp\left(-\tau_{\ell}/\tau_{c}\right)$ indicating that large events separated by long inter-event times ($\tau_{\ell}>\tau^{*}$) are independent and follow a Poissonian process. When we increase the energy threshold $E_{0}$ defining the large events (notice that $R_{E\geqslant E_{0}}\sim E_{0}^{1-\beta}$), we expect a linear relation between $\tau_{c}$ and $\tau^{*}$. However, we find that the increase of $\tau_{c}$ is slower than linear, and the best fit shows a power law with an exponent 0.86 $\pm$ 0.01 (Fig. S2b). This deviation may be caused by a lack of statistics concerning very large events and eventually may be used as an analytical tool to estimate biases in the obtained results due to insufficient statistical sets of data.

Omori: By defining the large acoustic emissions ($E_{\text{ac}}\geqslant 10^{8}$) as mainshocks, corresponding to about 4,500 events, we are able to reveal the existence of foreshocks and aftershocks following Omori laws as for real earthquakes Omori (1894): $n(t)=A/(c+t)^{p}+B$, where $A/c^{p}$ gives the rate increase associated to the mainshock, $B$ the background rate of the earthquakes, $c$ the time offset (positive and close to zero) Enescu et al. (2007) and $p$ the Omori exponent, around 1. On Fig. 2c we show the average of the AE rate of foreshocks and aftershocks around the mainshocks (where time of foreshocks is $t_{-}=t_{m}-t$ and time of aftershocks $t_{+}=t-t_{m}$, with $t_{m}$ the time of the mainshock). The foreshocks rate follows an Omori-like increase $n(t_{-})=9.05\leavevmode\nobreak\ \text{s}^{-1}+2.78/(0.22\leavevmode\nobreak\ \text{s}+t_{-})^{0.8}$ with a reduction of the activity in the last 0.1 s preceding the mainshock. Just after the mainshock, the aftershocks rate presents first a plateau, associated in real earthquakes to catalogues incompleteness (due to large quakes masking smaller ones), followed by a power-law decrease as $n(t_{+})=8.22\leavevmode\nobreak\ \text{s}^{-1}+7.81/(0.12\leavevmode\nobreak\ \text{s}+t_{+})^{1.0}$. When choosing larger threshold values defining the mainshocks (Fig. S2b) the number of events reduces; but the behavior is qualitatively the same.

Around mainshocks. Beyond analyzing the main laws of seismicity, we can also study the system’s behavior around mainshocks. We divide the time axis into windows of 200 ms duration, from 15 seconds before to 10 seconds after a mainshock, allowing computing a local probability distribution of acoustic energies found at a given time of any of the large emission. We find power law distributions, with a variable exponent represented on Fig. 3a. Far from a mainshock, we find a constant value close to 1.85, slightly bigger than the global $\beta$-value of 1.71. This difference is caused by the selection of time windows without extreme events since they are used as reference mainshocks. For about 4 seconds before the mainshock, the exponent shows a continuous but slight decrease that accelerates in the last second to reach a value of about 1.6 just before the mainshock. The first calculated $\beta$-value after the mainshock reaches 1.4 and then it jumps again to 1.85. These low $\beta$-values indicate an abundance of high-energy events, which are independent of the increase of the rate of the events (foreshocks and aftershocks). This subject is a source of controversy in the Statistical Seismology community Davidsen and Green (2011); Lippiello et al. (2012). Indeed, the analysis of local and relative fast variations of the $\beta$-value in real earthquakes is well documented, often associated to correlation between magnitudes Lippiello et al. (2012), but always affected by both the intrinsic lack of statistics and the incompleteness of the catalogues Davidsen and Green (2011). As our experiments rely on single measurements captured during a relative short period (one or a few days), the reported results are less affected from catalogue incompleteness than real earthquake data, making our statistical results very reliable.

We can focus on the two closest data points around the mainshock of Fig. 3a. The corresponding 2D histogram (Fig. 3b) shows that high-energy events cluster just after the mainshock. They have been detected –thanks to the wavelet analysis– inside the acoustic envelopes of the mainshocks, which can last for up to 30 ms. As large events correspond to rearrangements involving a large number of grains, it is expected to detect secondary peaks associated to the mainshocks. A priori the Omori law may not explain the statistics during these dynamic rearrangement, thus magnitude correlations seems the appropriated term to refer to this clustering of high energy events. Indeed, the activity seems to restart progressively from 40 ms having a maximum around 120 ms. These results, coherent with the aftershocks displayed on Fig. 2c, indicate that the flat part of the Omori law between 50 ms and 120 is not provoked by missing small events masked by previous large ones. We can also notice a clear quiescence in the 60 ms interval preceding the mainshock, also coherent with the decrease of the foreshock activity close to the mainshocks (Fig. 2c), a phenomenon that has been often reported prior to very large earthquakes Kanamori (1981). Longterm decrease of $\beta$-values has also been reported preceding very large earthquakes Nanjo et al. (2012), and also in controlled experiments Sammonds et al. (1992) and simulations Kun et al. (2013).

Discussions. A very large distribution of thresholds, required to achieve an earthquake-like dynamics, is directly related to the heterogeneous character of the system. In subcritical fracture experiments like those cited in the introduction Scholz (1968); Måløy et al. (2006); Baró et al. (2013); Stojanova et al. (2014); Barés et al. (2018), the heterogeneity is provided by a structural disorder, and the competition between the advancing crack and the fracture thresholds may result in a Gutenberg-Richter-like distribution of event sizes. However, most shear experiments present a main dynamics composed of quasi-periodic stick-slip events with a narrow distribution of sizes, which may be a consequence of a lack of disorder. In the case of similar experiments shearing a granular layer Rivière et al. (2018); Johnson et al. (2008); Scuderi et al. (2016), a very high number of particles and three-dimensional force chains may be responsible for an “averaging” effect that reduces the heterogeneity of the system, resulting in a regular stick-slip dynamics similar to the one obtained in solid flat interfaces.

In our system, the granular force network Miller et al. (1996); Howell et al. (1999); Majmudar and Behringer (2005) provides an evolving heterogeneity in terms of energy thresholds that is the key ingredient of the dynamics, and it is responsible for a distribution of events that resembles the Gutenberg-Richter law Gutenberg and Richter (1956); Choy and Boatwright (1995). The structure, dynamics and sizes of these heterogeneities in an actual fault remain as open questions. However, the two-dimensional nature of both our system and the inter-grain friction, and the low dimensionality of the force network (which may depend on the pressure between the plates Lherminier et al. (2014)) may serve as hints to eventually find them (see more details in Supp. Mat.). Our system is able to reproduce quantitatively the main statistical laws of seismicity, which indicates that both earthquakes and our experiment are governed by a similar physics, and opens a new pathway to the investigation of earthquake-like dynamics at a laboratory scale.