Correlated earthquakes in a self-organized model

Motivated by the fact that empirical time series of earthquakes exhibit long-range correlations in space and time and the Gutenberg-Richter distribution of magnitudes, we propose a simple fault model that can account for these types of scale-invariance. It is an avalanching process that displays power-laws in the event sizes, in the epicenter distances as well as in the waiting-time distributions, and also aftershock rates obeying a generalized Omori law. We thus confirm that there is a relation between temporal and spatial clustering of the activity in this kind of models. The fluctuating boundaries of possible slipping areas show that the size of the largest possible earthquake is not always maximal, and the average correlation length is a fraction of the system size. This suggests that there is a concrete alternative to the extreme interpretation of self-organized criticality as a process in which every small event can cascade to an arbitrary large one: the new picture includes fluctuating domains of coherent stress field as part of the global self-organization. Moreover, this picture can be more easily compared with other scenarios discussing fluctuating correlations lengths in seismicity.


Introduction
At the moment there is not a comprehensive explanation of the mechanisms giving rise to the complex phenomenology of earthquakes.
The magnitude of each earthquake is characterized by the Gutenberg-Richter (GR) law (Gutenberg and Richter, 1944), which is in fact a scaleinvariant distribution of energy release. Earthquakes are also long-range correlated with each other. It is indeed known that events are clustered in space and time (Turcotte, 1997;Scholz, 2002) and take place in complex fault patterns (Bonnet et al., 2001). The Omori law of aftershocks rate (Utsu et al., 1995) is an example of the temporal clustering of earthquakes, with a decay given by a scale-invariant law. The phenomenology of the distance between subsequent epicenters is also characterized by power-law distributions (Davidsen and Paczuski, 2005;Corral, 2006). Moreover, the values of magnitudes, waiting times and locations of earthquakes are part of a single scaling picture (Bak et al., 2002;Corral, 2003Corral, , 2004Corral, , 2005. Other examples are given by Mega et al. (2003) and Davidsen et al. (2006). Since seismicity is one of the most outstanding examples of a class of phenomena involving a wide range of energetic, spatial, and and temporal scales, it is expected that its modeling is problematic.
It is possible to build models based upon the phenomenology of earthquakes. For example, aftershock-sequence models require an assumed law of off-spring generation per event (Ogata, 1988;Helmstetter and Sornette, 2002;Turcotte et al., 2007;Lippiello et al., 2007). These models can yield realistic time-series, but by construction they use rather than explain laws like the GR one.
The scale-invariant distribution of earthquake sizes is reproduced by processes based on avalanches of stress redistribution, following the idea that there is self-organized criticality (SOC) (Bak, 1996;Sornette, 2000). The precursor of this concept in geophysics has been the slider-block model by Burridge and Knopoff (1967). It is evident from many models that the mechanism of avalanches of relaxations robustly leads to size-frequency power-laws. This behavior emerges from the collective organization of units that cooperate with very nonlinear rules, redistributing stress and typically dissipating it from open boundaries.
However, it has become also clear during the last years that the simplest SOC models cannot reproduce other important features of critical phenomena, usually involving correlations between events.
Models incorporating correlated events (Olami et al., 1992;Hainzl et al., 1999Hainzl et al., , 2000Hergarten and Neugebauer, 2002;Zöller et al., 2005;Huang et al., 1998;Lippiello et al., 2005;Baiesi and Maes, 2006;Lippiello et al., 2006;Abaimov et al., 2007) are a mi-nority within the literature on SOC. These few scattered results unfortunately have not constituted a large enough body for appropriately raising the issue of temporal organization to the attention of the scientific community.
In this paper we show that earthquakes phenomenology can guide us to build self-organized models with the appropriate features. In particular, we stress the importance of clustering events in space and time, an aspect leading us to develop a fault model that displays a full spectrum of powerlaw statistics (GR law, Omori law, waiting times and epicenter distances with broad distributions), not observed in previous models. Hence, the very basic idea of SOC is in fact achievable. In particular, the process self-organizes the epicenter locations, clustering them rather than spreading them randomly in space, as it is frequently imposed in other simple models.
A novel feature distinguishing the model we propose from previous ones is is the possibility to infer maximal areas of events from its configuration. It turns out that this model does not conform to the common picture associated with SOC in geophysics (Nature debate, 1999;Geller et al., 1997). The idea is that every tremor can in principle cascade in a large event, depending on minor details of the stress field. It is possible that the paradigm of sandpiles has been much influential in the consolidation of this view. Up to date, this interpretation has been a speculation, without any quantitative assessment of its validity. Below we show that we instead observe a mean correlation length limited to a given fraction of the whole fault, and a rich dynamical regime leading to complex patterns of possible slipping areas. The domains where avalanches can occur are not always maximal. Therefore, it is clear that in this model it is not possible to have a large earthquake at all times. We will come back to this point in the Section "Discussion". The next section contains the description of the model, while the numerical results are shown in section 3.

Model
The following model describes a one-dimensional fault with L units and with periodic boundary conditions. Each unit i represents the displacement h i of a plate with respect to a second one. Plates are sliding with respect to each other and thus the displacement h i corresponds to a slip accumulated with time. An external field σ i characterizes the speed of the strain accumulation in the unit: At each time step a unit i, chosen with probability p i ∼ exp(βσ i ), slips: If h i forms a high gradient with one of its neighbors j, in our case h i − h j ≥ 4, a local elastic instability occurs. This is relaxed by allowing the two nearest-neighbor units to get closer, If this process leads to the formation of new unstable couples (i, j), they are listed and processed into a random order until the list is empty, filling at the same time another list with eventual new unstable pairs. The new list is then processed, and so on. The iteration of this rule leads to a final state in which all bonds between units are stable again. The whole avalanche of relaxations represents an earthquake and is characterized by its size (the number of single relaxations, corresponding to the seismic moment), by its slipping area (the number of sites involved at least once), and by its epicenter (the unit where the avalanche started). It takes place by definition in one time step. The waiting time between avalanches is then measured by the number of time steps separating them. The aim of the field σ i is to reproduce some "external" tectonic loading, which should be originated by the crust portions that meet at the fault. Somewhat σ replaces the loading calculated explicitly with the laws of elasticity in other models (see for example (Ben-Zion, 1996;Ben-Zion et al., 2003)). Since earthquakes play the main role in reshaping the stress field in the crust, we let each σ i evolve with a rule that couples it with the activity in the system: Every time that a redistribution (2) occurs, the two corresponding fields are set equal to their average σ ij = (σ i +σ j )/2 plus a noise term δ drawn at random (for each site) from the interval [−1, 1] 1 : The evolution of the system is thus stochastic in many aspects. At the level of single redistributions involving (2) and (3), one has an update of σ's with random δ's. At the step (1) of forcing the system, the choice of i according to a probability p i is also stochastic. One can interpret the set of σ i as an array of local rates. Indeed, a micro-slip (h i → h i + 1) takes place with a rate proportional to exp(βσ i ).
A non-trivial regime emerges as long as β is sufficiently large to lead to a persistence of the earthquake activity in areas of the system. For β → ∞ one finds a choice of the position to apply (1) that corresponds to the site with the largest σ. This resembles an extremal dynamics for the field σ i . We rather chose β large but finite, such that many parts of the fault are likely to be active at the same time (if they are share similar values of σ i ). The evolution of the σ i guarantee a migration of active areas as well.
Despite the stochastic character of some of the microscopic updates, a rich phenomenology arises, with scale-free avalanches and with realistic interoccurrence statistics.

Results
We show results obtained by fixing β = 4, which is large enough to lead to clustering of epicenters. A preliminary check has shown qualitatively similar results in the range 2 ≤ β ≤ 6. For each L, initial configurations for simplicity have h i = 0 and σ i = 0. To be confident that the stationary regime has been reached, we first run a long transient of ≈ 10 8 ÷10 9 time steps without collecting statistics. From time step t = 0 we then collect time series composed by 2 ÷ 3 × 10 8 time steps. This constitutes a satisfactory statistics only if a large number of different profiles is sampled, which is the case for systems with 2000 units. We can thus collect data in a reasonable time for systems up to this size.
A first glance at the behavior of the model is proposed in Fig. 1, where we plot a sample of size and location of rupture areas as a function of time. One can see that the activity is an alternation of earthquakes of several sizes, with a persistence in active areas. This is confirmed by a plot of the increment of h i with respect to the values at time t = 0: Fig. 2(b) shows that the increments are concentrated in the active areas.
The statistics of several quantities turn out to be determined by power-laws. In order to display the frequency-size statistics, we adopt the following definition of magnitude: m = log 10 s Note that the usual prefactor 2/3 (Scholz, 2002) in the conversion from seismic moment to magnitude is not suitable for a one-dimensional model because the area of events is in fact a length. In Fig. 3 one can see that the number of events with magnitude ≥ m, denoted by N > (m), seems to follow a GR law, N > (m) ∼ 10 −bm , with b = 1.1 ± 0.1, though this distribution is most likely multiscaling, as it is often the case in one-dimensional automata (Kadanoff et al., 1989). We postpone the exact characterization of this distribution to future work. The distribution of slipping areas a instead has a clearer scaling: it develops a power-law tail ∼ a −τa for increasing L, with τ a = 1.5 (Fig 4), and obeys to standard finite-size scaling with D = 1 and where F is a scaling function, see inset of In addition to the avalanche size and area, in this model we can also measure metric properties characterizing the state of the system between two avalanches: one is the length of domains of units having constant sign in the slope of h i . Each profile h i is indeed an alternation of domains with increasing h and domains of decreasing h, forming in general a nontrivial landscape, see Fig. 2(a). This is also a result of the self-organization of the process, which includes the evolution of the σ i . Also domain lengths ℓ have a power-law distribution ∼ ℓ −τ ℓ with τ ℓ ≃ 1.9, see Fig. 5, which displays finite-size scaling also with D = 1 (inset of Fig. 5). Since τ a < τ ℓ , there is more chance to observe large areas than large domains. On the other hand, avalanches take place within domains. This suggests that avalanches are repetitive and appear more frequently in long domains.
Connected with the scale-invariance of domains, there is also a scaling of the correlation length of the stress f i = h i+1 − h i with the system size. The correlation length can be read from the shape of the correlation function where . . . means a statistical average over the sites and con- figurations 2 . It turns out that C L (r) conforms to a scaling function C L (r) ≃ C(r/L), with C(. . .) independent on L, as shown in Fig. 6. Hence, if we define the correlation length as the range where C L (r) > 0.1, we see (Fig. 6) that it has a value ≈ 10%L that diverges linearly with L, as one expects in critical systems. We will come back to this point is the Discussion.
Another quantity of interest is the jump between the position of grain addition at time t and the subsequent position of grain addition at t + 1. The jump distributions have also power-law tails, with exponent converging to ≈ −1, see Fig. 7. This distribution is thus similar to that of distances between subsequent earthquakes (Davidsen and Paczuski, 2005;Corral, 2006). Also the crossover to a background level for long jumps takes place at a length that is a fixed fraction the size of the catalogue (Davidsen and Paczuski, 2005;Corral, 2006 Fig. 8. Distribution of waiting times, for events larger than thresholds s (L = 2048). Two power-law fits are also shown for the two parts of the distribution relative to s = 30000.

Temporal correlations
During the last years part of the scientific debate on earthquake correlations has been focusing on the statistics of waiting times between events, see (Baiesi and Maes, 2006) for an overview. An issue was whether SOC models can have avalanches correlated with each other. Some models have waiting times between avalanches with an exponential distribution, suggesting that their events are completely uncorrelated. Clearly this is an unwanted feature in models of earthquakes. Recently Bak et al. (2002) and Corral (2003Corral ( , 2004Corral ( , 2005 have shown that waiting times have in general a nontrivial scaling form in their distributions. In Fig. 8 we plot some waiting time distributions that we observe in our model, for L = 2048 and for several minimum thresholds s of the size. These distributions have a shape with a double power-law form for high thresholds, as observed in catalogs of regional seismicity by Corral (2003) and in an aftershock-sequence model by Lippiello et al. (2007).
In Fig. 9 there is an attempt to collapse some of these distributions on a single curve, by rescaling the waiting times to scales in which their average value is 1, that is, by multiplying their values by the rate of events larger than the corresponding minimum thresholds. This procedure revealed an interesting scaling form for real earthquakes (Corral, 2003(Corral, , 2004(Corral, , 2005 (and also for solar flares, see ): in that case one observes a nice data collapse, with distributions being described by a single scaling function. The data collapse for this model is only approximate. We can conclude that the power-law tails in the distributions are a clear indication of a non-trivial organization and clustering in time of the avalanches, with some missing scale-invariance evidenced by the thresholding procedure.
It is also not trivial to observe aftershocks in simple models of seismicity. Indeed, one does not always observe Omori decay of aftershocks in synthetic catalogs. However, this L=2048, s=300 L=2048, s=1000 L=2048, s=3000 L=2048, s=10000 L=2048, s=30000 is a salient feature of seismicity, characterizing the occurrence of correlated events even for years (Utsu et al., 1995;Shcherbakov et al., 2004;Paczuski, 2004, 2005;Zaliapin et al., 2008). Our model does not yield time series with patterns clearly identifiable with aftershocks sequences, intended in the usual seismological sense. Nevertheless, an Omori-like decay can be detected, confirming the temporal clustering evidenced by waiting time statistics. To visualize the Omori decay, we use a simple definition of aftershocks, leaving more complicated spatio-temporal analysis (Shcherbakov et al., 2004;Paczuski, 2004, 2005;Baiesi, 2006;Zaliapin et al., 2008) for future works. Let us consider events with size s M as main shocks (to improve the statistics, we actually consider events in a range [0.9s M , 1.1s M ]). Each of these events collects aftershocks in a time-window following its occurrence time t M and including only events of smaller size. This time window t−t M thus ends if a new event of size at least 0.9s M occurs. The averaged statistics of the rate r(t − t M ) of avalanches after an main event of size s M is shown in Fig. 10 as a function of the time lag t − t M from the main shock, for several values of s M .
One can see that the aftershock decays depend on s M and follow a generalized Omori decay where A is a constant, t * is a characteristic time, and p is the exponent of the generalized decay (usually one observes p ≈ 1). As in real seismicity Paczuski, 2004, 2005), the onset of the power-law decay takes place at times t * that increase with the size of the main event. The same is true for the end of the Omori decay: data in Fig. 10 have an exponential decay after the Omori regime, as it was found for aftershocks (Baiesi and Paczuski, 2005). The exponent p takes values ranging from ≈ 1.3 for s M = 300, to ≈ 0.5 for s M = 10000. Its variability somewhat reflects the same lack of invariance for increasing thresholds manifested by waiting-time distributions.

Discussion
Some previous SOC models with realistic phenomenology are based on the mechanism of extremal dynamics (Olami et al., 1992;Hainzl et al., 1999Hainzl et al., , 2000Hergarten and Neugebauer, 2002;Zöller et al., 2005;Lippiello et al., 2005), in which an earthquake starts always from the weakest unit. Our stochastic model shows a more general mechanisms giving rise to correlated events within SOC, which involves activity suitably clustered in space and time, together with scale-free redistributions of energy in the form of avalanches. The random aspect cannot be excessive 3 : a load completely random in space has been for years the standard in several SOC cellular automata, maybe because it is the simplest protocol. In the field of seismicity this choice is not supported by phenomenological observations, as we know that epicenters are correlated and clustered. When a random load was imposed, avalanches were found to be uncorrelated (Baiesi and Maes, 2006). We thus argue that a (correct) clustering in space of events cannot be disentangled from the temporal clustering of events, both aspects being part of the same global organization in critical systems.
Regardless of the lack of dissipation from open boundaries, our process reaches a stationary critical regime. The reason is that its loading is not homogeneous and the evolution via avalanches generates the domains over which further large avalanches can occur. In the periodic system we have described, the minima of the accumulated slip profile are places where eventually avalanches must stop. These minima are not fixed but dynamic.
It is important to note that the dynamics of the accumulated slip profile, with domains that evolve in time, has nontrivial consequences. Each domain seems to represent what is normally observed in canonical SOC systems with open boundaries (Bak, 1996), the so called "sandpiles", which have a profile with a single slope, from the maximum at a closed boundary to a minimum at an open (dissipating) boundary. Eventually the whole process somewhat resembles a collection of smaller homogeneous SOC systems, whose number and position fluctuates in time. For each configuration, the maximum correlation length should be close to the length of the longest domain. Interestingly, this domain length is not always close to its possible maximum, which means that the system is often in a state incompatible with an earthquake spanning the whole fault. Moreover, we have seen that the range of the average correlation length is a fraction of the system size. On the one side, this says that we have to reconsider the typical value of correlation ranges upon change of scale of the whole system. Provided that we can meaningfully isolate an area from the rest of the crust, on the other hand, we can expect a finite mean correlation length within it. Hence, our model does not reproduce a popular picture associated with SOC, invoking a continuous state of "maximal" criticality in the crust due to an eventual infinite correlation length (Nature debate, 1999). According to this picture, earthquakes are inherently unpredictable in size, space and time because their cascade to large events depends on minor details of the stress field. This point has been used, for example, by Geller et al. (1997) to infer that earthquakes cannot be predicted. The validity of their argument can be limited by the lack of discussion about non-minor details. These major details in our models are those that are macroscopically visible when looking at the profile of the slip field h i , namely the different domains. Unfortunately patterns like these are not accessible in real measurements. Bak pointed out (Nature debate, 1999) that an earthquake does not "know how large it will become". This is not incompatible with our point that an earthquake "knows how large it cannot become". Perhaps both aspects should be taken into account in studies on earthquake prediction .
Therefore, according to our results, the following scenario is possible: The process of self-organization in seismicity, due to the slow load of the crust and its fast relaxation via earthquakes, converges to a dynamical SOC regime, with rise and fall of patterns of strongly correlated stress. These patterns may be associated with (local) fluctuating correlation lengths.
One could also have coexistence of SOC and other mechanisms (Sammis and . A previous SOC model with a heterogeneous fixed pattern of faults (Huang et al., 7 1998) has a behavior consistent with the hypothesis that the approach to large earthquakes is described by a criticalpoint picture (Jaumé and Sykes, 1999;Sammis and Sornette, 2002), with a finite-time singularity of Benioff strain release and a divergence of a correlation length (Zöller and Hainzl, 2002;Zaliapin et al., 2002). We have not investigated this point in our model yet, though it seems that its dynamics does not break all the correlations after a large earthquake. Indeed, a large slip along a domain lowers the total energy stored in the system, and eventually shifts the domain range of some units, but the domain itself should be ready for similar earthquakes without too much effort. However, an eventual merging with other coherent domains might lead to an increase of the correlation length in the area, with a possible connection with previous studies (Jaumé and Sykes, 1999;Sammis and Sornette, 2002;Zöller and Hainzl, 2002;Zaliapin et al., 2002). In any case, the stationary regime of our model appears to be different from that of intermittent criticality (Ben-Zion et al., 2003;Bowman and Sammis, 2004), in which every large event drives the system far from criticality, which is then slowly restored by the dynamics.

Conclusions
We have shown that it is possible to build stochastic processes with self-organized criticality that reproduce several power-laws found in earthquake statistics, like the GR law, the generalized Omori laws, the waiting-time distributions, and the distributions of distances between subsequent events. The robust scale-invariant statistics generated by avalanches is the leading principle of the study. We have stressed that it is important that the process generates activity clustered in space for eventually obtaining a clustering in time of events. In our model, contrary to previous examples, the clustering takes place even if the system is stochastic, showing that a moderate degree of randomness can be tolerated in SOC models with spatio-temporal correlation.
Our findings are contrary to a constant complete unpredictability of event sizes, even if SOC is one of the main mechanisms acting to generate the complexity of seismicity. The point is that every SOC system can have a finite size. We have described a system displaying domains that resemble a fluctuating collection of canonical SOC cellular automata. Interestingly, domains limit each other and their boundaries constitute the points where avalanches eventually must stop. The size of each domain quantifies locally the correlation length, which is thus a quantity fluctuating in space and time. As a result, the mean correlation length diverges with the system size, as expected, but it occupies only a finite fraction of the system. We have been able to visualize these features thanks to the simplicity of the one-dimensional model.
The oversimplified process that we have discussed is inspired by the phenomenology of earthquakes and tries to en-capsulate it, but clearly it is a geophysical model in an embryonic stage. Hopefully the results and discussion we have presented provide new ideas that will be useful for building models grounded on laws of geophysics and elasticity of solids, which still preserve the ability to reproduce earthquakes phenomenology. With models of this kind, for example, it would be interesting to see if creeping sections of faults can play the role of domain boundaries in the sense discussed in this paper.