# Common dependence on stress for the statistics of granular avalanches and earthquakes

Takahiro Hatano, Clément Narteau, Peter Shebalin

► To cite this version:

Takahiro Hatano, Clément Narteau, Peter Shebalin. Common dependence on stress for the statistics of granular avalanches and earthquakes. Scientific Reports, 2015, 5, p. 269-279. 10.1038/srep12280. insu-03579745

HAL Id: insu-03579745

https://insu.hal.science/insu-03579745

Submitted on 18 Feb 2022

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>
OPEN

# Common dependence on stress for the statistics of granular avalanches and earthquakes

Takahiro Hatano<fnref n="1" />, Clément Narteau² & Peter Shebalin³

Received: 19 November 2014

Accepted: 24 June 2015

Published: 21 July 2015

Both earthquake size-distributions and aftershock decay rates obey power laws. Recent studies have demonstrated the sensibility of their parameters to faulting properties such as focal mechanism, rupture speed or fault complexity. The faulting style dependence may be related to the magnitude of the differential stress, but no model so far has been able to reproduce this behaviour. Here we investigate the statistical properties of avalanches in a dissipative, bimodal particulate system under slow shear. We find that the event size-distribution obeys a power law only in the proximity of a critical volume fraction, whereas power-law aftershock decay rates are observed at all volume fractions accessible in the model. Then, we show that both the exponent of the event size-distribution and the time delay before the onset of the power-law aftershock decay rate are decreasing functions of the shear stress. These results are consistent with recent seismological observations of earthquake size-distribution and aftershock statistics.

Individual earthquakes may be regarded as large scale ruptures involving a wide range of structural and compositional heterogeneities in the crust. However, the statistical properties of a population of earthquakes are often described by simple power-laws. Among them, two laws are ubiquitous and occupy a central position in statistical seismology: the Gutenberg-Richter (GR) law<fnref n="1" /> and the Modified Omori law (MOL)²,³. The GR law describes the earthquake magnitude-frequency distribution

$$P(M_w) \propto 10^{-bM_w}, \qquad (1)$$

where $M_w$ is the moment magnitude and $b$ a constant with a value around 1 along active fault zones. The MOL describes the aftershock occurrence rate

$$n(t) \propto (t + c)^{-p}, \qquad (2)$$

where $t$ is the elapsed time from the triggering event (the so-called mainshock), $p$ a positive non-dimensional constant with a typical value of 1 and $c$ a time constant. The parameters in these two laws are believed to bear some information on the physical state of the crust. Indeed, *Schorlemmer et al.* find that the $b$-value of the GR law is decreasing going from normal (extension) over strike-slip (shear) to thrust (compression) earthquakes⁴. *Narteau et al.* also find that the time constant $c$ in the MOL has the same dependence on the faulting mechanism⁵. These two observations indicate that, under a simple assumption, $b$ and $c$ are decreasing functions of shear stress. Although the underlying mechanism needs further investigation, this may reflect a common time-dependent behaviour of fracturing in rocks during the propagation of earthquake ruptures and the nucleation of aftershocks.

Because the shear stress along an active fault is not directly measurable, a solution to address stress dependences in earthquake statistics is to analyse models that implement a restricted set of physical

<footnote n="1">¹Earthquake Research Institute, University of Tokyo, 113-0032 Tokyo, Japan. ²Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Univ Paris Diderot, UMR 7154 CNRS, 1 rue Jussieu, 75238 Paris, Cedex 05, France. ³Institute of Earthquake Prediction Theory and Mathematical Geophysics, Profsoyuznaya str., 84/32, Moscow 117997, Russia. Correspondence and requests for materials should be addressed to T.H. (email: hatano@eri.u-tokyo.ac.jp)</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** A conceptual model of a fault zone using a 3D granular system. We consider that the fault zone is a strongly damaged area that may be investigated through granular mechanics. The granular model represents a thin gouge layer. The constant grain size, the relatively high porosity and the low grain size to system size ratio of the model do not capture the natural structural and compositional complexities of fault zones in nature. Despite these strong differences, we explore similarities between the dynamics of avalanches in the model and earthquakes.</figcaption></figure>

processes. Such models are indeed numerous, ranging from rock fracture experiments$^{6-9}$ to computer simulations on cellular automata$^{10,11}$. Among them, sheared granular media$^{12-19}$ are simple representations of granular fault gouges (Fig. 1), which are commonly used in geophysics to analyse deformation of highly damaged rocks in fault zones$^{20-22}$. Additionally, both the energy and the stress can be easily defined in these models.

Here we use the simplest form of structural complexity - a bimodal distribution of particle size - to elucidate its fundamental control on the emergent scaling laws. We then perform numerical simulations showing that, as for real seismicity, avalanches in sheared granular matter obeys the GR law and the MOL with *b* and *c*-values being decreasing functions of the shear stress. The details of our numerical model are described in the Methods section.

## Results

**Characterisation of avalanches.** As observed in many previous works on amorphous particulate systems$^{14-17,20,21,23-30}$, Fig. 2 shows that the temporal fluctuations of the energy becomes volatile if the shear rate is sufficiently low and the volume fraction is sufficiently high. Under such condition the kinetic energy is negligible in comparison with the elastic energy $E(t)$, and therefore the elastic energy drop should approximate a transition from a local maximum to a local minimum in configurational energy. Then, an avalanche is defined as an abrupt drop of the elastic energy, $E(t_1) - E(t_2)$, where $t_1$ and $t_2$ denote the beginning and the end of an event, respectively. This is illustrated in the inset of Fig. 2b.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Typical time series of the elastic energy for a constant shear rate $\dot{\gamma} = 10^{-7}$ and different volume fraction $\phi$: (top) $\phi = 0.630$, (middle) $\phi = 0.644$, (bottom) $\phi = 0.650$. As shown by the normalization constant of the vertical axes, the total energy and the energy releases explore different ranges of magnitude according to volume fraction. Inset in the middle panel shows how we estimate the energy release associated with a single avalanche.</figcaption></figure>

Following the literature of earthquake studies³¹, we define the magnitude of avalanche as $M = \log_{10}[E(t_1) - E(t_2)] + m_0$, where $m_0 = 11$ is chosen so that the magnitude of the smallest avalanches is approximately zero. Note that, with this definition of magnitude, the GR law reads $P(M) \propto 10^{-2/3bM}$, as the moment magnitude for earthquake³¹ is defined as $M_w \approx 2/3 \log_{10}[E(t_1) - E(t_2)] + \text{const.}$ Hereafter, we use $\beta \equiv 2/3b$ instead of $b$. Other important quantities that characterise an avalanche are the initial stress $\sigma \equiv \hat{\sigma}(t_1)$ and the stress drop $\Delta\sigma \equiv \hat{\sigma}(t_1) - \hat{\sigma}(t_2)$. The definition of the shear stress and other technical remarks on avalanches are described in the Methods section.

**Magnitude-frequency distribution.** First, we discuss the nature of the avalanche magnitude-frequency distribution with respect to the two control parameters, the shear rate $\dot{\gamma}$ and the volume fraction $\phi$. Figure 3a,b show these distributions at several shear rates for $\phi = 0.644$ and $\phi = 0.650$, respectively. Both parameters control the shape of the magnitude-frequency distribution. For example, one can observe a break in scale-invariance for high shear rates at low volume fraction ($\phi = 0.644$ and $\dot{\gamma} = 10^{-5}$ in Fig. 3a). Similarly, characteristic-size distribution (i.e, peaked at a single magnitude) are observed for high volume fraction (Fig. 3b). However, the distribution is independent of the shear rate below a characteristic $\dot{\gamma}$-value, which may be interpreted as the inverse of the structural relaxation time. Not surprisingly, this threshold value is a decreasing function of the volume fraction. In the volume fraction range investigated here, it is approximately $10^{-6}$ for $\phi = 0.644$ (Fig. 3a) and $10^{-8}$ for $\phi = 0.650$ (Fig. 3b). Hereafter, we discuss such rate-independent behaviours by choosing sufficiently low shear rates.

Figure 3c shows rate-independent magnitude-frequency distributions at several volume fractions. They may be regarded as the GR law in the proximity of a critical volume fraction ($\phi \approx 0.644$), whereas they no longer obey the GR law at higher volume fractions. At much lower volume fractions ($\phi < 0.64$), we can hardly obtain a sufficient number of avalanches that ensures statistical significance. It should be noted that the $\beta$-value in the proximity of a critical volume fraction is sensitive to the minute changes in volume fraction. The $\beta$-value is a decreasing function of the volume fraction ranging from 0.47 to 0.78. We obtain the smallest value (0.47) at the largest volume fraction ($\phi = 0.645$) and the largest value (0.78).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Avalanche magnitude-frequency distributions: (a) for a volume fraction $\phi = 0.644$ and three shear rate values, $\dot{\gamma} \in [10^{-5}, 10^{-6}, 10^{-7}]$. (b) for a volume fraction $\phi = 0.650$ and three shear rate values, $\dot{\gamma} \in [10^{-7}, 10^{-8}, 10^{-9}]$. (c) for several volume fraction values at low shear rate, $\dot{\gamma} = 10^{-7}$ for $\phi \le 0.645$ and $\dot{\gamma} = 10^{-8}$ for $\phi > 0.645$.</figcaption></figure>

at the smallest volume fraction ($\phi = 0.642$). The $\beta$-value at intermediate volume fractions ($\phi = \{0.643, 0.644\}$) is approximately 0.64. This range of values is comparable to that of earthquakes⁴ for which the $\beta$-value varies from 0.50 along normal faults (i.e., extensional regime, low stress) to 0.73 along reverse faults (i.e., compressional regime, high stress).

We remark that distribution functions of other quantities also obey power laws. Among them, we find that the distribution functions of the avalanche duration, $t_2 - t_1$, exhibit power-law tails, the exponent of which is approximately 3.0 irrespective of the volume fraction. We also find that the exponent for the stress drop distribution is twice larger than that for energy drop.

As similar power law behaviours have been observed in many particulate systems, it is interesting to compare the present result with other studies on sheared granular matter. Actually, a range of $\beta$-value have been obtained: 0.36 to 0.95¹³, 0.82 to 0.89¹⁴, and 1.0¹⁵ have been reported. On the other hand, the GR law is not observed in an experiment that is conducted at lower volume fractions¹⁶. All these behaviours are consistent with the present result, where the $\beta$-value depends on both the volume fraction and the stress level. Thus, the $\beta$-value should not be an analogue for critical exponents. Some studies also report similar $\beta$-value variations¹⁷⁻¹⁹.

As an evidence of such non-universality, one can further illustrate a wide range of $\beta$-values obtained in various amorphous systems²³⁻³⁰. In a fracture experiment³² and a frictional system³³, the $\beta$-values are also found to be sensitive to control parameters. There may exist many unknown ingredients that control the $\beta$-value, and further investigation is required for the unified understanding of all these exponents in general amorphous systems.

To discuss the effect of shear stress on the $\beta$-value, we introduce the shear stress $\sigma$ at the beginning of an event as an additional argument to the magnitude-frequency distribution. Then, $P(M, \sigma)$ is the conditional probability of observing a magnitude $M$ avalanche under the (global) shear stress value $\sigma$. For convenience, $\sigma$ is integrated in a certain interval $S_i=[10^{-i/2}, 10^{-i/2+0.5}]$ ($i=\{1, 2, \dots, 9\}$). Thus, we obtain the following distribution function: $P_i(M) = \int_{S_i} d\sigma P(M, \sigma)$. Figure 4a shows the behaviours of $P_i(M)$ at $\phi = 0.643$. We can see that the probability of observing a larger avalanche increases as the shear stress increases. More importantly, the distribution function at each stress level develops a power-law tail with a $\beta$-value, which is a decreasing function of the shear stress (Fig. 4b). We confirm that this stress dependence is independent of the volume fraction in the range $0.642 \le \phi \le 0.645$. In addition, the shear stress dependence of the $\beta$-value is qualitatively the same as that in rock fracture experiments⁷.

**Aftershock statistics and the Modified Omori Law.** Second, we can also study time series of avalanches from our simulations to take them as analogues for mainshock-aftershocks sequences. Unlike the magnitude-frequency distribution, and despite the systematic occurrence of aftershocks in seismogenic areas, there are only few examples of such mainshock-aftershocks sequences in amorphous systems⁹,¹⁴. Triggered events may be difficult to distinguish from other events or simply too rare in individual sequences to exhibit a specific decay rate. Here, we use stacks of aftershocks to capture signals over
<figure><img src="image_5.png" /><figcaption>**Figure 4. Dependency of the avalanche magnitude-frequency distribution on the global shear stress. (a)** The avalanche magnitude-frequency distributions for different ranges of shear stress value, a volume fraction $\phi = 0.643$ and a shear rate $\dot{\gamma} = 10^{-7}$. Events $i$ are classified according to the value of the global shear stress at the initiation of the avalanches (See text for the exact ranges of shear stress value). (b) The slope of the magnitude-frequency distribution with respect to the global shear stress. The $\beta$-value is a decreasing function of the global shear stress.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5. The declustering method to select mainshock and aftershocks.** We consider two non-overlapping magnitude ranges for mainshocks [$M_{min}^M, M_{max}^M$] and aftershocks [$M_{min}^A, M_{max}^A$]. A magnitude $M \in [M_{min}^M, M_{max}^M]$ event occuring at time $t_M$ is selected as a mainshock if there is no larger event in the time interval [$t_M - \Delta T; t_M + \Delta T$]. All magnitude $M \in [M_{min}^A, M_{max}^A]$ events in the time window [$t_M, t_M + \Delta T$] are selected as aftershocks (see also the Methods section.)</figcaption></figure>

more than two decades in time. Following geophysical studies$^{5,34,35}$, the definitions of mainshocks and aftershocks as well as the declustering technique are described in Fig. 5 and in the Methods section.

Figure 6a shows the aftershock rates with respect to the time $\tau$ from the mainshock for several ranges of mainshock magnitude, which is denoted by $M^M$. The magnitude of aftershocks is denoted by $M^A$. The studied ranges of $M^M$ and $M^A$-values are chosen so that we have enough events in the numerical output to get sufficient statistics. We find that the MOL (Eq. 2) holds with the exponent $p \approx 1$ irrespective of the magnitude range of mainshocks. Figure 6b,c show the dependence of the aftershock rate on the volume fraction. We test several volume fractions ($\phi = \{0.642, 0.643, 0.644, 0.645\}$) and shear rate $\dot{\gamma} = \{10^{-6}, 10^{-7}, 10^{-8}\}$) to verify the existence of the MOL irrespective of these parameters. This relaxation process, characterised by an initial plateau followed by a power law decay, remains stable across the entire parameter space of the model. This is not the case for the GR law (Fig. 3). In addition, the time constant $c$ is insensitive to the volume fraction, as shown in Fig. 6b,c.

If aftershocks are triggered by external shear, aftershock statistics may be controlled by the shear strain that is applied after a mainshock, $\tau\dot{\gamma}$. In this case the aftershock decay rate $n(\tau)$ must be scaled as $n(\tau\dot{\gamma})$. Taking the MOL into account, this means that $c \propto 1/\dot{\gamma}$. However, we confirm that the $c$-value is independent of the shear rate. This indicates that aftershocks are not driven by the additional shear strain applied after a mainshock; rather, they are caused by the intrinsic relaxation dynamics. Thus, we do not expect aftershocks in the quasi-static shear deformation, in which a system fully relaxes to a stable configuration after an avalanche.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The aftershock decay rate. (a) Occurence rate of magnitude $M^A \in [1, 3]$ aftershocks for different mainshock magnitude ranges and a volume fraction $\phi = 0.644$. The dashed line is proportional to $1/t$ (i.e., the Omori Law). (b) Occurence rate of $M^A \in [1, 3]$ aftershocks for different volume fractions and magnitude $M^M \in [3, 4]$ mainshocks. (c) Occurence rate of $M^A \in [2, 4]$ aftershocks for different volume fractions and magnitude $M^M \in [4, 5]$. In all cases, the shear rate $\dot{\gamma} = 10^{-7}$ and individual aftershocks are stacked according to their main shock times to compensate for the small number of events in each sequence. There are at least 10 mainshocks and 1000 aftershocks in each sequence. Note the time delay before the onset of a power-law aftershock decay rate.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Dependency of the aftershock decay rate on the global shear stress. (a) Occurence rate of $M^A \in [1, 3]$ aftershocks for different ranges of shear stress value and a volume fraction $\phi = 0.644$. Aftershocks $i \in \{0, 1, 2, 3\}$ are classified according to the global shear stress value $\sigma_i \in [\exp(i-10), \exp(i-9)]$ at the initiation of the avalanches. The time delay before the onset of the power-law decay rate is systematically decreasing with the level of stress (i.e., an increasing $i$-value). (b) Negative dependence of the $c$-value on the global shear stress. Circles: $\phi = 0.644$, $M^M \in [3, 4]$, and $M^A \in [1, 3]$; Squares: $\phi = 0.645$, $M^M \in [3, 4]$, and $M^A \in [1, 3]$. Diamonds: $\phi = 0.644$, $M^M \in [2, 3]$, and $M^A \in [1, 2]$.</figcaption></figure>

Next, we analyse how the time constant $c$ depends on the magnitude of the shear stress. We define for aftershocks the stress range [$\sigma_{min}, \sigma_{max}$] and select only aftershocks that belong to this stress range and the magnitude range [$M^A_{min}, M^A_{max}$]. The aftershock rates are shown in Figure 7a, in which the MOL (with $p=1$) holds clearly and $c$ is a decreasing function of shear stress. We estimate the $c$-values by fitting the data with $A/(\tau+c)$ using the maximum-likelihood method. As shown in Fig. 7b, the $c$-value has a
negative dependence on the shear stress. We confirm this stress dependence for two volume fractions ($\phi = \{0.644, 0.645\}$) and for two magnitude ranges ($M^M \in [3, 4]$, $M^A \in [1, 3]$ and $M^M \in [4, 5]$, $M^A \in [2, 4]$). This negative shear-stress dependence of the c-value for aftershocks is consistent with the trend inferred from seismological observations$^{5,34,35}$.

## Discussion and conclusion

Many systems exhibit crackling noise behaviours that share common properties with earthquake statistics$^{36}$. However, none of them as yet been able to systematically capture the dependence of both the GR law and the MOL on the level of stress$^5$. Here, we have not only shown that these two fundamental laws of statistical seismology are relevant to avalanches in sheared granular matter, but also they have common dependence on the level of stress.

Concerning the GR law, which is commonly observed in a wide range of materials under different conditions from laboratory experiments to the field, we find scale invariance in the model only for a narrow range of volume fraction ($\phi \approx 0.644$). This does not contradict the previous studies, most of which involve different boundary conditions. For example, under constant pressure condition$^{13,15}$, a system can be self-organised to a critical volume fraction at sufficiently low shear rates. Additionally, in some experiments conducted at constant volume fractions, the GR law hardly approximates avalanche-size distributions if the volume fraction is lower than a critical value$^{16,17}$. In contrast, the GR law holds for earthquakes irrespective of the details: tectonic setting, depth, rock type, and loading rate. We do not have a clear answer about this difference. Nevertheless, the fact that our model can produce self-similar scalings over more than five orders of magnitude at the critical volume fraction is sufficient to differentiate between the GR law and other types of distribution. Hence, we can explore the sensibility of this power-law regime to different conditions, especially to shear stress. In addition, when the GR law is not respected, a characteristic avalanche size emerges, a behaviour that may be compared to the characteristic earthquake distribution$^{37,38}$.

Concerning aftershocks, we have shown that similar relaxation process may occur in sheared granular media and in the seismogenic crust. In addition, the negative stress dependences of the time delay before the onset of the power-law decay rate is the same as the one observed in real seismicity. This stress dependence is also coincident with the behaviour of some simple models based on a frictional constitutive law$^{39}$, damage mechanics$^{34}$, or static fatigue$^{40,41}$, despite the absence of such physical processes in our granular model. A quantitative result in the present study, that the time constant $c$ in MOL is an exponentially decreasing function of the applied shear stress, should be tested in more realistic and complex systems that are relevant to seismogenic zones as well as various industrial situations.

## Methods

**Simulation model.** Our three-dimensional granular system is made of frictionless spheres with diameters of $d$ and $0.7d$ (the ratio of populations is 1:1). For the sake of simplicity, we assume that the mass $m$ of these particles are the same. We limit ourselves to a rather small-size system ($N=1500$) for computational efficiency. Using the radius and the position of particle $i$, which are denoted by $R_i$ and $r_i$, respectively, the force between particles $i$ and $j$ is written as $\mathbf{F}_{ij} = [kh_{ij} - \zeta\mathbf{n}_{ij} \cdot \hat{\mathbf{r}}_{ij}]\mathbf{n}_{ij}$. Here $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$, $\mathbf{n}_{ij} = \mathbf{r}_{ij}/|\mathbf{r}_{ij}|$, and $h_{ij} = (R_i + R_j) - |\mathbf{r}_{ij}|$ is the overlap length. If $R_i + R_j < |\mathbf{r}_{ij}|$, particles $i$ and $j$ are not in contact so that the force vanishes.

Throughout this study, we adopt the units in which $d=1$, $m=1$, and $k=1$. We choose $\zeta=2.0$, which corresponds to the vanishing coefficient of restitution. A constant shear rate $\dot{\gamma}$ is applied to the system through the Lees-Edwards boundary conditions$^{42}$. Note that under these boundary conditions the system volume is constant. Thus, the important parameters are the shear rate $\dot{\gamma}$ and the packing fraction $\phi$.

To set-up the initial condition, particles are randomly distributed in a simulation box with zero shear stress. When the kinetic energy has relaxed to zero, a constant shear rate is applied. To avoid transient behaviours, we concentrate only on the data for which the total strain is greater than 100%. Simultaneously, we verify that a steady state of uniform shear rate and uniform volume fraction is achieved. Here, we investigate such uniform steady states only.

**Definition of avalanche.** The beginning of an avalanche is taken at time $t=t_1$ at which the total elastic energy $E(t)$ starts decreasing; i.e., $\dot{E}(t_1) = 0$ and $\ddot{E}(t_1) < 0$. Symmetrically, the end of an event is the time $t=t_2$ at which $\dot{E}(t_2) = 0$ and $\ddot{E}(t_2) > 0$. We also calculate the so-called global shear stress $\hat{\sigma}(t)$ using the virial$^{43}$. In analysing the time series, we did not set the noise threshold. Namely, any small increase in energy leads to the termination of an avalanche.

For simplicity, the spatial information of avalanches is discarded. In this case, one might overlook simultaneous avalanches occurring at different places. However, this is unlikely because the present system is small (i.e. the characteristic length is approximately of $9d$).

**Definitions of mainshocks and aftershocks.** Because aftershocks result from changes of stress induced by a mainshock, we disregard the smaller avalanches for which $M<0$ and consider ranges of magnitude for mainshocks [$M_{\min}^M$, $M_{\max}^M$] and aftershocks [$M_{\min}^A$, $M_{\max}^A$]. In practice, a magnitude
$M \in [M_{\min}^{M}, M_{\max}^{M}]$ event occurring at time $t_M$ is not a mainshock if there is at least one avalanche of the same or higher magnitude range in the time interval $[t_M - \Delta T, t_M + \Delta T]$. Thus, we only consider isolated mainshocks and, taking sufficiently large $\Delta T$-values, we also avoid overlapping aftershock sequences that belong to different mainshocks. Then, all $M < M_{\min}^{M}$ avalanches that follow a magnitude $M^{M}$ mainshock within the time interval $[t_M, t_M + \Delta T]$ are regarded as aftershocks. Finally, each after-shock is characterised by its magnitude $M^A$ and the elapsed time $\tau$ since the mainshock.

In order to reduce artefacts related to event detectability, we disregard larger magnitude range for mainshocks and smaller magnitude range for aftershocks. In addition, the magnitude range of main-shocks are chosen within the intermediate magnitude range in which the GR law holds. Using this strict methodology, we considerably reduce the number of events in our artificial catalogues. Therefore, in all ranges of magnitudes, aftershocks are stacked with respect to the time of their mainshocks to finally end-up with a single aftershock sequence.

## References

1. Gutenberg, B. & Richter, C. Frequency of earthquakes in California. *Bull. Seism. Soc. Am.* **34**, 185–188 (1944).

2. Omori, F. On after-shocks of earthquakes. *J. Coll. Sci. Imp. Univ. Tokyo* **7**, 111–200 (1894).

3. Utsu, T. Aftershocks and earthquake statistics (1): some parameters which characterize an aftershock sequence and their interrelations. *J. Facul. Sci. Hokkaido Univ. Ser. VII* **3**, 129–195 (1970).

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

5. Narteau, C., Byrdina, S., Shebalin, P. & Schorlemmer, D. Common dependence on stress for the two fundamental laws of statistical seismology. *Nature* **462**, 642–645 (2009).

6. Mogi, K. Magnitude-frequency relation for elastic shocks accompanying fractures of various materials and some related problems in earthquakes (2nd Paper). *Bull. Earthquake Res. Inst.* **40**, 831–853 (1962).

7. Scholz, C. H. The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes. *Bull. Seism. Soc. Am.* **58**, 399–415 (1968).

8. Hirata, T. Omori's Power Law aftershock sequences of microfracturing in rock fracture experiment. *J. Geophys. Res.* **92**, 6215–6221 (1987).

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

10. Turcotte, D. L. Self-organized criticality. *Rep. Prog. Phys.* **62**, 1377–1429 (1999).

11. Narteau, C. Classification of seismic patterns in a hierarchical model of rupture a new phase diagram for seismicity. *Geophys. Journ. Int.* **168**, 710–722 (2007).

12. Dalton, F. & Corcoran, D. Self-organized criticality in a sheared granular stick-slip system. *Phys. Rev. E* **63**, 061312 (2001).

13. Dalton, F. & Corcoran, D. Basin of attraction of a bounded self-organized critical state. *Phys. Rev. E* **65**, 031310 (2002).

14. Peña, A., McNamara, S., Lind, P. G. & Herrmann, H. J. Avalanches in anisotropic sheared granular media. *Granul. Matter* **11**, 243–252 (2009).

15. Bretz, M., Zaretzki, R., Field, S. B., Mitarai, N. & Nori, F. Broad distribution of stick-slip events in Slowly Sheared Granular Media: Table-top production of a Gutenberg-Richter-like distribution. *EPL* **74**, 1116–1122 (2006).

16. Twardos, M. & Dennin, M. Slow steady-shear of plastic bead rafts. *Granul. Matter* **7**, 91–96 (2005).

17. Behringer, R. P., Bi, D., Chakraborty, B. Henkes, S. & Hartley, R. R. Why Do Granular Materials Stiffen with Shear Rate? Test of Novel Stress-Based Statistics. *Phys. Rev. Lett.* **101**, 268301 (2008).

18. Kun, F., Varga, I., Lennartz-Sassinek, S. & Main, I. G. Approach to failure in porous granular materials under compression. *Phys. Rev. E* **88**, 062207 (2013).

19. Kun, F., Varga, I., Lennartz-Sassinek, S. & Main, I. G. Rupture cascades in a discrete element model of a porous sedimentary rock. *Phys. Rev. Lett.* **112**, 065501 (2014).

20. Mora, P. & Place, D. Simulation of the frictional stick-slip instability. *Pure Appl. Geophys.* **143**, 61–87 (1994).

21. Aharonov, E. & Sparks, D. Stick-slip motion in simulated granular layers. *J. Geophys. Res.* **109**, B09306 (2004).

22. Johnson, P. A. & Jia, X. Nonlinear dynamics, granular media and dynamic earthquake triggering. *Nature* **437**, 871–874 (2005).

23. Durian, D. J. Bubble-scale model of foam mechanics: Melting, nonlinear behavior, and avalanches. *Phys. Rev. E* **55**, 1739–1751 (1997).

24. Tewari, S. et al. Statistics of shear-induced rearrangements in a two-dimensional model foam. *Phys. Rev. E* **60**, 4385–4396 (1999).

25. Kabla, A., Scheibert, J. & Debregeas, G. Quasi-static rheology of foams. Part 2. Continuous shear flow. *J. Fluid Mech.* **587**, 45–72 (2007).

26. Maloney, C. & Lematre, A. Subextensive scaling in the athermal, quasistatic limit of amorphous matter in plastic shear flow. *Phys. Rev. Lett.* **93**, 16001 (2004).

27. Lerner, E. & Procaccia, I. Density scaling of avalanche statistics in amorphous solids. arXiv:1004.3193v1 (2010).

28. Heussinger, C., Chaudhuri, P. & Barrat, J.-L. Fluctuations and correlations during the shear flow of elastic particles near the jamming transition. *Soft Matter* **6**, 3050–3058 (2010).

29. Bailey, N. P., Schiotz, J., Lematre, A. & Jacobsen, K. W. Avalanche size scaling in sheared three-dimensional amorphous solid. *Phys. Rev. Lett.* **98**, 095501 (2007).

30. Salerno, K. M., Maloney, C. E. & Robbins, M. O. Avalanches in strained amorphous solids: Does inertia destroy critical behavior? *Phys. Rev. Lett.* **109**, 105703 (2012).

31. Hanks, T. C. & Kanamori, H. A moment magnitude scale. *J. Geophys. Res.* **84**, 2348–2350 (1979).

32. Ferenc, K., Lenkey, G., Takács, N. & Beke, D. Structure of magnetic noise in dynamic fracture. *Phys. Rev. Lett.* **93**, 227204 (2004).

33. Ciliberto, S. & Laroche, C. Experimental evidence of self organized criticality in the stick-slip dynamics of two rough elastic surfaces. *J. Phys. I France* **4**, 223–235 (1994).

34. Narteau, C., Shebalin, P. & Holschneider, M. Temporal limits of the power law aftershock decay rate. *J. Geophys. Res.* **107**, B122359 (2002).

35. Narteau, C., Shebalin, P. & Holschneider, M. Onset of power law aftershock decay rates in southern California. *Geophys. Res. Lett.* **32**, L22312 (2005).

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

37. Schwartz, D. & Coppersmith, K. Fault behavior and characteristic earthquake: examples from the Wasatch and San Andreas fault zones. *J. Geophys. Res.* **109**, 5681–5698 (1984).
38. Naylor, M., Greenhough, J., McCloskey, J., Bell, a. F. & Main, I. G. Statistical evaluation of characteristic earthquakes in the frequency-magnitude distributions of Sumatra and other subduction zone regions. *Geophys. Res. Lett.* **36**, L20303 (2009).

39. Dieterich, J. A constitutive law for rate of earthquake production and its application to earthquake clustering. *J. Geophys. Res.* **99**, 2601–2618 (1994).

40. Shcherbakov, R. & Turcotte, D. L. A damage mechanics model for aftershocks. *Pure Appl. Geophys.* **161**, 2379–2391 (2004).

41. Ben-Zion, Y. & Lyakhovsky, V. Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology. *Geophys. J. Int.* **165**, 197–210 (2006).

42. Allen, M. P. & Tildesley, D. J. Computer simulation of liquids (Oxford University Press, New York, 1987).

43. Zubarev, D. N. Nonequilibrium statistical thermodynamics (Consultants Bureau, New York, 1974).

**Acknowledgements**

We acknowledge support from JSPS Grand-in-Aid for Scientific Research, *New Perspective of the Great Subduction Zone Earthquake from the Super Deep Drilling (TH)* and the UnivEarthS LabEx program of Sorbonne Paris Cité, ANR-10-LABX-0023 and ANR-11-IDEX-0005-02 (CN).

**Author Contributions**

T.H. and C.N. planned the research. C.N. and P.S. designed the statistical data analysis for granular avalanches. T.H. did the numerical simulation. P.S. made the code for maximum likelihood method. T.H. and C.N. wrote the manuscript. All authors discussed the results.

**Additional Information**

**Competing financial interests:** The authors declare no competing financial interests.

**How to cite this article:** Hatano, T. et al. Common dependence on stress for the statistics of granular avalanches and earthquakes. *Sci. Rep.* **5**, 12280; doi: 10.1038/srep12280 (2015).

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

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/