Article Highlights

  • FORMOSAT-5/AIP and GIM TEC detect pre-earthquake ionospheric anomalies (PEIAs) and ionospheric storms

  • PEIA associated electric fields of the M7.3 earthquake and prompt penetration electric fields of the storm are 0.3 and 1.0–1.2 mV/m eastward, respectively

  • Dynamo electric fields can be estimated by the EIA crest locations, and computed by ion drift velocities probed by satellites

1 Introduction

Liu et al. (2001) first analyzed the total electron content (TEC) derived by ground‐based global positioning system (GPS) receivers in Taiwan to detect the temporal variation and the spatial distribution of anomalies possibly related to the 1999 M7.6 Chi‐Chi earthquake. Following that, many statistical analyses have been conducted to find the characteristics of temporal pre-earthquake ionospheric anomalies (PEIAs) from GPS TEC measured over Taiwan (Chen et al. 2004; Liu et al. 2004a, b, 2006a, b), Indonesia (Liu et al. 2010a), China (Liu et al. 2009, 2018a, b, c; Chen et al. 2015), and Japan (Kon et al. 2011; Liu et al. 2013a, b). These statistical analyses show that characteristics of the polarity (i.e., decrease or increase; negative or positive), appearance local time, duration, lead day, etc., of PEIAs could vary from place to place. Thus, to detect or identify PEIAs at a certain place, statistical analyses of TEC anomalies and earthquakes during a long-term period are required to first find its associated characteristics (Chen et al. 2015; Liu et al. 2018a). When observed anomalies meet the associated characteristics, it can be declared that a temporal PEIA has been detected. However, the observed PEIAs often suffer from global effects, such as solar disturbances, magnetic storms, etc.

To further discriminate the temporal detected PEIAs from global effects, the TEC in the global ionosphere map (GIM) (Mannucci et al. 1993, 1998) derived from the global navigation satellite system (GNSS) is ideally employed. Currently, the GIM has been routinely publishing with a 1-h time resolution and 1- or 2-day delay (cf. Sun et al. 2017). The spatial resolutions of the GIM covering ± 87.5° latitude and ± 180° longitude are 2.5° and 5°, respectively. Hence, each map consists of 5183 (= 71 × 73) grid points (lattices). After Liu et al. (2001), spatial PEIAs of GIM TEC associated with the 16 October 1999 Mw7.1 Hector Mine earthquake (Su et al. 2013), the 2004 M9.3 Sumatra–Andaman earthquake (Liu et al. 2010a), the 2008 M8.0 Wenchuan earthquake (Liu et al. 2009), the 2010 M7.0 Haiti earthquake (Liu et al. 2011a), and the 11 March 2011 M9.0 Tohoku earthquake (Liu et al. 2018b) have been reported. These studies show that anomalies similar to the detected temporal PEIAs frequently and persistently appear specifically around the epicenter of these earthquakes. Therefore, this global spatial search can be used to not only confirm the detected temporal PEIA but also locate possible forthcoming large earthquakes (Liu et al. 2009, 2010a, 2011a, 2018b; Su et al. 2013).

Satellite observations that provide a global and uniform coverage are ideally employed to monitor PEIAs and ionospheric weather as well as to discriminate local (such as earthquakes, etc.) from global (or space weather from above) effects. DEMETER (Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions) might be the first satellite designed specifically to find PEIAs and, therefore, its main scientific goals were to detect and characterize ionosphere electrical and magnetic anomalies in connection with seismic activities (Elie et al. 1999; Berthelier et al. 2006; Cussac et al. 2006; Lebreton et al. 2006). Many research articles have reported PEIAs in the ion/electron density and ion/electron temperature (Ho et al. 2013a, b, 2018; Ryu et al. 2014a, b, 2015a, b; Liu et al. 2015, 2016a; Shen et al. 2015, 2017; Tao et al. 2017; Yan et al. 2017) as well as amplitude changes of electromagnetic emissions over earthquake regions probed by DEMETER (http://smsc.cnes.fr/DEMETER/A_publications.htm). Besides DEMETER, some other satellites have also been used to find PEIAs of the ion/electron density, ion/electron temperature, magnetic field, and/or ion drift velocity, such as the HINOTORI satellite (Oyama et al. 2008), the DE-2 satellite (Oyama et al. 2011), ROCSAT-1 satellite (Lin et al. 2017; Liu and Chao 2017), and Swarm satellites (De Santis et al. 2017; Marchetti and Akhoondzadeh 2018; Akhoondzadeh et al. 2018).

On 12 November 2017, a severe magnitude M7.3 earthquake struck Iran near the border with Iraq (USGS, https://earthquake.usgs.gov/earthquakes/eventpage/us1000bjnz/executive). The seismic event occurred within the Zagros thrust and fold belt, a tectonically active region between the Arabian and Eurasian plates, with an epicenter near the Iranian town of Ezelgeh. Previous studies have shown that this area is characterized by moderate (M5-6) seismicity rupturing within the 8-to-13-km-thick sedimentary cover (Gombert et al. 2019; Nissen et al. 2019). Derived by the U.S. Geological Survey (USGS) from the most preferred data of the M7.3 earthquake, the focal mechanism indicates a fault plane with an oblique-reverse (oblique-thrust) slip on a 19° dipping angle.

The M7.3 earthquake struck on 12 November 2017, while two magnetic storms commenced on 7 and 21 November 2017, providing the science payload of Advanced Ionospheric Probe (AIP) onboard FORMOSAT-5 (F5) the first opportunity in studying PEIAs and ionospheric storm signatures in the ion density, the ion temperature, and the ion velocity. The GIM TEC routinely published by the Center for Orbit Determination in Europe (CODE) (http://aiuws.unibe.ch/ionosphere/) (cf. Schaer 1999) is used to validate the F5/AIP ion density, find the characteristic of temporal PEIAs at the studied area, as well as search the spatial distributions of anomalies similar to the temporal PEIAs and disturbances of ionospheric positive storm signatures. The F5/AIP ion density, ion temperature, and ion velocity are further examined to find the causal mechanisms of ionospheric weather induced by PEIAs of the M7.3 Iran–Iraq border earthquake and the two ionospheric storms.

2 Observations and Data Analyses

F5 satellite was launched at 18:51 UT (universal time) on 24 August 2017 with a Sun synchronous orbit at 720 km altitude and 98.28° inclination, which passes at almost the same local time of about 10:30 LT (local time) and 22:30 LT, and conducts observations with a 2-day revisit cycle (see Fig. 1) (Chang et al. 2017). The satellite carries a primary payload, the optical RSI (Remote Sensing Instrument), and a science payload, AIP (Advanced Ionospheric Probe), which are normally operating at the daytime of 10:30 LT and nighttime of 22:30 LT, respectively. The AIP mission is to continuously monitor PEIAs and observe ionospheric weathers in the nighttime (Liu and Chao 2017; Lin et al. 2017). F5/AIP has been recording ion quantities since the end of October 2017. Figure 1 displays the F5/AIP ion density, ion temperature, and ion downward/eastward velocities with 1-Hz time resolution in the nighttime during 2–3 November 2017. While the GIM TEC detects temporal and spatial PEIAs, F5/AIP observes concurrent and co-located ion density, ion temperature, and ion drift velocity. These allow us three-dimensionally examining plasma structures and dynamics of lithospheric activity/seismicity and ionospheric weather.

Fig. 1
figure 1

Global observations of F5/AIP a ion density (Ni), b ion temperature (Ti), c downward drift (VD), and d eastward ion drift (VE) probed at 22:30 LT on 2–3 November 2017. The red star denotes the epicenter of the 12 November 2017 M7.3 Iran–Iraq Border earthquake. The black rectangle represents the study area around the epicenter

Full size image

The magnitude M7.3 earthquake (34.9°N, 46.0°E) with a depth of 19 km struck Iran near the border with Iraq at 21:18 LT (18:18 UTC, Coordinated Universal Time) on 12 November 2017. Based on Dobrovolsky et al. (1979) and Su and Liu (2014), to study PEIA associated with the M7.3 Iran–Iraq earthquake, we examine GIM TEC and F5/AIP ion data within a rectangular area 20–50°N × 26–66°E centering on the round the epicenter in detail. Figure 2a from top to bottom displays the solar radio flux of F10.7 and the magnetic indices of AE, Kp, and Dst. The Dst index displays two magnetic storms: one moderate storm occurred at 05:00UT (08:00 LT, post-dawn) on 7 November and reached the maximum depression of − 70 nT on the same day (Storm 1), while the other small storm commenced at 18:00UT (21:00 LT, pre-midnight period) on 20 November and yielded the maximum depression of − 43 nT on 21 November (Storm 2). On 3–4 November, F10.7 is slightly greater than those on the two storm days; the maximum AE value exceeding 596 nT is smaller than those, 945 and 691 nT, on the storm days; the Kp index is about 3+; and the Dst index is about − 24 nT. Thus, it is relatively magnetic quiet on 3–4 November. TECs along the epicenter longitude and over the epicenter are isolated from a sequence of the GIM images in November of 2017 to study the pre-earthquake anomalies and ionospheric positive storm signatures. Figure 2b reveals the latitude–time–TEC plot, a time series of the TEC within ± 45°N along the M7.3 epicenter longitude of 46.0°E, in November of 2017. The EIA (equatorial ionization anomaly) crests move poleward and yield greater TEC values on 2–4, 7, and 21 November 2017, which indicates that the daily dynamo electric fields (Kelley 2009) increase on these days, and seismo- and storm-generated electric fields are in eastward (Lin et al. 2005; Liu et al. 2010b, 2013c; Cheng et al. 2022) during the PEIA and storm days. To identify abnormal signals in the TEC observation (O), a quartile-based process is applied. At each time point, we compute the median, \(\widetilde{M}\), based on the TECs in 15 days before the observation day as well as the first (or lower) and third (or upper) quartiles, denoted by LQ and UQ, respectively. We set the lower bound (LB) as,

$${\text{LB}} = \tilde{M} - \, k(\tilde{M} - {\text{LQ}})$$
(1a)

and upper bound (UB) as,

$${\text{UB}} = \tilde{M} + k({\text{UQ}} - \tilde{M}).$$
(1b)
Fig. 2
figure 2

The solar radio flux, magnetic condition, and ionospheric TEC variations in November 2017. a From top to bottom, the solar radio flux at 10.7 cm of F10.7, and the magnetic indices of AE, Kp, and Dst. The Dst index displays a moderate (Storm 1) and a small (Storm 2) magnetic storms with a maximum depression − 70 and − 43 nT (Kp 6+ and 50; the maximum AE 945 and 691) on 7 and 21 November 2017, respectively. The AE maximum 596 nT, Kp 3+, and Dst -24 nT, which show relatively magnetic quiet, on 2–4 November 2017. b A latitude-time-TEC plot along the epicenter longitude extracted from GIMs in November 2017. The white and black lines denote the epicenter latitude and the magnetic equator, respectively. c A time series of the GPS TEC over the epicenter (34.9°N, 46.0°E, around 35°N, 45.0°E) was extracted from GIMs in November 2017. The red, gray, and black curves demote the observation, median, and upper/lower bounds, respectively. The 12 November 2017 M7.3 Iran–Iraq Border earthquake is denoted by the vertical red solid line. The red/black dotted curves denote the deviation between the observed values and the computed median (∆TEC). The red/black shaded areas indicate the increase/decrease anomalous strength, which the observed GPS TECs exceed the associated upper/lower bound. The red, blue, and light blue dashed rectangular denote the study periods of the M7.3 earthquake PEIAs, Storm 1, and Storm 2, respectively

Full size image

Based on the TECs in the Iran–Iraq area during 1999–2016, the chance of observing a new TEC in the interval (LB, UB) is about 60% and 85% for k = 1.5 and k = 3.0, respectively. We focus on the deviation between the observed values and the computed median,

$$\Delta {\text{TEC}} = O \, - \tilde{M}.$$
(2)

Therefore, when an observed TEC is not in the associated (LB, UB), which is ∆TEC > k (UQ − \(\widetilde{M}\)) or ∆TEC <  − k (\(\widetilde{M}\)− LQ), we declare an upper (increase; positive) or lower (decrease; negative) abnormal TEC signal. For the M7.3 earthquake, we set a stringent criterion of k = 3.0. No significant negative TEC anomaly occurs at 34.9°N, 46.0°E in November 2017. However, as shown in Fig. 2c, positive TEC anomalies prominently appear in the three time periods: pre-midnight on 2 to post-midnight on 4 (i.e., 2–4 November) November (10–8 days before the earthquake), 7 November (Storm 1), and 21–22 November (Storm 2). Since it is relatively magnetic quiet, we would like to see whether the positive anomalies on 2–4 November are PEIAs of the M7.3 earthquake.

To see whether the anomalies detected above are possibly the temporal PEIA, the characteristics of the polarity (i.e., negative or positive), appearance local time, duration, lead day, etc., of anomalies associated with previous large earthquakes in the Iran–Iraq area have to be reached. If the observed anomalies well meet the characteristics, we then can consider the temporal PEIA being detected. To find the characteristics, we consider the GIM TEC and the associated 53 M ≥ 5.5 earthquakes in the Iran–Iraq area (25–40°N × 42.5–62.5°E), which are downloaded from USGS (https://earthquake.usgs.gov/earthquakes/browse/significant.php?year=2022) during 1999–2016 (Fig. 3 and Table 1). The one sample test (Neter et al. 1988; Chen et al. 2015; Liu et al. 2018c) is then used to find if the increase or decrease TEC is statistically significant for the earthquakes under study. Let π be the observed proportion of earthquake-related anomalies and π0 the background proportion of anomalies in the 18-year period of 6574 (1999 to 2016) days. The z value is then given by

$$z = \frac{{\pi - \pi_{0} }}{{\sqrt {\pi_{0} \left( {1 - \pi_{0} } \right)/n} }} ,$$
(3)

where n = 53 is the number of earthquakes. If z > 1.96, we claim, at significant level 0.05, that π > π0. By contrast, if z < − 1.96, we claim, at significant level 0.05, that π < π0. Note that the one sample test is conducted for negative and positive anomalies separately. Figure 4 displays four positive anomaly zones with statistical significance with level < 0.05 in Zone A, 00:00–04:00 UT 14–12 days; Zone B, 04:00–08:00 UT 12–9 days; Zone C, 08:00–1200 UT 14–11 days; and Zone D, 03:00–05:00 UT 7–6 days before the 53 M ≥ 5.5 earthquakes. Overall, positive anomalies tend to appear 04:00–16:00 UT day 6–14 days before the M ≥ 5.5 earthquakes, which indicates that the positive TEC anomaly appearing during 2–4 November (10–8 days before) is possibly associated with the M7.3 Iran–Iraq earthquake.

Fig. 3
figure 3

Locations of 53 M ≥ 5.5 earthquakes in the Iran–Iraq area (25–40°N × 42.5–62.5°E) during 1999–2016. The red circles stand for the earthquakes. The blue triangle (35°N, 45°E) denotes the location of the TEC value extracted from the GIM for the statistical analysis of finding the PEIA characteristics. The solid red star and open red star denote the epicenter of the 12 November 2017 M7.3 Iran–Iraq Border earthquake and the epicenter of the 1 December 2017 M6.1 earthquake, respectively

Full size image
Table 1 53 M > 5.5 Earthquakes in Iran–Iraq border area during 1999–2016
Full size table
Fig. 4
figure 4

Median values of ∆TEC (= TEC observation − \(\tilde{M}\)) at the fixed location (35°N, 45°E) 30 days before and after the 53 M ≥ 5.5 earthquakes. ∆TEC is in a TEC unit (TECu = 1 × 1016 #/m2). Red and black contours denote one sample test results of positive and negative anomalies with significant level < 0.05, respectively. Zones A, B, C, and D are TEC positive anomalies with significant level < 0.05 at 00:00–04:00 UT day 14–12 before, at 04:00–08:00 UT day 12–9 before, at 08:00–12:00 UT day 14–11 before, and at 03:00–05:00 UT day 7–6 before the earthquakes, respectively

Full size image

To verify whether the PEIAs are candidates for a reliable precursor, we treat the positive TEC anomalies as alarms for earthquakes in Zones A, B, C, and D (Fig. 4) and construct the statistical analysis of receiver operating characteristic (ROC) curve (Swets 1988) to evaluate the reliability of the earthquake alarming. Taking Zone B as an example, based on its PEIA characteristic for the threshold k = 1.5, when positive TEC anomalies appear more than one-third of the period of 0400–0800 UT, we issue an alarm for an earthquake with the magnitude of M ≥ 5.5 occurring in the following 9–12 days. Note that k = 1.5 is simply an empirical threshold in previous studies (Chen et al. 2015; Liu et al. 2018c). For small k values, many alarming days are issued, and hence, more false alarms are obtained. However, for large k values, limited alarming days are issued and both the false alarm and successful rates are drastically reduced. Hence, to test the preference of the positive anomalies, we consider constructing the ROC curve for various k values. For each k value, we examine four different conditions, an alarm day being followed by earthquakes or no earthquake, and a non-alarm day being followed by earthquakes or no earthquake within a certain lead day period. Let TP(k) and FN(k) stand for numbers of earthquake days with and without being led by alarm days, respectively, while FP(k) and TN(k) denote numbers of non-earthquake days with and without being alarmed, respectively. Then, we have a 2 × 2 contingency table and yield the true positive rate TPR(k) and false positive rate FPR(k) as given by

$${\text{TPR}}\left( k \right) = \frac{{{\text{TP}}\left( k \right)}}{{{\text{TP}}\left( k \right) + {\text{FN}}\left( k \right)}}$$
(4a)

and

$${\text{FPR}}\left( k \right) = \frac{{{\text{FP}}\left( k \right)}}{{{\text{FP}}\left( k \right) + {\text{TN}}\left( k \right)}}$$
(4b)

where TPR(k) is the probability that an earthquake is successfully alarmed, and FPR(k) is the probability to make a false alarm. Note that the value of k varies from 0 to 10 by increasing 0.05, and hence, there are 201 k values under investigation. Hence, the ROC curve with FPR(k) as the x-axis and TPR(k) as the y-axis can be constructed. Figure 5 shows the ROC curves (color curves) of the positive TEC anomalies in the four zones. Note that the dash line in Fig. 5 represents the null ROC curve where TPR(k) = FPR(k) for all k, an equal chance to alarm earthquake day and non-earthquake day. It means that the PEIAs are actually independent of the occurrence of earthquakes. The red curve is the 95% upper confidence bound for the null ROC curve (Sarlis and Christopoulos 2014). Since most part of the observed ROC curve in each of the four zones is above the red curve, the positive TEC anomalies are statistically associated with the earthquakes under study.

Fig. 5
figure 5

ROC curves for alarming M ≥ 5.5 earthquakes based on PEIA information from the four time zones, A-D, as indicated in Fig. 4. Four panels are simulations of random of the interevent time of the PEIA days. The color, gray, and red curves denote the ROC curves of the observations, 1000 simulations, and the 95% interval (Sarlis and Christopoulos 2014), respectively. The blue and red asterisks denote k = 1.5 and the best point yielding the maximum R score (= TPR − FPR), which is called the Youden index (Youden 1950), respectively. p value of the four zones are “zero”

Full size image

The area under the ROC curve (AUC) is used to assess the effectiveness of the PEIA as a precursor (Hanley and McNeil 1982). Note that, under the null ROC curve, AUC = 0.5. Therefore, a reliable precursor should have AUC > 0.5. To investigate whether the PEIAs are significantly related to the earthquakes under study, we perform a simulation-based statistical test where the PEIAs are treated as a random process, while the occurrence times of earthquakes are retained. In the statistical test, the null hypothesis H0 (AUC = 0.5) states that the occurrences of PEIAs and earthquakes are independent, and the alternative hypothesis H1 (AUC > 0.5) says that the occurrences of PEIAs and earthquakes are associated. Again, taking Zone B as an example, we find that with k = 1.5, there are 1892 PEIA days, more than one-third of positive TEC anomalies in the period of 04:00–08:00 UT being detected, which results in 1891 (= 1892 − 1) inter-PEIA times during 1999–2016. We then take 1000 bootstrap samples (Efron 1979) of size 1891 from the inter-PEIA times to form PEIA appearances for earthquake alarming. The same method is then applied to each of 201 k values, and hence, 1000 possible ROC curves are obtained. Figure 5 reveals 1000 simulated ROC curves (gray curves) with the associated TPR(k)*, FPR(k)*, and AUC*. The simulation results allow us further compute the p value, which is the proportion of the simulated AUC* larger than the observed AUC. Therefore, small p values (< 0.05) lead to the rejection of the null hypothesis H0. The resulted p values are, in fact, all zero. Thus, the statistical results of the one sample test and ROC-AUC analysis show that in the Iran–Iraq border area, the PEIA characteristic is the TEC significant increases (positive anomalies) at 00:00–12:00 UT day 14–6 before M ≥ 5.5 earthquakes. In fact, the observed ROC curve being above the 95% upper confidence bound curve and the AUCs being all significantly greater than 0.5 strongly suggest that PEIAs appearing on 2–4 November 2017 are the M7.3 earthquake-related.

Among the 53 earthquakes, there are 12 thrust, 22 thrust-oblique, 14 strike-slip, and 5 normal/normal-oblique earthquakes. Table 1 displays date, time, location, depth, magnitude, PEIA on day 9 before the earthquake, and focal mechanism of the 53 earthquakes. We first find odds (= X/(n − X), the ratio of earthquakes with-to-without the PEIA, of the thrust, thrust-oblique, strike-slip, and normal/normal-oblique earthquakes being 3.0(= 9/3), 2.67(= 16/6), 2.50(= 10/4), and infinite(= 5/0), respectively (Fig. 6a). The normal earthquakes yield the greatest odds, while the strike-slip earthquakes have the smallest one. The value of odds 2.36 and 3.17 meets significance level 0.1 and 0.05, respectively (Conover 1999). Therefore, the odds of the 53 earthquakes being 3.07 suggests that M ≥ 5.5 earthquakes in the Iran–Iraq area are more likely leaded by the PEIAs, regardless the focal mechanism.

Fig. 6
figure 6

a Odds (= X/(n − X), the ratio of earthquakes with-to-without the PEIA) of the thrust, thrust-oblique, strike-slip, and normal/normal-oblique earthquakes. White and open bars denote the number of earthquakes with and without the associated PEIAs, respectively. b The depth and magnitude of the earthquakes as well as the odds of PEIA. Solid and empty circles represent the earthquakes with and without the associated PEIAs, respectively. c Odds (every 10 earthquakes sliding by 1 from shallow to deep depth) versus the related median depth. d Odds of M ≥ 5.5, 6.0, 6.5, and 7.0 earthquakes. The dotted-dashed and dashed lines denote the value of odds 2.36 (significance level 0.1) and 3.17 (significance level 0.05), respectively

Full size image

To further understand the relationship between the occurrence of PEIAs and the parameters of the related earthquakes, the depths of the earthquakes with/without the associated PEIAs are plotted against the corresponding magnitudes (Fig. 6b). To find the earthquakes with certain depths that are more likely to experience PEIA, the odds of every 10 earthquakes sliding by 1 with PEIA are computed from small to large depths. Figure 6c depicts that the earthquakes with depth D ≤ 12 km are more likely to experience the PEIAs. To find the relationship between PEIA occurrence and earthquake magnitude, we further compute the odds of earthquakes equal/greater than a certain magnitude. Figure 6d displays that the odds of M ≥ 5.5, 6.0, 6.5, and 7.0 earthquakes are 3.07(= 40/13), 3.75(= 15/4), 1.33(= 4/3), and infinite (= 3/0), respectively. This suggests that the larger earthquakes generally have the better chance being leaded by the PEIAs. Note that the three M ≥ 7.0 earthquakes all have been leaded by PEIAs.

Similar to the procedure in identifying anomalies of Fig. 2c, Fig. 7a–c examines the global distribution of the percentage of positive anomaly occurrences (i.e., frequency counts) at the 5183 lattices of GIM during a 31-h (time-point) period starting from 22:00UT on 2 November and ending by 04:00UT on 4 November 2017 (i.e., pre-midnight of day 10 to post-midnight of day 8 before the earthquake), a 24-h period of Storm 1 on 7 November, and a 29-h period of Storm 2 on 21–22 November. Figure 7d, e, and f illustrates the distributions of the percentage exceeding 50% in Fig. 7a–c to have a better viewing on the distributions of positive anomalies associated with the PEIA, Storm 1 and Storm 2, respectively. Figure 7g displays that 12 lattices with significant TEC increases are mainly around the epicenter area, 0.23% (= 12/5183) of GIM, and frequently occur more than 27 time points (87% = 27/31) of the 31-h-time-point period (also see Fig. 7a and d). The significant TEC increases of the 12 lattices frequently occurring specifically around the epicenter as well as the agreement between the significant increase anomalies on 2–4 November and the characteristics confirm that the PEIAs in the GIM TEC associated with the M7.3 earthquake have been observed. Figure 7e displays that the significant TEC increases frequently (greater than 71% = 17/24) appear at 60–80°N, all longitudes, except 20–50°W; at 30–60°N, 80°W–15°E; and at 45–55°S, 75–140°E, which shows the positive storm signature of Storm 1 on 7 November. Similarly, Fig. 7f illustrates that the significant TEC increases frequently (72% = 21/29) appear at 50°N, 70°W; and 15–55°N, 0–80°E in the Northern hemisphere, as well as at 30–60°S, 80°W–115°W and 135°E–165°E in the Southern hemisphere. Note that the most frequent (100% = 29/29) TEC increases occurred inside the rectangular area during the 29-time-point period on 21–22 November 2017 (Fig. 7h). These indicate that positive storm signatures of Storm 2 at mid- and high latitudes have been mixed with some possible PEIAs. Figure 7h further reveals that the most frequent TEC increases with 100% specifically appear inside the rectangular area of 27 lattices, 0.52% (= 27/5183) of GIM, which strongly suggests the occurrence of forthcoming large earthquakes around the area. Again, based on the characteristic, the PEIA might precede forthcoming larger earthquakes by 6–14 days. It is surprising to find in the catalog published by USGS (https://earthquake.usgs.gov/earthquakes/eventpage/us1000bjnz/executive) report that a magnitude M6.1 earthquake (30.7°N, 57.3°E) with a depth of 9 km struck at 02:32 UTC on 1 December 2017, which is 9–10 days after the significant TEC increases on 21–22 November 2017. Nevertheless, the spatial analysis on significant TEC increases discriminates local effects of the PEIAs of the M6.1 earthquake and global effects of the positive storm signatures of Storm 2 during 21–22 November 2017.

Fig. 7
figure 7

Distributions of positive GIM TEC anomaly occurrence percentages during periods of a, d, and f, the M7.3 PEIAs from 22:00UT on 2 November to 04:00UT on 4 November 2017; b and e Storm 1 from 00:00UT on 7 November to 23:00UT on 7 November 2017; and c, f, and h Storm 2 from 00:00UT on 21 November to 04:00UT on 22 November 2017. a–c without any percentage threshold, d–f with the percentage threshold of 50%, g the top occurrences percentage of 77%, and h the top occurrences percentage of 100% The solid and open red stars denote the epicenter of the 12 November 2017 M7.3 earthquake and 1 December 2017 M6.1 earthquake, respectively

Full size image

The spatial coverage of GIM TEC with 5183 (71 × 73) lattices (locations) allows us further computing odds of positive anomalies on the globe. Following Fig. 7a, c, Fig. 8a, d illustrates those odds of positive anomalies on the globe during the 31-time-point M7.3 PEIA period (from 22:00UT on 2 November to 04:00UT on 4 November 2017) and during the 29-time-point Storm 2 plus M6.1 PEIA period (from 00:00UT on 21 November to 04:00UT on 22 November), respectively. The odds in 48 lattices around the M7.3 and M6.1 epicenters are 2.10–6.75 and 8.66–29.00, and however, the odds become very small with means about 0.12 and 0.46 on lattices far from the epicenters, respectively. Table 1 shows that no M > 5.5 earthquake occurs in the Iran–Iraq Border area during November 2014–November 2017. To contrast the globe odds between earthquakes and non-earthquakes, we also compute the globe odds of positive anomalies for non-earthquake time periods by random selecting 100 sets of 31 (29) time points during December 2014–October 2017. Figure 8b (8e) displays means of 100 odds on the globe and the overall mean of about 0.15 (0.19) for non-earthquakes, which demonstrate that odds are very small globally during non-earthquake periods. Figure 8c (8f) shows that the ratios of odds in Fig. 8a and b (8d and 8e) are 16.79–63.07 (27.37–156.25) in 48 (57) lattices around epicenter of the M7.3 (M6.1) earthquake, but the ratios become very small with the mean about 2.64 (3.54) on lattices far from the epicenter. The ratios of odds around the epicenters are significantly larger than those far from, which confirm that PEIAs are associated with the M7.3 and M6.1 earthquakes. White lattices in Fig. 8c and f denote ratios of odds being smaller than one, which indicates that positive anomalies appear less frequent during the M7.3 PEIA and Storm 2 plus M6.1 PEIA periods than their associated bases. In Fig. 8c, white lattices appearing in many large areas on the globe show that the magnetic condition during the M7.3 PEIA period is quieter than during its base period. By contrast, in Fig. 8f, much fewer white lattices appear on the globe, especially at mid-/high latitudes, which shows that the magnetic condition is rather disturbed, and Storm 2 effects have been observed.

Fig. 8
figure 8

Odds, odds bases constructed by 100 random simulations, and odds ratios during the M7.3 PEIA and Storm 2 plus M6.1 PEIA periods. a Odds in the M7.3 PEIA period, b Mean odds of the M7.3 reference/base, c Odds ratio of M7.3 earthquake PEIAs, (d) Odds in the M6.1 PEIA plus Storm 2 period, e Mean odds of M6.1 reference/base, and f Odds ratios of Storm 2 plus M6.1 earthquake PEIAs. White lattices denote ratios of odds being smaller than one

Full size image

For a cross-comparison, the GIM TEC observations along the same F5 orbits (i.e., co-located) and at about the same measurement time (i.e., concurrent) are extracted. We focus on the GIM TEC and F5/AIP data within the studied rectangular area (Fig. 1). Figure 9 displays observations and associated references constructed by the moving median 7 days before and after the observation day. From top to bottom, the observation and associated reference of GIM TEC, the F5/AIP ion density, ion temperature, downward velocity, and eastward velocity are illustrated, respectively. The reason why the moving window of 7 days before and after the observation day is used to construct the reference is that the F5/AIP starts measuring the ionospheric plasma since the end of October of 2017. Figure 9 illustrates that the ion density and the ion temperature are greater than their associated reference, while the downward velocity and eastward velocity, respectively, decrease and increase, during the PEIA of 2–3 November as well as the two storm days of 7 and 21–22 November, respectively. Due to the quasi-neutrality, the observations and references of GIM TEC are generally very similar to those of the ion density. This indicates that the plasma quantities probed by F5/AIP can also be used to study ionospheric disturbances.

Fig. 9
figure 9

Observations and References of GIM TEC and F5/AIP ion data at 22:30 LT in November 2017. The reference is the moving median 7 days before and after the observation day inside the rectangular area. From top to bottom, the observation and associated reference of the GIM TEC, F5/AIP ion density, ion temperature, downward velocity, and eastward velocity are illustrated, respectively. The red, blue, and light blue dashed lines denote the PEIAs on 2–3 November 2017, Storm 1 on 7 November 2017, and Storm 2 on 21–22 November 2017, respectively

Full size image

Ionospheric data are positive values, which inhabit a right-skewed and heavy-tailed distribution, and therefore, it is suitable to apply a median-based analysis. Note that the box-and-whisker procedure (Wilcox 2010), as a median-based analysis, has the advantage of visually observing the significant difference among multi-datasets simultaneously. Therefore, we employ box-and-whisker (box) plots (Fig. 10) to investigate the ion density, ion temperature, ion downward velocity, and ion eastward velocity anomalies inside the rectangular area on the observation and reference days shown in Fig. 9 during the PEIA of 2–3 November and two storm days of 7 and 21–22 November. The ends of the box in Fig. 10 are the upper and lower quartiles, where the lower (upper) quartile is the number such that at least 25% of observations are less (greater) than or equal to it. The horizontal line within the box denotes the median. If two boxes do not overlap with each other, we consider that there is a dramatic difference between the two boxes. However, when one shorter box with median is larger than the upper quartile or smaller than the lower quartiles of the other longer box, the two boxes might still be considerd to be different. Therefore, the observation and reference days might have different plasma parameter (or quantity) values when the two boxes are not overlapped.

Fig. 10
figure 10

Box plots of the observation (left-hand side color box) and the associated reference (right-hand side gray box) of the ion density (top row), ion temperature (second row), downward velocity (third row), and eastward velocity (bottom row) during the M7.3 earthquake PEIA (left column), Storm 1 (central column) and Storm 2 (right column). The horizontal line within the box denotes the median. The ends of the box are the first quartile (25% of the dataset, Q1) and third quartile (75% of the dataset, Q3), where the first (third) quartile is the middle value between the smallest (highest) and the median of the dataset. The difference between the Q1 and Q3 is called the inter-quartile range (IQR). If a value lower than Q1 − 1.5IQR (lower dot) and/or greater than Q3 + 1.5IQR (upper dot), it is declared the outlier (cross). The horizontal lines extending out from the box are the minimum and maximum value, which the minimum (maximum) is the smallest (largest) value within the range of outlier. The vertical lines out from the box to the minimum (maximum) are called the lower (upper) whiskers

Full size image

To have a more stringent investigation, we employ the Mann–Whitney U test (Corder and Foreman 2014) as a nonparametric test for possibly different plasma parameter values on the observation and reference days since the plasma parameter values may not be distributed according to the normal distribution. Let \({X}_{1}\), \({X}_{2}\),…, \({X}_{m}\) be the reference values and \({Y}_{1}\), \({Y}_{2}\),…, \({Y}_{n}\) be the observed values. Define I{\({Y}_{j}>{X}_{i}\)} = 1 if \({Y}_{j}>{X}_{i}\), = 0, otherwise for \(i=\mathrm{1,2},\dots ,m\) and \(j=\mathrm{1,2},\dots ,n.\) Then, the U statistic in the Mann–Whitney test is

$$U = \frac{{\mathop \sum \nolimits_{i = 1}^{m} \mathop \sum \nolimits_{j = 1}^{n} {\text{I}}\left\{ {Y_{j} > X_{i} } \right\} - mn/2}}{{\sqrt {mn\left( {m + n + 1} \right)/12} }},$$
(5)

and if the observation and reference days share the same plasma parameter value, the distribution of U can be well approximated by the standard normal distribution. Hence, under significance level of 0.05, we claim that the plasma parameter value during the observation days is larger (smaller) than that during the reference days if U > 1.96 (U < − 1.96). Results of the Mann–Whitney U test show that p values are all about of zero, except p = 0.60 for F5/AIP Ti during the Storm 1 period. These confirm that F5/AIP plasma parameters of the observation and the associated reference being significantly different, except Ti during the Storm 1 period. Therefore, the F5/AIP ion velocity can be used to derive the electric fields associated with the PEIA, Storm 1, and Storm 2.

Figure 10 shows that the two boxes of the observation (left-hand side box) and associated reference (right-hand side box) of the F5/AIP plasma quantities are either different or dramatically different. The top two rows of Fig. 10 show that the differences of the median values in the ion density (ion temperature) increase by 0.4 × 105, 0.6 × 105, and 0.4 × 105 #/cm3 (63, 128, and 40°K) during the PEIA, Storm 1, and Storm 2, respectively. Note that, again, Storm 2 convolved with the PEIAs of the M6.1 earthquake yields greater increases in the ion density and temperature, which well agrees with Figs. 2 and 9. The bottom two rows illustrate that the median values in the downward (eastward) velocity significantly decrease (increase) during the PEIA, Storm 1, and Storm 2 days. Based on the dynamo theory (cf. Kelley 2009), the electric field E can be expressed as,

$${\mathbf{E}} = - {\mathbf{V}} \times {\mathbf{B}},$$
(6)

where V is the ion velocity and the B is the Earth’s magnetic field. From the IGRF (International Geomagnetic Reference Field, https://wdc.kugi.kyoto-u.ac.jp/igrf/point/index.html) model, we find the B field at the satellite orbit height of 720 km altitude over the epicenter is 3.9 × 10–5 T with the magnetic dip of 68.30 degrees and the declination of 7.45° over the Iran–Iraq border area. By inserting the median value of the velocities of the two boxes into Eq. (1), we can calculate the electric fields on the observation and reference days, as well as their difference, subtracting the former from the latter. The bottom two rows in Fig. 10 show that the eastward (downward) electric field of 0.3, 1.2, and 1.0 mV/m (0.8, 0.5, and 0.8 mV/m) is generated during the PEIA, Storm 1, and Storm 2 period, respectively. Note that due to the mixture effect of Storm 2 and PEIA associated with the 1 December 2017 M6.1 earthquake, the eastward and downward electric fields are greater than those of Storm 1, respectively.

3 Discussion

Retrospective studies on the previous 53 M ≥ 5.5 earthquakes in the Iran–Iraq region during 1999–2016 have been conducted to find the characteristics of the temporal PEIAs (Figs. 3 and 4). Figure 4 shows that in the Iran–Iraq area, the PEIA characteristics have significant level 0.05 in one sample test that is the GIM TEC significant increases (positive anomalies) at 04:00–16:00 UT (07:00–19:00 LT (local time)) day 14–6 before M ≥ 5.5 earthquakes. The TEC significantly anomalously increases on day 10–8 before (2–4 November 2017) the 12 November 2017 M7.3 and on day 10–9 before (21–22 November) the 1 December 2017 M6.1 earthquake (Fig. 2c), well meeting the characteristics of positive anomalies appearing on day 14–6 before M ≥ 5.5 earthquakes reached by the one sample test (Fig. 4) and validated by the statistical results of the ROC-AUC analysis (Fig. 5). Therefore, we declare that the temporal PEIAs related to the two earthquakes have been detected.

By applying similar processes to Fig. 2c, the spatial analysis of the global 5183-lattice search is further used to confirm the temporal PEIAs being detected. Figure 7d and g shows that the positive TEC anomalies frequently, 87% (= 27/31) of the 31-h-time-point and 100% (= 28/28) of the 28-h-time-point period, and specifically, 0.23% (= 12/5183) and 0.52% (= 27/5183) of 5183 lattices on GIM, appear over the epicenters, which confirm that PEIAs of the 12 November 2017 M7.3 and the 1 December 2017 6.1 earthquake have been observed.

Figure 4 shows that PEIAs tend to appear at the same times, 04:00–16:00 UT (07:00–19:00 LT) day 6–14 before M ≥ 5.5 earthquakes, independent of the rupture time, which is rather random. The discrepancy might result from the rupture being seismological or mechanical processes, while PEIAs are related to the electromagnetic processes. Scientists observe traveling ionospheric disturbances almost right after large earthquakes or tsunamis, which induce by their ruptures via mechanical mechanisms of vertical motions of the Earth’s surface (e.g., Liu et al. 2006b, 2006c, 2010d, 2011b, 2012, 2016b, 2019, 2020, Liu and Sun 2011). By contrast, around the epicenter during the earthquake preparation period, electric fields near Earth’s atmosphere can be generated by underground seismo-electromagnetic processes, which are further mapped along the geomagnetic field line into the ionosphere and result in PEIAs. The mapping efficiency is mainly a function of the ionospheric conductivity, which yields the diurnal variation significantly (cf. Kelley 2009). Consequently, PEIAs generally appear at the same time. On the other hand, the characteristics of polarity, duration, and lead day might be related to underground structures, focal mechanisms, etc.

Ratcliffe (1972) and Kelley (2009) find that a stronger eastward electric field of the daily dynamo results in the EIA crest moving poleward. Liu et al. (2010c) examine the GPS TEC and M ≥ 5.0 earthquakes in Taiwan during 2001–2007, and find that the PEIA-associated electric fields can strongly perturb the daily dynamo electric fields and affect the EIA crest location few days before the earthquakes. On the other hand, the prompt penetration electric field can also superimpose with the daily dynamo electric field and affect the EIA crest location. The poleward motions of the EIA crest on 2–4, 7, and 21–22 November shown in Fig. 2b indicate that the PEIA-associated electric fields of the M7.3 earthquake and the prompt penetration electric fields are in the eastward direction.

For the temporal analyses, the significant TEC increases over the epicenter area on 2–4 November and on 21–22 November 2017 shown in Fig. 2c agree with the PEIA characteristics in the Iran–Iraq border area in Fig. 4, which indicates the temporal PEIAs of the two earthquakes being detected. Regarding the ionospheric storm, no significant decrease (negative) in TEC anomalies has been detected in November 2017. This suggests that the wind disturbance dynamo (Blanc and Richmond 1980; Lin et al. 2005; Kelley 2009; Fuller-Rowell 2011; Liu et al. 2013a, b) of the two storms are not prominent. By contrast, the positive storm signatures of significant TEC increase on 7 and 21–22 November 2017 result from the prominent prompt penetration electric fields in the eastward (Jaggi and Wolf 1973; Kelley 2009; Fuller-Rowell 2011; Liu et al. 2013a, b) of the two storms. Although Storm 1 is greater than Storm 2, the positive storm signatures of Storm 2 confounded by PEIAs of the M6.1 earthquake yield the greater TEC increase and the longer duration (Fig. 2).

Figure 6a shows that the normal earthquakes yield the greatest odds of infinite, and the strike-slip earthquakes have the smallest one of 2.50. The odds of the overall 53 earthquakes being 3.07 meets the significance level 0.1, which strongly suggests that M ≥ 5.5 earthquakes in the Iran–Iraq border area are more likely led by the PEIAs, regardless of the focal mechanism. Meanwhile, the odds studies show that the shallow earthquakes with depth D ≤ 12 km are more likely to experience the PEIAs (Fig. 6c), and the larger earthquakes generally have the better chance of being preceded by the PEIAs (Fig. 6d).

For the spatial analyses, the significant TEC increases frequently appearing specifically over the epicenter during 2–4 November confirm that PEIAs of the GIM TEC associated with the 2017 M7.3 Iran–Iraq border earthquake have been observed (Fig. 7a and d). The significant TEC increases frequently occur at worldwide mid- and high latitudes on 7 and 21–22 November 2017, which confirms that the positive storm signatures of Storm 1 and 2 have been detected (Fig. 7b–f). Figure 7c and f shows that in addition to the significant TEC increases at mid- and high latitudes, the most intense TEC increases appear inside the Iran–Iraq border area during 21–22 November 2017 (Fig. 7g), which was struck by the M6.1 earthquake (30.7°N, 57.3°E) on 1 December 2017. Therefore, the significant TEC increases on 21–22 November 2017 are the superposition of the positive storm signatures of Storm 2 and PEIAs on day 10–9 before the M6.1 earthquake. This explains that the TEC increase strength and duration of a small magnetic storm as Storm 2 are, respectively, greater and longer than those of a moderate one as Storm 1 (Fig. 2c).

Figure 8 displays the odds, odds bases constructed by 100 sets of random simulations, and ratios of odds during the M7.3 PEIA and Storm 2 plus M6.1 PEIA periods. Odds of about 0.46 far from the M6.1 epicenter are larger than those of 0.12 far from the M7.3 epicenter, which suggests that Storm 2 has the amplification factor of 3.83 (= 0.46/0.12). Figure 8c and 8f shows around the epicenters that the odds ratios of 16.79–63.07 during the M7.3 PEIA period are smaller those of 27.37–156.25 during the Storm 2 plus M6.1 PEIA period, respectively. Taking the amplification factor of 3.83 into consideration, the odds ratios of the M6.1 PEIA should be calibrated as 7.14–40.80, which are smaller than those of the M7.3 PEIA. This suggests that magnetic storms could affect PEIA occurrences, and larger earthquakes tend to experience more PEIAs. The 100 random simulations yield very small values of odds on the globe (Fig. 8b and e), while very large odds and odd ratios appear specifically in the small area of about 50 out of 5183 lattices around the epicenter (Fig. 8a, c, d, and f). These again show that PEIAs associated with the M7.3 and M6.1 earthquakes have been observed.

Figure 9 depicts the good agreements in observations and references between GIM TEC and F5/AIP ion density, which shows that F5/AIP can be useful to study PEIAs and ionospheric storms. The GIM TEC, F5/AIP ion density, ion temperature, and ion eastward (downward) velocity significantly increase (decreases) during the PEIA, Storm 1, and Storm 2 days. The box plots for the F5/AIP observations and references in Fig. 10 show that these increases during the PEIA and two storm periods are in difference and dramatic differences, respectively. The ion density and the ion temperature concurrently increasing during the three periods show that the cooling through Coulomb collisions (Kakinami et al. 2011) does not occur and, however, strongly suggests that some external mechanisms/forces might involve. The poleward motion of the EIA crests in Fig. 2b implies that the eastward electric fields on the PEIA, Storm 1, and Storm 2 days have been enhanced. In the third row of Fig. 10, the eastward electric field of 0.3 mV/m causes upward/northward E × B/B2 drift and results in the M7.3 PEIAs of the significant TEC and/or ion density increases, especially at the northward side of the epicenter (Fig. 7a), during the PEIA days of the 12 November 2017 M7.3 Iran–Iraq border Earthquake. Mozer and Serlin (1969) reported that the atmospheric field can be mapped without attenuation along the same magnetic field in the atmosphere, ionosphere, and magnetosphere. On the other hand, the eastward electric fields of 1.2 and 1.0 mV/m generally result from the prompt penetration electric fields of Storm 1 and Storm 2, respectively. Since the 1.0 mV/m electric field is contributed by Storm 2 and PEIAs of the M6.1 earthquake, the eastward electric field of either one of them should be smaller. Taking the amplification factor of 3.83 estimated by Fig. 8a and d into consideration, we find that the eastward electric fields of the M6.1 PEIAs would about 0.3 (= 1.0/3.83) mV/m, which is similar to that of the M7.3 PEIAs. This similarity implies that magnetic storms might affect PEIA-related electric fields. Likewise, the prompt penetration electric fields of Storm 2 would be about 0.7 (= 1.0–0.3) mV/m eastward. Nevertheless, the poleward motion of EIA crests in the GIM TEC shown in Fig. 2b and the upward motion of F5/AIP ion velocity confirm the eastward electric fields appearing on the PEIA, Storm 1, and Storm 2 days. In the bottom row, the downward electric field of 0.8 mV/m is related to the M7.3 earthquake PEIA, while those of 0.5 and 0.8 mV/m are due to the downward flow of Region 2 currents around pre-midnight (Kelley 2009) during Storm 1 and Storm 2 days, respectively. However, due to the high inclination angle of 68 degrees and the field aligned currents in Region 2, vertical electric fields might be rather difficult to be estimated correctly.

Akhoondzadeh et al. (2019) conducted Swarm satellites (Alpha, Bravo and Charlie) data analysis inside the Dobrovolsky area around the M7.3 Iran earthquake epicenter during the period from 1 August to 30 November 2017. They found that six Swarm measured parameters including electron density, electron temperature, and magnetic scalar and three vector components reveal irregular variations between 8 and 11 days prior to the earthquake, which generally agree with the result that PEIAs of GIM TEC and the F5/AIP ion density, ion temperature, and ion downward/eastward velocity (i.e., the eastward/downward electric field) appear day 9–8 before the M7.3 earthquake (Figs. 2, 7, 8, 9 and 10).

4 Conclusions

Six different statistical analyses of the quartile-based process, one sample test, spatial analyses, odds, box plot, and Mann–Whitney U test have been used to rigorously identify PEIAs and storm signatures in the ground-based remote sensing ionospheric GIM TECs and in site F5/AIP plasma quantities. The significant TEC increases appearing day 9–8 before the 12 November M7.3 earthquake and day 10–9 before the 1 December 2017 M6.1 earthquake agree well with the characteristics in the Iran–Iraq border area, which indicate that the temporal PEIAs of the two earthquakes have been observed. The spatial analyses together with odds studies show that the significant TEC increases frequently occur specifically over a small area (less 1% (= 50/5183) of GIM, the globe) of the two epicenters, which confirms that the PEIAs of the two earthquakes have been observed. The significant TEC increases in the high-latitude ionosphere are the positive storm signatures of Storm 1 and 2, which indicates that the penetration electric field is essential. The spatial analyses can be employed to discriminate local effects of earthquakes from the global ones of magnetic storms, etc. Similar tendencies in concurrent and co-located measurements of the GIM TECs and the F5/AIP ion density indicate that the two observations can be used to three-dimensionally detect PEIAs and to examine ionospheric storm signatures. F5/AIP ion density, ion temperature, and especially ion velocity can be employed to study PEIAs and ionospheric storms. In conclusion, the ion velocity leads having a better understanding of causal mechanisms of ionospheric disturbances. The M7.3 PEIA-associated electric field of 0.3 mV/m eastward and the prompt penetration electric field of 1.0–1.2 mV/m eastward for the first time are simultaneously estimated. This suggests that the ionospheric weather can be modulated by electric fields from above, from the magnetosphere/space, and from below, the atmosphere/lithosphere.