NOAA polar satellites use particle detectors which monitor fluxes of protons and electrons in polar orbits at altitudes between 807 and 854 km (Davis, 2007). The particle detectors (Space Environment Monitor SEM-2) consist of the total energy detector (TED) and the medium energy proton and electron detector (MEPED). The MEPED is composed of eight solid-state detectors measuring proton and electron fluxes from 30 keV to 200 MeV (Evans and Greer, 2004) which include the radiation belt populations, energetic solar particle events, and the low-energy portion of the galactic cosmic ray population. Data can be downloaded at the link http://www.ngdc.noaa.gov/stp/satellite/poes/dataaccess.html. As all of the sets of orbital parameters are provided every 8 s, this value was chosen as the basic time step for our study (Fidani and Battiston, 2008). Consequently, all other variables were defined with respect to the 8-s step. Thus, 8-s averages of the counting rates (CRs), latitude, longitude, MEPED, and omnidirectional data were calculated. Unreliable CRs with negative values were labeled and excluded from the analysis.

To systematically test the methodology proposed for NOAA data, a quantitative and rigorous definition of the concerned precursor (Console, 2001) was established. The daily averages of particle CRs exiting the entrapment in the geomagnetic field (precipitations) were calculated, and then the condition for which a CR fluctuation was not likely, due to possible statistical fluctuations, was set. This calculation was formulated with a probability larger than 99% (Fidani et al., 2010). The sudden increase in particle flux that satisfies this condition was named particle burst (Sgrigna et al., 2005), and for electrons, the same condition was named electron burst (EB). According to a previous work (Aleksandrin et al., 2003), the daily averages of CRs were calculated in the invariant coordinate space. Together with the L-shell and the pitch angle, it was necessary to take into account the CR amplitudes and their variations versus geomagnetic coordinates, since the spatial gradient of particle fluxes near the South Atlantic Anomaly (SAA) is too large (Fidani and Battiston, 2008). CR distributions inside invariant areas are compatible with a Poisson distribution. Being so, an amplitude threshold was introduced for the CRs to define the conditions for which a CR is a non-Poissonian fluctuation with 99% probability. Furthermore, NOAA satellites measure variations of ionospheric parameters not only due to EQs; indeed, they are principally due to solar activity (Sgrigna et al., 2005; Parrot, 2011). To reduce the effects of solar activity, both low values in Dst variations (http://wdc.kugi.kyoto-u.ac.jp/dst_final/index.html) and geomagnetic Ap indexes (https://www.ngdc.noaa.gov/geomag/data.shtml) were chosen to exclude CR data corresponding to the Sun’s influence (Fidani, 2015).

To conclude the quantification of the precursor (Console, 2001), the L-shell invariant parameter was considered to define the magnetic line where a physical interaction, whatever it is, can connect the seismic and ionospheric activities. Following the works by Aleksandrin et al. (2003) and Sgrigna et al. (2005), particle bursts were considered only when their L-shell values referred to M   ≥   6 EQ epicenter projections that occurred on the magnetic lines. This is equivalent to imagining that the physical interaction can occur in the same region near the vertical. The correlation was calculated by filling a histogram with the time differences Δ t between the EQs and EBs, which was indicated by ∑ { E Q ; E B } ( E Q   ×   E B )   following the work of Fidani (2015). This approach was performed by considering only EB on magnetic lines identified as projecting EQ coordinates at different altitudes with respect to EQ epicenters, from −600 km up to 3,200 km in increments of 100 km. A correlation peak at Δ t   =   2 - 3   h started to be significant only for 30–100-keV EBs when considering magnetic line altitudes above 1,400 km and was maximized for 2,200 km (Fidani, 2015) (see Figure 1). Correlations were maximized by using EQs with magnitudes M   ≥   6 downloaded at the link https://earthquake.usgs.gov/earthquakes/search/, located in both the Indonesian and the Philippine regions, having 90°–150° longitudes, with few events in South America. The EBs were detected high off the shore of the United States and the west coasts of South America, at longitudes between 200° and 280°. These different electron positions were associated in a causal way, due to the fact that electrons drift eastward and the EQ positions were located west of the EB detection positions (Fidani, 2018; Fidani, 2020). Being so, if the disturbances which caused electron precipitations from inner radiation belts occurred above the EQ epicenters in the ionosphere, they most likely anticipated the EQ times by 4–6.5 h.

However, to consider EB L-shells around the L-shells corresponding to the EQ epicenter projected at several altitudes constitutes an ambiguity in defining the phenomena that preceded EQs. In fact, the L-shell parameter to choose EB involves both EB and EQ events and each EQ event with a multiplicity of altitude projections. So, the request for the L-shell parameter, which is only related to EBs, is lost and EB cannot be chosen unambiguously. Furthermore, if the condition of L-shell similarity between considered EB detection and EQ projections is not satisfied, the correlation cannot be found. In conclusion, the steps used to calculate any statistical correlation between EBs and EQs are ambiguous, so they cannot be used to define a phenomenon that anticipates strong EQs. In particular, the step involving the L-shell of EBs needs to be modified for the EQ forecasting.

The NOAA-15 MEPED telescope used to monitor the electron flux coming from the zenith in three energy bands in the range of 30 keV–2.5 MeV (Evans and Greer, 2004) has been used. The energy detected for the electrons is a cumulative sum over three thresholds: E 1   =   30   keV , E 2   =   100   keV , and E 3   =   300   keV . Since different energies determine different behaviors in particle dynamics, new energy channels were derived from the difference of the energy thresholds to obtain electrons detected in the intervals 30–100 keV, 100–300 keV, and 300–2.5 MeV. CRs were then corrected for proton contamination (Rodger et al., 2010) from the lower energy range, based on both observations (Asikainen and Mursula, 2008) and simulations (Yando et al., 2011), and using software downloaded from the Virtual Radiation Belt Observatory (http://virbo.org/POES#Processing). Furthermore, the escaping conditions from trapped electrons were determined, thus selecting particles perturbed from the inner Van Allen Belts. These precipitating electrons were identified by calculating their minimum mirror point altitudes, h min , through the UNILIB libraries (Krunglanski, 2003). In fact, if h min   <   100   km along the drift period, the electrons having energies between 30 keV and 3 MeV are ensured to be absorbed in the residual atmosphere. This occurs at the SAA longitudes due to the geomagnetic field asymmetry. Then, electrons drifting eastward and escaping the trapped conditions can be found by enforcing the condition h min   <   100   km on detected CRs. Such electrons cross the NOAA altitudes and are thus able to be detected, up to the 80° longitude from the westward edge of the SAA.

The dynamics of electrons were described using adiabatic invariants such as the geomagnetic field at mirror points B m   =   B / cos ⁡ 2 ⁡ α , where B is the geomagnetic field, the pitch angle α is the difference between the electron velocity and geomagnetic field directions and the L-shell parameter. CRs of precipitating electrons were thus represented in a 4-dimensional matrix (t; L; α; B) including time. The introduction of B was useful for describing the strong spatial variability of the CRs when the satellite entered the SAA (Asikainen and Mursula, 2008). B covered the range of 16–47 µT which was divided into nine nonidentical intervals, shorter where the CR was higher, to better describe the CR spatial variations, and larger where the CR was less frequent, to have a greater number of samples (Fidani and Battiston, 2008). The considered intervals in B were as follows: 16.0–17.5, 17.5–19.0, 19.0–20.5, 20.5–22.0, 22.0–25.0, 25.0–28.0, 28.0–32.5, 32.5–37.0, and 37.0–47.0 µT. Having the SEM-2 detectors with a finite aperture of 30°, the α interval was chosen to be of 15°, dividing the complete excursion into 12 equal intervals. The B-field and the L-shell were re-evaluated on the NOAA-15 orbit using the International Geomagnetic Reference Field (IGRF-12) model (http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html) (Finlay et al., 2015). To concentrate the analysis in the inner Van Allen Belts, L was restricted at the interval of 1.0–2.2 with L-shell steps of 0.1 defining 12 equal intervals. To realize a measure of precipitating electrons which are disturbed from an action coming from the Earth’s surface, external Van Allen Belts are excluded by limiting L under 2.2, and the SAA is also excluded by a minimum value in B = 20.5 µT.

CRs were summed on 8-s time intervals and associated to each adiabatic interval in B, α, and L. CR histograms were created, and their distributions at all the adiabatic regions resulted as Poissonian. Therefore, to obtain less than 1% probability that a CR fluctuation was of a statistical origin, the condition Poisson (CR) < 0.01 had to be satisfied for the average value corresponding to the same adiabatic coordinates of that CR. Thus, such a CR was considered to be a significant fluctuation with a probability greater than 99%, corresponding to the same adiabatic intervals. The small geometric acceptance ( ∼ 0.1   cm 2 sr ) of NOAA detectors required a long time and large adiabatic intervals to obtain sufficient statistics for daily averages. However, in order to obtain a more accurate reading of the particle dynamics, small cell dimensions of adiabatic invariants should be preferred. Being so, an interpolated average value for each adiabatic interval was used to map L and B continuously from the centers of their intervals. In this way, cell sizes were not reduced, and daily averages were accepted only for those cells having at least 20 satellite passages a day. Furthermore, as CRs strongly increase near the SAA, a cubic nonlinear algorithm was used to better interpolate the averages. Starting from averages and variances, it was possible to verify if the NOAA MEPED detected any significant CR fluctuations along the entire satellite orbit.

Electron loss is primarily induced by solar activity; thus, time intervals, when the solar activity influences electron motion inside the internal Van Allen Belts, are excluded from the analysis. The exclusions are obtained by excluding time intervals when the Ap index overcame a threshold which is variable with seasons and years due to the solar cycle. The threshold was fixed by the year and the day of the year using the relation A p   =   11.1   +   0.8   sin [ 0.37 ( y e a r - 1996 ) ]   + { 2.1 - 0.1   sin [ 0.37 ( y e a r - 1996 ) ] }   cos [ 0.0172 ( d a y - 27 ) ] , where a clear phase shift was defined with the minimum of the Sun’s activity in 1996 and the 27-day modulation due to the Sun’s rotation. Moreover, being that the electron flux was related to substorm activity (Lam et al., 2010), CRs were not considered for the analysis when the Dst index was lower than −27 nT satisfying these conditions. The sudden increase of electron CRs were considered EBs influenced by the Earth’ surface, and more EBs detected along the same semi-orbit were considered as one EB.

The correlation between EBs and EQs was calculated after defining L-shells for an EQ  ( L E Q ) by projecting the EQ coordinates to different altitudes and then requiring the condition | L E Q   -   L E B |   ≤   0.1 . This was discussed in Reviewing NOAA Electrons’ Statistical Correlations as a problem for discriminating EBs unambiguously, and another criterion is needed. Moreover, the L-shell condition is equivalent to setting a position above the future EQ epicenter where the detected EB passed. So, it might be the position where electrons escaped the trapped conditions following some remote interaction with the EQ preparation zone in the ground (Fidani, 2020). Given that each L-shell is associated with a well-defined altitude at each geographical point and a physical link is reasonably able to reach a certain maximal altitude, it is plausible that for each L-shell, there corresponds a more-or-less defined interval of geographical coordinates. A plot of EQ latitudes with respect to the L E B indicates this correspondence in Figure 2 (left). The plot shows a quadratic dependence of L E B on the EQ latitude with a minimum around 10 ° . This depends on the shape of the internal Van Allen Belts above the EQ epicenters (see Figure 2; right) that crosses the altitude around 2,000 km with an increasing L-shell, as the latitude moves away from 10 ° . This asymmetry around the equator is produced by the inclination of the geomagnetic field with respect to the rotational axes. Being so, it is enough to select only EBs with L E B in a well-defined interval to guarantee that they correlate with strong EQs in Indonesia and the Philippines.

To verify the validity of the new EB condition, a correlation was recalculated between EQs and EBs, which now includes only EB parameters, with the L E B being in a restricted interval. After a complete study to maximize the correlation with respect to many EQ and EB parameters, the validity of the new condition was confirmed by choosing EBs with the following: 1.21   ≤   L E B   ≤   1.31 , pitch angles 56 °   ≤   α   ≤   74 ° and 108 °   ≤   α   ≤   126 ° , and positions − 35 °  to  15 ° in latitudes and 230 ° - 280 ° in longitudes. For EQs, the depth must be less than 200 km, the latitude in the − 6 °  to  26 ° interval, and the longitude in the 90 ° - 170 ° interval. The correlation was defined by filling the histogram ∑ { E Q ; E B } ( E Q   ×   E B ) with the differences Δ t   =   T E Q   -   T E B between the EQ time T E Q and the EB time T E B , only for those EBs with L in the interval 1.21–1.31. The correlation was also optimized according to both the time binning and the time shift. The optimization corresponded to the Δ t   =   1.5 - 3.5   h interval. Here, the time difference interval, used as the binning, was suggested as the time necessary for the EB to cover the EQ longitude interval of 80° for a 60-keV electron drifting eastward at L   =   1.26 , which was found to be about 2 h. After it, the number of correlation events is increased to 44, thus improving the correlation significance. The updated correlation is shown in Figure 3 (left), and the geographical distribution of correlated EQs is shown on the right. L E B and the EB geographic position were found to be the critical parameters to reveal EB true alarms, and the Δ t   =   1.5 - 3.5   h interval is used from here on out. It should be noted that the peaks around 48 h and around 0 h are consistent with the results of the study by Anagnostopoulos et al. (2012), even if with a low significance.

Following a work by Fidani (2018), the conditional probability of a strong EQ, given the EB measurement, can be calculated using the relation between the covariance and cross correlation (Billingsley, 1995). Moreover, this discussion is valid for binary events using any physical observation other than the EB, which is correlated with the EQ. If applied to the EQ and EB unitary events, the binary correlation is as follows: c o r r ( E Q , E B )   =   cov ( E Q , E B ) P ( E Q ) [ 1   -   P ( E Q ) ] P ( E B ) [ 1   -   P ( E B ) ] , (1) where the cov ( E Q , E B )   =   [ P ( E Q ∩ E B )   -   P ( E Q ) P ( E B ) ] can be explicated throughout the histogram of the EQ to EB coincidences ∑ { E Q ; E B } ( E Q   ×   E B ) in the following population formula: cov ( E Q , E B ) = N h ∑ { E Q ; E B } ( E Q × E B ) - N E Q N E B N h 2 , (2) where N E Q and N E B are the number of EQs and EBs that participated in the correlation, respectively, while N h is the number of total hours divided by two considered for the correlation. Being so, the probabilities of single events are P ( E Q )   =   N E Q / N h  and  P ( E B )   =   N E B / N h . The binary correlation histogram is then calculated as follows: c o r r ( E Q , E B ) = ∑ { E Q ; E B } ( E Q × E B ) / N h - P ( E Q ) P ( E B ) P ( E Q ) [ 1 - P ( E Q ) ] P ( E B ) [ 1 - P ( E B ) ] . (3) Being the joint probability by definition of the covariance and from Eq. 1, we have the following: P ( E Q ∩ E B ) = P ( E Q ) P ( E B ) + c o r r ( E Q , E B ) P ( E Q ) [ 1   -   P ( E Q ) ] P ( E B ) [ 1   -   P ( E B ) ] . (4) The conditional probability P ( E Q | E B )   =   P ( E Q ∩ E B ) / P ( E B ) can be rewritten as follows: P ( E Q | E B ) = P ( E Q ) + c o r r ( E Q , E B ) P ( E Q ) [ 1   -   P ( E Q ) ] [ 1   -   P ( E B ) ] P ( E B ) = ∑ { E Q ; E B } ( E Q   ×   E B ) N E B . (5) which means that if a correlation exists between EQs and EBs which is max Δ t [ ∑ { E Q ; E B } ( E Q   ×   E B ) ] and the time difference Δ t between EQ and EB is chosen to be that of correlations, the probability of an EQ with M ≥ 6 is increased by a term proportional to the correlation.

An equivalent approach to test the results obtained using NOAA particle data refers to the work by Console (2001). Here, a simple definition of an EQ forecasting hypothesis was suggested with a particular sub-volume of the total time–space volume, called the alarm volume, within which the probability of occurrence of strong EQs is higher than the average. Following the work by Console, the prediction related to the detection of a precursor consists in the occurrence of an EQ event of minimal magnitude in the alarm volume. In this framework, the performance of a specific method is carried out through the statistical parameters that can be evaluated in this example, such as the success rate N S / N A , the false alarm rate ( N A -   N S ) / N A , the alarm rate N S / N E , the failure rate ( N E   -   N S ) / N E , and the probability gain (Console, 2001), as shown below: G = N S N A V A V T N E , (6) where N S is the success number, N A is the alarm number, N E is the EQ number, V A is the alarm volume, and V T is the total volume. It should be noted that this description is completely equivalent to the previous being N S   =   max Δ t [ ∑ { E Q ; E B } ( E Q   ×   E B ) ] , N A   =   N E B ,   N E   =   N E Q ,   V A   =   A · N h , where A is the Indonesian and Philippine areas and V A   =   2 A is for the alarm duration of 2 h. V A based on NOAA particle data is constant for all the alarms. Thus, the success rate is exactly P ( E Q | E B ) , the false alarm rate is 1 −   P ( E Q | E B ) , the alarm rate is max Δ t [ ∑ { E Q ; E B } ( E Q   ×   E B ) ] / N E Q , the failure rate is 1   −   max Δ t [ ∑ { E Q ; E B } ( E Q   ×   E B ) ] / N E Q , and the probability gain is as follows: G = P Δ t ( E Q | E B ) P ( E Q ) = 1 + c o r r Δ t ( E Q , E B ) [ 1 / P ( E Q ) - 1 ] [ 1 / P ( E B ) - 1 ] , (7) where c o r r Δ t ( E Q , E B ) is the particular correlation and P Δ t ( E Q | E B ) is the particular conditional probability, both corresponding to a Δ t of 1.5–3.5 h.

A criterion for considering one model more valid than another can be made through the log-likelihood of observing that particular realization of the EQ process: under the hypothesis defining the probabilities of occurrence in P sub-volumes p _i , and c _i being the digital occurrence of at least one event in the sub-volume, the following is the case (Console, 2001): log ( L ) = ∑ i = 1 P { c i ⁡ log [ p i / ( 1 − p i ) ] − log ( 1 − p i ) } . (8) The geographical regions of Indonesia and the Philippines are considered, where strong ( M   ≥   6 ) EQs occurred over 16.5 years from July 1998, which were correlated with the NOAA EBs. The space–time alarm sub-volumes are in this case disjointed and separated in time only, each completely covering both areas for 2-h time intervals from 1.5 to 3.5 h after the EB observations. If so, the complete volume covers N _h hours, of which N _EQ are those where an EQ occurred with c i   =   1 , N E B are those where an alarm occurred with p i   =   P ( E Q | E B ) , N h − N E B are those where no alarm occurred with p i   =   P ( E Q ) , and ∑ { E Q ; E B } ( E Q   ×   E B ) are those where an EQ followed an alarm and where p i   =   P ( E Q | E B ) . The log-likelihood histogram is as follows: log ( L ) = log [ P ( E Q | E B ) 1 − P ( E Q | E B ) ] ∑ { E Q ; E B } ( E Q × E B ) + [ N E Q − ∑ { E Q ; E B } ( E Q × E B ) ] log [ P ( E Q ) 1 − P ( E Q ) ] +  N E B ⁡ log [ 1 − P ( E Q | E B ) ] + ( N h - N E B ) log [ 1 - P ( E Q ) ] , (9) and it is possible to compare this forecasting hypothesis with a simpler model, called the null hypothesis, that is, the Poisson hypothesis. The success rate of this model is constantly P ( E Q ) and the probability gain is always G   =   1.0 , so that the log-likelihood is calculated as follows: log ( L 0 ) = N E Q ⁡ log [ P ( E Q ) ] + ( N h − N E Q ) log [ 1 − P ( E Q ) ] . (10)

The log-likelihood ratio (Martin, 1971) can now be calculated to test the near-Earth space influence hypothesis against the null hypothesis for Δ t of the maximum correlation, evaluating the following: log ( L Δ t L 0 ) = max Δ t [ ∑ { E Q ; E B } ( E Q × E B ) ] log [ P Δ t ( E Q | E B ) P ( E Q ) ]                                             + { N E B − max Δ t [ ∑ { E Q ; E B } ( E Q × E B ) ] } × log [ 1   −   P Δ t ( E Q | E B ) 1 − P ( E Q ) ] . (11)

Moreover, the statistical link between EQs and EBs was tested for its significance, starting from their correlation distribution (Fidani, 2015). In the work by Fidani et al. (2010), it was reported that the EQ-to-EB correlation histogram (Eq. 3), obtained collecting ∑ { E Q ; E B } ( E Q × E B ) for many time differences, can satisfy a Poissonian process when only main shock EQs and semi-orbit EBs are used. Then, indicating the average correlation histogram with Ave, the standard deviation histogram is σ   =   A v e . Being so, the number of standard deviations N σ relative to Δ t , which is shown in Figure 4 of the work by Fidani (2015) for a Δ t lasting 2–3 h with respect to altitude projections, can be evaluated by calculating A v e Δ t and σ Δ t for the same Δ t as follows: N σ = max Δ t [ ∑ { E Q ; E B } ( E Q   ×   E B ) ] − A v e Δ t σ Δ t . (12)

The significance in terms of N σ is shown in Figure 4 for the entire interval of altitude projections corresponding to Figure 1. The significance α c o r r in terms of N σ can be obtained using tables of Poisson probabilities. The significance of the new correlation was also evaluated using N σ . The maximum obtained N σ   =   5.4 corresponded to a probability < 0 . 1 % , not being a statistical fluctuation. A summary of the correlation calculated using Eq. 3, time interval Δ t and probability gains calculated using Eq. 7, and number of events and N σ calculated using Eq. 12, corresponding to the altitude projections and to the new model based on L E B only, are reported in columns 2, 3, 4, 5, and 6 of Table 1, respectively. The values of significance calculated using tables starting from Eq. 12 are also reported in column 7 of Table 1.

Precursory information can be evaluated using Molchan’s error diagram (Molchan, 2003). The quantities needed to characterize the predictive properties of a strategy in an interval ( 0 , T ) are the relative number of failures to predict for an EQ magnitude greater than M, as follows: ν = 1 N ( T ) ∑ 0 < t < T [ 1 − D e ( t ) ] d N ( t ) , (13) and the relative alert time is as follows: τ   =   1 T ∑ 0 < t < T D e ( t ) t a , (14) where t a is the alert time, N ( t ) is the number of seismic events in the interval ( 0 , t ) , d N ( t )   =   N ( t   +   t a ) −   N ( t ) =   { 0   o r   1 } is the number of events in the alert time interval, and D e is a decision alert which can be { 0   o r   1 } in the interval ( t , t   +   t a ) . In the NOAA electron statistical results, t a   =   2   h , N ( T )   =   N E Q , the relative number of failures is the failure rate in the study of Console (2001), being ν E B   =   1   − max Δ t [ ∑ { E Q ; E B } ( E Q   ×   E B ) ] / N E Q , and the relative alert time is τ E B   =   P ( E B ) . The statistical significance α M o l c h of a given point ( ν E B , τ E B ) on the Molchan’s error diagram can be tested using the random probability of the Poisson model, which is the diagonal of the Molchan’s diagram, and is given by the binomial distribution (Kossobokov, 2006) as follows: B [ N E Q ( 1 − ν E B ) | N E Q , τ E B ] = C [ N E Q ; N E Q ( 1 − ν E B ) ] × τ E B N E Q ( 1 τ E B −   1 ) N E Q ν E B . (15) In light of this, the probability of obtaining N E Q ( 1   −   ν E B ) or more hits by chance, as there have been N E Q observed target EQs, is described by the following: α M o l c h = ∑ n = m N B [ n   | N E Q , τ E B ] , (16) which produces the confidence bounds and where the index m   =   max Δ t [ ∑ { E Q ; E B } ( E Q   ×   E B ) ] . G on the Molchan’s error diagram is the slope of the line connecting ( 0,1 ) to ( τ E B , ν E B ) (Zechar and Jordan, 2008), and it is simply calculated as G   =   ( 1   −   ν E B ) / τ E B (Molchan, 1991), which is identical to Eq. 7.

The scenario representing the model of EQ prediction needs to define volumes where EQs occur, where the target volume V _T is 2-d space + 1-d time–space. In this volume, the points of EQ occurrence can be identified, together with alarm volumes V _A , as success (S) and failure of prediction (F) events that are EQs occurring inside or outside V _A , respectively. In this case, a precursor volume V _P containing the alarm events must be defined, which is generally different from V _T ; V _P is the volume of the area where EB detection using NOAA satellites occurs, multiplied by the time of EB observations. An EB detection in V _P is an alarm event which defines V _A . With regard to the correlation mentioned above, for the Indonesian and Philippine latitudes and longitudes, V _T is obtained by multiplying this area by the time spanned by the EQ observations. In this scenario, the occurrence of an EQ event is considered only with M above a magnitude threshold M 0 , where M 0   =   6 . Unlike the models that consider EQs as precursors themselves, in this model, the EB precursor events are detected at different latitudes and longitudes, with respect to those of EQs. Correlations between EBs and EQs occurred for EB detection in the area to the west of the South Atlantic Anomaly. Thus, V _T concerns the longitude interval of 90°–170° and the latitude interval of −6°–26° multiplied by the time interval of the analysis, whereas V _P concerns the area of 230°–280° in longitude and the area of −35°–15° in latitude multiplied by the time interval of the analysis. V _A is generated by an observation of one EB in V _P . It has the same area as V _T multiplied by a duration T = 2 h, which results in the width of the correlation peak, occurring over the next 1.5–3.5 h, which is found to be the time position of the correlation peak. In this case, the alarm volume V _A is constant for all alarms. A success is added if an EQ with M   ≥   6 occurs in the V _A . A failure is added if an EQ occurs out of the V _A , which means not included in the time intervals, and any alarm EB not followed by EQ is classified as a false alarm. The described model can be represented by the three-dimensional space of geographical coordinates and time reported in Figure 5. Performance of the NOAA EB detection is conditioned by a V _P which is not continuous. In fact, as reported above, only a low number of days are magnetically calm enough to be used for the analysis. Moreover, the NOAA satellites go into the detection area west of the South Atlantic Anomaly intermittently, thus further reducing the total time of observations.

This is the cause of a noncontinuous V _T , where V _A appears to fill the same geographical area as V _T for a time interval of 2 h. A V _A within V _T is generated 1.5 h after an EB is observed in V _P . This is different from the model based exclusively on seismic activity, where the causal link between V _A and the precursor is near the vertical, given the seismic properties to cluster. The causal link between V _A and the EB observation event is represented near the horizontal. This causal link is the eastward electron drift, according to the electron energy, which is represented in Figure 5 by green thick arrows moving according to the longitude. The vertical distance between the starting point of arrows above the EQ epicenters and the base of V _A is the time anticipation of a possible physical interaction relating the EQ preparation volume to the ionosphere. In the analysis of 16.5 years of NOAA data, V _T concerned only the Indonesian and the Philippine areas multiplied by N h   =   6,953   h . It should be noted that this value and the following are different from those reported in past publications (Fidani, 2020), as the past reports were rough estimates. In this volume, the following occurred: a total number N E Q   =   N E   =   600 of EQs with M   ≥   6 , a total number N E B   =   N A   =   1,892 of alarms, and a total number N S   =   44 of success. Being so, the success rate of this model is N S / N A   =   0.023 , the false alarm rate is 1   −   ( N S / N A )   =   0.977 , the alarm rate is N S / N E   =   0.073 , the failure rate is 1   −   ( N S / N E )   =   0.927 , and the G   =   3.0 .

Here, the target volume is subdivided into nonoverlapping sub-volumes with time intervals of a day that fill V _T completely. For each day sub-volume, the probability of occurrence of at least one target event is estimated to be equal to P ( E Q ) with no EB observed and P ( E Q | E B ) with one EB. Analyzing the data, days with more than one burst can be found with a frequency of about 20%. These bursts can be far away in time when the time difference is more than 10,000 s (∼2.8 h) or neighbors when the time difference is 5,000–7,000 s (∼1.4–2.0 h); in the latter case, they belong to successive orbits according to NOAA POES orbit parameters. It should be remembered that all EBs in a semi-orbit were considered as only one EB to be counted for the correlation calculus (Fidani, 2015). When two bursts are far away, the time alarm of the first ends before the beginning of the second detection, so two disjoint V _A with the same increasing of the conditional probabilities occur. When detected bursts belong to successive orbits, the V _A time interval of 2 h is greater than the range of 5,000–7,000 s, and a partial overlapping between two consecutive alarms will also occur. The few cases of overlapping for NOAA alarms are not considered here and will be presented in a future publication.

Bayes’ theorem allows us to compute the probability that a hypothesis is true, provided that one knows the probable truth of all supporting arguments. It reverses the conditional probabilities and defines the probability of the hypothesis given the evidence. It shows that there is a significant probability gain in using precursors for prediction, even if a phenomenological occurrence is not the proof of a precursor. However, starting from the correlation results (Fidani, 2015), Bayes’ theorem, as shown below, P ( E B | E Q ) = P ( E B )   P ( E Q | E B ) P ( E Q ) , (17) can be employed using statistical bases (Eq. 5). Being so, the alarm rate can be rewritten as P ( E B | E Q ) , the failure rate as 1   −   P ( E B | E Q ) , and the probability gain as follows: G E B = P ( E B | E Q ) P ( E B ) = G . (18) Bayes theorem allows us to compute the probability of an EB given the measurement of an EQ. If calculated, it appears surprisingly high, equal to 0.8 for the correlation. However, this result must not be misinterpreted. In fact, the probability gain remains the same, suggesting that the high probability of detecting an EB when an EQ is observed is due to its greater frequency of EB occurrence.

Bayes’ theorem tests hypotheses and can be updated on the basis of new information. If used under the condition of many independent precursors EB, EC, ED …, the EQ conditional probability in a small time interval of a given area after the simultaneous detection of one EB, one EC, one ED ... can be approximated by (Aki, 1981) the following: P ( E Q | E B ∩ E C ∩ E D ∩ ... ) ≃ P ( E Q ) G E B G E C G E D · · · . (19)

It is the product among the unconditional probability and the probability gains of each precursor; those are the ratios between the conditional probability of each precursor and unconditional probability. However, the condition of independent precursors is difficult to prove, and from the studies reported in the Introduction, a set of physical links for only a part of them is suspected (see, for example, the study by Pulinets et al. (2015)). Therefore, it seems that dependent candidate precursors represent the most common occurrence. Therein, the conditional probability of an E Q with a magnitude greater than M 0 is not increased by a product of a further probability gain of another detected precursor, if this has a certain degree of dependence on the first to be detected. Thus, the conditional probability cannot be approximated by Eq. 19; it must be recalculated. Starting with only two dependent precursors E B and E C that generate alarms in the same V T , the conditional probability on E Q , given the observations of both precursors E B and E C , can be expressed using the relations (Eq. 7) and (Eq. 1) as follows: G E B   ∩   E C = 1 + cov ( E Q , E B ∩ E C ) P ( E Q ) P ( E B ∩ E C ) , (20) where the covariance can be explicated throughout the histogram of E Q to E B ∩ E C coincidences ∑ { E Q ; E B ∩ E C } [ E Q   ×   ( E B ∧ E C ) ] , and by considering the total number N E B ∩ E C of correlated precursor events. Finally, P ( E B | E C ) can be calculated through the c o r r ( E B , E C ) of the relation (Eq. 5), as with all the other conditional probabilities. However, a more interesting question might arise upon using two observation networks whose observables are dependent: what is the overall probability gain upon using observations without differentiating them? In this case, the warning corresponds to a detection from the set of dependent observables E B ∪ E C , which means that the probability, as shown below, P ( E Q | E B ∪ E C ) = P ( E Q ∩ ( E B ∪ E C ) ) P ( E B ∪ E C ) , (21) is conditioned by the observation of an EB or an EC or of both an EB and an EC, indifferently. Then, using simple algebra with P ( E Q ∩ ( E B ∪ E C ) )   =   P ( ( E Q ∩ E B ) ) ∪ ( E Q ∩ E C ) ) , we obtain the following: P ( E Q | E B ∪ E C ) = P ( E Q ∩ E B ) + P ( E Q ∩ E C ) − P ( E Q ∩ E B ∩ E C ) P ( E B ) + P ( E C ) − P ( E B ∩ E C ) , (22) where P ( E Q   ∩   E B )   =   P ( E Q | E B ) P ( E B ) , P ( E Q ∩ E C )   =   P ( E Q | E C ) P ( E C ) , P ( E B    ∩    E C )   =   P ( E B | E C ) P ( E C ) , P ( E Q    ∩    E B    ∩    E C ) =   P ( E Q | E B    ∩    E C ) P ( E B    ∩    E C ) . Considering the relation (Eq. 7), the probability gain due to the observation of an event in E B ∪ E C is as follows: G E B   ∪   E C = G E B P ( E B ) + [ G E C − G E B   ∩   E C P ( E B | E C ) ] P ( E C ) P ( E B ) + [ 1 − P ( E B | E C ) ] P ( E C ) , (23) where G E B and G E C are the probability gains of the single precursor E B and the single precursor E C , respectively. G E B ∩ E C is the probability gain (Eq. 20) due to the observation of both E B and E C events that are correlated between them. Alternatively, using the relations (Eq. 7) and (Eq. 1), we obtain the following: G E B ∪ E C = 1 + cov ( E Q , E B ∪ E C ) P ( E Q ) P ( E B ∪ E C ) , (24) where the covariance can be explicated throughout the histogram of E Q to E B ∪ E C coincidences ∑ { E Q ; E B ∪ E C } [ E Q   ×   ( E B ∨ E C ) ] , and by considering the total number N E B ∪ E C   =   N E B   +   N E C   −   N E B ∩ E C .

It needs to be highlighted that a statistical correlation between the two time series does not generally mean that they are physically related (Aldrich, 1995). A causal link has been hypothesized between EB measurements and EQ occurrence (Fidani, 2018; 2020), which supports the validity of the hypothesis. Although the Δ t   =   T E Q   −   T E B time difference of the correlation is in agreement with the physical migration of electrons eastward, this migration has not yet been observed for EBs correlated with EQs. Furthermore, a physical link between the EQ preparation zone and the ionosphere above the future epicenter, separated by about 2,000 km, has not been demonstrated. Finally, the existence of some physical phenomena occurring at the future EQ epicenter, which is enough to influence the ionosphere, remains only a hypothesis until all of these passages have been fully demonstrated. To discover EBs with the correct times at different longitudes, which would satisfy their physical migration in the ionosphere, more satellites are needed. This verification is currently possible for EBs, as different NOAA satellites fly together, even though cases having suitable satellite positions must be found. Till date, this has not been calculated, as no correlations have been found with EBs selected using other satellite databases during the same periods.

Regarding a physical link, able to cover more than 2,000 km between the Earth’s crust and the ionosphere throughout the atmosphere, magnetic pulses have been hypothesized to influence the high-energy charged particles’ motion by pitch angle diffusion (Galper et al., 1995). Even if other processes have been proposed more recently, such as the injection of radioactive substances and charged aerosols into the atmosphere, leading to a change in the vertical electric current in the atmosphere and to a modification of the electrical field in the ionosphere (Sorokin et al., 2001), from a solid-state physics perspective (Freund, 2011), and radioactive ionization to model the lithosphere–atmosphere–ionosphere–magnetosphere coupling by the latent heat flux (Pulinets et al., 2015), by AGW (Yang et al., 2019; Piersanti et al., 2020). Recent measurements of magnetic pulses on the occasions of strong EQs (Bleier et al., 2009; Orsini, 2011; Nenovski, 2015) have suggested an efficient process which may influence the ionosphere (Fidani et al., 2020). In this latter work, the magnetic data analysis at the Chieti Station of the Central Italy Electromagnetic Network, performed using two independent sample systems of the same signal, showed that a large number of pulses were recorded in the ELF band below 10 Hz with amplitudes mostly in the range of 2.5–80 nT. Specifically, the model proposed for analyzing magnetic pulses consisted of diffused underground electrical currents throughout a conductive strip between the Apennines and the Adriatic Sea. The current required to induce detected pulses is greater than i = 40 kA for the strongest pulses. Unipolar magnetic pulses can be, for example, generated deep in the rock column by peroxy defects when rocks are subjected to increasing deviatoric stresses (Freund et al., 2021). The proposed strip of diffuse current constitutes the electromagnetic source of ULF waves, which is able to produce low intensities of magnetic inductions on the Earth’s surface, even if it is measurable both near and far from the EQ epicenter. The strip of diffuse current is able to produce significant fields, also far from the ground by integration. To demonstrate this and calculate magnetic induction in the ionosphere, two simple models can be considered. An infinitely long strip of about 6-km-thick and 150-km-large conductive soil has been considered to calculate the magnetic induction very close to the larger strip surface. In a real case, a finite-length strip should be utilized, which gives a correct result even for a coil magnetometer very close to the strip. Referring to Figure 6, on the left side, where the lines of induction are generated by a not-finite strip and wire lengths are compared, it is possible to see that the B intensities near the strip are lower than those near the wire. This is due to the current density differences flowing near the observational point. Relations (Eq. 2) and (C4) of the study by Fidani et al. (2020) were used for the wire and the strip, respectively. Moving away from the current density, magnetic inductions generated by the wire and the strip become equal, when the distance overcomes the strip width. However, currents which are not of finite length deviate too much from the real case at large distances. So, a finite length l of the wire is considered and the magnetic induction is simply calculated as follows: B = μ o i l 4 π d ( d 2 + l 2 ) , (25) where d is the distance. Figure 6, on the right side, depicts a comparison of B intensities of the two models with the distance. The B intensity results are limited to 0.2 nT at distances of 1,000 km with a total current of 40 kA. This is in agreement with the intensities obtained from the theoretical calculations, which have shown that only a magnetic type source with frequency < 10–20 Hz can be effective for the penetration of fields into the upper ionosphere and magnetosphere (Molchanov et al., 1995). Moreover, magnetic disturbances were observed above moderate EQ epicenters for the frequency band of 0.1–10 Hz (Strakhov et al., 1994). From the upper ionosphere, these waves travel as Alfven waves further along the geomagnetic field line and reach the inner boundary of the inner Van Allen Belt. Alfven waves are thought to resonantly interact with trapped charged particles; this process is most intensive in the equatorial part of the magnetosphere, for L-shell values equal to or less than 2 (Molchanov and Mazhaeva, 1993). In fact, the measured pulse frequency under 10 Hz is around the bouncing resonance of electrons with energies between 30 and 100 keV (Walt, 1994); it is exactly as measured on board NOAA satellites. This means that energy can be efficiently transferred from magnetic pulses producing electron pitch angle disturbances (Galper et al., 1995). Furthermore, the B intensity due to an Earth surface strip of current is a very low value, compared to an example of geomagnetic activity. However, the frequency range of geomagnetic activity is about 0.0001–0.01 Hz (Francia and Villante, 1997), so that the B rate is out of resonance. Finally, currents produced by lightning are generally lower in intensity, around i = 10 kA, with frequency emissions in the upper part of the ELF band and the VLF band still out of bouncing resonance (Inan et al., 2010).

Indeed, magnetic signals have been correlated with strong EQs in different regions of the world such as Japan (Ohta et al., 2013), Kamchactka (Schekotov et al., 2019), and California (Kappler et al., 2019). The Japanese and Kamchactka studies obtained correlations with time differences of 2–5 days before seismic events, and probability gains of about 1.6 (Han et al., 2014) to 2.6 (Hayakawa et al., 2019) were reported. Magnetic pulses, from here on identified with ECs, have been hypothesized to induce EBs, and therefore, they may be considered dependent events for hypothesis. The relation (Eq. 23) can be used to calculate the probability gain of EQ probability due to possible observations of both EBs and ECs. This possibility can be useful as the observation of ECs on the Earth’s surface can occur when a NOAA satellite is not in a suitable position to reveal the possibly induced EBs. Moreover, magnetic detectors cannot be installed at some positions offshore or are not able to detect ECs where there are EQs. In other words, EBs and ECs can compensate for the gaps on each other in a forecasting experiment, where V _T is very large, as for the West Pacific. The magnetic pulse precursor can be represented by the red volume of Figure 5. Being so, the mutual interrelation of the sets of EBs, ECs, and EQs is a completely real case, which can be studied for a total probability gain with the expression (Eq. 23). To show the advantage in using two dependent observables, a forecasting experiment can be imagined where both EB and EC data are available. Here, one can imagine the presence of a magnetometer network distributed on several islands of the West Pacific. Although not still existent, this network is currently achievable. We suppose that it exists and is able to obtain a G E C   =   1.6 for EQs within a certain distance around the stations, with a time advance of ECs of 4–6 h with respect to EQs. P ( E C ) is necessary for Eq. 23 and is supposed to be 0.05, which means N E C   =   324 magnetic pulses, or sets of magnetic pulses, considered as magnetic alarms on the same time interval as for EBs. For what concerns the NOAA observable, a G E B   =   3 was found from the correlation analysis. However, even if ECs from days with the Ap index above the previously defined threshold were excluded as for EBs, the daily NOAA observation time was always half a day (Fidani, 2020) due to the satellite’s orbit crossing the EB detection region. Being so, the probability gain should be evaluated on a double time interval and being able to detect the same number of EBs. The probability gain is thus approximately G E B   =   1.5 . With regard to P ( E B | E C ) , the correlation between EBs and ECs would need to be calculated using observational data. P ( E B | E C ) depends on the histogram maximum max Δ t ∗ [ ∑ { E B ; E C } ( E B   ×   E C ) ] , where Δ t ∗   =   4 ( 6 )   −   Δ t is such that the total time in advance of ECs with respect to EQs is 4–6 h. This could be chosen in the range from 0, no correlation, to 324, complete correlation, the N E B being a much higher number, in order to consider a range of possibilities. Concerning G E B ∩ E C , the correlation between E Q and E B ∩ E C would need to be calculated utilizing observational data for max Δ t [ ∑ { E q ; E B ∩ E C } ( E Q   ×   ( E B ∧ E C ) ) ] . If E B ∩ E C at Δ t are the set of events EB and EC related by chance, ∑ { E B ; E C } ( E B   ×   E C ) should be a low value with G E B ∪ E C   ≃   [ G E B P ( E B )   +   G E C P ( E C ) ] / [ P ( E B )   +   P ( E C ) ] . The dependence between EBs and ECs can be introduced in Eq. 23 by fixing a time shift in ECs of Δ t ∗ so that the number of correlation events is max   [ ∑ { E B ; E C } ( E B   ×   E C ) ]   =   max Δ t ∗ [ ∑ { E B ; E C } ( E B   ×   E C ) ] . It would also be necessary in this case to consider a range of possibilities from 0, no correlation between E Q and E B ∩ E C , to 10 common events. The probability gain improvements that are obtainable thanks to a network of EC measurements added to the NOAA satellite EB detection are reported in the contour plot of Figure 7.

The maximum value of G E B ∪ E C   =   2.05 is obtained for 314 ECs which are correlated with EBs, but none of these EBs are correlated with EQs, shown by point (314, 0) of Figure 7, which means that 10 ECs are correlated with EQs and another part of EBs not correlated with ECs is correlated with EQs; thus, the total number of different observations correlated with EQs is increased by 10. G E B ∪ E C is a minimum when the correlation events between EBs and ECs are max Δ t ∗ [ ∑ { E B ; E C } ( E B   ×   E C ) ]   =   0 , shown by the point (0, 0) of Figure 7. That is, when there are no ECs correlated with EQs nor with EBs, and G E B ∪ E C is slightly less than 1.6 because it is a weighted average between G _EB and G _EC . When the EB-to-EC correlation increases to 10 events, a maximum of 10 correlations with EQs can happen, as shown by point (10, 10) of Figure 7, where G E B ∪ E C is still lower due to redundancy. G E B ∪ E C does not exist to the left of the white dotted line and on the right of the black dotted line. When the EB-to-EC correlation is at the maximum and the correlation of such ECs with EQs reaches 10 events, as shown by point (324, 10) of Figure 7, G E B ∪ E C returns as near 1.5 due to the low number of correlated EQ events.