Fault heterogeneity and the connection between aftershocks and afterslip

TL;DR

This paper presents an analytical and numerical study linking aftershock rates to afterslip dynamics, showing that fault heterogeneity influences seismic activity patterns and reproduces key empirical laws.

Contribution

It introduces a model incorporating fault heterogeneity and friction variability, demonstrating how aftershock rates relate to afterslip and reproducing observed seismic laws.

Findings

Aftershock rate is proportional to afterslip rate with a correction factor.

Omori law arises from logarithmic afterslip evolution.

Synthetic catalogs reproduce Gutenberg-Richter law with realistic b-value.

Abstract

Whether aftershocks originate directly from the mainshock and surrounding stress environment or from afterslip dynamics is crucial to the understanding of the nature of aftershocks. We build on a classical description of the fault and creeping regions as two blocks connected elastically, subject to different friction laws. We show analytically that, upon introduction of variability in the fault plane's static friction threshold, a non trivial stick-slip dynamics ensues. In particular we support the hypothesis \citep{PA04} that the aftershock occurrence rate is proportional to the afterslip rate, up to a corrective factor that is also computed. Thus, the Omori law originates from the afterslip's logarithmic evolution in the velocity-strengthening region. We confirm these analytical results with numerical simulations, generating synthetic catalogs with statistical features in good agreement with instrumental catalogs. In particular we recover the Gutenberg-Richter law with a realistic $b$-value ($b\simeq 1$) when Coulomb stress thresholds obey a power law distribution.

Figures

Equations

$$\lambda (t)=\frac{K}{t+c},$$

$${\tau }_u(t)={\sigma }_N\left({\mu }_c+A\log \left(\frac{\dot{u}(t)}{V_c}\right)\right)$$

$${\tau }_u(t)=k\left(h(t)-u(t)\right)+k_0\left(V_{0}t-u(t)\right),$$

$$\dot{u}(t)=V_c\exp \left(\frac{h(t)}{\rho }+\frac{t}{t_R}-\frac{{\mu }_c}{A}-\frac{u(t)}{{\rho }_0}\right)$$

$$u(t)$$

$$F(t)$$

$$h(0)\to h(0)+\Delta h$$

$$u(0)\to u(0)+q\Delta h.$$

$$F(t)$$

$$u(t)$$

$$u(t)=u(0)+nq\Delta h+{\rho }_0\log \left(1+\frac{e^{\frac{n\Delta h}{\rho }}}{t_R}t\right).$$

$$t_{AS}=\frac{t_R}{e^{\frac{n\Delta h}{\rho }}}\left(e^{\frac{\delta F_h}{k{\rho }_0}}-1\right)\text{mbox}with{t}_{AS}\ll t_R.$$

$$u(t)=V_{c}t+\text{mbox}const.$$

$$\lambda (t)=Q\left(k(u(t)-h(t))\right)k\frac{\partial (u(t)-h(t))}{\partial t}$$

$$Q(x)=\frac{g\left(x\right)}{{\int }_x^{\infty }g(f)df}.$$

$$\lambda (t)=Q\left(k(u(t)-h(t))\right)k\dot{u(t)}.$$

$$t_M=\frac{1}{k{V}_c}\delta F_h=t_R\frac{\delta F_h}{k{\rho }_0}$$

Full text

I Introduction

A mainshock is followed by the increase of seismic activity caused by aftershocks as well as by a significant time-dependent postseismic deformation known as afterslip. The common interpretation is that afterslip is activated by the stress increase due to the mainshock’s coseismic instantaneous deformation and mostly occurs in regions with a velocity strengthening rheology. This interpretation dates back to the seminal paper by Marone et al. (1991) who analytically obtained a hyperbolic time decay of the post-seismic deformation rate $\dot{u}(t)\sim 1/t$, for a velocity strengthening region . The hyperbolic temporal decay is similar to the one usually observed for the aftershock occurrence rate, $\lambda(t)$:

[TABLE]

and known as the Omori law.

The proportionality between aftershock occurrence rate, $\lambda(t)$, and stress or strain rate, $\dot{u}(t)$, has been documented by the postseismic deformation measured after several large earthquakes (Perfettini and Avouac, 2004; Perfettini et al., 2005; Hsu et al., 2006; Perfettini and Avouac, 2007; Perfettini and Ampuero, 2008; Savage and Yu, 2007; Savage and Langbein, 2008; Savage, 2010; Canitano et al., 2018; Perfettini et al., 2018). Perfettini and Avouac (2004) explained the observed proportionality under the assumption that aftershocks are induced by afterslip. This result was supported by the analytical solution of a single spring-slider model under velocity strengthening friction which models brittle creeping. The proportionality has been also used to obtain the friction parameters as well as the stressing rate, in a given velocity strengthening region, from the temporal behavior of the recorded aftershocks (Frank et al., 2017).

In this letter we present a minimal model of the fault as a sliding block connected to the afterslip region treated, as in Marone et al. (1991); Perfettini and Avouac (2004), as a second block with velocity strengthening rheology (see Fig. 1). In this two-block model the slip of the fault block induces the afterslip relaxation of the second block which in turns promotes further failures of the fault. The central assumption is that instability thresholds $f^{th}$, on the fault plane, is not uniform but random with a distribution $g(f^{th})$. This makes more realistic the description of the fault as a single-slider block. Previous studies (Kaneko and Lapusta, 2007; Ader et al., 2014), indeed, have shown that the absence of heterogeneities is responsible for important differences between single-slider models and two-dimensional continuum models. The presence of random thresholds allows us to analytically demonstrate that, without any assumption on the initial stress distribution, the proportionality $\lambda(t)\propto\dot{u}(t)$ is a stable feature of earthquake triggering. We present analytical and numerical results of the model evolution.

II The model

We describe the fault as a block $H$ coupled by means of a spring of elastic coupling $k$ to a second block $U$ which represents the creeping region of the crust, with velocity strengthening friction. Marone et al. (1991) identified the region where afterslip occurs in a zone of unconsolidated sediments above the fault whereas Perfettini and Avouac (2004); Perfettini et al. (2005); Perfettini and Avouac (2007) assumed the creeping region zone to be deeper because of the transition towards ductile behavior caused by the increase of temperature with depth. In other studies afterslip has been observed very close (Crescentini et al., 1999) or even on the fault within the seismogenic zone (Miyazaki et al., 2004; Johnson et al., 2006; Freed, 2007). In our approach, the precise position of the velocity strengthening region is not relevant and we assume that the block $U$ is embedded in a more extensive region creeping at a constant velocity $V_{0}$. We model this interaction with a second spring of elastic constant $k_{0}$ whose free end moves with velocity $V_{0}$ (Fig. 1). To simplify the model we assume that both blocks move in the $V_{0}$ direction and, therefore, only scalar quantities can be considered. Indicating with $u(t)$ and $h(t)$ the positions of the block $U$ and the block $H$ at the time $t$, respectively, the shear stress acting on $U$ is $k\left(h(t)-u(t)\right)+k_{0}\left(V_{0}t-u(t)\right)$ and the friction force $\tau_{u}(t)$, under steady state condition, can be written as

[TABLE]

where $\dot{u}(t)$ is the block velocity, $\sigma_{N}$ is the effective normal stress, $\mu_{c}$ is the static friction coefficient when the block $U$ slides at the steady velocity $V_{c}$ and $A>0$ for a velocity strengthening material. In the overdamped limit, the constitutive equation for $U$ reads:

[TABLE]

which admits a stationary solution for $\dot{u}(t)=V_{c}=\frac{k_{0}}{k_{0}+k}V_{0}$. In general one has

[TABLE]

where $\rho_{0}=\frac{A\sigma_{N}}{k+k_{0}}$ and $\rho=\frac{A\sigma_{N}}{k}$ are two characteristic length whereas $t_{R}=\frac{A\sigma_{N}}{k_{0}V_{0}}=\frac{\rho_{0}}{V_{c}}$ represents the long time scale associated with the extended creeping region at velocity $V_{0}$.

We next consider the evolution of the block $U$ starting from the stationary solution $\dot{u}(0)=V_{c}$ at the time $t=0$. The solution for $u(t)$ was written in Perfettini and Avouac (2004) in term of the evolution of $h(t)$:

[TABLE]

Concerning the dynamics of the block $H$ we assume the Coulomb failure criterion (CFC): The instability is only controlled by the Coulomb stress $S_{H}(t)$ acting at time $t$ on the block $H$, with $S_{H}(t)=k(u(t)-h(t))-\mu_{H}\sigma_{H}(t)$. Here $k(u(t)-h(t))$ is the shear stress, $\mu_{H}$ is a friction coefficient and $\sigma_{H}(t)$ the effective normal stress represented by the normal stress reduced by the pore pressure. According to the CFC, under the assumption that $\sigma_{H}(t)$ is constant, the block $H$ is unstable as soon as the shear stress overcomes a reference frictional stress $f^{th}$. Treating each slip as instantaneous implies that $H$ either is at rest or it slips in an infinitely short time. This approximation corresponds to a vanishing nucleation size coherently with the steady-state approximation in the constitutive equation (Eq.(3)) (Rice and Ben-Zion, 1996). Within our hypothesis, the dynamics of the two blocks can be written in terms of the Eqs.(5) and split in three regimes:

• •

Slip regime. If at time $t=0$, $k(u(0)-h(0))>f^{th}$, the position of the block $H$ becomes unstable and jumps of $\Delta h$ inducing the coslip of the block $U$ which depends on the precise dynamics of the block $H$. For example in Fig. 1 we present an explicit solution when $h(t)$ moves steadily at velocity $V_{s}\gg V_{0}$ for a short time $t_{s}$. In all cases the block $U$ slips coseismically with the block $H$ of an amount $q\Delta h$ with $q\in[0,1]$, so that the slip is given by:

[TABLE]

The precise value of $q$ is a complicated function of model parameters and of the specific slip profile. In our minimal model we assume that $q$ is a fixed value equal for all slips.

At the end of the slip, the block $H$ experiences a stress drop $k(1-q)\Delta h$, and if $k(u(t_{s})-h(t_{s}))<f^{th}$ the fault is stuck again $h(t)=h(t_{s})$. In the hypothesis of a uniform frictional stress $f^{th}$, the fault exhibits a trivial dynamics characterized by jumps of equal size ($\Delta h$) separated by a constant time interval. In this case the pre-stress value $f^{th}-k(u(0)-h(0))$ is peaked at a characteristic value and stress and seismic rate are not proportional. Nevertheless, different sources of heterogeneity should affect the local value of the friction coefficient: The presence of asperities which induces fluctuations in the normal and pore pressure values, the variation of the basal friction coefficient or, again, the fluctuations of $A$. These parameters are temperature and water dependent and are likely to vary along the fault. Thus it is reasonable that the frictional stress experienced after a slip $\Delta h$ is typically different from the previous one (see the variable intensity of pinning points – red dots – in Fig. 1). Concretely, after each slip it is reasonable that the block $H$ comes up against a new value of $f^{th}=f_{1}^{th}$ extracted from a distribution $g(f^{th})$ with a finite probability that $k(u(t_{s})-h(t_{s}))>f_{1}^{th}$. Hence, the block is still unstable and accordingly can perform $n$ subsequent slips until the distance to failure $\delta F_{h}=f_{n}^{th}-k\left(u(0)-h(0)\right)-n(1-q)\Delta h>0$. In this case, the slip instability corresponds to a single slip of size $n\Delta h$ (a single earthquake) and the dynamics is characterized by earthquakes of different sizes. Large earthquakes can be, therefore, viewed as a succession of smaller seismic ruptures. After the slip the block $H$ is stuck in the novel position $h(t)=h(0)+n\Delta h$ and the evolution of $u(t)$ is given by the explicit solution of Eqs.(5):

[TABLE]

Two distinct regimes can be identified.

• •

Afterslip regime. At times shorter than $t_{R}$ one can replace $t_{R}\left(e^{\frac{t}{t_{R}}}-1\right)\simeq t$ so that:

[TABLE]

In this regime the motion of $U$ increases the stress on $H$ inducing further slips of the block $H$, i.e. the aftershocks. This occurs if $k\rho_{0}\log\left(1+\frac{e^{\frac{n\Delta h}{\rho}}}{t_{R}}t_{AS}\right)=\delta F_{h}$, namely

[TABLE]

• •

Interseismic regime. At $t>t_{R}$ the motion of the block $U$ is dominated by the creeping velocity $V_{0}$. In this regime one can assume $\left(e^{\frac{t}{t_{R}}}-1\right)\simeq e^{\frac{t}{t_{R}}}$ and neglect the one in the logarithm of Eq.(10). We obtain that the block $U$ slides at the steady velocity $V_{c}=\frac{k_{0}}{k+k_{0}}V_{0}$:

[TABLE]

III Statistics of slip events: analytical results

In the following we explore the statistical features of the evolution of a single fault represented by the block $H$ which, because of the coupling with the block $U$, can experience multiple slip events. The statistics of the single fault is expected to be representative of the statistics of a population of faults with random initial stress conditions.

The block $H$ exhibits a stick-slip dynamics with non trivial temporal correlations for a sufficiently broad distribution $g(f)$ which leads to a broad distributed distance to the next instability. In particular, we indicate with $t_{0}=0$ the time of the last slip and consider the probability $P_{0}(t)$ of no slip up to time $t$. Since $u(t)$ is monotonically increasing with time (Eq. (5)), $P_{0}(t)$ corresponds to the probability to extract a friction threshold $f_{1}^{th}$ larger than $k(u(t)-h(t))$, $P_{0}(t)=\int_{k(u(t)-h(t))}^{\infty}g(f)df$. In the hypothesis that $g(f)$ does not present sharp discontinuities, the evolution of the block $H$ can be described as a time-dependent non-homogeneous Poisson process with $P_{0}(t)=\exp{\left(-\int_{0}^{t}\lambda(t^{\prime})dt^{\prime}\right)}$ where $\lambda(t)$ is the seismicity rate, i.e. the number of earthquake for unit time triggered by the mainshock. As a consequence, $\lambda(t)=-\frac{\partial\log\left(P_{0}(t)\right)}{\partial t}$ leading to

[TABLE]

with

[TABLE]

Taking into account that $h(t)=const$, outside the very short slip regime, we find

[TABLE]

Because of the presence of the time-dependent coefficient $Q\left(k(u(t)-h(t))\right)$, Eq.(15) shows that in our model the seismic rate $\lambda(t)$ is not exactly proportional to the stress rate. However, we argue that in the afterslip regime $Q(k(u-h))$ can only depend logarithmically on time and the proportionality between aftershock and stress rate is, at first order, satisfied. More precisely, we first focus on three types of distributions for the friction thresholds: a Gaussian distribution $g(f)=\sqrt{\frac{\alpha}{\pi}}\exp{\left(-\alpha\left(f-f_{0}\right)^{2}\right)}$ restricted to positive $f$, a power law distribution $g(f)=\frac{(\beta-1)f_{0}^{\beta-1}}{(f+f_{0})^{\beta}}$ and an exponential distribution $g(f)=\frac{1}{\gamma}\exp\left(-\gamma f\right)$. In the case of the exponential distribution $Q(k(u-h))$ is exactly a constant whereas it displays logarithmic time dependence for the power law ($Q(k(u-h))\sim 1/(k(u-h))\sim 1/\log(t)$) and for the Gaussian distribution ($Q(k(u-h))\sim k(u-h)\sim\log(t)$). Therefore, using $\frac{\partial k(u-h)}{\partial t}=k\dot{u}(t)\sim 1/t$ in the afterslip regime (Eq. (10)), we always obtain the Omori hyperbolic decay $\lambda(t)\sim\dot{u}\sim 1/t$, with possible logarithmic corrections coming from $Q(t)$. An important exception to this behaviour is represented by the distributions with an upper cut-off $f_{\max}$. In this case, one can easily derive111for example by considering the uniform distribution $[0,f_{\max}]$ that $Q(k(u-h))\propto 1/(f_{\max}-k(u(t)-h(t)))$ which can diverge at finite time. Our results then shows that the key ingredients for the Omori law are represented by heterogeneities in the frictional stress values combined with the logarithmic stress relaxation induced by the velocity strengthening rheology. This behaviour is very general provided that the maximal values of the frictional stress are large compared to the increase of the stress that can be experienced during the afterslip phase.

The presence of frictional heterogeneities also produces a non trivial slip size distribution $p(n)$, which corresponds to the probability that the fault $H$ performs a slip of size $n\Delta h$. Assuming a small stress drop $k(1-q)\Delta h$ and the independence between subsequent jumps we can us the mapping to a record problem. More precisely we neglect the stress change after each slip assuming $k(u(t)-h(t))-k(1-q)\Delta h\simeq k(u(t_{0})-h(t_{0}))$, and calculate $p(n)$ as the probability to draw $n+1$ independent and identically distributed random variables $f_{0}^{th},f_{1}^{th},...,f_{n}^{th}$, such that $f_{j}^{th}<f_{0}^{th}$ for $j\in[1,n-1]$ and $f_{n}^{th}>f_{0}^{th}$. In this case we find $p(n)\propto 1/\left(n^{2}+n\right)$ for any probability density function $g(f_{0})$ and independently of its domain $(f_{\min},f_{\max})$ (here we assumed $f_{\min}\geq 0$ and $f_{\max}=\infty$). To show it, we remark: $p(n)=\int_{f_{\min}}^{f_{\max}}df_{0}g(f_{0})\left(P_{<}(f_{0})\right)^{n-1}\times P_{>}(f_{0})$, where $P_{<}(f_{0})=\int_{f_{\min}}^{f_{0}}df^{\prime}g(f^{\prime})$ is the probability to draw a threshold smaller than $f_{0}$ and, similarly, $P_{>}(f_{0})=1-P_{<}(f_{0})$ which gives $g(f_{0})=\frac{dP_{<}(f_{0})}{df_{0}}$. As a consequence, $p(n)=\int_{f_{\min}}^{f_{\max}}df_{0}\frac{dP_{<}(f_{0}))}{df}$ $\left(\left(P_{<}(f_{0})\right)^{n-1}-\left(P_{<}(f_{0})\right)^{n}\right)$ $=\frac{1}{n}-\frac{1}{n+1}$. For sufficiently large $n$ we obtain the power law behavior $p(n)\sim n^{-\eta}$ with $\eta=2$. This is in qualitative agreement with the Gutenberg-Richter (GR) law for the magnitude distribution. Indeed, taking into account that the size of a slip $n\Delta h$ is proportional to the seismic moment released during a slip instability, the GR law combined with the logarithmic relation between magnitude and seismic moment corresponds to a power law behavior $p(n)\sim n^{-\eta}$ with $\eta=1+(2/3)b$, where $b\simeq 1$ is the coefficient of the GR law.

IV Numerical Results

In numerical simulations, the block $H$ is at rest at time $t_{0}$. We consider $k(u(t_{0})-h(t_{0}))<f^{th}$ and $u(t)$ evolving according to Eq.(12) which corresponds to the interseismic regime. The evolution of the block $U$ increases the shear stress $\tau_{h}(t)$ and the first slip occurs at the time $t_{0}+t_{M}$ with

[TABLE]

obtained by the inversion of Eq.(12). The block $H$, then performs $n$ subsequent slips before reaching the stable condition $\tau_{h}^{(n)}<f_{n}^{th}$ which corresponds to an earthquake of size $n\Delta h$ occurred at the time $t_{0}+t_{M}$. The subsequent slip will occur at a time $t_{0}+t_{M}+\delta t$ where $\delta t$ is the minimum between $t_{AS}$ and $t_{M}$ given by Eq.(11) and Eq.(16), respectively. If $\delta t=t_{AS}$ the new slip is correlated to the first one and they are considered belonging to the same seismic sequence. The process is iterated and therefore the sequence can contain many correlated earthquakes. A new sequence starts as soon as $\delta t=t_{M}$. In numerical simulations we fix the value $\Delta h=10\rho$ in order to have a clear time separation between $t_{AS}$ and $t_{M}$ (Eq.s(11),(16)) and vary the only the $g$-related parameters $\alpha,\beta,\gamma,f_{0}$.

A typical numerical sequence is plotted in Fig. 2a, where each point corresponds to an event with occurrence time $t$ and size $n\Delta h$. The presence of correlated sequences, i.e. the clusters, is manifest. In the majority of cases (more than $70\%$) the first event of the cluster is also the largest event belonging to the same sequence (Fig. 2c). In the remaining sequences the largest event, the mainshock, is anticipated by few smaller events, the foreshocks (Fig. 2b).

In Eq.(11) we define aftershocks as those events following the first event of the sequence. In the following we adopt the more standard definition as those events following the mainshock. Since the number of foreshocks is always much smaller than the aftershock one, statistical features of aftershocks do not depend on the specific definition. Fig. 3 shows the rate of aftershocks $\lambda(t)$ as a function of time since the mainshock which is in good agreement with the hyperbolic time decay predicted by the Omori law, for any distribution $g(f)$ and for different choices of parameters.

In Fig. 4 we plot the slip size distribution $p(n)$ in the synthetic catalog. Results of numerical simulations (Fig. 4) show that for all choices of the distribution $g(f)$ and for different values of parameters, $p(n)$ exhibits a power law decay. We find the exponent $\eta\simeq 2$ for the Gaussian and exponential $g(f)$. Very interestingly, for a power law distributed $g(f)$ and for a wide range of $\beta$-values ($\beta\in[3,8]$), we find $\eta\in[1.6,1.7]$ corresponding to a $b$-value close to $1$ which is in quantitative agreement with the GR law of instrumental catalogs.

V Conclusions

We have presented a minimal model for a seismic fault described as a system of two elastically interconnected blocks, with a block representing the fault plane being subject to a random frictional force as opposed to another block having a velocity strengthening rheology. The evolution of the model presents the typical stick-slip behavior of real fault systems with mainshocks followed by aftershocks distributed in time according to the Omori law. Furthermore, the model reproduces the GR with a realistic $b$-value when the friction distribution $g(f)$ is a power law. This supports the hypothesis of power law distributed Coulomb stress thresholds in agreement with other indications of the self-similar nature of seismic occurrence (Scholz, 2002).

The main difference with previous results is the presence of randomness in the frictional thresholds. This makes the distance to the failure in our model always broadly distributed and allows us to avoid any assumption on the pre-stress distribution. In our study, indeed, the stress conditions before each slip instability originate from the previous stage of earthquake occurrence. This is a novel result with respect to other scenarios where the Omori decay is obtained only starting from a population of fault patches in which pre-stresses are uniformly distributed. Assuming the CFC, for instance, the number of faults which slip in the time interval $\Delta t$ are those with $k(u(0)-h(0))\in[f^{th}-\dot{\tau}(0)\Delta t,f^{th}]$ and are proportional to the stress rate $\dot{\tau}(0)$ for a uniformly distributed $f^{th}-k(u(0)-h(0))$. In particular, under stationary conditions $\dot{\tau}(t)=const$ this hypothesis gives a steady seismicity rate, whereas the Omori decay is obtained when the stress on the fault increases because of afterslip relaxation (Perfettini et al., 2005; Perfettini and Avouac, 2007; Perfettini et al., 2018). In an alternative interpretation (Dieterich, 1994), still assuming that pre-stresses are uniformly distributed, aftershocks are produced by the response of a population of rate weakening patches to a coseismic stress change. In presence of a rate weakening friction, indeed, the relationship between $\lambda(t)$ and $\dot{u}(t)$ is no further linear (Helmstetter and Shaw, 2006) and $\lambda(t)$ decays consistently with the Omori law in response to a stress drop. Savage (2010) has shown that the combination of this relationship, which holds for a velocity weakening fault (Dieterich, 1994), with the hypothesis that the driving stress is generated by afterslip evolution (Perfettini and Avouac, 2004) leads to an improved description of post-seismic relaxation. A rate weakening description of the fault block $H$ can be incorporated in our two-block model instead of the more simple CFC. In our approach we do not include a rate weakening description of the block $H$ and keep the more simple CFC since it allows us to perform analytical calculations and, at the same time, makes numerical simulations much more simple. Nevertheless, because of heterogeneities in the frictional thresholds and since slip instabilities occur on the instantaneous time scale $t_{s}$, details of the friction law acting on the $H$-block are not expected to be relevant. The main effects of the friction laws should reflect in changes of the fault slip $\Delta h$ which, however, is not a relevant parameter of our model. Indeed no significant differences are observed in numerical simulations where $\Delta h$ is not constant but is randomly extracted from a Gaussian distribution.

The comparison with instrumental aftershock sequences would be the subsequent step to support our main conclusions. A further aspect, not investigated in this paper, concerns the understanding of the mechanism leading to the presence of small earthquakes, i.e. the foreshocks, which in our numerical simulations anticipate the occurrence of large earthquakes. This could provide new insights on the outstanding question of the nature of foreshocks (de Arcangelis et al., 2016; Lippiello et al., 2017). As a further step, the two-block model can be natural generalized to the many-block Burridge-Knopoff (BK) model (Burrige and Knopoff, 1967) coupled with a velocity strengthening region. Recent studies (Jagla, 2010; Jagla and Kolton, 2010; Jagla et al., 2014; Lippiello et al., 2015; de Arcangelis et al., 2016; Landes and Lippiello, 2016) have shown that, introducing an intermediate time scale for relaxation in a cellular automata version of the BK model, one recovers statistical features of instrumental catalogs, such as the Omori law and the GR law with a realistic $b$-value. Our results represent a justification for this class of models and provide insights in the mechanisms leading to realistic aftershock features.

Data and Resources No data were used in this paper.

Acknowledgments We thank Eduardo Jagla and Hugo Perfettini for useful discussions. This work was partially supported by the grant from the Simons Foundation, for FL (# 454935 Giulio Biroli, # 327939 Andrea Liu, # 454951 David Reichman)