Suspected Seismo-Ionospheric Anomalies before Three Major Earthquakes Detected by GIMs and GPS TEC of Permanent Stations

Many studies have reported that there is a coupling mechanism between ionosphere and earthquake (EQ). Ionospheric anomalies in the form of abnormal increases and decreases of ionospheric Total Electron Content (TEC) are even regarded as precursors to EQs. In this paper, TEC anomalies associated with three major EQs were investigated by Global Ionospheric Maps (GIMs) and GPS-TEC, including Kumamoto-shi, Japan—EQ occurred on 15 April 2016 with Mw = 7.0; Jinghe, China—EQ occurred on 8 August 2017 with Mw = 6.3; and Lagunas, Peru—EQ occurred on 26 May 2019 with Mw = 8.0. It was found that the negative ionospheric anomalies linger above or near the epicenter for 4–10 h on the day of the EQ. For each EQ, the 10-min sampling interval of TEC was extracted from three permanent GPS stations around the epicenter within 10 days before and after the EQ. Variations of TEC manifest that the negative ionospheric anomalies first appear 10 days before the EQ. From 5 days before to 2 days after the main shock, the negative ionospheric anomalies were more prominent than the other days, with the amplitude of negative ionospheric anomaly reaching −3 TECu and the relative ionospheric anomaly exceeding 20%. In case of Kumamoto-shi EQ, the solar-geomagnetic conditions were not quiet (Dst < −30 nT, Kp > 4, and F10.7 > 100 SFU) on the suspected EQ days. We discussed the differences between ionospheric anomalies caused by active solar-geomagnetic conditions and EQ. Combining the analysis results of Jinghe EQ and Lagunas EQ, under quiet solar-geomagnetic conditions (Dst > −30 nT, Kp < 4, and F10.7 < 100 SFU), it can be found that TEC responds to various solar-geomagnetic conditions and EQ differently. The negative ionospheric anomalies could be considered as significant signals of upcoming EQs. These anomalies under different solar-geomagnetic conditions may be effective to link the lithosphere and ionosphere in severe seismic zones to detect EQ precursors before future EQs.


Introduction
Since the study of Leonard and Barnes [1] discovered that there is a potential coupling mechanism between earthquake (EQ) and abnormal ionospheric disturbance from the Alaska EQ in 1964, several researches have been devoted to detect short-and long-term seismo-ionospheric anomalies from ionospheric remote sensing measures. Variations in electrical resistivity; the greatest plasma frequency in the ionosphere, foF2; electromagnetism (EM); and ionospheric Total Electron Content (TEC) were correlated with large amounts of energy released during the EQ incubation period, which are probably signs of a impending main shock [2][3][4][5]. TEC can be obtained by dual-frequency, satellite-based measurements from Global Positioning System (GPS) [6]. A large number of published studies have reported that anomalous TEC variations appear either as precursory effects from a few days to a few weeks before an EQ or as late effect around the main shock [7,8]. For example, during the 27-day period before the EQ, negative and positive anomalies appeared in the ionosphere, which were not related to active solar-geomagnetic conditions [9]. Statistical analysis of 10 days of TEC data before global Mw ≥ 6.0 EQs during 1998-2014 showed that within 5 days before the EQs, the positive and negative ionospheric anomalies were the most prominent [10]. Similarly, significant ionospheric TEC depletions were observed 3-4 days before the Chi-Chi EQ-Mw = 7.6-occurred in Taiwan and were different from the equatorial positive anomalies, which were placed on each side of the equator [11]. However, local conditions of the ionosphere are subject to numerous influences such as solar radiation, geomagnetic activity, meteorological events, anthropogenic effects, atmospheric gravity waves (AGW), and traveling ionospheric disturbances (TID) [12][13][14]. When the surroundings of EQ suffer from abnormal space weather effects, with GPS it is difficult to accurately detect the ionospheric TEC anomaly information [12,15]. Thomas et al. [16] even strongly stated that it is impossible to take pre-earthquake ionospheric anomalies as EQ precursor information according to long-term statistical analysis.
Nevertheless, several published researches have reported that solar-geomagnetic activity and other inducements of ionospheric anomalies lead to variability in the ionosphere. There are temporal and spatial characteristics of ionospheric anomalies induced by EQ, which provides an opportunity for us to study ionospheric precursors [15,17,18]. To explain the characteristics, the process of EQ-affected ionosphere has been widely reported. Tectonic forces inside the Earth are the main cause of EQ. The highly active free-electron carrier is activated in the squeezed rocks within seismogenic zones, and they propagate upwards. The rise of conductivity in the air induced by the active free electron, eventually causing ionospheric perturbations [19]. Similarly, Kuo et al. [20] proposed that the deformation of the earth's surface will cause a change in the magnetic field, and the resulting electric field will have a significant impact on the ionosphere. In other studies, Pulinets et al. [12] proposed a Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) model to explain variations in the atmosphere and ionosphere parameters during the seismic incubation period, showing that the decrease or increase in the thermal conductivity of atmosphere is caused by continuous Radon gas emission from the EQ breed zone.
In this paper, the research of seismo-ionospheric anomalies before three earthquake cases in different regions, with different focal depths and under different solar-geomagnetic conditions was implemented. We analyzed the temporal and spatial distribution of ionospheric anomalies associated with three Mw > 6.0 EQs on the day of each EQ. The temporal series of Vertical Total Electron Content (VTEC) was retrieved from permanent GPS stations within the EQ breed zone from 10 days before to 10 days after the EQ. Sliding interquartile range method was used to distinguish seismo-ionospheric anomalies specifically triggered by the impending main shock. The ionospheric TEC variation related to solar-geomagnetic conditions were discussed in detail with enough evidences. Section 2 introduces the materials and methods; then, Section 3 shows our results and description. In Section 4, we summarize the characteristics of seismo-ionospheric anomalies. The process of EQ-induced ionospheric anomalies is also discussed in detail with previous studies. Finally, in Section 5, our conclusions are provided.

Earthquake Data
In this paper, we studied the temporal and spatial TEC anomalies related to three Mw > 6.0 EQs, including the Kumamoto-shi, Japan EQ that occurred on 15 April 2016 with Mw = 7.0; the Jinghe, China EQ that occurred on 8 August 2017 with Mw = 6.3; and the Lagunas, Peru EQ that occurred on 26 May 2019 with Mw = 8.0. They occurred in different locations around the world and under various geomagnetic conditions. EQ occurrence time is converted to Coordinated Universal Time (UTC). The detailed information of these three EQs retrieved from the U.S. Geological Survey (USGS) via http://www.earthquake. usgs.gov/earthquakes/search (accessed on 5 November 2021) is shown in Table 1.

GIMs and GPS-TEC Data
Many scientists have utilized GIM containing grid data of the Vertical Total Electron Content (VTEC) to investigate ionospheric phenomena. The GIM is generated by data collected from almost 400 stations distributed around the world, with a temporal resolution of 1 h and spatial resolution of 2.5 • × 5 • (lat. × lon.) [23]. Preseismic ionospheric anomalies were investigated by GIMs acquired from the Center for Orbit Determination in Europe (CODE) via http://ftp.aiub.unibe.ch/CODE/ (accessed on 5 November 2021), which is one of the most reliable Data Products when monitoring ionosphere [24]. The VTEC retrieved from permanent GPS ground-stations has higher accuracy and temporal resolution. For each EQ, three available GPS stations were selected by the method proposed by Dobrovolsky et al. [25]. The EQ breed zones were estimated by the following equation: where M w is EQ magnitude and R is the EQ breed zone's radius. In this study, all the stations operated within the EQ breed zone are shown in Figure 1. RINEX data of permanent GPS observatories in China (XJDS, XJXY, XJYN) were obtained from the Crustal Movement Observation Network of China (CMONOC), which includes GNSS stations across mainland China and enables the continuous monitoring of the ionosphere over China as accurately as possible [26]. The rest (AIRA, DAEJ, SMST, BOGT, POVE, RIOP) were obtained from International GNSS services (IGS) observatories via http://www.igs.gnsswhu.cn/index.php/Home/DataProduct/igs.html (accessed on 5 November 2021). We extracted the ionospheric VTEC of each site with 10 min sampling intervals. According to the research of Liu et al. [11] and Shah and Jin [10], the 10 days before and after the EQ were selected as EQ-related days. The main part of the total ionospheric delay or Slant Total Electron Content (STEC) in the propagation direction of the GPS signal can be extracted from dual-frequency GPS measurements [6] and expressed as the following equation [27]: where f 1 , f 2 stand for the frequency of the carrier; P 4,sm can be obtained by using the carrier phase observations to smooth the pseudorange observations; c is speed of light in a vacuum; and DCB i , DCB j stand for differential code biases of the satellites and differential code biases of the receivers, respectively. STEC is calculated along the ray path of the observed signal in the unit of TECu (1 TECu = 10 16 electron/m 2 ; Li et al., 2012) and retrieved from GPS stations in a one-square-meter tube. Then, STEC is transformed into VTEC by the following equation [28]: in this equation, R denotes the Earth radius and H is the elevation of ionosphere's upper limit. Similarly, Z is the satellite elevation angle [29].

Anomaly Analysis Method
The TEC was extracted at the same time each day and sorted in ascending order (the TECs mentioned below all stand for VTEC). To identify the anomalous signals in TEC measurements, the Sliding Interquartile Range (IQR) analysis method was adopted on daily TEC processing, which was proposed and implemented by Liu et al. [3]. The median (TEC median ) and the associated interquartile range (IQR) of every successive 10-day value were calculated to construct TEC anomaly bounds. Under the assumption of a normal distribution with mean µ and standard deviation σ for the TEC, the expected value of TEC median and IQR are µ and 1.34σ, respectively [3]. In this paper, we constructed the upper bound TEC median + 1.5IQR and lower bound TEC median − 1.5IQR, presented in Equation (4). These bounds are expected to be µ and 2.0σ, respectively. If the original TEC occurs outside of either bound, the variations are anomalous with almost 95% confidence that a positive (moved up beyond upper bound) or negative (moved down below lower bound) signal is detected [30]. The IQR method has advantages in eliminating gross errors.
The ionospheric abnormal value is calculated by the following equation: where ∆TEC p represents a positive ionospheric anomaly, ∆TEC n represents a negative ionospheric anomaly, ∆TEC 0 means no anomaly occurred. To determine the extent of the ionospheric anomaly, the Relative TEC (RTEC) is calculated by the following equation: When the ionospheric TEC is beyond the upper bound or below the lower bound, it is judged as an abnormal ionospheric TEC enhancement or depletion. The extent of ionospheric anomaly is measured by RTEC. However, we only calculated the RTEC of negative ionospheric anomalies to observe the extent and distribution characteristics of TEC depletion related to the main shock.

Results and Description
In this paper, the ionospheric anomalies associated with three Mw > 6.0 EQs were detected by GIMs and GPS-TEC. We calculated the ionospheric TEC anomalies at each grid point in the bihourly GIMs to show the drift trend and spatial distribution of global ionospheric anomalous clouds over time. The GPS-TEC was extracted from 3 permanent GPS stations within the EQ breed zone and the time span was from 10 days before to 10 days after these EQs. For each EQ, we examined solar-geomagnetic conditions. Generally, when establishing the ionospheric anomaly detection method, it is important to confirm whether the method produces incorrect detection under normal conditions and if the method fails to produce a correct detection under abnormal conditions.

Kumamoto-shi Earthquake
From 8 to 17 April, the solar-geomagnetic conditions were not quiet. It can be seen from Figure 2 that the Kp and Dst index are unstable during this period. On 7, 12, 14, 16 April, a sharp increase appears on the Kp index and the value of Kp even exceeds 4. At the same time, a sharp decrease appears in the Dst index, which means active geomagnetic conditions. In addition, during these 10 days, the solar radiation index F10.7 exceeds 100 SFU. On 15 April 2016, a Mw = 7.0 EQ occurred as a result of strike-slip faulting at shallow depth, with a focal depth of 10 km. As shown in Figure 3, the bihourly GIMs reveal ionospheric anomalies in different locations around the world. At UT = 0 h of the EQ day, significant negative ionospheric anomalies appear over the northeast of the epicenter, which linger until UT = 4 h, and the amplitude of TEC anomaly is about -5 TECu. Then, the anomalous clouds gradually drift southwestward. From UT = 0-8 h, there is a significant positive ionospheric anomaly near the equator, which appears to be more in line with the intensification of Equatorial Ionization Anomaly (EIA). Starting from UT = 12 h, 12 h before the main shock, there are continuous negative ionospheric anomalies lingering in the northeast of the epicenter until UT = 22 h, and the anomaly amplitude is about −3 TECu. The TEC anomalies detected by the correlation analysis before the EQ have some peculiar characteristics. Firstly, the ionospheric anomalies always lingered over or near the epicenter, unlike diurnal changes in the ionosphere or TID, which drift westward at a speed of 100-1000 m/s. According to the statistical analysis of TIDs, the occurrence rate of TID strongly depends on the season and Local Time (LT), and in these regionsdawn (05:00-07:00 LT) and dusk (17:00-20:00 LT) in summer (May-August), and daytime (08:00-12:00 LT) in winter (November-February)-TID has the highest probability of being discovered. [31,32] Secondly, the range and extent of ionospheric anomalies are related to the magnitude of the EQ and the depth of the focal point. The earthquake occurred at 1:25 am (LT) on 15 April (spring), which is not in the high-occurrence period of TID. In addition, the ionospheric cloud over the epicenter that appeared on 15 April did not have the drift characteristics of TID. So, the persistent, slightly negative ionospheric anomaly northwest of the epicenter that lingered from UT = 12-22 h may be related to the impending main shock. We compared the ionospheric anomalies presented by GIMs on 13-15 April and analyzed the similarities and differences among them to verify that the occurrence of ionospheric anomalies was related to EQ rather than active solar-geomagnetic conditions (Dst < −30 nT, Kp > 4, and F10.7 > 100 SFU) on 15 April. The solar-geomagnetic conditions on 13 April were very similar to the main shock day. As can be seen in Figure 4, from UT = 0-4 h, significant positive ionospheric anomalies appear globally, without spatial distribution characteristics related to the impending main shock. Starting from UT = 6 h, the global ionospheric anomalies begin to fade away. At UT = 12 h, significant positive and negative ionospheric anomalies appear again on both sides of the equator. The positive and negative ionospheric anomalies are distributed globally, and their characteristics are similar to those induced by active solar-geomagnetic conditions. There is a time lag between the ionospheric anomalous response and the solar-geomagnetic conditions, as they are not exactly simultaneous [33]. Referring to Figure 5, a similar phenomenon appears on 14 April. At UT = 0 h, a slightly positive ionospheric anomaly appears near the epicenter and disappears soon after. Starting from UT = 12 h, significant positive ionospheric anomalies begin to appear over the epicenter and drift to the southwest. Suddenly intensified positive ionospheric anomalous clouds of ionosphere appear on both sides of the equator at UT = 18-20 h.
In general, the ionospheric anomalies caused by active solar-geomagnetic conditions are different in magnitude and scope from the ionospheric anomalies caused by the EQ. Otsuka et al. proposed that the drift trend from northeast to southwest and the shape of ionospheric anomalous clouds are consistent with the previous study of TID over Japan induced by active solar-geomagnetic conditions [34]. The significant positive ionospheric anomalies that appeared on 13-14 April are attributed to active solar-geomagnetic conditions. By comparing Figures 4 and 5 with Figure 3, we can find that the ionospheric anomalies in the GIMs are completely different. As a whole, the anomalous clouds induced by EQ on 15 April stay near the epicenter without drifting. The TEC anomalies before large EQs do not propagate as a circular TID, because the mechanism is radically different from the mechanism of TIDs [18].  In order to further explore the preseismic anomaly of the ionosphere, the ionospheric anomalies were studied by 10 min-sampling-interval VTEC from three GPS stations. The VTEC time-series and its associated lower and upper bounds are shown in Figure 6. It can be seen that intermittent ionospheric TEC enhancement and depletion are observed. The pre-ionospheric enhancement on 14 April, reaching 5 TECu, is more likely to be caused by the active solar-geomagnetic conditions. This is consistent with the phenomena presented in the GIMs (Figure 5), which provide evidence to justify our statistical bounds manifesta-tion. The earliest negative ionospheric anomaly appeared 8 days before the EQ, observed by three GPS stations. Only from two days before the EQ to the EQ day, there was a prominent eruption of negative ionospheric anomalies, reaching −3 TECu. Especially, on 15 April 2016-the main shock day-significant ionospheric negative anomalies were observed by three GPS stations. After the EQ, intermittent positive and negative ionospheric anomalies were also discovered.  Table 1 for EQ details). The upper bound (UB) and lower bound (LB) are shown in gray, and TEC is shown in red. The blue curve represents ∆TEC.
Moreover, the extent of relative negative ionospheric anomalies calculated by Equation (6) are shown in Figure 7. It is found that the negative ionospheric anomalies are particularly obvious on the day of the EQ. The abnormal ionospheric TEC depletions before the EQ are likely to be caused by the process of EQ incubation. This phenomenon is consistent with Figures 6 and 7, showing that significant negative ionospheric anomalies were detected on 15 April (earthquake day), including GIMs and VTEC from GPS stations.

Jinghe, China Earthquake
Variations in the solar-geomagnetic conditions are obstacles to accurately identify ionospheric anomalies before EQs. We selected the EQ that occurred on 8 August 2017 in Jinghe, China as the research object under quiet solar-geomagnetic conditions (Dst > −30 nT, Kp < 4, and F10.7 < 100 SFU), as shown in Figure 8. Kp 07/29 07/30 07/31 08/01 08/02 08/03 08/04 08/05 08/06 08/07 08/08 08/09 08/10 08/11 08/12 08/13 08/14 08/15 08/16 08/17 08/18  The epicenter is far from the ocean and at a higher latitude than the Kumamoto-shi EQ; so, the ionosphere above the epicenter is not easily affected by the EIA. The bihourly GIMs (Figure 9) show that at UT = 4 h, a slightly positive ionospheric anomaly appears close to the equator. The abnormal value gradually increases, reaching a peak of 5 TECu at UT = 6 h and drifting westward. At UT = 12 h, a positive conjugate ionospheric anomaly appears on the south side of the equator. The characteristics of such anomalies are consistent with the TID. At the same time, 9 h before the EQ, a slightly negative ionospheric anomaly begins to appear near the epicenter and far from the equator. At UT = 14 h, it hits a peak, which is lower than −3 TECu. These negative ionospheric anomalous clouds linger until UT = 20 h, 3 h before the EQ. There is no significant drift and the clouds remain above the epicenter. We investigated the solar-geomagnetic conditions when the EQ occurred, and attributed these anomalous clouds of TEC to the impending main shock under quiet solar-geomagnetic conditions, as it is different with characteristics of two-crest equatorial region. The temporal series of GPS TEC extracted from three permanent GPS stations is shown in Figure 10. It can be found that the slightly negative ionospheric anomalies first appear on 29 July, 10 days before the EQ. Then, from two days before to the day of the EQ, significant negative ionospheric anomalies are observed by three GPS stations, with amplitudes of −2 TECu. On 10 August, two days after the EQ, the negative ionospheric anomaly appears again. What cannot be ignored is that during 3-5 August, positive ionospheric anomalies erupt intensely, especially peaking with 10 TECu on 4 August, 4 days before the main shock. We can see from Figure 8 that the Kp exceeds 4 on 4 and 6 August, and the Dst shows a positive sudden change, then drops sharply on 4 August. Therefore, the positive ionospheric anomalies during these days are mainly related to active geomagnetic activities (Kp > 4, Dst < −30 nT). It further proves that ionospheric TEC will respond to active geomagnetic conditions. In the case of quiet solar-geomagnetic conditions, the occurrence of these negative ionospheric anomalies are mainly related to the main shock.

Lagunas, Peru Earthquake
The previous two cases show that negative ionospheric anomalies are found near the epicenter on EQ-related days and the extent of ionospheric anomalies caused by EQs is different. However, there are no exact characteristics about seismo-ionospheric anomalies. Shah and Jin [10] claimed that GPS TEC anomalies depend on the magnitude and hypocentral depth of future EQs and can also be attributed to different fault systems. We selected the Lagunas, Peru EQ that occurred on 26 May 2019 with Mw = 8.0 for research, which is the result of normal faulting at an intermediate depth, approximately 110 km beneath the Earth's surface. During the EQ breed period, the Solar-geomagnetic conditions are quiet, with Kp < 4, −30 nT < Dst < 30 nT, and F10.7 < 100 SFU, as shown in Figure 12. Similarly, we analyzed the temporal and spatial distribution of ionospheric anomalies in the GIMs. As can be seen from Figure 13, at UT = 0 h, a slightly negative ionospheric anomaly appears above the epicenter and the peak value is −3 TECu; then, the anomalous clouds gradually drift westward. At UT = 6 h, a slightly negative ionospheric anomaly appears again over the north of the epicenter. Until UT = 10 h, after the main shock occurs, the slightly negative ionospheric anomaly begins to fade. The range and extent of negative ionospheric anomalies are small, which may be related to the deep focal depth of the EQ. At UT = 12 h, significant conjugate negative ionospheric anomalies appear on both sides of the equator and drift from east to west. The slightly negative ionospheric anomalous clouds at UT = 0-10 h may be caused by the process of EQ incubation, according to the characteristics of the seismo-ionospheric anomaly.
The TEC data from 3 GPS stations near the epicenter were extracted for research. It can be seen from Figure 14 that as early as 10 days before the EQ, slightly negative ionospheric anomalies were observed by the BOGT station. From 5 days before the EQ to the day of the main shock, the negative ionospheric anomalies begins to appear frequently, and the amplitude reaches −3 TECu. During the 10 days after the EQ, negative ionospheric anomalies appear occasionally, not as significantly as before the EQ. During the earthquake incubation period, the ionosphere could be continuously affected by the energy released from epicenter, resulting in changes in the scale of the negative ionospheric anomaly. No significant positive ionospheric anomalies were found by three GPS stations in these 21 days under quiet solar-geomagnetic conditions (Dst > −30 nT, Kp < 4, and F10.7 < 100 SFU).   Table 1 for EQ details). The upper bound (UB) and lower bound (LB) are shown in gray, and TEC is shown in red. The blue curve represents ∆TEC.
We studied the RTEC of this EQ. As shown in Figure 15, the statistical results of the three stations show that from 5 days before the EQ to the day of the EQ, significant negative ionospheric anomalies are observed. The negative ionospheric anomalies before the EQ show some temporal correlation with the main shock. Therefore, we believe that the ionospheric anomalies above the epicenter are caused by EQ. In the analysis of these three cases, we found that negative ionospheric anomalous clouds appeared 4-8 h before the main shock and lingered for 4-10 h. During the 21 EQ-related days, negative ionospheric anomalies were observed as early as 10 days before the EQ and as late as 8 days after the EQ. The specific pre/post-TEC precursor information is shown in Table 2. Table 2. Associated pre/post-TEC precursors observed by GPS stations (the negative sign represents the number of days before the EQ, and the positive sign represents the number of days after the EQ).

EQ Location GNSS Station Pre EQ Day Post EQ Day
Kumamoto-shi, Japan

Discussion
Ionospheric anomalies related to three major Mw > 6.0 EQs were analyzed in this paper. In the bihourly GIMs, we found that before the main shock, there are significant negative ionospheric anomalies over or near the epicenter, which is different from ionospheric anomalies caused by active solar-geomagnetic conditions or TID. For each EQ, during the 21 EQ-related days, we discovered some negative ionospheric anomalies related to the process of EQ incubation. Through the study of the RTEC of the three cases, a common phenomenon is found that the negative ionospheric anomalies during 5 days before to 2 days after the EQs are more prominent than the other days, under various solar-geomagnetic conditions. There are some similarities with the research results proposed by Pulinets and Davidenko [18] on the temporal characteristics of seismo-ionospheric anomalies.
In case of the Kumamoto's Mw = 7.0 EQ, which occurred on 15 April 2016, active solar-geomagnetic conditions are recorded from 8-17 April. In previous researches, the active solar-geomagnetic conditions and EQ are both considered to be main causes of ionospheric anomalies [35]. Although not all geomagnetic activities can cause anomalies in the ionosphere [4], it is difficult to conclude that the ionospheric anomalies above the Kumamoto-shi on 15 April are induced by the upcoming EQ. We analyzed the ionospheric anomalies on 13-14 April, under almost the same active solar-geomagnetic conditions (Dst < −30 nT, Kp > 4, and F10.7 > 100 SFU). It is found that a wide range of positive ionospheric anomalies appear near the equator in these two days and the characteristics of anomalous clouds are very consistent with the mesoscale ionospheric anomalies caused by active geomagnetic conditions. The shape and drift trend of ionospheric anomalies over the epicenter on 15 April is significantly different. Furthermore, only the significant negative ionospheric anomalies on 15 April were observed by three permanent GPS stations during 21 EQ-related days. Therefore, the negative ionospheric anomalies on this day were mainly caused by the EQ. Iwata and Umeno [31] proposed the same conclusion by studying the characteristics of abnormal ionospheric cloud drift over Japan.
For the Jinghe EQ, ionospheric TEC responds to variations in solar-geomagnetic conditions differently. It is found that on the day of the EQ, there was a significant negative ionospheric anomaly over the epicenter, under quiet solar-geomagnetic conditions (Dst > −30 nT, Kp < 4, and F10.7 < 100 SFU). The ionospheric TEC responding to active solar-geomagnetic conditions further prove that under quiet conditions, EQ is the cause of negative ionospheric anomalies. Under absolutely quiet solar-geomagnetic conditions, similar temporal distribution of negative ionospheric anomalies were also found in the Lagunas EQ. However, the ionospheric negative anomalies are so weak that they cannot remain above the epicenter for a long time due to the fact that energy released by the seismic fault zone needs to travel a long distance to reach the surface.
In order to further explain how EQs cause ionospheric anomalies, some studies have explained them from the perspective of physics. Friedemann and Freund [2] proposed that the collision between rocks generates high-speed charge carriers, which form a fluctuating charge cloud; this can explain the electrical signals and electromagnetic radiation associated with EQs. Similarly, the investigation of temporal and spatial measurements suggested that abrupt anomalies observed within 10 days before the 2007 Mw = 7.7 Chile EQ were considered to be the fragments of the plasma intensification process in the southern hemisphere activated by the rock compression in the quake region [36]. Kuo et al. [20] proposed a coupling model for the stressed-rock-Earth surface-charge-atmosphere system. The stressed-rock acts as the dynamo to provide the currents for the coupling system. The current in the EQ fault zone has an impact on the electric fields and currents in the atmosphere and the lower boundary of ionosphere as well as the formation of a plasma bubble. The constant current flow from the EQs epicenter can further stimulate the horizontal electric field at the bottom of the ionosphere to initiate the vertical or zonal drift of ionospheric plasma following the magnetic field direction. On the other hand, Pulinets et al. [12] demonstrated that the observed variations of the atmosphere and ionosphere parameters have a common cause-the air ionization by radon. A significant depletion in ionospheric TEC was observed in the month of earthquake occurrence by performing long-term statistics. This is consistent with the findings of this paper.
The coupling mechanism between the EQ and ionosphere is very complex and difficult to measure. It still needs more observation by ionospheric measurements from GPS and other techniques. In other studies, the observations of different ionospheric parameters is available; for example, the variation in foF2, observed by the ionosonde, has a similar tendency and is highly correlated with overhead TEC recorded by GPS before EQ. Namely, within 1-5 days before the EQ, significant depletions on foF2 and TEC were observed [11].
This has some similarities with the results of this paper. The electron temperature measurements from Langmuir Probe (ISL) of Detection of Electromagnetic Emissions Transmitted from EQ Region (DEMETER) can also be used as an indicator to judge the ionospheric anomalies. The disturbances in the electron density observed by CHAMP satellite and the concurrent electric field and Ne changes suggest that EIA intensification could be triggered by the E field disturbances over the epicenter [37]. In terms of ionospheric TEC, the ambiguity-fixed carrier-phase ionospheric observable is more accurate than the carrierphase leveled-code one [38]. The distribution and changes of electrons in the ionosphere can be observed by ionospheric tomography method in more detail [39].

Conclusions
In this paper, the seismo-ionospheric TEC anomalies associated with three Mw > 6.0 EQs were detected by GIMs and GPS-TEC. We analyzed the spatial and temporal variations of the ionospheric TEC under various solar-geomagnetic conditions. The main findings of this research are as follows: 1.
The negative ionospheric anomalies observed within 10 days before the main shock could be considered significant signals of upcoming EQs. In the analysis of the three EQ cases, we found significant negative ionospheric anomalies with ∆TEC reaching −3 TECu before the main shock. On the day of Jinghe EQ and Lagunas EQ, the negative ionospheric anomaly clouds linger over or near the epicenter for 4-10 h, under quiet solar-geomagnetic conditions (Kp < 4, Dst > −30 nT, and F10.7 < 100 SFU). In the case study of the Kumamoto-shi, Japan EQ, on the day of the EQ, the negative ionospheric anomalies are observed under active solar-geomagnetic conditions and are significantly different from those caused by the same solar-geomagnetic conditions. 2.
The negative ionospheric anomalies are more prominent during 5 days before to 2 days after the EQ than the other days. In the analysis of three EQ cases, negative ionospheric anomalies were found with RTEC exceeding 20% during this period. These ionospheric anomalies manifest a significant temporal correlation with the main shock.

3.
Ionospheric TEC responds differently to various solar-geomagnetic conditions. In the case study of Kumamoto-shi, Japan EQ, abnormal ionospheric TEC enhancement appears under active solar-geomagnetic conditions (Kp > 4, Dst < −30 nT, and F10.7 > 100 SFU) 2 days before the main shock and the abnormal amplitude reaches 5 TECu. Similar TEC variations also appear in the research of Jinghe EQ on the 4th day before the main shock with amplitude reaching 10 TECu, under active geomagnetic conditions (Kp > 4 and Dst < −30 nT). On 15 April 2016 and 17-18 August 2017, the ionospheric TEC enhancement did not occur, even though the solar-geomagnetic conditions are still active.
Finally, we discussed the other parameters that reflect the characteristics of the ionosphere and the merging process of seismo-ionospheric TEC anomalies. More detailed studies need to be performed to distinguish the ionospheric anomalies caused by EQs.