Universal Earthquake-Occurrence Jumps, Correlations with Time, and Anomalous Diffusion

Spatiotemporal properties of seismicity are investigated for a worldwide (WW) catalog and for Southern California in the stationary case (SC), showing a nearly universal scaling behavior. Distributions of distances between consecutive earthquakes (jumps) are magnitude independent and show two power-law regimes, separated by jump values about 200 km (WW) and 15 km (SC). Distributions of waiting times conditioned to the value of jumps show that both variables are correlated in general, but turn out to be independent when only short or long jumps are considered. Finally, diffusion profiles reflect the shape of the jump distribution.

Earthquake statistics provides important insights about seismic processes [1]. Simple laws of seismicity such as the Gutenberg-Richter frequency-magnitude relation [2,3,4] and the Omori law for the temporal decay of the seismic activity [4] allow for instance the prediction of the size of a characteristic earthquake, the evaluation of the spatial distribution of stress in a fault, or the estimation of seismic risk after a strong event [5,6,7]. It is natural to expect that more sophisticated, multidimensional statistical studies can be very valuable for hazard assessment as well as for understanding fundamental properties of earthquakes. In this way, the structure of earthquake occurrence in magnitude, space, and time should reflect all the complex interactions in the lithosphere.
The first systematic analysis of seismicity taking into account its multidimensional nature was performed by Bak et al., who studied the dependence of waiting-time distributions on the magnitude range and the size of the spatial region selected for analysis, Their study revealed, by means of a scaling law, the complexity of a self-similar dynamical process with a hierarchy of structures over a wide range of scales [8,9,10,11]. Waiting times (also referred to as return times, interevent times, etc.) are just the time intervals between consecutive earthquakes in a region, and can be studied in two ways: i) by using single-region waiting-time distributions, in which an arbitrary region is characterized by its own distribution [12], or ii) the original approach [8], in which a large region is divided into smaller, equally sized areas and waiting times are measured for the smaller areas but are included together into a unique mixed distribution for the whole region. (This procedure can be interpreted as the measurement of first return times of seismic activity to different points in a large, spatially heterogeneous system, see Ref. [9].) The outcomes of both approaches are clearly different, but in any case scaling turns out to be a fundamental tool of analysis, reducing the multidimensional dependence of the waiting-time distributions to simple univariate functions.
Further, it is possible to understand this scaling approach in terms of a renormalization-group transformation [13].
Equally important for risk estimations and forecasting, though much less studied [14,15], should be the statistics of the distances between consecutive earthquakes, which we can identify with jumps (or flights) of earthquake occurrence. Very recently, Davidsen and Paczuski have provided a coherent picture using Bak et al.'s mixed-distributions procedure [16]; in contrast, our paper undertakes the study of the earthquake-distance problem considering the simpler approach of single-region distributions, and for the case of stationary seismicity. The results will lead us to examine the distribution of waiting times conditioned to different values of the jumps; finally, we measure diffusion profiles for seismic occurrence, turning out to be that the profiles are largely determined by the distribution of jumps.
Consider that an arbitrary spatial region and a range of magnitudes have been selected for analysis. The unit vector locating the epicenter of the i−th earthquake on the portion of the Earth surface selected is given byr where the angles ϕ i and θ i denote longitude and latitude, respectively. The spatial distance, or jump, between the i−th event and the immediately previous in time event, i − 1, can be obtained from the angle α i defined by the two vectors, In this way one can measure distances as angles, in degrees; the distance in km is obtained by multiplying α i in radians by the Earth radius (about 6370 km, assuming a spherical Earth).
Given a set of values of the jumps, their probability density D(α) is defined as the probability per unit distance that the distance is in a small interval containing α. In order to avoid boundary problems it is convenient to start the analysis of seismicity on a global scale, which has also the advantage of stationarity (or, more properly, homogeneity in time, at least for the last 30 years).
Stationarity means that any short time period is described by (roughly) the same seismic rate, R (defined as the number of earthquakes per unit time); in such a case, a linear increase of the cumulative number of earthquakes versus time must be observed. Notice that stationarity does not mean that aftershock sequences are not present in the data, rather, many sequences can be interwinned but without a predominant one in the spatial scale of observation.  [16], implying that the spatial occurrence of large earthquakes is no different than the occurrence of small ones, in contrast to claims by some authors [18]. There is, nevertheless, an exception for small distances, less than about 0.1 • (≃ 10 km) for M c from 5.5 to 6, comparable with the size of rupture [19]. The distributions are characterized by a decreasing power-law regime from about 0.1 • to 2 • , with an exponent about 1.6, and a possible increasing second power law for α > 2 • , with exponent 0.3, decaying abruptly close to the maximum distance. Note the difference between these exponents and the value 0.6, obtained in Ref. [16] mixing many areas in Southern California and including non-stationary periods.
However, the measurement of D(α) contains an inconsistency when it is used for distances over a sphere; for instance, after any event, there are many ways to obtain a jump of α = 90 • (360 "ways" in intervals of 1 • ) but there is only one way to get α = 180 • (which is the maximum possible distance, corresponding to the antipodes of the initiating event). Therefore, instead of working with the probability density defined per unit distance, i.e., per unit angle, it is preferable to use (and easier to interpret) the probability density defined per unit solid angle (taking into account all the possible orientations for a given distance α). A straightforward way to estimate this quantity is by means of the transformation due to the fact that 2π sin α dα is the element of solid angle defined by the points at distance α.
The results of the transformation are clearly seen in Fig. 1(a) (top curves). In addition to the power law for the range 0.1 • − 2 • , whose exponent turns out to be ∼ 2.6, the second power law, now decreasing, becomes more apparent from about 2 • to the maximum distance, 180 • , with an exponent ∼ 0.7. This second power law seems to imply a dependence of events at large spatial scales. Indeed, for a given direction, it is more likely a jump with α = 20 • than one of 160 • (independence would imply that the probabilities were the same); however, this effect is due to the fact that earthquake occurrence is not uniform over the Earth but fractal, and the power law reflects the fractal structure of the epicenters. If for long distances there is no causal relation between events, the distribution of earthquake jumps is equivalent to the distribution of distances between any pair of earthquakes (not necessarily consecutive); and as the density corresponding to this distance (per unit distance) is the derivative of the well-known correlation integral, widely used to obtain fractal dimensions, independence at long distances implies where d f is the correlation dimension. We have measured the probability density of the distance between any two earthquakes (per unit solid angle), obtaining a behavior proportional to 1/α 0.85 , for events with M ≥ 6, which implies a correlation dimension d f ≃ 1.15 in reasonable agreement withD(α) ∼ 1/α 0.7 , and confirming the independence of events for α > 2 • .
Turning back to the short-jump regime, the excess of probability given by the corresponding power-law with respect the long-jump power law (associated to the uncorrelated regime) is a clear sign of spatial clustering in earthquake occurrence, extending up to distances of about 2 • ≃ 200 km. It is likely that this clustering extends beyond this limit, but it is not detectable with this procedure as it is hidden below the uncorrelated, long-distance regime. In any case, the power law behavior implies the no existence of a finite correlation length, at least up to 200 km, at variance with the findings of Ref. [20,21]. In consequence we can consider the earthquake jumps as Lévy flights, which are one of the main causes of anomalous diffusion.
In contrast to worldwide seismicity, regional seismicity turns out to be nonstationary, in general, and in the same way as for waiting-time distributions [12], the distributions of earthquake jumps depend on the time window selected for analysis. This problem will be avoided here by considering specific time windows characterized by stationary seismicity; an important realization then is that stationary seismicity is characterized by stationary distributions. Additionally, stationary seismicity has the advantage of a lower magnitude threshold of completeness, due to the absence of large events.
We consider Southern-California seismicity, from the waveform cross-correlation catalog by Shearer et al. [22]. The analysis has been performed on 11 stationary periods, comprised between the years 1986 and 2002, yielding a total time span of 9.25 years and containing 6072 events for M ≥ 2.5. It is very striking, as Fig. 1(b) shows, that the distributions of jumps resemble very much those of the worldwide case, with two power-law regimes, but in a different scale (the separation is at about 0.1 • ). In principle, for a smaller area such as California one would expect that the distribution of jumps is just a truncation of the global one, but Fig. 1(b) clearly refutes this fact, providing a clear illustration of earthquake-occurrence self-similarity in space. In fact, the figure shows these distributions rescaled by a factor L, with L the maximum distance for the region, which is L ≃ 6.5 • for Southern California. and L = 180 • for the worldwide case (also included in the plot). The reasonable data collapse is a signature of the existence of a scaling law for the jump distribution,

D(α) ≃ g(α/L)/L
(forD(α) we need to include an extra factor C of normalization). The scaling law can only be given the status of approximated, as the exponent for short jumps for Southern California seems to be different than in the worldwide case, 2.15 versus 2.6 for the distributionD. Nevertheless, this variation could be due to artifacts in the short-distance properties of the catalogs. On the other hand, the long-jump regime is well described by the worldwide exponent. We recall that the exponents are far from the value of Davidsen and Paczuski for the nonstationary California case [16]. The main difference with the worldwide case is the abrupt decay ofD(α) close to the maximum α. Clearly, this is a boundary effect. We have tried several corrections for this effect (in particular the factor 2π sin α is not right when the ring defined by the set of points at a distance α is not totally included in the region under study) but we have not found a substantial improvement and then we have preferred to keep the things simple.
The distributions of jumps we have obtained, together with previous work on the distribution of waiting times [12], provide a simple picture of earthquake occurrence as a continuous-time random walk. This means that we can understand seismic activity as a random walk over the Earth surface, where one event comes after another at a distance that follows D(α) and after a waiting time τ given by the waiting-time probability density, D(τ). However, an important ingredient is missing in this picture, as we need to take into account the correlations between distances and waiting times.
With this goal in mind, we introduce the conditional waiting-time probability density, D(τ | X ), where | X means that only the cases where X is verified are taken into account; in our case X will be a set of values of the jumps. It turns out that if D(τ | X ) does not depend on X , then τ and X are independent, whereas when D(τ | X ) changes with X , then τ and X are correlated (nonlinearly, in general).
From the behavior of D(α), a natural threshold for the jumps in the worldwide case is α ≃ 2 • . Figure 2(a) compares D(τ | α short ) with D(τ | α long ) for worldwide seismicity, where short and long refer to sets of distances below and above 2 • , respectively. The differences are clear: for short jumps the waiting-time distribution is a decreasing power law for several decades, with an exponent close to 1, and ends in an exponential decay. We can identify these distributions with highly correlated events, i.e., aftershock sequences. On the other hand, for long α the waitingtime distribution seems to be exponential for its full range, compatible with a Poisson process and therefore with independent occurrence. But more surprising than the differences between short and long jumps are perhaps the similarities for short jumps. In fact, there seems to be a radical change of behavior separated by α ≃ 2 • , in the sense that if we are above, or below, this value, the distributions do not change, in other etc., see Fig. 2(a) (for the range between 1 • and 2 • the behavior is not clear as the statistics is low). This means that for short jumps the waiting-time distribution is independent on the value of the jump, and the same happens for long jumps, but when the whole range of jumps is considered this is no longer true and both variables become dependent, in contrast to Ref. [16]. For each set of curves we can fit a gamma distribution, D(τ|α) ∝ e −τ/a /τ 1−γ , turning out to be γ ≃ 0.17 for short α and γ ≃ 0.9 for long α. The latter value is very close to 1, the characteristic value of a Poisson process; in fact, the difference between 1 and 0.9 is not significant within the uncertainty of our data, but if the hypothesis of a Poisson process for long distances could be rejected at any reasonable significance level (for which much more data would be necessary) and a value γ < 1 could be significantly established, this would constitute a support for the existence of long-range earthquake triggering. In any case, our findings are in disagreement with the hypothesis put forward in Ref. [20], which considers seismicity as a process uncorrelated in time but correlated in space.
This result confirms that the universal scaling law found for stationary seismicity in Ref. [12] is in fact a mixture of aftershocks and independent events, turning out to be very striking that this mixture leads to a universal behavior. The reason could be a universal proportion of aftershocks versus mainshocks in stationary seismicity [23]. Further, the conditional time distributions verify a scaling law, where the indices + and − refer to long and short jumps, respectively, and R −1 (R, M c ) is the mean waiting time of the unconditional distribution for the region R and M ≥ M c (obviously, a scaling with the mean of the conditional distribution also holds). Figure 2(b) shows this behavior. The next step is to measure earthquake diffusion directly from data. The fundamental quantity is p(α,t), which we call diffusion profile and gives the probability density that two earthquakes (not consecutive) are at a distance α when they are separated by a time t, i.e., In practice, the single value of t is replaced by a range of values and, as above, we introducẽ p(α,t) ≡ p(α,t)/(2π sin α). In the case of normal diffusionp(α,t) is given by a Gaussian distribution (more precisely, a semi-Gaussian, as α > 0), with a second moment scaling as α 2 ∼ t, whereas p(α,t) is given by a Rayleigh distribution (or a Maxwell distribution in three dimensions).
In contrast, our measurements yield to results far from normality, see Fig. 3. For short times, the profile resembles the distribution of jumps with 2 regimes; as time evolves events migrate farther from the origin, towards the long-range part of the curve, which becomes dominant for all α for long times, with an exponent of value 0.8 forp.
A scaling law also holds for the diffusion profile, i.e.,

p(α,t) = h(α/L,t)/L.
Finally, the process turns out to be clearly subdiffusive, as α 2 grows more slowly than linearly, although we cannot provide a definite value of a exponent as in Refs. [24,25].
This paper would not have been possible without the effort of the researchers who have made earthquake catalogs publicly available. The Ramón y Cajal program and the projects BFM2003-06033 (MCyT) and 2001-SGR-00186 (DGRGC) are acknowledged.