Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 1137177 10.3389/feart.2023.1137177 Earth Science Original Research Statistical analysis of ionospheric vertical total electron content anomalies before global Mw≥6.0 shallow earthquakes during 2000–2020 Ma et al. 10.3389/feart.2023.1137177 Ma Yifang 1 Zhou Xiaohui 2 * Yang Yilin 3 Hu Leyin 1 Dong Hongyan 1 Yan Rui 1 1 Beijing Earthquake Agency, Beijing, China 2 School of Geodesy and Geomatics, Wuhan University, Wuhan, China 3 Nordic Volcanological Center, Institute of Earth Sciences, University of Iceland, Reykjavík, Iceland

Edited by: Giovanni Martinelli, National Institute of Geophysics and Volcanology, Italy

Reviewed by: Patricl Taylor, National Aeronautics and Space Administration (NASA), United States

Angelo De Santis, Istituto Nazionale di Geofisica e Vulcanologia (INGV), Italy

Funchun Chen, Chinese Academy of Sciences (CAS), China

 *Correspondence: Xiaohui Zhou, xhzhou@sgg.whu.edu.cn 14 04 2023 2023 11 1137177 04 01 2023 03 04 2023 Copyright © 2023 Ma, Zhou, Yang, Hu, Dong and Yan. 2023 Ma, Zhou, Yang, Hu, Dong and Yan

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

To quantitatively investigate the relationship between earthquakes and ionospheric anomalies, this paper presents a statistical study of pre-earthquake vertical total electron content (VTEC) variations. A total of 1522 shallow (≤60 km) strong (Mw≥6.0) earthquakes in the global area during 2000-2020 are selected, and classified according to different magnitudes, latitudes and focal depths. A quartile-based process with different lengths of sliding windows, equaling 10 days, 15 days and 27 days, respectively, has been utilized to detect VTEC anomalies. The abnormal level is first defined, and then VTEC anomalies occurrence probabilities ( P o ) and occurrence rates ( P E ) within 1-10 days before 1522 earthquakes have been calculated. Besides, VTEC anomalies occurrence rates of the background days (PN) are also calculated. The results show that the significant correlation between Po and epicentral latitudinal locations could be observed within 1-10 days before earthquakes. The values of Po increase with larger magnitudes in the equatorial and low-latitude regions, but decrease with greater magnitudes in the mid- and high-latitude regions to some degree. Within 1-5 days before earthquakes, the overall trend of P E shows an increase with larger magnitudes, but the correlation between the values of P E and magnitudes is relatively weak in the southern mid- and high-latitude regions. There is no evident causality between PN and the magnitude, and most of the values of P E / P N are larger than 1, indicating that VTEC anomalies within a few days before earthquakes are probably related with the forthcoming earthquakes. Moreover, when the abnormal level exceeds 60%, different sliding window lengths have a significant impact on the values of P o and P E in the mid- and high-latitude regions. In particular, there are obvious systematic deviations between the values of P o obtained from different sliding windows in the southern mid- and high-latitude regions. However, the selection of the optimal sliding window needs to be further studied.

 shallow strong earthquakes seismo-ionospheric anomaly VTEC statistical analysis different lengths of sliding window 1 Introduction

The seismic process is not only confined to the lithosphere, but also has impacts on the troposphere, ionosphere and even magnetosphere through the electromagnetic fields effect. The ionospheric anomalies within a few days before the earthquakes are relatively stable at short time scales, and have been studied widely in the field of earthquake prediction (Pulinets and Boyarchuk, 2004; Liperovsky et al., 2008; Pulinets and Ouzounov, 2011). A large number of studies have shown that ionospheric perturbations before many earthquakes could be identified (e.g., Liu et al., 2001; Le et al., 2011; Zhu et al., 2014; Tang et al., 2015; Sun et al., 2016; Parrot and Li, 2018; Pulinets et al., 2021). The possible earthquake-related NmF2 (F2 layer peak electron density) and TEC anomalies have been widely discussed in recent decades (e.g., Nishihashi et al., 2009; Ma et al., 2014). Especially, because of the development of the Global Navigation Satellite system (GNSS), GNSS VTEC have attracted more and more attention in the investigation of the ionospheric variations prior to large earthquakes (Liu et al., 2004; Shah and Jin, 2015). For the first time, Liu et al. (2001) used GNSS(GPS) VTEC to study the ionospheric disturbance before the Chi-Chi earthquake, and found that the VTEC over the epicenter decreased significantly 1, 3 and 4 days before the earthquake. Later on, more and more scientists began to focus on GNSS VTEC variations before earthquakes, aiming to detect the potential ionospheric anomalies related to the forthcoming earthquakes. For example, Yao et al. (2012) analyzed ionospheric variations prior to the 2011 Mw9.0 Japan earthquake, and indicated that ionospheric anomalies occurring on 8 March might be a precursor of the earthquake; Ho et al. (2013) showed that TEC increased 9–19 days before the 2010 M8.8 Chile earthquake and specifically over the epicenter; Su et al. (2013) investigated ionospheric TEC variations before the Hector Mine earthquake, and found that ionospheric disturbance appeared just above the epicenter 5 days before the earthquake. These studies show that the GNSS TEC anomalies appear a few days before the earthquake with different magnitude and focal depth.

Over about 50 years of research, no consensus in the scientific community has been formed on the existence of ionospheric earthquake precursors (Rishbeth, 2006; Dautermann et al., 2007; Masci, 2012; Ovalle et al., 2013; Masci and Thomas, 2014; Zolotov et al., 2019). Dautermann et al. (2007) indicated that there was no statistically significant correlation between TEC anomalies and earthquakes in Southern California during 2003–2004; Kon et al. (2011) selected M≥6.0 earthquakes in Japan during 1998–2010, and found that significant positive TEC anomalies within 1000 km above the epicenter appeared within 1–5 days before the earthquakes. According to Masci (2012), the analysis of Kon et al. (2011) was not reliable because of the influence of global geomagnetic events. Ovalle et al. (2013) concluded that it remained controversial whether the observed NmF2 and TEC anomalies were unambiguously related to the 2010 M8.8 Chile earthquake. Background geomagnetic events may impact revealing the relationship between TEC anomalies and earthquakes.

To validate the relationship between ionospheric anomalies and earthquakes in response to the controversy, many scientists have undertaken a multitude of studies on the physical mechanism of generating ionospheric anomalies. Firstly, the morphological characteristics of ionospheric anomalies before a large number of earthquakes are summarized. For instance, Liu et al. (2004) analyzed Ms5.0+ earthquakes in Taiwan from 1999 to 2002, and found that obvious negative TEC anomalies occurred 5 days before the earthquakes. Le et al. (2011) made a statistical study of global 736 M≥6.0 earthquakes during 2002–2010, and proposed that occurrence rates of abnormal days are larger for earthquakes with greater magnitude and lower depth. De Santis et al. (2019) analyzed the electron density and magnetic field data from 3 Swarm satellites to detect possible anomalies associated with 1312 M≥5.5 shallow earthquakes from January 2014 to August 2018, and the results showed that anomalies occurred between a few days and 80 days before the earthquakes with larger peaks at around 10, 20 and 80 days, and supported the Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) with clear statistical significance. Shah et al. (2020) studied the ionospheric anomalies before the global Mw≥5.0 earthquakes from 1998 to 2019, and the results revealed that prominent ionospheric anomalies appeared within 5 days before and after the earthquakes. Based on the characteristics of ionospheric anomalies prior to a large number of earthquakes, the physical mechanisms of seismic LAIC have been extensively studied (Freund, 2011; Klimenko et al., 2012; Pulinets, 2012; Zolotov et al., 2012). For example, Freund (2011) proposed that positive holes released by stressed rocks are highly mobile and can reach the Earth’s surface, and then ionize the atmosphere and change the vertical electric field between the ground and the lower edge of the ionosphere. In addition, a large number of studies have shown that anomalous atmospheric electric field variations in the earthquake preparation zone are likely to be the main cause of ionospheric disturbance (Zhang et al., 2014; Jiang et al., 2017; Pulinets and Davidenko, 2018; Davidenko and Pulinets, 2019). For example, Namgaladze et al. (2009) proposed that vertical plasma motion in the ionospheric F2 region under the action of the zonal electric field is the main disturbance formation factor, and ionospheric anomalies before strong earthquakes at middle and low latitudes verified this mechanism. Liu et al. (2010) studied the crest of equatorial ionization anomaly (EIA) variations before 150 M≥5.0 earthquakes in Taiwan, and the results implied that the weak atmospheric electric field a few days before the earthquakes may cause the EIA crest anomalies.

Both the characteristics of ionospheric anomalies and the physical mechanisms of lithosphere-atmosphere-ionosphere coupling present diversity and complexity, and the influence factors may include magnitudes, focal depths, latitude and longitude of the epicenter, focal mechanisms, the weather, the season, solar and geomagnetic activity and so on. However, many studies focused on a single large earthquake, and there are relatively few statistical results based on a large number of earthquakes. The corresponding characteristics of the ionospheric anomalies are still not fully understood. In addition, a quartile-based process is the most common method for extracting ionospheric VTEC anomalies, but different authors use different numbers of days as the lengths of the sliding windows, such as, 10 days (e.g., Zhou et al., 2009; Zhu et al., 2014), 15 days (e.g., Liu et al., 2010; Ke et al., 2016; Liu and Xu, 2017), 27 days (e.g., Xu et al., 2011; Guo et al., 2015), etc. It should be noted that the effects of different sliding windows on the ionospheric anomalies are rarely studied (Zolotov et al., 2019). To solve the above problems, this study uses GIM VTEC to carry out a statistical analysis by studying the VTEC anomalies within 1–10 days before 1522 global shallow (≤60 km) Mw≥6.0 earthquakes during 2000–2020. The factors, including the magnitude, focal depth and the latitude of the epicenter are considered for all the earthquakes, and the effects of different sliding windows on ionospheric anomalies are investigated. Moreover, ionospheric anomalies during background days are also analyzed to compare with those prior to the earthquakes. This study aims at helping the research on the physical mechanisms of lithosphere-atmosphere-ionosphere coupling by summarizing the characteristics of ionospheric anomalies comprehensively, and determining whether different lengths of sliding windows affect ionospheric anomalies features.

 2 Data and method 2.1 Data source

The worldwide Mw≥6.0 earthquakes during 2000–2020 are selected to analyze ionospheric anomalies in this study. The data are retrieved from the Global Centroid Moment Tensor (CMT) Project (http://www.globalcmt.org/). The earthquakes selected in this study are declustered from aftershocks following the method of Michael (2011), and the earthquakes occurring at the similar location but with the short interval (<10 days) from the previous ones are also excluded to avoid possible confounded effects from adjacent earthquakes. Finally, 1522 shallow (≤60 km) earthquakes are selected, and Figure 1 illustrates epicenter locations of these earthquakes.

Locations of global 1522 Mw≥6.0 earthquakes during 2000–2020.

The GIM VTEC is derived using the observations from hundreds of global GNSS stations (Hernández-Pajares et al., 2009). The GIM covers ± 87 . 5 ∘ latitude and ± 180 ∘ longitude ranges with spatial resolutions of 2 . 5 ∘ and 5 ∘ , respectively, and the time interval of the VTEC is 2 h. For each earthquake, the cell including the epicenter was selected as the point to analyze ionospheric VTEC anomalies. According to Dobrovolsky, (1979), the radius of the M6.0 earthquake preparation zone is about 380 km, corresponding to 3 . 5 ∘ . Therefore, the spatial resolution of GIM VTEC is sufficient to extract ionospheric anomalies using the nearest grid to the epicenter.

The equatorial geomagnetic activity index (Dst) data provided by the World Data Center for Geomagnetism, Kyoto (https://wdc.kugi.kyoto-u.ac.jp/) are used to represent the geomagnetic activity.

 2.2 Statistical method

In this study, a quartile-based process is performed to detect ionospheric VTEC anomalies within 1–10 days prior to each earthquake. As the length of sliding window is limited by the seasonal variability of the ionosphere at longer timescales, 10, 15, and 27 days are chosen as the candidate lengths of sliding windows based on the previous studies. At each time point on any day, the median x ¯ is computed using the VTEC at the same time point within 10, 15, and 27 days before this day as the background value, respectively, and the associated inter-quartile range I Q R is also obtained to construct the upper or lower bound x ¯ ± 1.5 I Q R . If VTEC continuously exceeds the associated upper or lower bounds for at least 6 h during a day, this day would be considered as an anomalous day. Moreover, the abnormal level (AL) is defined as the percentage of the largest deviation from the median (Le et al., 2011). Ionospheric anomalies with AL<20% are regarded as the daily effects of solar activities (Le et al., 2011, Personal Communication), therefore, only ionospheric anomalies with AL>20% are analyzed in this study. In addition, the ionospheric anomalies with AL>40% and AL>60% (i.e., n × A L , n = 1,2,3 with A L = 20 % ) are also checked for the purpose of studying whether the characteristics of ionospheric anomalies are similar with different AL. If a day with Dst≤−40 nT or Dst≥40nT, this day and the following 3 days are excluded to avoid the interference of the magnetic disturbed activity. After removing the effects of daily solar activity and geomagnetic disturbance, ionospheric anomalies occurring on 1 day are recorded as the seismo-ionospheric anomalies.

After the analysis of seismo-ionospheric anomalies for each event, 1522 earthquakes are divided into three different latitudinal zones (as shown in Figure 1), and VTEC anomalies occurrence probabilities ( P o ) and occurrence rates ( P E ) will be investigated respectively. For each zone, the earthquakes are firstly classified by magnitudes in increments of 0.1 or by depths in increments of 20km, and then we calculate P o and P E of each group. P o can be computed as the ratio of the number of earthquakes with seismo-ionospheric anomalies and the total number of earthquakes (Fujiwara et al., 2004), as shown in Eq. 1. In this equation, N o A L and N o T o t a l are the number of earthquakes with seismo-ionospheric anomalies and the total number of all the earthquakes in each group, respectively. For example, in the group (Mw ≥ 6.5 in Zone A with AL>20%), P o is the number of earthquakes with seismo-ionospheric anomalies divided by the total number of earthquakes. P o = N o A L N o T o t a l

The occurrence rates for the n th earthquake P E n can be calculated as the ratio of the number of the seismo-ionospheric abnormal days and the total quiet days: N A L , T n / T − Δ S n , and P E is defined as the mean of P E n , as shown in Eq. 2, derived from Le et al. (2011). In this equation, K is the number of earthquakes in each group (for example, Mw ≥ 6.5 in Zone A with AL>20%); N A L , T n is the number of seismo-ionospheric days with different AL (for example, AL>20%) within the T days before the n th earthquake, and T = 1,2,3 , . . . , 10 ; Δ S n is the number of magnetic disturbed days during the 1–10 days before the n th earthquake; N A L , T n / T − Δ S n is the number of the seismo-ionospheric abnormal days divided by the number of the total quiet days before the n th earthquake. P E = 1 K ∑ n = 1 K N A L , T n T − Δ S n × 100 %

3 Results and discussions 3.1 Seismo-ionospheric anomalies occurrence probabilities

According to the method described above, we calculated P o for earthquakes with different magnitudes and depths. Figure 2 shows P o with AL > 20%, 40%, and 60% within 1–10 days before different magnitude earthquakes with the depth ≤20, ≤40, and ≤60 km, respectively. It can be seen from Figure 2 (left), there is no significant correlation between the values of P o and the magnitude for the earthquakes of 6.0≤Mw < 6.5 in all three zones. There are larger values of P o in Zone A for larger magnitude earthquakes of Mw≥6.5, but the values of P o decrease with the magnitudes increasing in Zone B and C for Mw≥6.5 earthquakes to some extent. The values of P o in Zone A are higher than those in the other zones for Mw≥6.7 earthquakes, and all the results in Zone A are larger than those in the other zones with AL>20% and AL>40%. Ionospheric enhancements in equatorial regions may be the main reason (Liu et al., 2010; Shah et al., 2020). However, the results in Zone C for 6.0≤Mw < 6.7 earthquakes are higher than those in the other two zones with AL>60%. It needs to note that for the earthquakes in the mid- and high-latitude regions, other magnetic indices are not considered except the Dst index. So we compared the variations of the Kp index and Dst index, and found that if Kp=4 is chosen as the threshold value, the impact of not considering Kp can be ignored. But for the earthquakes in the high-latitude regions (as shown in Table 1), the number of disturbed day based on the AE index (>500 nT) is larger than that based on Dst. That is, external magnetic fields contamination is still probably not excluded in this study for these 51 earthquakes in Zone B and Zone C. However, the most of values of P o in Zone A are still higher than those in other zones. It indicates that the latitude of the epicenter has a significant influence on the P o , that is, the results of the low-latitude and equatorial regions are higher than those of the mid- and high-latitude regions. To a certain extent, P o increases with the magnitude increasing in the low-latitude and equatorial region, while P o decreases with the magnitude increasing in the mid- and high-latitude region.

Seismo-ionospheric anomalies occurrence probabilities of AL>20%, 40%, and 60% with different magnitudes (A, C, E) and depths (B, D, F), respectively (15-day sliding window). (A, C) are the results in Zone A. (C, D) are the results in Zone B. (E, F) are the results in Zone C.

The difference of disturbed days based on Dst and AE before 51 Earthquakes at high latitudes.

No Date Location of epicenter Magnitude (Mw) Place Disturbed days based on Dst (Doy) Disturbed days based on AE > 500 nT (Doy)
Year Doy Latitude (°) Longitude (°)
1 2014 198 60.42 −140.31 6 Zone B 190–191;193
2 2010 120 60.45 −177.71 6.5 Zone B 112–114;119
3 2006 142 60.86 165.81 6.6 Zone B 132–134;138–141
4 2006 110 60.89 167.05 7.6 Zone B 100–109 100;103–108
5 2018 334 61.49 −150.02 7.1 Zone B 329;331
6 2020 9 62.27 171 6.4 Zone B 3;8
7 2002 307 63.23 −144.89 7.8 Zone B 297–306 297–304;306
8 2008 150 63.92 −21.17 6.3 Zone B 142–144;146;149
9 2000 169 63.99 −20.47 6.5 Zone B 160–165 159–160;162–168
10 2020 173 66.46 −18.72 6 Zone B
11 2013 45 67.65 142.51 6.7 Zone B 35;38–39;43–44
12 2008 174 67.71 141.43 6.1 Zone B 167–170 166–170;172
13 2018 224 69.74 −144.78 6.4 Zone B 215;219;223
14 2011 29 70.99 −6.65 6.2 Zone B 19–20;28
15 2012 243 71.44 −9.84 6.7 Zone B 233;235–239
16 2018 313 71.51 −10.81 6.8 Zone B 309–312 308–312
17 2009 232 72.22 0.84 6 Zone B 231
18 2012 145 73.01 5.59 6.3 Zone B 137–144 135;137;139;141;143–144
19 2017 8 74.44 −92.06 6.1 Zone B 366;1–7
20 2009 188 75.33 −72.49 6 Zone B 179–181
21 2008 52 77.02 19.28 6.1 Zone B 42–47;49–50
22 2009 65 80.33 −2.32 6.5 Zone B 55;58
23 2005 65 84.93 98.69 6.3 Zone B 55 55–61;64
24 2004 214 −63.65 −166.92 6 Zone C 204–213 204–211;213
25 2016 31 −63.14 169.7 6 Zone C 21–25 21–24;28
26 2014 107 −62.65 155.43 6.2 Zone C 97;101–103
27 2013 15 −62.6 −161.94 6.1 Zone C 13
28 2020 336 −61.97 154.9 6.1 Zone C 326–328;330–334
29 2007 102 −61.72 161.2 6 Zone C 92–94 92–94;100
30 2017 281 −61.56 154.32 6.2 Zone C 271–274 271–276;278–279
31 2019 204 −61.31 154.26 6.1 Zone C 194;196;198;202–203
32 2006 232 −61.27 −34.52 7 Zone C 222 224;229–231
33 2011 196 −61.12 −22.85 6 Zone C 186–189 186–188;190–193;195
34 2006 2 −61.12 −21.39 7.4 Zone C 357–361 360–363;365
35 2013 196 −61.05 −23.51 7.3 Zone C 187–195 186–195
36 2008 41 −61.05 −25.01 6.5 Zone C 33–36 31–35;38
37 2007 118 −61.04 −20.12 6.1 Zone C 113–117
38 2014 306 −61.03 153.88 6 Zone C 296 296–301
39 2009 59 −61.03 −24.39 6.3 Zone C 54–55;58
40 2014 70 −61.01 −19.92 6.4 Zone C 60–62 60;69
41 2020 208 −60.97 −25.01 6.3 Zone C 198–199 198;206–207
42 2003 216 −60.8 −43.21 7.6 Zone C 207–215 206–215
43 2009 106 −60.71 −26.55 6.7 Zone C 99–101
44 2012 283 −60.65 153.39 6.6 Zone C 274–279;282 274–275;280;282
45 2012 15 −60.62 −56.47 6.6 Zone C 5;10
46 2019 239 −60.54 −25.82 6.6 Zone C 238
47 2013 321 −60.49 −45.32 7.8 Zone C 311–320 311;313–315;319–320
48 2018 58 −60.23 150.18 6 Zone C 48–50;52–55;57
49 2000 64 −60.2 150.21 6.3 Zone C 61–63 54–59;61–62
50 2009 300 −60.05 −65.54 6 Zone C 296–299 295–298
51 2001 103 −60.04 −24.37 6.2 Zone C 93–102 94–102

According to Figure 2 (right), one can find that the values of P o decrease in Zone A and C, and increase in Zone B slightly, with depths increasing. The maximum difference of P o between different depths in Zone A, B and C was 0.0075, 0.0372 and 0.0291, respectively. Therefore, it reveals that for shallow earthquakes (≤60 km), the influence of depths on P o could be ignored. For the case of AL>20% or AL>40%, the results of Zone A are the largest, followed by those of Zone B. For the case of AL>60%, the results of Zone C are the largest, while the rest two zones have no significant difference. Therefore, selecting different AL has an obvious impact on P o in different zones with different magnitudes and depths. Considering that all the earthquakes in this study are shallow earthquakes, the impact of different depths will not be further investigated in the following sections.

To study the effects of different lengths of the sliding windows on P o , Figure 3 shows the differences between the results of different sliding windows, indicating the effect of sliding window length on latitudinal zones, that is, the most affected zones are Zone C, B and A in decreasing order. For the earthquakes in the southern mid- and high-latitude region (Zone C), compared with the results of 15-day sliding window: 1) the results of 10-day or 27-day sliding window are systematically increased, and the differences are larger using 27-day sliding window; 2) For the cases of AL>20% or AL>40%, most of the relative changes are less than 20%, while the maximum difference and relative change can reach 0.16% and 80%, respectively with AL>60%. Moreover, the differences increase with the increase of the magnitude selecting AL>60%.

The differences and relative changes of seismo-ionospheric anomalies occurrence probabilities between the results of different sliding windows. (A, C, E) show the difference between the results of 10-day and 15-day sliding windows in Zone A, Zone B and Zone C, respectively. (B, D, F) show the difference between the results of 27-day and 15-day sliding windows in Zone A, Zone B and Zone C, respectively.

For the earthquakes in the northern mid- and high-latitude region (Zone B), compared with the results of 15-day sliding window: 1) most differences are random, and are also greater using 27-day sliding window; 2) For the cases of AL>20% or AL>40%, most relative changes are between −6% and 10%. For the cases of AL>60%, the 27-day results increase systematically, and the maximum difference and relative change can reach 0.092% and 35.71%, respectively.

For the earthquakes in the low-latitude and equatorial region (Zone A), compared with the results of 15-day sliding window: 1) the results vary systematically using the 27-day sliding window; 2) For the cases of AL> 20% or AL>40%, most relative changes are between −8% and 8%. For the cases of AL>60%, the 27-day results increase systematically and the differences decrease with the magnitude increasing, with the maximum difference and relative change reaching 0.039% and 19.59%, respectively.

When AL > 20% or AL > 40% is selected to extract seismo-ionospheric anomalies, most of the relative changes of different sliding window are under 20%. In this study, GIM produced by the Jet Propulsion Laboratory (JPL) is used, and the error of VTEC estimation during the creating process of global maps by JPL is about 10%–17% (Zakharenkova et al., 2008). Therefore, for the case of AL > 20% or AL > 40%, different lengths of sliding windows have a statistically insignificant influence on P o . When AL>60% is selected, different lengths of sliding windows have non-negligible impacts on P o in the mid- and high-latitude region, particularly in the southern Hemisphere where there are obvious systematic deviations between the results of different sliding windows. The reasons for the observed systematic deviations may include: 1) the real effects of different lengths of sliding windows; 2) the number of GNSS stations in the southern Hemisphere is smaller than those in the northern Hemisphere, resulting in larger error of VTEC in the southern Hemisphere; 3) the number of earthquakes in Zone C is the smallest. 4) external magnetic fields contamination is still probably not excluded in this study for some earthquakes in Zone B and Zone C (Table 1). However, it can be noted that in Zone A, there are also obvious systematic deviations between the results of 15-day and 27-day sliding windows, indicating the influence of different sliding window lengths cannot be ignored.

Figure 4 shows the Index Number in percentage of the seismo-ionospheric anomalies within 1–10 days before the earthquakes. The Index Number is calculated as the ratio of the cumulative number of seismo-ionospheric anomalies in a single day and the total number of seismo-ionospheric anomalies (Shah et al., 2020). The Index Number enhances within 5 days before the earthquakes in all the three zones, especially in Zone A and B. Figures 5, 6 show the differences and relative changes between the results of different sliding windows. It can be seen that the differences between the results of different sliding windows are random, and are larger between the results of 15-day and 27-day sliding windows. The differences are smallest in Zone A, followed by Zone B. The differences raise with AL increasing. For the case of AL>20% or AL>40%, compared with the results of 15-day sliding window, the most of relative changes are between −20% and 25%. For the case of AL > 60%, the differences are more significant, as the relative changes can exceed 50% in all the three zones. So when AL>60% is selected, the impacts of different sliding windows on the Index Number can not be neglected.

The Index Numbers of the seismo-ionospheric anomalies within 1-10 days before the global earthquakes during 2000-2020 (15-day sliding window). (A, C, E) are the results of 6.0≤Mw<6.5 in Zone A, Zone B and Zone C, respectively. (B, D, F) are the results of Mw≥6.5 in Zone A, Zone B and Zone C, respectively.

The differences and relative changes of the Index Numbers of the seismo-ionospheric anomalies between the results of 10-day and 15-day sliding windows. (A, C, E) represent the results of 6.0≤Mw<6.5 in Zone A, Zone B and Zone C, respectively. (B, D, F) represent the results of Mw≥6.5 in Zone A, Zone B and Zone C, respectively.

The differences and relative changes of the Index Numbers of the seismo-ionospheric anomalies between the results of 27-day and 15-day sliding windows. (A, C, E) represent the results of 6.0≤Mw<6.5 in Zone A, Zone B and Zone C, respectively. (B, D, F) represent the results of Mw≥6.5 in Zone A, Zone B and Zone C, respectively.

 3.2 Seismo-ionospheric anomalies occurrence rates

According to Eq. 2, seismo-ionospheric anomalies occurrence rates P E are obtained for all the three zones. Figure 7 shows P E with AL>20%, AL>40% and AL>60% within 1–10 days before different magnitude earthquakes. We found that the larger magnitude of earthquakes and the closer prior to the earthquake occurrence, the larger values of P E . For instance, the value of P E increases from 26.82% for Mw≥6.0 earthquakes to 31.29% for Mw≥7.0 earthquakes 1 day before the earthquakes. It can be seen that the correlation is not obvious between the values of P E in Zone C and the magnitude with AL>60%. It needs to be taken into account with the time (for example, 1–5 days prior to the earthquakes) and the latitude of the epicenter (for example, in the low-latitude and equatorial region) that the values of P E increase with the magnitude increasing. In addition, as mentioned in the previous section, the number of earthquakes is smaller and the errors of VTEC are larger in Zone C, so the results in Zone C could be biased to some extent. Besides, the effects of external magnetic fields contamination for some earthquakes in Zone B and Zone C may not be excluded, which would also bias the results.

Seismo-ionospheric anomalies occurrence rates of AL>20%, 40%, and 60% within 1-10 days before earthquakes of Zone A (A, D, G), Zone B (B, E, H), and Zone C (C, F, I), respectively (15-day sliding window).

Figures 8, 9 represent the P E differences of the results using different sliding windows. For Zone A, different sliding window lengths have little influence. For Zone B, compared with the results of the 15-day sliding window, the values of P E significantly increase for Mw≥6.7 earthquakes using 27-day sliding window. For Zone C, the differences of the results of different sliding windows are larger than those in the other zones without systematic pattern. Therefore, it indicates that the impacts of sliding window lengths are related to the latitude of the epicenter. Different sliding windows should be selected to investigate the earthquakes in different zones. For Zone A, 10-day, 15-day and 27-day sliding windows may be selected arbitrarily. For Zone B, 27-day sliding window may be a better choice. But the selection of the optimal sliding window for Zone C still needs further research in the future.

The differences of seismo-ionospheric anomalies occurrence rates between the results of 10-day and 15-day sliding windows. (A, D, G) represent the results in Zone A. (B, E, H) represent the results in Zone B. (C, F, I) represent the results in Zone C.

The differences of seismo-ionospheric anomalies occurrence rates between the results of 27-day and 15-day sliding windows. (A, D, G) represent the results in Zone A. (B, E, H) represent the results in Zone B. (C, F, I) represent the results in Zone C.

To compare the difference between the seimo-ionospheric anomalies a few days before the earthquakes and the day-to-day ionospheric variation, the ionospheric anomalies occurrence rates during background days P N , as shown in Eq. 3, derived from Le et al. (2011), are also calculated. For each earthquake, 61–300 days before this earthquake are selected as the background days, and disturbed days Δ W by geomagnetic storms and by the Mw≥6.0 earthquakes at the adjacent places are also excluded. P N = ∑ n = 1 K N A L n K × 240 − Δ W

In this equation, K is the number of earthquakes in each group (for example, Mw ≥ 6.5 in Zone A with AL>20%); N A L n is the number of seismo-ionospheric days with different AL (for example, AL>20%) during the 61–300 days before the n th earthquake; Δ W is the number of disturbed days.

Figure 10A shows the values of P N for different magnitude earthquakes. The results are smaller than those of P E , and seem to decrease very slightly with the magnitude increasing. The effects of different latitudinal locations and AL on the values of P N is insignificant. Figures 10B, C) presents the differences of P N derived from different sliding window lengths are obviously systematic, and are related to AL. For example, when AL > 20% is assumed, all the results of 27-day sliding window are reduced compared with those of 15-day sliding window. Therefore, different sliding window lengths have a systematic effect on the values of P N .

(A) Ionospheric anomalies occurrence rates of A >20%, 40%, and 6% during background days; (B) The differences of ionospheric anomalies occurrence rates between the results of 10-day and 15-day sliding windows; (C) The differences of ionospheric anomalies occurrence rates between the results of 27-day and 15-day sliding windows.

Figure 11 shows the ratio of P E and P N with different magnitude using 15-day sliding window. The results show that most of the values of P E / P N are larger than 1, indicating large differences of the occurrence rates between the days prior to earthquakes and the background days.

The ratio of the seimo-ionospheric anomalies occurrence rates within 1-10 days before the earthquakes to the ionospheric anomalies occurrence rates during background days (15-day sliding window). (A, D, G) represent the results in Zone A. (B, E, H) represent the results in Zone B. (C, F, I) represent the results in Zone C.

 4 Conclusion

Both the latitude of the epicenter and the magnitude can affect the characteristics of seismo-ionospheric anomalies before shallow earthquakes. In terms of seismo-ionospheric anomalies occurrence probabilities P o , ionospheric enhancement in the low-latitude and equatorial region is more significant (Zone A), and their P o are larger with the magnitude increasing, especially for Mw≥6.6 earthquakes; in the mid- and high-latitude region (Zone B and C), a slight negative correlation is presented between the values of P o and the magnitude. In terms of seismo-ionospheric anomalies occurrence rates P E , the values of P E increase with the magnitude increasing in all the three zones, but the correlation between the values of P E and the magnitude is faint in Zone C. Because of the small number of earthquakes and the low accuracy of VTEC in Zone C, the reliability of the results still needs to be further studied and confirmed.

Both the number of days before earthquakes and the AL can affect the characteristics of seismo-ionospheric anomalies before shallow earthquakes. The number of seismo-ionospheric anomalies within 1–5 days before the earthquakes increases significantly in all the three zones, and the positive correlation between the values of P E and the magnitude is more strong within 1–5 days before the earthquakes. For the case of AL>20% and AL>40%, the values of P o in the low-latitude and equatorial region are higher than those in the mid- and high-latitude region; For the case of AL>60%, the values of P o in the southern mid- and high-latitude region are higher than those in the other zones before 6.0≤Mw < 6.7 earthquakes. Therefore, when seismo-ionospheric anomalies in different regions are analyzed, the choice of AL may affect the outcome significantly.

For the mid- and high-latitude region, the effects of different lengths of sliding windows on seismo-ionospheric anomalies can not be ignored. For the case of AL>60%, different sliding windows have a significant impact on the values of P o and P E , and there are systematic deviations between the values of P o using different sliding windows in the southern mid- and high-latitude region. However, the selection of the optimal sliding window needs to be further studied, especially for the Zone C.

Moreover, there are large differences between the seismo-ionospheric anomalies occurrence rates P E 1–10 days prior to the earthquakes and ionospheric anomalies occurrence rates P N during the background days, indicating that the seismo-ionospheric anomalies within just a few days before the earthquakes are probably related with the forthcoming earthquakes. The 1522 earthquakes in this study are divided into multiple groups, resulting in the small number of earthquakes in some groups, and the accuracy of VTEC in different regions is uneven, which may bias the results in this study. Although the earthquakes occurring at the similar location but with the short interval (<10 days) from the previous ones are excluded, some earthquakes occurring close to their boundaries in two confining zones may also have some confounded effects on the results. Besides, whether a pixel of 2 . 5 ∘ × 5 ∘ is appropriate for the earthquakes larger than Mw6.0 needs to be further studied.

 Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

 Author contributions

MY and ZX contributed to conception of the study. Data were collected by HL, DH, and YR. MY, ZX, and YY performed the experimental analysis. MY wrote the first draft of the manuscript. ZX revised the manuscript. All authors read and approved the final manuscript.

 Funding

This work was supported by Beijing Natural Science Foundation Grant Numbers 8224094.

We are grateful to the World Data Center for Geomagnetism, Kyoto providing Dst index, the Global CMT Project providing parameters of earthquakes, and IGS providing GIM.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Dautermann T. Calais E. Haase J. Garrison J. (2007). Investigation of ionospheric electron content variations before earthquakes in southern California, 2003-2004. J. Geophys. Res. 112, 1074–1086. 10.1029/2006JB004447 Davidenko D. V. Pulinets S. A. (2019). Deterministic variability of the ionosphere on the eve of strong (M≥6) earthquakes in the regions of Greece and Italy according to long-term measurements data. Geomagn. Aeron. 59 (4), 493–508. 10.1134/S001679321904008X De Santis A. Marchetti D. Pavón-Carrasco F. J. Cianchini G. Perrone L. Abbattista C. (2019). Precursory worldwide signatures of earthquake occurrences on Swarm satellite data. Sci. Rep. 9, 20287. 10.1038/s41598-019-56599-1 Dobrovolsky I. P. Zubkov S. I. Miachkin V. I. (1979). Estimation of the size of earthquake preparation zones. Pure Appl. Geophys. 117 (5), 1025–1044. 10.1007/BF00876083 Freund F. (2011). Pre-earthquake signals: Underlying physical processes. J. Asian Earth Sci. 41 (4-5), 383–400. 10.1016/j.jseaes.2010.03.009 Fujiwara H. Kamogawa M. Ikeda M. Liu J. Y. Sakata H. Chen Y. I. (2004). Atmospheric anomalies observed during earthquake occurrences. Geophys. Res. Lett. 31, L17110. 10.1029/2004GL019865 Guo J. Li W. Yu H. Liu Z. Zhao C. Kong Q. (2015). Impending ionospheric anomaly preceding the Iquique Mw8.2 earthquake in Chile on 2014 April 1. Geophys. J. Int. 203 (3), 1461–1470. 10.1093/gji/ggv376 Hernández-Pajares M. Juan J. Sanz J. Orus R. Garcia-Rigo A. Feltens J. (2009). The IGS VTEC maps: A reliable source of ionospheric information since 1998. J. Geod. 83 (3-4), 263–275. 10.1007/s00190-008-0266-1 Ho Y. Y. Jhuang H. K. Su Y. C. Liu J. Y. (2013). Seismo-ionospheric anomalies in total electron content of the GIM and electron density of DEMETER before the 27 February 2010 M8.8 Chile earthquake. Adv. Space Res. 51 (12), 2309–2315. 10.1016/j.asr.2013.02.006 Jiang W. Ma Y. Zhou X. Li Z. An X. Wang K. (2017). Analysis of ionospheric vertical total electron content before the 1 April 2014 Mw 8.2 Chile earthquake. J. Seismol. 21 (6), 1599–1612. 10.1007/s10950-017-9684-y Ke F. Wang Y. Wang X. Qian H. Shi C. (2016). Statistical analysis of seismo-ionospheric anomalies related to Ms>5.0 earthquakes in China by GPS TEC. J. Seismol. 20 (1), 137–149. 10.1007/s10950-015-9516-x Klimenko M. V. Klimenko V. V. Zakharenkova I. E. Pulinets S. A. (2012). Variations of equatorial electrojet as possible seismo-ionospheric precursor at the occurrence of TEC anomalies before strong earthquake. Adv. Space Res. 49 (3), 509–517. 10.1016/j.asr.2011.10.017 Kon S. Nishihashi M. Hattori K. (2011). Ionospheric anomalies possibly associated with M ≥ 6.0 earthquakes in the Japan area during 1998–2010: Case studies and statistical study. J. Asian Earth Sci. 41 (4), 410–420. 10.1016/j.jseaes.2010.10.005 Le H. Liu J. Y. Liu L. (2011). A statistical analysis of ionospheric anomalies before 736 M6.0+ earthquakes during 2002-2010. J. Geophys. Res. 116 (A2), A02303. 10.1029/2010JA015781 Liperovsky V. A. Pokhotelov O. A. Meister C. V. Liperovskaya E. V. (2008). Physical models of coupling in the lithosphere-atmosphere-ionosphere system before earthquakes. Geomagn. Aeron. 48, 795–806. 10.1134/S0016793208060133 Liu J. Y. Chen C. H. Chen Y. I. Yang W. H. Oyama K. I. Kuo K. W. (2010). A statistical study of ionospheric earthquake precursors monitored by using equatorial ionization anomaly of GPS TEC in Taiwan during 2001–2007. J. Asian Earth Sci. 39 (1), 76–80. 10.1016/j.jseaes.2010.02.012 Liu J. Y. Chen Y. I. Chuo Y. J. Tsai H. F. (2001). Variations of ionospheric total electron content during the Chi-Chi earthquake. Geophys. Res. Lett. 28 (7), 1383–1386. 10.1029/2000GL012511 Liu J. Y. Chuo Y. J. Shan S. J. Tsai Y. B. Chen Y. I. Pulinets S. A. (2004). Pre-earthquake ionospheric anomalies registered by continuous GPS TEC measurements. Ann. Geophys. 22 (5), 1585–1593. 10.5194/angeo-22-1585-2004 Liu W. Xu L. (2017). Statistical analysis of ionospheric TEC anomalies before global Mw≥7.0 earthquakes using data of CODE GIM. J. Seismol. 21 (4), 759–775. 10.1007/s10950-016-9634-0 Ma X. Lin Z. Chen H. (2014). Analysis on ionospheric perturbation of TEC and NmF2 based on GPS and COSMIC data before and after the Wenchuan earthquake. Chin. J. Geophys 57 (8), 2415–2422. 10.6038/cjg20140803 Masci F. (2012). The study of ionospheric anomalies in Japan area during 1998–2010 by Kon et al.: An inaccurate claim of earthquake-related signatures? J. Asian Earth Sci. 57, 1–5. 10.1016/j.jseaes.2012.06.009 Masci F. Thomas J. (2014). Comment on “Temporal and spatial precursors in ionospheric total electron content of the 16 October 1999 M w  7.1 Hector Mine earthquake” by Su et al. (2013). J. Geophys. Res. 119, 6994–6997. 10.1002/2014JA019896 Michael A. J. (2011). Random variability explains apparent global clustering of large earthquakes. Geophys. Res. Lett. 38, L21301. 10.1029/2011GL049443 Namgaladze A. A. Klimenko M. V. Klimenko V. V. Zakharenkova I. E. (2009). Physical mechanism and mathematical modeling of earthquake ionospheric precursors registered in total electron content. Geomagn. Aeron. 49 (2), 252–262. 10.1134/S0016793209020169 Nishihashi M. Hattori K. Jhuang H. K. Liu J. Y. (2009). Possible spatial extent of ionospheric GPS-TEC and NmF2 anomalies related to the 1999 Chi-Chi and Chia-Yi earthquakes in Taiwan. Terr. Atmos. Ocean. Sci. 20, 779–789. (T). 10.3319/tao.2009.01.22.01(t) Ovalle E. Villalobos C. Bravo M. Foppiano A. (2013). Maximum electron concentration and total electron content of the ionosphere over Concepción, Chile, prior to the 27 February 2010 earthquake. Adv. Space Res. 52, 1274–1288. 10.1016/j.asr.2013.07.005 Parrot M. Li M. (2018). “Statistical analysis of the ionospheric density recorded by DEMETER during seismic activity,” In pre-earthquake processes: a multidisciplinary approach to earthquake prediction studies (geophysical monograph 234). Editors Ouzounov D. Pulinets S. Hattori K. Taylor P. (Washington, DC: American Geophysical Union John Wiley and Sons, Inc.). Pulinets S. A. Boyarchuk K. (2004). Ionospheric precursors of earthquakes. Berlin: Springer Press. Pulinets S. A. Davidenko D. V. (2018). The nocturnal positive ionospheric anomaly of electron density as a short-term earthquake precursor and the possible physical mechanism of its formation. Geomagn. Aeron. 58, 559–570. 10.1134/s0016793218040126 Pulinets S. A. Krankowski A. Hernandez-Pajares M. Marra S. Budnikov P. Zakharenkova I. (2021). Ionosphere sounding for pre-seismic anomalies identification (INSPIRE): Results of the Project and perspectives for the short-term earthquake forecast. Front. Earth Sci. 9, 610193. 10.3389/feart.2021.610193 Pulinets S. (2012). Low-latitude atmosphere-ionosphere effects initiated by strong earthquakes preparation process. Int. J. Geophys 2012, 1–14. 10.1155/2012/131842 Pulinets S. Ouzounov D. (2011). Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) model-An unified concept for earthquake precursors validation. J. Asian Earth Sci. 41 (4), 371–382. 10.1016/j.jseaes.2010.03.005 Rishbeth H. (2006). Ionoquakes: Earthquake precursors in the ionosphere? Eos. Trans. Am. Geophys. Union 87, 316. 10.1029/2006EO320008 Shah M. Ahmed A. Ehsan M. Khan M. Tariq M. A. Calabia A. (2020). Total electron content anomalies associated with earthquakes occurred during 1998-2019. Acta Astronaut. 175, 268–276. 10.1016/j.actaastro.2020.06.005 Shah M. Jin S. (2015). Statistical characteristics of seismo-ionospheric GPS TEC disturbances prior to global Mw≥5.0 earthquakes (1998–2014). J. Geodyn. 92, 42–49. 10.1016/j.jog.2015.10.002 Su Y. C. Liu J. Y. Chen S. P. Tsai H. F. Chen M. Q. (2013). Temporal and spatial precursors in ionospheric total electron content of the 16 october 1999M w7.1 hector mine earthquake: Ionoquack OFM7.1 hector mine earthquake. J. Geophys. Res. 118 (10), 6511–6517. 10.1002/jgra.50586 Sun Y. Liu J. Lin C. Tsai H. F. Chang L. C. Chen C. Y. (2016). Ionospheric f2 region perturbed by the 25 april 2015 Nepal earthquake: Ionosphere perturbed by Nepal earthquake. J. Geophys. Res. 121, 5778–5784. 10.1002/2015JA022280 Tang J. Yao Y. Zhang L. (2015). Temporal and spatial ionospheric variations of 20 April 2013 earthquake in Yaan, China. IEEE Geosci. Remote Sens. Lett. 12, 2242–2246. 10.1109/LGRS.2015.2463081 Xu T. Chen Z. Li C. Wu J. Hu Y. Wu Z. (2011). GPS total electron content and surface latent heat flux variations before the 11 March 2011 M9.0 Sendai earthquake. Adv. Space Res. 48 (8), 1311–1317. 10.1016/j.asr.2011.06.024 Yao Y. Chen P. Wu H. Zhang S. Peng W. (2012). Analysis of ionospheric anomalies before the 2011 Mw 9.0 Japan earthquake. Chin. Sci. Bull. 57, 500–510. 10.1007/s11434-011-4851-y Zakharenkova I. E. Shagimuratov I. I. Tepenitzina N. Y. Krankowski A. (2008). Anomalous modification of the ionospheric total electron content prior to the 26 September 2005 Peru earthquake. J. Atmos. Sol-terr. Phy. 70 (15), 1919–1928. 10.1016/j.jastp.2008.06.003 Zhang X. Shen X. Zhao S. Yao L. Ouyang X. Qian J. (2014). The characteristics of quasistatic electric field perturbations observed by DEMETER satellite before large earthquakes. J. Asian Earth Sci. 79 (2), 42–52. 10.1016/j.jseaes.2013.08.026 Zhou Y. Wu Y. Qiao X. Zhang X. (2009). Ionospheric anomalies detected by ground-based GPS before the Mw7.9 Wenchuan earthquake of May 12, 2008, China. J. Atmos. Sol-terr. Phy. 71, 959–966. 10.1016/j.jastp.2009.03.024 Zhu F. Zhou Y. Lin J. Su F. (2014). A statistical study on the temporal distribution of ionospheric TEC anomalies prior to M7.0+ earthquakes during 2003-2012. Astrophys. Space Sci. 350, 449–457. 10.1007/s10509-014-1777-2 Zolotov O. V. Knyazeva M. A. Romanovskaya Y. V. (2019). Computer analysis of total electron content in the Earth’s ionosphere in problems of searching for and detection of earthquake precursors: Current problems and challenges. Russ. J. Phys. Chem. B 13 (4), 681–684. 10.1134/S1990793119040146 Zolotov O. V. Namgaladze A. A. Zakharenkova I. E. Martynenko O. V. Shagimuratov I. I. (2012). Physical interpretation and mathematical simulation of ionospheric precursors of earthquakes at midlatitudes. Geomagn. Aeron. 52 (3), 390–397. 10.1134/S0016793212030152