Natural Time Analysis: The Area under the Receiver Operating Characteristic Curve of the Order Parameter Fluctuations Minima Preceding Major Earthquakes.

It has been reported that major earthquakes are preceded by Seismic Electric Signals (SES). Observations show that in the natural time analysis of an earthquake (EQ) catalog, an SES activity starts when the fluctuations of the order parameter of seismicity exhibit a minimum. Fifteen distinct minima-observed simultaneously at two different natural time scales and deeper than a certain threshold-are found on analyzing the seismicity of Japan from 1 January 1984 to 11 March 2011 (the time of the M9 Tohoku EQ occurrence) 1 to 3 months before large EQs. Six (out of 15) of these minima preceded all shallow EQs of magnitude 7.6 or larger, while nine are followed by smaller EQs. The latter false positives can be excluded by a proper procedure (J. Geophys. Res. Space Physics 2014, 119, 9192-9206) that considers aspects of EQ networks based on similar activity patterns. These results are studied here by means of the receiver operating characteristics (ROC) technique by focusing on the area under the ROC curve (AUC). If this area, which is currently considered an effective way to summarize the overall diagnostic accuracy of a test, has the value 1, it corresponds to a perfectly accurate test. Here, we find that the AUC is around 0.95 which is evaluated as outstanding.

1. Introduction

Earthquakes (EQs) exhibit complex correlations in time, space and magnitude, e.g., [1,2,3,4,5,6,7,8,9,10]. The observed EQ scaling laws point to [11,12,13] the existence of phenomena closely associated with the proximity of the system to a critical point (the mainshock is the new phase [13,14]). In this frame, the order parameter of seismicity is the quantity by means of which one can identify the approach of the dynamical system to the critical point. Such a parameter has been introduced [15] upon analyzing the seismicity, in a new time domain, termed natural time [13,14,16,17,18] (see below) which has found useful applications in diverse fields [13,19,20]. This analysis unveils hidden properties in time series of complex systems [13] and has recently been also used by Turcotte and coworkers as basis of a new methodology to estimate the current seismic risk level [21,22,23,24,25,26,27,28].

In a time series comprising N EQs, the natural time for the occurrence of the k-th EQ of energy is defined as , i.e., we ignore the time intervals between consecutive events, but preserve their order as well their energy . In natural time analysis, the evolution of the pair is studied, where is the normalized energy and is estimated by means of the relation [29] , where stands for the EQ magnitude. It has been argued [15] that the variance of natural time weighted for , namely may serve as an order parameter of seismicity. The entropy S in natural time is defined [18] by where the brackets denote averages with respect to the distribution , i.e., . Upon considering time reversal , i.e., , the value S changes to a value :

The physical meaning of the entropy change in natural time under time reversal has been discussed in References [13,19].

The study of the fluctuations of this order parameter of seismicity reveals challenging results. To compute fluctuations, we follow the procedure described in detail in References [30,31] by using a sliding natural time window of constant length, i.e., consisting of a number W of EQs that would occur on average within the crucial scale [32] of a few months or so, which is the lead time of Seismic Electric Signals (SES) activities. These are series of low frequency transient changes of the electric field of the Earth [33,34] that are detected before major EQs (both in Japan [35] and Greece [36,37,38]). We then compute the average value and the standard deviation of the ensemble of obtained. The quantity is termed [13] variability of . The time evolution of the value can then be pursued by sliding the natural time 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 in the EQ catalog. The corresponding minimum value is labeled . Please note that of Equation (4) is reminiscent [39] of the square root of the Ginzburg criterion introduced within the frame of the mean field theories of the Ginzburg-Landau type, e.g., see p. 175 of Goldenfeld [40], discussed also by Holliday et al. [12].

The following results were revealed: Using the Japan Meteorological Agency (JMA) seismic catalog [30] and considering all EQs of magnitude from 1 January 1984 to 11 March 2011 (the time of the M9 Tohoku EQ) within the area 25–46 E, 125–148 E (Figure 1), fifteen distinct minima—observed simultaneously (see, e.g., Appendix A of Reference [31]) at and with a ratio in the range 0.95 to 1.08 and —of the fluctuations of the order parameter of seismicity were found 1 to around 3 months before large EQs. All shallow EQs of magnitude 7.6 or larger during this 27 year period (see Figure 1) were preceded by six (out of 15) of these minima (see Figure 2), the spatiotemporal variations of which also reveal an estimate of the candidate epicentral area [41]. Remarkably, among the minima, the deepest minimum was observed around 5 January 2011 (which is almost two weeks after the minimization on 22 December 2010 of the entropy change under time reversal [42]), i.e., around two months before the M9 Tohoku EQ.

Figure 1

(color online) Epicenters (red stars) of all six shallow EQs with M ≥ 7.6 within the area N E since 1 January 1984 until the M9 Tohoku EQ. The smaller green stars indicate the epicenters of all M ≥ 3.5 EQs.

Figure 2

(color online) Variability versus the conventional time depicted in consecutive 6 year periods in the upper graphs of the panels (a–e), respectively. The six minima preceded the M ≥ 7.6 EQs are marked by red circles at the curve while the nine minima followed by smaller EQs by green squares again at the curve; the values of are also written in red and green, respectively (see, e.g., Tables 1 and 2 of Reference [30]). No data are plotted in (e) after M9 Tohoku EQ. The horizontal red line corresponds to the shallowest minimum that preceded a M ≥ 7.6 EQ and the EQs are marked as black arrows whose magnitude can be read in the right scale. In addition, below the variability graph in each panel, we depict the time dependent seismicity rate of the temporal epidemic-type aftershock sequence (ETAS) model [47,48,49] according to Equation (1) of Ogata et al. [50] together with the seismicity (black vertical lines ending at circles whose magnitude can be read in the right scale). The ETAS model parameters are the same as those presented in Figure 2a in Reference [50].

Using Monte Carlo calculations, Sarlis et al. [43] showed that the probability to achieve the above results by chance is of the order of (we shall return to this point later). This conclusion was also recently strengthened by Christopoulos et al. [44] using the event coincidence analysis (ECA) [45]. It is the scope of the present paper to investigate the diagnostic accuracy of the precursory minima before major EQs in Japan during 1984–2011 by employing the area under the receiver operating characteristic (ROC) [46] curve (AUC), which is a very recent technique for judging the quality of binary predictions.

2. Receiver Operating Characteristics Technique

The performance of a diagnostic test in the case of a binary predictor can be evaluated using the measures of sensitivity [51] (where sensitivity = (True positives)/(True positives + False negatives)) and specificity (where specificity = (True negatives)/(True negatives + False positives)) (e.g., Reference [52]). This is achieved by a ROC curve that includes all the possible decision thresholds from a diagnostic test results. Simply defined, a ROC curve is a plot of the sensitivity versus the quantity ‘1-specificity’, i.e., the False positive rate=(False positives)/(True negatives + False positives), of a diagnostic test. Hence, a ROC diagram depicts the sensitivity or hit rate (or True positive rate) versus the False positive rate (or false alarm rate) thus showing the trade-off between hits and false alarms [46]. The AUC is an effective way to summarize the overall diagnostic accuracy of the test. It takes a value of 1 for a perfectly accurate test. An AUC of 0.5 suggests no discrimination (i.e., ROC curve falls on the diagonal), 0.7 to 0.8 is considered acceptable, 0.8 to 0.9 is considered excellent, and more than 0.9 is considered outstanding [51]; see also Chapter 5, pp. 160–164 of Reference [53]. Moreover, as shown by Mason and Graham [54] AUC can be used to estimate the statistical significance of the prediction scheme. Recently, a method was proposed [55] that can estimate the AUC—and hence the statistical significance—corresponding to an operating point in the ROC plane using the so-called k-ellipses which are constructed on the basis of confidence ellipses. According to this method, when only one operating point of a prediction method is known, then the AUC of the corresponding ROC curve can be approximated by the AUC of the k-ellipse that passes through this point, e.g., see Figure 3 in Reference [55].

By means of the latter method, the AUC corresponding to the precursory variations of the Earth’s electric and magnetic field that precede EQs in Greece in the 1980s, see, e.g., Dologlou [56], (i.e., the two points depicted in Figure 4 of Reference [57]) and EQs in Japan in the 2000s, see, e.g., Han et al. [58], (i.e., the two points in Figure 8a,b of Reference [57]) can be estimated. They are 0.625 and 0.658 for the two operating points in Figure 4 of Reference [57] and 0.869 and 0.943 for the the two operating points shown in Figure 8a,b of Reference [57], respectively. These four AUC exceed the one shown in Figure 11 of Reference [59] which has been calculated for the relationship between magnetic field pulses and EQs in California.

3. Data Analyzed

Here, as in References [30,41], the JMA seismic catalog was used and the energy of each EQ was obtained from after converting [60] to the moment magnitude defined by Hanks and Kanamori [61]. We consider all EQs—by setting a threshold = 3.5 in order to assure data completeness (see Reference [30]; see also Figure 6 in Nanjo et al. [62])—from 1 January 1984 to the M9 Tohoku EQ occurrence on 11 March 2011, within the area N E, which covers the whole Japanese region (see Figure 1). Since 47,204 EQs occurred in this period of about 326 months, we have on average EQs per month, thus we choose the natural time window lengths W = 200 and 300 for the calculation of that would correspond on average to a few months period in accordance with the SES activities observations, as already mentioned.

4. Results

In Figure 2, we plot the variability for W=200 (red) and W=300 (blue) (left scale) along with all ( in the right scale) versus the conventional time in four consecutive 6 year periods (a), (b), (c), (d), respectively, and one period (e) of almost 3 years. The six minima preceeding the EQs (see Table 1 of Reference [30]) are marked with red circles, while the nine followed by the EQs (see Table 2 of Reference [30]) with green squares. Other minima below the threshold (i.e., ) that have not been marked by red circles or green squares, did not obey any of the following criteria (selected for reasons explained in detail by Sarlis et al. [30] and Varotsos et al. [31]): and should appear simultaneously (see, e.g., Appendix A of Reference [31]) with a ratio in the range 0.95 to 1.08.

Figure 3 depicts the ROC graph, i.e., the plot of True positive rate versus the False positive rate, as a function of the total rate of alarms with an alarm period of 3 months, which is tuned by a threshold in the predictor. In particular, we vary the value of the shallowest from 0.295 down to 0.157 and we obtain the ROC depicted by the red solid line which has AUC=0.951. Similarly, upon increasing the value of the lower threshold of the ratio / from 0.95 to 1.08, we obtain the green dotted ROC with AUC = 0.965. Finally, the dashed blue ROC with AUC = 0.943 corresponds to the case that we decrease the maximum value of the ratio / from 1.08 down to 0.95. In the same figure, we also depict the k-ellipses that correspond [55] to the 95% and 99% confidence intervals. We observe that in all these three cases the observed ROCs exceed the 99% confidence interval exhibiting statistical significance (see also below).

Figure 3

ROC curves resulting from the study of the minima of the fluctuations of the order parameter of seismicity in Japan. For all the three ROCs, the corresponding AUCs have a mean value 0.95 reflecting an outstanding discrimination (see p. 162 of Reference [53]). The cyan and magenta lines are the k-ellipses that correspond [55] to the 95% and 99% confidence intervals. They indicate how far away from the diagonal the ROC curve of a random predictor may scatter with probability 1/20 or 1/100.

5. Discussion

AUC is currently considered [51] an effective way to summarize the overall accuracy of a diagnostic test, as already mentioned in Section 2. It is characterized outstanding when AUC is more than 0.9 [51], which is of course the present case of since the AUC computed here is around 0.95. This result corroborates with the fact that the statistical significance of the precursory nature of has been shown by two other independent procedures, in particular Monte Carlo calculations in Reference [43] and ECA in Reference [44].

The above result (AUC ≈ 0.95) can be also used to to estimate the overall diagnostic accuracy of SES and the associated Earth’s magnetic field variations, in view of the following experimental fact emerged from measurements of independent research groups: Observations show that there exists simultaneous appearance of the minima of the fluctuations of the order parameter of seismicity with the initiation of SES activities, before major EQs in Japan and Greece. This was [63] the first time in the literature that before major EQs anomalous changes are found to appear simultaneously (as well as they are also linked in space) in two independent data sets of different geophysical observables (geoelectrical measurements and seismicity). Focusing on measurements in Japan for example, let us consider the volcanic seismic swarm activity in 2000 in the Izu Island region, which was then characterized by JMA as being the largest EQ swarm ever recorded [64], and the M9 Tohoku EQ on 11 March 2011 which is the strongest EQ ever recorded in Japan. Specifically, a straightforward analysis of the JMA seismic catalog in natural time, by employing a sliding natural time window comprising the number of events that would occur in a few months, it was observed [63] that the fluctuations of the order parameter of seismicity exhibit a clearly detectable minimum at the time of the initiation of the pronounced SES activity identified by Uyeda et al. [35,65] almost two months before the onset of the volcanic-seismic swarm activity in 2000 in the Izu Island region, Japan. Concerning the before the M9 Tohoku EQ, this was observed [30] on 5 January 2011 being the deepest minimum during the period from 1 January 1984 until the M9 Tohoku EQ occurrence, as mentioned. This date almost coincides with the detection of anomalous magnetic field variations on the z component during the period from 4 to 14 January 2011 at two measuring sites lying at epicentral distances of around 130 km [66,67,68] pointing to the initiation of an SES activity (as observed for SES activities in Greece, see, e.g., pp. 8–9 of Reference [37], in which the associated magnetic field variations are clearly detectable at epicentral distances up to around 150–200 km for EQs of magnitude 6.0 or larger). This occurred almost two weeks after the observation [42] on 22 December 2010 of the minimum of the entropy change of seismicity under time reversal.

The probability to obtain such a minimum by chance was shown to be approximately 3%, thus it is statistically significant. Such a minimum is of precursory nature signaling that a large EQ is impending according to the conclusions deduced from the natural time analysis of the Olami-Feder-Christensen (OFC) model for EQs [69], which is probably [70] the most studied non-conservative self-organized criticality (SOC) model. In particular, it has been shown that exhibits a clear minimum [13] (or maximum if we define [71] instead of ) before a large avalanche in the OFC model, which corresponds to a large EQ. The minimum observed on 22 December 2010 has been recently shown to be of profound importance in identifying the occurrence time of the M9 Tohoku EQ [72,73]. It was also found [74] (see p. 321) that it corresponds to the first stage of the physical model for the SES generation according to which an excess stress disturbance starts gradually increasing until reaching the critical stress. This minimum was accompanied by an abrupt increase of the order parameter of seismicity [75] to which we now turn in view of the interesting property that it exhibits a natural time scale dependence that has a functional form suggested by Penrose and coworkers [76] consistent with the phase transition kinetics of Lifshitz and Slyozov [77]. Beyond the six (out of 15) minima preceding all shallow EQs of magnitude 7.6 or larger, an inspection of Figure 2 reveals that a increase appears upon the occurrence of a major EQ [75]. In particular, a careful study of Figure 2 shows that there exist six prominent fluctuations increases of (higher than 1.2) before the M9 Tohoku EQ. These are accompanied by major EQs and specifically the increase observed in the beginning of 1993 after the M7.5 EQ on 15 January 1993 at 42.92 N 144.35 E and the increases of associated with all shallow EQs of magnitude 7.6 or larger that occurred from 1 January 1984 until the M9 Tohoku EQ occurrence, which are the following: The M7.8 EQ on 12 July 1993, the M8.2 EQ on 4 October 1994, the M7.6 EQ on 28 December 1994, the M8.0 EQ in 2003, and the M7.8 EQ in 2010. The two larger are in 2003 and in 2010: First, the large increase of in 2003 appears upon the occurrence of the M8 EQ on 26 September 2003 (attaining its maximum on 28 September 2003, see Figure 2b of Reference [75]). Second, the increase of on 22 December 2010 is observed upon the occurrence on the same day of the M7.8 near Chichi-jima EQ. The latter, i.e., the one in 2010, appears almost simultaneously with the minimum the entropy change of seismicity under time reversal, which occurs also on 22 December 2010 as it was found upon applying the procedure of Reference [42] to the entire Japanese region seismic data.

In order to assure the simultaneity of the minima with the initiation of SES activities as well as their statistical significance, attention is drawn to a few misunderstandings explained below:

First, only the minima of the order parameter fluctuations defined in Equation (4), but not the minima of (as claimed in Reference [78]), are of precursory nature [30,43,44].

Second, as mentioned in the introduction, the probability to achieve the results concerning of Reference [30] by chance is of the order . This value was obtained both by Monte Carlo calculations as well as with the ROC technique. For example, in the latter technique, the calculation was made as follows: Sarlis et al. [30] found 15 distinct minima during the 27 year period from 1 January 1984 to 11 March 2011 and six of these minima preceded (by 1–3 months) all shallow EQs of magnitude 7.6 or larger. In general, the ROC curves exhibit fluctuations which depend on the positive P cases (the number of significant events, i.e., the 6 EQs of magnitude 7.6 or larger in the present case) and the negative Q cases (the number of non-significant events, i.e., M < 7.6) to be predicted. Thus, by dividing the 27 year period by almost 3 months, we have 109 three-month periods (i.e., P + Q = 109) out of which only six included significant events (P = 6). Hence, all six significant events were successfully predicted, thus the hit rate is 100%. The  nine minima that were followed within three months by smaller EQs may be considered false alarms, thus the false alarm rate (False positive rate) is 9/103 ≈ 8.74%. A straightforward calculation by using the FORTRAN code VISROC.f, provided in Reference [55], gives that the probability p to obtain this operation point (8.74%, 100%) by chance based on k-ellipses is of the order of . In other words, this result () demonstrates that it is incorrect to claim [78] that since there exists an abundance of false positives (i.e., 6 true positives while 9 false positives in the present case) the method of identifying that are simultaneous with SES activities is unusable without, however, paying attention to its statistical significance (cf. This can be also shown analytically by applying the formulation given in the Appendix of Reference [79]). Such a line of thinking (abundance of false positives) is in contrast with the proper use of ROC curves as discussed in detail in Section 6 of Reference [46], see Equation (1) there, according to which the selection of the operating point should be selected on the basis of the trade of the losses in case of false negative and the cost of a false positive.

Third, the aforementioned nine false alarms (false positives) could have been avoided by following the procedure developed in Reference [31] in which the analysis of seismicity of Japan and the computation of was made in two different areas NE and NE (instead of the single area NE in Reference [30]) which are compatible with the results obtained in Reference [8] by means of EQ networks based on similar activity patterns. Comparing the results between these two areas all false alarms have been excluded since they did not simultaneously obey the criteria applied to both areas. This way, as true positives were selected only the minima preceding the strongest EQs in the smaller area, while the remaining minima preceding EQs of smaller magnitude were excluded.

Fourth, concerning the detection of true SES activities, we clarify that electrical disturbances that may look to be similar to SES activities, but not obeying the criteria developed in References [17,18,36] to distinguish them from noise, are not classified as SES activities.

6. Conclusions

Upon analysing the seismicity of Japan during an almost 27 year period, i.e., from 1 January 1984 to 11 March 2011 (the time of the M9 Tohoku EQ occurrence), we identify the minima of the fluctuations of the order parameter of seismicity (which are in general accompanied by simultaneous SES activities). Focusing on the area under the receiver operating characteristic curve, which is currently considered an effective way to summarize the overall diagnostic accuracy of a test, we found a value close to AUC = 0.95. This shows that the precursory nature of these minima is outstanding, after taking into account that AUC takes the value of 1 for a perfectly accurate diagnostic test.