# Seismic quiescence in a frictional earthquake model

Seismic quiescence in a frictional earthquake model Summary We investigate the origin of seismic quiescence with a generalized version of the Burridge–Knopoff model for earthquakes and show that it can be generated by a multipeaked probability distribution of the thresholds at which contacts break. Such a distribution is not assumed a priori but naturally results from the aging of the contacts. We show that the model can exhibit quiescence as well as enhanced foreshock activity, depending on the value of some parameters. This provides a generic understanding for seismic quiescence, which encompasses earlier specific explanations and could provide a pathway for a classification of faults. Statistical methods, Earthquake dynamics, Theoretical seismology 1 INTRODUCTION While two main seismological laws, the Gutenberg–Richter (GR) law for the magnitude distribution of earthquakes (Gutenberg & Richter 1954) and the Omori law for the time evolution of the frequency of aftershocks (Omori 1894; Utsu et al.1995), are well-established empirically and even may be justified theoretically (Burridge & Knopoff 1967; Olami et al.1992; Socolar et al.1993; Grassberger 1994; Rundle & Klein 1995; Rice & Ben-Zion 1996; Stein 1999; Hainzl et al.2000a,b; Helmstetter et al.2004; Hergarten & Krenn 2011; Serino et al.2011; Jagla 2010; Jagla & Kolton 2010; Jagla 2013; Braun & Peyrard 2013; Mega et al.2003; Kawamura et al.2012; Pelletier 2000), our knowledge about foreshocks is still limited. If a clear pattern could be distinguished in the foreshock activity it might provide an early warning for earthquakes, but we are far from this stage because the observations reveal very different pictures. Many earthquakes are preceded by foreshocks, but their frequency can nevertheless vary widely depending on the type of earthquake (Bouchon et al.2013). Moreover, some earthquakes are preceded by an unexpected calm period, lasting for several hours or more (Jones & Molnar 1979). It is such a period of quiescence, viewed as a characteristic feature of the imminence of the main shock that allowed the only successful prediction of an earthquake (Raleigh et al.1977), which saved a large number of lives in China in 1975. Although the mechanisms that generate foreshocks are still under discussion (Bouchon et al.2013), the existence of some activity before a major earthquake (MEQ) does not seem surprising. In recent experiments with slider pushed on the side, Fineberg et. al. (Rubinstein et al.2004, 2007, 2011) showed that the global sliding, which corresponds to an MEQ, is preceded by propagation of several cracks called precursors that may be used as indicators of the upcoming MEQ. This is also consistent with some ideas that relate earthquakes to self-organized criticality or cascade mechanisms (Olami et al.1992; Socolar et al.1993; Grassberger 1994; Helmstetter et al.2004; Hergarten & Krenn 2011). The fact that a quiescence period could be a characteristic announcement of an MEQ looks more surprising. Various mechanisms have been considered to explain seismic quiescence (Hergarten & Krenn 2011; Main & Meredith 1991). Here, we propose a generic mechanism related to the distribution of the threshold for the breaking of the contacts along a fault and illustrate it by simulations of a simple earthquake model. The Burridge and Knopoff (BK) spring-block model (Burridge & Knopoff 1967) further developed in a number of works (Olami et al.1992; Hainzl et al.2000a,b; Helmstetter et al.2004; Hergarten & Krenn 2011; Serino et al.2011; Jagla 2010; Jagla & Kolton 2010; Jagla 2013; Braun & Peyrard 2013; Braun & Tosatti 2014; Kazemian et al.2015) [see also reviews (Pelletier 2000; Kawamura et al.2012), and references therein] has been used as a generic model which reproduces many features of earthquakes. Simulations showed that some generalized versions of the BK model may demonstrate the GR law and even the Omori law, but the existence of foreshocks in this model is not clearly demonstrated and explained yet (e.g. see Hainzl et al.2000a,b; Helmstetter et al.2004; Hergarten & Krenn 2011; Jagla 2010; Jagla & Kolton 2010; Jagla 2013; Kazemian et al.2015; Pelletier 2000). In BK-type models, the top block (the slider) is coupled with the bottom block (the base assumed to be fixed) by a set of frictional contacts. When the slider moves, a frictional contact i, modeled as an elastic spring, elongates so that the force at the contact point, fi = kxi (we assume that all contacts have the same rigidity k for the sake of simplicity), increases until it reaches a threshold value fsi for which the contact breaks. Then it forms again, with zero stretching xi ∼ 0, and the process can repeat. The evolution of such a system is determined by the distribution of thresholds Pc(fsi). The general features of this distribution also determines the nature of the foreshock activity, as one can realize by considering two simple examples. If Pc(fsi) is one-peaked (Braun & Peyrard 2010), e.g. Gaussian with the centre at fs and a standard deviation Δfs, then the total frictional force on N contacts, $$F = \sum _{i=1}^N f_i$$ linearly increases with the slider displacement X until the contact forces reach values fi ∼ fs − Δfs. Then some contacts start to break, one by one, giving rise to small earthquakes. The growth of F(X) slows down and turns into a decay for larger values of X, and under certain conditions (Braun 2015) an elastic instability may occur. The top block slides for a distance ΔX ≳ Δfs/k, while most of the contacts break. This event corresponds to a large earthquake. In this case, the earthquake activity increases just before the MEQ. Such a scenario is common in natural earthquakes (Kazemian et al.2015; Bowman & King 2001; Jordan & Jones 2010). If the threshold distribution is two-peaked, e.g. when the interface consists of two types of contacts, weak and strong ones with thresholds fs1 and fs2, respectively, as recently considered in Kazemian et al. (2015), the scenario is qualitatively different. To see how it works let us consider a simple case. We assume that, when a contact breaks it forms again as a weak contact with probability p and strong contact with probability 1 − p, and that the breaking of a contact reduces the friction force by φ. At time t0 = 0, we start with $$n_1^{(0)} = p N$$ weak contacts and $$n_2^{(0)} = (1-p)N$$ strong contacts. If the slider moves at velocity v the driving force grows as Nkvt. When it reaches fs1 at time t1, all weak contacts break. The friction force drops by $$\varphi n_1^{(0)}$$ and the broken contacts reform, leading to $$n_1^{(1)} = p n_1^{(0)} = p^2 N$$ weak contacts and a growing number of strong contacts. The same process can repeat again, leading to an infinite series of foreshocks, while the number of weak contacts decreases as $$n_1^{(j)} = p^{j+1} N$$. During the time interval tj − tj − 1, the force has to grow enough to reach fs1 again after the breaking $$n_1^{(j-1)}$$ weak contacts. It requires a time tj − tj − 1 = φpjN/(Nkv), so that the total time for the infinite series of foreshocks stays finite $$\sum _{j=1}^{\infty } (t_{j} - t_{j-1}) = (\varphi / kv) p / (1 - p)$$. After that time only strong contacts persist. As the slider keeps moving, there is a calm period without any foreshock until the strong threshold fs2 is reached. Then all strong contacts break suddenly, leading to an MEQ. This picture of contacts with a two-peaked distribution is oversimplified, and would not show all features of real earthquakes, such as the GR law or the Omori law for aftershocks, but it shows how the properties of the contact thresholds are sufficient to lead to seismic quiescence. Let us now consider a more realistic case and investigate its properties. 2 MODEL We consider the earthquake model previously developed in Braun & Tosatti (2014) and Braun & Scheibert (2014). It is based on the BK-type model, i.e. a set of frictional contacts between a slider and a base, but with the following two important ingredients. First, by contact we mean a macrocontact consisting of a large number of microcontacts on the area $$\lambda _c^2$$, where λc ≈ a2E/k is the elastic correlation length (Caroli & Nozieres 1998; Braun et al.2012) (here a is the average distance between the microcontacts and E is the Young modulus of the slider). For a frictional metal/metal interface λc is typically of the order of μm (Caroli & Nozieres 1998; Braun et al.2012), but for a seismic fault it may be much larger. On the distance λc, the slider may be treated as rigid, while at larger distances we have to account for the deformation of the block. Thus, in our model a contact is a macroscopic object. It has its own ‘degrees of freedom’, which may lead to its aging resulting in a growth of its threshold value as the lifetime of the contact increases. Contact aging is a complex stochastic process which can have various physical or chemical origins. We model it by a simple Langevin equation (Braun & Peyrard 2013; Braun & Tosatti 2014)   $$df_{si} (t)/dt = B(f_{si}) + G \, \xi (t) \,,$$ (1)where B(f) and G are the so-called drift and stochastic forces, respectively (Gardiner 1985), and ξ(t) is a Gaussian random force. If we chose   $$B(f) = \left( \frac{2\pi f_s}{t_0} \right) \beta ^2 \frac{\left( 1 - f/f_s \right) (f/f_s)^\nu }{1 + \varepsilon (f/f_s)^{\nu +2}} \,,$$ (2)and   $$G = (4\pi /t_0)^{1/2} \beta \, \Delta f_s \,,$$ (3)where β, ε and ν are dimensionless parameters, the threshold values grow with time as shown in Fig. 1. For ε = 0 and ν = 0, the stationary solution of the Langevin equation leads to a distribution of thresholds Pc(fsi) which is a Gaussian centred at fs with half-width Δfs, while for ε > 0, the stationary distribution of thresholds has a power-law tail. The factor (f/fs)ν introduces a ‘delay’ in contact formation as demonstrated in Fig. 1 (inset). A newborn contact is very weak and initially its threshold grows faster than in the asymptotic limit. Figure 1. View largeDownload slide Aging of macrocontacts: growth of the average contact threshold 〈fsi〉 with its lifetime for β = 100, $$\varepsilon =2/\Delta f_{s}^{2}$$, ν = 0 (blue dashed) and ν = 2 (solid curve), fs = 1, Δfs = 0.3 and N = 3001. The inset shows the short-time behaviour in log–log scale (N = 100 001). Figure 1. View largeDownload slide Aging of macrocontacts: growth of the average contact threshold 〈fsi〉 with its lifetime for β = 100, $$\varepsilon =2/\Delta f_{s}^{2}$$, ν = 0 (blue dashed) and ν = 2 (solid curve), fs = 1, Δfs = 0.3 and N = 3001. The inset shows the short-time behaviour in log–log scale (N = 100 001). The choice ε(Δfs/fs)2 ∼ 1 leads to the GR law with the exponent b ∼ 1 (Braun & Peyrard 2013). Note that an aging of the thresholds is a necessary condition for the existence of a stick-slip behaviour of a sliding contact (Braun 2015). For spring-block models, the hypothesis made on the properties of the contacts are crucial to determine the behaviour of the model. While few studies were devoted to foreshocks, many model developments were derived from the goal of reaching a proper description of aftershocks. It was soon realized that contact aging is important. Since the pioneering work of Dieterich (1972), this is one of the most accepted mechanisms of aftershocks. The first approach was to explicitly introduce a time-dependent friction coefficient using an expression deduced from experimental observations (Dieterich 1979a,b). This approach further evolved in a rate- and state-dependent representation of the fault constitutive properties (Ruina 1983; Marone 1988; Dieterich 1994). In these approaches, the friction coefficient depends on a state variable θ and the model is completed by a deterministic differential equation which determines the time dependence of the state variable. The meaning of this state variable depends on the model. It can be the time of contact or velocity (Marone et al.1991; Dieterich 1994), but it can also be the local slip (Ruina 1983). Recent studies testing these two approaches (Helmstetter & Shaw 2009) have shown that both evolution laws produce qualitatively similar behaviours (afterslips, slow earthquakes and aftershocks) but that the slip law is more unstable than the aging law. Velocity strengthening has also been involved to explain some features of post-seismic slip (Perfettini & Ampuero 2008). Introducing a state variable appeared therefore as an essential step to describe the experimental observations on earthquakes, in particular the aftershocks and their distribution. However, the diversity of the proposals for the evolution of the state variable shows that the choice of the best description is still an open problem. Our approach proceeds along the same idea that a contact has internal degrees of freedom, with two major differences (Braun & Peyrard 2013). First, the internal state of a contact evolves according to a stochastic equation (eq. 1), instead of a deterministic differential equation. This takes into account the influence of external phenomena, such as for instance some tremors in the earth crust that could come from far away, and it models the fluctuations at the scale of the local contacts. Second, the description of a contact is based from a underlying physical model at a smaller scale. The idea is that, as explained above, at the scale of a fault a contact is actually a macrocontact which is the result of a large number of local contacts. Within this viewpoint, the laws of friction can be established from a master equation describing the breaking and reforming of the many microcontacts (Braun & Peyrard 2010). The interest of this approach that it provides a basis to include the physical properties of the local contacts, for instance the formation of chemical bonds, or local plasticity, in the properties of the macrocontacts which enter in the spring-block model, instead of postulating an equation for the state variable. Second, we incorporate the elasticity of the slider. The standard BK model assumes that, when a contact breaks, the released stress is arbitrarily redistributed among neighbouring contacts. Actually, this redistribution is due to the elastic interaction between the contacts and depends on the 3-D deformation of the slider, which cannot be rigorously treated, neither analytically nor numerically, unless the full geometry of the fault is taken into account. Instead, here we use the 1-D model, where the slider is divided in two layers as proposed in Braun & Scheibert (2014) and Braun & Tosatti (2014), which gives a reasonable approximation of the exact behaviour. In this model shown schematically in Fig. 2, the bottom layer is divided into rigid $$\lambda _c^3$$-cubes coupled together by springs of rigidity K = Eλc and coupled by frictional springs of rigidity k with the base. It describes the ‘interface’ layer (IL). The other part of the slider, the upper layer (UL), is divided into a chain of parallelepipeds of height NLλc coupled by springs KL = NLEλc. The UL and IL are coupled by the set of N transverse springs KT = Eλc/[2(1 + σP)NL], where σP is the slider Poisson ratio. Note that the geometry of the blocks is introduced to allow us to make a link between the elastic constants of the model (K, KT and KL) and the elastic parameters of the sliding medium. It is not an essential feature of the model. In this model, if we push the most-left block of the UL until the most-left frictional contact breaks, a crack emerges and propagates through the interface for a finite distance Λ ≫ λc till it is arrested in qualitative agreement with experiments (Rubinstein et al.2004, 2007, 2011). The UL plays the role of a reservoir, where the elastic energy is stored and partially released at the main shock, while the unreleased part of the elastic energy results in aftershocks satisfying the Omori law (Braun & Tosatti 2014). Figure 2. View largeDownload slide The model. The upper layer (UL) is split in rigid blocks of size λc × λc × NLλc connected by springs of elastic constant KL. The interface layer (IL) is split in rigid blocks of size λc × λc × λc connected by springs of elastic constant K. The UL and IL are coupled by springs of elastic constant KT. The IL is connected with the rigid bottom block (the base) by contacts, represented by ‘frictional’ springs of elastic constant k, which break when the local stress exceeds a threshold value. The UL is driven with the velocity v through springs of elastic constant Kd. Figure 2. View largeDownload slide The model. The upper layer (UL) is split in rigid blocks of size λc × λc × NLλc connected by springs of elastic constant KL. The interface layer (IL) is split in rigid blocks of size λc × λc × λc connected by springs of elastic constant K. The UL and IL are coupled by springs of elastic constant KT. The IL is connected with the rigid bottom block (the base) by contacts, represented by ‘frictional’ springs of elastic constant k, which break when the local stress exceeds a threshold value. The UL is driven with the velocity v through springs of elastic constant Kd. Finally, we attach springs Kd to the top boundary of the UL and drive the ends of these springs with a velocity v to simulate the driving of the seismic fault. For the simulations, we use periodic boundary conditions in the direction of the chain of blocks. 3 SIMULATION RESULTS The simulation results are presented in Figs 3–6. For the chosen set of parameters, we ran 20 independent calculations using different sets of random numbers, and with the protocol described in the Material and method section we obtained 10.8 × 106 shocks above the background level and analysed 7012 main shocks which occurred during the time interval T = 0.67 tmax ≈ 1.34 × 106 after the first 33 per cent of data were discarded in each simulation. Figure 3. View largeDownload slide Time evolution of the system: (a) the frictional force F(t), and (b) the (global) amplitude of EQs $${\cal A} (t)$$ versus time. Figure 3. View largeDownload slide Time evolution of the system: (a) the frictional force F(t), and (b) the (global) amplitude of EQs $${\cal A} (t)$$ versus time. Figure 4. View largeDownload slide Statistics of EQ magnitudes presented in Fig. 3 showing the Gutenberg–Richter power-law behaviour; the dashed line corresponds to the exponent b = 1. Figure 4. View largeDownload slide Statistics of EQ magnitudes presented in Fig. 3 showing the Gutenberg–Richter power-law behaviour; the dashed line corresponds to the exponent b = 1. Figure 5. View largeDownload slide Typical earthquakes in the model. Top panel: the EQ amplitude $${\cal A} (t)$$ versus time for the last 5 per cent of simulation time. Middle panel: detailed $${\cal A} (t)$$ dependences for two time intervals. Bottom panel: the corresponding colour maps of the EQ amplitude on the (t, x) plane, using the colour scale indicated on the right. F—foreshocks, C—calm, M—main shock and A—aftershocks. Figure 5. View largeDownload slide Typical earthquakes in the model. Top panel: the EQ amplitude $${\cal A} (t)$$ versus time for the last 5 per cent of simulation time. Middle panel: detailed $${\cal A} (t)$$ dependences for two time intervals. Bottom panel: the corresponding colour maps of the EQ amplitude on the (t, x) plane, using the colour scale indicated on the right. F—foreshocks, C—calm, M—main shock and A—aftershocks. Figure 6. View largeDownload slide Foreshocks and aftershocks statistics in the presence of contact aging (β = 100): the rate of fore- and aftershocks n(τ) (red solid symbols) and their associated standard deviations (red bars), and the shock amplitudes $${\cal A} (\tau )$$ (blue squares) relative to the corresponding main shock. The horizontal dash lines show the average magnitudes of the foreshocks and aftershocks. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. The results plotted on this figure have been obtained from the analysis of 7012 main shocks, collected over 20 independent simulations. Figure 6. View largeDownload slide Foreshocks and aftershocks statistics in the presence of contact aging (β = 100): the rate of fore- and aftershocks n(τ) (red solid symbols) and their associated standard deviations (red bars), and the shock amplitudes $${\cal A} (\tau )$$ (blue squares) relative to the corresponding main shock. The horizontal dash lines show the average magnitudes of the foreshocks and aftershocks. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. The results plotted on this figure have been obtained from the analysis of 7012 main shocks, collected over 20 independent simulations. The dependences of the driving force F(t) (Fig. 3a) and the amplitude of shocks $${\cal A}(t)$$ (Fig. 3b) on a large timescale both look rather stochastic. As we observed with shock visualization during simulation, and as can also be guessed from comparison of Figs 1 and 3(a), the frictional force is mainly determined by one or few large thresholds with fsi ≫ fs that survived for a long enough time; the breaking of these contacts results in MEQs. The statistics of the shock amplitudes follows the GR power law (see Fig. 2) for an interval of magnitudes $$\Delta {\cal {M}} \gtrsim 2$$, which is however limited by the largest possible magnitude determined by the ratio of the aging rate β and the driving velocity v, $${\cal {M}}_{\rm max} \propto {\rm log}_{10} (\beta /\sqrt{v})$$ (Braun & Peyrard 2013) (as we use here a rather small quasi-1-D system, we cannot expect a very large value of $$\Delta {\cal {M}}$$). Fig. 5 shows the $${\cal A}(t)$$ dependence for the last 5 per cent of the simulation time (top panel) and more detailed pictures at a finer timescale of two typical MEQs (middle panel); the bottom panel of Fig. 5 shows the colour maps of the earthquake amplitude. Fig. 6 presents the main result of the data analysis, the average number of shocks n(τ) in an interval δτ ( δτ = 0.0648), versus τ the rescaled time interval from (τ < 0), or after (τ > 0), the corresponding main shock. We show both the mean value n(τ), computed for the 7012 analysed shocks, and its fluctuations for different MEQs, measured by its standard deviation over all analysed MEQs. The variation of the fluctuations versus τ are the most interesting. Before the MEQ, they are large, and vary widely with τ in the range τ ≲ −0.7, showing a significant and random foreshock activity. Then suddenly they drop to a small value and stay so for the last period before the MEQ, −0.7 ≲ τ < 0. In this range, none of the 7012 analysed MEQs shows any notable foreshock activity. As expected, after the MEQ, Fig. 6 shows a large aftershock activity. Fig. 7 shows that the sharp drop of foreshock activity before large events is clearly due to the aging of the contacts, described by the parameter β in eq. (2). When β is reduced from β = 100 to 3 (β = 0 corresponding to no aging at all), all other parameters being preserved, the calm period disappears and a slight increase of activity immediately before the MEQ is even observed. We can also note that, in the absence of aging, the fluctuations of the foreshock and aftershock activity with time, are much smaller. Similarly the standard deviations of n(τ) between the different main shocks are smaller and almost time-independent. The joint evolution of the aftershocks and foreshocks when β varies shows some similarity with the correlation between the rates of foreshocks and aftershocks found in earthquakes catalogues in different geographic regions (Lippiello et al.2017). Figure 7. View largeDownload slide Same as Fig. 6 with β = 3, that is, with very little aging of the contacts. In this case, the magnitudes of the main events are not as large as in the presence of aging, so that we define the main shocks by $${\cal A}_{\rm MEQ} = \kappa \, {\cal A}_{\rm max}$$, with κ = 0.4. A set of 1205 MEQs was used to draw this figure. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. Figure 7. View largeDownload slide Same as Fig. 6 with β = 3, that is, with very little aging of the contacts. In this case, the magnitudes of the main events are not as large as in the presence of aging, so that we define the main shocks by $${\cal A}_{\rm MEQ} = \kappa \, {\cal A}_{\rm max}$$, with κ = 0.4. A set of 1205 MEQs was used to draw this figure. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. Relation with the actual scales of earthquakes. To couple our dimensionless units with real ones, we use the following estimation (Helmstetter et al.2004): let us calculate the average time between the MEQs as trec = T/nMEQ ≈ 2.8 × 103. Then, if we identify this time span with, e.g. 100 d (which is the recurrence time of $${\cal M} > {\cal M}_0 = 5$$ earthquakes in Southern California), we obtain the correspondence between time units in our model and in real seismicity:   \begin{equation*} t_0 = 1 \longleftrightarrow \tau _0 \equiv t_0/\alpha \approx 3.8 \times 10^{-3} \longleftrightarrow t_{\rm real} \sim 1 \, {\rm hr}. \end{equation*} Thus, the time span in Fig. 3 corresponds to about 230 yr and the average calm period before the MEQs, τcalm ≲ 0.7 observed in Fig. 6, corresponds to about 7 d. However, this estimation depends on the value of $${\cal M}_0$$ we choose. For example, at the Parkfield segment along the San Andreas fault, California, interplate earthquakes of magnitude about $${\cal M}_0 = 6$$ have occurred at recurrence intervals of 23 ± 9 yr (Kawamura et al.2012), while along the Nankai trough, where the Philippine Sea plate subducts beneath southwestern Japan, great earthquakes of magnitude $${\cal M}_0 = 8$$ have repeatedly occurred every 100 yr (Kawamura et al.2012). If earthquakes in a given region follow the GR law, then the recurrence period depends on $${\cal M}_0$$ as $$t_{\rm rec} \propto \exp ({\cal M}_0)$$. In a simple model, we cannot provide an unquestionable value of the magnitude that could be compared with actual magnitudes, but the main difficulty for a comparison with actual data is that the limited size of the model restricts the range of the observable events. Fig. 6 shows events with a range that spans less than three orders of magnitude, while actual foreshocks can be much smaller with respect to an MEQ. The model cannot show the smallest events and tends to overestimate the duration of the calm period. In actual earthquakes, it may sometimes only last for a few hours or less than a day as in Raleigh et al. (1977), but may extend to much longer periods reaching several months (Main & Meredith 1991). 4 DISCUSSION AND CONCLUSION In order to correctly describe various features of earthquakes, the model that we used in this work is more complex than the simplified case that we presented in the Introduction, and the probability distribution of the breaking thresholds is not imposed a priori but instead results from the dynamics of the model through eq. (1). Fig. 8 shows that the qualitative mechanism for the origin of a calm period that we presented in the introduction is nevertheless also present in this more elaborate case. While the initial probability distribution of the thresholds for contact breaking Pc(fsi) was a Gaussian, the evolution of the thresholds according to eq. (1) leads to a very different shape as time evolves. The maximum value of fsi grows, and moreover the distribution tends to split into a large cluster of contacts which break easily (the ‘weak contacts’ of the qualitative model) and a few very strong contacts. The contacts with intermediate thresholds tend to disappear. It is the few strong contacts that prevent the sliding and are responsible for the MEQs when they finally break. The gap between these strong contacts and the many weak contacts is responsible for the calm period. This evolution is mostly governed by the aging of the contacts. Figure 8. View largeDownload slide Evolution of the probability distribution of the contact thresholds as time grows (from top to bottom figure, as indicated in each frame. Figure 8. View largeDownload slide Evolution of the probability distribution of the contact thresholds as time grows (from top to bottom figure, as indicated in each frame. In natural seismicity quiescence is observed in about 40 per cent of events (Turcotte et al.2000; Kawamura et al.2012): the frequency of small events is gradually enhanced preceding the main shock, whereas, just before the main shock, it is suppressed in a close vicinity of the epicentre of the upcoming event. In about 60 per cent of natural earthquakes another scenario is observed. Large earthquakes are preceded by a period during which the surrounding region experiences a phase of accelerated seismic release (Ben-Zion & Lyakhovsky 2002; McGuire et al.2005). The rate of these foreshocks also follows the Omori law (Kagan & Knopoff 1978; Jones & Molnar 1979), usually with a smaller value of tO. Our model may also demonstrate such a behaviour for smaller values of the ratio β2/v (Braun & Tosatti 2014), when the distribution of large thresholds is more uniform. However, in the model as in reality, earthquakes exhibit a broad dispersion of their features from event to event. This is why seismic quiescence is notable on the standard deviations of the foreshock activity much more than on the mean value of n(τ) averaged over a large number of MEQs (inset of Fig. 6). The intensity of the foreshock activity prior to the calm period varies strongly from event to event and can even stay rather low for some MEQs. This means that using quiescence as a warning signal may be unreliable even for a fault which have the characteristics which lead to quiescence. In nature large earthquakes are almost always followed by a very large seismicity rate (aftershocks). In our model aftershocks are less likely for parameters that lead to a calm period before big events (see Figs 6 and 7). This is consistent with the picture shown in Fig. 8 because the quiescence is associated to an exhaustion of the weak contacts before the MEQ so that the large event is more likely to fully release the stress in the fault. It would be interesting to check whether such an effect is actually observed for natural earthquakes showing quiescence. However, the decrease of aftershock activity may be exaggerated in our calculations due to the restricted size of the model. In an actual earthquake a complete fault line does not break at once, and instead the partial breaking tends to accumulate stress in some regions, as our model does for some parameter sets. This study points out the crucial role of the distribution of thresholds to control the pattern of foreshocks. The existence of a few strong contacts tends to induce a period of seismic quiescence before a large earthquake. Such a distribution could of course come from the geometry of a fault. Some studies correlated different foreshock patterns to general features of faults, such as interplate or intraplate earthquakes (Bouchon et al.2013). However, one important result of our study is that the aging of the contacts can also be a determining factor, in slowly building a multipeaked distribution of contact thresholds. Several mechanisms may be responsible for aging of the sliding interface, such as removing of dislocations from asperities, growth of contact sizes due to plastic deformation, squeezing of a ‘lubricant’ (e.g. water or powdered rock) from the interface or its solidification due to high pressure, coalescence of contacts, etc. These features may be characteristic of a particular fault or a type of rock found in some area, so that statistics of quiescence could differ from place to place. As the timescale of quiescence is short compared to the periods generally involved in foreshock studies, there are very few available data which could support this idea. However, a recent paper (Lippiello et al.2017) presents (in fig. 3) statistics of foreshocks which show some drop of activity before a main shock for the Northern California Earthquake Catalog (Waldhauser & Schaff 2008; particularly for magnitude 2) which does not appear for other catalogues. This could suggest that such local differences exist although this is to be taken with caution because the study did not pay a particular attention to quiescence. A recent study has stressed the dominant role of the formation of interfacial chemical bonds (Qunyang Li et al.2011). It is interesting that the chemistry is also proposed as a cause for seismic quiescence (Main & Meredith 1991). It would be interesting if geologists could correlate the chemical properties of the rocks involved in the faults and the pattern of foreshocks. Real earthquakes, as well as the model that we have presented when its parameters are changed, can exhibit different patterns of foreshocks. An acceleration of the events prior to a large earthquake is not always the rule. Seismic quiescence may have deep implications by showing that one should not take too strictly the viewpoint that earthquakes are associated with criticality. The earthquake may not be the divergence of some accelerating events that culminate in a main shock, but instead appear after some unusually quiet period. Our model of aging shows that such a behaviour may emerge from what could, at first glance, be considered as a secondary effect, the aging of the contacts, which in turn may strongly alter the distribution of the thresholds at which contacts break. Of course one has to be careful in extending conclusions drawn from a simple, 1-D model to the complexity of natural earthquakes. As discussed above, we hope that such a model can nevertheless give useful hints for further studies of real earthquakes. But it is also obvious that developments of the model are clearly needed, in particular, generalizations to 2-D models of the interface and 3-D models of the elastic slider which would improve the GR statistics and provide a power law for earthquake spatial distribution. The process describing the aging of the contacts also certainly needs further investigations as it is presently largely arbitrary. We use a stochastic process that goes to a Gaussian distribution in one limit and a power law in another, but how this process can be generated by the dynamics of ‘macrocontacts’ is still open. Acknowledgements This work was partially supported by CNRS-Ukraine PICS grant no. 151539. OB was supported in part by COST Action MP1303. REFERENCES Ben-Zion Y., Lyakhovsky V., 2002. Accelerated seismic release and related aspects of seismicity patterns on earthquake faults, Pure appl. Geophys . 159, 2385– 2412. Google Scholar CrossRef Search ADS   Bouchon M., Durand V., Marsan D., Karabulut H., Schmittbuhl J., 2013. The long precursory phase of most large interplate earthquakes, Nature Geosci ., 6, 299– 302. Google Scholar CrossRef Search ADS   Bowman D.D., King G.C.P., 2001. Accelerating seismicity and stress accumulation before large earthquakes, Geophys. Res. Lett ., 28, 4039– 4042. Google Scholar CrossRef Search ADS   Braun O.M., 2015. Stick-slip vs. smooth sliding in the multicontact interface, Europhys. Lett ., 109, 48004-1– 48004-6. Google Scholar CrossRef Search ADS   Braun O.M., Peyrard M., 2013 Role of aging in a minimal model of earthquakes, Phys. Rev. E , 87, 032808-1–032808-7. Braun O.M., Scheibert J., 2014. Propagation length of self-healing slip pulses at the onset of sliding: a toy model, Tribol. Lett ., 56, 553– 562. Google Scholar CrossRef Search ADS   Braun O.M., Tosatti E., 2014. Aftershocks in a frictional earthquake model, Phys. Rev. E , 90, 032403-1– 032403-8. Google Scholar CrossRef Search ADS   Braun O.M., Peyrard M., Stryzheus D.V., Tosatti E., 2012. Collective effects at frictional interfaces, Tribol. Lett ., 48, 11– 24. Google Scholar CrossRef Search ADS   Braun O.M., Peyrard M., 2010. Master equation approach to friction at the mesoscale, Phys. Rev. E , 82, 036117-1– 036117-19. Google Scholar CrossRef Search ADS   Burridge R., Knopoff L., 1967. Model and theoretical seismicity, Bull. seism. Soc. Am ., 57, 341– 371. Caroli C., Nozieres Ph., 1998. Hysteresis and elastic interactions of microasperities in dry friction, Eur. Phys. J. B , 4, 233– 246. Google Scholar CrossRef Search ADS   Dieterich J.H., 1972. Time dependent friction as a possible mechanism for afterschocks, J. geophys. Res ., 77, 3771– 3781. Google Scholar CrossRef Search ADS   Dieterich J.H., 1979a. Modeling of rock friction. 1. Experimental results and constitutive equations, J. geophys. Res ., 84, 2161– 2168. Google Scholar CrossRef Search ADS   Dieterich J.H., 1979b. Modeling of rock friction. 2. Simulation of preseismic slip, J. geophys. Res ., 84, 2169– 2175. Google Scholar CrossRef Search ADS   Dieterich J.H., 1994. A constitutive law for the rate of earthquaque production and its application to earthquake clustering, J. geophys. Res ., 90, 2601– 2618. Google Scholar CrossRef Search ADS   Gardiner C.W., 1985. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences , Springer-Verlag, Berlin. Google Scholar CrossRef Search ADS   Grassberger P., 1994. Efficient large-scale simulations of a uniformly driven system, Phys. Rev. E , 49, 2436– 2444. Google Scholar CrossRef Search ADS   Gutenberg B., Richter C.F., 1954. Earthquake magnitude, intensity, energy and acceleration, Bull. seism. Soc. Am ., 46, 105– 146. Hainzl S., Zöller G., Kurths J., 2000a. Similar power laws for foreshock and aftershock sequences in a spring block model for earthquakes, J. geophys. Res ., 104, 7243– 7253. Google Scholar CrossRef Search ADS   Hainzl S., Zöller G., Kurths J., 2000b. Self-organization of spatio-temporal earthquake clusters, Nonlin. Proces. Geophys ., 7, 21– 29. Google Scholar CrossRef Search ADS   Helmstetter A., Hergarten S., Sornette D., 2004. Properties of foreshocks and aftershocks of the nonconservative self-organized critical Olami-Feder-Christensen model, Phys. Rev. E , 70, 046120-1– 046120-13. Google Scholar CrossRef Search ADS   Helmstetter A., Shaw B.E., 2009. Afterslip and aftershocks in the rate-and-state friction law, J. geophys. Res. , 114, B01308-1– B01308-24. Google Scholar CrossRef Search ADS   Hergarten S., Krenn R., 2011. Synchronization and desynchronization in the Olami-Feder-Christensen earthquake model and potential implications for real seismicity, Nonlin. Process. Geophys ., 18, 635– 642. Google Scholar CrossRef Search ADS   Jagla E.A., 2010. Realistic spatial and temporal earthquake distributions in a modified Olami-Feder-Christensen model, Phys. Rev. E , 81, 046117-1– 046117-8. Google Scholar CrossRef Search ADS   Jagla E.A., Kolton A.B., 2010. The mechanisms of spatial and temporal earthquake clustering, J. geophys. Res ., 115, B05312-1– B05312-8. Google Scholar CrossRef Search ADS   Jagla E.A., 2013. Forest-fire analogy to explain the b value of the Gutenberg-Richter law for earthquakes, Phys. Rev. Lett ., 111, 238501-1– 238501-5. Google Scholar CrossRef Search ADS   Jones L.M., Molnar P., 1979. Some characteristics of foreshocks and their possible relationship to earthquake prediction and premonitory slip on faults, J. geophys. Res ., 84, 3596– 3608. Google Scholar CrossRef Search ADS   Jordan T.H., Jones L.M., 2010. Operational earthquake forecasting: Some thoughts on why and how, Seismol. Res. Lett ., 81, 571– 574. Google Scholar CrossRef Search ADS   Kagan Y.Y., Knopoff L., 1978. Statistical study of the occurrence of shallow earthquakes, Geophys. J. R. astr. Soc ., 55, 67– 86. Google Scholar CrossRef Search ADS   Kawamura H., Hatano T., Kato N., Biswas S., Chakrabarti B.K., 2012. Statistical physics of fracture, friction, and earthquakes, Rev. Mod. Phys ., 84, 839– 884. Google Scholar CrossRef Search ADS   Kazemian J., Tiampo K.F., Klein W., Dominguez R., 2015. Foreshock and aftershocks in simple earthquake models, Phys. Rev. Lett ., 114, 088501-1– 088501-6. Google Scholar CrossRef Search ADS   Lippiello E., Giacco F., Marzocchi W., Godano G., de Archangelis L., 2017. Statistical features of foreshocks in instrumental and ETAS catalogs, Pure appl. Geophys ., 174, 1679– 1697. Google Scholar CrossRef Search ADS   Main I.G., Meredith P.G., 1991. Stress corrosion constitutive laws as a possible mechanism for intermediate-term and short-term seismic quiescence, Geophys. J. Int ., 107, 363– 372. Google Scholar CrossRef Search ADS   Marone C., 1988. Laboratory-derived friction laws and their application to seismic faulting, Ann. Rev. Earth. Planet. Sci ., 26, 643– 696. Google Scholar CrossRef Search ADS   Marone C., Scholz C.H., Bilham R., 1991. On the mechanics of earthquake afterslip, J. geophys. Res ., 96 8441– 8452. Google Scholar CrossRef Search ADS   McGuire J.J., Boettcher M.S., Jordan T.H., 2005. Foreshock sequences and short-term earthquake predictability on East Pacific Rise transform faults, Nature , 434, 457– 461. Google Scholar CrossRef Search ADS PubMed  Mega M.S., Allegrini P., Grigolini P., Latora V., Palatella L., Rapisarda A., Vinciguerra S., 2003. Power-law time distribution of large earthquakes, Phys. Rev. Lett ., 90, 188501-1– 188501-4. Google Scholar CrossRef Search ADS   Olami Z., Feder H.J.S., Christensen K., 1992. Self-organized criticality in a continuous, nonconservative cellular automaton modeling earthquakes, Phys. Rev. Lett ., 68, 1244– 1247. Google Scholar CrossRef Search ADS PubMed  Omori F., 1894. On after-shocks of earthquakes, J. Coll. Sci. Imp. Univ. Tokyo  7, 111– 120. Pelletier J.D., 2000. Spring-block models of seismicity: review and analysis of a structurally-heterogeneous model coupled to a viscous asthenosphere, Geophys. Monogr ., 120, 27– 42. Perfettini H., Ampuero J.-P., 2008. Dynamics of a velocity strengthening fault region: implications for slow earthquakes and postseismic slip, J. geophys. Res.  113, B09411-1– B09411-22. Google Scholar CrossRef Search ADS   Qunyang L., Tullis T.E., Goldsby D., Carpick R.W., 2011. Frictional ageing from interfacial bonding and the origins of rate and state friction, Nature , 480, 233– 237. Google Scholar CrossRef Search ADS PubMed  Raleigh C.B. et al.  , 1977. Prediction of the Haicheng earthquake, EOS, Trans. Am. geophys. Un ., 58, 236– 272. Google Scholar CrossRef Search ADS   Rice J.R., Ben-Zion Y., 1996. Slip complexity in earthquake fault models, Proc. Natl. Acad. Sci ., 93, 3811– 3818. Google Scholar CrossRef Search ADS   Rubinstein S.M., Cohen G., Fineberg J., 2004. Detachment fronts and the onset of dynamic friction, Nature , 430, 1005– 1009. Google Scholar CrossRef Search ADS PubMed  Rubinstein S.M., Cohen G., Fineberg J., 2007. Dynamics of precursors to frictional sliding, Phys. Rev. Lett.  98, 226103-1– 226103-4. Google Scholar CrossRef Search ADS   Rubinstein S.M., Barel I., Reches Z., Braun O.M., Urbakh M., Fineberg J., 2011. Slip sequences in laboratory experiments resulting from inhomogeneous shear as analogs of earthquakes associated with a fault edge, Pure appl. Geophys ., 168, 2151– 2166. Google Scholar CrossRef Search ADS   Ruina A., 1983. Slip instability and state variable friction laws, J. geophys. Res ., 88, 10 359–10 370. Rundle J.B., Klein W., 1995. Dynamical segmentation and rupture patterns in a “toy” slider-block model for earthquakes, Nonlin. Process. Geophys ., 2, 61– 79. Google Scholar CrossRef Search ADS   Serino C.A., Tiampo K.F., Klein W., 2011. New approach to Gutenberg-Richter scaling, Phys. Rev. Lett ., 106, 108501-1– 108501-4. Google Scholar CrossRef Search ADS   Socolar J.E.S., Grinstein G., Jayaprakash C., 1993. On self-organized criticality in nonconserving systems, Phys. Rev. E , 47, 2366– 2376. Google Scholar CrossRef Search ADS   Stein R.S., 1999. The role of stress transfer in earthquake occurrence, Nature , 402, 605– 609. Google Scholar CrossRef Search ADS   Turcotte D.L., Shcherbakov R., Rundle J.B., 2000. Complexity and earthquakes, in Treatise on Geophysics , 4 Chap. 23, pp. 675– 700, ed. Kanamori H., Elsevier. Utsu T., Ogata Y., Matsu’ura R.S., 1995. The centenary of the Omori formula for a decay law of aftershock activity, J. Phys. Earth , 43, 1– 33. Google Scholar CrossRef Search ADS   Waldhauser F., Schaff D.P., 2008. Large-scale relocation of two decades of Northern California seismicity using cross-correlation and double-difference methods, J. Geophys. Res.: Solid Earth , 113 (B8) B08311, doi:10.1029/2007JB005479. APPENDIX: METHODS A1 Parameters We use dimensionless units, where four model parameters are set to 1, k = 1, λc = 1, t0 = 1 and fs = 1 so that the other parameters should be compared with those. The IL consists of N = 101 $$\lambda _c^3$$-cubes, each of mass m = 10−4, as chosen in our analysis of stick-slip phenomena (Braun 2015). The springs between the cubes have the elastic constant K = 100. The UL consists of parallelepipeds of height NL = 25 (so that their masses are NLm); for the Poisson ratio, we took σP = 0.3. With these parameters, the crack propagation distance (where aftershocks are expected) is Λ ∼ 30 (Braun & Scheibert 2014). For the rate of aging, we took β = 100, and for the parameter ε we used $$\varepsilon =2 / \Delta f_s^2$$. With these parameters, the statistics of the magnitudes of earthquakes follow the GR law for an interval of magnitudes $$\Delta {\cal M} \sim 2$$, which is of course limited by the size of the model. For the dispersion of thresholds, we took Δfs = 0.3, and used ν = 2 to mimic a delay in contact formation. The initial stretching of the contacts as well as the initial thresholds of the newborn contacts, when they form again after breaking, are taken from a Gaussian distribution with zero mean and dispersion Δfs. The UL is moved through the springs with elastic constant Kd = 0.5 attached to its top and driven with the velocity v = 0.1. The equations of motion are solved by a Runge–Kutta method with a viscous damping η = 0.3 ω0, where ω0 = (K/m)1/2, so that the dynamics of the system is underdamped. With these parameters, the motion of a single $$\lambda _c^3$$-block exhibits a stick-slip behaviour [recall that stick-slip exists only for a certain interval of driving velocities (Braun 2015)]. The duration of the runs is tmax = 2 × 106. The first 33 per cent of a trajectory are discarded to allow the system to lose memory from the initial configuration and reach a steady state (dashed vertical line in Fig. 3). A2 Protocol and data analysis In a simulation, we calculate and store the driving force F(t) and the dimensionless quake amplitude defined as the sum of the force drops due to broken contacts at every time step,   $${\cal A} (t) = \sum _{i=1}^N \Delta F_i (t)/( N f_s ) \,.$$ (A1)The rate of shocks, n(t), is defined as the number of shocks per time unit regardless of their amplitude. The dimensionless earthquake magnitude is defined as   $${\cal M} = \frac{2}{3} \log _{10} \left( {\cal A}/{\cal A}_0 \right) ,$$ (A2)where $${\cal A}_0$$ is a constant which defines the magnitude scale. As for the study of real earthquakes (Bouchon et al.2013), the analysis of foreshocks is delicate because one has to distinguish ‘main shocks’ in the middle of many events, occurring randomly and with a broad range of amplitudes. To be qualified as ‘foreshocks’ events must be related to the main shock both in time and space, and moreover a large foreshock should not be considered as an independent MEQ. In analysing the results, we used the following protocol (Braun & Tosatti 2014). First we remove small ‘background’ earthquakes with amplitudes below some level from the analysis. We only retain $${\cal A} > 2 \, \langle {\cal A} (t) \rangle$$ (see broken red line in Fig. 3). Next, we single out the main shocks above some level $${\cal A}_{\rm MEQ}$$. We took here $${\cal A}_{\rm MEQ} = \kappa \, {\cal A}_{\rm max}$$, with κ = 0.03 (solid red line in Fig. 3). We also used the following rescaling procedure to detect large events which could be related. If nMEQ is the total number of MEQs in a simulation, let us call σ = S/nMEQ, the average area occupied by a single MEQ in the (t, x) plane, with S = Nλctmax. We then rescale the time coordinate t → τ = t/α with α = σ/λ2, where λ ∼ Λ is some distance chosen in such a way that the distribution of MEQs on the (τ, x) plane becomes isotropic (λ ≈ 30 and α ≈ 300 for the parameters used in Fig. 3). In this rescaled space, events which are related to a particular MEQ lie within a circle around this MEQ so that they can be easily identified. We can now scan all MEQ coordinates on the (τ, x) plane and, if the distance ρij = [(τi − τj)2 + (xi − xj)2]1/2 between two MEQs i and j is smaller than some value ρcut then only the larger of these two MEQs remains as the MEQ, while the lower one is removed from the list of MEQs. We chose ρcut = 0.4λ ∼ 0.4Λ, that is, to be considered as related to each other two events must not be separated by more than about 0.4 × the average distance along which a crack propagates. With this protocol, we have obtained a set of well-separated MEQs isotropically occupying the (τ, x) plane, and we may calculate the temporal and spatial distributions of all earthquakes within some area around each MEQ. We consider that all the events separated from the corresponding MEQ by less than ρ0 = ρcut/3 are related to this particular MEQ. They can be either foreshocks or aftershocks. However, we have to make sure that some events that occur before one of the MEQ do not belong to the tail of aftershocks following a violent MEQ which occurred earlier. This is done by examining only MEQs which are separated by a time interval long enough to allow the tail of aftershocks to die out. In practice, only MEQs separated by a τ interval greater than 0.4 × the full τ interval that we scan around MEQs are analysed. This lead us to discard more than 30 per cent of the events that would otherwise qualify as MEQs. Then we collapse all data together, designating τ = 0 for every main shock and normalizing shocks amplitudes on the corresponding main shock value. We compute n(τ) the number of secondary events (foreshocks or aftershocks) in a given δτ = 0.0648 interval averaged over all the investigated MEQs as well as the standard deviations of their values for the different MEQs, which provides an estimate of the error bars for n(τ) and gives us an information on the fluctuations of the earthquake activity in a given τ interval around an MEQ. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Seismic quiescence in a frictional earthquake model

, Volume 213 (1) – Apr 1, 2018
8 pages

/lp/ou_press/seismic-quiescence-in-a-frictional-earthquake-model-1q10PilFvm
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy008
Publisher site
See Article on Publisher Site

### Abstract

Summary We investigate the origin of seismic quiescence with a generalized version of the Burridge–Knopoff model for earthquakes and show that it can be generated by a multipeaked probability distribution of the thresholds at which contacts break. Such a distribution is not assumed a priori but naturally results from the aging of the contacts. We show that the model can exhibit quiescence as well as enhanced foreshock activity, depending on the value of some parameters. This provides a generic understanding for seismic quiescence, which encompasses earlier specific explanations and could provide a pathway for a classification of faults. Statistical methods, Earthquake dynamics, Theoretical seismology 1 INTRODUCTION While two main seismological laws, the Gutenberg–Richter (GR) law for the magnitude distribution of earthquakes (Gutenberg & Richter 1954) and the Omori law for the time evolution of the frequency of aftershocks (Omori 1894; Utsu et al.1995), are well-established empirically and even may be justified theoretically (Burridge & Knopoff 1967; Olami et al.1992; Socolar et al.1993; Grassberger 1994; Rundle & Klein 1995; Rice & Ben-Zion 1996; Stein 1999; Hainzl et al.2000a,b; Helmstetter et al.2004; Hergarten & Krenn 2011; Serino et al.2011; Jagla 2010; Jagla & Kolton 2010; Jagla 2013; Braun & Peyrard 2013; Mega et al.2003; Kawamura et al.2012; Pelletier 2000), our knowledge about foreshocks is still limited. If a clear pattern could be distinguished in the foreshock activity it might provide an early warning for earthquakes, but we are far from this stage because the observations reveal very different pictures. Many earthquakes are preceded by foreshocks, but their frequency can nevertheless vary widely depending on the type of earthquake (Bouchon et al.2013). Moreover, some earthquakes are preceded by an unexpected calm period, lasting for several hours or more (Jones & Molnar 1979). It is such a period of quiescence, viewed as a characteristic feature of the imminence of the main shock that allowed the only successful prediction of an earthquake (Raleigh et al.1977), which saved a large number of lives in China in 1975. Although the mechanisms that generate foreshocks are still under discussion (Bouchon et al.2013), the existence of some activity before a major earthquake (MEQ) does not seem surprising. In recent experiments with slider pushed on the side, Fineberg et. al. (Rubinstein et al.2004, 2007, 2011) showed that the global sliding, which corresponds to an MEQ, is preceded by propagation of several cracks called precursors that may be used as indicators of the upcoming MEQ. This is also consistent with some ideas that relate earthquakes to self-organized criticality or cascade mechanisms (Olami et al.1992; Socolar et al.1993; Grassberger 1994; Helmstetter et al.2004; Hergarten & Krenn 2011). The fact that a quiescence period could be a characteristic announcement of an MEQ looks more surprising. Various mechanisms have been considered to explain seismic quiescence (Hergarten & Krenn 2011; Main & Meredith 1991). Here, we propose a generic mechanism related to the distribution of the threshold for the breaking of the contacts along a fault and illustrate it by simulations of a simple earthquake model. The Burridge and Knopoff (BK) spring-block model (Burridge & Knopoff 1967) further developed in a number of works (Olami et al.1992; Hainzl et al.2000a,b; Helmstetter et al.2004; Hergarten & Krenn 2011; Serino et al.2011; Jagla 2010; Jagla & Kolton 2010; Jagla 2013; Braun & Peyrard 2013; Braun & Tosatti 2014; Kazemian et al.2015) [see also reviews (Pelletier 2000; Kawamura et al.2012), and references therein] has been used as a generic model which reproduces many features of earthquakes. Simulations showed that some generalized versions of the BK model may demonstrate the GR law and even the Omori law, but the existence of foreshocks in this model is not clearly demonstrated and explained yet (e.g. see Hainzl et al.2000a,b; Helmstetter et al.2004; Hergarten & Krenn 2011; Jagla 2010; Jagla & Kolton 2010; Jagla 2013; Kazemian et al.2015; Pelletier 2000). In BK-type models, the top block (the slider) is coupled with the bottom block (the base assumed to be fixed) by a set of frictional contacts. When the slider moves, a frictional contact i, modeled as an elastic spring, elongates so that the force at the contact point, fi = kxi (we assume that all contacts have the same rigidity k for the sake of simplicity), increases until it reaches a threshold value fsi for which the contact breaks. Then it forms again, with zero stretching xi ∼ 0, and the process can repeat. The evolution of such a system is determined by the distribution of thresholds Pc(fsi). The general features of this distribution also determines the nature of the foreshock activity, as one can realize by considering two simple examples. If Pc(fsi) is one-peaked (Braun & Peyrard 2010), e.g. Gaussian with the centre at fs and a standard deviation Δfs, then the total frictional force on N contacts, $$F = \sum _{i=1}^N f_i$$ linearly increases with the slider displacement X until the contact forces reach values fi ∼ fs − Δfs. Then some contacts start to break, one by one, giving rise to small earthquakes. The growth of F(X) slows down and turns into a decay for larger values of X, and under certain conditions (Braun 2015) an elastic instability may occur. The top block slides for a distance ΔX ≳ Δfs/k, while most of the contacts break. This event corresponds to a large earthquake. In this case, the earthquake activity increases just before the MEQ. Such a scenario is common in natural earthquakes (Kazemian et al.2015; Bowman & King 2001; Jordan & Jones 2010). If the threshold distribution is two-peaked, e.g. when the interface consists of two types of contacts, weak and strong ones with thresholds fs1 and fs2, respectively, as recently considered in Kazemian et al. (2015), the scenario is qualitatively different. To see how it works let us consider a simple case. We assume that, when a contact breaks it forms again as a weak contact with probability p and strong contact with probability 1 − p, and that the breaking of a contact reduces the friction force by φ. At time t0 = 0, we start with $$n_1^{(0)} = p N$$ weak contacts and $$n_2^{(0)} = (1-p)N$$ strong contacts. If the slider moves at velocity v the driving force grows as Nkvt. When it reaches fs1 at time t1, all weak contacts break. The friction force drops by $$\varphi n_1^{(0)}$$ and the broken contacts reform, leading to $$n_1^{(1)} = p n_1^{(0)} = p^2 N$$ weak contacts and a growing number of strong contacts. The same process can repeat again, leading to an infinite series of foreshocks, while the number of weak contacts decreases as $$n_1^{(j)} = p^{j+1} N$$. During the time interval tj − tj − 1, the force has to grow enough to reach fs1 again after the breaking $$n_1^{(j-1)}$$ weak contacts. It requires a time tj − tj − 1 = φpjN/(Nkv), so that the total time for the infinite series of foreshocks stays finite $$\sum _{j=1}^{\infty } (t_{j} - t_{j-1}) = (\varphi / kv) p / (1 - p)$$. After that time only strong contacts persist. As the slider keeps moving, there is a calm period without any foreshock until the strong threshold fs2 is reached. Then all strong contacts break suddenly, leading to an MEQ. This picture of contacts with a two-peaked distribution is oversimplified, and would not show all features of real earthquakes, such as the GR law or the Omori law for aftershocks, but it shows how the properties of the contact thresholds are sufficient to lead to seismic quiescence. Let us now consider a more realistic case and investigate its properties. 2 MODEL We consider the earthquake model previously developed in Braun & Tosatti (2014) and Braun & Scheibert (2014). It is based on the BK-type model, i.e. a set of frictional contacts between a slider and a base, but with the following two important ingredients. First, by contact we mean a macrocontact consisting of a large number of microcontacts on the area $$\lambda _c^2$$, where λc ≈ a2E/k is the elastic correlation length (Caroli & Nozieres 1998; Braun et al.2012) (here a is the average distance between the microcontacts and E is the Young modulus of the slider). For a frictional metal/metal interface λc is typically of the order of μm (Caroli & Nozieres 1998; Braun et al.2012), but for a seismic fault it may be much larger. On the distance λc, the slider may be treated as rigid, while at larger distances we have to account for the deformation of the block. Thus, in our model a contact is a macroscopic object. It has its own ‘degrees of freedom’, which may lead to its aging resulting in a growth of its threshold value as the lifetime of the contact increases. Contact aging is a complex stochastic process which can have various physical or chemical origins. We model it by a simple Langevin equation (Braun & Peyrard 2013; Braun & Tosatti 2014)   $$df_{si} (t)/dt = B(f_{si}) + G \, \xi (t) \,,$$ (1)where B(f) and G are the so-called drift and stochastic forces, respectively (Gardiner 1985), and ξ(t) is a Gaussian random force. If we chose   $$B(f) = \left( \frac{2\pi f_s}{t_0} \right) \beta ^2 \frac{\left( 1 - f/f_s \right) (f/f_s)^\nu }{1 + \varepsilon (f/f_s)^{\nu +2}} \,,$$ (2)and   $$G = (4\pi /t_0)^{1/2} \beta \, \Delta f_s \,,$$ (3)where β, ε and ν are dimensionless parameters, the threshold values grow with time as shown in Fig. 1. For ε = 0 and ν = 0, the stationary solution of the Langevin equation leads to a distribution of thresholds Pc(fsi) which is a Gaussian centred at fs with half-width Δfs, while for ε > 0, the stationary distribution of thresholds has a power-law tail. The factor (f/fs)ν introduces a ‘delay’ in contact formation as demonstrated in Fig. 1 (inset). A newborn contact is very weak and initially its threshold grows faster than in the asymptotic limit. Figure 1. View largeDownload slide Aging of macrocontacts: growth of the average contact threshold 〈fsi〉 with its lifetime for β = 100, $$\varepsilon =2/\Delta f_{s}^{2}$$, ν = 0 (blue dashed) and ν = 2 (solid curve), fs = 1, Δfs = 0.3 and N = 3001. The inset shows the short-time behaviour in log–log scale (N = 100 001). Figure 1. View largeDownload slide Aging of macrocontacts: growth of the average contact threshold 〈fsi〉 with its lifetime for β = 100, $$\varepsilon =2/\Delta f_{s}^{2}$$, ν = 0 (blue dashed) and ν = 2 (solid curve), fs = 1, Δfs = 0.3 and N = 3001. The inset shows the short-time behaviour in log–log scale (N = 100 001). The choice ε(Δfs/fs)2 ∼ 1 leads to the GR law with the exponent b ∼ 1 (Braun & Peyrard 2013). Note that an aging of the thresholds is a necessary condition for the existence of a stick-slip behaviour of a sliding contact (Braun 2015). For spring-block models, the hypothesis made on the properties of the contacts are crucial to determine the behaviour of the model. While few studies were devoted to foreshocks, many model developments were derived from the goal of reaching a proper description of aftershocks. It was soon realized that contact aging is important. Since the pioneering work of Dieterich (1972), this is one of the most accepted mechanisms of aftershocks. The first approach was to explicitly introduce a time-dependent friction coefficient using an expression deduced from experimental observations (Dieterich 1979a,b). This approach further evolved in a rate- and state-dependent representation of the fault constitutive properties (Ruina 1983; Marone 1988; Dieterich 1994). In these approaches, the friction coefficient depends on a state variable θ and the model is completed by a deterministic differential equation which determines the time dependence of the state variable. The meaning of this state variable depends on the model. It can be the time of contact or velocity (Marone et al.1991; Dieterich 1994), but it can also be the local slip (Ruina 1983). Recent studies testing these two approaches (Helmstetter & Shaw 2009) have shown that both evolution laws produce qualitatively similar behaviours (afterslips, slow earthquakes and aftershocks) but that the slip law is more unstable than the aging law. Velocity strengthening has also been involved to explain some features of post-seismic slip (Perfettini & Ampuero 2008). Introducing a state variable appeared therefore as an essential step to describe the experimental observations on earthquakes, in particular the aftershocks and their distribution. However, the diversity of the proposals for the evolution of the state variable shows that the choice of the best description is still an open problem. Our approach proceeds along the same idea that a contact has internal degrees of freedom, with two major differences (Braun & Peyrard 2013). First, the internal state of a contact evolves according to a stochastic equation (eq. 1), instead of a deterministic differential equation. This takes into account the influence of external phenomena, such as for instance some tremors in the earth crust that could come from far away, and it models the fluctuations at the scale of the local contacts. Second, the description of a contact is based from a underlying physical model at a smaller scale. The idea is that, as explained above, at the scale of a fault a contact is actually a macrocontact which is the result of a large number of local contacts. Within this viewpoint, the laws of friction can be established from a master equation describing the breaking and reforming of the many microcontacts (Braun & Peyrard 2010). The interest of this approach that it provides a basis to include the physical properties of the local contacts, for instance the formation of chemical bonds, or local plasticity, in the properties of the macrocontacts which enter in the spring-block model, instead of postulating an equation for the state variable. Second, we incorporate the elasticity of the slider. The standard BK model assumes that, when a contact breaks, the released stress is arbitrarily redistributed among neighbouring contacts. Actually, this redistribution is due to the elastic interaction between the contacts and depends on the 3-D deformation of the slider, which cannot be rigorously treated, neither analytically nor numerically, unless the full geometry of the fault is taken into account. Instead, here we use the 1-D model, where the slider is divided in two layers as proposed in Braun & Scheibert (2014) and Braun & Tosatti (2014), which gives a reasonable approximation of the exact behaviour. In this model shown schematically in Fig. 2, the bottom layer is divided into rigid $$\lambda _c^3$$-cubes coupled together by springs of rigidity K = Eλc and coupled by frictional springs of rigidity k with the base. It describes the ‘interface’ layer (IL). The other part of the slider, the upper layer (UL), is divided into a chain of parallelepipeds of height NLλc coupled by springs KL = NLEλc. The UL and IL are coupled by the set of N transverse springs KT = Eλc/[2(1 + σP)NL], where σP is the slider Poisson ratio. Note that the geometry of the blocks is introduced to allow us to make a link between the elastic constants of the model (K, KT and KL) and the elastic parameters of the sliding medium. It is not an essential feature of the model. In this model, if we push the most-left block of the UL until the most-left frictional contact breaks, a crack emerges and propagates through the interface for a finite distance Λ ≫ λc till it is arrested in qualitative agreement with experiments (Rubinstein et al.2004, 2007, 2011). The UL plays the role of a reservoir, where the elastic energy is stored and partially released at the main shock, while the unreleased part of the elastic energy results in aftershocks satisfying the Omori law (Braun & Tosatti 2014). Figure 2. View largeDownload slide The model. The upper layer (UL) is split in rigid blocks of size λc × λc × NLλc connected by springs of elastic constant KL. The interface layer (IL) is split in rigid blocks of size λc × λc × λc connected by springs of elastic constant K. The UL and IL are coupled by springs of elastic constant KT. The IL is connected with the rigid bottom block (the base) by contacts, represented by ‘frictional’ springs of elastic constant k, which break when the local stress exceeds a threshold value. The UL is driven with the velocity v through springs of elastic constant Kd. Figure 2. View largeDownload slide The model. The upper layer (UL) is split in rigid blocks of size λc × λc × NLλc connected by springs of elastic constant KL. The interface layer (IL) is split in rigid blocks of size λc × λc × λc connected by springs of elastic constant K. The UL and IL are coupled by springs of elastic constant KT. The IL is connected with the rigid bottom block (the base) by contacts, represented by ‘frictional’ springs of elastic constant k, which break when the local stress exceeds a threshold value. The UL is driven with the velocity v through springs of elastic constant Kd. Finally, we attach springs Kd to the top boundary of the UL and drive the ends of these springs with a velocity v to simulate the driving of the seismic fault. For the simulations, we use periodic boundary conditions in the direction of the chain of blocks. 3 SIMULATION RESULTS The simulation results are presented in Figs 3–6. For the chosen set of parameters, we ran 20 independent calculations using different sets of random numbers, and with the protocol described in the Material and method section we obtained 10.8 × 106 shocks above the background level and analysed 7012 main shocks which occurred during the time interval T = 0.67 tmax ≈ 1.34 × 106 after the first 33 per cent of data were discarded in each simulation. Figure 3. View largeDownload slide Time evolution of the system: (a) the frictional force F(t), and (b) the (global) amplitude of EQs $${\cal A} (t)$$ versus time. Figure 3. View largeDownload slide Time evolution of the system: (a) the frictional force F(t), and (b) the (global) amplitude of EQs $${\cal A} (t)$$ versus time. Figure 4. View largeDownload slide Statistics of EQ magnitudes presented in Fig. 3 showing the Gutenberg–Richter power-law behaviour; the dashed line corresponds to the exponent b = 1. Figure 4. View largeDownload slide Statistics of EQ magnitudes presented in Fig. 3 showing the Gutenberg–Richter power-law behaviour; the dashed line corresponds to the exponent b = 1. Figure 5. View largeDownload slide Typical earthquakes in the model. Top panel: the EQ amplitude $${\cal A} (t)$$ versus time for the last 5 per cent of simulation time. Middle panel: detailed $${\cal A} (t)$$ dependences for two time intervals. Bottom panel: the corresponding colour maps of the EQ amplitude on the (t, x) plane, using the colour scale indicated on the right. F—foreshocks, C—calm, M—main shock and A—aftershocks. Figure 5. View largeDownload slide Typical earthquakes in the model. Top panel: the EQ amplitude $${\cal A} (t)$$ versus time for the last 5 per cent of simulation time. Middle panel: detailed $${\cal A} (t)$$ dependences for two time intervals. Bottom panel: the corresponding colour maps of the EQ amplitude on the (t, x) plane, using the colour scale indicated on the right. F—foreshocks, C—calm, M—main shock and A—aftershocks. Figure 6. View largeDownload slide Foreshocks and aftershocks statistics in the presence of contact aging (β = 100): the rate of fore- and aftershocks n(τ) (red solid symbols) and their associated standard deviations (red bars), and the shock amplitudes $${\cal A} (\tau )$$ (blue squares) relative to the corresponding main shock. The horizontal dash lines show the average magnitudes of the foreshocks and aftershocks. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. The results plotted on this figure have been obtained from the analysis of 7012 main shocks, collected over 20 independent simulations. Figure 6. View largeDownload slide Foreshocks and aftershocks statistics in the presence of contact aging (β = 100): the rate of fore- and aftershocks n(τ) (red solid symbols) and their associated standard deviations (red bars), and the shock amplitudes $${\cal A} (\tau )$$ (blue squares) relative to the corresponding main shock. The horizontal dash lines show the average magnitudes of the foreshocks and aftershocks. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. The results plotted on this figure have been obtained from the analysis of 7012 main shocks, collected over 20 independent simulations. The dependences of the driving force F(t) (Fig. 3a) and the amplitude of shocks $${\cal A}(t)$$ (Fig. 3b) on a large timescale both look rather stochastic. As we observed with shock visualization during simulation, and as can also be guessed from comparison of Figs 1 and 3(a), the frictional force is mainly determined by one or few large thresholds with fsi ≫ fs that survived for a long enough time; the breaking of these contacts results in MEQs. The statistics of the shock amplitudes follows the GR power law (see Fig. 2) for an interval of magnitudes $$\Delta {\cal {M}} \gtrsim 2$$, which is however limited by the largest possible magnitude determined by the ratio of the aging rate β and the driving velocity v, $${\cal {M}}_{\rm max} \propto {\rm log}_{10} (\beta /\sqrt{v})$$ (Braun & Peyrard 2013) (as we use here a rather small quasi-1-D system, we cannot expect a very large value of $$\Delta {\cal {M}}$$). Fig. 5 shows the $${\cal A}(t)$$ dependence for the last 5 per cent of the simulation time (top panel) and more detailed pictures at a finer timescale of two typical MEQs (middle panel); the bottom panel of Fig. 5 shows the colour maps of the earthquake amplitude. Fig. 6 presents the main result of the data analysis, the average number of shocks n(τ) in an interval δτ ( δτ = 0.0648), versus τ the rescaled time interval from (τ < 0), or after (τ > 0), the corresponding main shock. We show both the mean value n(τ), computed for the 7012 analysed shocks, and its fluctuations for different MEQs, measured by its standard deviation over all analysed MEQs. The variation of the fluctuations versus τ are the most interesting. Before the MEQ, they are large, and vary widely with τ in the range τ ≲ −0.7, showing a significant and random foreshock activity. Then suddenly they drop to a small value and stay so for the last period before the MEQ, −0.7 ≲ τ < 0. In this range, none of the 7012 analysed MEQs shows any notable foreshock activity. As expected, after the MEQ, Fig. 6 shows a large aftershock activity. Fig. 7 shows that the sharp drop of foreshock activity before large events is clearly due to the aging of the contacts, described by the parameter β in eq. (2). When β is reduced from β = 100 to 3 (β = 0 corresponding to no aging at all), all other parameters being preserved, the calm period disappears and a slight increase of activity immediately before the MEQ is even observed. We can also note that, in the absence of aging, the fluctuations of the foreshock and aftershock activity with time, are much smaller. Similarly the standard deviations of n(τ) between the different main shocks are smaller and almost time-independent. The joint evolution of the aftershocks and foreshocks when β varies shows some similarity with the correlation between the rates of foreshocks and aftershocks found in earthquakes catalogues in different geographic regions (Lippiello et al.2017). Figure 7. View largeDownload slide Same as Fig. 6 with β = 3, that is, with very little aging of the contacts. In this case, the magnitudes of the main events are not as large as in the presence of aging, so that we define the main shocks by $${\cal A}_{\rm MEQ} = \kappa \, {\cal A}_{\rm max}$$, with κ = 0.4. A set of 1205 MEQs was used to draw this figure. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. Figure 7. View largeDownload slide Same as Fig. 6 with β = 3, that is, with very little aging of the contacts. In this case, the magnitudes of the main events are not as large as in the presence of aging, so that we define the main shocks by $${\cal A}_{\rm MEQ} = \kappa \, {\cal A}_{\rm max}$$, with κ = 0.4. A set of 1205 MEQs was used to draw this figure. The inset shows the average number of foreshocks n(τ) using a log–log magnified scale. Relation with the actual scales of earthquakes. To couple our dimensionless units with real ones, we use the following estimation (Helmstetter et al.2004): let us calculate the average time between the MEQs as trec = T/nMEQ ≈ 2.8 × 103. Then, if we identify this time span with, e.g. 100 d (which is the recurrence time of $${\cal M} > {\cal M}_0 = 5$$ earthquakes in Southern California), we obtain the correspondence between time units in our model and in real seismicity:   \begin{equation*} t_0 = 1 \longleftrightarrow \tau _0 \equiv t_0/\alpha \approx 3.8 \times 10^{-3} \longleftrightarrow t_{\rm real} \sim 1 \, {\rm hr}. \end{equation*} Thus, the time span in Fig. 3 corresponds to about 230 yr and the average calm period before the MEQs, τcalm ≲ 0.7 observed in Fig. 6, corresponds to about 7 d. However, this estimation depends on the value of $${\cal M}_0$$ we choose. For example, at the Parkfield segment along the San Andreas fault, California, interplate earthquakes of magnitude about $${\cal M}_0 = 6$$ have occurred at recurrence intervals of 23 ± 9 yr (Kawamura et al.2012), while along the Nankai trough, where the Philippine Sea plate subducts beneath southwestern Japan, great earthquakes of magnitude $${\cal M}_0 = 8$$ have repeatedly occurred every 100 yr (Kawamura et al.2012). If earthquakes in a given region follow the GR law, then the recurrence period depends on $${\cal M}_0$$ as $$t_{\rm rec} \propto \exp ({\cal M}_0)$$. In a simple model, we cannot provide an unquestionable value of the magnitude that could be compared with actual magnitudes, but the main difficulty for a comparison with actual data is that the limited size of the model restricts the range of the observable events. Fig. 6 shows events with a range that spans less than three orders of magnitude, while actual foreshocks can be much smaller with respect to an MEQ. The model cannot show the smallest events and tends to overestimate the duration of the calm period. In actual earthquakes, it may sometimes only last for a few hours or less than a day as in Raleigh et al. (1977), but may extend to much longer periods reaching several months (Main & Meredith 1991). 4 DISCUSSION AND CONCLUSION In order to correctly describe various features of earthquakes, the model that we used in this work is more complex than the simplified case that we presented in the Introduction, and the probability distribution of the breaking thresholds is not imposed a priori but instead results from the dynamics of the model through eq. (1). Fig. 8 shows that the qualitative mechanism for the origin of a calm period that we presented in the introduction is nevertheless also present in this more elaborate case. While the initial probability distribution of the thresholds for contact breaking Pc(fsi) was a Gaussian, the evolution of the thresholds according to eq. (1) leads to a very different shape as time evolves. The maximum value of fsi grows, and moreover the distribution tends to split into a large cluster of contacts which break easily (the ‘weak contacts’ of the qualitative model) and a few very strong contacts. The contacts with intermediate thresholds tend to disappear. It is the few strong contacts that prevent the sliding and are responsible for the MEQs when they finally break. The gap between these strong contacts and the many weak contacts is responsible for the calm period. This evolution is mostly governed by the aging of the contacts. Figure 8. View largeDownload slide Evolution of the probability distribution of the contact thresholds as time grows (from top to bottom figure, as indicated in each frame. Figure 8. View largeDownload slide Evolution of the probability distribution of the contact thresholds as time grows (from top to bottom figure, as indicated in each frame. In natural seismicity quiescence is observed in about 40 per cent of events (Turcotte et al.2000; Kawamura et al.2012): the frequency of small events is gradually enhanced preceding the main shock, whereas, just before the main shock, it is suppressed in a close vicinity of the epicentre of the upcoming event. In about 60 per cent of natural earthquakes another scenario is observed. Large earthquakes are preceded by a period during which the surrounding region experiences a phase of accelerated seismic release (Ben-Zion & Lyakhovsky 2002; McGuire et al.2005). The rate of these foreshocks also follows the Omori law (Kagan & Knopoff 1978; Jones & Molnar 1979), usually with a smaller value of tO. Our model may also demonstrate such a behaviour for smaller values of the ratio β2/v (Braun & Tosatti 2014), when the distribution of large thresholds is more uniform. However, in the model as in reality, earthquakes exhibit a broad dispersion of their features from event to event. This is why seismic quiescence is notable on the standard deviations of the foreshock activity much more than on the mean value of n(τ) averaged over a large number of MEQs (inset of Fig. 6). The intensity of the foreshock activity prior to the calm period varies strongly from event to event and can even stay rather low for some MEQs. This means that using quiescence as a warning signal may be unreliable even for a fault which have the characteristics which lead to quiescence. In nature large earthquakes are almost always followed by a very large seismicity rate (aftershocks). In our model aftershocks are less likely for parameters that lead to a calm period before big events (see Figs 6 and 7). This is consistent with the picture shown in Fig. 8 because the quiescence is associated to an exhaustion of the weak contacts before the MEQ so that the large event is more likely to fully release the stress in the fault. It would be interesting to check whether such an effect is actually observed for natural earthquakes showing quiescence. However, the decrease of aftershock activity may be exaggerated in our calculations due to the restricted size of the model. In an actual earthquake a complete fault line does not break at once, and instead the partial breaking tends to accumulate stress in some regions, as our model does for some parameter sets. This study points out the crucial role of the distribution of thresholds to control the pattern of foreshocks. The existence of a few strong contacts tends to induce a period of seismic quiescence before a large earthquake. Such a distribution could of course come from the geometry of a fault. Some studies correlated different foreshock patterns to general features of faults, such as interplate or intraplate earthquakes (Bouchon et al.2013). However, one important result of our study is that the aging of the contacts can also be a determining factor, in slowly building a multipeaked distribution of contact thresholds. Several mechanisms may be responsible for aging of the sliding interface, such as removing of dislocations from asperities, growth of contact sizes due to plastic deformation, squeezing of a ‘lubricant’ (e.g. water or powdered rock) from the interface or its solidification due to high pressure, coalescence of contacts, etc. These features may be characteristic of a particular fault or a type of rock found in some area, so that statistics of quiescence could differ from place to place. As the timescale of quiescence is short compared to the periods generally involved in foreshock studies, there are very few available data which could support this idea. However, a recent paper (Lippiello et al.2017) presents (in fig. 3) statistics of foreshocks which show some drop of activity before a main shock for the Northern California Earthquake Catalog (Waldhauser & Schaff 2008; particularly for magnitude 2) which does not appear for other catalogues. This could suggest that such local differences exist although this is to be taken with caution because the study did not pay a particular attention to quiescence. A recent study has stressed the dominant role of the formation of interfacial chemical bonds (Qunyang Li et al.2011). It is interesting that the chemistry is also proposed as a cause for seismic quiescence (Main & Meredith 1991). It would be interesting if geologists could correlate the chemical properties of the rocks involved in the faults and the pattern of foreshocks. Real earthquakes, as well as the model that we have presented when its parameters are changed, can exhibit different patterns of foreshocks. An acceleration of the events prior to a large earthquake is not always the rule. Seismic quiescence may have deep implications by showing that one should not take too strictly the viewpoint that earthquakes are associated with criticality. The earthquake may not be the divergence of some accelerating events that culminate in a main shock, but instead appear after some unusually quiet period. Our model of aging shows that such a behaviour may emerge from what could, at first glance, be considered as a secondary effect, the aging of the contacts, which in turn may strongly alter the distribution of the thresholds at which contacts break. Of course one has to be careful in extending conclusions drawn from a simple, 1-D model to the complexity of natural earthquakes. As discussed above, we hope that such a model can nevertheless give useful hints for further studies of real earthquakes. But it is also obvious that developments of the model are clearly needed, in particular, generalizations to 2-D models of the interface and 3-D models of the elastic slider which would improve the GR statistics and provide a power law for earthquake spatial distribution. The process describing the aging of the contacts also certainly needs further investigations as it is presently largely arbitrary. We use a stochastic process that goes to a Gaussian distribution in one limit and a power law in another, but how this process can be generated by the dynamics of ‘macrocontacts’ is still open. Acknowledgements This work was partially supported by CNRS-Ukraine PICS grant no. 151539. OB was supported in part by COST Action MP1303. REFERENCES Ben-Zion Y., Lyakhovsky V., 2002. Accelerated seismic release and related aspects of seismicity patterns on earthquake faults, Pure appl. Geophys . 159, 2385– 2412. Google Scholar CrossRef Search ADS   Bouchon M., Durand V., Marsan D., Karabulut H., Schmittbuhl J., 2013. The long precursory phase of most large interplate earthquakes, Nature Geosci ., 6, 299– 302. Google Scholar CrossRef Search ADS   Bowman D.D., King G.C.P., 2001. Accelerating seismicity and stress accumulation before large earthquakes, Geophys. Res. Lett ., 28, 4039– 4042. Google Scholar CrossRef Search ADS   Braun O.M., 2015. Stick-slip vs. smooth sliding in the multicontact interface, Europhys. Lett ., 109, 48004-1– 48004-6. Google Scholar CrossRef Search ADS   Braun O.M., Peyrard M., 2013 Role of aging in a minimal model of earthquakes, Phys. Rev. E , 87, 032808-1–032808-7. Braun O.M., Scheibert J., 2014. Propagation length of self-healing slip pulses at the onset of sliding: a toy model, Tribol. Lett ., 56, 553– 562. Google Scholar CrossRef Search ADS   Braun O.M., Tosatti E., 2014. Aftershocks in a frictional earthquake model, Phys. Rev. E , 90, 032403-1– 032403-8. Google Scholar CrossRef Search ADS   Braun O.M., Peyrard M., Stryzheus D.V., Tosatti E., 2012. Collective effects at frictional interfaces, Tribol. Lett ., 48, 11– 24. Google Scholar CrossRef Search ADS   Braun O.M., Peyrard M., 2010. Master equation approach to friction at the mesoscale, Phys. Rev. E , 82, 036117-1– 036117-19. Google Scholar CrossRef Search ADS   Burridge R., Knopoff L., 1967. Model and theoretical seismicity, Bull. seism. Soc. Am ., 57, 341– 371. Caroli C., Nozieres Ph., 1998. Hysteresis and elastic interactions of microasperities in dry friction, Eur. Phys. J. B , 4, 233– 246. Google Scholar CrossRef Search ADS   Dieterich J.H., 1972. Time dependent friction as a possible mechanism for afterschocks, J. geophys. Res ., 77, 3771– 3781. Google Scholar CrossRef Search ADS   Dieterich J.H., 1979a. Modeling of rock friction. 1. Experimental results and constitutive equations, J. geophys. Res ., 84, 2161– 2168. Google Scholar CrossRef Search ADS   Dieterich J.H., 1979b. Modeling of rock friction. 2. Simulation of preseismic slip, J. geophys. Res ., 84, 2169– 2175. Google Scholar CrossRef Search ADS   Dieterich J.H., 1994. A constitutive law for the rate of earthquaque production and its application to earthquake clustering, J. geophys. Res ., 90, 2601– 2618. Google Scholar CrossRef Search ADS   Gardiner C.W., 1985. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences , Springer-Verlag, Berlin. Google Scholar CrossRef Search ADS   Grassberger P., 1994. Efficient large-scale simulations of a uniformly driven system, Phys. Rev. E , 49, 2436– 2444. Google Scholar CrossRef Search ADS   Gutenberg B., Richter C.F., 1954. Earthquake magnitude, intensity, energy and acceleration, Bull. seism. Soc. Am ., 46, 105– 146. Hainzl S., Zöller G., Kurths J., 2000a. Similar power laws for foreshock and aftershock sequences in a spring block model for earthquakes, J. geophys. Res ., 104, 7243– 7253. Google Scholar CrossRef Search ADS   Hainzl S., Zöller G., Kurths J., 2000b. Self-organization of spatio-temporal earthquake clusters, Nonlin. Proces. Geophys ., 7, 21– 29. Google Scholar CrossRef Search ADS   Helmstetter A., Hergarten S., Sornette D., 2004. Properties of foreshocks and aftershocks of the nonconservative self-organized critical Olami-Feder-Christensen model, Phys. Rev. E , 70, 046120-1– 046120-13. Google Scholar CrossRef Search ADS   Helmstetter A., Shaw B.E., 2009. Afterslip and aftershocks in the rate-and-state friction law, J. geophys. Res. , 114, B01308-1– B01308-24. Google Scholar CrossRef Search ADS   Hergarten S., Krenn R., 2011. Synchronization and desynchronization in the Olami-Feder-Christensen earthquake model and potential implications for real seismicity, Nonlin. Process. Geophys ., 18, 635– 642. Google Scholar CrossRef Search ADS   Jagla E.A., 2010. Realistic spatial and temporal earthquake distributions in a modified Olami-Feder-Christensen model, Phys. Rev. E , 81, 046117-1– 046117-8. Google Scholar CrossRef Search ADS   Jagla E.A., Kolton A.B., 2010. The mechanisms of spatial and temporal earthquake clustering, J. geophys. Res ., 115, B05312-1– B05312-8. Google Scholar CrossRef Search ADS   Jagla E.A., 2013. Forest-fire analogy to explain the b value of the Gutenberg-Richter law for earthquakes, Phys. Rev. Lett ., 111, 238501-1– 238501-5. Google Scholar CrossRef Search ADS   Jones L.M., Molnar P., 1979. Some characteristics of foreshocks and their possible relationship to earthquake prediction and premonitory slip on faults, J. geophys. Res ., 84, 3596– 3608. Google Scholar CrossRef Search ADS   Jordan T.H., Jones L.M., 2010. Operational earthquake forecasting: Some thoughts on why and how, Seismol. Res. Lett ., 81, 571– 574. Google Scholar CrossRef Search ADS   Kagan Y.Y., Knopoff L., 1978. Statistical study of the occurrence of shallow earthquakes, Geophys. J. R. astr. Soc ., 55, 67– 86. Google Scholar CrossRef Search ADS   Kawamura H., Hatano T., Kato N., Biswas S., Chakrabarti B.K., 2012. Statistical physics of fracture, friction, and earthquakes, Rev. Mod. Phys ., 84, 839– 884. Google Scholar CrossRef Search ADS   Kazemian J., Tiampo K.F., Klein W., Dominguez R., 2015. Foreshock and aftershocks in simple earthquake models, Phys. Rev. Lett ., 114, 088501-1– 088501-6. Google Scholar CrossRef Search ADS   Lippiello E., Giacco F., Marzocchi W., Godano G., de Archangelis L., 2017. Statistical features of foreshocks in instrumental and ETAS catalogs, Pure appl. Geophys ., 174, 1679– 1697. Google Scholar CrossRef Search ADS   Main I.G., Meredith P.G., 1991. Stress corrosion constitutive laws as a possible mechanism for intermediate-term and short-term seismic quiescence, Geophys. J. Int ., 107, 363– 372. Google Scholar CrossRef Search ADS   Marone C., 1988. Laboratory-derived friction laws and their application to seismic faulting, Ann. Rev. Earth. Planet. Sci ., 26, 643– 696. Google Scholar CrossRef Search ADS   Marone C., Scholz C.H., Bilham R., 1991. On the mechanics of earthquake afterslip, J. geophys. Res ., 96 8441– 8452. Google Scholar CrossRef Search ADS   McGuire J.J., Boettcher M.S., Jordan T.H., 2005. Foreshock sequences and short-term earthquake predictability on East Pacific Rise transform faults, Nature , 434, 457– 461. Google Scholar CrossRef Search ADS PubMed  Mega M.S., Allegrini P., Grigolini P., Latora V., Palatella L., Rapisarda A., Vinciguerra S., 2003. Power-law time distribution of large earthquakes, Phys. Rev. Lett ., 90, 188501-1– 188501-4. Google Scholar CrossRef Search ADS   Olami Z., Feder H.J.S., Christensen K., 1992. Self-organized criticality in a continuous, nonconservative cellular automaton modeling earthquakes, Phys. Rev. Lett ., 68, 1244– 1247. Google Scholar CrossRef Search ADS PubMed  Omori F., 1894. On after-shocks of earthquakes, J. Coll. Sci. Imp. Univ. Tokyo  7, 111– 120. Pelletier J.D., 2000. Spring-block models of seismicity: review and analysis of a structurally-heterogeneous model coupled to a viscous asthenosphere, Geophys. Monogr ., 120, 27– 42. Perfettini H., Ampuero J.-P., 2008. Dynamics of a velocity strengthening fault region: implications for slow earthquakes and postseismic slip, J. geophys. Res.  113, B09411-1– B09411-22. Google Scholar CrossRef Search ADS   Qunyang L., Tullis T.E., Goldsby D., Carpick R.W., 2011. Frictional ageing from interfacial bonding and the origins of rate and state friction, Nature , 480, 233– 237. Google Scholar CrossRef Search ADS PubMed  Raleigh C.B. et al.  , 1977. Prediction of the Haicheng earthquake, EOS, Trans. Am. geophys. Un ., 58, 236– 272. Google Scholar CrossRef Search ADS   Rice J.R., Ben-Zion Y., 1996. Slip complexity in earthquake fault models, Proc. Natl. Acad. Sci ., 93, 3811– 3818. Google Scholar CrossRef Search ADS   Rubinstein S.M., Cohen G., Fineberg J., 2004. Detachment fronts and the onset of dynamic friction, Nature , 430, 1005– 1009. Google Scholar CrossRef Search ADS PubMed  Rubinstein S.M., Cohen G., Fineberg J., 2007. Dynamics of precursors to frictional sliding, Phys. Rev. Lett.  98, 226103-1– 226103-4. Google Scholar CrossRef Search ADS   Rubinstein S.M., Barel I., Reches Z., Braun O.M., Urbakh M., Fineberg J., 2011. Slip sequences in laboratory experiments resulting from inhomogeneous shear as analogs of earthquakes associated with a fault edge, Pure appl. Geophys ., 168, 2151– 2166. Google Scholar CrossRef Search ADS   Ruina A., 1983. Slip instability and state variable friction laws, J. geophys. Res ., 88, 10 359–10 370. Rundle J.B., Klein W., 1995. Dynamical segmentation and rupture patterns in a “toy” slider-block model for earthquakes, Nonlin. Process. Geophys ., 2, 61– 79. Google Scholar CrossRef Search ADS   Serino C.A., Tiampo K.F., Klein W., 2011. New approach to Gutenberg-Richter scaling, Phys. Rev. Lett ., 106, 108501-1– 108501-4. Google Scholar CrossRef Search ADS   Socolar J.E.S., Grinstein G., Jayaprakash C., 1993. On self-organized criticality in nonconserving systems, Phys. Rev. E , 47, 2366– 2376. Google Scholar CrossRef Search ADS   Stein R.S., 1999. The role of stress transfer in earthquake occurrence, Nature , 402, 605– 609. Google Scholar CrossRef Search ADS   Turcotte D.L., Shcherbakov R., Rundle J.B., 2000. Complexity and earthquakes, in Treatise on Geophysics , 4 Chap. 23, pp. 675– 700, ed. Kanamori H., Elsevier. Utsu T., Ogata Y., Matsu’ura R.S., 1995. The centenary of the Omori formula for a decay law of aftershock activity, J. Phys. Earth , 43, 1– 33. Google Scholar CrossRef Search ADS   Waldhauser F., Schaff D.P., 2008. Large-scale relocation of two decades of Northern California seismicity using cross-correlation and double-difference methods, J. Geophys. Res.: Solid Earth , 113 (B8) B08311, doi:10.1029/2007JB005479. APPENDIX: METHODS A1 Parameters We use dimensionless units, where four model parameters are set to 1, k = 1, λc = 1, t0 = 1 and fs = 1 so that the other parameters should be compared with those. The IL consists of N = 101 $$\lambda _c^3$$-cubes, each of mass m = 10−4, as chosen in our analysis of stick-slip phenomena (Braun 2015). The springs between the cubes have the elastic constant K = 100. The UL consists of parallelepipeds of height NL = 25 (so that their masses are NLm); for the Poisson ratio, we took σP = 0.3. With these parameters, the crack propagation distance (where aftershocks are expected) is Λ ∼ 30 (Braun & Scheibert 2014). For the rate of aging, we took β = 100, and for the parameter ε we used $$\varepsilon =2 / \Delta f_s^2$$. With these parameters, the statistics of the magnitudes of earthquakes follow the GR law for an interval of magnitudes $$\Delta {\cal M} \sim 2$$, which is of course limited by the size of the model. For the dispersion of thresholds, we took Δfs = 0.3, and used ν = 2 to mimic a delay in contact formation. The initial stretching of the contacts as well as the initial thresholds of the newborn contacts, when they form again after breaking, are taken from a Gaussian distribution with zero mean and dispersion Δfs. The UL is moved through the springs with elastic constant Kd = 0.5 attached to its top and driven with the velocity v = 0.1. The equations of motion are solved by a Runge–Kutta method with a viscous damping η = 0.3 ω0, where ω0 = (K/m)1/2, so that the dynamics of the system is underdamped. With these parameters, the motion of a single $$\lambda _c^3$$-block exhibits a stick-slip behaviour [recall that stick-slip exists only for a certain interval of driving velocities (Braun 2015)]. The duration of the runs is tmax = 2 × 106. The first 33 per cent of a trajectory are discarded to allow the system to lose memory from the initial configuration and reach a steady state (dashed vertical line in Fig. 3). A2 Protocol and data analysis In a simulation, we calculate and store the driving force F(t) and the dimensionless quake amplitude defined as the sum of the force drops due to broken contacts at every time step,   $${\cal A} (t) = \sum _{i=1}^N \Delta F_i (t)/( N f_s ) \,.$$ (A1)The rate of shocks, n(t), is defined as the number of shocks per time unit regardless of their amplitude. The dimensionless earthquake magnitude is defined as   $${\cal M} = \frac{2}{3} \log _{10} \left( {\cal A}/{\cal A}_0 \right) ,$$ (A2)where $${\cal A}_0$$ is a constant which defines the magnitude scale. As for the study of real earthquakes (Bouchon et al.2013), the analysis of foreshocks is delicate because one has to distinguish ‘main shocks’ in the middle of many events, occurring randomly and with a broad range of amplitudes. To be qualified as ‘foreshocks’ events must be related to the main shock both in time and space, and moreover a large foreshock should not be considered as an independent MEQ. In analysing the results, we used the following protocol (Braun & Tosatti 2014). First we remove small ‘background’ earthquakes with amplitudes below some level from the analysis. We only retain $${\cal A} > 2 \, \langle {\cal A} (t) \rangle$$ (see broken red line in Fig. 3). Next, we single out the main shocks above some level $${\cal A}_{\rm MEQ}$$. We took here $${\cal A}_{\rm MEQ} = \kappa \, {\cal A}_{\rm max}$$, with κ = 0.03 (solid red line in Fig. 3). We also used the following rescaling procedure to detect large events which could be related. If nMEQ is the total number of MEQs in a simulation, let us call σ = S/nMEQ, the average area occupied by a single MEQ in the (t, x) plane, with S = Nλctmax. We then rescale the time coordinate t → τ = t/α with α = σ/λ2, where λ ∼ Λ is some distance chosen in such a way that the distribution of MEQs on the (τ, x) plane becomes isotropic (λ ≈ 30 and α ≈ 300 for the parameters used in Fig. 3). In this rescaled space, events which are related to a particular MEQ lie within a circle around this MEQ so that they can be easily identified. We can now scan all MEQ coordinates on the (τ, x) plane and, if the distance ρij = [(τi − τj)2 + (xi − xj)2]1/2 between two MEQs i and j is smaller than some value ρcut then only the larger of these two MEQs remains as the MEQ, while the lower one is removed from the list of MEQs. We chose ρcut = 0.4λ ∼ 0.4Λ, that is, to be considered as related to each other two events must not be separated by more than about 0.4 × the average distance along which a crack propagates. With this protocol, we have obtained a set of well-separated MEQs isotropically occupying the (τ, x) plane, and we may calculate the temporal and spatial distributions of all earthquakes within some area around each MEQ. We consider that all the events separated from the corresponding MEQ by less than ρ0 = ρcut/3 are related to this particular MEQ. They can be either foreshocks or aftershocks. However, we have to make sure that some events that occur before one of the MEQ do not belong to the tail of aftershocks following a violent MEQ which occurred earlier. This is done by examining only MEQs which are separated by a time interval long enough to allow the tail of aftershocks to die out. In practice, only MEQs separated by a τ interval greater than 0.4 × the full τ interval that we scan around MEQs are analysed. This lead us to discard more than 30 per cent of the events that would otherwise qualify as MEQs. Then we collapse all data together, designating τ = 0 for every main shock and normalizing shocks amplitudes on the corresponding main shock value. We compute n(τ) the number of secondary events (foreshocks or aftershocks) in a given δτ = 0.0648 interval averaged over all the investigated MEQs as well as the standard deviations of their values for the different MEQs, which provides an estimate of the error bars for n(τ) and gives us an information on the fluctuations of the earthquake activity in a given τ interval around an MEQ. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial