*Article* **A Study of Possible Correlations between Seismo-Ionospheric Anomalies of GNSS Total Electron Content and Earthquake Energy**

**Yung-Chih Su 1,2,3 and Jinming Sha 1,2,3,\***


**Abstract:** In this study, we conduct a correlation analysis between the daily occurrence times of the increase and decrease anomalies in the global total electron content (TEC) in the ionosphere, and the daily earthquake energy release within 110–130◦E longitude over the following three latitude regions: A: 13◦S–0.5◦S (22.3◦S–10◦S geomagnetic), B: 0.5◦S–19.5◦N (10◦S–10◦N geomagnetic), and C: 19.5◦N–32.1◦N (10◦N–22.5◦N geomagnetic). The TEC data from global ionosphere maps (GIMs) during earthquake events of *M* ≥ 2.5 that occurred in 2015–2018 are used in this study. The time series of daily seismic wave energy releases within the three regions and the daily occurrence times of the TEC anomalies in each GIM grid are computed. By time-shifting the time series, the correlations are calculated and compared globally, and the temporal characteristics are also examined. The disturbance storm time (Dst) index, planetary geomagnetic index Kp, and daily observed 10.7 cm solar flux (F10.7) are used to remove data associated with space weather variations. Although the seismo-ionospheric precursor is not confirmed by the statistical investigations, the greater occurrence times of TEC decrease anomalies are observed in the southeast in Region A, and the conjugate point 13 days prior to a *M*6.9 earthquake in Region A, which occurred on 5 August 2018, in accordance with the statistical results. Therefore, it is required to apply more parameters to understand the causes of the ionospheric TEC variations and investigate whether ionospheric variations are caused by earthquakes.

**Keywords:** earthquake energy; total electron content; global ionosphere maps; seismo-ionospheric anomaly

#### **1. Introduction**

Earthquakes are one of the most devastating types of natural disasters, and large earthquakes can cause great damage to life and property. Therefore, predicting the occurrence of earthquakes in advance has been of interest to researchers and the general public. Scientific attempts related to earthquake prediction include monitoring variations in the geomagnetic field, surface displacements, groundwater levels, etc. [1–5]. In addition, many studies have shown that the electron density in the Earth's ionosphere might vary anomalously within approximately 7 days prior to earthquakes [6–11]. Liu et al. [6,7] conducted regional ionospheric observations using the Chung-Li ionosonde (25.0◦N, 121.2◦E) and the Global Positioning System (GPS) network in the Taiwan area to observe the maximum plasma frequency foF2 in the ionospheric F2 region and the total electron content (TEC), respectively, and detected ionospheric variations before the 1999 *M*7.7 Chi-Chi earthquake, which occurred in central Taiwan. The results show that the ionospheric electron density over the epicenter anomalously decreased 1, 3, and 4 days before the Chi-Chi earthquake.

**Citation:** Su, Y.-C.; Sha, J. A Study of Possible Correlations between Seismo-Ionospheric Anomalies of GNSS Total Electron Content and Earthquake Energy. *Remote Sens.* **2022**, *14*, 1155. https://doi.org/ 10.3390/rs14051155

Academic Editors: Shuanggen Jin and Gino Dardanelli

Received: 3 December 2021 Accepted: 22 February 2022 Published: 26 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

To understand the ionospheric precursor signatures, it is important to examine global ionospheric variations prior to earthquakes. Global ionosphere maps (GIMs) provided by the Center for Orbit Determination in Europe (CODE) are based on time-continuous observations with global coverage and can be used to investigate the global distributions of ionospheric anomalies. Liu et al. [8–10] employed GIMs to investigate the seismoionospheric precursors (SIPs) for the 2008 *M* 8.0 Wenchuan, China, the 2004 *M* 9.1 Sumatra, and the 2010 *M* 7.0 Haiti earthquakes, respectively. They found that ionospheric TEC anomalies appeared around the epicenters 1–6 days prior to the three earthquakes. Furthermore, satellite observations have detected ionospheric variations before earthquakes. For example, Liu et al. [11] utilized data obtained by the French satellite Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions (DEMETER) to analyze the differences between the observed ionospheric parameters before and after the 2008 Wenchuan earthquake and found that the ionospheric nighttime electron density and ion density decreased significantly over the epicenter 1–6 days before the earthquake. On 2 February 2018, the China Seismo-Electromagnetic Satellite (CSES) was launched. It measures the electric field, magnetic field, electron density, ion density, and other important parameters. Yan et al. [12] showed that the ionospheric parameters observed by CSES around the epicenters of four earthquakes with a *M* > 7.0 during August 2018 were perturbed prior to the earthquakes. Statistical results also support pre-earthquake anomalous variations in the ionospheric electron density [6,8,9,13–15]. Several mechanisms have been proposed to explain the occurrence of SIPs, including atmospheric conductivity changes in the global electric circuit (GEC) related to ionization caused by an increase in the radon emanations before earthquakes [16], amplification of internal gravity waves (IGWs) produced due to surface motion through interaction with planetary waves and their subsequent propagation to the ionosphere where they modify the dynamo electric fields [17], and large-scale electric field-induced ionospheric perturbations [18,19].

There are also studies showing the opposite findings. For example, Masci et al. [20] analyzed ionospheric TEC during the 6 April 2009 *M*6.1 L'Aquila (Italy) earthquake and showed that the hump-like shape in the TEC difference (DTEC) time series calculated from two GPS receivers in central Italy, which was interpreted as a possible earthquake effect by Nenovski et al. [21], appears to be a diurnal variation and cannot be considered an earthquake-related effect. Dautermann et al. [22] and Thomas et al. [23] investigated long periods of data; the former used 2 years (2003–2004) of TEC data from the Southern California Integrated GPS Network (SCIGN) for earthquakes in southern California, and the latter used GIMs for the 1279 *M* ≥ 6.0 earthquakes globally during the years 2000–2014. Both studies showed no evidence to support ionospheric perturbations prior to earthquakes. The supporting and opposite results might be because the variability in SIPs might be large, making them sometimes difficult to detect.

Many studies have shown that the ionospheric anomalies in GIM TEC associated with earthquakes might have frequent or prolonged occurrence times (e.g., Liu et al. [10]). In this study, we would like to investigate the possibility of establishing an earthquake prediction model by using GIM TEC, considering the frequent or prolonged feature of SIPs. This study focuses on earthquakes that occurred during 2015–2018 within three selected low geomagnetic latitude regions in the 110–130◦E longitude sector. The day is used as a unit, and the daily seismic energy released within the three latitude regions is computed. The correlation coefficients between the seismic energy and the daily occurrence times of GIM TEC anomalies within each GIM grid for ±15 days from the energy release are calculated. Thus, we could compare the correlations globally and examine the temporal characteristics. The data that are associated with geomagnetic disturbances by using the disturbance storm time (Dst) index are removed. In addition, the results of using declustered seismic energy and removing data by using Dst and considering the planetary geomagnetic index Kp and daily observed 10.7 cm solar flux (F10.7) are investigated.

#### **2. Materials and Methods**

#### *2.1. Earthquake Data*

In this study, the earthquake data from U.S. Geological Survey (USGS, https:// earthquake.usgs.gov/earthquakes/search/, last accessed on 26 January 2022) and GIM TEC from CODE (ftp://ftp.aiub.unibe.ch/CODE/, last accessed on 29 August 2021) are used. All earthquakes that occurred within a longitude range of 110–130◦E and during 1 January 2015, to 31 December 2018, with *M* ≥ 2.5 from the USGS are included in this study. We presumed that the uncertainties, which might be associated with the pre-earthquake effects near the Earth's surface and their coupling to the ionosphere, such as magnetic dip, and the geographically related conditions (for example, ocean, urban area, etc.), might produce different ionospheric variations, and the selected longitude sector is subdivided into three latitude regions, which are A: 13◦S–0.5◦S (22.3◦S–10◦S geomagnetic), B: 0.5◦S–19.5◦N (10◦S–10◦N geomagnetic), and C: 19.5◦N–32.1◦N (10◦N–22.5◦N geomagnetic). The three regions are separated mainly according to the magnetic dip angle. The distribution of the earthquakes used is shown in Figure 1.

**Figure 1.** The distribution of earthquakes of magnitude ≥ 2.5 occurring within 110–130◦E longitude from 1 January 2015 to 31 December 2018 (white dots), and the global distribution of GIM TEC observed at 1800 UT on 21 March 2018 in TEC units (1 TECU = 10 <sup>16</sup> electron/m2). Two vertical white lines denote longitudes of 110 and 130◦E. The three study Regions A (13◦S–0.5◦S), B (0.5◦S–19.5◦N), and C (19.5◦N–32.1◦N) are divided by four red lines. The equatorial ionization anomaly (EIA) can be seen between approximately −150 and 30◦E longitudes.

It is shown that SIPs might be related to earthquake energy [13,14]. One of the variables used in computing the correlation coefficient is the daily total energy of seismic waves radiated by earthquakes in a given latitude region. The seismic wave energy released by an earthquake of magnitude *M* was estimated using the following relationship developed by Gutenberg and Richter [24]:

$$
\log E = 11.8 + 1.5 \, M \tag{1}
$$

where *E* is seismic wave energy in erg (10−<sup>7</sup> J). The total seismic wave energy released within a region by all the earthquakes was then summed for each day during 2015–2018. The time series of daily seismic energy and the temporal distributions of earthquakes with different magnitude ranges within 110–130◦E and in Regions A, B, and C are provided in Figure S1 in Supplementary Material S1. There were 123,505 seismic events with *<sup>M</sup>* ≥ 2.5 in total worldwide, releasing a total seismic wave energy of 1.28 × <sup>10</sup> <sup>25</sup> erg during 2015–2018. Table 1 shows that there were 8080 earthquakes with *M* ≥ 2.5 in 110–130◦E during 2015–2018, with four *M* ≥ 7 earthquakes among them, among which one occurred in Region A and three in Region B. The total seismic wave energy of the 8080 earthquakes was 3.43% of the net seismic energy released worldwide during 2015–2018. Among the three Regions A–C, Region B spans the widest latitude range and has the largest earthquake

number and total seismic energy, while Region C has the smallest earthquake number and energy.


**Table 1.** Earthquake number and energy of *M* ≥ 2.5 within 110–130◦E in 2015–2018.

Note: The earthquake catalog: USGS (https://earthquake.usgs.gov/earthquakes/search/, last accessed on 26 January 2022).

#### *2.2. GIM TEC*

TEC is the integral of electron density along the signal path from the global navigation satellite system (GNSS) satellite to the receiver. Since the signal path to the ground-based receiver is tilted, the following measured TEC is called slant TEC (STEC):

$$STEC = \int\_{\text{Sat}}^{\text{Rx}} N\_{\text{c}} ds\_{\text{s}} \tag{2}$$

where *Ne* is the electron density and Sat and Rx are the locations of the GNSS satellite and ground-based receiver, respectively. As the satellite's zenith angle changes, the length of the path of the signal through the ionosphere changes, and hence, STECs have to be converted into vertical TECs (VTECs) to eliminate the zenith angle effect. The GIM TECs provided by CODE are generated by using data from approximately 300 GNSS sites worldwide using spherical harmonics expansion. The GIM covers a range of ±87.5◦N latitude and ±180◦E longitude, with spatial resolutions of 2.5◦ and 5◦, respectively. During the study period, the CODE GIM provides 24-h global TEC distributions (from 00 to 23 UT) each day [25]. Please note that the accuracy of GIM is lower over the ocean areas. However, the GIM data are suitable statistical data for global analysis and to investigate the possible obvious and prolonged SIPs. The global distribution of GIM TEC at 1800 UT on 21 March 2018, is displayed in Figure 1.

#### *2.3. TEC Anomalies*

We use the method used in Liu et al. [8] and Chen et al. [13] to define anomalies in the GIM TEC variations. For an observation (*O*) at a certain time and location (GIM grid), the median TEC and the corresponding first and third quartiles, denoted by *M*, *LQ,* and *UQ*, respectively, are computed based on the TECs of the previous 15 days at the same time and location. The upper bound (*UB*) and lower bound (*LB*) are defined as the following:

$$\begin{array}{l} LB = M + k(MQ - M), \\ LB = M - k(M - LQ). \end{array} \tag{3}$$

We set *k* to be 1.5 [8–10,13]. If *O* > *UB* (*O* < *LB*), then an anomalous increase (decrease) in the TEC observation at this time and location is declared. Under the assumption of a normal distribution, the probability of an observed TEC in the interval (*LB*, *UB*) is approximately 69%. The occurrence times of the increase anomaly and decrease anomaly and the sum of the two are counted separately for each day. The daily occurrence times of TEC anomalies (the maximum is 24 times), together with the daily seismic wave energy release, are used to compute the correlation coefficients.

#### *2.4. Removing Storm Effects and Other Experimental Designs*

When magnetic disturbances or magnetic storms occur, the ionospheric electron density is perturbed and could anomalously increase and/or decrease relative to that in a quiet time period. Liu et al. [26] showed that storm-related ionospheric perturbations can last as long as 4 days after storm onset. The Dst index, obtained from the average of the horizontal component of the low-latitude magnetic field, can be used to identify storm onset. The Dst temporal variations from the World Data Center (WDC) for Geomagnetism in Kyoto (http://wdc.kugi.kyoto-u.ac.jp/dstdir/index.html, last accessed on 24 January 2022) during the study period were employed. In this study, we refer to Zhu et al. [27], and if Dst < −40 nT or the absolute value of the Dst difference from the previous 1-h value > 40 nT, the occurrence times of TEC anomalies from 1 day before to 4 days after this UT day (6 days in total) were not considered. Figure S1 in Supplementary Material S1 shows the days that were not considered in this study.

To investigate the temporal correlation, the time series of occurrence times of TEC increase anomalies, decrease anomalies, and both anomalies together at each GIM grid were shifted ±15 days forward and backward relative to the time series of daily seismic energy release. Therefore, the correlation coefficients between the daily seismic wave energy release within 110–130◦E in Regions A, B, and C and the occurrence times of TEC anomalies in each GIM grid prior to and after energy release could be computed. During the study period, the average sample size used in calculating correlation coefficients is 905 days (with a minimum of 893 and a maximum of 917 days depending on the day shifts) after removing the storm effects by using the Dst index.

To test whether the two variables are correlated (population correlation coefficient = 0), the *t*-test was applied. The *t*-value (*t*) is given by the following:

$$t = r\sqrt{\frac{n-2}{1-r^2}}\tag{4}$$

where *r* is the correlation coefficient and *n* is the sample size. The null hypothesis is that is not significantly different from zero, i.e., there is not a significant linear relationship between the occurrence times of TEC anomalies and the earthquake energy in the population [28]. Setting a significance level (α) of 0.01, the critical value *t*<sup>α</sup> = 0.01 for the two-tailed test with a sample size of 905 is approximately 2.581, and the associated *r* is 0.086. If *r* > 0.086 (*r* < −0.086), we consider the occurrence times of TEC anomalies to be positively (negatively) correlated with the earthquake energy. In addition, the *p* value (*p*), which is defined as the probability that *t* is greater than the observed *t* with a degree of freedom of *n* − 2, was also computed. The minimum *p* was found to further examine the spatial relationships between the positive correlations and earthquakes.

Since larger earthquakes are more likely to produce SIPs [13,14], for Regions A, B, and C, the sample pairs with daily seismic wave energy < 1.1220 × <sup>10</sup> <sup>20</sup> erg (equivalent to an *M*5.5 earthquake) were removed in the following analyses. In addition, a 1-unit difference in *M* corresponds to a difference of approximately 32 times in earthquake energy, and 2 units correspond to 1000 times (Equation (1)). Thus, there could be effects from the outliers. In the analyses of Regions A, B, and C, the outliers were found from seismic wave energies that were ≥ *M*5.5 by using a 1.5 interquartile range above or below the associated two quartiles, and the associated sample pairs were further removed to eliminate the possible effects from outliers. The outliers were removed, and the ranges of seismic energy were fixed before eliminating the storm effects. Table 2 lists the parameters of the samples with daily seismic wave energy release ≥ *M*5.5 with/without outliers in Regions A, B, and C. The threshold of the correlation coefficient is α = 0.01 (*r*<sup>α</sup> = 0.01), and maximum daily seismic wave energy releases for samples with energy release ≥ *M*5.5 and with/without outliers for Regions A, B, and C are computed. The parameters vary depending on whether the storm effect is significant when the days are shifted.


**Table 2.** Parameters for daily energy release equivalent to *M* ≥ 5.5 earthquakes.

Note: The *r*<sup>α</sup> = 0.01 represents the threshold value of *r* at significance level α = 0.01.

#### **3. Results**

#### *3.1. 110–130*◦*E*

In this section, we briefly describe the results of the correlation analyses. The results are presented in Supplementary Material S2. The figures which are of more importance or could be used as examples are also shown in the paper. Figures 2–4 (Figures S1–S3 in Supplementary Material S2) show the temporal and spatial distributions of *r* between the daily released seismic wave energy within 110–130◦E (all the earthquakes in Figure 1 were considered) and its associated occurrence times of TEC anomalies 15 days before and after the energy release. Panels from left to right denote the days from daily energy release (D); for example, D = −7 represents the 7 days before the energy release. The main focus of the study is over 120◦E longitude. However, we also wanted to examine whether similar variations existed in other longitude sectors for the same selected earthquakes. From top to bottom, the four longitude sectors of 120 (the study area of earthquakes), 30, −60, and −150◦E, with a width of ±20◦ and separated by 90◦, are displayed. Each panel displays an area of a 40◦ longitude range (x-axis) and a ±90◦N latitude range (y-axis). Note that the results of the −150◦E sector have a lower accuracy since most areas in this sector are ocean areas, and the results are provided for comparisons and references. For the decrease anomaly (Figure 2), a larger area of positive correlation is seen in the 120◦E sector from approximately 30–60◦N at D = −13, and many positive correlations appear in northern higher latitudes and other latitudes at approximately D = −14 to −7. The positive correlations are more apparent from D = −6 to −1 over the 120 and 30◦E sectors with changing latitudes. After the energy release, there were fewer positive correlations in the four longitudinal sectors. For the increase anomaly (Figure 3), positive correlations mainly appear at approximately D = 3 to 12, with D = 3 to 7 in the southern higher latitudes and D = 8 and 9 in all four sectors. For the occurrence times of both increase and decrease anomalies (Figure 4), positive correlations are relatively rare and sporadic in time and space. There are only a few negative correlations for all three types of anomalies in Figures 2–4.

#### *3.2. Region A*

Figure 5 shows the temporal and spatial distributions of *r* between the daily seismic wave energy released within Region A and its associated occurrence times of TEC decrease anomalies in ±15 days. Positive correlations are evident mainly in the northern mid- and high latitudes, and roughly during D = −14 to −10, −7 to −5, and −3 to 0 in Figure 5. There are only a few correlations for D > 0. In contrast, the distributions of the positive correlation for the increase anomalies are relatively irregular, and the possible precursory information is not as clear as shown in Figure S5, which is provided in Supplementary Material S2. In addition, no specific feature in the temporal and spatial distribution of the correlation for both increase and decrease anomalies counted together could be identified. The results for using all energy releases available within Region A are provided as Figures S4–S6 in Supplementary Material S2.

**Figure 2.** Temporal and spatial distributions of correlation coefficient between daily seismic wave energy release within 110–130◦E and occurrence times of TEC decrease anomalies. Panels from left to right represent −15 to 15 days from the seismic wave energy release (D). Top to bottom panels display four selected longitude sectors centered over 120, 30, −60, and −150◦E, respectively, each with a width of ±20◦. The x-axis of each panel denotes the 40◦ longitude range, and the y-axis denotes latitude. The contour denotes *r*<sup>α</sup> = 0.01 = ±0.086.

**Figure 3.** Same as Figure 2 for the results of occurrence times of TEC increase anomalies.

**Figure 4.** Same as Figure 2 for the results of occurrence times of both the TEC increase and decrease anomalies counted together.

**Figure 5.** Temporal and spatial distributions of correlation coefficient between the daily seismic wave energy release within Region A and the occurrence times of the TEC decrease anomalies. Panels from left to right represent −15 to 15 days from the seismic wave energy release (D). Top to bottom panels display four selected longitude sectors centered over 120, 30, −60, and −150◦E, respectively, each with a width of ± 20◦. The contour denotes *r*<sup>α</sup> = 0.01 = ±0.086. The black box on the top first panel indicates Region A.

Figure 6 shows the temporal and spatial distributions of *r* for samples with daily seismic wave energy ≥ *<sup>M</sup>*5.5 (1.1220 × <sup>10</sup> <sup>20</sup> erg) released within Region A. As the sample size is reduced, the required correlation coefficients for significance at the 99% level increase. On D = −13, there is a larger area of positive correlation for the occurrence times of decreasing anomalies around the northern mid-latitudes of the 120◦E sector. Positive correlations in the northern mid- and high latitudes can be seen when D = −7 to 0. They also appear at other latitudes and longitudes and after energy release. Figure 7 shows the results of the decrease anomaly when outliers of daily seismic wave energy ≥ *M*5.5 were removed. Positive correlations can be found in the mid- and high latitudes of both hemispheres, as well as in the low latitudes, during D = −2 to 5. There are fewer correlations and no obvious pattern for the TEC increase anomaly for energy release ≥ *M*5.5. After removing outliers of larger energy releases and for increase anomalies, positive correlations appear on D = −12 to −11 and −8 to −4 in the four sectors. The results for energy releases ≥ *M*5.5 with/without outliers within Region A are provided as Figures S7–S10.

**Figure 6.** Same as Figure 5 for daily seismic wave energy release ≥ *M*5.5 with outliers within Region A and for the TEC decrease anomalies. The contour denotes *r*<sup>α</sup> = 0.01 according to the sample size at each D.

#### *3.3. Region B*

The temporal and spatial distributions of *r* between the seismic wave energy released within Region B and the associated occurrence times of the TEC decrease anomalies and increase anomalies in a period of ±15 days are similar to those obtained from all the earthquakes occurring within 110–130◦E during the four years, as shown in Figures 2 and 3. The correlation features between the occurrence times of both increase and decrease anomalies counted together and earthquake energy are fewer. The results are provided in Figures S11–S13 in Supplementary Material S2.

Figure 8 shows the results for samples with an energy release ≥ *M*5.5 and for the decrease anomalies. It is revealed that apparent positive correlations of the decrease anomaly appear between D = −6 and 0. Correlations are also seen up toD=6 after energy release. After removing outliers with large energy releases, the apparent correlations from D = −6 to 0 in Figure 8 disappear. Positive correlations appear during D = −14 to −11, and the possible

precursory information becomes less clear. The distributions of positive correlations of increase anomalies with energy release ≥ *M*5.5 are more similar to those in Figure 3 (apparent patterns between approximately D = 3 and 12). After the outliers were removed, there were positive correlations on D = −10 to −3, but the precursory information is less clear. The positive correlations obtained with outliers shown in Figures 8 and S16 could be caused by the outlier with a large energy release. The results for energy releases ≥ *M*5.5 with/without outliers within Region B are shown in Figures S14–S17.

**Figure 7.** Same as Figure 5 for the daily seismic wave energy release ≥ *M*5.5 without outliers within Region A and for the TEC decrease anomalies. The contour denotes *r*<sup>α</sup> = 0.01 according to sample size at each D.

#### *3.4. Region C*

The temporal and spatial distributions of *r* for the seismic wave energy released within Region C using all the samples available during the study period show that there is no clear precursory information that could be directly identified for the TEC decrease anomalies and increase anomalies since the positive correlations appear in the four sectors and after energy release. There are also fewer precursory features associated with both the increase and decrease anomalies counted together based on the correlation coefficients (Figures S18–S20).

The results for samples with an energy release ≥ *M*5.5 with outliers and for decrease anomalies and increase anomalies are presented in Figures 9 and 10, respectively, because we also reference them in the following analyses. The sample size for an energy release ≥ *M*5.5 in Region C is relatively small (14–20, see Table 2), and the correlation coefficients could be very large. Figure 9 shows that positive correlations for the decrease anomalies appear in large areas over the globe at approximately D = −13 to −11 and D = −5. There are also positive correlations at approximately D = 9 to 15. On the other hand, Figure 10 shows that more apparent positive correlations for the increase anomalies appear at approximately D = −6 to 3 and on D = 7. After the outliers with large energy releases are removed, for the decrease anomalies, positive correlations can be seen at approximately D = −11 to −7 in the four sectors. For the increase anomalies, there are apparent correlations during D = −15 to −14 and −6 to −2 in the four sectors. The results for energy releases ≥ *M*5.5 with/without outliers within Region C are provided in Figures S21–S24.

**Figure 8.** Temporal and spatial distributions of correlation coefficient between the daily seismic wave energy release ≥ *M*5.5 with the outliers within Region B and the occurrence times of the TEC decrease anomalies. Panels from left to right represent −15 to 15 days from the seismic wave energy release (D). Top to bottom panels display the four selected longitude sectors centered over 120, 30, −60, and −150◦E, respectively, each with a width of ±20◦. The contour denotes *r*<sup>α</sup> = 0.01 according to sample size at each D. The black box on the top first panel indicates Region B.

**Figure 9.** Temporal and spatial distributions of correlation coefficient between the daily seismic wave energy release ≥ *M*5.5 with the outliers within Region B and the occurrence times of the TEC decrease anomalies. Panels from left to right represent −15 to 15 days from the seismic wave energy release (D). Top to bottom panels display the four selected longitude sectors centered over 120, 30, −60, and −150 ◦E, respectively, each with a width of ±20◦. The contour denotes *r*<sup>α</sup> = 0.01 according to sample size at each D. The black box on the top first panel indicates Region C.

**Figure 10.** Same as Figure 9 for the results of occurrence times of the TEC increase anomalies.

#### *3.5. Geographical Dependence*

To further investigate the possible precursory information, the *p* value (*p*), which is defined as the probability that *t* is greater than the observed *t* (Equation (4)) with a degree of freedom of sample size–2 under the hypothesis of no correlation, was computed from the *r* values of the decrease and increase anomalies shown in Supplementary Material S2, and the minimum *p* values during D = −15 to 15 and over the global 71 × 73 (5183) GIM grids were extracted. The D, location, and *r* of the minimum *p* were also extracted in order to examine the possible geographical dependence (Table 3). The results reveal that the minimum *p* of the decrease anomaly is D = −13 in the 120◦E sector and approximately 30–35◦N, when all samples available from 110–130◦E and Region A and those with energy releases ≥ *M*5.5 with outliers in Region A were used. In addition, the minimum *p* of the decrease anomaly appears at D = −6 in the 120◦E sector for energy releases ≥ *M*5.5 with outliers in Region B. For energy releases in Regions A and B, the minimum *p* values, which appear before the energy release and in the 120◦E sector, are not above the latitude ranges of the two regions. We search for the minimum *p* and wish to find the time and locations where the SIPs are most likely to appear, and also expect that it would be more efficient in earthquake prediction by using GIM.


**Table 3.** The minimum *p*, associated D, locations, and *r* in analyses of this study shown in Supplementary Material S2.

Note: The minimum *p* during D = −15 to 15 and over the 5183 GIM grids were extracted for occurrence times of the decrease anomalies and increase anomalies, respectively. The bold font denotes the minimum *p* appearing before the energy release and in the 120◦E sector. The asterisk near *r* denotes that the associated *r* is the maximum *r* during D = −15 to 15 and over the 5183 GIM grids.

#### **4. Discussions**

First, the positive correlations that appear before the earthquake energy release are summarized. For all the samples of energy released within 110–130◦E longitude during the years 2015–2018, the most apparent correlations for decrease anomalies are around D = −6 to −1 over 120 and 30◦E sectors, and this pattern is similar to the results of Region B. The minimum *p* of the decrease anomalies is found at D = −13 in the 120◦E sector. One of the features of Region A is that positive correlations of decrease anomalies occur in the northern, higher latitudes. The minimum *p* of the decrease anomalies is also found at D = −13 in the 120◦E sector when using all samples and those with energy releases ≥ *M*5.5 with outliers. In Region B, the apparent correlations of the decrease anomalies appear at approximately D = −6 to −1 in the longitude sectors of 120 and 30◦E, as obtained with all samples and those with energy releases ≥ *M*5.5. These should be produced by the outliers with large energy releases. The minimum *p* of decrease anomalies is found at D = −6 in the 120◦E sector for energy releases ≥ *M*5.5 with outliers. For Region C, apparent correlations of the decrease anomalies are seen in all four sectors for energy releases ≥*M*5.5 with and without outliers. Positive correlations of the increase anomalies appear a few days before energy releases for energy releases ≥ *M*5.5 with and without outliers.

Since positive correlations appear in the four sectors and after energy release, the possible precursory information cannot be confirmed. In addition, the minimum *p* values occurring before energy release and in the 120◦E sector are characterized by samples that have not been declustered and have outliers. In Supplementary Material S3, we first examine the results of the declustered samples with energy releases of ≥*M*5.5 and without outliers. The advantage of the declustering process is that is avoids the influence of consecutive larger energy releases during a space weather event. The declustering method was as follows: if there were two energy releases of ≥ *M*5.5 separated by ≤ 6 days, the larger energy release was retained. Moreover, two methods were used to eliminate space weather effects. The first one involved the Dst index, which was identical to that used in the paper. The second method involved considering Kp and F10.7. The condition that for days with at least one Kp > 3 (3o), the occurrence times of the TEC anomalies were not considered as set. In addition, five periods of strengthening solar activity were found by using F10.7 and excluded (please see Supplementary Material S3 for details).

In Supplementary Material S3, Figures S1–S12 show that, for Region A (Figures S1–S4), one obvious characteristic of the decrease anomalies is that a large area of positive correlation appears at D = −2 around Region A when using Dst (the odd figures) and using Kp and F10.7 (the even figures), which are also provided in Figures 11 and 12, respectively. The apparent correlations of the increase anomalies mainly appear at D ≤ −4. When applying Kp and F10.7, the minimum *p* of the decrease anomalies appears at D = −2 and in the 120◦E sector at −20◦N (Table 4; Table S2). For Region B, the precursory features of the decrease and increase anomalies seem difficult to identify when using Dst and using Kp and F10.7, since the correlations occur in the four sectors and after the energy release (Figures S5–S8), and no minimum *p* of the decrease or increase anomalies appears before the energy release in the 120◦E sector. For Region C, some of the temporal characteristics of the positive correlation of decrease and increase anomalies are closer to the results than the samples without the declustering process and including the outliers were used (Figures S9–S12), as shown in Figures 9 and 10. When applying Kp and F10.7, the minimum *p* of the increase anomalies appears at D = −1 in the 120◦E sector; however, at D = −1, larger areas of correlation are also observed in the four sectors in the distributions of Figure S12. The parameters for the declustered samples of energy release ≥ *M*5.5 without outliers are shown in Table S1 (similar to Table 2). The results of the declustered samples of energy release ≥ *M*5.5 and with outliers are provided in Supplementary Material S4 (Figures S1–S12 and Tables S1–S2). For daily seismic wave energy released within Regions A, B, and C with the use of the Dst index to eliminate the space weather effects and within Region C with the use of Kp and F10.7, the results for the decrease and increase anomalies are similar to those of energy releases ≥ *M*5.5 with outliers in the paper and Supplementary Material S2. When applying Kp and F10.7, the correlations of the decrease anomalies improve between D = −13 and −6 in Region A (Figure S2), while in Region B, frequent positive correlations can be seen in the northern higher latitudes at approximately D ≤ 2 (Figure S6). On the other hand, there is no explicit precursory information in the correlation coefficients of the increase anomalies for Regions A and B (Figures S4 and S8). Only in Region A with the use of the Dst index does the minimum *p* of the decrease anomalies occur on D = −13 in the 120◦E sector (Tables 5 and S2).

Furthermore, the results of the declustered samples with energy releases ≥ *M*5.5 over the entire area of Regions A–C (13◦S–32.1◦N) are also examined in Supplementary Material S3. The results of retaining and excluding the outliers, as well as using the same two methods to remove the space weather effects, are presented in Figures S13–S20. In addition, the results might show a large variability in the occurrence times of the TEC anomalies when using *k* = 1.5 in the definition of TEC anomalies, and *k* = 3 (a probability of 0.043 for an increase or decrease anomaly under a normal distribution) was used for the energy release within the entire area of Regions A–C. Figures S13–S16 show the results of retaining the outliers. The distributions of the decrease anomalies and increase anomalies are similar to or resemble some of the features found in Region B, as shown in Supplementary Materials S2 and S4. After removing the outliers, more sporadic correlations are shown in Figures S17–S20, and the minimum *p* values of the increase anomalies are found at D = −4 and −15 in the 120◦E sector at −47.5◦N and −72.5◦N when using Dst and using Kp and F10.7, respectively (Table S4). Nevertheless, in the associated distributions of Figures S19 and S20, correlations also appear at higher latitudes in other sectors. The parameters for the declustered samples with energy releases ≥ *M*5.5 over the entire area of Regions A–C are shown in Table S3.

**Figure 11.** Distributions of the correlation coefficient of the TEC decrease anomaly for the declustered samples with energy releases ≥ *M*5.5 within Region A and without outliers. Data are removed using the Dst index.

**Figure 12.** Same as Figure 11, but data are removed using Kp and F10.7.


**Table 4.** The minimum *p*, associated D, locations, and *r* for the results of the declustered samples with energy releases ≥ *M*5.5 without outliers shown in Supplementary Material S3.

Note: The minimum *p* during D = −15 to 15 and over the 5183 GIM grids were extracted for occurrence times of decrease anomalies and increase anomalies, respectively. The bold font denotes the minimum *p* appearing before energy release and in the 120◦E sector. The asterisk near *r* denotes that the associated *r* is the maximum *r* during D = −15 to 15 and over the 5183 GIM grids.

**Table 5.** Same as Table 4 for the results of the declustered samples with energy releases ≥ *M*5.5 with outliers shown in Supplementary Material S4.


Note: The minimum *p* during D = −15 to 15 and over the 5183 GIM grids were extracted for occurrence times of decrease anomalies and increase anomalies, respectively. The bold font denotes the minimum *p* appearing before energy release and in the 120◦E sector. The asterisk near *r* denotes that the associated *r* is the maximum *r* during D = −15 to 15 and over the 5183 GIM grids.

Since the correlations appear frequently, to examine whether the possible precursory information, such as the minimum *p* appearing before an energy release in the 120◦E sector, could be accidental, we performed three times of random tests for Region A. The time series of daily seismic wave energy releases for 2015–2018 within Region A were permuted randomly. Then, the dates with energy release ≥ *M*5.5 were selected, and the declustering process and removal of outliers were implemented. The results show that regardless of using Dst or using Kp and F10.7, there are positive correlations of TEC anomalies in the four sectors and after energy release. The minimum *p* values are all quite small, the same as those shown in Table 4, and the minimum *p* values could also appear in the 120◦E sector and before energy release. These indicate that the precursory information in this study might be accidental. The results of random permutations are shown in Figures S1–S12 and Tables S1–S2 in Supplementary Material S5.

Although the precursory information could be coincidental, some interesting results might be seen from this study. For example, the results of Region B show that larger areas of correlation appear before energy release when including the outliers, and the minimum *p* of decrease anomaly for energy releases ≥ *M*5.5 with outliers occurs on D = −6 and in the 120◦E sector (Table 3). This possibly indicates an SIP event has appeared with other ionospheric disturbances. Another interesting result is that for Region A, the minimum *p* values of the decrease anomaly are on D = −13 and in the 120◦E sector when including the outliers and using the Dst index, and the minimum *p* appears on D = −2 and in the 120◦E sector for the declustered samples with energy releases ≥ *M*5.5 without outliers when using Kp and F10.7. This possibly indicates that the SIPs might not only appear closer to the epicenter region, but also occur along the earthquake longitude.

Linear regression models were fitted to samples for which the minimum *p* values appear in the 120◦E sector and before energy release in Tables 3–5 and Table S4 in Supplementary Material S3. The results are shown in Figure 13. The top panels, from left to right, are the results of using all samples, energy releases ≥ *M*5.5 with outliers, and declustered samples with energy releases ≥ *M*5.5 with outliers for Region A, respectively. The time and locations of samples, anomaly types, and the methods used to eliminate the space weather effects are indicated in the title. When using all samples, the R<sup>2</sup> is merely 0.0244 (= *r*2), because when the seismic energies are smaller, there are still larger occurrence times of the TEC decrease anomaly. After removing samples with energy releases < *M*5.5, the R<sup>2</sup> values increase in the other two panels. There are four energy releases ≥ *M*6.9 with the occurrence times of the decrease anomaly ≥ 12 times on D = −13 and at (32.5◦N, 125◦E). The result of declustered samples with energy releases ≥ *M*5.5 without outliers within Region A, and using Kp and F10.7, shows that the R2 is the greatest (0.7168) among six panels on D = −2 and at (–20◦N, 135◦E) for the decrease anomaly (bottom left panel). For energy releases ≥ *M*5.5 with outliers within Region B, there is an energy release of approximately *M*7.3 that has 20 occurrence times of the decrease anomaly (bottom middle panel). The regression line is affected by the sample, with an energy release of approximately *M*7.3. However, since the possible precursory information shown in this study cannot be confirmed to be related to the earthquakes, the regression models can only be provided as references.

Because the CODE GIM has provided TEC maps of 1-h resolution since 19 October 2014 (day of year 292), the period around 2015–2018 was selected as the study period. Earlier and later years of data might be included in the analysis to include a greater number of earthquakes and larger earthquakes. Moreover, the earthquake depth is one of the factors that possibly influences SIPs. The earthquake depth could be considered in future studies. Since there are positive correlations at higher latitudes, removing possible high-latitude effects might also be considered.

The TEC derived from the Jet Propulsion Laboratory (JPL) GIM was used to examine the daily occurrence times of TEC anomalies during a 30-day period from 21 July to 19 August 2018. The time resolution of JPL GIM is 2 h. During this period, there were six *M* ≥ 6.0 earthquakes occurring on the four days, and two events with *M*6.9 among the six earthquakes, which occurred on 5 August and 19 August, respectively. Figure 14 shows the temporal variations of Dst, Kp, and F10.7 between 16 July and 23 August 2018. The Kp index reveals that the geomagnetic disturbances with Kp ≥ 3o occurred around 17, 21, and 24–25 July and 7, 11, and 15–20 August. The Dst variations show that a weaker magnetic storm occurred on 14 August. Figures 15 and 16 show the daily occurrence times of the TEC decrease anomaly and increase anomaly from JPL with *k* = 1.5, respectively. Frequent ionospheric positive effects are observed in the associated periods with Kp ≥ 3o at higher latitudes and over the globe, which are followed by the ionospheric negative effects. These effects should be associated with the magnetic disturbances. On the other hand, despite the ionospheric negative effects over the northern higher latitudes on 23 July (D = −13 for the *M*6.9 earthquake on 5 August), larger occurrence times of the decrease anomaly also appear in the southeast of Region A and the associated conjugate point. In addition, on 6 August (D = 1 and D = −13 for another *M*6.9 event on 19 August), greater occurrence times of the decrease anomaly appear between approximately 120 and 180◦E in the northern mid-latitudes. It is noted that the greater occurrence times of the increase anomaly are seen between 29 July and 4 August (D = −7 to −1) around Region A. In fact, the TEC increase and decrease anomalies with *k* = 1.5 tend to appear globally. The larger occurrence times of TEC anomalies observed prior to the earthquakes might be only the natural variability; however, there is also the possibility of existence of the earthquake effects.

**Figure 13.** Linear regression models for samples with the minimum *p* values appearing in the 120◦E sector and before the energy release. The time and locations of samples, anomaly types, and the methods used to eliminate the space weather effects are indicated in title. The sample size (N), *r*, equation of best fit line (x: occurrence times of TEC anomaly, y: daily seismic wave energy release in erg), RMSE (in erg and magnitude), and R2 are shown in each panel. The unit of y-axis is erg, and the values and ticks are displayed with the corresponding magnitudes.

**Figure 14.** Temporal variations of Dst, Kp, and F10.7 during 16 July to 23 August 2018. The blue lines denote the time of six earthquakes during the period (1. *M*6 at 17:07:23 on 28 July, 7.10◦S, 122.73◦E. 2. *M*6.4 at 22:47:38 on 28 July, 8.24◦S, 116.51◦E. 3. *M*6.9 at 11:46:38 on 5 August, 8.26◦S, 116.44◦E. 4. *M*6.5 at 15:35:01 on 17 August, 7.37◦S, 119.80◦E. 5. *M*6.3 at 04:10:22 on 19 August, and 8.34◦S, 116.60◦E. 6. *M*6.9 at 14:56:27 on 19 August, 8.32◦S, 116.63◦E.). The time of the earthquake catalog is in UTC.

**Figure 15.** Daily occurrence times of the TEC decrease anomaly from JPL. The time resolution of JPL GIM is 2 h. The UT date is noted on top of each panel. The white box denotes Region A, and the stars denote the earthquakes occurring on the four days.

**Figure 16.** Same as Figure 15 for the daily occurrence times of the TEC increase anomaly.

Due to the conjugate effects of the ionospheric variations observed prior to the earthquakes, the electric field perturbations could be important. The origin of the electric field perturbations might be associated with the neutral wind modification. For example, Oyama et al. [17] proposed that the internal gravity waves (IGWs) of extremely small amplitude produced by small-amplitude ground motion before a large earthquake could be one of the possible origins. The IGWs interact with planetary-scale waves below 10 km and are amplified. When the amplified IGWs propagate to the dynamo region, they modify the

wind system or conductivity and, in turn, modify the electric field. Due to the frequent occurrences of ionospheric anomalies, the earthquake signals must constantly interact with the ionosphere. In addition, the statistical results in this study show a possibility of earthquake effects at higher latitudes; therefore, the area of neutral wind modification should be quite large. To understand the ionospheric variability, it is required to apply more parameters to the analyses and investigate the possibility of connections with earthquakes.

#### **5. Conclusions**

In conclusion, the temporal and spatial correlation between the seismic wave energy releases within the three low geomagnetic latitude regions in 110–130◦E and the occurrence times of the ionospheric TEC anomalies have been examined. Four longitude sectors of 120, 30, −60, and −150◦E are investigated, and there is no clear difference in the distributions of correlation coefficients between 120◦E and the other sectors. Although the statistical results do not show sufficient evidence for the specific SIPs for the three study regions, the greater occurrence times of the TEC decrease anomaly are found on D = −13 prior to two *M*6.9 earthquakes in August 2018 in the southeast of Region A and the conjugate point by using JPL GIM, according to the results of statistical investigations. Therefore, it is required to apply more parameters to analyze the causes of the ionospheric TEC variations and investigate the possible connections with earthquakes.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/rs14051155/s1, Supplementary Material S1 (Figure S1): Dst variations, temporal distributions of earthquake events and variations in daily total seismic wave energy release during 2015–2018 for earthquakes and energy releases within 110–130◦E, Region A, Region B, and Region C, Supplementary Material S2 (Figures S1–S24): results of correlation analyses shown in Section 3, Supplementary Material S3 (Figures S1–S20, Tables S1–S4): results of correlation analyses for the declustered samples with energy releases ≥ *M*5.5 without outliers for Regions A, B, and C, and for the declustered samples with energy releases ≥ *M*5.5 from the entire area of Regions A–C, Supplementary Material S4 (Figures S1–S12, Tables S1–S2): results of correlation analyses for the declustered samples with energy releases ≥ *M*5.5 with outliers for Regions A, B, and C. Supplementary Material S5 (Figures S1–S12, Tables S1–S2): results of correlation analyses for the declustered samples with energy releases ≥ *M*5.5 without outliers from three times of random permutations of energy releases within Region A.

**Author Contributions:** Conceptualization, Y.-C.S. and J.S.; methodology, Y.-C.S.; software, Y.-C.S.; validation, Y.-C.S. and J.S.; formal analysis, Y.-C.S.; investigation, Y.-C.S.; resources, J.S.; data curation, Y.-C.S.; writing—original draft preparation, Y.-C.S.; writing—review and editing, Y.-C.S.; visualization, Y.-C.S.; supervision, J.S.; project administration, J.S.; funding acquisition, J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Multigovernment International Science and Technology Innovation Cooperation Key Project of the National Key Research and Development Program of China (grant number 2018YFE0184300), EU Erasmus+ Project "Innovation on Remote Sensing Education and Learning (IRSEL)" (586037-EPP-1-2017-1-HU-EPPKA2-CBHE-JP), EU Erasmus+ Project "GIS and Remote Sensing for Sustainable Forestry and Ecology (SUFOGIS)" (598838-EPP-1-2018-EL-EPPKA2-CBHE-JP), Fujian Provincial Department of Science and Technology (2018I0005), and Fujian Provincial Department of Human Resources and Social Security (2017 strait postdoctoral exchange funding project).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The GIM data is obtained from the Center for Orbit Determination in Europe (CODE) (http://www.aiub.unibe.ch/download/CODE, last accessed on 24 February 2022). The earthquake catalog is obtained from the U.S. Geological Survey (https://earthquake.usgs.gov/ earthquakes/search/, last accessed on 26 January 2022). The Dst index is obtained from the World Data Center for Geomagnetism in Kyoto (http://wdc.kugi.kyoto-u.ac.jp/dstdir/, last accessed on 24 January 2022). The results presented in this paper rely on the Kp index, calculated and made available by ISGI Collaborating Institutes from data collected at magnetic observatories. We thank

the involved national institutes, the INTERMAGNET network, and ISGI (isgi.unistra.fr, last accessed on 30 January 2022). The F10.7 is obtained from Natural Resources Canada (https://spaceweather.gc. ca/solarflux/sx-5-en.php, last accessed on 30 January 2022). The GIM TEC data published by JPL is downloaded from the Crustal Dynamics Data Information System (ftp://gdc.cddis.eosdis.nasa.gov/, last accessed on 26 January 2022).

**Acknowledgments:** The authors are grateful to J.Y. Liu (Department of Space Sciences & Engineering, National Central University) for his guidance. The authors are grateful to P.K. Rajesh (Department of Earth Sciences, National Cheng Kung University) for his help in correcting the manuscript. The authors are also deeply grateful to the anonymous reviewers for their critical comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

