Observing Spatio-temporal Clustering and Separation Using Interevent Distributions of Regional Earthquakes

Past studies that attempted to quantify the spatio-temporal organization of seismicity have defined the conditions by which an event and those that follow it can be related in space and/or time. In this work, we use the simplest measures of spatio-temporal separation: the interevent distances R and interevent times T between pairs of successive events. We observe that after a characteristic value R * , the distributions of R begin to follow that of a randomly shuffled sequence, suggesting that events separated by R > R * are more likely to be uncorrelated events generated independent of one another. Interestingly, the conditional T distributions for short-distance (long-distance) events, R ≤ R * (R > R *), peak at correspondingly short (long) T values, signifying the spatio-temporal clustering (separation) of correlated (independent) events. By considering different threshold magnitudes within a range that ensures substantial catalogue completeness , invariant quantities related to the spatial and temporal spacing of correlated events and the rate of generation of independent events emerge naturally.


Introduction
A better understanding of the processes governing seismicity, coupled with the advancements in instrumentation, resulted in the ability to measure and record the time of arrival, hypocentre location, and magnitude of an earthquake event with a substantially high level of sensitivity.Because of this, decades-long catalogues of earthquake events for different regions and even for a global scale are now available for study.Despite the current limitation in accurate short-term earthquake prediction, much information about the underlying properties of the processes generating earthquakes can be drawn from these historical records.In particular, these records can reveal how earthquake events cluster in space and time (Utsu and Ogata, 1995;Kagan and Knopoff, 1980).
Previous works have put forward different metrics by which this spatio-temporal organization in regional seismicity can be quantified.Davidsen et al. (2008) incorporated the spatial information in the criteria for identifying the recurrence of a particular event.When applied to seismicity, they observed scaling laws for both the spatial and temporal separation distances between recurrences (Davidsen et al., 2006).Baiesi and Paczuski (2004) included the magnitude information via the Gutenberg-Richter law in creating a spacetime window for aftershock collection, revealing a scale-free network of correlated earthquakes.Zaliapin et al. (2008) extended the analysis and revealed bimodal distributions of the generalized spatio-temporal separation of earthquake events.All of these works establish relationships between an event and one or more of other succeeding events, and use the data from the southern California seismic region, where extensive records exist for several decades of observation.
Here, we look at the simplest measures of separation distances between earthquake events.The interevent distances and times are the simple separation between two events that follow each other in time.Very recently, such approaches have been used to identify clusters of correlated earthquakes for various seismogenic regions (Anderson and Nanjo, 2013).Understandably, there are apparent limitations in using these simplest quantities to reveal spatio-temporal clustering.It relies heavily on the completeness of the catalogue under study (i.e. to ensure that no events are left out, thereby affecting the order of the recorded events) and fails to account for long-range correlations (i.e. it does not go beyond the succeeding event, and there is the possibility that non-related events may get in between correlated ones in the sequence).For these reasons, many previous Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
Perhaps the closest work to have shown clustering (separation) of correlated (uncorrelated) events using only the interevent times is one by Touati et al. (2009).They report observable differences between the return time distributions of regional and global earthquake catalogues: the histogram of interevent times of southern California earthquakes shows two distinct peaks, signifying the difference in characteristic waiting times between correlated (same aftershock sequence) and independent (different sequences) events, while global statistics reveals a single characteristic peak due to overlapping sequences from various locations.They explain the results using the epidemic-type aftershock sequence (ETAS) model (Kagan and Knopoff, 1981;Ogata, 1988;Sornette and Helmstetter, 2002), a five-parameter model wherein the main-shock generation rate µ is used as a proxy for spatial extent.While the work clearly demonstrated the different timescales involved for correlated and independent events similar to the bimodal behaviour observed by Zaliapin et al. (2008) for the same region, it does not incorporate the actual measure of spatial separation between earthquake epicentres in the statistical analysis.
In this work, we show that, despite the apparent limitations, it is possible to use the simple pairwise interevent distances and times to show spatio-temporal clustering of earthquakes.The advantage of using only the pairwise interevent separations, apart from the simplicity and the robustness to addition of data, is that there are no a priori assumptions regarding earthquake relationships and no regional factors need to be accounted for.Thus, we are able to show analyses of different regions with different catalogue completeness levels.Our analyses are guided by the fact that spatio-temporal clustering is a well-established phenomenon in seismicity (Utsu and Ogata, 1995;Kagan and Knopoff, 1980) and must therefore manifest even in the simplest measures of spatiotemporal separations.

Regional data sets and parameters
In our analyses, we used three catalogues from different regions: the Philippines, PH (1973PH ( -2012)), taken from the subset (4-24 • N and 115-130 • E) of the global Preliminary Determination of Epicentres catalogue (PDE, 2012); Japan, JP (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998), from the Japan University Network Earthquake catalogue (JUNEC, 1998);andsouthern California, SC (1982-2012), obtained from the Southern California Earthquake Data Center (SCEDC, 2012) and filtered to remove any man-made seismic events.As noted earlier, we need to ensure catalogue completeness for analyses that are based on simple interevent separations.For this purpose, we consider only events with magnitudes m above a certain threshold magnitude M, m ≥ M. We repeated the analyses for different M values in increments of M = 0.1 for a range of threshold magnitudes that ensure substantial data completeness.
To gauge the completeness levels of each data sets, we plot the cumulative distributions of magnitudes in the left panels of Fig. 1.Despite the expected power-law behaviour given by the Gutenberg-Richter (GR) law, we observe plateaulike "rollover" regions for smallest magnitudes resulting from limited resolution and other factors that hinder accurate recording of very weak events.For our purposes, therefore, the regimes following the GR law are deemed to be substantially complete as they are beyond these "rollover" regions.In addition to ensuring that the threshold magnitudes chosen are within the power-law regime, we ensure that the number of data points are enough to get statistically sound results in all our subsequent analyses.The ranges of threshold magnitudes satisfying these criteria are M ∈ [4.5, 4.8] for PH, M ∈ [2.5, 3.8] for JP, and M ∈ [2.5, 3.5] for SC.These are shown as the shaded regimes in the left panels of Fig. 1.
For each threshold magnitude in increments of M = 0.1, we obtain the time series of pairwise interevent times T i between events i and i + 1, denoted by where t is the time of occurrence in minutes relative to the first recorded event.To complement the temporal analysis, the corresponding interevent distance R i is defined as where the spatial coordinates are based on the latitude (φ) and longitude (θ) coordinates (in radians) and R E = 6371 km is the approximate radius of the earth.This definition of the interevent distance based on epicentres and assuming a spherical surface is a special case of the general hypocentre separation distance used by Kagan and Knopoff (1980).
Clustering behaviour is easily apparent from looking at the scatter plots of R and T .In the right panels of Fig. 1, we plot all pairs of R and T values for the highest threshold magnitude considered (plots for highest magnitudes are the least dense and are presented for better visualization, but the behaviour is the same for other threshold magnitudes).All scatter plots show a generally increasing trend, and closer inspection reveals dense concentration of points at both lower left and upper right regions.These scatter plots not only show the spatio-temporal clustering (concentration of points at short R and T region) and separation (at long R and T region) of earthquake events, but they also give a glimpse of the regional differences in the earthquake-generation mechanisms.
In the succeeding sections, we discuss the individual statistical distributions of both R and T .Using a characteristic R * obtained for each region, conditional distributions of T subject to the value of the corresponding R further highlight the clustering and separation behaviour, and the regional variability of these statistics.

Basic R distributions
We present the distributions of R and T both in the form of unnormalized histograms h(R) and h(T ) and normalized probability densities p(R) and p(T ).All distributions are presented in logarithmic binning for better visualization in double-logarithmic plots.The inclusion of unnormalized histograms is made in view of previous observations that some features of the distributions may not be noticeable upon normalization (Touati et al., 2009).We observe two regimes in both h(R) (left panels of Fig. 2) and p(R) (right panels of Fig. 2), similar to those observed in a previous work (Corral, 2006).These two regimes are visually discernible from the shape of the distribution, which is consistent even for different threshold magnitudes.Inspired by the similar result where a bimodal distribution of the combined spatio-temporal distance is observed for both the SC data and ETAS model results, we interpret the resulting distributions as a crossover between two component distributions having different characteristic lengths, where the component distribution having the longer characteristic length originates from a random process (Zaliapin et al., 2008;Zaliapin and Ben-Zion, 2013).To check whether this is the case, we generate a random sequence by reshuffling the events from the original series and comparing the R distributions of the original and shuffled data.Randomization is done by dividing the entire time of observation from the original data into time slices t = 0.1 s (for simultaneous events, one is moved to the next empty time slice) and randomly reshuffling these slices.In Fig. 2, the interevent distance distributions from the shuffled sequences are plotted as connected symbols alongside the corresponding distributions from the original sequences.Clearly, the shuffled distributions are unimodal, and the peak values correspond to the long-R peak of the original distribution.
It is interesting to note that all distributions collapse under a single curve by simple normalization (i.e. by number of events for the histograms and by number of events and bin widths for the probability density plots), hence eliminating the need for a region-dependent characteristic scaling factor to describe the data.Another interesting feature of the result is in the values of R where the original and shuffled sequences intersect, as denoted by the arrows in the left panels of Fig. 2. In Fig. 2g, we observe that these values have no apparent variation within the range of M considered.The average of these values is denoted as R * and marked by a broken line in the right panels of Fig. 2. The obtained values for the different regions are R * PH = 125 ± 1 km, R * JP = 164 ± 7 km, and R * SC = 79 ± 6 km, where the error denotes the maximum absolute difference between the values and the mean.
The emergence of R * by simply comparing the original and shuffled data is an indication that it represents a characteristic length scale in the system.The fact that its value differs per region suggests that it might be related to the underlying earthquake-generating mechanisms at work.As it represents a crossover separation distance between the strongly clustered events and the randomly occurring events, we believe that R * is an approximate measure of the spatial extent of triggering of related events.In the following, we use R * as a boundary between "short" and "long" separation distances.This allows us to study the relationship between the spatial and temporal interevent properties without imposing an artificially selected boundary.

Conditional T distributions
We divided the set of all interevent times into two groups based on the value of their corresponding R relative to R * -T in = {T |R ≤ R * } and T out = {T |R > R * } -and compare the conditional distributions of T in and T out with those of all T .Conditional distributions are important indicators of independence: if T is independent of R, both conditional distributions of T in and T out should follow the same behaviour, and collapse under the same curve upon normalization (Livina et al., 2005).Our results, however, point to a strong de-pendence between these two properties.We present in Fig. 3 the conditional histograms, h(T in ) and h(T out ) plotted with the total histogram h(T ) (left panels) and the corresponding conditional density functions p(T in ) and p(T out ) and the total interevent time probability density function p(T ) (right panels).In this case, we plot only the results for the smallest M considered, but similar behaviours are observed for all the other M values.
Both histograms and probability density plots show a strong dependence between spatial and temporal interevent properties.Without spatio-temporal clustering, the distributions of both T in and T out should just follow that of T .Instead, in Fig. 3 we see that the crossover between these two conditional distributions gives rise to the total distribution.For events with R ≤ R * , the h(T in ) shows a peak at a shorter characteristic value of the interevent time.Additionally, p(T in ) > p(T ) for short interevent times.This is a clear indication of the clustering of correlated events both in space and time.On the other hand, for events with R > R * , the h(T out ) shows a peak at a longer characteristic interevent time, and p(T out ) < p(T ) for shorter interevent times.Thus, T out represents the independently generated events, which are more likely to be separated by longer distances and times.
The difference between the conditional distributions and the total distribution can be viewed as a manifestation of the disparity in the timescales involved in the driving and relaxation mechanisms of earthquake events.The former, which we believe is responsible for the conditional histograms of long-range events, involves longer timescales, as it is driven by the slow process of tectonic motion (in the order of several cm yr −1 ), and results in significantly long waiting times before the generation of a new independent event.The latter may explain the origin of shorter waiting time durations for nearby events, as individual earthquakes in the same aftershock sequence happen in minutes, and entire sequences happen over a duration of several days or weeks.
The crossover between the mechanisms of clustering and separation is further highlighted upon looking at the behaviour of the distributions obtained from the shuffling procedure described earlier.As expected, shuffling the sequence will reduce the occurrence of both the very short and very long interevent times.The randomization would result in a significant decrease in the occurrence of temporal clustering of correlated events.Similarly, the shuffling procedure will remove the long temporal separation required to generate independent events.In Fig. 3, we observe this both in the histograms and probability density functions from the shuffled sequences, shown as broken lines in the plots.In the left panels, we observe that the histograms from the shuffled sequences peak around the weighted average of the peak values of the component distributions.In the right panels, we show the behaviour of the resulting distributions upon normalization.Though not readily apparent in a log-log scale, the distribution obtained from the shuffled series is found to be exponential, as expected of random, memoryless events.and (e)-(f) SC, M = 2.5.Hollow symbols are total distributions, thin broken lines are the distributions obtained from the shuffling procedure, and connected filled symbols are the conditional distributions.Using R * to separate the interevent times resulted in the the decomposition of the histograms into two component histograms with different characteristic times, as observed by Touati et al. (2009) in the ETAS model.The normalized probability densities further show that events separated by short (long) distances are also more likely to be separated by short (long) time intervals.The shuffled sequences show random (Poisson) statistics, with reduced occurrences of very short and very long interevent times.(g) The fraction of T out , γ , is roughly constant within the range of M considered: γ PH = 0.7; γ JP = 0.8; and γ SC = 0.4.
What is easily discernible, however, is that the probability of occurrence of very short and very long interevent times is significantly decreased upon shuffling.
The relationship between the total and conditional histograms in the left panels of Fig. 3 is reminiscent of the ETAS model results and empirical data analysis that show bimodal behaviour resulting from the crossover of the distributions of correlated (same aftershock sequence) and independent (different aftershock sequence) events (Zaliapin et al., 2008;Touati et al., 2009;Zaliapin and Ben-Zion, 2013).Touati et al. (2009) observed that the parameter that most influences the separation between the peaks of the correlated and independent histograms is the rate of the generation of independent events µ: lower µ results in a long separation of peaks while higher µ produces component distributions with almost overlapping peaks.In the left panels of Fig. 3, we observe that the disparity in the short characteristic peak of h(T in ) and the long characteristic peak of h(T out ) is more pronounced for the case of the PH (Fig. 3a) and SC (Fig. 3c).On the other hand, the peaks of the conditional and total histograms for JP (Fig. 3b) show almost overlapping peaks, suggesting a relatively higher level of seismic activity in the region.
While it is difficult to provide quantitative measures of their association, we observe that the level of seismic activity as measured by the overlap of the T in and T out histograms is qualitatively consistent with what is known about the regions considered.In particular, we expect the almost overlapping T in and T out histograms for Japan, as it is situated in a very active seismic region where at least four major plates are interacting (i.e. the North American, Pacific, Eurasian, and Philippine Sea plates) (Hirata, 1989).Previous works have indicated that the highly complex network structure of the fault system in Japan differs significantly from the almost one-dimensional strike-slip fault system in the Philippines (Matsumoto et al., 1992), which is also indicated from the differences in the histograms in Fig. 3a and c.On the other hand, previous studies have suggested similarities in the fault movements and structures of the Philippine and San Andreas faults (Rutland, 1967;Acharya and Aggarwal, 1980), which can also be seen in corresponding histograms Fig. 3a and e.
Here, we observe another emergent property of the system that is almost invariant within the range of magnitudes considered.Because the crossover separation distance R * is the same for all magnitudes, the fraction of events with R ≤ R * and with R > R * is preserved.We denote by γ the fraction of "long" distance events, γ = count(R > R * )/N, where N is the total number of events.In Fig. 3g, we obtain the following average values for this ratio: γ PH = 0.7, γ JP = 0.8, and γ SC = 0.4.Factoring in the threshold magnitude M and the total length of time considered, γ may in fact be related to the rate of background activity µ (Hainzl et al., 2006), which would explain the extent of the crossover between the component histograms as observed in the model (Touati et al., 2009).

Effect of higher threshold magnitudes
The spatio-temporal clustering of correlated events and the separation of independent events is still observed even for higher threshold magnitudes.Upon considering higher threshold magnitudes, weaker events are neglected from the analysis, thereby lengthening the mean waiting time between the occurrence of two "successive" events.Despite this, we still observe the separation of the temporal histogram into two conditional histograms based on separation distance.In Fig. 4, we show h(T in ) and h(T out ) side by side to compare their behaviour upon increasing M.
As shown in the left panels of Fig. 4, the effect of increasing M is to broaden the tails of the T in distributions.This is even more apparent when we look at the normalized p(T in ) in the insets, where the distributions all follow the same trend for several orders of magnitude before having distinct tails that grow longer with increasing M.This indicates that the T in distributions are more likely to be generated by correlated mechanisms; regardless of threshold magnitude, the shorter characteristic waiting times between correlated events should follow the same behaviour.Conversely, in the right panels of Fig. 4, we observe that the distributions of T out are completely shifted to longer waiting times by increasing M.This strengthens the fact that T out distributions are due to independent events, which are guided by the same random mechanisms but have increasing characteristic values for increasing M.
We plot in Fig. 4g-i the peak of the histograms for different regions.For increasing M, the peak of the T in histograms shows negligible change, in contrast with the fast rate of increase in the peak of the T out histograms.The average value of the h(T in ) peaks, τ , is of particular interest to us, as it represents another quantity that emerged that is very slowly varying for all the magnitude thresholds considered.The average values (τ PH = 131 min; τ JP = 89 min for JP; and τ SC = 9 min for SC) may be related to the characteristic waiting time between correlated aftershocks.

Conclusions
Spatio-temporal clustering of earthquakes in the form of foreshocks and aftershocks is not entirely new and is in fact evident even from experience.Therefore, this should manifest even in the simplest metrics that are readily available from the data.As shown here, within a range of threshold magnitudes that ensure substantial data completeness, simply taking the events pairwise based on their arrival times is enough to reveal these features of seismicity.The analysis relies heavily on the completeness of the catalogues, but this is easily addressed by substantially complete records for different regions and relatively long observation times.The main advantage of the method is the lack of predefined conditions to characterize the relationships between successive events.Comparison with the distributions obtained from randomly shuffled sequences reveals important characteristics of the system in space and time.The characteristic separation distance R * obtained from the intersection of the interevent distance distributions of the original and shuffled sequences has been shown to be constant for the range of magnitudes considered.In the temporal domain, the interevent time distributions of the shuffled sequences show a decrease in occurrence of both the very short waiting times (due to correlated earthquakes) and the very long waiting times (the energy accumulation period to generate independent events) compared to those of the original sequences.
Grouping successive events based on their separation distance relative to R * have resulted in the decomposition of the interevent time distribution into two components with different characteristic times: events separated by "short" ("long") distances are also more likely to be separated by "short" ("long") waiting times, suggesting spatio-temporal clustering (separation) between correlated (independent) events.We believe that the observed spatio-temporal clustering is strongly dominated by aftershock sequences, particularly because we are looking at significantly high-magnitude events that are more likely to trigger aftershocks, but we do not discount other possible sources of clustered events (e.g.earthquake swarms; Weaver and Hill, 1978).
From the simple analysis, we obtained characteristic values that may be related to the physical properties of the regions: (1) the R * may be related to the spatial extent of correlated earthquake triggering; (2) the peak of the T in histogram, τ , may be related to the corresponding characteristic waiting times between correlated events; and (3) the fraction of events with R > R * , γ , may be related to the rate of background activity µ.These values, which emerged naturally from the analyses, may have resulted from the confluence of many factors pertaining to the earthquake-generating mechanisms specific to each region.However, in the absence of the capability to disentangle the contribution from each of these factors, these metrics provide useful approximations for characterizing and comparing these systems.The method presented here, with its simplicity and robustness to addition of new events, can easily be adapted in understanding other catalogues to reveal simple estimates of spatial and temporal scales involved in regional seismicity.

Figure 1 .
Figure 1.Cumulative magnitude distributions and representative R-T plots for the (a)-(b) Philippines, PH; (c)-(d) Japan, JP; and (e)-(f) southern California, SC.The threshold magnitudes considered for analyses are shaded in (a), (c),and (e).These are within the power-law regimes of the magnitude distributions and have a sufficient number of events to ensure substantial completeness.The R-T scatter plots in (b), (d), and (f), obtained for the highest threshold magnitudes considered for each of the catalogues, reveal a visually discernible separation into two clusters situated at different regimes: for short-R-short-T and long-R-long-T .

Figure 2 .
Figure 2. Interevent distance R histograms (left panels) and probability density (right panels) plots for (a)-(b) PH, (c)-(d) JP, and (e)-(f) SC, for the corresponding magnitude ranges shown in Fig. 1.Symbols are distributions obtained from the original sequences, while connected symbols are from the shuffled sequences.The original and shuffled distributions begin to follow the same trend after characteristic R * values, as indicated by the arrows in (a), (c), and (e).In (g), the value of R * is shown not to vary significantly within the threshold magnitude ranges considered.The average values of R * are shown as broken lines in the right panels: (b) R * PH = 125 ± 1 km, (d) R * JP = 164 ± 7 km, and (f) and R * SC = 79 ± 6 km.These R * values are used to separate the "short" and "long" R regimes.

Figure 3 .
Figure 3. Representative interevent time T histograms (left panels) and probability density (right panels) plots superimposed with the conditional distributions of T in and T out , for the case of the smallest M considered: (a)-(b) PH, M = 4.5; (c)-(d) JP, M = 2.5; and (e)-(f) SC, M = 2.5.Hollow symbols are total distributions, thin broken lines are the distributions obtained from the shuffling procedure, and connected filled symbols are the conditional distributions.Using R * to separate the interevent times resulted in the the decomposition of the histograms into two component histograms with different characteristic times, as observed byTouati et al. (2009) in the ETAS model.The normalized probability densities further show that events separated by short (long) distances are also more likely to be separated by short (long) time intervals.The shuffled sequences show random (Poisson) statistics, with reduced occurrences of very short and very long interevent times.(g) The fraction of T out , γ , is roughly constant within the range of M considered: γ PH = 0.7; γ JP = 0.8; and γ SC = 0.4.

Figure 4 .
Figure 4. Conditional histograms of T in (left panels, hollow symbols) and T out (right panels, filled symbols), with the corresponding probability density plots (inset): (a)-(b) PH, (c)-(d) JP, and (e)-(f) SC.For higher M, events with smaller magnitudes are neglected; this results in the broadening of the tail of T in histograms (left panels) and the complete shifting of the T out distributions to longer values (right panels).In (g)-(h), the peak interevent time for all histograms are tracked, showing the almost negligible shift of the peaks of T in and the continuous increase in the peaks of T out .