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

Authors: Stavros-Richard G. Christopoulos; Nicholas V. Sarlis
Journal: Complexity
 (2017)

---

## Abstract

Recently, the study of the coherent noise model has led to a simple (binary) prediction algorithm for the forthcoming earthquake magnitude in aftershock sequences. This algorithm is based on the concept of natural time and exploits the complexity exhibited by the coherent noise model. Here, using the relocated catalogue from Southern California Seismic Network for 1981 to June 2011, we evaluate the application of this algorithm for the aftershocks of strong earthquakes of magnitudeM≥6. The study is also extended by using the Global Centroid Moment Tensor Project catalogue to the case of the six strongest earthquakes in the Earth during the last almost forty years. The predictor time series exhibits the ubiquitous 1/f noise behavior.


---

## Body

## 1. 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–4] and natural time [5–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–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.In a time series comprisingN avalanches (or EQs), natural time χk=k/N serves as an index for the occurrence of the kth event [5–7]. Natural time focuses on the sequential order of the events and the analysis is usually made by using the pair (χk,Qk), where Qk is a quantity proportional to the size (and hence to the energy) of the kth event. Natural time analysis has found useful applications in a variety of fields: Statistical Physics (e.g., [17–23]), Cardiology (e.g., [24–27]), Geophysics (e.g., [6, 28–32]), Atmospheric Sciences (e.g., [33, 34]), Seismology (e.g., [35–40]), Physics of EQs (e.g., [41–45]), EQ prediction (e.g., [5, 46–52]), and so on (for a recent review see [7]). For the case of EQs [36], Qk∝10cMk, where Mk is the magnitude of the kth EQ and c≈1.5. The two quantities k(=Nχk) and Mk form the basis of the aftershock magnitude prediction algorithm, which is discussed in the next section.


## 2. Materials and Methods

### 2.1. 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 M≥6.0 events that occurred within this area are shown in Table 1. As seen in this table, there are 4 stronger EQs with M≥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]).Table 1
All EQs with magnitudeM≥6.0 during the period 1981 to June 2011 within the reporting area for local events of SCSN (see the polygon in Figure 1 of Hauksson et al. [8]). The dimensions of the rectangular aftershock area, centered at each epicenter, used in this study are also shown.
EQ name M Origin time (UTC) Lat (°) Lon (°) Area considered Landers 7.30 1992/06/28 11:57:33 34.202 - 116.435 1.1°  ×  1.1° El Mayor-Cucapah 7.20 2010/04/04 22:40:42 32.264 - 115.295 1.0°  ×  1.0° Hector Mine 7.10 1999/10/16 09:46:43 34.595 - 116.271 0.9°  ×  0.9° Northridge 6.70 1994/01/17 12:30:54 34.206 - 118.549 0.6°  ×  0.6° 2nd Superstition Hills 6.60 1987/11/24 13:15:56 33.017 - 115.843 0.5°  ×  0.5° Big Bear aftershock 6.30 1992/06/28 15:05:30 34.202 - 116.828 0.3°  ×  0.3° 1st Superstition Hills 6.20 1987/11/24 01:54:14 33.082 - 115.779 0.3°  ×  0.3° Joshua Tree 6.10 1992/04/23 04:50:22 33.968 - 116.313 0.3°  ×  0.3°Figure 1
Map of the aftershocks considered (see Table1) 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.
(a) (b) (c) (d)Figure 2
Map of the aftershocks considered (see Table1) 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.
(a) (b) (c) (d)Following Shcherbakov et al. [13], we consider as aftershocks all reported events that occurred during a period T 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 100.5Mm 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 Mm 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 T examined was one year for the four stronger mainshocks while it varies for the four smaller EQs (see Results and Discussion).


### 2.2. 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 Na agents; for example, the points of contact in a subterranean fault, for each agent i we associate a threshold xi, i=1,2,…,Na that represents the amount of stress that an agent withstands before it breaks. Without loss of generality [3], xi may come from the uniform distribution in the interval 0≤x<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 pstress(η) and replace all xi that are smaller than η with new values resulting from the uniform distribution in the interval 0≤x<1. The number of agents whose thresholds are renewed is the size S of the avalanche. During the “aging” step, a fixed fraction f 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≤x<1. If we assume that Na→∞, the thresholds xi are represented by a threshold distribution function pthres(x) which initially (k=0) is considered uniform in the interval 0≤x<1; that is, pthres0x=1 (see Figure 3). The size S1 of the first avalanche (k=1) is just the probability Probx<η1=∫0η1pthres0xdx=η1 which represents the “mass” of the agents that had thresholds smaller than η1. After the subsequent “aging” step the threshold distribution becomes pthres1x. Upon repeating the two steps for the second time using η2 we can obtain the size S2 of the second avalanche and the new threshold distribution pthres2x and so on.Figure 3
Schematic diagram of the threshold distribution functionpthreskx after the kth avalanche in the coherent noise model with Na→∞: for k=0 (solid red), k=1 (dashed green) exhibiting a single step at x=0.2, k=2 (dashed blue) with two steps at x=0.1, 0.2, and k=3 (dotted magenta) with a single step at x=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 f=0.1 has been used.In the case of the coherent noise model forNa→∞, the threshold distribution pthreskx after the kth avalanche is a piecewise constant function having [1] the following general form:(1)pthreskx=∑n=0nkdnkΘx-bnk,where nk is the number ofsteps present in the threshold distribution function after the kth avalanche (e.g., see Figure 3), b0k=0, bnk, and dnk are positive real parameters, and Θ(x) is the Heaviside (unit) step function; that is, Θ(x)=0 if x<0 and Θ(x)=1 if x≥0. The analysis of numerical results has shown (e.g., see Figure 5 of [1]) that the integer number nk may serve as apredictor for the size Sk+1 of the next avalanche. It has been also found [1] that nk as a function of k (and hence versus natural time) exhibits the ubiquitous 1/f noise [59–63] (together with a logarithmic trend, e.g., see Eqs. (31) and (32) of [1]).The number of stepsnk changes when a “stress” level η is applied: (a) eliminating at least two smaller “stress” levels previously applied (cf. the case k=3 in Figure 3) or (b) being smaller than the smaller nonzero bnk (cf. the case k=2 in Figure 3). When all ηk<1, nk coincides with the cardinality |Ek| of the sets Ek of successive extrema formed when studying the time series of the random stresses ηk. The sets Ek ofsuccessive extrema are defined [25] as follows: E0 equals the empty set ∅. Each subsequent Ek is obtained by the procedure described below for k times: select a random number ηk from a given probability density function f(η) and compare it with all the members of Ek-1. In order to construct the set Ek, we discard from the set Ek-1 all its members that are smaller than ηk and furthermore include ηk. Thus, Ek≠∅ forallk>0 and Ek is a finite set of real numbers whose members are always larger than or equal to ηk. Moreover, max⁡[Ek]≥max⁡[Ek-1]. The increase of the cardinality |Ek| of these sets is at the most 1, but its decrease may be as large as |Ek|-1, thus exhibiting a characteristic asymmetry as k increases.The aforementionedpredictability of the coherent noise model on the basis of nk could be also seen as follows: when a relatively large random stress ηk 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 x=0 (cf. the case k=3 for x<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 Mk (Mk≥Mthres, where Mthres is a threshold magnitude) as a function of their order of occurrence k, that is, natural time, by assuming ηk=Mk/10. Setting k=0 at the occurrence time of the mainshock, as shown in the upper left rectangular box of Figure 4, the corresponding sets Mk of the successive extrema (as defined above) in the aftershock magnitude time series are constructed and the cardinality ek≡|Mk| (e.g., see Figure 4) plays a role similar to that of nk in the case of the coherent noise model. A simple computer program that calculates ek is presented in Appendix.Figure 4
How the setsMk of the successive extrema of aftershocks are formed in the case of the Landers EQ. In each rectangular box the aftershocks with Mk≥Mthres are arranged vertically according to their order of occurrence. Each time a new aftershock Mk (shown in italics) takes place, the members of Mk (shown in bold) are obtained by discarding all the previous members, that is, those of Mk-1, that are smaller than Mk. The number of elements of Mk is the cardinality |Mk| denoted by ek.Figure5(a) depicts the aftershock magnitude time series of the Landers EQ (see Table 1) versus conventional time together with the time series of ek. The argument of the latter time series has been intentionally shifted by unity in order to compare the value of ek-1, that is, the number of successive extrema before the occurrence of the kth aftershock, with the magnitude Mk of the kth aftershock. We observe a characteristic pattern which reflects the already mentioned asymmetry when ek is considered as time series with respect to thenatural number k, but no clear correlation between Mk and ek-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 ek-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 ek-1 was equal to 4. Thus, we observe that for Landers EQ the suggestion that ek-1 might be a useful predictor could be considered plausible. The subject of the next section is to examine statistically whether ek-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–66].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 ofek-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.
(a) (b) (c)


## 2.1. 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 M≥6.0 events that occurred within this area are shown in Table 1. As seen in this table, there are 4 stronger EQs with M≥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]).Table 1
All EQs with magnitudeM≥6.0 during the period 1981 to June 2011 within the reporting area for local events of SCSN (see the polygon in Figure 1 of Hauksson et al. [8]). The dimensions of the rectangular aftershock area, centered at each epicenter, used in this study are also shown.
EQ name M Origin time (UTC) Lat (°) Lon (°) Area considered Landers 7.30 1992/06/28 11:57:33 34.202 - 116.435 1.1°  ×  1.1° El Mayor-Cucapah 7.20 2010/04/04 22:40:42 32.264 - 115.295 1.0°  ×  1.0° Hector Mine 7.10 1999/10/16 09:46:43 34.595 - 116.271 0.9°  ×  0.9° Northridge 6.70 1994/01/17 12:30:54 34.206 - 118.549 0.6°  ×  0.6° 2nd Superstition Hills 6.60 1987/11/24 13:15:56 33.017 - 115.843 0.5°  ×  0.5° Big Bear aftershock 6.30 1992/06/28 15:05:30 34.202 - 116.828 0.3°  ×  0.3° 1st Superstition Hills 6.20 1987/11/24 01:54:14 33.082 - 115.779 0.3°  ×  0.3° Joshua Tree 6.10 1992/04/23 04:50:22 33.968 - 116.313 0.3°  ×  0.3°Figure 1
Map of the aftershocks considered (see Table1) 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.
(a) (b) (c) (d)Figure 2
Map of the aftershocks considered (see Table1) 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.
(a) (b) (c) (d)Following Shcherbakov et al. [13], we consider as aftershocks all reported events that occurred during a period T 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 100.5Mm 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 Mm 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 T examined was one year for the four stronger mainshocks while it varies for the four smaller EQs (see Results and Discussion).


## 2.2. 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 Na agents; for example, the points of contact in a subterranean fault, for each agent i we associate a threshold xi, i=1,2,…,Na that represents the amount of stress that an agent withstands before it breaks. Without loss of generality [3], xi may come from the uniform distribution in the interval 0≤x<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 pstress(η) and replace all xi that are smaller than η with new values resulting from the uniform distribution in the interval 0≤x<1. The number of agents whose thresholds are renewed is the size S of the avalanche. During the “aging” step, a fixed fraction f 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≤x<1. If we assume that Na→∞, the thresholds xi are represented by a threshold distribution function pthres(x) which initially (k=0) is considered uniform in the interval 0≤x<1; that is, pthres0x=1 (see Figure 3). The size S1 of the first avalanche (k=1) is just the probability Probx<η1=∫0η1pthres0xdx=η1 which represents the “mass” of the agents that had thresholds smaller than η1. After the subsequent “aging” step the threshold distribution becomes pthres1x. Upon repeating the two steps for the second time using η2 we can obtain the size S2 of the second avalanche and the new threshold distribution pthres2x and so on.Figure 3
Schematic diagram of the threshold distribution functionpthreskx after the kth avalanche in the coherent noise model with Na→∞: for k=0 (solid red), k=1 (dashed green) exhibiting a single step at x=0.2, k=2 (dashed blue) with two steps at x=0.1, 0.2, and k=3 (dotted magenta) with a single step at x=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 f=0.1 has been used.In the case of the coherent noise model forNa→∞, the threshold distribution pthreskx after the kth avalanche is a piecewise constant function having [1] the following general form:(1)pthreskx=∑n=0nkdnkΘx-bnk,where nk is the number ofsteps present in the threshold distribution function after the kth avalanche (e.g., see Figure 3), b0k=0, bnk, and dnk are positive real parameters, and Θ(x) is the Heaviside (unit) step function; that is, Θ(x)=0 if x<0 and Θ(x)=1 if x≥0. The analysis of numerical results has shown (e.g., see Figure 5 of [1]) that the integer number nk may serve as apredictor for the size Sk+1 of the next avalanche. It has been also found [1] that nk as a function of k (and hence versus natural time) exhibits the ubiquitous 1/f noise [59–63] (together with a logarithmic trend, e.g., see Eqs. (31) and (32) of [1]).The number of stepsnk changes when a “stress” level η is applied: (a) eliminating at least two smaller “stress” levels previously applied (cf. the case k=3 in Figure 3) or (b) being smaller than the smaller nonzero bnk (cf. the case k=2 in Figure 3). When all ηk<1, nk coincides with the cardinality |Ek| of the sets Ek of successive extrema formed when studying the time series of the random stresses ηk. The sets Ek ofsuccessive extrema are defined [25] as follows: E0 equals the empty set ∅. Each subsequent Ek is obtained by the procedure described below for k times: select a random number ηk from a given probability density function f(η) and compare it with all the members of Ek-1. In order to construct the set Ek, we discard from the set Ek-1 all its members that are smaller than ηk and furthermore include ηk. Thus, Ek≠∅ forallk>0 and Ek is a finite set of real numbers whose members are always larger than or equal to ηk. Moreover, max⁡[Ek]≥max⁡[Ek-1]. The increase of the cardinality |Ek| of these sets is at the most 1, but its decrease may be as large as |Ek|-1, thus exhibiting a characteristic asymmetry as k increases.The aforementionedpredictability of the coherent noise model on the basis of nk could be also seen as follows: when a relatively large random stress ηk 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 x=0 (cf. the case k=3 for x<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 Mk (Mk≥Mthres, where Mthres is a threshold magnitude) as a function of their order of occurrence k, that is, natural time, by assuming ηk=Mk/10. Setting k=0 at the occurrence time of the mainshock, as shown in the upper left rectangular box of Figure 4, the corresponding sets Mk of the successive extrema (as defined above) in the aftershock magnitude time series are constructed and the cardinality ek≡|Mk| (e.g., see Figure 4) plays a role similar to that of nk in the case of the coherent noise model. A simple computer program that calculates ek is presented in Appendix.Figure 4
How the setsMk of the successive extrema of aftershocks are formed in the case of the Landers EQ. In each rectangular box the aftershocks with Mk≥Mthres are arranged vertically according to their order of occurrence. Each time a new aftershock Mk (shown in italics) takes place, the members of Mk (shown in bold) are obtained by discarding all the previous members, that is, those of Mk-1, that are smaller than Mk. The number of elements of Mk is the cardinality |Mk| denoted by ek.Figure5(a) depicts the aftershock magnitude time series of the Landers EQ (see Table 1) versus conventional time together with the time series of ek. The argument of the latter time series has been intentionally shifted by unity in order to compare the value of ek-1, that is, the number of successive extrema before the occurrence of the kth aftershock, with the magnitude Mk of the kth aftershock. We observe a characteristic pattern which reflects the already mentioned asymmetry when ek is considered as time series with respect to thenatural number k, but no clear correlation between Mk and ek-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 ek-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 ek-1 was equal to 4. Thus, we observe that for Landers EQ the suggestion that ek-1 might be a useful predictor could be considered plausible. The subject of the next section is to examine statistically whether ek-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–66].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 ofek-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.
(a) (b) (c)


## 3. Results and Discussion

We first examine whetherek as it results from the aftershock time series exhibits 1/f noise. To this end, we use the Detrended Fluctuation Analysis (DFA) [67–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–72] and heart dynamics [68, 73] to meteorology [74, 75], atmospheric physics [76, 77], geophysics [28, 40, 78, 79], and economics [80–85]. The major advantage of DFA is the systematic elimination of polynomial trends of different order [86–89] (and hence of the much slower aforementioned logarithmic trend in nk). In DFA, a fluctuation measure FDFA(l) is constructed, based on the residuals of a piecewise polynomial fitting to the profile time series, and studied versus the scale l used for the polynomial detrending (e.g., for more details see Section 1.4.3 of [7]). If FDFA(l)∝lα the time series exhibits power-law correlations and the power spectrum S(f)(∝1/fβ) exponent β is related [70] to the DFA exponent α through β=2α-1. Figure 6 depicts the results of the DFA for the ek 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 ek indeed exhibits 1/f noise.Figure 6
DFA results for theek time series after the four stronger EQs (see Table 1) in Southern California during the period 1981 to June 2011. For reasons of clarity, log10⁡[FDFA(l)] 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 α.We now turn to the prediction algorithm based onek 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 (k-1)th aftershock when ek-1 is equal to or below some threshold et and lasts until the occurrence of the (next) kth aftershock with Mk≥Mthres. Based on Figure 1 of Shcherbakov et al. [13], we set Mthres=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 et. The hit rate is the ratio of the cases for which TIP was on and Mk≥Mtarget over the total number of cases that Mk≥Mtarget. The false alarm rate is the ratio of the cases for which TIP was on and Mk<Mtarget over the total number of cases for which Mk<Mtarget. 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 positiveP cases (when Mk≥Mtarget) and the negative Q cases (when Mk<Mtarget) to be predicted. The statistical significance of an ROC curve depends [93] on the area A under the curve in the ROC plane (cf. this area is also usually termed AUC). Mason and Graham [93] have shown that (2)A=1-UPQ,where U follows the Mann–Whitney U-statistics [94]. Recently, a visualization scheme for the statistical significance of ROC curves has been proposed [95]. It is based on k-ellipses which are the envelopes of the confidence ellipses (cf. a point lies outside a confidence ellipse with probability exp⁡(-k/2)) obtained when using a random predictor and vary the prediction threshold. These k-ellipses cover the whole ROC plane and using their A we can have a measure [95] of the probability p 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 p 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 p value gives the probability of obtaining by chance aprediction scheme that leads to an AUC equal to the observed one for given values of P and Q. 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] short-term aftershock incompleteness (STAI) which mainly arises because directly after a large EQ the overlapping arrivals of waves from different events can be masked [97–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.


### 3.1. The Case of the Four Stronger EQs(M≥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 Table1 centered at the epicenter of each mainshock and for a period T equal to one year. The aftershocks thus obtained (see Figure 1) have been sorted according to their origin time and the time series Mk versus the order of occurrence k has been constructed. Figure 5 shows an example of the time series ek obtained from the successive extrema of the aftershock magnitude time series Mk of the 1992 Landers EQ. We then apply the prediction algorithm as described in the previous section for various target magnitudes Mtarget: when ek-1 is smaller than or equal to et the TIP turns on: if the magnitude of the (next) aftershock is Mk≥Mtarget we have a true positive successful prediction and if Mk<Mtarget we have a false positive unsuccessful prediction. When TIP is off, that is, ek-1>et, and Mk<Mtarget we have a true negative successful prediction; if Mk≥Mtarget we have false negative unsuccessful prediction, a miss. This way we can obtain ROC curves for various Mtarget.We first focus on a highMtarget, selected on the first decimal digit, that is determined by the formula Mtarget=Mm-2, where Mm is—as mentioned—the mainshock magnitude. The results for these Mtarget’s for the four strong EQs under discussion are shown in Figure 7. The corresponding p 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 Mtarget=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 p values when considering a wide range of Mtarget’s from 4.0 to 5.5.Figure 7
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=Mm-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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown.
(a) (b) (c) (d)Figure 8
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown.
(a) (b) (c) (d)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 versusMtarget while (b) and (d) depict the corresponding p values.
(a) (b) (c) (d)


### 3.2. The Case of the Four Smaller EQs(6.6≥M≥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 Mk are shown in the last column of Table 1, but now the period T 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 M6.2 EQ, we considered as T the period until the occurrence of the second Superstition Hills M6.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 T 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 T 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 theMk time series, and from the latter we obtained the predictor time series ek. This time series has been used for the prediction of the aftershock magnitude time series Mk and the ROC curves obtained are shown in Figures 10 and 11 for Mtarget=Mm-2 and Mtarget=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 p=0.7%). If we focus our attention on the high Mthres’s (see Figure 10), larger p 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 p 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 p values when considering a wide range of Mtarget’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 p values are below 5% pointing to the statistical significance of the method.Figure 10
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=Mm-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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown. The red arrows indicate the points on the ROC diagram that correspond to et=11 and et=7 that coincide with the values of ek-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.
(a) (b) (c) (d)Figure 11
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown.
(a) (b) (c) (d)At this point we have to comment on the values ofek 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 ek-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 et 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 M=2.36 EQ occurred at 34.046 N 116.287 W, 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 et=11 ROC point leads to an extremely high false alarm rate. Hence, the selection of et 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.


### 3.3. 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 smallMtarget (e.g., see Figures 8 and 11), which also pertains when focusing on Mtarget=Mm-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 p 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 different from California, we proceeded, by employing the same procedure as that described in Materials and Methods, to the analysis of allM≥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 Mthres=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 p values obtained for the wide range of Mtarget 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 102 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 ek-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).Table 2
All EQs with magnitudeM≥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.
EQ name M Origin time (UTC) Lat (°) Lon (°) Area considered Sumatra-Andaman 9.0 2004/12/26 00:58:50 3.30 95.78 7.8° × 7.8° Sumatra-Nias 8.6 2005/03/28 16:09:36 2.09 97.11 4.9° × 4.9° Sumatra, Indonesia 8.5 2007/09/12 11:10:26 - 4.44 101.37 4.4° × 4.4° Chile 8.8 2010/02/27 06:34:15 - 35.85 - 72.71 6.2° × 6.2° Tohoku, Japan 9.1 2011/03/11 05:46:23 38.32 142.37 8.7° × 8.7° Indian Ocean 8.6 2012/04/11 08:38:36 2.33 93.06 4.9° × 4.9°Figure 12
Results obtained for allM≥8.5 EQs according to the GCMT catalogue (see Table 2) by ROC diagrams. Panels (a) and (b) depict the AUC and the corresponding p values, respectively.
(a) (b)An inspection of Figure12(b) also shows that the p values are much larger than those of Figure 9. This is related to the fact that Mthres 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 Mthres does not allow ek to vary significantly since there are no small EQs to increase ek. For example, the maximum value of ek 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 et has to be modified accordingly. Figure 13 depicts two selections of operating points, corresponding to fixed et for each catalogue, for the totality of the 14 EQs studied when focusing on the prediction of aftershocks with M≥Mm-2 and could be seen as a summary evaluation for various regions and two different catalogues (we will return to this point later).Figure 13
Operating points on the ROC diagram when usingek as a predictor for Mtarget=Mm-2 for all mainshocks studied (see Tables 1 and 2) (a) when adopting for Southern California (SC) etSC=5 and for GCMT etGCMT=3. Panel (b) corresponds to the selection etSC=7 and etGCMT=4.
(a) (b)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 Mthres given by(3)MthresMm,t=max⁡Mm-4.5-0.75log10⁡t,2.0,where t 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 p values only 6 exceed 7% whereas the vast majority, that is, 107 cases, are below 5%.Figure 14
Results obtained when focusing on the first 24 hours after the mainshock with time variable thresholdMthres(Mm,t) from the analysis of the four stronger ((a), (b)) and the four smaller ((c), (d)) EQs by ROC diagrams. Panels (a) and (c) depict AUC versus Mtarget while (b) and (d) depict the corresponding p values.
(a) (b) (c) (d)As an additional check, we examined whether the validity of the proposed method originates from STAI by the following test: each time the magnitudeMk of a reported aftershock satisfied Mk≥Mthres(Mm,t), Mk was considered as the new element of a time series labeled Mlorig (original). At the same time, a magnitude Mk′[≥Mthres(Mm,t)] has been randomly selected from the original distribution of Mk and considered as the new element of a time series labeled Mlsynth (synthetic). The two time series Mlorig and Mlsynth have been analyzed by ROC using the prediction algorithm and the corresponding AUCs have been compared 104 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 M6.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 should have been below 10% with probability 96.2% (≈0.98+8·0.1·0.97+4·7·0.12·0.96) due to the binomial distribution).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 10Mm with the mainshock magnitude Mm (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 forb = 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 MJMA=6.2 at the northern Miyagi-Ken, northern Japan, and the aftershock data for almost 6 days of the 16 August 2005 EQ of MJMA=7.2 offshore of western Fukuoka-Ken, Kyushu, Japan. Both the real test data and the simulating b=1 time series have been considered with Mthres=2.0 and the algorithm has been applied by using the code presented in Appendix. Figure 15 depicts the AUCs and the corresponding p values. An inspection of this figure shows that the results obtained for the real data clearly outperform the simulated ones since in only 3 cases (out of 48) the AUC of the ETAS simulated data exceeded that of the real data. Thus, the experimentally observed predictability cannot be reproduced. These results also show that the present algorithm conveys additional information and for this reason it may be used together with standard and well-accepted methods [109, 110] in which the forecasting of aftershocks in time and magnitude domain is carried out by employing the modified Omori law and the Gutenberg-Richter law.Figure 15
Comparison of (a) AUC and (b)p values obtained when using ek as a predictor for various Mtarget for the real aftershock magnitude time series of Fukuoka and Miyagi EQs and an ETAS simulation with b=1 produced by etasim.f of SASeis2006 [104]. In all cases, Mthres=2.0 has been used.
(a) (b)The models discussed in the previous paragraphs do not impose any direct correlation between successive magnitudes apart from those inherited from the combination of the Gutenberg-Richter law and Omori law, the fact that the rate of triggered events increases with the magnitude of the triggering event, that is, the ETAS model and STAI. It is a fact, however, that magnitude correlations exist in both regional [17, 36, 40, 111] and global [38, 39] seismic catalogues. Thus, one should also investigate whether such correlations are responsible for the predictability observed in Figure 9. Particularly, for California these correlations have been extensively studied by Davidsen and Green [96] and Lippiello et al. [112]. Davidsen and Green [96] showed that they cannot be fully accounted for by STAI as was also found in Lippiello et al. [112] (e.g., see their Figure 2). In both of these studies, based on the investigation of the magnitude correlations firstly introduced by Lippiello et al. [113], the common method of investigating such correlations was the construction of the difference δP(m0)=P(Δm<m0)-P(Δm∗<m0) between the cumulative distributions of the observed differences between successive EQs Δm≡Mk+1-Mk and those obtained Δm∗≡Mk+1∗-Mk∗ from randomly shuffled copies {Mk∗} 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 Mthres(Mm,t), and investigated the presence of such correlations. The cumulative distributions P(Δm<m0) and P(Δm∗<m0) 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 δP(m0) (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 Mthres(Mm,t)) 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 103 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.Figure 16
The cumulative distributions of the event incrementsΔm≡Mk+1-Mk 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).
(a) (b) (c) (d) (e) (f)  (g) (h)Figure 17
The quantityδP(m0) for the observed aftershocks magnitude time series (solid red) together with the corresponding quantity resulting from two example surrogate time series (dashed green and dashed blue) and the average δP(m0) resulting of the whole ensemble of the surrogate time series (thick black curve).
(a) (b) (c) (d) (e) (f)  (g) (h)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 thresholdMthres(Mm,t), and examined the predictability in the [1 day, 1 month] time-interval. The corresponding AUCs and p values are shown in Figure 18. These values (apart from those concerning the second Superstition Hills M6.6 EQ) certainly point to the statistical significance of the present method (86 cases out of the 129 lead to p values smaller than 10% with 79 of them below 5%), although they are higher than the p 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 kth 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 q-exponential [117], almost inversely proportional to k. 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 off-equilibrium behavior which can be best described by the present algorithm as shown in Figure 14. As the order k of the 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.Figure 18
Results obtained when considering as target period the [1 day, 1 month] time-interval after the mainshock and time varying thresholdMthres(Mm,t) from the analysis by ROC diagrams. Panel (a) depicts AUC versus Mtarget while (b) the corresponding p values.
(a) (b)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 predictorek-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 etSC=7 and etGCMT=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 etSC=5 and etGCMT=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 observeek≤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 Mk≥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 et=8 we can successfully predict (more than) 80% of these strong aftershocks (cf. only the cases labeled L4, E4, H3b, and N4 exhibit ek-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 et=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 (et=6) to 3-4 days (et=8).Table 3
The 5 stronger aftershocks as determined by truncating magnitude to the first decimal digit for the four stronger EQs(M≥6.7) in Southern California studied (see also Table 1) together with the corresponding values of the predictor ek-1.
EQ name M e k - 1 Origin time (UTC) Landers 7.30  1992/06/28 11:57:33 L1 6.30 4 1992/06/28 15:05:30 L2 5.77 0 1992/06/28 12:00:44 L3 5.70 1 1992/06/28 12:01:15 L4 5.69 9 1992/06/29 14:08:37 L5 5.53 8 1992/06/28 14:43:21 El Mayor-Cucapah 7.20  2010/04/04 22:40:42 E1 5.70 3 2010/04/04 22:50:17 E2 5.43 3 2010/04/04 23:15:13 E3 5.38 6 2010/04/04 23:25:06 E4 5.29 11 2010/04/08 16:44:25 E5 5.26 0 2010/04/04 22:43:01 Hector Mine 7.10  1999/10/16 09:46:43 H1 5.77 4 1999/10/16 09:59:34 H2 5.37 4 1999/10/16 12:57:20 H3a 4.99 8 1999/10/21 01:54:06 H3b 4.99 11 1999/10/22 16:08:47 H5a 4.98 4 1999/10/21 01:54:33 H5b 4.98 3 1999/10/21 01:54:28 H5c 4.98 2 1999/10/16 09:52:53 Northridge 6.70  1994/01/17 12:30:54 N1 5.89 0 1994/01/17 12:31:56 N2 5.58 8 1994/01/17 23:33:30 N3 5.24 4 1994/01/18 00:43:08 N4 5.24 15 1994/03/20 21:20:12 N5 5.20 3 1994/01/17 12:40:34Figure 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 valueet=4 (a), et=6 (b), and et=8 (c).
(a) (b) (c)


## 3.1. The Case of the Four Stronger EQs(M≥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 Table1 centered at the epicenter of each mainshock and for a period T equal to one year. The aftershocks thus obtained (see Figure 1) have been sorted according to their origin time and the time series Mk versus the order of occurrence k has been constructed. Figure 5 shows an example of the time series ek obtained from the successive extrema of the aftershock magnitude time series Mk of the 1992 Landers EQ. We then apply the prediction algorithm as described in the previous section for various target magnitudes Mtarget: when ek-1 is smaller than or equal to et the TIP turns on: if the magnitude of the (next) aftershock is Mk≥Mtarget we have a true positive successful prediction and if Mk<Mtarget we have a false positive unsuccessful prediction. When TIP is off, that is, ek-1>et, and Mk<Mtarget we have a true negative successful prediction; if Mk≥Mtarget we have false negative unsuccessful prediction, a miss. This way we can obtain ROC curves for various Mtarget.We first focus on a highMtarget, selected on the first decimal digit, that is determined by the formula Mtarget=Mm-2, where Mm is—as mentioned—the mainshock magnitude. The results for these Mtarget’s for the four strong EQs under discussion are shown in Figure 7. The corresponding p 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 Mtarget=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 p values when considering a wide range of Mtarget’s from 4.0 to 5.5.Figure 7
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=Mm-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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown.
(a) (b) (c) (d)Figure 8
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown.
(a) (b) (c) (d)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 versusMtarget while (b) and (d) depict the corresponding p values.
(a) (b) (c) (d)


## 3.2. The Case of the Four Smaller EQs(6.6≥M≥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 Mk are shown in the last column of Table 1, but now the period T 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 M6.2 EQ, we considered as T the period until the occurrence of the second Superstition Hills M6.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 T 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 T 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 theMk time series, and from the latter we obtained the predictor time series ek. This time series has been used for the prediction of the aftershock magnitude time series Mk and the ROC curves obtained are shown in Figures 10 and 11 for Mtarget=Mm-2 and Mtarget=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 p=0.7%). If we focus our attention on the high Mthres’s (see Figure 10), larger p 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 p 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 p values when considering a wide range of Mtarget’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 p values are below 5% pointing to the statistical significance of the method.Figure 10
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=Mm-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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown. The red arrows indicate the points on the ROC diagram that correspond to et=11 and et=7 that coincide with the values of ek-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.
(a) (b) (c) (d)Figure 11
Receiver Operating Characteristics (red squares) when usingek as a predictor for Mtarget=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 p value to obtain by chance an ROC curve based on k-ellipses [95]; the k-ellipses with p=10%, 5%, and 1% are also shown.
(a) (b) (c) (d)At this point we have to comment on the values ofek 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 ek-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 et 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 M=2.36 EQ occurred at 34.046 N 116.287 W, 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 et=11 ROC point leads to an extremely high false alarm rate. Hence, the selection of et 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.


## 3.3. 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 smallMtarget (e.g., see Figures 8 and 11), which also pertains when focusing on Mtarget=Mm-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 p 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 different from California, we proceeded, by employing the same procedure as that described in Materials and Methods, to the analysis of allM≥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 Mthres=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 p values obtained for the wide range of Mtarget 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 102 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 ek-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).Table 2
All EQs with magnitudeM≥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.
EQ name M Origin time (UTC) Lat (°) Lon (°) Area considered Sumatra-Andaman 9.0 2004/12/26 00:58:50 3.30 95.78 7.8° × 7.8° Sumatra-Nias 8.6 2005/03/28 16:09:36 2.09 97.11 4.9° × 4.9° Sumatra, Indonesia 8.5 2007/09/12 11:10:26 - 4.44 101.37 4.4° × 4.4° Chile 8.8 2010/02/27 06:34:15 - 35.85 - 72.71 6.2° × 6.2° Tohoku, Japan 9.1 2011/03/11 05:46:23 38.32 142.37 8.7° × 8.7° Indian Ocean 8.6 2012/04/11 08:38:36 2.33 93.06 4.9° × 4.9°Figure 12
Results obtained for allM≥8.5 EQs according to the GCMT catalogue (see Table 2) by ROC diagrams. Panels (a) and (b) depict the AUC and the corresponding p values, respectively.
(a) (b)An inspection of Figure12(b) also shows that the p values are much larger than those of Figure 9. This is related to the fact that Mthres 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 Mthres does not allow ek to vary significantly since there are no small EQs to increase ek. For example, the maximum value of ek 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 et has to be modified accordingly. Figure 13 depicts two selections of operating points, corresponding to fixed et for each catalogue, for the totality of the 14 EQs studied when focusing on the prediction of aftershocks with M≥Mm-2 and could be seen as a summary evaluation for various regions and two different catalogues (we will return to this point later).Figure 13
Operating points on the ROC diagram when usingek as a predictor for Mtarget=Mm-2 for all mainshocks studied (see Tables 1 and 2) (a) when adopting for Southern California (SC) etSC=5 and for GCMT etGCMT=3. Panel (b) corresponds to the selection etSC=7 and etGCMT=4.
(a) (b)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 Mthres given by(3)MthresMm,t=max⁡Mm-4.5-0.75log10⁡t,2.0,where t 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 p values only 6 exceed 7% whereas the vast majority, that is, 107 cases, are below 5%.Figure 14
Results obtained when focusing on the first 24 hours after the mainshock with time variable thresholdMthres(Mm,t) from the analysis of the four stronger ((a), (b)) and the four smaller ((c), (d)) EQs by ROC diagrams. Panels (a) and (c) depict AUC versus Mtarget while (b) and (d) depict the corresponding p values.
(a) (b) (c) (d)As an additional check, we examined whether the validity of the proposed method originates from STAI by the following test: each time the magnitudeMk of a reported aftershock satisfied Mk≥Mthres(Mm,t), Mk was considered as the new element of a time series labeled Mlorig (original). At the same time, a magnitude Mk′[≥Mthres(Mm,t)] has been randomly selected from the original distribution of Mk and considered as the new element of a time series labeled Mlsynth (synthetic). The two time series Mlorig and Mlsynth have been analyzed by ROC using the prediction algorithm and the corresponding AUCs have been compared 104 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 M6.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 should have been below 10% with probability 96.2% (≈0.98+8·0.1·0.97+4·7·0.12·0.96) due to the binomial distribution).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 10Mm with the mainshock magnitude Mm (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 forb = 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 MJMA=6.2 at the northern Miyagi-Ken, northern Japan, and the aftershock data for almost 6 days of the 16 August 2005 EQ of MJMA=7.2 offshore of western Fukuoka-Ken, Kyushu, Japan. Both the real test data and the simulating b=1 time series have been considered with Mthres=2.0 and the algorithm has been applied by using the code presented in Appendix. Figure 15 depicts the AUCs and the corresponding p values. An inspection of this figure shows that the results obtained for the real data clearly outperform the simulated ones since in only 3 cases (out of 48) the AUC of the ETAS simulated data exceeded that of the real data. Thus, the experimentally observed predictability cannot be reproduced. These results also show that the present algorithm conveys additional information and for this reason it may be used together with standard and well-accepted methods [109, 110] in which the forecasting of aftershocks in time and magnitude domain is carried out by employing the modified Omori law and the Gutenberg-Richter law.Figure 15
Comparison of (a) AUC and (b)p values obtained when using ek as a predictor for various Mtarget for the real aftershock magnitude time series of Fukuoka and Miyagi EQs and an ETAS simulation with b=1 produced by etasim.f of SASeis2006 [104]. In all cases, Mthres=2.0 has been used.
(a) (b)The models discussed in the previous paragraphs do not impose any direct correlation between successive magnitudes apart from those inherited from the combination of the Gutenberg-Richter law and Omori law, the fact that the rate of triggered events increases with the magnitude of the triggering event, that is, the ETAS model and STAI. It is a fact, however, that magnitude correlations exist in both regional [17, 36, 40, 111] and global [38, 39] seismic catalogues. Thus, one should also investigate whether such correlations are responsible for the predictability observed in Figure 9. Particularly, for California these correlations have been extensively studied by Davidsen and Green [96] and Lippiello et al. [112]. Davidsen and Green [96] showed that they cannot be fully accounted for by STAI as was also found in Lippiello et al. [112] (e.g., see their Figure 2). In both of these studies, based on the investigation of the magnitude correlations firstly introduced by Lippiello et al. [113], the common method of investigating such correlations was the construction of the difference δP(m0)=P(Δm<m0)-P(Δm∗<m0) between the cumulative distributions of the observed differences between successive EQs Δm≡Mk+1-Mk and those obtained Δm∗≡Mk+1∗-Mk∗ from randomly shuffled copies {Mk∗} 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 Mthres(Mm,t), and investigated the presence of such correlations. The cumulative distributions P(Δm<m0) and P(Δm∗<m0) 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 δP(m0) (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 Mthres(Mm,t)) 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 103 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.Figure 16
The cumulative distributions of the event incrementsΔm≡Mk+1-Mk 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).
(a) (b) (c) (d) (e) (f)  (g) (h)Figure 17
The quantityδP(m0) for the observed aftershocks magnitude time series (solid red) together with the corresponding quantity resulting from two example surrogate time series (dashed green and dashed blue) and the average δP(m0) resulting of the whole ensemble of the surrogate time series (thick black curve).
(a) (b) (c) (d) (e) (f)  (g) (h)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 thresholdMthres(Mm,t), and examined the predictability in the [1 day, 1 month] time-interval. The corresponding AUCs and p values are shown in Figure 18. These values (apart from those concerning the second Superstition Hills M6.6 EQ) certainly point to the statistical significance of the present method (86 cases out of the 129 lead to p values smaller than 10% with 79 of them below 5%), although they are higher than the p 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 kth 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 q-exponential [117], almost inversely proportional to k. 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 off-equilibrium behavior which can be best described by the present algorithm as shown in Figure 14. As the order k of the 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.Figure 18
Results obtained when considering as target period the [1 day, 1 month] time-interval after the mainshock and time varying thresholdMthres(Mm,t) from the analysis by ROC diagrams. Panel (a) depicts AUC versus Mtarget while (b) the corresponding p values.
(a) (b)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 predictorek-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 etSC=7 and etGCMT=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 etSC=5 and etGCMT=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 observeek≤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 Mk≥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 et=8 we can successfully predict (more than) 80% of these strong aftershocks (cf. only the cases labeled L4, E4, H3b, and N4 exhibit ek-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 et=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 (et=6) to 3-4 days (et=8).Table 3
The 5 stronger aftershocks as determined by truncating magnitude to the first decimal digit for the four stronger EQs(M≥6.7) in Southern California studied (see also Table 1) together with the corresponding values of the predictor ek-1.
EQ name M e k - 1 Origin time (UTC) Landers 7.30  1992/06/28 11:57:33 L1 6.30 4 1992/06/28 15:05:30 L2 5.77 0 1992/06/28 12:00:44 L3 5.70 1 1992/06/28 12:01:15 L4 5.69 9 1992/06/29 14:08:37 L5 5.53 8 1992/06/28 14:43:21 El Mayor-Cucapah 7.20  2010/04/04 22:40:42 E1 5.70 3 2010/04/04 22:50:17 E2 5.43 3 2010/04/04 23:15:13 E3 5.38 6 2010/04/04 23:25:06 E4 5.29 11 2010/04/08 16:44:25 E5 5.26 0 2010/04/04 22:43:01 Hector Mine 7.10  1999/10/16 09:46:43 H1 5.77 4 1999/10/16 09:59:34 H2 5.37 4 1999/10/16 12:57:20 H3a 4.99 8 1999/10/21 01:54:06 H3b 4.99 11 1999/10/22 16:08:47 H5a 4.98 4 1999/10/21 01:54:33 H5b 4.98 3 1999/10/21 01:54:28 H5c 4.98 2 1999/10/16 09:52:53 Northridge 6.70  1994/01/17 12:30:54 N1 5.89 0 1994/01/17 12:31:56 N2 5.58 8 1994/01/17 23:33:30 N3 5.24 4 1994/01/18 00:43:08 N4 5.24 15 1994/03/20 21:20:12 N5 5.20 3 1994/01/17 12:40:34Figure 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 valueet=4 (a), et=6 (b), and et=8 (c).
(a) (b) (c)


## 4. 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 ubiquitous1/f 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 withM≥Mm-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 withMm≥6.7 could have been successfully predicted with a number of alarms amounting only 2% of the total number of aftershocks (with M≥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 Figure13(b) and leads to an average behavior with a hit rate 77% and a false alarm rate of 34%.



---

*Source: 6853892-2017-02-20.xml*