# Bulletin of the Seismological Society of America

- Copyright © 2006 Bulletin of the Seismological Society of America

## Abstract

The statistics of time delays between successive earthquakes has recently been claimed to be universal and to show the existence of clustering beyond the duration of aftershock bursts. We demonstrate that these claims are unjustified. Stochastic simulations with Poissonian background activity and triggered Omori- type aftershock sequences are shown to reproduce the interevent-time distributions observed on different spatial and magnitude scales in California. Thus the empirical distribution can be explained without any additional long-term clustering. Furthermore, we find that the shape of the interevent-time distribution, which can be approximated by the gamma distribution, is determined by the percentage of mainshocks in the catalog. This percentage can be calculated by the mean and variance of the interevent times and varies between 5% and 90% for different regions in California. Our investigation of stochastic simulations indicates that the interevent- time distribution provides a nonparametric reconstruction of the mainshock magnitude-frequency distribution that is superior to standard declustering algorithm.

## Introduction

In the past, earthquake statistics were successful in revealing some stable characteristics of seismicity such as the Gutenberg–Richter relation for the magnitude-frequency distribution (Gutenberg and Richter, 1956) and the Omori law for the decay of aftershock activity (Utsu *et al.*, 1995). Recently, the statistics of the waiting times between consecutive earthquakes (so-called interevent times) have also increasingly become the focus of research (Bak *et al.*, 2002; Corral, 2004; Davidsen and Goltz, 2004). Analyzing a number of different seismic catalogs, Corral (2004) found that a unique probability density function can describe the observed interevent times, namely the gamma distribution 1 with constants *C* = 0.5 ± 0.1, *γ* = 0.67 ± 0.05, and *β* = 1.58 ± 0.15. Here, *τ* is the normalized interevent time that is obtained by multiplying the interevent time Δ*t* with the earthquake rate *λ*, that is, *τ* = *λ*Δ*t*. This distribution was claimed to be universal for stationary seismicity; that is, it should hold from worldwide to local scales and for all magnitude ranges. However, the condition of stationarity implies a selection of regions where aftershock activity is not dominant.

Based on some general assumptions, Molchan (2005) could theoretically show that, in agreement with equation (1), the distribution decays exponentially for large interevent times and that the value 1/*β* is the fraction of mainshocks among all seismic events. His only assumption is that the seismicity consists of a Poissonian background activity and triggered aftershocks that are supposed to follow the Omori law. Thus, equation (1) can only be universal if the fraction of mainshocks within the activity is constant and close to 60% (for *β* = 1.58). However, the existing estimations of the mainshock fraction vary between 10% and over 90% for different regions and seem to be inconsistent with a single, universal number (e.g., Reasenberg, 1985; Kagan, 1991). If the fraction of mainshocks is not universal, then the constant 1/*β* in equation (1) should vary for different regions. In this case, analyzing the interevent-time distribution could yield a nonparametric estimate of the mainshock rate. In time- independent seismic-hazard studies, the determination of mainshock activity usually involves applying a declustering procedure such as the Reasenberg algorithm (Reasenberg, 1985), which requires some more or less arbitrary definitions of space and time windows. Thus, a nonparametric estimate could be of great importance.

In this study, we systematically investigate the correlation of the interevent-time distribution and the level of background activity in the western United States and in stochastic simulations where temporally uncorrelated background events trigger Omori-type aftershock sequences. We show in the Analysis of Interevent-Time Distribution Section that the simulations are in agreement with the empirical seismicity from the viewpoint of interevent-time statistics and that the shape of the interevent-time distribution is strongly correlated with the fraction of background seismicity in the catalog. In the Reconstruction of the Mainshock Magnitude-Frequency Distribution Section, we demonstrate that the interevent-time distribution can be used to reconstruct the magnitude-frequency distribution of the mainshock activity.

## Analyzed Data

### Earthquakes in California/Nevada

We analyze the Advanced National Seismic System (anss) catalog of earthquakes having occurred in the time interval between 1 January 1980 and 1 January 2004 within 32.5° and 43° latitude and −113° and −123° longitude (Advanced National Seismic System, 2005). These particular spatial and temporal intervals are chosen in order to guarantee a high-quality and homogeneous data set. The frequency-magnitude distribution of this data set follows the Gutenberg–Richter relation for magnitudes above 2, indicating that the catalog is complete for *M* >2.

The whole region is divided into squares of size *L* × *L*, where all earthquakes with magnitudes *M* ≥*M*_{min} are taken into account. The standard choice is *L* = 100 km and *M*_{min} = 3; however, we additionally analyze all combinations of *L* = 50, 100, 200, and 500 km, and *M*_{min} = 2.5, 3, and 4.

### Earthquake Simulations

Additionally, we analyze simulations of the epidemic- type aftershock sequences (etas) model, which is a stochastic point process, where each earthquake has some magnitude-dependent ability to trigger its own Omori-law- type aftershocks (Ogata, 1988; Helmstetter and Sornette, 2002). In particular, the rate of aftershocks induced by an earthquake that occurred at time *t _{i}* with magnitude

*M*is given by 2 for time

_{i}*t*>

*t*. The parameters

_{i}*K*,

*α*,

*c*, and

*p*are constant for all earthquakes of a given area. The total occurrence rate is the sum of the rate of all preceding earthquakes and a constant background rate,

*λ*(

*t*) =

*λ*

_{0}+ . The background rate

*λ*

_{0}is usually assumed to result from stress accumulation at tectonic-plate boundaries due to tectonic- plate motion.

We produce Monte Carlo simulations of the etas model using the inverse transform method proposed by Felzer *et al.* (2002). Earthquake magnitudes are drawn from a Gutenberg–Richter magnitude-frequency distribution. The background events that occur with a constant rate *λ*_{0} are uniformly located in a 100 × 100 km square. The timing of the aftershocks is calculated from a nonstationary Poissonian function based on the modified Omori law (equation 2). The aftershock locations are chosen according to an isotropic probability density distribution 3 where *q* = 1.69 is chosen according to the inverted value given by Zhuang *et al.* (2004) and *R _{m}* = 0.011 · 10

^{0.5αMm}(km) is the assumed extension of a mainshock rupture with magnitude

*M*(Reasenberg, 1985).

_{m}We performed 1000 simulations of 50 years' duration, where the minimum magnitude is chosen to be *M*_{min} = 3 and the maximum magnitude *M*_{max} = 7. To account for variable conditions, for each simulation, the parameters of the magnitude-frequency distribution are randomly selected from the uniform distributions *a* ∈ [3.0, 5.0] and *b* ∈ [0.8, 1.2]. The etas parameters are chosen from the intervals *α* ∈ [0.7, 1.0], *c* ∈ [1*min*, 1*h*], *p* ∈ [1.05, 1.2], and *n* ∈ [0.4, 0.95], where *n* is the branching parameter of the etas model, *n* = *K* · *f*(*α, b, c, p, M*_{min}, *M*_{max}) = *Kbc*^{1−p}(1 − 10^{(α−b)(Mmax−Mmin)})/[(*b − α*)(*p* − 1)(1 − 10^{−b(Mmax−Mmin)})] (Helmstetter and Sornette, 2002). The branching parameter is the average number of aftershocks per mainshock and defines the fraction of aftershocks in infinite long simulations. Setting these parameters, the parameter *K* is fixed as *K* = *n*/*f*(*α, b, c, p, M*_{min}, *M*_{max}). We fix *n* instead of *K* to get reasonable levels of aftershock activity and to avoid supercritical branching parameters (*n* ≥ 1).

In addition to the etas model simulations, we analyzed model simulations in which only background events can trigger aftershocks according to equation (2) and where the aftershock magnitudes cannot exceed the mainshock magnitude. We will call this model the simple-type aftershock sequence (stas) model. We also performed 1000 stas model simulations with the same parameters as for the etas model.

## Analysis of Interevent-Time Distribution

The interevent times Δ*t _{i}* =

*t*

_{i}− t_{i−1}are calculated from the occurrence times

*t*of the earthquakes and normalized by multiplying with the total earthquake rate

_{i}*λ*=

*N*/

*T*(

*T*= catalog time span), that is,

*τ*=

_{i}*λ*Δ

*t*. We determine the interevent-time distribution in all cases where the number

_{i}*N*of earthquakes with

*M*≥

*M*

_{min}is larger than 50. Figure 1a shows the stacked probability distributions of normalized intervent times for the cells in California with

*L*= 100 km and

*M*

_{min}= 3. The analogous plot for the etas simulations is shown in Figure 1b. The result for both observed and simulated data is almost identical. The distributions clearly deviate from an exponential distribution that would be expected in the case of a Poissonian process, that is, in the case of temporarily uncorrelated activity. For small interevent times, the distribution decays according to

*τ*

^{−1}. This can be explained by Omori-like aftershock activity with a

*p*-value close to 1. For a single aftershock sequence according to the Omori law

*t*

^{−p}the interevent times are known to be distributed according to a power law

*τ*

^{−p̃}, where the relation between

*p*and

*p̃*is given by

*p̃*= 2 −

*p*

^{−1}(Utsu

*et al.*, 1995).

Now we fit the data by the normalized gamma distribution. In this case, the constant *C* in equation (1) is given by *C* = (*β ^{γ}*Γ(

*γ*))

^{−1}, where Γ(

*x*) is the gamma function. The parameter

*β*is determined simply by calculating the mean and the variance of the interevent times, namely,

*β*= and

*γ*= /

*β*(e.g., Denker and Woyczynski, 1998). Due to the normalization of the interevent times, the mean is in our case simply = 1 and we have and

*γ*= 1/

*β*. The interevent-time distributions can be fitted quite well by the gamma distribution. The result is shown for two examples from California in Figure 2a, where 1/

*β*is determined to be 0.11 and 0.47. Following Molchan (2005), this corresponds to relative background activity of 11% and 47%. In Figure 2b, the same data sets are compared with the resulting interevent-time distributions for long simulations of the etas model (200,000 events) with background activity of 11% and 47%. The agreement between the simulated and observed data is almost perfect and shows once more that the etas model captures the mechanisms responsible for the shape of the interevent-time distribution. Furthermore, it is a first indication that 1/

*β*really defines the mainshock fraction.

To investigate the correlation between the shape of the interevent-time distribution and the background level *μ* in a systematic way, we calculate *μ*_{inter} = 1/*β* = for all cells in California and compare these values with the estimations resulting from a standard declustering procedure. For this purpose, we apply the declustering algorithm of Reasenberg (1985) to the anss earthquake catalog, with the declustering parameters set to standard values for California (*Q* = 10, *P* = 0.99, *τ*_{0} = 1 day, *τ*_{max} = 10 days). The algorithm identifies 55% of the earthquakes as aftershocks in the catalog with *M* ≥2.5; respectively 50% for *M* ≥3 and 41% for *M* ≥4. Analyzing the declustered catalog, for each cell we get an independent estimate of the background level by dividing the number of mainshocks *N*_{main} by the total number of events in the cell, *μ*_{decl} = *N*_{main}/*N*. The dependence between the estimation based on the interevent-time distribution and that resulting from the declustering procedure is shown in Figure 3a for all cells and all combinations of *L* and *M*_{min}. A clear correlation between both estimations is observed, where the linear regression gives a dependence of *μ*_{inter} = −(0.06 ± 0.01) + (0.64 ± 0.01) · *μ*_{decl}.

Figure 3b shows the same analysis for 1000 simulations of the etas model and 1000 simulations of the stas model with parameters given in the Earthquake Simulations section. The correlation is found to be independent of the model type and very similar to that for the earthquakes in California. The scattering is in both cases similar, and the regression line calculated for the California earthquakes fits those of the stochastic simulations. Only for high mainshock levels do the simulations show a deviation from the linear trend that is less pronounced in the observed data.

The results in Figure 3 indicate that the analysis of the interevent-time distribution systematically leads to lower values of the mainshock rate than those expected from the declustering procedure. For the stochastic simulations, we are able to judge the quality of both estimations, because for each earthquake it is known whether it belongs to the Poissonian background rate or is triggered by a mainshock. Thus the real fraction of background events is known for each simulation and can be compared with the estimation based on the interevent times and that of the declustering process. This comparison is shown in Figure 4. It is found that the estimates based on the the interevent-time distribution scatter around the true value whereas the values based on the Reasenberg decluster algorithm are systematically overestimating the real level. Thus the chosen time and space windows of the declustering algorithm for California (*Q* = 10, *P* = 0.99, *τ*_{0} = 1 day, *τ*_{max} = 10 days) seem to be insufficient for declustering all aftershocks in the synthetic catalog.

### Robustness of the Estimations

The synthetic catalogs analyzed so far have been adapted in their length to those from California in order to enable a comparison. The number of events in each catalog varies between 50 and over 5000. To study the quality of the estimates in dependence on the number of earthquakes in the catalog, we separately analyze catalogs with 100, 1000, and 10,000 earthquakes. The results are shown in Figure 5. The scatter is rather large for 100 events, but is already small for 1000 events. However, we find that the scatter does not totally vanish for very large catalogs and the estimation is slightly biased toward underestimating the true value. This indicates that the gamma distribution is only an approximation of the true, more complicated interevent-time distribution. However, the systematic shift Δ can be corrected because of its approximate parabolic shape. Least-square fitting yields the term 4 which has to be added to 1/*β* in order to get an unbiased estimate (see bold line in Fig. 5). In this way, the fraction of mainshocks within the seismicity can be estimated by means of the gamma distribution with an accuracy better than 0.1 for catalogs with about 1000 events.

It is important to note that only the rate of mainshocks is really constant in the etas model. The mainshock fraction is only constant if it is measured on long time intervals. Sornette and Werner (2005) showed that this long-term average value should be less than 45% due to physical constraints. However, the mainshock fraction varies largely on time intervals that are small compared to the recurrence time of the largest earthquake. In time periods consisting of a large event, the aftershock fraction is usually very large, but small on periods with only small-magnitude background events. We now analyze different subintervals of the California catalog in order to demonstrate this effect as well as to quantify the stability of our estimations. In particular, we use five different data sets consisting of the earthquakes that occurred within (1) 1980–2004 (full catalog), (2) 1980– 2000, (3) 1980–1990, (4) 1990–2000, and (5) 1990–2004, respectively. We put a grid on the complete region of the catalog with a spacing of 50 km and chose those grid points where each of the five subcatalogs consists of at least 50 events with magnitude *M* ≥3 in the square region of size 100 × 100 km. Figure 6a shows that, in some regions, the average earthquake rate varies largely in the different time periods. At the same time, we find also that the estimated fraction of mainshocks is highly variable in many regions (Fig. 6b) but that the estimated rate of mainshocks is rather stable (Fig. 6c). This is in good agreement with the etas model. Furthermore, the estimated mainshock fraction is in most places less than 45%, as theoretically predicted by Sornette and Werner (2005). It is no contradiction that a number of places show larger fractions because all values are estimated on relatively short time intervals of maximum 24 years.

In order to demonstrate that the estimates are not only stable for different time periods but also coherent in space, we also analyze the complete catalog on a fine grid of 5 km spacing. For visual reasons, we use a circular region of 50 km radius around each grid point instead of a square region. Figure 6d shows the resulting map of the estimated mainshock rate where the largest values are found for the volcanic region of Long Valley caldera. The map shows a rather smooth behavior demonstrating that our estimations are robust.

## Reconstruction of the Mainshock Magnitude- Frequency Distribution

The interevent-time distribution cannot only be used to estimate the mainshock rate for a given magnitude cutoff, but also to reconstruct the full magnitude-frequency distribution of mainshocks.

The procedure is simple. Let us assume that we want to analyze an earthquake catalog of time length *T* and with minimum magnitude of completeness *M*_{min}. The first step is to choose the magnitude level *M* = *M*_{min} and to consider only those *N*(*M*) earthquakes with magnitude larger than *M*. In the second step, the mean and the variance of interevent times between the earthquakes with magnitude larger than *M* are determined. Together with the small correction term Δ (see equation 4), this gives an estimation of the mainshock fraction (*μ*_{inter} = + Δ) and thus of the number of mainshocks with magnitude larger than *M* per year: _{main}(*M*) = *μ*_{inter} · *N*(*M*)/*T*. Now, increase the magnitude level by a certain increment and go back to the first step and repeat the steps iteratively. In this way, we get an estimation of the magnitude-frequency distribution of mainshocks _{main}(*M*).

We perform this procedure for simulations, where the branching parameter *n* is fixed to 0.6 and the *a*-value to 4.3, whereas *α*-, *c*- and *p*-values randomly vary from simulation to simulation in the intervals defined in The Earthquake Simulations section. The catalog length is chosen so that the total number of earthquakes is equal to 1000 or 10,000. In these simulations, we choose the *b*-value of the mainshocks (*b*_{main}) independently of the *b*-values of aftershocks (*b*_{after}). For each simulation, we reconstruct the magnitude-frequency distribution of mainshocks in the described way. The resulting semilogarithmic distributions log((*M*)) are fitted by the Gutenberg–Richter relation log((*M*)) = *a* − *bM*. For that we perform a weighted least-square fit, where the inverse of the error was used as a weighting factor. The error interval of the activity rate is determined by the 10% and 90% quantiles of calculated for samples of *N* values randomly taken from a gamma distribution with parameter *β*. Additionally, we analyzed the mainshock distribution resulting from declustering and the real distribution of background events. In this case, the weights for the least-square fit are determined from the 90% confidence interval of observing an activity rate under the assumption of a Poissonian distribution. As an alternative, we have used the maximum likelihood method to calculate the *a*- and *b*-values (Aki, 1965). The results are almost identical.

In Figure 7, the estimated parameters of the mainshock magnitude-frequency distribution inverted from the interevent times are shown for two cases: (1) *b*_{main} = 0.8, *b*_{after} = 1.2, and (2) *b*_{main} = *b*_{after} = 1. These estimations are compared with those for the whole activity, the declustered catalog, and the real background events. In all cases, the interevent-time-based estimations are close to the real parameters. In particular, the estimated rates are unbiased and, in the first case, the small *b*-value of the background events is correctly identified. Only in the second case are the estimated *b*-values slightly too low. However, the same shift is observed for the estimations based on the declustered catalog and points to a more general problem of the *b*-value estimation.

The *a*-value can be directly determined from the reconstructed rate of mainshocks (*M*_{min}) and the *b*-value according to *a* = log((*M*_{min}) + *bM*_{min}. Thus, the described procedure yields an estimation of both parameters of the Gutenberg–Richter relation.

## Conclusions

The interevent-time distribution of earthquakes has been recently claimed to be universal and to indicate the existence of clustering beyond aftershock occurrence. By analyzing observed and synthetic earthquake catalogs, we find that both statements are invalid. Firstly, simulations of the etas model are shown to reproduce the observed interevent-time distributions in California without any fine-tuning of the model parameters. Thus the combination of a Poissonian background activity and triggered aftershocks according to the Omori law explain the observed shape of the interevent- time distribution. There is no need for assuming any additional mechanisms of long-term clustering such as, for example, accelerated or deaccelerated seismic activity preceding mainshocks (Wyss, 1997; Jaumé and Sykes, 1999). Secondly, we find that the shape of the interevent-time distribution is correlated with the number of earthquakes that remain after applying Reasenberg’s declustering algorithm. Thus the interevent-time distribution is not universal; in particular, it is correlated with the percentage of mainshocks in the catalog. However, exactly this nonuniversality of the interevent-time distribution offers the opportunity to extract local characteristics. Particularly, the analysis of the interevent times yields a new, independent, and nonparametric estimate of the mainshock rate and the whole mainshock magnitude-frequency distribution. The investigations presented here show that for the synthetic catalogs where the parameters are known, this estimation is better than the application of a standard declustering procedure. For the California catalog this would suggest that the Reasenberg algorithm seems to overestimate the background level, at least if standard parameters are used. Although the declustering procedure can probably be improved by fine-tuning the involved space and time parameters and adapting them to the specific seismic region, the interevent-time-based method seems to be superior because no parameters have to be adjusted. At least, it yields an important boundary condition for any applied declustering procedure.

## Acknowledgments

The authors wish to thank the editor Jeanne Hardebeck and two anonymous reviewers for valuable suggestions. This work was supported by the Deutsche Forschungsgemeinschaft (SCHE280/14-1).

- Manuscript received 8 March 2005.