Investigation of the Spatial Clustering Properties of Seismic Time Series : A Comparative Study from Shallow to Intermediate-Depth Earthquakes

In this paper, a size-independent modification of the general detrended fluctuation analysis (DFA) method is introduced. With this modified DFA, seismic time series (m ≥ 4.5) pertaining to most seismically active regions of the world from the year 1972 up to the year 2016 are comparatively analyzed. An eminent homogeneity of spatial clustering behaviors in worldwide range is detected and DFA scaling exponents coincide with previous results for local regions. Furthermore, universal nontrivial spatial clustering behaviors are revealed fromshallow to intermediate-depth earthquakes by varying the depth of the hypocenter: in shallow earthquakes, the depth range corresponding to the minimum spatial clustering is first verified to be spatiotemporal and identical; in intermediate-depth earthquakes, strong weakening long-range correlations compared to shallow earthquakes are unveiled to strengthen the recent findings in seismology. Our work aims to suggest a possible statistical approach to explore the dynamic mechanisms implied in the seismicity as well as in other analogous dynamic evolution processes.


Introduction
Among multifarious natural hazards, seismicity is a representative complex spatiotemporal phenomenon, in which underlying mechanisms are not yet fully understood [1,2].Statistics analyzing provides important insights into fundamental properties of seismic processes, such as the well-known Gutenberg-Richter laws [3] for the magnitude distribution of earthquakes and the Omori laws [4] for the time evolution of the frequency of aftershocks.Inspired by these fractal analogies, Bak has invoked the self-organized criticality (SOC) concept [5] into seismicity.For reviews on this topic, see [6][7][8] and references therein.
Usually, in statistical analysis, earthquakes can be treated as spatiotemporal stochastic point processes [9,10] due to their nonhomogeneous occurrence time.In such processes, spatial and temporal clustering behaviors featured by correlations come to light and develop into important properties of seismic activity [11][12][13][14][15][16][17][18][19][20].In seismology, such correlations originate from multifarious microscopic dynamics, governing complex interactions of many components over a worldwide range of time or space scales.It is known that an earthquake can trigger the subsequent one even more than 1000 km far away [21,22], implying that the long-range correlations exist.In such a situation, identifying the properties of long-range correlations may shed light on the fundamental dynamic seismic mechanisms.This approach has been continuously attracting seismologists, geologists, and physicists' attentions [22][23][24][25][26][27][28][29][30][31].
However, identifying long-range correlations in seismicity is usually still a tricky business mainly due to the irregular fluctuations and severe nonstationarity embedded in the emitted seismic signals [25,32,33].Such nonstationarity is of a much larger order than that of any other known natural signal.Peng et al. [34] propose the DFA as a robust method which can provide a quantitative representation of the correlation properties of a series of signals.The advantage of DFA over many other methods is that it permits the detection of long-range correlations embedded in seemingly nonstationary time series and also avoids the spurious 2 Complexity detection of apparent long-range correlation which is an artifact of nonstationarity.In seismology, the earthquakes can be described as points represented by the geophysical coordinates (longitude and latitude) in space and a sequence of points distributed for time axis where such temporal points are determined by their occurrence time.Accordingly, the seismic time series can be initialized by the spatial or temporal interval between subsequent earthquakes for the analysis of seismic mechanisms [13,27,29,[35][36][37][38].With the aid of DFA analysis on these series, scaling behaviors of fluctuations can be explored to qualify the earthquakes correlation characteristics.Already, there have been many works performed in this approach: Telesca has analyzed the spatiotemporal clustering behavior in earthquakes, especially the dependence of the time-clustering behavior on the depth (0-30km) of the earthquakes in southern California (SC) [36]; Z. Y Zheng et al. have found DFA scaling exponent of interval time series has no dependence on hypocenter depth or epicentral distance in Japan earthquakes [25]; P. A Varotsos and N. V Sarlis have discovered that the minima earthquake preceding violent earthquake is followed by characteristic changes of temporal correlations between earthquake magnitudes [39].Other related results can be found in [40][41][42][43][44].
Recently, the studies of the global earthquake network strengthen the idea that the Earth behaves as a self-organized critical system with long spatial correlations between earthquakes [45,46].This means the earthquake spatial clustering may be not only regional but also global.Nevertheless, to the best of our knowledge, most of the current works which are probing earthquakes spatial clustering with DFA method have been carried out in the absence of worldwide comparative study.In this work, a comparative DFA analysis will be performed to explore the possible consistency of spatial clustering behaviors across the seismic regions around the world.This work seems to contribute to the idea of a possible worldwide spatial long-range correlation as well as a long-range temporal correlation, i.e., "memory" [22,25,[47][48][49][50] between earthquakes temporally apart from each other.In addition, the spatial long-range correlation can be quantified as a scaling exponent through DFA analysis on the time series of earthquake separation distance (spatial interval), whose fluctuation is obviously governed by the seismic region dimension.Such effect may lead to nonuniformity in comparative study across different regions of the world.Here, we develop a size-independent modification of the general DFA method to remove system dimension effect from the fluctuation scaling behavior.This modification is evoked from our analysis of correlations in the Bak-Sneppen model [51], which is one of the simplest SOC models for earthquakes, forest fires, biological evolution, etc.
There is also another theme that appeals to us: is there any universal difference in spatial clustering behavior between shallow and intermediate-depth earthquakes?This doubt is "echoed" with one of the current hot debates in seismology: what constitutes an intermediate-depth earthquake [52][53][54][55]?It is well-known that the shallow earthquake stems from brittle and ductile behavior of the crust [56].However, the physical mechanisms behind intermediate-depth earthquakes remain uncertain (two plausible hypothesis: dehydration embrittlement [57] and thermal shear runaway instability [52,58]).In particular, the depth cut-off to distinguish between shallow and intermediate-depth earthquakes is still difficult to define [53].On this account, a depthdependent spatial clustering analysis is implemented with the depth range extending from the shallow seismic region to the intermediate-depth seismic region.One expects that such results may provide valuable clues on detection of the intermediate-depth earthquake mechanisms as well as help explain differences in their seismological behaviors compared to shallow earthquakes.
This paper is organized as follows: In the next section, we review the general DFA method and introduce our modified DFA.In Section 3, the earthquake catalogs and regions are presented.In Section 4, we comparatively analyze the spatial clustering and its depth-dependence with the modified DFA.The last section is devoted to a summary.

Description of the Method
DFA is a scaling analysis method providing a quantitative representation of the correlation properties of a series of signals.With this method, the long-range correlation can be identified from a nonstationary time series.Such superiority has been confirmed in many research fields [59][60][61][62][63].The general DFA procedure consists of steps as follows.Let us suppose that () is a length series , and this series is of compact support, i.e., () = 0, for an nonsignificant fraction of the values only.First, we assume that () are increments of a random walk process around the average  V and integrate the time series to acquire the "profile" or "trajectory" of the signal: where Then () can be viewed as the "displacement" from the origin due to a set of random "walks" () −  V , because of the integration that the correlations in () get involved.Next, the integrated time series is divided into boxes of equal length .In each box of length, a least-squares line is fit to the data, representing the trend in that box.In DFA, one fits in each box the integrated time series to polynomial function which is the integration of the supposed trend.Here,   () is the fitting polynomial in box .Linear, quadratic, cubic, or higher order polynomials can be used in the fitting procedure (usually called DFA1, DFA2, DFA3, . ..).The root-mean-square (rms) fluctuation of this integrated and detrended time series is calculated by Repeating this calculation over all box sizes, we obtain a relationship between () representing the average fluctuation as a function of box size and the box size  representing the number of spikes in a box which is the size of the watch window.If () behaves as a power-law function of , data presents scaling: where the scaling exponent  provides information about the type of correlation of the signal.If there is no correlation (white noise),  = 0.5.Any long-range correlation in the signal will give  ̸ = 0.5.An exponent  > 0.5 corresponds to positive correlations in the signals, while  < 0.5 represents negative correlations.
As mentioned above, earthquakes spatial clustering properties can be studied by analyzing the long-range correlations in earthquakes sequences, which are defined in terms of the time series of the separation distance [14].In this work, we regard the sequence of the earthquakes as a stochastic point process, displaying earthquake which occurs at some random location (hypocenter) in three-dimensional (3) space.Such process can be described by Dirac's delta function finite sum to the location ( 1  ,  2  ,  3  ): presents the number of earthquake events.The three-dimensional separation distance between two consecutive earthquakes  and  is measured as Euclidean distance between their hypocenters: where the Earth radius is  0 = 6.36 × 10 6 .(  ,   ) are the latitude and longitude, in radians.ℎ  is the depth of the hypocenter of th earthquake in the earthquake data catalog.
The range of the location distance in each dimension is then defined as As pointed above, earthquakes sequence can be featured by 3 time series of the separation distance   .The generalization of the DFA to 3 records is straigthforward as one needs only to change the basic equation ( 3), which is rewritten as where y() =  1   1 +  2   2 +  3   3 , similarly for the fitting polynomial vector y  .So far one gets the rms fluctuation of the hypocenter separation distance in general DFA method.Considering the range of the location distance, y() ("displacement") is the order of .The fitted local trend y  () is expected to have a similar order of magnitude.If the lℎ order polynomial fitted curve in each box can cancel the  order part in the "displacement", F() will more likely scale with  1/2 .Therefore, the domination of seismic region dimension on the rms fluctuation could be removed by the scaled F().Another considerable scaling component with  is the appearance of earthquake hypocenters, with fractal dimension   [12,31,64].Thus, more reliable scaling renormalization of F() is and, for convenience, we chose the single   = 2.6 in our calculation [12], although, precisely,   varies slightly for different seismic areas and depths [64].Consequently, one can employ this modified DFA to uniformly compare the scaling behaviors of the   fluctuation in different earthquake catalogs.Also, it should be noted that DFA analysis is inappropriate regardless of the order of the polynomial fit applied, because of the presence of nonstationary in time series of the earthquakes separation distance [32].The detrending of the time series is carried out by the subtraction of the polynomial fits from the "profile"; different order DFA differ in their ability of eliminating trends in the series.For instance, DFA0 can remove just constant trends in the data, DFA1 wipes off constant and linear trends, DFA2 eliminates trends up to second order, and so forth.Thus a comparison of the results obtained from different orders of DFA provides information on the order of trend in the data.In this work, we will present DFA analysis up to the third-order DFA3.

The Regions and Catalogs
One of the two earthquakes catalogs used in our analysis is derived from the Global Earthquake Catalog, provided by US (United States of America) Geological Survey (USGS), exactly at Advanced National Seismic System [65], in which seismicity from the entire Globe has been recorded.The catalog contains all seismicity completely from January 1, 1972, to December 31, 2016.However, such a limitation in this catalog-that it is not consistent in all regions of the world-merely comprises earthquakes of full magnitudes for the US but only events whose magnitude ≥ 4.5 (of the Richter scale) for the other areas in the Globe (only if they are informed that the event is clearly felt or causes damage).In seismology, the homogeneity of a global catalog is indispensable with a view of the long-range seismic activity in the world.Thus, in order to obtain homogeneous results from our comparative analysis for different seismic regions over the world, we have screened out only earthquakes with magnitude ≥ 4.5.The magnitudes Mb, ML, Ms, and Mw are taken into account accordingly and the artificial earthquake events (nuclear tests and quarry blasts) have been excluded from the catalog.Finally, 209560 events remain, and 80% of them happen near the surface of the Earth (depth ≤ 100 km).Besides, to compare our results with previous studies, we adopt the other earthquakes catalog: southern California (SC) earthquakes catalog relocated by Richards-Dinger and Shearer [66], which is widely employed for its high quality relocation.
The regions chosen in our study are all regions of high seismic potential.The Himalayas are among the most seismically frequent regions in the Globe.Most of the earthquakes there are contributed by the Eurasian plates and north-south convergence of the Indian Eurasian plates, the east-west convergence of the Indo-Burmese mountain, and the Indian plate underthrusting below the Eurasian plate [67].Japan is entangled with huge seismic risk, due to the Pacific-Philippine-Eurasian triple plate junction, where the interactions of three tectonic plates are complex and potential activity is involved [68].South Asia strides over subduction zone on the eastern side of the India plate, the strike slip structural belt at the back edge of the subduction zone, and the compressed structural zone of the India ocean plate and the pacific plate [69].The large multicentury earthquakes result from a broad collisional zone where earthquakes with magnitude ≥ 7.0 are common inland.Most of the shallow major earthquakes and all the other intermediate-depth and deep earthquakes, outside of the Circum-Pacific seismic belt, are distributed here.The Americas (North and South) are surrounded by the Circum-Pacific seismic belt, from the south of Alaska, along the eastern coast of North America, through Canada, California, and western Mexico to Columbia, Peru, and Chile in South America.This wide region is also the most seismically active zone in the world [70].The geological range of the above regions and corresponding number of earthquake events in our analysis are illustrated in Table 1 and the geographic regional maps are shown in Figure 1.

Results
Firstly, the spatial clustering properties are analyzed from the separation distance time series of the catalogs.To see whether our modification is valid, we plot () as a function of box size  in Figures 2 and 3, by employing general and modified DFA up to the third order, respectively.The obvious long-range correlation of separation distances can be found with scaling exponents  greater than 0.5 in both Figures 2 and 3. Furthermore, the scaled rms fluctuations in Figure 3 have weak dimensional dependence in contrast with rms fluctuations in Figure 2, which not only verifies our modification satisfactorily but also pretty indicates such long-range correlations come from the earthquake dynamical fluctuations rather than the finite dimension effect.In Figure 3,  probably equals 0.70 − 0.76, consistent with the previous results of the local seismic regions [41,71].The differences of  between the orders of DFA are so small with the interval of uncertainty ±0.05, confirming the precision of our algorithm.Hence, it may be possible for the persistent existence of long-range correlations of separation distances across the Globe as well as in the local seismic regions.
Besides, in Figure 3, scaling exponents  in different catalogs are similar: South Asia: 0.73±0.02,Himalayas: 0.70± 0.03, America: 0.71 ± 0.01, Globe: 0.72 ± 0.05, and Japan: 0.76 ± 0.03.Such almost identical power-law behaviors of () mean some homogeneity of spatial clustering behavior in worldwide range seems already apparent, at least with the magnitude ≥ 4.5.Seismologically, this homogeneity of hypocenters spatial clustering could be considered as a manifestation of the fault network structure and the nature of earthquakes long-range interactions across the Globe.Thus, our results are inclined to have a better agreement with the conjecture that the Earth behaves as SOC system with possible spatial long-range relationships between distantly located earthquakes [5,45,46].Meanwhile, what needs to be explained here is that the local spatial clustering [11,12] cannot be effectively excluded from our  exponent due to the absence of the catalogs declustering procedure as a multitype algorithm widely used in seismology [72].Specifically, aftershock sequences responsible for local spatial clustering can be identified and removed in declustered catalogs.This means a sufficient work is needed in the future to fully discriminate between long-range or local correlation.
Next, we investigate depth-dependent spatial clustering properties from the catalogs same as above, with the depth ranging from 0 to 230 km.We first divide the depth range into 23 classes of 10 km each; then we apply the modified DFA analysis on the series of the separation distances   belonging to each class.To ensure the efficiency of modified DFA, we check up the adequacy of the 23 depth classes.In Figure 4, we show  as a function of hypocenter depth for different catalogs.From these figures, one can catch sight of the fact that not all the depth classes are filledin for each catalog, due to the inherent event (earthquake) undersampling of the corresponding deep regions in catalogs (Japan and Himalayas).After that, in every catalog, we calculate  for each depth class, associating the value of  with the average depth of the class.As a result, phased nontrivial behaviors of  in Figure 4 can be observed.In the earlier depth range of 0 − 40 km (0 − 30 km for Himalayas),  shows a singular behavior irrespective of the catalogs.In this valley-like structure,  plummets from about 0.9 to the assumed minimum value approximately 0.7 in depth: 15.68 km (South Asia), 15.23 km (Himalayas), 15.28 km (America), 15.26 km (Globe), and 15.59 km (Japan), respectively.It should be noted that these depths are the average depth of the earthquakes in the class 10 − 20 km.In order to obtain a possible more precise depth range corresponding to the minimum , the modified DFA analysis as above is repeated with the depth range focused on 0 − 40 km and 8 classes of 5 km each.The higher depth resolution is not considered due to the depth measurement error which results from the fact that recording stations are confined to the Earth's surface; the exact depth of the hypocenter is very hard to accurately resolve.Nonetheless, Figure 5 shows exactly similar valley-like structures of  for different catalogs.The minimum  occurs at depth: 17.85 km (South Asia), 17.67 km (Himalayas), 17.76 km (America), 17.74 km (Globe), and 17.77 km (Japan), respectively.Such extremely unanimous depth range (15 − 20 km) for minimum  would remind us of   whether a determined universal depth range corresponding to the minimum  exists.To further prove this, we apply the same calculation to the SC catalog, relocated by Richards-Dinger (2000).The small depth error (741 m) in this catalog enables us to analyze  in 1 km class.In Figure 6, the valleylike structure exhibits again, associating a minimum  at the depth 15.49 km.More interestingly, in Telesca's results from SC catalog [36], the scaling exponent of temporal clustering shows a minimum corresponding to the depth range of 15−16 km.Worthy of note is that similar range of  has been used to calculate  in Figures 4-6 and each () as a function of  presents a good scaling.Although an exhaustive study should be further performed, it seems that the spatial clustering, as well as the temporal clustering of shallow earthquakes, shows a minimum  corresponding to the depth of 15 − 20 km.This minimum could be derived from the brittle and ductile behavior of the crust versus depth, as presented in the synoptic shear zone model [73].Therefore, we argue that the clustering behavior as a typical seismic dynamic evolution may be characterized with spatiotemporal consistency.One can further find in Figure 4 that  sustains falling with the increasing depth in all catlogs, until the next nontrivial depth region is approaching, and in each catalog  almost gradually flattens after depth 50 km.It is worth noting that 50 km is the depth at which well-developed aftershock sequences become infrequent [74] and also the maximum depth of  interplate subduction zone earthquakes [75].Either of the above analogies accounts for the weakening long-term correlations in deep earthquakes.After transition depth (50km), exactly beyond the depths 54.98 km (South Asia), 55.09 km (America), 64.88 km (Japan), 65.19 km (Himalayas), and 74.94 km (Globe),  presents smooth fluctuations ranging within 0.52 − 0.62.Such "smoothness" indicates the presence of tendency to flicker-noise type dynamics in the earthquake process [36,71].On the other hand, these closely related but not exactly the same depths may correspond to the fact that it is hard to define a depth cut-off to distinguish between shallow and intermediate-depth earthquakes.For example, Wadati (1929) defines shallow earthquake as ≥ 60 km and intermediate-depth earthquake as 100 − 200 [76]; Gutenberg and Richter (1949) define intermediate-depth earthquake as 70 − 300 [77]; Frohlich (2006) and Houston (2007) define intermediate-depth earthquake as 50 − 300 km [74,78].These inconsistent definitions as well as our results imply that shallow and intermediate-depth earthquake dynamic mechanisms may not be differentiable at least in the depth range 50 − 70 km.This puzzle, as one can see, remains, to us, regrettable but predictable for two conceivable reasons: one is the authentic existence of strong interactions between the occurrence of intermediate-depth and shallow earthquakes in recent seismological research [79]; the other is the disparate seismotectonics between seismic regions, e.g., Himalayas dominated by the shallow earthquakes but the South Asia frequently invaded by intermediate-depth and deep earthquakes.Nevertheless, the significant and Complexity synchronized suppression of  in different catalogs shows that the spatial clustering of intermediate-depth earthquakes is widely different from that of shallow earthquakes.These low and stable spatial clustering properties in intermediate-depth earthquakes may originate from two plausible hypotheses: dehydration embrittlement and thermal shear runaway instability.Based on them, a particular kind of intermediate-depth earthquakes clustering could be regarded as a earthquake nest [80] which is isolated from nearby seismic activity [53].Such nest mode would lead to very small rupture regions and/or slow rupture velocities compared with shallow earthquakes [54,55].In a word, the spatial long-range correlations (perhaps the temporal ones as well) may become strong weakening in intermediate-depth earthquakes, relative to a cascading thermal shear runaway, with a leading inefficient radiation process followed by a quick and active efficient rupture [81,82].Meanwhile, in Figure 4, the error and instability of  are distinct beyond the depth about 170 km, associating with the event undersampling in catalogs.Even so, any uniform earthquakes long-range dynamic mechanisms in such deeper regions could not be precluded on account of our results, because some possible earthquake nests may be found or confirmed as the extensive regional and global networks [53].Although our results are not sufficient, it may open up a intriguing prospect for the exploration of the dynamic correlation mechanisms in seismicity.

Summary
Analyzing the earthquakes clustering behaviors gains insight into the underlying dynamics of seismicity.In the present paper, we introduced a modified version of general DFA method to demonstrate the extremely similar spatial clustering behaviors for different wide seismic regions over the world.Some kind of homogeneity involved worldwide in the spatial dynamic seismic evolution is shown.
Moreover, the depth-dependent spatial clustering properties have been identified on the comparison of different catalogs.A phased nontrivial clustering behavior has been revealed from shallow to intermediate-depth earthquakes.In the shallow regions,  has been identified with a minimum corresponding to the depth range of 15 − 20 km.Comparing with the previous results in the temporal clustering analysis, we argue that a consistent characteristic depth range may exist in both spatial and temporal clustering behaviors of shallow earthquakes.In the intermediatedepth regions, earthquakes spatial clustering has been found faded with the much lower and stable value of , which means the strong weakening spatial long-range correlations in intermediate-depth earthquakes compared to shallow earthquakes, probably due to the earthquake nest mechanisms.
We should admit that our results could not preclude other clustering behaviors in intermediate-depth earthquakes, because the extensive regional and global networks of earthquake nests may be truly natural.Such evidence may be supplied by more comprehensive statistical correlation analysis of seismic time series.We hope to report our progress in this regard in the future.

Figure 2 :
Figure2: log 10 ()-log 10  plots for the time series of hypocenter separation distance in catalogs: South Asia, Himalayas, America, Globe, and Japan.The employed method for calculations is the general DFA.

Figure 3 :Figure 4 :Figure 5 :
Figure3: log 10 ()-log 10  plots for the time series of hypocenter separation distance in catalogs: South Asia, Himalayas, America, Globe, and Japan.The employed method for calculations is the modified DFA.

Figure 6 :
Figure 6: Variation with the depth of , calculated by the modified DFA for the time series of hypocenter separation distance in SC catalog relocated by Richards-Dinger and Shearer.

Table 1 :
Geological range (latitude and longitude) and earthquake event number of the regions.