Minima of the fluctuations of the order parameter of seismicity worldwide

TL;DR

This study identifies precursory minima in seismicity order parameter fluctuations, occurring about five and a half months before major earthquakes, and proposes a statistically significant prediction method based on these minima.

Contribution

It demonstrates that minima in seismicity fluctuations reliably precede large earthquakes and introduces a new prediction approach with high accuracy.

Findings

Minima occur ~5.5 months before earthquakes >8.4 magnitude

Minima are consistent across different major earthquakes

Proposed prediction method has 100% hit rate and 6.67% false alarm rate

Abstract

It has been shown [Phys. Rev. E 84, 022101 (2011); Chaos 22, 023123 (2012)] that earthquakes of magnitude $M$ greater or equal to 7 are globally correlated. Such correlations were identified by studying the variance $\kappa_1$ of natural time which has been proposed as an order parameter for seismicity. Here, we study the fluctuations of this order parameter using the Global Centroid Moment Tensor catalog for a magnitude threshold $M_{thres}=5.0$ and focus on its behavior before major earthquakes. Natural time analysis reveals that distinct minima of the fluctuations of the order parameter of seismicity appear {within} almost five and a half months{-on average-} before all major earthquakes of magnitude larger than 8.4. This phenomenon corroborates the finding [Proc. Natl. Acad. Sci. U.S.A. 110, 13734 (2013)] that similar minima of the seismicity order parameter fluctuations had preceded all major shallow earthquakes in Japan. Moreover, on the basis of these minima a statistically significant binary prediction method for earthquakes of magnitude larger than 8.4 with hit rate 100\% and false alarm rate 6.67\% is suggested.

Tables

Table 1. Table 1: All EQs with M ≥ 8.5 𝑀 8.5 M\geq 8.5 during the period 1 January 1976 to 1 October 2014 together with the outcome of the binary prediction scheme discussed in Section V .

EQ date	area	Lat. (oN)	Long. (oE)	$M_w$	outcome
2004-12-26	Sumatra-Andaman	3.30	95.78	9.0	TP
2005-03-28	Sumatra-Nias	2.09	97.11	8.6	TP
2007-09-12	Sumatra, Indonesia	-4.44	101.37	8.5	TP
2010-02-27	Chile	-35.85	-72.71	8.8	TP
2011-03-11	Tohoku, Japan	38.32	142.37	9.1	TP
2012-04-11	Indian Ocean	2.33	93.06	8.6	TP

Table 2. Table 2: The nine variability minima that were found during the whole study period together with the strongest EQ that followed min ⁡ ( β 160 ) subscript 𝛽 160 \min(\beta_{160}) within nine months ( Δ ​ t 100 Δ subscript 𝑡 100 \Delta t_{100} and Δ ​ t 160 Δ subscript 𝑡 160 \Delta t_{160} denote the corresponding lead times in months for W = 100 𝑊 100 W=100 and W ′ = 160 superscript 𝑊 ′ 160 W^{\prime}=160 , respectively).

$\min \left({\beta }_{100}\right)$ (date)	$\min \left({\beta }_{160}\right)$ (date)	$\frac{\min \left({\beta }_{160}\right)}{\min \left({\beta }_{100}\right)}$	EQ date	EQ area	$M$	$\Delta t_{100}$	$\Delta t_{160}$
0.2270(2004-04-05)	0.2431(2004-04-05)	1.071	2004-12-26	Sumatra-Andaman	9.0	8.8	8.8
0.1605(2005-01-28)	0.1702(2005-02-02)	1.060	2005-03-28	Sumatra-Nias	8.6	2.0	1.8
0.2770(2006-12-02)	0.2971(2006-12-20)	1.073	2007-09-12	Sumatra, Indonesia	8.5	9.5	8.9
0.2789(2008-08-25)	0.3049(2008-08-25)	1.088	2009-01-03	Papua, Indonesia	7.7	4.4	4.4
0.2317(2010-02-01)	0.2462(2010-02-16)	1.063	2010-02-27	Chile	8.8	0.9	0.4
0.2368(2010-11-29)	0.2638(2010-11-30)	1.114	2011-03-11	Tohoku, Japan	9.1	3.4	3.4
0.2849(2011-07-27)	0.3232(2011-08-04)	1.134	2012-04-11	Indian Ocean	8.6	8.6	8.4
0.2789(2012-05-20)	0.3053(2012-06-03)	1.095	2013-02-06	Solomon Islands	7.9	8.7	8.3
0.2612(2013-12-28)	0.2767(2014-01-13)	1.059	2014-04-01	Iquique, Chile	8.1	3.1	2.6

Equations

$$N(\ge M)\sim 10^{a-bM}$$

$$r(t,M)=\frac{1}{{\tau }_0{\left[1+t/c(M)\right]}^p},$$

$$p_k=Q_k/\sum _{n=1}^N{Q}_n$$

$${\kappa }_1=\sum _{k=1}^N{p}_k{\chi }_k^2-{\left(\sum _{k=1}^N{p}_k{\chi }_k\right)}^2\equiv ⟨{\chi }^2⟩-⟨\chi {⟩}^2.$$

$$Q_k\propto 10^{1.5M_k},$$

Full text

Keywords: seismicity, complex systems, order parameter, magnitude correlations

**Earthquake (EQ) time series manifest complexity in many aspects showing complex correlations between space, time and magnitude $M$. Recently, it has been found that there exist correlations between the successive EQs of magnitude 7.0 or larger in a global scale. This suggests that the whole of the solid Earth crust should be considered as a single complex system. Adopting this view, we study here these correlations focusing on the prediction of the time of occurrence of the strongest (magnitude class 9) EQs in the globe. We make use of a methodology recently applied to the Japanese region that is based on natural time analysis which allows the introduction of an order parameter for seismicity. The study of the fluctuations of this order parameter has revealed distinct minima one to three months before the occurrence of all major EQs around Japan. This kind of analysis is now generalized for the period 1 January 1976 to 1 October 2014 considering all EQs with magnitude greater than or equal to 5.0 in the Global Centroid Moment Tensor catalog. Considering a sliding natural time window of length equal to the number of EQs that on average occur within a few months, we show that similar minima can be also identified for the global seismicity. In particular, it is shown that all magnitude class 9 EQs, i.e, with magnitude greater than 8.4, have been preceded by such minima with an average lead time of five and a half months. This lead time almost coincides with the maximum lead time observed for the low frequency precursory electric signals termed Seismic Electric Signals (SES) activities, thus supporting the recent finding that the initiation of an SES activity marks the time when the system enters a critical stage. The properties of these minima are investigated and a simple prediction method based on these minima is suggested. Furthermore, the present results strengthen the view that the solid Earth crust behaves as a single complex system. **

I Introduction

Earthquakes (EQs) exhibit complex correlations in space, time and magnitude $M$ (e.g. see Refs. Sornette, 2000; Corral, 2004; Davidsen and Paczuski, 2005; Saichev and Sornette, 2006; Holliday et al., 2006; Eichner et al., 2007; Lennartz et al., 2008; Lippiello et al., 2009; Telesca and Lovallo, 2009; Telesca, 2010; Bottiglieri et al., 2010; Sarlis et al., 2010a; Lennartz et al., 2011; Sarlis, 2011; Sarlis and Christopoulos, 2012) which are inevitably interrelated with the fundamental power laws of seismicity. Principal among the latter is the Gutenberg-Richter law for the distribution of earthquake magnitudes which states that the (cumulative) number of earthquakes with magnitude greater than $M$ occurring in a specified area and time is given by

[TABLE]

where $b$ is a constant, which varies only slightly from region to region being generally in the range $0.8\leq b\leq 1.2$ (see Ref. Rundle et al., 2003 and references therein) and $a$ the logarithm of the earthquakes with magnitude greater than zeroShcherbakov et al. (2004) measuring the intensity of regional seismicityTurcotte et al. (2002) (note that Eq.(1) holds both regionally and globally). It is also well known that EQs are clustered in time and the Omori (power) law describes the temporal decay of aftershock activityUtsu (1961) (see also Ref. Shcherbakov et al., 2004):

[TABLE]

where $r(t,M)$ denotes the rate of occurrence of aftershocks with magnitudes greater than $M$ per day, $t$ is the time that has elapsed since the mainshock and $\tau_{0}$ and $c(M)$ are characteristic times. Note that $p$ is close to unity for large earthquakes (e.g., see Ref. Saichev and Sornette, 2007). There are also power-laws referring to the distributionSornette and Davy (1991) $\sim 1/L^{2}$ of fault lengths $L$, the fractal structure of fault networksDavy et al. (1990) as well as the universal law for the distribution of waiting times and seismic rates derivedBak et al. (2002) from the analysis of space-time windows. Considering seismicity as a (non-equilibrium) critical phenomenon, the observedTurcotte (1997) EQ power laws point to the proximity of the system to a critical point.Holliday et al. (2006) Adopting this view, an order parameter for seismicity labeled $\kappa_{1}$ has been introduced (Varotsos et al., 2005) in the frame of the analysis in a new time domain termed natural time $\chi$ (see Section II) which has been found to reveal novel dynamical features hidden in the time series of complex systems. (Varotsos et al., 2011a)

The study of the order parameter $\kappa_{1}$ in global seismicity has revealedSarlis (2011); Sarlis and Christopoulos (2012) that there exist correlations between successive EQ magnitudes when $M\geq 7.0$. These correlations have been identified by comparing the statistics of the experimentally observed $\kappa_{1}$ with those obtained when destroying by shuffling the order of occurrence of these EQs. The present paper focuses on the global fluctuations of $\kappa_{1}$ in view of the recent findingSarlis et al. (2013) that in Japan such fluctuations exhibit characteristic minima before the occurrence of EQs of magnitude $M\geq 7.6$. In short, the distribution of $\kappa_{1}$ for the number of EQs that on average occur within a few months has been studiedSarlis et al. (2013) by means of the ratio $\beta$ (calledSarlis et al. (2010b) variability, see Section II) of the standard deviation over the mean of the corresponding $\kappa_{1}$ values and $\beta$ minima have been identified before the stronger EQs.

Here, we analyze the Global Centroid Moment Tensor catalog, described in Section III, in order to examine the existence of similar $\beta$ minima. The results, which are presented in Section IV, show that such $\beta$ minima do exist in global seismicity and precede magnitude class 9 EQs. These results are further discussed in Section V and the corresponding conclusions are summarized in Section VI.

II Natural time analysis. Background

Natural time analysis has been shownAbe et al. (2005) to be optimal for enhancing the signals in time-frequency space. For a time series comprising $N$ events, e.g. EQs, the natural time for the occurrence of the ${k-\rm{th}}$ event of energy $Q_{k}$ is defined(Varotsos et al., 2001; Varotsos et al., 2002a, b) by $\chi_{k}=k/N$. We then study the evolution of the pair $(\chi_{k},p_{k})$ where

[TABLE]

is the normalized energy and construct the quantity $\kappa_{1}$ which is the variance of $\chi$ weighted by $p_{k}$

[TABLE]

In the case of EQs, the quantity $Q_{k}$ of Eq.(3) is the radiated seismic energy which is proportionalKanamori (1978) to the seismic scalar moment $M_{0}$ and can be estimated by means of the usual relationKanamori (1978)

[TABLE]

where $M_{k}$ is the moment magnitude of the ${k-\rm{th}}$ EQ.Varotsos et al. (2005); Sarlis et al. (2009, 2010b); Varotsos et al. (2011b, 2012a, 2012b, 2012c); Ramírez-Rojas and Flores-Márquez (2013); Varotsos et al. (2013); Flores-Márquez et al. (2014) The quantity $\kappa_{1}$ has been proposedVarotsos et al. (2005) as an order parameter for seismicity since it changes abruptly to zero upon the occurrence of a strong event and exhibitsVarotsos et al. (2005, 2006); Sarlis (2011); Sarlis et al. (2011); Sarlis and Christopoulos (2012) a behavior similar to that of the order parameter for various equilibrium and non-equilibrium critical systems.

As mentioned in Section I, the fluctuations of this order parameter versus conventional time can be studiedSarlis et al. (2010b); Varotsos et al. (2011b, 2012a, 2012b); Sarlis et al. (2013, 2015a) by means of its variability $\beta$ within excerpts of the EQ catalog ending just before a target EQ. The calculated $\beta$ value is assigned to this target EQ and hence to its (conventional) occurrence time. Specifically,Sarlis et al. (2015a) let us consider an EQ catalog $\{M_{p}\},p=1,2,\ldots,L$ and its excerpts, labeled by $j=0,1,\ldots,L-W$, comprising $W$ consecutive EQs, e.g., $\{M_{j+m}\},m=1,2,\ldots,W$. For such an excerpt, we form its subexcerpts consisting of the $n^{th}$ to $(n+5)^{th}$ EQs, $(n=1,2,\ldots W-5)$ (since at least $l=6$ EQs are neededVarotsos et al. (2005) for obtaining reliable $\kappa_{1}$) and compute $\kappa_{1}$ through Eq.(4) for each of them by assigning $\chi_{k}=k/6$ and $p_{k}=10^{1.5M_{j+n+k-1}}/\sum_{k^{\prime}=1}^{6}10^{1.5M_{j+n+k^{\prime}-1}}$, $k=1,2,\ldots 6$, to the $k^{th}$ member of the subexcerpt. This procedure leads to $(W-5)$ values of $\kappa_{1}$. We iterate the same process for new subexcerpts comprising $l=7$ members (giving rise to $(W-6)$ values of $\kappa_{1}$ ), 8 members …and finally $W$ members. Then, we compute the average $\mu(\kappa_{1})$ and the standard deviation $\sigma(\kappa_{1})$ of the thus-obtained $(W-5)+(W-6)+\ldots+1[=(W-4)(W-5)/2]$ $\kappa_{1}$ values. The quantity $\beta_{W}(j)\equiv\sigma(\kappa_{1})/\mu(\kappa_{1})$ is the variability of the order parameter $\kappa_{1}$ for this excerpt of length $W$ labeled by $j$. This $\beta_{W}(j)$ value is assigned to the $j+W+1$-th EQ in the catalog, the target EQ that took place at (conventional) time $t_{j+W+1}$. Hence, for the $\beta_{W}$ value of a target EQ only its past EQs are used in the calculation. The time evolution of the $\beta_{W}$ value is then pursued by sliding the excerpt $W$ through the EQ catalog, i.e., by considering the pairs $(t_{j+W+1},\beta_{W}(j)),j=0,1,\ldots,L-W$, for examples see Fig. 1.

III The Data Analyzed

The data analyzed come from the Global Centroid Moment Tensor (CMT) ProjectDziewoński et al. (1981); Ekström et al. (2012) that covers global seismicity since 1 January 1976. For EQs that took place before 2011, the 1976 to 2010 CMT catalog was used, whereas for EQs since 1 January 2011 to 1 October 2014 the monthly CMT catalogs have been employed (all these catalogs are publicly available from [http:

www.globalcmt.org/CMTfiles.html](http:%5C%5C%0Awww.globalcmt.org/CMTfiles.html)). The CMT catalog provides the scalar seismic moment $M_{0}$ for each listed event and thus this value was directly used as $Q_{k}$ in Eq.(3) for the calculation of the $\kappa_{1}$ values. In accordance with a previous study,Sarlis (2011) we considered all EQs of magnitude greater than or equal to $M=5.0$. This resulted in 38,006 EQs, whose epicenters are shown in Fig. 1 of Ref. Sarlis et al., 2015b, during the concerned period of 38 years and 9 months (1 January 1976 to 1 October 2014). Hence, we have a monthly rate of approximately 80 EQs. Table 1 summarizes all the magnitude class 9 EQs in CMT.

IV Results

The purpose of the present study is, as mentioned in Section I, to investigate whether there exist variability minima in the fluctuations of the order parameter $\kappa_{1}$ of global seismicity similar to those observedSarlis et al. (2013); Varotsos et al. (2014) in Japan. For this purpose, we examine $\beta_{W}$ when using as $W$ the number of events that on average occur within a few months. The reasoning behind this selection is that Seismic Electric Signals (SES), which are low frequency ($\leq 1$Hz) transient changes of the electric field of the Earth (the study of which has been motivated by solid state physics aspects on point defects, see Ref. Varotsos, 2008 and references therein) that have been found to precede EQs, exhibit a lead time within this range of a few months. The latter has been verified by numerous observations of SES activities in GreeceVarotsos and Lazaridou (1991); Varotsos et al. (1993, 1996, 2008); Varotsos et al. (2011a) which agree with later observations in Japan (Uyeda et al., 2000, 2002, 2009; Orihara et al., 2012; Xu et al., 2013) and China(Huang, 2011), as well as in California (see Bernardi et al. (1991); Fraser-Smith et al. (1990) where magnetic field variations similar to those associated with the SES activities in GreeceVarotsos et al. (2009) have been reported) and Mexico (see Ramírez-Rojas et al. (2011) and references therein), which also showed that within a few months after the SES activity observation a strong main shock takes place. Thus, this time scale of a few months provides an estimate of the time in which the system enters a critical stage before a strong EQ.

Along these lines, $\beta_{100}$, $\beta_{160}$, $\beta_{240}$, $\beta_{320}$ and $\beta_{400}$ have been plot in Fig.2 of Ref. Sarlis et al., 2015b versus the conventional time. It is observed that as $W$ decreases the corresponding $\beta_{W}$ fluctuates rapidly becoming more difficult to interpret. Indeed, when $W_{2}>W_{1}$, the quantity $\beta_{W_{2}}$ seems to smoothen the results obtained for $\beta_{W_{1}}$. This is so because the $\kappa_{1}$ values used for the calculation of various $\beta_{W_{1}}$s are included in the ensemble of $\kappa_{1}$ values from which $\beta_{W_{2}}$ is obtained (see Section II). Thus, one might say that when increasing the window length from $W_{1}$ to $W_{2}$ apparently a (non-linear) smoothing is performed. Along these lines, we start our analysis from Fig.2(e) of Ref. Sarlis et al., 2015b, we observe two clear precursory minima of $\beta_{400}$ with $\beta_{400}<0.26$: one before the 2004 Sumatra-Andaman $M$9 EQ and another one before the 2007 Sumatra, Indonesia $M$8.5 EQ (see Table 1) (cf. another minimum also appears but after the 2012 Indian Ocean $M$8.6 EQ, see the brown arrows in Fig.2(e) of Ref. Sarlis et al., 2015b which mark the four deepest minima of $\beta_{400}$). As the scale $W$ decreases, we observe, for example in Fig.2(c) of Ref. Sarlis et al., 2015b, that another more pronounced minimum of $\beta_{W}$ for $W\leq 240$ appears before the 2005 Sumatra-Nias $M$8.6 EQ (see Table 1), but the brown arrows in Fig.2(c) of Ref. Sarlis et al., 2015b show that there are still two deep minima of $\beta_{240}$ which are either postseismic or precursory to smaller EQs. For even smaller $W$, i.e., $W=100$ and $W=160$, two more minima become discernable before the 2010 Chile $M$8.8 EQ and the 2011 Tohoku $M$9.1 EQ. Thus, a comparison of Fig.2 of Ref. Sarlis et al., 2015b with Table 1 reveals that before all the six $M\geq 8.5$ EQs minima of the variability $\beta$ of the order parameter $\kappa_{1}$ of global seismicity are observed. Moreover, an inspection of the four deepest minima indicated by the arrows in Fig.2 Sarlis et al., 2015b shows that as $W$ decreases these deepest minima concentrate before the strongest EQs.

In order to clearly identify these minima, we focus on the smaller scales, i.e., $\beta_{100}$ and $\beta_{160}$, which, though more complex than those of larger scales, exhibit their deeper minima before all the $M\geq 8.5$ EQs. The situation then is similar to that in Japan where two scales $W=200$ and $W^{\prime}=300$ have been used.Sarlis et al. (2013) Interestingly, the seismic rate in the Japanese catalog wasSarlis et al. (2013) 148 EQs/month leading to the fact that $W=200$ and $W^{\prime}=300$ correspond to the average number of EQs in 1.35 and 2.03 months, respectively. The latter values do not differ much from 1.25 and 2.0 months calculated for $\beta_{100}$ and $\beta_{160}$ in the global CMT catalog (see Sec. III), thus supporting our hypothesis that the lead time of SES is a characteristic time scale for the solid Earth crust. The selection of two scales $W$ and $W^{\prime}$ is made because upon considering an almost simultaneous appearance of their local minima and examine if their ratio $\min(\beta_{W^{\prime}})/\min(\beta_{W})$ lies in the region around unitySarlis et al. (2013); Varotsos et al. (2014) a variability minimum can be unambiguously identified.Sarlis et al. (2013); Varotsos et al. (2014); Sarlis et al. (2015a) This reflects the premiseSarlis et al. (2013) that there exists a *single *(critical) process of consecutive EQs -charactererized by the ratio $\min(\beta_{W^{\prime}})/\min(\beta_{W})$- that leads to the variability minimum. After this variability minimum, the system tries to retain its equilibrium (cf. the curves in Fig.1 do not exhibit any significant trend although they consist of observables, i.e., $\beta_{W}$’s, based on a few months’ data for an almost 39 years period) with a preseismic process which eventually culminates with the occurrence of the strong EQ. Here, in order to identify a local minimum of either $\beta_{100}$ or $\beta_{160}$ we require that this should be a minimum for at least its 15 previous and 15 future values. Thus, such a local minimum is at least a local minimum for almost two weeks on average. Moreover, $\beta_{100}$ and $\beta_{160}$ should exhibit simultaneous minima. To secure this conditionVarotsos et al. (2014) we require at least 90% of the EQs that are included in the calculation of the local $\beta_{100}$ minimum are also included in the calculation of the local $\beta_{160}$ minimum. Once simultaneous local $\beta_{100}$ and $\beta_{160}$ minima are observed, we examine whether their ratio $\min(\beta_{160})/\min(\beta_{100})$ is close to unity as in the case of Japan.Sarlis et al. (2013); Varotsos et al. (2014) When this occurs, we identify the so called variability minimum. Selecting the range $(r_{1}=)1.05<\min(\beta_{160})/\min(\beta_{100})<1.15(=r_{2})$, which was determined on the basis of the minima observed before the $M\geq 8.5$ EQs, we obtain the nine variability minima shown in Table 2. For their determination, any value of $\beta_{100}$ below $(\beta_{0}=)$0.285 (which is the shallowest minimum $\beta_{100}$ value observed before an $M\geq 8.5$ EQ, see Table 2) was tested in order to examine whether it constitutes a local $\beta_{100}$ minimum. When such a local minimum was found, we investigated whether a value of $\beta_{160}$ that includes at least 90 EQs of the local $\beta_{100}$ minimum was a local $\beta_{160}$ minimum. If such a simultaneous local $\beta_{160}$ minimum was found, we examined if the ratio $\min(\beta_{160})/\min(\beta_{100})$ lied in the aforementioned range. If this was the case then we have a variability minimum and the corresponding $\min(\beta_{100})$ and $\min(\beta_{160})$ values are indicated with the red and green circles in Fig.1, respectively.

V Discussion

An inspection of Table 2 and Fig.1 shows that the variability minima in the global CMT catalog appear nine months at the most before the occurrence of magnitude class 9 EQs. Such a lead time is approximately three times larger than the corresponding lead time observedSarlis et al. (2013); Varotsos et al. (2014) in Japan for EQs of $M\geq 7.8$. This indicates that the preparation stage in global scale is presumably of larger duration than that in regional scale.

Moreover, a comparison of the present results with those obtainedSarlis et al. (2013); Varotsos et al. (2014) from the regional study of Japan before the strongest EQ, i.e., the 2011 Tohoku $M$9.1 EQ, leads to a compatible explanation: The variability minimum observed for this EQ in the global CMT catalog was on 29-30 November 2010. On the same date (see Table 1 of Ref. Sarlis et al., 2013 where local time is used) a minimum was also observed in the regional Japanese catalog that was associatedSarlis et al. (2013); Varotsos et al. (2014) with the 2010 Near Chichi-jima EQ that occurred 22 days later. The deepest variability minimum for the 27 year period studiedSarlis et al. (2013); Varotsos et al. (2014) in Japan was observed on 5 January 2011 being almost simultaneousSkordas and Sarlis (2014) with the anomalous geomagnetic diurnal variations recordedXu et al. (2013) at Esashi station at about 135km away from the Tohoku epicenter. Thus one may assumeSkordas and Sarlis (2014) that on the date of the deepest variability minimum in Japan an SES activity was recorded (whose magnetic field componentSarlis and Varotsos (2002); Varotsos et al. (2003) was recorded at Esashi). Then, the natural time analysis of seismicitySarlis et al. (2008); Varotsos et al. (2010) around the Tohoku epicenter leadsSkordas and Sarlis (2014) to the conclusion that the criticality conditionVarotsos et al. (2011c) $\kappa_{1}=0.07$ was fulfilled 31 hours before the occurrence of the 2011 Tohoku $M$9.1 EQ. These results strengthen the view that the (regional) variability minimum is observed simultaneously with the emission of an SES activity as found in Ref. Varotsos et al., 2013. If we now take into account the observationOrihara et al. (2014) that the level and temperature of confined groundwater in north eastern Honshu (Japan’s largest island) showed anomalous changes three months before the 2011 Tohoku $M$9.1 EQ (starting approximately on 10 and 15 December 2010, respectively, see their Fig.1) and such changes are precursory to large subduction EQs, we may have the view that the variability minimum in the global CMT catalog indicated the beginning of a preparatory process in a global scale, presumably related with long-range spatial correlations between EQ activity patterns (cf. such patternsTenenbaum et al. (2012) have been foundVarotsos et al. (2013) very important for the appearrence of regional variability minima), that was followed within 10-15 days by the aforementioned changes in confined groundwater which was later followed by the simultaneous appearance of the regional variability minimum and the SES emission, and finally by the 2011 Tohoku $M$9.1 EQ which is the largest EQ both regionally and globally. The fact that the minimum before Tohoku EQ was not found to be the deepest when using CMT, might be understood in the following context: In global scale, the Tohoku EQ magnitude does not differ much from the other $M\geq 8.5$ EQs that occurred in other regions of the world, but it does so in Japan where the second larger EQ was the 1994 East-Off Hokkaido $M$8.3 EQ and the relation between the value of $\min(\beta_{100})$ and the magnitude of the impending EQ is not very sharp (for example $\min(\beta_{100})$ before the 2014 Sumatra-Andaman $M$9.0 EQ is shallower than that before the 2005 Sumatra-Nias $M$8.6EQ and a similar picture holds also for Japan, see Table 1 of Ref. Sarlis et al., 2013, where the aforementioned 1994 East-Off Hokkaido $M$8.3 EQ exhibits the shallowest of the variability minima identified there). It is also noteworthy, that if the study of the global CMT catalog was made until the occurrence of the the 2011 Tohoku $M$9.1 EQ, then an inspection of Table 2 shows that one could have selected instead of 0.285 the value 0.277 (that preceded the 2007 Sumatra, Indonesia $M$8.5 EQ) to search for $\beta_{100}$ minima and in that case all EQs with $M\geq 8.5$ would have been preceded by variability minima (including the 2011 Tohoku $M$9.1 EQ) without any false alarm. By the same token, if such a selection was made until 1 October 2014, only the 2012 Indian Ocean $M$8.6 would have been missed leading to a hit rate of 80% with no false alarms.

We now turn to the statistical significance of the variability minima in the global CMT catalog reported in Table 2. As it was doneSarlis et al. (2016) in the case of the Japanese catalog, we follow the Receiver Operating Characteristics (ROC) method.Fawcett (2006) ROC is a technique to depict the quality of binary predictions. It is a plot of the hit rate (or True positive rate) versus the false alarm rate (or False positive rate), as a function of the total rate of alarms, which is tuned by a threshold in the predictor. The hit rate is the ratio of the cases for which the alarm was on and a significant event occurred over the total number of significant events. The false alarm rate is the ratio of the cases for which the alarm was on and no significant event occurred over the total number of non-significant events. Only if the hit rate exceeds the false alarm rate, a predictor is useful (for example, the ROC analysis has been recently usedTelesca et al. (2014) to disciminate between seismograms of tsunamigenic and non-tsunamigenic EQs). Random predictions generate on average equal hit and false alarm rate (thus, falling on the blue diagonal in Fig. 4 of Ref. Sarlis et al., 2015b), and the corresponding ROC curves exhibit fluctuations which depend on the positive $P$ cases (i.e., the number of significant events) and the negative $Q$ cases (i.e., the number of non-significant events) to be predicted. The statistical significance of an ROC curve dependsMason and Graham (2002) on the area $A$ under the curve in the ROC plane. It has been shownMason and Graham (2002) that $A=1-{U}/(PQ),$ where $U$ follows the Mann-Whitney U-statistics.Mann and Whitney (1947) Recently, a visualization scheme for the statistical significance of ROC curves has been proposed.Sarlis and Christopoulos (2014) It is based on k-ellipses which are the envelopes of the confidence ellipses -cf. a point lies outside a confidence ellipse with probability $\exp(-{\rm k}/2)$- obtained when using a random predictor and vary the prediction threshold. These k-ellipses cover the whole ROC plane and upon using their $A$ we can have a measureSarlis and Christopoulos (2014) of the probability $p$ to obtain by chance (i.e., using a random predictor) an ROC curve passing through each point of the ROC plane. In order to apply the ROC method, we divide the whole period covering 465 months into 51 nine-month periods (i.e., $P+Q=51$) out of which only 6 included significant events $(P=6)$, i.e., EQs with $M\geq 8.5$, see Table 2. Hence, the hit rate is 100%. On the other hand, the 3 minima which were followed by smaller EQs (see the fourth, eighth and ninth line of Table 1) may be considered as false alarms giving rise to a false alarm rate of $3/45\approx 6.67\%$. By using the FORTRAN code VISROC.f provided by Ref. Sarlis and Christopoulos, 2014 we obtain: (a) the ROC diagram of Fig. 4 of Ref. Sarlis et al., 2015b in which we depict by the red circle the operation point that corresponds to the results presented in Table 2 and (b) the probability $p$ to obtain this point by chance based on k-ellipses which results in $p=0.0047\%$. Interestingly, this $p$-value is of the order $10^{-5}$, thus being similar with that foundSarlis et al. (2016) for the variability minima in the regional studySarlis et al. (2013) of Japan. Moreover, when we vary any of the three parameters $(\beta_{0},r_{1},r_{2})$ that determine the operation point in order to obtain an ROC curve the resulting $p$-values estimated on the basis of the corresponding $A$-values lead to $p$-values ranging from $1.3\times 10^{-4}$ to $6.5\times 10^{-5}$, see the red, blue and cyan lines in Fig.4 of Ref. Sarlis et al., 2015b. Here, it is worthwhile to comment also on fact that the average value of the lower and the upper limit of the ratio of $\beta$ minima, $(r_{1}+r_{2})/2(=1.1)$, almost coincides with the value of the same ratio observed before the 2011 Tohoku $M$9.1 EQ which is the strongest EQ. This property -which is also observedSarlis et al. (2013); Varotsos et al. (2014) in the case of Japan- certainly points to the existence of a single (critical) process which is best monitored for the strongest EQ. It is also of interest to investigate whether such a prediction scheme could be extended to smaller target magnitude thresholds. Such an attempt is made in the Appendix.

In general, we have to stress the following: First, from the present analysis we may infer that a major EQ is likely approaching, but we cannot gain information on the likely epicenter location. This is clearly something to further consider in future work, probably by employing a procedure similar to that usedSarlis et al. (2015a); Huang (2015) to determine the future EQ epicentral location in the regional study of Japan. Similarly, the exact timing of the EQ cannot be predicted (just the average lead time of the obtained dynamical signatures is known from the present analysis, which varies from case to case). Second, the parameters $(\beta_{0},r_{1},r_{2})$ used for identifying precursory signatures have been optimized a posteriori based on a few major events (see also the Appendix). Detailed long-term observations will be required to further constrain these parameters (based on a larger number of strong events) and assess whether the observed signatures are sufficiently stable among different regions and EQ types.

Finally, we discuss the following two important points emerged in the present study. First, the scales $W$ that gave rise to the minima identified both globally and regionally correspond to the number of EQs that on average occur within one to two months. Such time scales have been first identifiedVarotsos and Lazaridou (1991); Varotsos et al. (1993) as lead times for SES activities. Second, the average lead time of the variability minima in the global CMT catalog is five and a half months which strikingly agrees with the maximum lead time observed for SES activities.Varotsos et al. (2011a) As for the regional study of Japan,Sarlis et al. (2013); Varotsos et al. (2014) the corresponding average lead time of the minima is 2.25 months that compares favorably with the average lead time of SES activities.Varotsos et al. (2011a)

VI Conclusions

The Global Moment Tensor Catalog in the period 1 January 1976 to 1 October 2014 was investigated here for the existence of precursory variability $\beta$ minima of the order parameter $\kappa_{1}$ of seismicity in natural time as it was previously foundSarlis et al. (2013) for the strong EQs in Japan. If we consider appropriate natural time scales $W$, corresponding to the number of events that on average occur within one to two months, we were able to identify minima of $\beta_{100}$ and $\beta_{160}$ that precede all the six magnitude class 9 EQs ($M\geq 8.5)$. Their lead time varies from 1 to 9 months before the main shock and using their properties we separated them from all other similar local minima of $\beta_{100}$ and $\beta_{160}$ except of 3 cases in 2008, 2012 and late 2013. The probability to obtain such a selection of precursory phenomena by chance is of the order of 10*-5*. The present results strengthen the view that the solid Earth crust behaves as a single complex system.

Appendix A Possible extension to smaller target magnitude thresholds

In the main text we presented a prediction scheme focused on the occurrence time of magnitude class 9 EQs, here we attempt an extension of this method to smaller target magnitude thresoholds. As a first step, we focus on the prediction of the occurrence time of the single $M=8.4$ EQ in global CMT which is the 2001 Peru $M$8.4 EQ, see Fig.5(b) of Ref. Sarlis et al., 2015b. This EQ exhibits a $\min(\beta_{100})$ value equal to 0.352 (see the inset of Fig.5(b) of Ref. Sarlis et al., 2015b) whereas the corresponding ratio ${\min(\beta_{160})}/{\min(\beta_{100})}\approx 1.09$ lies within the narrow range $[1.05,1.15]$ used for the prediction of $M\geq 8.5$ EQs (see Table 2) when using $\beta_{0}=0.285$. Increasing $\beta_{0}$ to 0.353 leads to 17 additional minima followed by smaller EQs (false alarms) giving rise to the uppermost point of the ROC diagram depicted in Fig.6 of Ref. Sarlis et al., 2015b (cf. a selection of $r_{1}=1.060$ and $r_{2}=1.135$, which is also compatible with the results of Table 2 but has not been selected based on the fact that one should also allow for a plausible experimental error of 0.01 in the ${\min(\beta_{160})}/{\min(\beta_{100})}$ ratio, leads in this case to only 12 false alarms). An inspection of Fig.6 of Ref. Sarlis et al., 2015b also reveals that such a scheme remains statistically significant with a $p$-value $1.4\times 10^{-4}$, comparable to the one presented in the main text.

When attempting to lower the target magnitude threshold below 8.4, e.g., to 8.3 (thus attempting to predict 12 EQs), a problem arises since the 1977 Sumbawa (Indonesia) $M$8.3 EQ occurred almost one and a half year after the beginning of the CMT catalog and the calculated $\beta_{160}$ values are available for only six months before this EQ, see Fig.5(a) of Ref. Sarlis et al., 2015b. This constrains $\min(\beta_{100})$ to $\min(\beta_{100})=$0.484 (see the dashed square in Fig.5(a) of Ref. Sarlis et al., 2015b) with a ratio ${\min(\beta_{160})}/{\min(\beta_{100})}\approx 1.04$, whereas a smaller value $\min(\beta_{100})$=0.459 observed early January 1997 cannot be used since there are no $\beta_{160}$ values to be paired with it (cf. there exist only 108 EQs in the CMT catalog for 1976). We mention that although the value of the ratio is close to the range $[1.05,1.15]$ used for the prediction of $M\geq 8.5$ EQs the value of $\min(\beta_{100})=$0.484 is very high (cf. the mean value and the standard deviation of $\beta_{100}$ for the whole period studied is $\mu=$0.593 and $\sigma=$0.220, respectively). The use of $(\beta_{0},r_{1},r_{2})=(0.484,1.04,1.15)$ in this case results in very frequent alarms amounting to $\tau$=47% of the total time (when assuming that the alarm is stopped upon the occurrence of an $M\geq 8.3$ EQ) and a hit rate of 83.3%(=10/12). However, by increasing only $r_{2}$ to 1.19 -thus selecting $(\beta_{0},r_{1},r_{2})=(0.353,1.05,1.19)$- one can obtain the same hit rate with $\tau=$32%, while when $(\beta_{0},r_{1},r_{2})=(0.353,1.05,1.32)$ a hit rate of 91.7%(=11/12) is achieved with $\tau=$37%. These results indicate that the prediction scheme becomes less efficient upon focusing on lower target magnitude thresholds. This is not incompatible with the view that EQs with $M<8.5$ may or may not lead the whole Earth crust system to criticality.

If we further decrease the target magnitude threshold down to 8.0, the number of EQs to be predicted within the almost 39 years period of our study increases to 30. We cannot any more employ a method similar to the one presented above in which the alarm was stopped upon the occurrence of an $M\geq 8.3$ EQ, because for example three days before the 2004 Sumatra-Andaman $M$9.0 EQ an $M$8.1 EQ took place, see Fig.1(a). For this reason, we considered the case that the alarm lasts exactly nine months after the observation of $\min(\beta_{160})$ indepent of whether a strong EQ takes place or not. In this case, the aforementioned selection of $(\beta_{0},r_{1},r_{2})=(0.353,1.05,1.19)$ leads to a hit rate of 66.7% with $\tau=38\%$. Alternatively, by decreasing $r_{1}$ to 1.04 -and hence using $(\beta_{0},r_{1},r_{2})=(0.353,1.04,1.19)$- we can obtain a hit rate of 73.3% with $\tau=41$% (see Fig.5 of Ref. Sarlis et al., 2015b). Interestingly, in the latter case, $(r_{1}+r_{2})/2$ lies very close to the value of the ratio ${\min(\beta_{160})}/{\min(\beta_{100})}$ observed before the 2011 Tohoku $M$9.1 EQ, see Table 2. An observation of Figs. 5(a) and 5(b) of Ref. Sarlis et al., 2015b reveals that during the period from 1976 to 2002, when the total number of EQs with $M\geq 8.0$ was only thirteen, the total alarm time amounted to $\tau=$22% with a hit rate 7/13$\approx 54$%. However, Fig.5(c) of Ref. Sarlis et al., 2015b, that covers the period 2002 to 2014 comprising of seventeen $M\geq 8.0$ EQs, reflects that the corresponding quantities are $\tau=$75% and 15/17$\approx 88$%. These numbers reveal that even for EQs with $M\geq 8.0$ such a scheme performs fairly well in both quiet and seismic active periods since the obtained hit rates are above 50% and larger than their corresponding $\tau$ values.