An Application of the Coherent Noise Model for the Prediction of Aftershock Magnitude Time Series

1Faculty of Engineering, Environment and Computing, Coventry University, Priory Street, Coventry CV1 5FB, UK 2Solid Earth Physics Institute, Department of Physics, School of Science, National and Kapodistrian University of Athens, Panepistimiopolis, Zografos, 157 84 Athens, Greece 3Section of Solid State Physics, Department of Physics, School of Science, National and Kapodistrian University of Athens, Panepistimiopolis, Zografos, 157 84 Athens, Greece


Introduction
The prediction of the magnitudes and occurrence times of aftershocks is of crucial importance for restricting the losses caused by strong earthquakes (EQs, hereafter) because buildings already damaged by the mainshock may collapse upon the occurrence of a strong aftershock.Recently, an algorithm has been suggested [1] on the basis of the coherent noise model [2][3][4] and natural time [5][6][7] that might be useful for the determination of both magnitudes and occurrence times of aftershocks.It is the main scope of the present paper to investigate the applicability of such an algorithm to aftershock time series in Southern California.This area has been selected in view of the publication [8] of an accurate waveform relocated EQ catalogue for Southern California from 1981 to June 2011 [9] that exhibits tighter spatial clustering of seismicity than the routinely generated catalogue does.
The coherent noise model [2][3][4] is a model that shows reorganization events (avalanches) whose size distribution follows a power law over many decades and displays aftershock events.These events have been shown [3,10,11] to exhibit a behavior similar to that of the Omori-Utsu law [12]; see also [13,14], for real EQ aftershocks.Moreover, it has been recently shown [15] that it is compatible with the unified scaling law [16] of waiting times between EQs.

Data Analyzed.
As mentioned in the introduction, we consider the waveform relocated catalogue [8] for Southern California from 1981 to June 2011 [9].In order to filter out distant poorly constrained events, we use the geographic polygon shown in Figure 1 of Hauksson et al. [8] that covers the Southern California Seismic Network (SCSN) reporting area for local events.The  ≥ 6.0 events that occurred within this area are shown in Table 1.As seen in this table, there are 4 stronger EQs with  ≥ 6.7 which are the 1992 Landers, the 2010 El Mayor-Cucapah, the 1999 Hector Mine, and the 1994 Northridge EQs (see Figure 1), and 4 smaller EQs which are the two Superstition Hills large events [53] in 1987 and the 1992 Joshua Tree and Big Bear events (see Figure 2).The latter two EQs are related [54,55] to the Landers EQ which is the strongest one in California for the last sixty years (e.g., see Table 1 of [56]).Following Shcherbakov et al. [13], we consider as aftershocks all reported events that occurred during a period  within a region centered at the epicenter of the each strong EQ (see the last column of Table 1).The linear dimension of each region scales with 10 0.5  starting from the region of 1.1 ∘ × 1.1 ∘ selected by Shcherbakov et al. [13] for Landers EQ (cf. the scaling of the aftershock zone with mainshock magnitude   was firstly introduced by Utsu and Seki [57]).This scaling has been employed since it allows the determination of the "aftershock" spatial window upon the occurrence of the mainshock and does not make use of any a posteriori information concerning the aftershocks.Hence, the spatial range of the space-time window usually called Omori regime-which is identified [58] by examining the best fits of the Omori-Utsu law to the aftershock data-may differ from this window.The period  examined was one year for the four stronger mainshocks while it varies for the four smaller EQs (see Results and Discussion).

Coherent Noise
Model and the Algorithm.The events, to be called avalanches, in the coherent noise model result from the following procedure [2,3].Consider a system of   agents; for example, the points of contact in a subterranean fault, for each agent  we associate a threshold   ,  = 1, 2, . . .,   that represents the amount of stress that an agent withstands before it breaks.Without loss of generality [3],   may come from the uniform distribution in the interval 0 ≤  < 1.The dynamics of the model consists of two steps, a "stress" step, which is more important and sufficient to produce large avalanches, and an "aging" step.During the "stress" step, we select a random number (or "stress" level)  from some probability distribution function  stress () and replace all   that are smaller than  with new values resulting from the uniform distribution in the interval 0 ≤  < 1.The number of agents whose thresholds are renewed is the size  of the avalanche.During the "aging" step, a fixed fraction  of agents is selected at random and their thresholds are also replaced with new thresholds obtained again from the uniform distribution in the interval 0 ≤  < 1.If we assume that   → ∞, the thresholds   are represented by a threshold distribution function  thres () which initially ( = 0) is considered uniform in the interval 0 ≤  < 1; that is,  (0)  thres () = 1 (see Figure 3).The size  1 of the first avalanche thres () =  1 which represents the "mass" of the agents that had thresholds smaller than  1 .After the subsequent "aging" step the threshold distribution becomes  (1)  thres ().Upon repeating the two steps for the second time using  2 we can obtain the size  2 of the second avalanche and the new threshold distribution  (2)  thres () and so on.In the case of the coherent noise model for   → ∞, the threshold distribution  ()  thres () after the th avalanche is a piecewise constant function having [1] the following general form: where   is the number of steps present in the threshold distribution function after the th avalanche (e.g., see Figure 3),  () 0 = 0,  ()  , and  ()  are positive real parameters, and Θ() is the Heaviside (unit) step function; that is, Θ() = 0 if  < 0 and Θ() = 1 if  ≥ 0. The analysis of numerical results has  shown (e.g., see Figure 5 of [1]) that the integer number   may serve as a predictor for the size  +1 of the next avalanche.It has been also found [1] that   as a function of  (and hence versus natural time) exhibits the ubiquitous 1/ noise [59][60][61][62][63] (together with a logarithmic trend, e.g., see Eqs. (31) and ( 32) of [1]).
The number of steps   changes when a "stress" level  is applied: (a) eliminating at least two smaller "stress" levels previously applied (cf. the case  = 3 in Figure 3) or (b) being smaller than the smaller nonzero  ()   (cf. the case  = 2 in Figure 3).When all   < 1,   coincides with the cardinality |  | of the sets   of successive extrema formed when    studying the time series of the random stresses   .The sets   of successive extrema are defined [25] as follows:  0 equals the empty set 0. Each subsequent   is obtained by the procedure described below for  times: select a random number   from a given probability density function () and compare it with all the members of  −1 .In order to construct the set   , we discard from the set  −1 all its members that are smaller than   and furthermore include   .Thus,   ̸ = 0 for all  > 0 and   is a finite set of real numbers whose members are always larger than or equal to   .Moreover, max[  ] ≥ max[ −1 ]. thres () after the th avalanche in the coherent noise model with   → ∞: for  = 0 (solid red),  = 1 (dashed green) exhibiting a single step at  = 0.2,  = 2 (dashed blue) with two steps at  = 0.1, 0.2, and  = 3 (dotted magenta) with a single step at  = 0.3.The stress levels  1 = 0.2,  2 = 0.1, and  3 = 0.3 have been applied during the first, second, and third "stress" step, respectively.The value of  = 0.1 has been used.
The increase of the cardinality |  | of these sets is at the most 1, but its decrease may be as large as |  | − 1, thus exhibiting a characteristic asymmetry as  increases.
The aforementioned predictability of the coherent noise model on the basis of   could be also seen as follows: when a relatively large random stress   is applied, the coherent noise model almost "forgets" the effect of previously applied smaller stresses acquiring relatively high threshold density for a relatively wide threshold range starting from  = 0 (cf. the case  = 3 for  < 0.3 in Figure 3).Since the archetypal physical system of the coherent noise model assumed (e.g., see the last few lines of the first page of [2]) that the agents correspond to the points of contact in a subterranean fault, it has been suggested [1] that the coherent noise model predictability can be exploited in the case of EQ aftershocks.In the latter case, one studies the successive extrema formed in the time series of the aftershock magnitudes   (  ≥  thres , where  thres is a threshold magnitude) as a function of their order of occurrence , that is, natural time, by assuming   =   /10.Setting  = 0 at the occurrence time of the mainshock, as shown in the upper left rectangular box of Figure 4, the corresponding sets M  of the successive extrema (as defined above) in the aftershock magnitude time series are constructed and the cardinality   ≡ |M  | (e.g., see Figure 4) plays a role similar to that of   in the case of the coherent noise model.A simple computer program that calculates   is presented in Appendix.
Figure 5(a) depicts the aftershock magnitude time series of the Landers EQ (see Table 1) versus conventional time together with the time series of   .The argument of the latter time series has been intentionally shifted by unity in order to  compare the value of  −1 , that is, the number of successive extrema before the occurrence of the th aftershock, with the magnitude   of the th aftershock.We observe a characteristic pattern which reflects the already mentioned asymmetry when   is considered as time series with respect to the natural number , but no clear correlation between   and  −1 can be seen.For this reason, we plot in Figure 5(b) the excerpt corresponding to the first three and a half hours after the mainshock occurrence.During this period, the strongest aftershock of Landers EQ (see the sixth line of Table 1), the Big Bear aftershock [54], took place.We may observe that most of the time strong aftershocks occur when  −1 is below 5.In order to further clarify what happens, it is most appropriate to depict the results in natural time.This is done in Figure 5(c) where we observe that the Big Bear aftershock occurs when  −1 was equal to 4. Thus, we observe that for Landers EQ the suggestion that  −1 might be a useful predictor could be considered plausible.The subject of the next section is to examine statistically whether  −1 can be used for the prediction of the magnitude of the forthcoming aftershock, in the sense that this magnitude exceeds a threshold.Here, we note that the coherent noise model method might be able to predict the magnitude of the next aftershock; however it cannot say when it will occur.This makes impossible the use of other verification methods as the Molchan diagram [64][65][66].

Results and Discussion
We first examine whether   as it results from the aftershock time series exhibits 1/ noise.To this end, we use the Detrended Fluctuation Analysis (DFA) [67][68][69] which is a standard robust method suitable for detecting long-range power-law correlations and has been applied to diverse fields where scale-invariant behavior emerges, for example, from DNA [70][71][72] and heart dynamics [68,73] to meteorology [74,75], atmospheric physics [76,77], geophysics [28,40,78,79], and economics [80][81][82][83][84][85].The major advantage of DFA is the systematic elimination of polynomial trends of different order [86][87][88][89] (and hence of the much slower aforementioned logarithmic trend in   ).In DFA, a fluctuation measure  DFA () is constructed, based on the residuals of a piecewise polynomial fitting to the profile time series, and studied versus the scale  used for the polynomial detrending (e.g., for more details see Section 1.4.3 of [7]).If  DFA () ∝   the time series exhibits power-law correlations and the power spectrum ()(∝ 1/  ) exponent  is related [70] to the DFA exponent  through  = 2 − 1. Figure 6 depicts the results of the DFA for the   time series deduced from the aftershock magnitude time series of the four stronger EQs of Table 1.We observe that the resulting DFA exponents  are very close to unity, thus showing that   indeed exhibits 1/ noise.
We now turn to the prediction algorithm based on   discussed in the previous section.This is applied in natural time.Thus, the time increased probability (TIP) [90,91] is turned on after the occurrence of the ( − 1)th aftershock when  −1 is equal to or below some threshold   and lasts until the occurrence of the (next) th aftershock with   ≥  thres .Based on Figure 1 of Shcherbakov et al. [13], we set  thres = 2.0.Under these assumptions, the aftershock magnitude prediction reduces to a binary prediction.We can therefore use the Receiver Operating Characteristics (ROC) graph [92] to depict the prediction quality.This is a plot of the  hit rate (or True Positive rate, TPr) versus the false alarm rate (or False Positive rate, FPr), as a function of the total rate of alarms, which is tuned by the threshold   .The hit rate is the ratio of the cases for which TIP was on and   ≥  target over the total number of cases that   ≥  target .The false alarm rate is the ratio of the cases for which TIP was on and   <  target over the total number of cases for which   <  target .Random predictions generate equal hit and false alarm rate, and hence they lead to the diagonal in ROC plot.Thus, only when the points lie above this diagonal a predictor is useful.Of course, random predictions lead to ROC curves that exhibit fluctuations which depend on the positive  cases (when   ≥  target ) and the negative  cases (when   <  target ) to be predicted.The statistical significance of an ROC curve depends [93] on the area  under the curve in the ROC plane (cf.this area is also usually termed AUC).Mason and Graham [93] have shown that where  follows the Mann-Whitney -statistics [94].Recently, a visualization scheme for the statistical significance of ROC curves has been proposed [95].It is based on ellipses which are the envelopes of the confidence ellipses (cf. a point lies outside a confidence ellipse with probability exp(−/2)) obtained when using a random predictor and vary the prediction threshold.These -ellipses cover the whole ROC plane and using their  we can have a measure [95] of the probability  to obtain by chance (i.e., using a random predictor) an ROC curve passing through each point on the ROC plane.In the present case, the  values that will be reported below correspond to the whole ROC curve and are calculated on the basis of its AUC using Eq. ( 2) by means of the computer program VISROC.f available in Sarlis and Christopoulos [95].This  value gives the probability of obtaining by chance a prediction scheme that leads to an AUC equal to the observed one for given values of  and .It is worthwhile to mention that in the standard method, by using the Gutenberg-Richter law which is a skewed distribution for which ROC graphs are [92] especially useful, the magnitudes are assumed to be independent, thus unpredictable, and this remains so even if we consider the so-called [96] shortterm aftershock incompleteness (STAI) which mainly arises because directly after a large EQ the overlapping arrivals of waves from different events can be masked [97][98][99] thus temporarily increasing the completeness threshold of the EQ catalogue (see also the discussion below).On the other hand, the present method tries to predict large aftershocks from the correlation between magnitudes.Below in this paper, we focus on evaluating the statistical significance of this method.Other methods may be also needed to evaluate its practical utility.

The Case of the Four Stronger EQs
( ≥ 6.7) in Southern California.These are the cases of the 1992 Landers, the 2010 El Mayor-Cucapah, the 1999 Hector Mine, and the 1994 Northridge EQ.As mentioned the analysis is made by considering the areas shown in Table 1 centered at the epicenter of each mainshock and for a period  equal to one year.The aftershocks thus obtained (see Figure 1) have been sorted according to their origin time and the time series   versus the order of occurrence  has been constructed.Figure 5 shows an example of the time series   obtained from the successive extrema of the aftershock magnitude time series   of the 1992 Landers EQ.We then apply the prediction algorithm as described in the previous section for various target magnitudes  target : when  −1 is smaller than or equal to   the TIP turns on: if the magnitude of the (next) aftershock is   ≥  target we have a true positive successful prediction and if   <  target we have a false positive unsuccessful prediction.When TIP is off, that is,  −1 >   , and (%) Figure 7: Receiver Operating Characteristics (red squares) when using   as a predictor for  target =   − 2 for the aftershock time series of Landers, El Mayor-Cucapah, Hector Mine, and Northridge EQs (see Table 1).In each case, the colored contours present the  value to obtain by chance an ROC curve based on -ellipses [95]; the -ellipses with  = 10%, 5%, and 1% are also shown.
<  target we have a true negative successful prediction; if   ≥  target we have false negative unsuccessful prediction, a miss.This way we can obtain ROC curves for various  target .We first focus on a high  target , selected on the first decimal digit, that is determined by the formula  target =   − 2, where   is-as mentioned-the mainshock magnitude.The results for these  target 's for the four strong EQs under discussion are shown in Figure 7.The corresponding  values to obtain such ROC curves by chance is 0.012%, 0.036%, 0.81%, and 0.24% for the Landers, El Mayor-Cucapah, Hector Mine, and Northridge EQ, respectively.Figure 8 depicts the ROC curves obtained when selecting a lower target magnitude  target = 4.0.We observe that in this case also we obtain ROC curves of high statistical significance.Finally, in Figures 9(a) and 9(b), we present the results obtained for AUC and the corresponding  values when considering a wide range of  target 's from 4.0 to 5.5.

3.2.
The Case of the Four Smaller EQs (6.6 ≥  ≥ 6.0).As already mentioned, this is the case of the two Superstition Hills large events [53] in 1987 and the 1992 Joshua Tree and Big Bear events which are related [54,55] to the Landers EQ.The areas used for the selection of the aftershock magnitude time series   are shown in the last column of Table 1, but now the period  for the study of aftershocks may vary from one year since these EQs cannot be considered as clearly independent mainshocks.Specifically, in the case of the first Superstition Hills 6.2 EQ, we considered as  the period until the occurrence of the second Superstition Hills 6.6 EQ, that is, roughly 11 hours and 23 minutes.Moreover, since Joshua Tree EQ is clearly related [54,55] to the Landers EQ, we considered as  the period until the occurrence time of the Landers EQ, that is, roughly 66 days and 7 hours.In the remaining two cases treated here, a period  of one year has been used.The aftershocks considered in each case are shown in Figure 2.
Following the procedure used in the previous subsection, we sorted the aftershocks according to their occurrence time and constructed the   time series, and from the latter we obtained the predictor time series   .This time series has been used for the prediction of the aftershock magnitude time series   and the ROC curves obtained are shown in Figures 10 and 11 for  target =   − 2 and  target = 3.5, respectively, for each EQ.On the latter figure, we clearly observe that the predictions made on the basis of the proposed algorithm are  far beyond chance (cf.for the ROC related to the Joshua Tree EQ in Figure 11 we have  = 0.7%).If we focus our attention on the high  thres 's (see Figure 10), larger  values are obtained which are 5%, 0.1%, 0.9%, and 11% for the second Superstition Hills, Big Bear, first Superstition Hills, and Joshua Tree EQs, respectively.These  values are also larger than those obtained in the previous subsection when studying the stronger EQs in Southern California.In Figures 9(c) and 9(d), we present the results obtained for AUC and the corresponding  values when considering a wide range of  target 's from 3.5 to 5.0. Figure 9(d) shows that in the vast majority (45 out) of the 59 examined cases the obtained  values are below 5% pointing to the statistical significance of the method.
At this point we have to comment on the values of   obtained after the last aftershock considered in the case of the first Superstition Hills and the 1992 Joshua Tree EQs.As mentioned our study above was terminated just before the occurrence of the second Superstition Hills and the Landers EQ, respectively, in each case.Hence, these values coincide with the  −1 values just before the occurrence of the latter two EQs and are 11 and 7, respectively.The operating points on the ROC diagram corresponding to   equal to these values are shown with the red arrows in Figure 10.An inspection of this figure shows that depending of the selection of the operating point of the proposed algorithm an alarm could have been probably on almost 23 hours before the occurrence of the Landers EQ (actually after 12:38:20 UTC on June 27, 1992, when an  = 2.36 EQ occurred at 34.046 N 116.287W, i.e., 22 km away from the Landers EQ epicenter) whereas this is rather improbable for the case of the second Superstition Hills EQ since the corresponding to   = 11 ROC point leads to an extremely high false alarm rate.Hence, the selection Joshua Tree (%) Figure 10: Receiver Operating Characteristics (red squares) when using   as a predictor for  target =   − 2 for the aftershock time series of the second Superstition Hills, the Big Bear aftershock, the first Superstition Hills, and the Joshua Tree EQs (see Table 1).In each case, the colored contours present the  value to obtain by chance an ROC curve based on -ellipses [95]; the -ellipses with  = 10%, 5%, and 1% are also shown.The red arrows indicate the points on the ROC diagram that correspond to   = 11 and   = 7 that coincide with the values of  −1 just before the occurrence of the second Superstition Hills and the Landers EQ obtained from the analysis of the aftershocks of the first Superstition Hills and Joshua Tree EQs, respectively.
of   is also very important.The latter point together with its implications for a practical application of the proposed algorithm will be discussed in the next subsection.

Further Analysis of the Results
. We first focus on the fact that the prediction algorithm presented here is statistically significant.This statistical significance is very high when considering a relatively small  target (e.g., see Figures 8 and   11), which also pertains when focusing on  target =   − 2 (e.g., see Figures 7 and 10).A comparison of Figures 9(b) and 9(d) shows that this statistical significance may decrease when we consider the smallest in magnitude EQ, that is, the Joshua Tree EQ.Even in this case, however, we obtain  values which exclude the possibility of a random result.In order to evaluate the statistical significance of the proposed algorithm out of sample and for EQs in regions Figure 11: Receiver Operating Characteristics (red squares) when using   as a predictor for  target = 3.5 for the aftershock time series of the second Superstition Hills, the Big Bear aftershock, the first Superstition Hills, and the Joshua Tree EQs, respectively.In each case, the colored contours present the  value to obtain by chance an ROC curve based on -ellipses [95]; the -ellipses with  = 10%, 5%, and 1% are also shown.
different from California, we proceeded, by employing the same procedure as that described in Materials and Methods, to the analysis of all  ≥ 8.5 EQs in the Earth since 1976 (see Table 2).For this study the data analyzed come from the Global Centroid Moment Tensor (CMT) Project [100,101] that covers global seismicity since 1 January 1976 and we considered  thres = 5.0, as in Sarlis et al. [51].For EQs that took place before 2011, the 1976 to 2010 CMT catalogue was used, whereas for EQs since 1 January 2011 to 1 July 2015 the monthly CMT catalogues have been employed (all these catalogues are publicly available from http://www.globalcmt.org/CMTfiles.html).For reasons of brevity, we will refer to these data as the GCMT catalogue.The AUC and  values obtained for the wide range of  target from 6.0 to 7.5 are depicted in Figure 12 in a fashion similar to that of Figure 9. From the six cases studied, only the results obtained for the 2005 Sumatra-Nias EQ do not appear statistically significant.This fact has been also verified by producing 10 2 randomly shuffled copies of the original time series and comparing whether the AUCs obtained for a shuffled copy time series were higher than those depicted in the majority of cases presented in Figure 12(a) for each EQ: For the 2005 Sumatra-Nias EQ in 54% of the cases the AUCs for the shuffled copies were higher than those depicted in Figure 12(a), while this percentage falls to 9% and 8% for the Sumatra-Andaman and Chile cases.For the other three cases studied the obtained percentages were smaller than 2%.We note, however, that the 2005 Sumatra-Nias EQ occurred almost three months after the Sumatra-Andaman EQ and is included within the rectangular region of the latter EQ (see Table 2).The corresponding  −1 value (almost three and a half days) before the 2005 Sumatra-Nias EQ was 4, which is the same as that before the Big Bear aftershock; see Figure 5(c).An inspection of Figure 12(b) also shows that the  values are much larger than those of Figure 9.This is related to the fact that  thres used for GCMT is much larger than that used in Southern California where an accurate waveform relocated regional catalogue [8,9] has been employed.Using a higher  thres does not allow   to vary significantly since there are no small EQs to increase   .For example, the maximum value of   is 10, 9, 10, 11, 17, and 6 for the six EQs of GCMT presented in the first column of Table 2, respectively, compared to 22, 22, 22, 24, 15, 19, 14, and 15 for the EQs presented in the first column of Table 1 in Southern California.Hence, the threshold   has to be modified accordingly.Figure 13 depicts two selections of operating points, corresponding to fixed   for each catalogue, for the totality of the 14 EQs studied when focusing on the prediction of aftershocks with  ≥   − 2 and could be seen as a summary evaluation for various regions and two different catalogues (we will return to this point later).
Let us now return to the cases in Southern California and discuss the effect of STAI on the statistical significance of the method.In order to correct from this effect, we followed Helmstetter et al. [102,103], considered for each EQ of Table 1 a time varying  thres given by  thres (  , ) = max [  − 4.5 − 0.75 log 10 () , 2.0] , (3) where  is the elapsed time measured in days from the mainshock, that is, Eq. ( 15) of Helmstetter et al. [103], and repeated the calculations focusing only on the aftershocks during the first 24 hours after each mainshock where STAI is stronger.The results are shown in Figure 14 and show that out of the 115 obtained  values only 6 exceed 7% whereas the vast majority, that is, 107 cases, are below 5%.
As an additional check, we examined whether the validity of the proposed method originates from STAI by the following test: each time the magnitude   of a reported aftershock satisfied   ≥  thres (  , ),   was considered as the new element of a time series labeled  orig  (original).At the same time, a magnitude    [≥  thres (  , )] has been randomly selected from the original distribution of   and considered as the new element of a time series labeled  synth  (synthetic).The two time series  orig  and  synth  have been analyzed by ROC using the prediction algorithm and the corresponding AUCs have been compared 10 4 times.The results showed that the percentage of cases for which the AUC of a synthetic time series was larger in the majority of cases presented in Figure 9 than that of the original one is 4.2%, 8.2%, 0.6%, 9.9%, 59%, 0.005%, 7.2%, and 9.5% for the 8 EQs presented in the first column of Table 1, respectively.These results (apart from those concerning the second Superstition Hills 6.6 EQ that will be discussed below) indicate that the observed predictability cannot be fully accounted for by STAI (cf.if that was the case, only up to two of the above 8 numbers Complexity Figure 13: Operating points on the ROC diagram when using   as a predictor for  target =   − 2 for all mainshocks studied (see Tables 1  and 2) (a) when adopting for Southern California (SC)   SC = 5 and for GCMT   GCMT = 3. Panel (b) corresponds to the selection   SC = 7 and   GCMT = 4.
One might also argue that the present results could be alternatively reproduced by synthetic models generated using the Epidemic Type Aftershock Sequence (ETAS) model [105], coupling the main statistical laws describing the occurrence of earthquakes in time, space in magnitude, that is, the Gutenberg-Richter magnitude distribution and the Omori-Utsu law as well as the increase of the number of aftershocks as about 10   with the mainshock magnitude   (or Bath law).Following Zhuang and Touati [106], we used the code etasim.f of SASeis2006 [104]; see also [105,107,108], for ETAS.By simulating magnitude by Gutenberg-Richter law for b = 1 with cutoff magnitude 2.0, we compared the results thus obtained (blue asterisks in Figure 15) with those of the test data provided in SASeis2006 (http://www.ism.ac.jp/ ∼ogata/Ssg/softwares.html).These test data come [104] from the Japan Meteorogical Agency (JMA) EQ catalogue and concern the aftershock data for almost 18 days of the 26 July 2003 EQ of  JMA = 6.2 at the northern Miyagi-Ken, northern Japan, and the aftershock data for almost 6 days of the 16 August 2005 EQ of  JMA = 7.2 offshore of western Fukuoka-Ken, Kyushu, Japan.Both the real test data and the simulating  = 1 time series have been considered with  thres = 2.0 and Figure 15: Comparison of (a) AUC and (b)  values obtained when using   as a predictor for various  target for the real aftershock magnitude time series of Fukuoka and Miyagi EQs and an ETAS simulation with  = 1 produced by etasim.f of SASeis2006 [104].In all cases,  thres = 2.0 has been used.the cumulative distributions of the observed differences between successive EQs Δ ≡  +1 −  and those obtained Δ * ≡  * +1 − *  from randomly shuffled copies { *  } of the original catalogue.Since such correlations are expected [112] to be stronger for EQs which are closer in time and increases with the magnitude threshold of the catalogue used, we focused our attention on the first 24 hours after the mainshock, applied a time-dependent threshold  thres (  , ), and investigated the presence of such correlations.The cumulative distributions (Δ <  0 ) and (Δ * <  0 ) are depicted by the red (solid) and the green (dashed) curves, respectively, in Figure 16 and they actually differ as it is also shown by the ( 0 ) (red solid) curve in Figure 17.Hence, it is of crucial importance to confirm that the predictability of large aftershocks presented in Figure 14 (for the same time period and the same  thres (  , )) is not caused by the phenomenon observed in Figures 16 and 17.For this reason, we employed the method of surrogates proposed by Schreiber and Schmitz [114] as implemented in the TISEAN package [115] and applied it to the magnitude time series.Such a method creates surrogate data with the same Fourier amplitudes and the same distribution of values.This procedure reproduces the behavior found in Figures 16 and 17 (see the blue dashed curve in Figure 16 and the green, blue, and thick black curves in Figure 17).We generated 10 3 such surrogate magnitude time series, analyzed them by the prediction algorithm, and compared their AUCs with the AUCs of the observed magnitude time series that led to Figure 14.The percentage the AUCs of a surrogate time series was greater than the AUC of the observed magnitude time series in the majority of cases shown in Figure 14 are 9.0%, 0.2%, 0.6%, 43%, 0.1%, 0.1%, 0.7%, and 54% for the cases of the Landers, El Mayor-Cucapah, Hector Mine, Northridge, 2nd Superstition Hills, Big Bear, 1st Superstition Hills, and Joshua Tree EQs, respectively.Such percentages (cf.if we repeat the procedure for the whole study period the percentages for Northridge and Joshua Tree EQs become 14% and 32%, resp.)show that the behavior found in Figures 16 and 17 which is, as mentioned, related to the magnitude correlations found by Davidsen and Green [96] and Lippiello et al. [112] cannot fully account for the observed predictability.Moreover, these results show that the predictability of the aftershock time series is a nonlinear property [116] of the magnitude time series of aftershocks.The latter fact is compatible with the high nonlinearity exhibited by the coherent noise model.
Moreover, as a more direct analysis to avoid the effect of data incompleteness, we excluded the first one day period after the mainshock from the target period of the forecasting, employed time varying threshold  thres (  , ), and examined the predictability in the [1 day, 1 month] time-interval.The corresponding AUCs and  values are shown in Figure 18.These values (apart from those concerning the second Superstition Hills 6.6 EQ) certainly point to the statistical significance of the present method (86 cases out of the 129 lead to  values smaller than 10% with 79 of them below 5%), although they are higher than the  values corresponding to Figure 9.This effect is not incompatible with the behavior exhibited by the coherent noise model.It has been recently shown [21,22] that the expected avalanche size of the th avalanche of the coherent noise model in a statistical ensemble of initially identical systems relaxes to the steady state value by following a power-law decay, described by the Tsallis -exponential [117], almost inversely proportional to .When we started our study in the second day after the mainshock, we excluded almost 10% of the first month aftershocks.These events, however, express the most important part of the offequilibrium behavior which can be best described by the present algorithm as shown in Figure 14.As the order  of the Figure 16: The cumulative distributions of the event increments Δ ≡  +1 −   for the real aftershocks magnitude time series (solid red) together with those corresponding to randomly shuffled magnitude time series (dashed green) and for surrogate magnitude time series (dashed blue).
avalanche increases, the system deviates less from the steady state and hence it is plausible to assume that its predictability decreases.The fact that the predictability decreases when we exclude the first 10% of the avalanches produced and focus on the latter 90% has been also verified by simulations of the coherent noise model.Thus, in summary, we have shown that the observed predictability of the present algorithm cannot be attributed solely to STAI nor is a result included in up to date ETAS modelling nor is a direct consequence of the magnitude correlations studied by Davidsen and Green [96] and Lippiello et al. [112] in Southern California.The present method tries to predict large aftershocks from the correlation between magnitudes which appear to mimic the behavior of the coherent noise model.According to our view, this is the reason behind the predictability observed in Figure 9 that corresponds to a large time period of study as well as in Figure 14 that corresponds to the first day.
If we focus on the strongest aftershock (which may be the most destructive one), the value of the predictor  −1 before its occurrence for all the 14 EQs studied is 4, 3, 4, 0, 2, 6, 1, and 7 for the Landers, El Mayor-Cucapah, Hector Mine, Northridge, second Superstition Hills, Big Bear, first Superstition Hills, and Joshua Tree EQs in California, whereas it is 4, 3, 0, 0, 0, and 0 for the Sumatra-Andaman, Sumatra-Nias, Sumatra Indonesia, Chile, Tohoku Japan, and Indian Ocean EQs in GCMT.Thus, when selecting for Southern California the threshold   SC = 7 and   GCMT = 4 for GCMT, as in Figure 13(b), one could have successfully predicted the strongest aftershock for all the 14 EQs studied with the corresponding operating points above the diagonal in the ROC diagram.Moreover, in 12 out of the 14 cases studied, that is, if we discard the cases of Sumatra-Nias (which, as mentioned, could have been predicted based on the analysis of Sumatra-Andaman) and Indian Ocean EQ, these operating points lie in the left uppermost quadrant of the ROC diagram with an average hit rate 76% and an average false alarm rate 28% (cf. the corresponding average values for all 14 EQs are 77% and 34%, resp.).Similarly, when selecting   SC = 5 and   GCMT = 3 the 12 out the 14 operating points lie in the region FPr < 0.3 and TPr > 0.3; see Figure 13(a).The 14 operating points in this case lead to an average hit rate 61% and an average false alarm rate 16%.
Returning now to the results for Southern California, we have to note that the probability to observe   ≤ 4 in the examined predictor time series is 0.8%, 0.4%, 1.2%, and 2.0% for the four stronger EQs, respectively, thus leading to a number of alarms smaller than 2% of the observed aftershocks (with   ≥ 2.0) in each case (cf. the corresponding total alarm duration in conventional time, which is depicted in Figure 19(a), varies from four hours and forty minutes for the Landers EQ to one hour for the El Mayor-Cucapah EQ).Moreover, when considering the five stronger aftershocks in each case (see Table 3), we observe that by setting   = 8 we can successfully predict (more than) 80% of these strong aftershocks (cf.only the cases labeled L4, E4, H3b, and N4 exhibit  −1 > 8).As an example, we also plot in Figures 19(b) and 19(c) the total alarm duration in conventional time as a function of the time elapsed from the mainshock for the predictor thresholds   = 6 and 8.We observe that, during the first month from the mainshock, which is the most highly prone period to strong aftershocks, this quantity varies from 9-30 hours (  = 6) to 3-4 days (  = 8).

Conclusions
In conclusion, in the present study we presented a simple algorithm for predicting aftershock magnitudes based on the analysis of the coherent noise model and examined its performance in Southern California, where a waveform relocated catalogue [8] from 1981 to June 2011 [9] is available.This study has been extended, for comparison purposes, to the six strongest EQs in the Earth during almost the last forty years according to the GCMT catalogue [100,101].The main conclusions are the following: (1) The predictor time series exhibits the ubiquitous 1/ noise.
(2) The algorithm leads to statistically significant results for 13 out of the 14 EQs studied.Only the case of Sumatra-Nias EQ could be considered by chance, but this EQ might have been anticipated (three and a half days before its occurrence) by the similar analysis for Sumatra-Andaman EQ.
(3) When focusing on predicting aftershocks with  ≥   −2 for all the 14 cases studied, an average behavior with a hit rate 61% and a false alarm rate of 16% can  be achieved since the corresponding percentages may fluctuate for each particular case; see Figure 13(a).
(4) The strongest aftershock of each EQ in Southern California with   ≥ 6.7 could have been successfully predicted with a number of alarms amounting only 2% of the total number of aftershocks (with  ≥ 2.0) and a total alarm duration of less than five hours in each case.
The 3 cases when a stronger EQ occurs within or close to the rectangular area studied during the study period of another EQ have revealed the following.The second Superstition Hills EQ would have been missed, whereas Landers and Sumatra-Nias EQs might have been anticipated.The corresponding selection of operating points on the ROC diagram can be seen in Figure 13(b) and leads to an average behavior with a hit rate 77% and a false alarm rate of 34%.

Figure 1 :
Figure 1: Map of the aftershocks considered (see Table 1) in the case of (a) Landers, (b) El Mayor-Cucapah, (c) Hector Mine, and (d) Northridge EQs whose epicenters are depicted by the largest (black) stars, respectively.
Figure 1: Map of the aftershocks considered (see Table 1) in the case of (a) Landers, (b) El Mayor-Cucapah, (c) Hector Mine, and (d) Northridge EQs whose epicenters are depicted by the largest (black) stars, respectively.

5 4. 5 >Figure 2 :
Figure 2: Map of the aftershocks considered (see Table 1) in the case of (a) 2nd Superstition Hills, (b) Big Bear aftershock, (c) 1st Superstition Hills, and (d) Joshua Tree EQs whose epicenters are depicted by the largest (black) stars, respectively.
Figure 2: Map of the aftershocks considered (see Table 1) in the case of (a) 2nd Superstition Hills, (b) Big Bear aftershock, (c) 1st Superstition Hills, and (d) Joshua Tree EQs whose epicenters are depicted by the largest (black) stars, respectively.

Figure 5 :
Figure 5: (a) The magnitude time series (right scale) of the Landers EQ (open circle) and its aftershocks (red vertical bars) together with the time series of  −1 (blue broken line, left scale).(b) is an excerpt of (a) for the first three and a half hours after the mainshock.(c) The results of (b) depicted in natural time.The green arrows in (b) and (c) indicate the Big Bear aftershock.

98 Figure 6 :
Figure6: DFA results for the   time series after the four stronger EQs (see Table1) in Southern California during the period 1981 to June 2011.For reasons of clarity, log 10 [F DFA ()] have been vertically displaced by one unit starting from the bottom to the top.The least squares straight lines in each case are also shown together with the corresponding DFA exponent .

Figure 8 :
Figure8: Receiver Operating Characteristics (red squares) when using   as a predictor for  target = 4.0 for the aftershock time series of Landers, El Mayor-Cucapah, Hector Mine, and Northridge EQs, respectively.In each case, the colored contours present the  value to obtain by chance an ROC curve based on -ellipses[95]; the -ellipses with  = 10%, 5%, and 1% are also shown.

Figure 9 :
Figure 9: Results obtained from the analysis of the four stronger ((a), (b)) and the four smaller ((c), (d)) EQs in Southern California by ROC diagrams.Panels (a) and (c) depict AUC versus  target while (b) and (d) depict the corresponding  values.

Figure 19 :
Figure 19: Total alarm duration in conventional time as a function of the time elapsed from the occurrence of the mainshock for the Landers (solid red), El Mayor-Cucapah (dashed green), Hector Mine (dashed blue), and Northridge (dotted magenta) EQs when considering as predictor threshold the value   = 4 (a),   = 6 (b), and   = 8 (c).

.77 5.77 5.77 5.77 5.77 5.77 5.77 5.77 5.77 5.77 5.77 5.70 5.70 5.70 5.70 5.70 5.70 5.70 5.70 5.70 5.70 5.00 5.00 5.00 5.00 5.00 4.00 4.00 2.20 4.20 5.49 5.49 5.49 5.41 5.41 4.44 5.00 3.50 3.50 3.00
of the successive extrema of aftershocks are formed in the case of the Landers EQ.In each rectangular box the aftershocks with   ≥  thres are arranged vertically according to their order of occurrence.Each time a new aftershock   (shown in italics) takes place, the members of M  (shown in bold) are obtained by discarding all the previous members, that is, those of M −1 , that are smaller than   .The number of elements of M  is the cardinality |M  | denoted by   .

Table 1
) in Southern California during the period 1981 to June 2011.For reasons of clarity, log 10 [F DFA ()] have been vertically displaced by one unit starting from the bottom to the top.The least squares straight lines in each case are also shown together with the corresponding DFA exponent .

Table 2 :
All EQs with magnitude  ≥ 8.5 during the period 1 January 1976 to 1 July 2015 reported in the GCMT catalogue.The dimensions of the rectangular aftershock area, centered at each epicenter, used here are also shown.

Table 3 :
The 5 stronger aftershocks as determined by truncating magnitude to the first decimal digit for the four stronger EQs ( ≥ 6.7) in Southern California studied (see also Table1) together with the corresponding values of the predictor  −1 .