1 Introduction

Natural time analysis, introduced in the beginning of 2000s (Varotsos et al. 2001b, 2002a, b, 2003a, b), uncovers unique dynamic features hidden behind the time series of complex systems. It has found applications in a large variety of diverse fields, and the relevant results have been compiled in a monograph (Varotsos et al. 2011c). In general, in the investigation of time series associated with complex systems, the following two types of complexity measures are considered (Varotsos et al. 2011a): scale-specific and scale-independent measures. Upon employing natural time analysis, it has been also found (Varotsos et al. 2011a, c) that complexity measures that employ fluctuations on fixed timescales seem to complement those measures that take into account fluctuations on different timescales. For example, when analyzing electrocardiograms in natural time aiming at the distinction between sudden cardiac death individuals (SD) and healthy (H) individuals (e.g., Varotsos et al. 2007), we make use of the entropy fluctuations and the aforementioned complementarity holds in the following sense (Varotsos et al. (2011a); see also pp. 410–413 of Varotsos et al. (2011c)): If in the frame of the latter complexity measures (i.e., different natural timescales) an ambiguity emerges in the distinction between SD and H, the former complexity measures (i.e., fixed natural timescales) give a clear answer. Such a complexity measure for fixed natural timescales rather than the order parameter fluctuations has been applied by Varotsos et al. (2011a, 2013) to the complex time series of seismicity to which we now turn.

In the case of seismicity, in a time series comprising N earthquakes, the natural time χ k  = k/N serves as an index for the occurrence of the kth earthquake. This index together with the energy Q k released during the kth earthquake of magnitude M k , i.e., the pair (χ k ,Q k ), is studied in natural time analysis. Alternatively, one studies the pair (χ k ,p k ), where

$$p_{k} = \frac{{Q_{k} }}{{\sum\nolimits_{n = 1}^{N} {Q_{n} } }}$$
(1)

stands for the normalized energy released during the kth earthquake. The variance of χ weighted for p k , designated by κ 1, is given by (Varotsos et al. 2001b, 2003a, b, 2005, 2011a, c)

$$\kappa_{1} = \sum\limits_{k = 1}^{N} {p_{k} } \left( {\chi_{k} } \right)^{2} - \left( {\sum\limits_{k = 1}^{N} {p_{k} } \chi_{k} } \right)^{2}$$
(2)

where Q k , and hence p k , for earthquakes is estimated through the usual relation (Kanamori 1978).

$$Q_k \propto 10^{1.5M_k} .$$
(3)

Earthquakes exhibit complex correlations in time, space and magnitude (e.g., Huang 2008, 2011; Lennartz et al. 2008, 2011; Rundle et al. 2012, 2016). Taking the view that earthquakes are (non-equilibrium) critical phenomena (e.g., Holliday et al. 2006; Varotsos et al. 2011b), we employ the analysis in natural time χ, because in this analysis an order parameter for seismicity has been proposed. In particular, it has been explained (Varotsos et al. 2005) (see also pp. 249–253 of Varotsos et al. 2011c) that the quantity κ 1 given by Eq. (2)—or the normalized power spectrum in natural time Π(ω) as defined by Varotsos et al. (2001b, 2002a) for natural angular frequency ω → 0—can be considered as an order parameter for seismicity since its value changes abruptly when a mainshock (the new phase) occurs, and in addition, the statistical properties of its fluctuations resemble those in other non-equilibrium and equilibrium critical systems. Note that at least 6 earthquakes are needed for obtaining reliable κ 1 (Varotsos et al. 2005).

Let us consider a sliding natural time window of fixed length comprising W consecutive events. After sliding, event by event, through the earthquake catalog, the computed κ 1 values enable the calculation of their average value µ(κ 1) and their standard deviation σ(κ 1). We then determine the variability β of κ 1, i.e., the quantity β W (Sarlis et al. 2010)

$$\beta_{W} = \frac{{\sigma \left( {\kappa_{1} } \right)}}{{\mu \left( {\kappa_{1} } \right)}}$$
(4)

that corresponds to this natural time window of length W. To compute the time evolution of β W , we apply the following procedure (Varotsos et al. 2011a, c; Sarlis et al. 2013, 2015): First, take an excerpt comprising W(≥100) successive earthquakes from the seismic catalog. We call it excerpt W. Second, since at least 6 earthquakes are needed for calculating reliable κ 1 (Varotsos et al. 2005), as mentioned above, we form a window of length 6 (consisting of the 1st to the 6th earthquake in the excerpt W) and compute κ 1 for this window. We perform the same calculation of successively sliding this window through the whole excerpt W. Then, we iterate the same process for windows with length 7, 8 and so on up to W. [Alternatively, one may use windows with length 6, 7, 8 and so on up to l, where l is markedly smaller than W, e.g., l ≈ 40 (Varotsos et al. 2011a, c, 2013)]. We then calculate the average value µ(κ 1) and the standard deviation σ(κ 1) of the ensemble of κ 1 is thus obtained. The quantity β W  ≡ σ(κ 1)/µ(κ 1) defined in Eq. (4) for this excerpt of length W is assigned to the (W + 1)th earthquake in the catalog, the target earthquake. (Hence, for the β W value of a target earthquake only its past earthquakes are used in the calculation.) The time evolution of the β W value can then be pursued by sliding the excerpt W through the earthquake catalog, and the corresponding minimum value (for at least W values before and W values after) is labeled β W,min.

Varotsos et al. (2011a, 2013) and Sarlis et al. (2013) found that the results become exciting upon using a crucial scale, i.e., when these W consecutive events extend to a time period comparable to the lead time of the seismic electric signals (SES) activities, i.e., of the order of a few months (Varotsos and Lazaridou 1991; Varotsos et al. 1993, 2011c), which are series of low frequency (≤0.1 Hz) electric signals that precede major earthquakes (e.g., Varotsos et al. 2008 and references therein). This originates from an effect that has been recently uncovered and explained in the next section. In particular, in Sect. 2 we indicate that, if geoelectrical data are available, an SES activity initiates when a minimum β W,min of the fluctuations of the order parameter of seismicity is observed. An unprecedented minimum β W,min has been identified in the beginning of January 2011, i.e., almost two months before the Tohoku earthquake occurrence in Japan, as it will be explained in Sect. 3. A discussion follows in Sect. 4, and finally in Sect. 5, we summarize our conclusions.

2 An SES activity initiates when the fluctuations of the order parameter of seismicity minimize

This phenomenon has been reported by Varotsos et al. (2013) who used as an example the pronounced SES activity with innumerable signals observed by Uyeda et al. (2002, 2009) almost two months before the onset of the volcanic-seismic swarm activity in 2000 in the Izu Island region, Japan. This swarm was then characterized by Japan Meteorological Agency (JMA) as being the largest earthquake swarm ever recorded (Japan Meteorological Agency 2000). In particular, the date of the initiation of SES activity was reported to be on April 26, 2000, and the first major earthquake of magnitude 6.5 occurred on July 1, 2000, close to Niijima Island (Fig. 1).

Fig. 1
figure 1

Epicenters (stars) of all mainshocks with magnitude 7.6 or larger (depth < 400 km) within the area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\) (black rectangle) studied by Sarlis et al. (2013) since January 1, 1984, until the M W9 Tohoku EQ on March 11, 2011. The smaller area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\) studied by Varotsos et al. (2013) is also shown with yellow rectangle. The small yellow square indicates the location of the Niijima Island where the SES activity of the volcanic-seismic swarm activity in 2000 in the Izu Island region has been recorded (Uyeda et al. 2002, 2009)

Full size image

Varotsos et al. (2013) analyzed in natural time the series of the earthquakes reported during this period by the JMA seismic catalog and considered all EQs within the area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\) (the yellow rectangle in Fig. 1) which almost covers the whole Japanese region. The seismic moment M 0 that is proportional to the energy released during an EQ (which is the quantity Q k of the kth event used in natural time analysis) was obtained from the magnitude M JMA in the JMA catalog by using the approximate formulae of Tanaka et al. (2004) that interconnect M JMA with M w. They set a threshold M JMA > 3.4 to assure data completeness, leading to 52,718 EQs from 1967 to the time of Tohoku EQ, thus having on the average of the order of 102 EQs per month. Since we are interested in timescales comparable to that of the lead time of the SES activity, as mentioned, Varotsos et al. (2013) employed natural time window lengths W = 100, 200, 300 and 400 earthquakes. Furthermore, for each earthquake e i in the seismic catalog, they calculated the κ 1 values upon using for example the previous 6 to 40 consecutive earthquakes. They found that β W minimum was observed on the following date(s): April 23, 2000, for W = 100, 26 April for W = 200, 21–23 April for W = 300 and 23–24 April for W = 400 (see their Fig. 1b, which is here reproduced in Fig. 2c for the sake of comparison). These results show that the date of the β W minimum of seismicity, after considering that the experimental error in its determination is around a few days, approximately coincides with the date (≈April 26, 2000) of the initiation of the SES activity. In other words, Varotsos et al. (2013) concluded that an SES initiates almost when the fluctuations of the order parameter of seismicity minimize.

Fig. 2
figure 2

Variability β W of κ 1 of seismicity versus the conventional time from March 1, 2000, to July 1, 2000, which is the date of occurrence of the M6.5 EQ close to Niijima Island during the volcanic-seismic swarm activity in the Izu Island region in 2000. a In the smaller area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\), b in the larger area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\). For the sake of comparison, we also plot in (c) the β W values of κ 1 of seismicity in the smaller area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\) for the values W = 100, 200, 300 and 400 as computed in Varotsos et al. (2013) and shown in their Fig. 1(b)

Full size image

In a later study (Varotsos et al. 2014), a set of criteria has been developed to distinguish truly precursory minima β W,min of the order parameter of seismicity from non-precursory ones. One of these criteria demands that a truly precursory minimum β W,min should exhibit spatial invariance in the following sense: It should appear almost on the same date (within the aforementioned experimental error of course) when changing the boundaries of the area studied. Let us now investigate here whether this criterion is obeyed in the aforementioned case studied by Varotsos et al. (2013). Along these lines, Varotsos et al. (2014) investigated also the area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\) which is somewhat larger than the area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\) studied by Varotsos et al. (2013). These two areas are shown with black rectangle (\({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\)) and yellow rectangle (\({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\)), respectively, in Fig. 1 along with a small yellow square indicating the location of the Niijima Island where the pronounced SES activity of the volcanic-seismic swarm activity in 2000 in the Izu Island region was recorded by Uyeda et al. (2002, 2009). In addition, Varotsos et al. (2014) increased the magnitude threshold from M JMA = 3.4 to M JMA = 3.5 to assure a more safe data completeness. Setting this threshold, and considering the period from January 1, 1984, to March 11, 2011 (the date of M W9 Tohoku earthquake), we are left with 47,204 EQs and 41,277 earthquakes in the concerned period of about 326 months in the larger (black rectangle) and smaller (yellow rectangle) area, respectively. Thus, we have on the average ∼145 and ∼125 earthquakes per month for the larger and smaller area, respectively. For the sake of brevity in the calculation of β W values for both areas, we use the values W = 200 and W = 300, which would cover a period of around a few months before each target earthquake. Furthermore, for the sake of comparison between the two areas, we also investigate the case of W = 250 since this value in the larger area roughly corresponds to the case W = 200 in the smaller area. Figure 2 provides an overview of these results. In particular, the following quantities are plotted versus the conventional time during the period from March 1, 2000, to July 1, 2000 (the date of the occurrence of the M6.5 earthquake close to Niijima Island): In Fig. 2a, we show the quantities β 200 and β 300 (in red and blue, respectively) for the smaller area, and in Fig. 2b, we depict β 200, β 250 and β 300 (in red, green and blue, respectively) for the larger area. In these figures, we then observe that after the last days of March the variability β W exhibits a gradual decrease and a minimum β W,min appears on a date around the date of the initiation of the SES activity. In particular, in Fig. 2b, the relevant curve (green) for β 250 in the larger area minimizes on April 25, 2000, which is approximately the date of the initiation of the SES activity reported by Uyeda et al. (2002, 2009), lying also very close to the date (21 April)—if the experimental error is considered—at which in the smaller area the β 200 curve (red) in Fig. 2a minimizes, thus verifying the main conclusion of Varotsos et al. (2013).

The two different geophysical observables under discussion, i.e., SES activity and minimum of fluctuations of the order parameter of seismicity, were shown by Varotsos et al. (2013) to be also linked in space as follows: Using a narrow 5° × 5° spatial window sliding through the (21° × 21°) area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\), we investigated the earthquake events that would occur in an almost two months period. We searched for the small 5° × 5° areas that result in a clearly observable (i.e., the deepest) β minimum the date of which differs by no more than a few days from the date April 26, 2000, of the β minimum observed by Varotsos et al. (2013) in the 21° × 21° area for W = 200. This occurs for the seismicity in the small areas \({\text{N}}_{33}^{38} {\text{E}}_{135}^{140}\) and \({\text{N}}_{31}^{36} {\text{E}}_{135}^{140}\) that include the Izu Island region which became active a few months later. In other words, these two 5° × 5° areas surrounding the forthcoming earthquakes that occurred later on July 1 and 30, 2000, result in the deeper β minima that are observed almost simultaneously with the initiation of the SES activity on April 26, 2000. This inspired the determination of the epicentral area of a forthcoming earthquake on the basis of the β minimum of a sliding spatial window discussed later in Sect. 4.

3 An unprecedented minimum of the fluctuations of the order parameter of seismicity almost two months before the Tohoku earthquake occurrence

For the vast majority of major earthquakes in Japan, the simultaneous appearance of the minima of the fluctuations of the order parameter of seismicity with the initiation of SES activities discussed in Sect. 2 could not be directly verified due to the lack of geoelectrical data. In view of this lack of data, an investigation was made by Sarlis et al. (2013) that was focused on the question whether minima β W,min of the fluctuations β W of the order parameter of seismicity are observed before major earthquakes from January 1, 1984, to March 11, 2011 (the date of the M W9 Tohoku earthquake occurrence), in Japanese area. Actually such minima were identified to appear a few months before large earthquakes; see Fig. 3 (cf. in this figure, beyond the mean and the standard deviation of the β W values in panels (e) and (f), we also observe β W maxima just after the occurrence of large earthquakes. This is understood in the context that, as explained by Varotsos et al. (2005), the κ 1 value abruptly decreases to zero upon the occurrence of a large earthquake, thus leading to a large κ 1 fluctuation and hence to a β W maximum). Among these minima, the minimum observed on January 5, 2011, i.e., almost two months before the M W9 Tohoku earthquake, was the lowest. No minimum of β W with this depth was observed at any other time of the whole period studied. This result was in agreement with the results obtained by Varotsos et al. (2011a, b) in California and Greece. In particular, the California seismic catalog was investigated over a 25-year period (1979–2003) and the Greek seismic catalog over a 10-year period (1999–2008). In both cases, the fluctuations of the order parameter of seismicity exhibited a global minimum before the strongest earthquake occurred during the period investigated.

Fig. 3
figure 3

Variability β W (left scale) versus the conventional time given in three consecutive 10-year periods in (a), (c) and (e) for the smaller area and in (b), (d) and (f) for the larger area. In addition, all M JMA ≥ 6.5 earthquakes (in black, M JMA in the right scale) are plotted. In panels (e) and (f), the descriptive statistics of β W in each case are shown by the error bars that correspond to the mean ± one standard deviation interval. The maxima of β W are related either to the appearance of a bimodal feature in the distribution of κ 1; e.g., see Sarlis et al. (2010); Varotsos et al. (2012a) or due to a unimodal distribution maximized at κ 1 = 0 (e.g., see Fig. 2 of Varotsos et al. 2012c). The latter behavior is usually observed when a strong earthquake is included in the corresponding catalog excerpt

Full size image

Sarlis et al. (2013) also concluded that it appears that there are two kinds of β W minima, namely precursory and non-precursory to large earthquakes. In particular, Sarlis et al. (2013) analyzed in natural time the seismicity of Japan from January 1, 1984, to March 11, 2011, using sliding natural time window of lengths W comprising the number of events that would occur in a few months within the area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\). They set a threshold M JMA = 3.5 to assure a safe data completeness which leads on the average to a number of the order of ∼102 earthquakes per month. They chose the values W = 200, 300 and 400, which would cover a period of a few months before each target earthquake and obtained the following results: Distinct minima of β were found one month to three months before the occurrence of all other Japanese major mainshocks (M JMA ≥ 7.6, depth < 400 km) during 1984–2011. With less certitude, nine other minima of β may have also been precursory to large earthquakes. The minima of β which seem to be precursory to sizable earthquakes commonly show the β 300,min /β 200,min ratio close to unity in the range of 0.95–1.08, whereas the other minima show the ratio outside this range.

In other words, Sarlis et al. (2013) found that the fluctuations of the order parameter of seismicity exhibited 15 distinct minima—deeper than a certain threshold—one to a few months before the occurrence of large earthquakes that occurred in Japan from January 1, 1984, to March 11, 2011. Six (out of 15) of these minima were followed by all the shallow mainshocks of magnitude 7.6 or larger during the whole period studied. Sarlis et al. (2016), in a more recent study, demonstrated the statistical significance of the findings by Sarlis et al. (2013) since they showed that the probability to achieve the latter result by chance is of the order of 10−5. This conclusion was strengthened by employing also the receiver operating characteristic technique.

It is important to note that the minima β W,min become almost invisible for greater W that would correspond to time intervals appreciably larger than a few months, e.g., W = 2000 and 3000 (see Fig. 2C of Sarlis et al. 2013), thus emphasizing that it is crucial to use a timescale comparable to that of the lead time of SES activities.

Let us now investigate whether the deepest β W,min before the M W9 Tohoku earthquake obeys spatial invariance in a similar fashion as discussed in Sect. 2, i.e., by considering both the larger area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\) and the smaller area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\). (Recall that Sarlis et al. (2013) considered only the larger area in their study.) We find that the results of this investigation, which are depicted in Fig. 4 and summarized in Table 1, remain unchanged upon varying the area studied, thus verifying the criterion of spatial invariance.

Fig. 4
figure 4

Excerpt of Fig. 3 plotted in expanded timescale before the 2011 M W9.0 Tohoku earthquake for the smaller area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\) (a) and the larger area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\) (b), respectively

Full size image
Table 1 Values of β W,min, in the larger area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{148}\) and in the smaller area \({\text{N}}_{25}^{46} {\text{E}}_{125}^{146}\). The dates in parentheses correspond to the occurrence time of the last earthquake that led to the appearance of the corresponding β W,min with the format YYYY-MM-DD
Full size table

4 Discussion

The following two important comments are now in order:

First, it is remarkable that the minimum β W,min before the M W9 Tohoku earthquake was observed during a period (from January 4, 2011, to January 14, 2011, see the blue curves in Fig. 4) in which unique anomalous geomagnetic field variations appeared on the z component, at two stations lying at epicentral distances of around 130 km (Xu et al. 2013; Han et al. 2015, 2016), which are characteristic of an accompanying strong SES activity (Varotsos et al. 2001a; Sarlis and Varotsos 2002; Varotsos et al. 2003c).

Second, generalizing a procedure employed by Varotsos et al. (2013) to show that the two phenomena, SES activity and minimum of the fluctuations of the order parameter of seismicity, are also linked in space, we have recently ascertained (Sarlis et al. 2015) that a spatiotemporal study of the minima β W,min observed before the six M JMA ≥ 7.6 mainshocks during 1984–2011 in Japan leads to an estimation of their epicentral areas (see also Huang 2015, on this finding). The extent of the usefulness of this result for the case of the M W9 Tohoku EQ in 2011 and in particular for the determination of its occurrence time has been recently investigated by Varotsos et al. (2017). They reviewed, in the light of the recent advances, the procedure to identify the occurrence time which in general is achieved by computing the κ 1 values of the seismicity occurring in the candidate epicentral area after the initiation of the SES activity. Along these lines, for the case of the M W9 Tohoku EQ, Varotsos et al. (2017) computed the κ 1 values of seismicity in the epicentral area estimated in Sarlis et al. (2015) by starting this computation from January 5, 2011, since it is the date of the appearance of the deepest minimum β W,min (Sarlis et al. 2013), which was found to be simultaneous (Skordas and Sarlis 2014)—as mentioned above—with the anomalous geomagnetic field variations and hence with strong SES activity. The results of this computation showed that the condition κ 1 = 0.070 (which signals that the mainshock is going to occur within a few days or so) was fulfilled almost one day before the Tohoku EQ (Varotsos et al. 2017).

Thus, in other words, Varotsos et al. (2013) reported for the first time the physical interconnection between SES activity and the minimum of the fluctuations of the order parameter of seismicity, i.e., that these two phenomena are linked not only in time (since the SES activity before the Izu Island volcanic-seismic activity in 2000 initiates when β W,min of seismicity appears) but also in space. Here, this physical interconnection is further strengthened since it is shown to be valid also for another case, i.e., the M W9 Tohoku EQ in the following sense: By combining the results obtained from the natural time analysis of the seismicity data (Skordas and Sarlis 2014; Sarlis et al. 2013, 2015) and the anomalous geomagnetic field variations reported by Xu et al. (2013) and Han et al. (2015, 2016), we find that β W,min of seismicity is almost simultaneous with the latter anomalous variations that are characteristic of a strong SES activity; in addition, these two phenomena are also linked in space since a spatiotemporal study of β W,min of seismicity by Sarlis et al. (2015) led to a region that lies within a few hundred kilometers from the M W9 Tohoku EQ epicenter.

5 Conclusions

When we adopt natural time analysis and use a sliding natural time window of length W so that these W consecutive events correspond to a time period comparable to the lead time of SES activities, i.e., of the order of a few months, the following precursory phenomena are unveiled:

An unprecedented minimum β W,min of the fluctuations of the order parameter of seismicity appears approximately on January 5, 2011, i.e., almost two months before the M W9 Tohoku earthquake that occurred on March 11, 2011. This becomes almost invisible if we use considerably larger values of W.

The aforementioned unprecedented minimum β W,min occurs almost simultaneously with the appearance of unique anomalous variations of the geomagnetic field that are characteristic of a strong SES activity since they are observed mainly on the z component.

The above simultaneous appearance is strikingly similar with the one identified by Varotsos et al. (2013) between the minimum β W,min of the fluctuations of the order parameter of seismicity and the pronounced SES activity observed by Uyeda et al. (2002, 2009) almost two months before the volcanic-seismic activity in 2000 in the Izu Island region, Japan.