On the Statistical Signiﬁcance of the Variability Minima of the Order Parameter of Seismicity by Means of Event Coincidence Analysis

: Natural time analysis has led to the introduction of an order parameter for seismicity when considering earthquakes as critical phenomena. The study of the ﬂuctuations of this order parameter has shown that its variability exhibits minima before strong earthquakes. In this paper, we evaluate the statistical signiﬁcance of such minima by using the recent method of event coincidence analysis. Our study includes the variability minima identiﬁed before major earthquakes in Japan and Eastern Mediterranean as well as in global seismicity.


Introduction
In many fields of science, and especially in the study of complex systems, a very important problem is the possible interconnection between two point processes. In simple words, when we have two event time-series how can we decide that they influence one another? To answer this crucial question, Donges et al. [1] recently introduced the event coincidence analysis (ECA) which takes the general viewpoint of possibly coupled point processes and quantifies the strength, directionality, and time lag of the statistical relations between these two event series. This is done by comparing the empirical co-occurrence frequencies with those one may expect for a null hypothesis model, e.g., uncoupled Poisson processes [2]. ECA has already found useful applications in various disciplines from epidemiology [1] and biogeography [3] to environmental and sustainability science [4], hydrology [5,6], geophysics [7,8], or even complex network physics [2]. Here, we employ ECA together with natural time analysis (NTA) to attack a problem related with the physics behind the preparation of strong earthquakes (EQs) which are one of the deadliest natural disasters.
The present paper is organized as follows. In the next Section, we briefly summarize the background of NTA (Section 2.1) and ECA (Section 2.2). Section 3 presents the results obtained for both regional (Sections 3.1 and 3.2) and global (Section 3.3) seismicity. These results will be discussed in Section 4 and our conclusions follow in Section 5.

Natural Time Analysis
For the NTA of a time-series comprising N EQs, the natural time is defined by χ k = k/N and serves as an index for the occurrence of the k-th EQ. The quantity χ k together with the energy E k released during the k-th EQ of magnitude M k , i.e., the pair (χ k , E k ), is studied in NTA. Alternatively, we study the pair (χ k , p k ), where denotes the normalized energy released during the k-th EQ. The variance of χ weighted for p k , labeled κ 1 , is given by [9,11,26,63]: where E k , and therefore p k , for EQs is estimated through the usual relation suggested by Kanamori [64]: and M k corresponds to the moment magnitude M w [65]. The magnitude reported in each EQ catalog should be converted to M w ; for example, in the case of Japan the magnitudes reported from the Japan Meteorological Agency (JMA) M J MA were converted to M w using the formulae suggested by Tanaka et al. [25]. As already mentioned, the quantity κ 1 , i.e., the variance of natural time χ k , plays an important role in seismicity. Varotsos et al. [26] proposed that κ 1 given by Equation (2) can be considered as an order parameter for seismicity since its value changes abruptly when a strong EQ occurs and its fluctuations have statistical properties similar to other nonequilibrium and equilibrium critical systems [16,26]. When considering an excerpt of an EQ catalog comprising W consecutive events, we can estimate various κ 1 values that correspond to subexcerpts of consecutive EQs, e.g., by using the first 6 EQs, the second 6 EQ, the first seven EQs, etc. This multitude of κ 1 values enable the calculation of their average value µ(κ 1 ) and their standard deviation σ(κ 1 ). We then determine the variability β of κ 1 [56]: that corresponds to this natural time window of length W and quantifies the intensity of the fluctuations of the order parameter of seismicity. Note that β W of Equation (4) could be identified [62] as effectively the square root of the Ginzburg criterion, e.g., see Equation (6.25) on p. 175 of Goldenfeld [66], the importance of which in EQ processes has been discussed by Holliday et al. [53]. The time evolution of β W is followed by sliding the window of W consecutive EQs, event by event, through the EQ catalog and assigning to its value the occurrence time of the EQ which follows the last EQ of the window studied in the EQ catalog. The fluctuations of the order parameter of seismicity as quantified by β W have been studied in various seismically active regions, like Japan [30,31,59] or Eastern Mediterranean [62], and in the global seismicity as well [60,61]. These studies have revealed that before strong EQs β W s exhibit characteristic minima labeled β min and were focused, so far, on a small number of the stronger EQs in each geographic region (see Tables 1-5) for which their preceding minima exhibited a clear characteristic behavior. Moreover, in the β W computations in each region only EQs with magnitude exceeding a certain value have been considered to ensure data completeness [30,31,[59][60][61][62].

Event Coincidence Analysis
In ECA is based on counting (possibly lagged) coincidences between events of the two different types but it should be stressed that B type events are considered as possibly influencing the timings of A type events, and not vice versa [1,3]. The assumption to be quantified and tested by ECA is that the events of B precede the events of A. This is made by introducing an instantaneous coincidence if two events with timings t A i and t B j (t A i ≥ t B j ) are closer in time than a coincidence interval ∆T, that means t A i − t B j ≤ ∆T, and generalizing it to a lagged coincidence if holds, where τ ≥ 0 is the time lag parameter. For the quantification of the strength of the interrelations between the two event series, two variants of the coincidence rate addressing B type events as either precursors or trigger have been introduced [1]: The precursor coincidence rate where Θ(x) is the Heaviside unit step function (equal to 0 for x ≤ 0 and 1 for x > 0) and 1 I (x) is the indicator function for the interval I (equal to 1 for x ∈ I, and 0 otherwise), which measures the fraction of A type events that are preceded at least by one B type event (i.e., multiple B type events within ∆T are counted once) and the trigger coincidence rate which measures the fraction of B type events that are followed by at least one A type event (i.e., multiple A type events within ∆T are counted once). The distinction between precursor and trigger coincidence rates allows the introduction of a notion of directionality while the parameter τ explicitly takes into account lagged interrelations between B type and A type events. Based on appropriate assumptions for example for the inter-event time distibutions [1] various statistical tests can be made to examine whether B type events are precursors to A type events for a risk enhancement test [4] or whether B type events are triggers for A type events trigger test [4]. Typical examples are the climate-related disasters as risk enhancement factor for armed conflicts in ethnically fractionalized countries [4] or the role of flood events as triggers of epidemic outbreaks [1]. Here, we made use of the CC.eca.es function of the CoinCalc package [3] for R [67] that implements ECA and compared against the (default) Poisson test, though it is well known that EQs appear in sequences in view of aftershocks (cf. CC.eca.es function allows two additional tests called "shuffle" and "surrogate". For the first case, the significance test is based on an empirical cumulative distribution function of the coincidence rates estimated by a large number of randomly shuffled time series having the same number of events like the original time series-thus numerically simulating a Poisson process-whereas for the second case, the significance test is based on a similar calculation for a large number of surrogate time series having the same waiting time distributions as the original data). The Poisson test has been selected because the magnitude range of the EQs (that define the A type events) considered here (see the first column of Tables 1-5) is such that it barely exceeds one magnitude unit while according to Båth law [11,33,[68][69][70] the difference in magnitude between the mainshock and its largest aftershock is approximately 1.2 magnitude units. Therefore, there are no aftershocks considered in our study and such EQs are expected to follow a Poisson process [71]. The function CC.eca.es provides, versus the time lag τ and the coincidence interval ∆T, the probability (p-value) to obtain by chance the observed precursor coincidence rate r p (∆T, τ) according the test model, the p-value to obtain by chance the observed trigger coincidence rate r t (∆T, τ) as well as the rates r p (∆T, τ) and r t (∆T, τ) themselves.
As an additional check, we also employed the "surrogate" test which could be in principle used for any empirical time-series and found p-values smaller than the values p J , p EM , p G 1 , p G 2 , and p G 3 , which will be reported below.

Results
In this section, we present the results of ECA for the interrelation between the strong (target) EQs and the order parameter fluctuation minima β min that have been found in the regional studies of Japan [59], Eastern Mediterranean [62], and in the global seismicity [60,61], see Figure 1. For these reasons, the time-series of target EQs will serve as event time-series A and the time-series of β min as event time-series B, whereas the unit of time for ∆T and τ is one day. Following a previous study [7], which employed ECA to estimate the statistical significance of Seismic Electric Signals (SES) [72,73] as EQ precursors in Greece and Japan, we focus our study in the cases where the precursory p-value is smaller than or equal to the trigger p-value. Moreover, we also focus on the cases where the precursor coincidence rate equals to unity, as in our studies there have been always observed β min before the target EQs. World map depicting the areas of Japan (green rectangle) and Eastern Mediterranean (yellow rectangle) at which regional studies based on the variability β have been made. The locations of the target EQs are also shown for both the global (red stars) and regional studies (green and yellow stars), see Tables 1, 2, and 4, respectively. The star corresponding to the M9 Tohoku EQ that occurred on 11 March 2011 in Japan is half red and half green since it is included both in global and regional studies.

Japan
Varotsos et al. [59] presented an analysis of seismicity of Japan based on the EQ catalog of JMA in natural time in which they studied the variability β W of the order parameter κ 1 [74] by means of EQ networks based on similar activity patterns (ENBOSAP). In this sense, they generalized the variability study of Sarlis et al. [30], which was made in the area N 46 o 25 o E 148 o 125 o , and by comparing the results between the two areas they managed to eliminate false alarms. In other words, the minima selected on the basis of criteria applied simultaneously in both areas are only the ones that precede the strongest EQs in the smaller area and there are no other β min that precede EQs of smaller magnitude there. Table 1 Figure 1) together with their precursory minima. When using the last two columns of Table 1 and apply ECA, we obtain the results presented by the color contours in Figure 2. Specifically, Figure 2A depicts r p (∆T, τ) as a function of the ∆T and τ, which are measured-as mentioned-in days, whereas in Figure 2B displays the corresponding p-values. According to the discussion of the previous paragraph, our region of interest is the one with r p (∆T, τ) = 1 that exhibits the lowest p-values. This is marked by the rectangle in Figures 2A,B, whereas in Figure 2C, we depict in expanded scale the p-values in this region. An inspection of the latter panel reveals that for τ = 22d and ∆T = 52d we have the minimum p-value, which turns out to be p J = 1.88 × 10 −7 .   Table 4 and Section 5 of [59]).

Eastern Mediterranean
The interconnection of ENBOSAP with the NTA of seismicity has been further developed in a later study [62] focusing in the area N 52 o 23 o E 50 o 5 o around Eastern Mediterranean which was based on the EQ catalog of the United States National Earthquake Information Center. In a fashion similar to that in Japan, the simultaneous application of criteria to two different areas selected on the basis of the network properties of ENBOSAP has led to characteristic minima β min that precede the strongest EQs of M ≥ 7.1 in the area N 52 o 28 o E 44 o 7 o , see the yellow dashed rectangle in Figure 1. Table 2 compiles the target EQs and their preceding β min . The results obtained by ECA for the last two columns of this table are shown in Figure 3. We observe that for τ = 19d ∆T = 164d a very small p-value is found, which turns out to be p EM = 5.80 × 10 −7 .   Table 2 and Section 5 of [62]) (cf. in the EQ catalog the reported magnitude is either mb or M w with mb being reported [75] for the smaller M(< 5.5) EQs when there is no authoritative M w available. For small magnitude EQs, however, mb practically equals M w , see Equation (31.II) of [76] and Equation (4) of [65]

Global Seismicity
The aforementioned regional studies of the variability of the order parameter κ 1 of seismicity have been also generalized to global seismicity by analyzing in natural time the Global Centroid Moment Tensor Catalog [77,78] of the period 1 January 1976 to 1 October 2014 [60,61]. Two different target EQ magnitude thresholds M thres have been used, i.e., M thres = 8.5 and M thres = 8.4, and for the latter threshold the method has been also applied when using the mid-scale time-series of seismicity [79] (cf. as discussed in detail in Appendix A in [60] the EQ prediction scheme becomes less efficient upon considering lower target magnitude thresholds). As such studies cannot inherently include two different areas, there also appear minima β min , which correspond to false alarms (FA), as they have been followed by EQs of magnitude smaller than M thres .

Focusing on EQs with M ≥ 8.5
During the study period, there exist six strong EQs with M ≥ 8.5 as shown in Table 3. They have been preceded within nine months by an equal number of β min (see the first six rows of the last column of Table 3), but the relevant study [60] has also led to three β min , which should be considered as FA in the sense discussed in the previous paragraph. In the present case, the six EQ occurrence dates constitute the A type event time-series in ECA, while the nine β min dates correspond to the B type event time-series. Figure 4 depicts the results of ECA and shows that when τ = 11d and ∆T = 255d a very small p-value arises that is p G 1 = 1.19 × 10 −5 . Table 3. The epicenters and dates of occurrence of all EQs with M ≥ 8.5 during the study period 1 January 1976 to 1 October 2014 as reported by the Global Centroid Moment Tensor Catalog together with the dates of their precursory variability minima. In the last three rows, we also include the dates of three more variability minima that could not be separated from the precursory ones (see Table 2 in [60]

Focusing on EQs with M ≥ 8.4
Sarlis et al. [60] and later Sarlis et al. [61] extended the methodology for the identification of β min to predict smaller EQs, i.e., the ones with M ≥ 8.4. This led to the results shown in Table 4. An inspection of this Table, shows that the seven EQs (in the first three columns) are preceded within nine months from eight β min , whereas there are eleven more FA β min . Using the seven EQ occurrence dates as A type event time-series and the nineteen β min dates as B type event time-series, we obtain the ECA results shown in Figure 5. Restricting ourselves to the cases when r t (∆T, τ) = 1, we find that the smallest p-value is p G 2 = 1.82 × 10 −4 that was achieved for the same values of τ and ∆T as before, i.e., for τ = 11d and ∆T = 255d.  It has been shown [61] that the EQ magnitude time-series of global seismicity can be separated into three different time-series, namely, the micro-, the mid-, and macro-scale time-series, which have drastically different multifractal properties. This is made possible by using empirical mode decomposition [80][81][82][83][84] and multifractal Detrended Fluctuation Analysis [85], the latter being a generalization of the well known Detrended Fluctuation Analysis (DFA) [86][87][88]. Out of the three component time-series, the mid-scale EQ magnitude time-series has been found to be the most useful one for EQ prediction purposes [61]. Table 5 summarizes the results found for the minima β min of the fluctuations of κ 1 estimated using only the mid-scale EQ magnitude time-series. We observe that now there are seven β min that precede within eight months the 7 EQs with M ≥ 8.4 and only 10 FA β min . Using the fourth and fifth columns of Table 5 as A type and B type event time-series in ECA, we find the results shown in Figure 6. In particular, one obtains a p-value p G 3 = 1.37 × 10 −5 for τ = 49d and ∆T = 186d.

Discussion
It has been shown by Varotsos et al. [89] that, for regional studies, the precursory β min occur upon the initiation of a series of precursory low frequency (≤0.1 Hz) electric signals, termed SES activity [90,91], as for example the one recorded by Uyeda et al. [27,92]. Briefly, once an SES activity has been recorded, a major EQ takes place a few weeks to 5 1 2 months later within a candidate epicentral area that can be determined on the basis of the so-called "selectivity" map of the recording station as well as by considering the ratio of the two horizontal electric field SES components (see, e.g., [90,91,93]). Considering the case of the strongest EQ in Japan which is tabulated in the last row of Table 1, the minimum β min of the fluctuations of the order parameter of seismicity before the M9 Tohoku EQ that occurred on 11 March 2011 was observed on 5 January 2011 being the deepest (β min ever observed in Japan, [30,59]) and is almost simultaneous with the initiation of anomalous magnetic field variations mainly in the z component. The latter were measured [94][95][96] at two geomagnetic stations (at Esashi and Mizusawa) lying approximately 130 km from the epicenter of the Tohoku EQ (cf. detectable magnetic field variations in the vertical component are expected to accompany SES activities before strong EQs, see, e.g., [97]). In view of this interrelation between SES activities and β min , we employed a range (170d) of around 5 1 2 months in the ECA of the regional studies depicted in Figures 2 and 3. The result that the lowest p-values, i.e., p J and p EM , were found for τ = 22d, ∆T = 52d (τ + ∆T < 6 months) and τ =19d, ∆T = 164d (τ + ∆T ≈6 months), respectively, strengthens the aforementioned simultaneity of the initiation of SES activities and β min . Additionally, the τ values found in both cases lie in the range 18 to 24d, which was revealed by ECA as the shortest period between the observation of an SES activity and the impending EQ occurrence in Greece [7] (as the cases studied in the latter work involve EQs of considerably lower target magnitudes, it is worthwhile in a future work to compare by means of ECA the resulting lead times for both precursors, i.e., β min and SES activities, and discuss the results against their simultaneity found in [89]). Finally, the fact that both p J (= 1.88 × 10 −7 ) and p EM (= 5.80 × 10 −7 ) are almost two orders of magnitude smaller than the p-values estimated [98] for the original variability study of Japan [30], shows that the use of appropriate EQ networks (i.e., the ENBOSAP suggested by [74]) may significantly optimize the predictive power of the minima of the fluctuations of the order parameter of seismicity in natural time.
We now turn to ECA for the β min of the global seismicity depicted in Figures 4-6. First, we have to mention that in these cases, where we do not employ ENBOSAP, we inevitably have β min that are followed by EQs of a magnitude smaller than the target M thres , which as mentioned correspond to FA. This fact explains why p G 1 , p G 2 , and p G 3 are almost two orders magnitude larger than p J and p EM of the regional cases (i.e., they are comparable to those found for the Japanese case by [98]). Second, in all three cases presented, p G 1 , p G 2 , and p G 3 are almost one order of magnitude smaller than the corresponding trigger p-values, thus ECA reveals the precursory character (risk enhancement factor, [4]) of β min for the strongest EQs. Third, the values of τ = 11d, ∆T = 255d and τ = 49d, ∆T = 186d sum up to a maximum lead time of nine or eight months, respectively, pointing to a different mechanism behind these minima compared to the regional ones. For example, the study of the situation before the aforementioned 2011 Tohoku EQ, as described in Tables 3-5, shows that the global β min was observed around 30 November to 3 December 2010, i.e., at dates markedly different than those of the regional β min and the simultaneous initiation of the anomalous variations in the vertical component of the geomagnetic field at Esashi and Mizusawa. Thus, it is highly likely that a preparation stage has already started even before the observation of the local β min at 5 January 2011. Interestingly, recent findings [8,37,99,100] based on the analysis of seismicity in natural time and non-extensive statistical mechanics [101][102][103] have revealed signs of critical behavior related with the preparation of a strong EQ around 22 December 2010 which lies between the date of the global and the regional β min . Fourth, p G 3 , that corresponds to the β min observed when analyzing the mid-scale global seismicity in natural time, is one order of magnitude smaller than p G 2 strengthening the importance of the mid-scale for EQ prediction (cf. mid-scale focuses, see Figure 7 of [61] to scales that comprise a number of seismic events that on average occur within a period of around a few months or so thus it compares favorably with the SES activities time lag and lead time). Finally, the fact that the three p-values (p G 1 , p G 2 , and p G 3 ) are by two orders of magnitude smaller than 5% indicates the statistical significance of the appearance of β min as precursors to strong EQs in a global scale.
It would be fruitful to investigate, in a future work, whether such minima appear in other systems of rich complexity (see, e.g., [104][105][106]) that exhibit avalanches or large sudden events, such as stock market crashes or turbulence (cf. the similarity between the latter two systems has been discussed in [107] while stock markets are known to exhibit EQ like phenomena [108]). Furthermore, we plan to study β min as precursors to strong avalanches observed in EQ models which can produce a large amount of data such as the Carlson-Langer model [109] or the Olami-Feder-Christensen model [110] to see if they show a similar drop in the variability. Concerning the latter model, it has been studied by means of the entropy in natural time [17] revealing minima of the entropy change ∆S under time reversal which have been also observed before strong EQs [36,37].

Conclusions
The statistical significance of the minima β min of the fluctuations of the order parameter of seismicity in natural time has been examined by event coincidence analysis as a possible precursor to strong earthquakes in both regional and global level. The ECA results show that β min are indeed statistically significant EQ precursors in both cases. The time lag τ and the coincidence interval ∆T found in the regional studies are compatible with the view that regional β min are simultaneous with the initiation of SES activities. At global scale, τ and ∆T point to the existence of a preparatory stage starting even earlier than the SES activity which is compatible with recent findings in the literature and may support the view that the whole solid Earth crust behaves as a single complex system.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript.

ECA
Event coincidence analysis ENBOSAP Earthquake networks based on similar activity patterns EQ Earthquake FA False alarm JMA Japan Meteorological Agency NTA Natural time analysis SES Seismic electric signals