A stress test to evaluate the usefulness of Akaike information criterion in short-term earthquake prediction

Akaike information criterion (AIC) has been recently adopted to identify possible earthquake precursors in ionospheric total electron content (TEC). According to the authors of this methodology, their technique allows finding abrupt increases (positive breaks) in vertical TEC rate of change 25–80 min before the occurrence of large earthquakes, highlighting a promising implication of AIC method in Mw > 8 earthquakes alert strategies. Due to the relevance of this matter, a lively scientific debate ensued from these results. In this study, we carefully evaluate AIC method potentiality in searching earthquake TEC precursory signatures. We first investigate the dependence of the detected breaks number on the adjustable AIC method parameters. Then, we show that breaks occurrence clusters around specific local times and around moderate and high solar and geomagnetic activity. The outcome of this study is that AIC method is not concretely usable for issuing large earthquakes alerts.

It can be noted that, a few minutes after the earthquakes, sTEC and vTEC time series reported in the two papers mentioned above show evident co-seismic ionospheric disturbances (CIDs) lasting for hours (see Heki 3 Fig. 2; Heki and Enomoto 4 Figs. 2 and 3). CIDs are usually observed in TEC data to appear shortly after very large earthquakes as ionosphere's response to atmospheric pressure waves propagating upward created by the sudden piston-like motion of the ground/ocean surface 6,7 .
According to Masci et al. 2 , the pre-earthquake TEC increase reported in Heki 3 and Heki and Enomoto 4 is identified by using a flawed reference curve as the normal TEC background. To be more precise, Masci et al. 2 put in evidence that: i) In both papers, the normal TEC background is defined by a polynomial fitting of TEC time series excluding 1 hour period of data, from -40 to +20 minutes with respect to the earthquake time (see, e.g., Fig. 8 in Heki and Enomoto 4 ).
ii) CID mainly results as a sudden depletion of the TEC lasting for hours for large earthquakes, to which resonant atmospheric oscillations are superimposed. Kakinami et al. 8 named "ionospheric hole" the sudden long-lasting post-seismic TEC depletion. Astafyeva et al. 9 interpreted the TEC hole as the negative phase on an N-like shape wave, a compression-rarefaction wave generated by the earthquake.
iii) The amplitude and duration of a CID mainly depend on the earthquake magnitude. That is, stronger earthquakes cause larger and longer TEC perturbations 9 . Usually, CIDs generated by large earthquakes are observed lasting from 1 to a few hours before TEC recovers the normal state. iv) Considering the 2011 Tohoku-Oki earthquake case studies, sTEC and vTEC time series reported by Heki 3 and Heki and Enomoto 4 show the onset of a CID effect at about 10 minutes after the shock. This is followed by a significant post-seismic depletion (the hole) that lasts at least 2 hours after the shock, more than the 20 minutes period of data excluded in the Heki's fitting procedure. TEC curves also show resonant atmospheric oscillations superimposed to the depletion. Saito et al. 10 show that the area of the TEC depletion after the Tohoku-Oki earthquake had horizontal size of about 500 km.
Masci et al. 2 concluded that it is evident that such post-seismic long-lasting depletion strongly influences the definition of the reference curve. This because, Heki 3 and Heki and Enomoto 4 used in the fitting procedure TEC data with a deep post-seismic depletion. This generates an apparent pre-seismic TEC increase that results in the identification of a flawed pre-earthquake TEC increase.
We would like to put in evidence, that Heki and Enomoto 1 , while trying to refute the Masci's criticism, at page 7007 state that: • The reference curve method is essentially impractical for short-term earthquake prediction. Even if we could constrain the onset of the anomalies, we need the TEC data after earthquakes to pin down the reference curves (extrapolation from the preseismic part is hardly satisfactory).
• As the response to Criticism #1, we will propose a new approach to identify "breaks" (abrupt increase in rate) in absolute vTEC time series as a substitute for the reference curves (section 2.3).
Still at page 7008: • In both Figures 1 and S2, we could intuitively recognize preseismic vTEC increase and postseismic recovery. There, we drew "reference curves" as we did in Heki and Enomoto (2013). Although it became easier to identify the onset of the anomaly by the sTEC to vTEC conversion, we still need the data after earthquakes to draw such reference curves.
Note that, Heki and Enomoto 1 , in a very unusual way for a rebuttal, agree with Masci et al. 2 that the reference curve method is a flawed method for identifying pre-earthquake increases in TEC time-series because data after the earthquake time are required. Therefore, they explore a new method in order to identify anomalous pre-earthquake TEC changes.
In conclusion, we do not see in Heki and Enomoto 1 any rebuttal of the criticism #1.

Mw dependence and reference curve method
Heki 3 (see his Fig. 4b) shows a dependence between the size of the TEC anomaly and the earthquake magnitude Mw, that is, stronger earthquakes are preceded by larger pre-earthquake TEC increases. Masci et al. 2 pointed out that the dependence of the anomaly size with Mw is seemingly induced by the reference curve method. To be more precise, Astafyeva et al. 9 have shown that the depth of the post-seismic TEC depletion mainly depends on the earthquake magnitude. That is, stronger earthquakes trigger deeper post-seismic TEC depletions. Therefore, by using the Heki's reference curve method, a deeper post-seismic TEC depletion results in the identification of a larger pre-earthquake TEC increase.
Note that at page 7010 Heki and Enomoto 1 state that Heki 3 , while claiming a direct dependence between the earthquake magnitude and precursor size, in its Fig. 4b, compared sizes of the precursors of four earthquakes as the cumulative departure of the data from the reference curves at their occurrence times. This quantity, however, depends on the definition of reference curves.

2/14
Once again, unexpectedly for a rebuttal, Heki and Enomoto 1 agree with Masci's remark maintaining that the reference curve method identifies pre-earthquake increases whose amplitude in induced by analysis procedure and also leads to an apparent direct relationship between Mw and the pre-earthquake anomaly size.
Criticism #2: natural variability of TEC As a further investigation of the TEC anomalous pre-earthquake increase shown in Heki 3 before the 2011 Tohoku-Oki earthquake, Masci et al. 2 performed a superposed epoch analysis of sTEC curves (± 30 days around the earthquake time). They show that the TEC increase reported by Heki 3 was not particularly anomalous, but it may be explained in terms of normal global-scale ionospheric variability.
Concerning to this criticism, Heki and Enomoto 1 maintain that: • page 7007: Rebuttal to Criticism #2 is not straightforward because we agree that the natural variability overwhelms the precursors in terms of amplitudes especially when geomagnetic activity is high.
• page 7018: We did not simply rebut to Criticism #2.  2 strongly stressed that if the 40 minutes exceptional repeatability had seismic origin, it should be related to a geological and tectonic context that must be necessarily common to the areas where the earthquakes occurred. That is, strong earthquakes occurred in different regions of the Earth were able to activate a mechanism that generated pre-earthquake TEC anomalies at the same lead-time of 40 minutes. However, this cannot be considered a realistic scenario because the earthquakes investigated by Heki (2011) and Heki and Enomoto 4 had very inhomogeneous seismological properties, and were generated in completely different tectonic and geodynamic scenarios.
Anyway, while not considering the geological context, by excluding TEC data from -40 to +20 minutes with respect to the earthquake time in defining the reference curve, is evident that Heki 3 and Heki and Enomoto 4 a priori considered TEC data outside this range as the normal TEC background. Consequently, they a-priori make a decision regarding the onset time, the observation period, and the amplitude of the pre-earthquake TEC anomaly they are looking for.
Note also that Heki and Enomoto 1 by using the AIC method for the same cases studies of Heki 3 and Heki and Enomoto 4 , where the curve fit method is used, identified different onset times for pre-earthquake TEC increases. For instance, for the 2004 Sumatra-Andaman earthquake, the onset time changes from -40 to -80 minutes, while, for the 2014 Iquique it goes from -40 to -25 minutes. Only for Tohoku-Oki earthquake, Heki and Enomoto 1 confirmed the break detection 40 minutes prior to the shock.
In summary, in Heki and Enomoto 1 there is no rebuttal to the Masci's criticism #3. Instead, it is evident that the Heki and Enomoto's 1 AIC analysis does not support the 40 minutes repeatability of the onset of pre-earthquake anomalies claimed by Heki 3

and Heki and Enomoto 4 .
Criticism #4: the anomaly in the magnetic declination at the time of the 2011 Tohoku earthquake Heki and Enomoto 4 (Figs. 4 and A3) compared magnetic data of Kakioka and Kanoya station at the time of the Tohoku-Oki earthquake. They show a pre-earthquake anomaly in the magnetic declination D at Kakioka. The D anomaly seems to arise almost simultaneously with the TEC anomaly 40 minutes before the earthquake. No anomalies seem to be presents in components H and Z, and in the total field F. According to them, this suggests that the earthquake-related perturbing magnetic field is dominantly in east-west direction.
At page 7007 Heki and Enomoto 1 state that Masci et al. 2 commented on the geomagnetic field and thought it unlikely that the pre-seismic anomaly ⇠ 40 min before the Tohoku-Oki earthquake is seen only in the declination time series (Criticism #4). They also maintain on that this criticism seem to come from misunderstandings by Masci et al. (2015).
The main remark by Masci et al. 2 concerned the idea that the perturbing magnetic field was dominantly (i.e., almost polarized) in the east-west direction. According to Masci et al. 2 , it is hard to accept that: i) The existence of an earthquake-related effect able to generate polarized magnetic disturbances.
ii) The idea that alleged earthquake-related polarized disturbances may preserve polarization before being observed.

3/14
Finally, Masci et al. 2 (Fig. 5) put in evidence that Heki and Enomoto 4 identified the D anomaly by using a curve fit method similar to that they used for the TEC anomaly. Therefore, it is very likely that also the D anomaly should be an artifact due to the fitting procedure, where, as for TEC, both the onset and the end time of the anomaly are imposed by the authors.
Note that, Heki and Enomoto 1 at page 7018 agree with the Masci's conclusion stating that: It is true that only declination showed the "clear" changes with the reference curve method. Thus, once again, we do not see any rebuttal of the criticism #4.

A further Heki's remark on Masci et al. (2015): sTEC or vTEC?
Heki and Enomoto 1 (page 7006) criticized the choice of Masci et al. 2 which analyze sTEC time series stating that they did not give a reason why they did not use absolute vTEC, but use sTEC.
According to Heki and Enomoto 4 , vTEC time series allow a better identification of pre-earthquake anomalies because are free from TEC U-shaped long-term changes seen in sTEC. We would like to point out that Masci et al. 2 after testing the results by Heki 3 did not consider necessary to replicate their analysis using vTEC data. This because, by using the curve fit method, both Heki 3 and Heki and Enomoto 4 show the same flawed 40 minutes repeatability for the onset time of pre-earthquake TEC increases.
Note that Heki and Enomoto 1 at page 7008 state: In both Figures 1 and S2, we could intuitively recognize preseismic VTEC increase and postseismic recovery. There, we drew "reference curves" as we did in Heki and Enomoto 4 . Although it became easier to identify the onset of the anomaly by the sTEC to vTEC conversion, we still need the data after earthquakes to draw such reference curves.
Concerning this point, we would like to stress that is surprising that in Kelley et al. 11 , a paper co-authored by Heki, the authors: i) Analyzed sTEC and not vTEC data, which according to Heki and Enomoto (2013) are more suitable for identifying pre-earthquake TEC changes.
ii) In addition, concerning Criticism #1, they authors still use the (flawed) reference curve method despite that previously Heki and Enomoto 1 have hypothesized the AIC method is more efficient in the search of pre-earthquake increases in vTEC time series.
Therefore, we can conclude that the conversion from sTEC to vTEC cannot be considered a key point in the study of the occurrence of pre-earthquake increases in TEC time series, but this is just a reorganization of the same data. This because, both using sTEC and vTEC, to define the reference curve in the Heki's papers authors needs of post-seismic TEC data perturbed by a sudden long-lasting post-seismic depletion.

Conclusions
Here we addressed the alleged Heki and Enomoto 1 rebuttal to criticism by Masci et al. 2 . The aim of a rebuttal is to refute what has been reported in a previous published paper. However, we have shown that Heki and Enomoto 1 , while suggesting a new numerical approach in the search of pre-earthquake anomalies in TEC data, cannot be actually considered a direct rebuttal of Masci et al. 2 criticisms on Heki 3 and Heki and Enomoto 1 , neither Masci's criticisms can be considered to come from authors misunderstandings, as claimed by Heki and Enomoto 1 .

Supplementary Discussion S2: Description and processing of GPS-TEC data
Receiver-independent exchange (RINEX) files containing GPS code and carrier phase observables acquired every 30 seconds are used to obtain calibrated TEC values by applying the Ciraolo et al. 12 method, improved several times in the years (details at htt ps : //drive.google.com/ f ile/d/1tDn0kqtJSZllJbrHmZC9T LiJUEb3bo8 /view?usp = sharing). To avoid problems related to the use of leveled observations, this method tries to estimate the phase offsets b arc affecting the observations of the differential phase delay of each arc, i.e. S arc S arc = sTEC arc + b arc (S1) To accomplish this task, sTEC values are mapped as a two-dimensional (2D) surface by means of the classical thin shell method (at an altitude that is typically set to 350 km) where vTEC, which is the equivalent vertical TEC, is a 2D unknown function over the thin shell and c the angle formed by the line of sight receiver-satellite and the perpendicular to the shell at the ionospheric pierce point. It is worth highlighting that

4/14
the fact that vTEC is uniquely defined only for observations from a single station limits the applicability of the method to a multi-day, single-station solution. vTEC is expanded as a polynomial, linear in local time and of the fourth-order in Modip, the modified dip latitude proposed by Rawer et al. 13 . In virtue of the aforementioned considerations, the observations S arc can be expressed as S arc = secS n c n p n (LT, where p n is the term of the polynomial of order n and c n the corresponding coefficient. Equation (S3) is linear in the unknown coefficients c n and phase offsets b arc , and so can be solved via standard or more sophisticated least squares methods. Fig. 3 (a-1) by Heki and Enomoto 1 . We have analyzed the same GPS data from GEONET Satellite 15 and receiver 3009 by using the AIC algorithm we developed and the same AIC parameters, i.e., Dt = ±30 and ±40 min, T h R =75%, T h A =3 TECU/hour.  Table 2 in the main text. Ntot indicates the total number of breaks found by AIC method (Dt = 60 minutes, T h A =3 TECU/hour and T h R =75%) in the specified vTEC time series. Step #2