Seismo-Ionospheric Effects Prior to Two Earthquakes in Taiwan Detected by the China Seismo-Electromagnetic Satellite

Abstract

In this paper, we focused on the characteristics of the seismo-ionospheric effects related to two successive earthquakes, namely, the earthquakes in 2022 in Taitung Sea, Taiwan, China, with magnitudes (M) of 6.7 and 6.3, at 23.45° N, 121.55° E and 23.39° N, 121.52° E and with the same focal depth of 20 km, which were detected by the China Seismo-Electromagnetic Satellite (CSES). By applying the sliding interquartile range method to electron density (Ne) data acquired by the Langmuir probe (LAP) onboard the CSES and the grid total electron content (TEC) data obtained from the Center for Orbit Determination in Europe (CODE), positive anomalies were found under quiet geomagnetic conditions on 2–3 March and 8–9 March 2022—that is, 19–20 and 13–14 d before the earthquakes, respectively, and the global ionospheric mapping (GIM) TEC data suggested that anomalies may also have been triggered in the magnetic conjugate area 13–14 d prior to the earthquakes occurrences. In addition, the CSES Ne data showed enhancements 3 and 5 d before the earthquakes occurred. Furthermore, 138 earthquakes with M ≥ 5.0 that occurred in Taiwan and the surrounding region during the period February 2019 to March 2022 were statistically analyzed using the CSES Ne data. The results show that most of the Ne anomalies were positive. Moreover, the greater the earthquake magnitude, the greater the frequency of the anomalies; however, the amplitude of the anomalies did not increase with the earthquake magnitude. The anomalies were concentrated during the period of 10 d before to 5 d after the earthquakes. No increase in the amplitude of anomalies was observed as the time of the earthquakes approached. Finally, based on evidence relating to earthquake precursor anomalies, we conclude that it is possible that earthquakes in Taiwan and the surrounding region affect the ionosphere through the geochemical, acoustic, and electromagnetic channels, as described by the lithosphere–atmosphere–ionosphere coupling (LAIC) model, and that the two studied earthquakes in Taiwan may have induced ionospheric effects through the geochemical channel.

1. Introduction

Seismo-ionospheric detection is one of the new technologies used for the short-term monitoring and prediction of earthquakes, as well as an important means of studying the sphere-coupling mechanism [1]. Since the ionospheric perturbations that occurred before the 1964 M 9.2 Alaska earthquake were observed by Leonard and Barnes (1965) [2], many studies [3,4,5,6,7] have demonstrated that these ionospheric perturbations do exist and that there is a high probability (>70%) of detecting the ionospheric perturbations related to M > 6.0 earthquakes [8].

Seismo-ionospheric monitoring technologies can be divided into two types: ground-based and space-borne. Ground-based ionospheric monitoring technology includes vertical/slant ionosounding, total electron content (TEC), very low frequency (VLF) electromagnetic wave observations, and Schumann resonance detection [1]. Using vertical/slant ionosounding, Pulinets and Boyarchuk (2004) [9] presented the first comprehensive study of F2 layer critical frequency (foF2) anomalies observed before strong earthquakes. Liu et al. (2016) [10] carried out a comparative study of foF2 anomalies in northern China using an oblique sounding network and a vertical ionosonde. It was shown that these anomalies exhibited similar characteristics and could be used to investigate seismo-ionospheric effects. Ahmed et al. (2022) [11] used foF2 data to analyze the ionospheric anomalies that occurred before the 2019 M 5.6 Mirpur earthquake and reported that the anomalies appeared 10–20 d before the earthquake. The first abnormal TEC signal was detected by Calais and Minister (1995) [12] after the 1994 M 6.7 California earthquake. Zhao et al. (2008) [13] suggested that an unusually large increase in the TEC appeared 3 d before the 2008 great Wenchuan earthquake. A statistical analysis of the global ionospheric mapping (GIM) TEC data related to 56 M ≥ 6.0 earthquakes in China was performed by Chen et al. (2015) [14]. Their results demonstrated that seismic activity was responsible for the TEC anomalies. Song et al. (2018) [15] built a regional prediction model for the ionospheric TEC over China by using genetic algorithm-based neural networks (GA-NNs). Simha et al. (2020) [16] studied atmospheric and ionospheric anomalies that preceded the Taiwan earthquakes that occurred on 4 February (M 6.1) and 6 February (M 6.4), 2018. The results indicated that TEC anomalies occurred 1 d before the M 6.1 earthquake and 1, 3, and 5 d before the M 6.4 earthquake and had an absolute amplitude of 5–8 TECU. Observations of the VLF electromagnetic waves were used by Hayakawa et al. (2010) [17] to investigate the correlation between M > 6.0 earthquakes in Japan and lower ionospheric perturbations as observed in the propagation of subionospheric VLF/low frequency (LF) waves, and Zeren et al. (2014) [18] observed VLF radio signal anomalies associated with the 2010 M 7.1 Yushu earthquake and the 2013 M 7.0 Lushan earthquake. Schumann resonance detection was first used to detect anomalous phenomena in Japan before the Taiwan Chi-chi earthquake by Hayakawa et al. (2005) [19], and Ouyang et al. (2015) [20] analyzed the variability of Schumann resonance parameters observed at low-latitude stations in China.

Space-borne ionospheric monitoring technology is based on satellites and space stations. With the development of satellite observation technology, many countries have launched satellites specifically for monitoring space environment and natural disasters. Examples include the QuakeSat satellite launched by the United States, Ukraine’s SICH-1M satellite, Russia’s COMPASS-1 satellite, the European Space Agency’s SWARM satellites, the French DEMETER satellite, and China’s China Seismo-Electromagnetic Satellite (CSES). Important results have been achieved in studying the ionospheric perturbations that occur before earthquakes by using data from a range of satellites [21,22,23,24]. Until now, there has been little research based on CSES electron density (Ne) data; however, some earlier studies confirmed that these data do exhibit a clear response to large earthquakes. In 2018, Yan et al. (2018) [25] used multiple types of CSES data to analyze the seismo-ionospheric effects observed before four M > 7.0 earthquakes that occurred around the globe in August 2018. Zhang et al. [26] took the 2018 M 6.9 Indonesia earthquake as a case study for investigating VLF electromagnetic wave anomalies as well as anomalies in multiple plasma parameters and in the electromagnetic field that were recorded by the CSES. In 2020, Song et al. [27] used CSES Ne data and GIM TEC data to study the seismo-ionospheric effects observed before four M > 5.8 earthquakes that occurred in Indonesia in July and August 2018 and concluded that CSES Ne anomalies appeared up to 5 d before the earthquakes. These authors also statistically analyzed GIM TEC data for the periods before and after a total of 35 M ≥ 5.8 earthquakes that occurred during the period 2007–2017 in Indonesia and the surrounding region. The results showed that positive TEC anomalies occurred 1–7 d prior to the earthquakes. Li et al. [8] studied the seismic influence on ionospheric parameters—for example, Ne and ion density (Ni)—as recorded by the CSES and the DEMETER satellite prior to M ≥ 4.8 earthquakes. In 2021, Xie et al. [28] used multiple sources of data, including TEC data as well as CSES and SWARM Ne and electron temperature (Te) data, to detect possible ionospheric perturbations related to M 6.4 and M 7.1 earthquakes in California. It was shown that precursor anomalies occurred 5–8, 12, and 14 d before the earthquakes. Zhang et al. [29] employed CSES Ne data, TEC data from Global Positioning System (GPS) observations made on the Chinese mainland, and TEC mapping data published by the Jet Propulsion Laboratory (JPL) to analyze the ionospheric perturbations associated with the M 6.4 Yutian (China) earthquake that occurred on 26 June 2020, and suggested that the seismo-ionospheric perturbations were concentrated one week before the earthquake, mainly with increasing features. In 2022, Du and Zhang [30] used CSES Ne, Te, and ion composition data acquired before the 2021 M 6.4 Yangbi and M 7.4 Maduo earthquakes that occurred in China. It was shown that Ne and Te anomalies appeared about 40 d before these earthquakes.

The aims of this study included identifying the characteristics of seismo-ionospheric effects in Taiwan and the surrounding region, developing new methods for short-term earthquake prediction, and studying the sphere-coupling mechanism. To do this, we investigated two successive earthquakes that occurred in Taiwan (the 2022 M 6.7 and M 6.3 Taitung Sea earthquakes) using CSES Ne data. The datasets used are described in Section 2. The detailed results of the data analysis are described in Section 3. Section 4 includes a statistical analysis of the CSES Ne data associated with 138 M ≥ 5.0 earthquakes (corresponding to 65 groups) in Taiwan and the surrounding region. In Section 5, the lithosphere–atmosphere–ionosphere coupling (LAIC) mechanism is discussed in relation to earthquakes in Taiwan and the surrounding region, including the two earthquakes that were investigated in detail. The conclusion follows in Section 6.

2. Datasets

2.1. Seismic Events

Taiwan is located on the eastern edge of the Asian continent and in the subduction zone between the Eurasian plate and the Philippine Sea plate. To the east of Taiwan, along the Ryukyu trench, the oceanic crust of the Philippine Sea plate is subducting northwestward underneath the continental crust of the Eurasian plate. Meanwhile, to the south of the island, along the Manila trench, the oceanic crust of the South China Sea is subducting eastward underneath the oceanic crust of the Philippine Sea plate [31]. As a result of this complex subduction and convergence, Taiwan and the surrounding region constitute one of the most seismically active regions in the world. This seismic activity provides a natural window for observing the characteristics of seismo-ionospheric anomalies and the mechanism behind sphere coupling. Figure 1A,B illustrate the subduction of the Eurasian plate and the Philippine Sea plate in this area.

Two consecutive earthquakes (the 2022 M 6.7 and M 6.3 Taitung Sea earthquakes) occurred in Taiwan on 23 March 2022, at 23.45° N, 121.55° E and 23.39° N, 121.52° E and at 01:41 BJT (Beijing time) and 04:29 BJT, respectively. In this paper, these two earthquakes are referred to as EQ1 and EQ2, respectively. According to the China Earthquake Networks Center (CENC) catalog, the two earthquakes had the same focal depth of 20 km. EQ2 occurred only 168 min after EQ1, and the earthquake epicenters were only 7.3 km apart. Thus, in this study, we considered that, together, EQ1 and EQ2 constitute the same event. Figure 1C shows the locations of the earthquake epicenters.

Table 1. Detailed information about EQ1 and EQ2 in Taitung Sea.

No.	BJT	UT	Lat./°N	Lon./°E	Depth/km	Magnitude
EQ1	23 March 2022 01:41:38	22 March 2022 17:41:38	23.45	121.55	20	6.7
EQ2	23 March 2022 04:29:58	22 March 2022 20:29:58	23.39	121.52	20	6.3

gives details of the two earthquakes, including the time of occurrence according to Beijing time (BJT) and universal time (UT), as well as their locations, depths, and magnitudes.

2.2. Ne Data Acquired by the LAP Onboard the CSES

The CSES was launched on 2 February 2018, and follows a sun-synchronous orbit, with a height of 507 km and an inclination of 97.4°. Between latitudes 65° N and 65° S, the satellite passes over the same region on the Earth’s surface every 5 d and has an ascending node time of 2 a.m. LT (local time) and a descending node time of 2 p.m. LT [27]. The CSES carries nine scientific instruments. These include a Langmuir probe (LAP), which can measure the plasma parameters Ne and Te, with a sweep period of 3 s in survey mode and 1.5 s in burst mode [8]. In this study, we collected and analyzed Ne data supplied as the LAP Level 02 product (https://www.leos.ac.cn, accessed on 27 June 2022). Given that the daytime plasma density is strongly influenced by solar activity [32], we used only nighttime (ascending orbit) Ne data for analysis.

2.3. CODE TEC Data

The Centre for Orbit Determination in Europe (CODE) provides GIM TEC data (http://ftp.aiub.unibe.ch/CODE/, accessed on 27 June 2022), which have been widely used by researchers around the world. These data have a spatial resolution of 2.5° in latitude and 5° in longitude; the temporal resolution is 1 h [29].

2.4. Solar-Geomagnetic Data

Many factors can induce ionospheric anomalies, which include earthquakes, geomagnetic perturbations, and solar activity. Therefore, to extract information about earthquake-related ionospheric perturbations, the influence of geomagnetic fluctuations and solar activity has to be excluded. To do this, we used the Kp, Dst, and F_10.7 indices (http://wdc.kugi.kyoto-u.ac.jp for Kp and Dst and http://www.sepc.ac.cn/ for F_10.7, both accessed on 12 July 2022). In general, space weather-related disturbances correspond to values of Kp ≥ 3, Dst ≤ −30 nT, or F_10.7 ≥ 150 sfu.

Figure 2 shows the variations in Kp, Dst, and F_10.7 from 21 February to 22 March 2022—that is, for the 30 d preceding EQ1 and EQ2. The dates with poor space weather are as follows: (a) Kp ≥ 3—21–22 February, 24–25 February, 27 February, 4–7 March, 10–15 March, 20 March, and 22 March; (b) Dst ≤ −30 nT—22 February, 5–7 March, and 11–15 March; and (c) F_10.7 ≥ 150—none. In short, the ionosphere was rarely seriously affected by solar activity but was affected by geomagnetic perturbations within the 30 d prior to EQ1 and EQ2.

3. Analysis and Results

3.1. Time-Series Analysis of CSES Ne Data

To investigate the seismo-ionospheric effects before EQ1 and EQ2, we analyzed CSES Ne data using the sliding interquartile range method. By selecting the daily orbit that fell within the study area (Figure 1A), resampling the original data to a resolution of 0.5° in the latitudinal direction, using the data from the previous 5 revisiting orbits (during 25-d observations) to calculate the background value, and computing the median value and 1.5 IQR (IQR = upper quartile–lower quartile), the relative change in Ne was calculated using Equation (1) [33]:

$$dObs=\left\{\begin{array}{c}\frac{Obs-(Median+kIQR)}{Median+kIQR}*100\%\\ 0\\ \frac{Obs-(Median-kIQR)}{Median-kIQR}*100\%\end{array}\right.\begin{array}{c}Obs>Median+kIQR\\ Obs-kIQR\le Obs\le Obs+kIQR\\ Obs<Median-kIQR\end{array}$$

(1)

where Obs is the observed Ne value, Median is the median value, k is 1.5, and IQR is the interquartile range.

The blue line in Figure 3 shows the nighttime Ne value from 21 February to 22 March 2022—that is, for the 30-d period before the occurrence of EQ1 and EQ2. The black line represents the median value, and the pink lines represent the upper and lower boundaries of the IQR. The dotted red lines represent the thresholds (±60%) used to judge whether a Ne anomaly exists. After excluding those dates that were seriously affected by space weather, the relative change in Ne was found to exceed the 60% threshold on 3, 9, 17, and 19 March—that is, 19, 13, 5, and 3 d prior to EQ1 and EQ2; the corresponding anomalies had amplitudes of 145%, 200%, 85%, and 65%, respectively. For days when the space weather was calm, these positive Ne anomalies can provisionally be considered as seismo-ionospheric anomalies related to EQ1 and EQ2.

3.2. Analysis of CODE TEC Mapping Data

To verify whether the positive anomalies in the CSES Ne data were related to EQ1 and EQ2, we used GIM TEC data. The sliding interquartile range method was also applied to these data. The median value of the ionospheric TEC for the 15 d prior to the day being considered was taken as the background value; the 1.5 IQR was also calculated. The absolute change value of TEC would be computed when it exceeded thresholds using Equation (2) [29]:

$$\Delta TEC=\left\{\begin{array}{c}TEC-(TE{C}_m+kIQR)\\ 0\\ TEC-(TE{C}_m-kIQR)\end{array}\right.\begin{array}{c}TEC>TE{C}_m+kIQR\\ TE{C}_m-kIQR\le TEC\le TE{C}_m+kIQR\\ TEC<TE{C}_m-kIQR\end{array}$$

(2)

where TEC is the observed TEC value, TEC_m is the median value, k is 1.5, and IQR is the interquartile range.

Dobrovolsky et al. [34] proposed a theoretical size for the active precursor manifestation zone, which could be calculated using Equation (3):

$$\rho ={10}^{0.43M}km$$

(3)

where $\rho$ is the radius of the earthquake preparation zone, and M is the earthquake magnitude. According to this formula, the radius of a M 6.7 earthquake (EQ1) preparation zone is 760 km.

GIM TEC global ionospheric anomaly maps are shown in Figure 4. After excluding the dates that were significantly influenced by space weather, a day was regarded as being an abnormal day if an anomaly (|∆TEC| ≥ 2 TECU, 1 TECU = 10¹⁶ el/m²) occurred and then persisted within the radius of the earthquake preparation zone for more than 2 h. It can be seen that TEC anomalies occurred from 14:00 to 16:00 UT on 2 March; from 04:00 to 06:00 and at 14:00 and 18:00 UT on 3 March; from 08:00 to 12:00 UT on 8 March; and from 08:00 to 10:00 UT on 9 March for the 30 d period before the occurrence of EQ1 and EQ2. These dates correspond to 20, 19, 14, and 13 d before EQ1 and EQ2. It is worth noting that prominent positive anomalies also occurred in the southern hemisphere conjugate zone from 10:00 to 12:00 UT on 8 March and 08:00 to 10:00 UT on 9 March, indicating that EQ1 and EQ2 may have triggered anomalies in the magnetic conjugate area.

Under calm space weather conditions, the TEC anomalies may be induced by the medium scale traveling ionospheric disturbances (MSTIDs). Based on the TEC data from GPS, many studies [35,36,37] have demonstrated that MSTIDs propagate southwestward with a propagation velocity of 100–400 m/s in the northern hemisphere. Lee et al. [38] used the TEC data from GPS to analyze the characteristics of MSTIDs in Taiwan. The results indicated that MSTIDs propagate southwestward with a propagation velocity of 100–160 m/s in Taiwan. According to the propagation velocity of MSTIDs, the shortest propagation distance of MSTIDs is 720 km within 2 h in Taiwan and the surrounding region, which is close to the radius of the earthquake preparation zone (760 km). In other words, the TEC anomalies will move out of the earthquake preparation zone after 2 h if the TEC anomalies are related to MSTIDs. The rules that we made to identify a TEC anomaly above are “After excluding the dates that were significantly influenced by space weather, a day was regarded as being an abnormal day if an anomaly (|∆TEC| ≥ 2 TECU) occurred and then persisted within the radius of the earthquake preparation zone for more than 2 h”. Thus, we conclude that the TEC anomalies in this study were not affected by MSTIDs.

Taiwan is located under the northern crest of equatorial ionization anomaly (EIA). Based on previous studies [39,40,41,42], the plasma fountain electric field should be within ±7° N magnetic. The latitude of EQ1 is about 23° N geographic, corresponding to 12° N magnetic. Thus, the distance between EQ1 and the northern boundary of the fountain electric field is about 5° in latitude which is about 555 km. Again, on the base of Dobrovolsky et al. [34], the preparation zone of EQ1 (the radius of a M 6.7 earthquake preparation zone is 760 km) could overlap with the fountain electric field region. According to the fact that the GIM TEC anomalies are located within the preparation zone of EQ1, it can be inferred that the northern crest of EIA moves poleward before the occurrence of EQ1, otherwise it will move out of the preparation zone of EQ1. Thus, we can reasonably infer that the whole process might be like this: a local electric field was produced in the earthquake preparation zone 19–20 and 13–14 d prior to the occurrence of EQ1. The local electric field led to the enhancement of the eastward fountain electric field. Driven by the enhanced E × B drift, the northern crest of EIA moved poleward, which in turn resulted in the occurrence of TEC positive anomalies in the earthquake preparation zone. Liu et al. [39,43,44,45,46] statistically analyzed the GPS TEC anomalies related to M ≥ 5.0 earthquakes in Taiwan. The statistical results demonstrated that most of these anomalies were negative and the EIA crest usually moved equatorward, which does not agree with the results of this study. The reasons behind have yet to be explored.

The results of the analysis of the GIM TEC anomalies were consistent with those for the CSES Ne anomalies: positive anomalies were found in both cases and, moreover, were concentrated on the dates 2–3 March and 8–9 March; that is, 19–20 d and 13–14 d before EQ1 and EQ2. Therefore, the analysis of the GIM TEC data verifies the analysis of the CSES Ne anomalies and indicates that the positive anomalies in the CSES Ne and GIM TEC values are related to seismo-ionospheric effects that occurred in advance of EQ1 and EQ2. Furthermore, the occurrence of CSES Ne anomalies on 17 and 19 March, corresponding to 5 and 3 d before the two earthquakes, possibly demonstrates that the CSES Ne anomalies lasted longer than the GIM TEC anomalies; alternatively, this may be due to the lack of GPS sites in the study area, which results in the GIM TEC data having a lower resolution than the CSES Ne data.

4. Statistical Analysis

To understand the general characteristics of the ionospheric anomalies associated with earthquakes in Taiwan and the surrounding region, a statistical analysis of CSES Ne anomalies related to seismic events was performed using the sliding interquartile range method. A total of 138 M ≥ 5.0 earthquakes that occurred between February 2019 and March 2022 within the range of 14.2° N–34.2° N and 111.8° E–131.8° E were investigated using the information in the CENC catalog. Figure 1A shows the locations of the earthquake epicenters. After setting a threshold of ±60%, 65 Ne anomalies that occurred during the period 2 January 2019 to 10 April 2022, were found. Most (89%) of these anomalies were positive.

To make the statistical results more convincing, the 138 earthquakes were grouped based on the following rules [47]. For the 15 d prior to or following an earthquake: (a) define the magnitude as that of the earthquake with the greatest magnitude within the radius of the earthquake preparation zone, and then group all the earthquakes within this zone as one group; (b) define the magnitude of the main shock as being the magnitude and group the earthquakes in the foreshock–main shock–aftershock sequence as one group (the difference in magnitude between the main shock and the other earthquakes must not be less than 0.7); and (c) for a double shock-type or swarm-type sequence, calculate the magnitude using Equation (4) [48] and group the earthquakes in the sequence as one group:

$$\lg E=1.5M+11.8$$

(4)

Here, E is the earthquake energy and M is the surface wave magnitude. According to the above grouping rules, the 138 earthquakes were divided into 65 groups of earthquakes.

Table 2. Relationship between the earthquake magnitude and the number and frequency of Ne anomalies (from 30 d before to 10 d after the earthquakes).

Magnitude	Number of Seismic Groups	Total Anomalies		Anomalies before Earthquakes		Anomalies after Earthquakes
		Number	Frequency	Number	Frequency	Number	Frequency
5.0–5.4	36	65	1.81	42	1.17	23	0.64
5.5–5.9	14	41	2.93	30	2.14	11	0.79
6.0–6.4	7	15	2.14	11	1.57	4	0.57
6.5–6.9	8	28	3.50	19	2.38	9	1.13
Total	65	149	2.30	102	1.57	47	0.72

lists the number and frequency of Ne anomalies corresponding to different earthquake magnitudes for the period of 30 d before to 10 d after the earthquakes belonging to each group. In general, the frequencies (the frequency of all anomalies, the frequency of anomalies before earthquakes, and the frequency of anomalies after earthquakes) gradually increase as the earthquake magnitude increases. The exception is for earthquakes in the M range of 6.0–6.4: the lack of earthquakes in this range leads to the frequency of the earthquakes in this group being lower than that of the M 5.5–5.9 earthquakes. Du [47] statistically analyzed the anomalies in CSES Ne data related to 116 M ≥ 6.0 earthquakes from around the globe and found that the frequency of the anomalies increased as the earthquake magnitude increased, which agrees well with the results of this study.

Figure 5A–D show the number and the accumulation of the Ne anomalies that occurred on different days during the period of 30 d before to 10 d after the earthquakes with magnitudes of 5.0–5.4, 5.5–5.9, 6.0–6.4, 6.5–6.9, respectively. The number of anomalies increased 8 d before the earthquakes, but dropped dramatically 7 d after the earthquakes with magnitudes of 5.5–5.9; the number of anomalies increased significantly 13 d before the earthquakes, but decreased 8 d after the earthquakes with magnitudes of 6.0–6.4; the number of anomalies increased remarkably 6 d before the earthquakes, but dropped notably 3 d after the earthquakes with magnitudes of 6.5–6.9. Overall, the Ne anomalies were concentrated during the period of 10 d before to 5 d after the earthquakes. This agrees well with the results of previous studies [8,47].

Figure 5E–H show the amplitude of the Ne anomalies that occurred on different days during the period of 30 d before to 10 d after the earthquakes with magnitudes of 5.0–5.4, 5.5–5.9, 6.0–6.4, 6.5–6.9, respectively. The largest abnormal amplitudes (≥200%) were observed 22, 11, 9, and 1 d before and 3 d after the earthquakes. In addition, the amplitudes of anomalies corresponding to the M 6.5–6.9 earthquakes are higher than those corresponding to the M 6.0–6.4 earthquakes. These results are not exactly inconsistent with the idea that the amplitude of anomalies gradually increases as the time of the earthquake occurrence approaches or as the earthquake magnitude increases.

Some of the results described above agree with the results of previous studies, whereas others do not. This illustrates that the relationship between earthquakes and ionospheric perturbations is complex. The ionospheric perturbations can be affected by many factors, such as the focal depth, tectonic background, and latitude, or by the topography close to the earthquake epicenter [8,24,27]. In turn, this can affect the number and amplitude of the detected anomalies.

5. Discussion

At present, there are three main models that describe the seismo-ionospheric coupling mechanism [49]. The first model, the lithosphere–atmosphere–ionosphere multichannel coupling model, was proposed by Hayakawa [50] in 2004. This model was further summarized by Ding et al. [51] in 2006. The second model, the unified lithosphere–atmosphere–ionosphere coupling model, was proposed by Pulinets and Boyarchuk [9] in 2004. According to this model, precursor anomalies are caused by radon, inert gas, and greenhouse gas emissions, as well as air turbulence. The advantage of this model is that it explains the mechanisms behind the generation and propagation of abnormal LF, high frequency (HF), and very high frequency (VHF) radio waves, which the first model does not. In 2011, using the further developed model, which they named the unified sphere-coupling model, Pulinets and Ouzounov [52] explained the process by which the local electric field anomaly is formed—possibly through the decay of radon and its progeny through a series of physical and chemical reactions during the earthquake preparation. The third model, the lithosphere–coversphere–atmosphere coupling model, was proposed by Wu et al. [53] in 2012. This model is mainly used to explain the mechanism behind the formation of thermal infrared anomalies. The first model involves the interactions between three spheres and has been widely used. Therefore, in this study, we based our considerations of the seismo-ionospheric coupling mechanism on the first model and, to some extent, on the second model.

The LAIC model includes three channels: the geochemical channel, the acoustic channel, and the electromagnetic channel. Based on the collected data related to earthquake precursor anomalies in Taiwan and the surrounding region, we believe that all three of these channels may constitute pathways by which pre-earthquake ionospheric anomalies form in the study area. The evidence for this includes the following. (1) Geochemical channel: Liao et al. [54] statistically analyzed underground fluid data from Fujian Province acquired before nine M ≥ 7.0 Taiwan earthquakes that occurred during the period of 1990–2004 and detected water radon, F⁻, Cl⁻, water-soluble SiO₃²⁻, and water conductivity anomalies before the earthquakes. King and Chia [55] suggested anomalous streamflow and groundwater-level changes before the 1999 Chi-chi earthquake in Taiwan. Groundwater-levels recorded at a well 1.5 km east of the seismogenic fault showed an anomalous rise 2 d before the earthquake, and then a unique 4 cm drop beginning 3 h before the earthquake. In other words, the radon anomalies and changes in water elevation that occurred in the geochemical channel before the earthquakes may have induced ionospheric perturbations. (2) Acoustic channel: In most cases, brightness temperature increases with the growing terrestrial heat flow. Genzano et al. [56] monitored thermal infrared anomalies that occurred before M > 4.0 earthquakes from 1995 to 2002 using the Robust Satellite Techniques (RST). Taking the 1999 Chi-chi earthquake as an example, a thermal infrared anomaly appeared two weeks in advance of the earthquake, close to the epicenter and along the associated tectonic lineaments. Based on Moderate Resolution Imaging Spectroradiometer (MODIS) thermal infrared satellite remote sensing data, Yuan et al. [57] reported that a thermal anomaly near the epicenter reached its maximum size and amplitude 6 d before the 2016 M 6.2 Taiwan earthquake. Furthermore, Lin [58] studied the precursor anomalies characteristics of underground fluid data from Fujian Province that preceded the 2003 M 7.0 Taitung Sea earthquake. The results indicated that the water temperature anomalies with increasing features at the Wuyazui well occurred 180 d before the earthquake. It was concluded that abnormal heat sources in the acoustic channel before the earthquake may have caused ionospheric perturbations. (3) Electromagnetic channel: Qian et al. [59] observed ultra-low frequency (ULF) electromagnetic anomalies that preceded the 1999 Chi-chi earthquake. After performing simulated experiments, the authors concluded that the strong stress that occurred before the earthquake in the source area may have led to the development of rock microcracks and that the ULF signals were then produced by piezoelectric or piezomagnetic effects. It is clear that the abnormal ULF signals in the electromagnetic channel before the earthquake may have induced ionospheric perturbations.

By combining tectonic background information with the results of previous studies [9,49,50,51,52,54,55,56,57,58,59], we developed a LAIC model of Taiwan and the surrounding region. A three-dimensional (3D) representation of this is shown in Figure 6. Earthquakes frequently occur in Taiwan and the surrounding region due to subduction and the convergence between the Eurasian plate and the Philippine Sea plate. Before an earthquake occurs, micro-fractures, liquid diffusion, and pressure variations are concentrated near the earthquake source. A series of geochemical interactions near the epicenter, including radon decay, and the changes in water elevation lead to changes in the composition of the air and its conductivity. This, in turn, produces a local electric field anomaly. Finally, this field influences the plasma density in the ionosphere. The anomalous heat sources open the acoustic channel, causing acoustic gravity wave (AGW) turbulence, and enhancing gravity and planetary waves; this then leads to synchronous oscillations in the ionospheric plasma density and temperature when these waves are transmitted to the ionosphere. The electromagnetic channel is opened by the anomalous ULF signals. These signals penetrate the ionosphere and interact with energetic particles in the magnetosphere, causing the particles to precipitate into the ionosphere and produce density modulation.

In order to study the characteristics of the LAIC mechanism associated with EQ1 and EQ2, we collected and analyzed underground fluid data from the Fujian Province near Taiwan at southeast China within the radius of the earthquake preparation zone (760 km). Because of the high quality and prediction efficiency of the data, F⁻ data from the Xiamen Dongfu well, gas radon data from the Ningde well, and water radon, F⁻, Cl⁻, HCO₃⁻, and water conductivity data from the Huaan Tainei well have been used as mid-short term earthquake prediction indices [60]. In this study, we focused on whether precursor anomalies existed in these underground fluid indices during the three-month period before EQ1 and EQ2 (23 December 2021–22 March 2022). Three such anomalies were found;

Table 3. Details of underground fluid precursor anomalies before EQ1 and EQ2 (23 December 2021–22 March 2022).

Basic Information				Description of Anomaly		Relationship between the Anomaly and EQ1 and EQ2
Well	Index	Lat./°N	Lon./°E	Start and End Time (BJT)	Character Description	Number of Days between the Occurrence of the Anomaly and EQ1 and EQ2/d	Distance/km
Xiamen Dongfu	F−	24.5	117.9	5 March 2022	>1.6 σ (standard deviation)	18	387
Ningde	Gas radon	26.7	119.5	31 January–24 February 2022	The gas radon concentration occurred, but the temperature did not change.	51	411
Huaan Tainei	Water radon	24.7	117.6	11 February 2022	The water radon value was higher than 90 Bq/L and oscillated strongly. The temperature was normal.	40	430

lists detailed information about these. The F⁻ anomaly at the Xiamen Dongfu well, the gas radon anomaly at the Ningde well, and the water radon anomaly at the Huaan Tainei well appeared 18, 51, and 40 d before the earthquakes, respectively. On this basis, we suggest that it is the geochemical channel that provides the pathway for ionospheric perturbations. Due to the lack of relevant data, it cannot be said whether the acoustic channel and electromagnetic channel were open or not during these two cases of EQ1 and EQ2.

6. Conclusions

The ionospheric anomalies associated with two successive earthquakes that occurred in 2022 in Taitung Sea, Taiwan, China, with magnitudes (M) of 6.7 and 6.3, at 23.45° N, 121.55° E and 23.39° N, 121.52° E and with the same focal depth of 20 km, were investigated in detail in this study. Under calm space weather conditions, the CSES Ne data and CODE TEC data verified each other and both exhibited anomalies on 2–3 March and 8–9 March—that is, 19–20 and 13–14 d before the earthquakes. In addition, the GIM TEC data suggested that positive anomalies may have been triggered in the magnetic conjugate area 13–14 d prior to the earthquakes, and the CSES Ne data exhibited positive anomalies 3 and 5 d before the earthquakes.

Next, a statistical analysis of 138 M ≥ 5.0 earthquakes that occurred in Taiwan and the surrounding region during the period February 2019 to March 2022 showed that most of the associated anomalies in the CSES Ne data were positive. It was found that, as the earthquake magnitude increased, the frequency of the anomalies increased, but their amplitudes did not. The statistical results also demonstrated that the anomalies were concentrated during the period of 10 d before to 5 d after the earthquakes. Moreover, as the dates of the earthquakes approached, the amplitude of the anomalies did not increase, whereas we had predicted that they would gradually increase. In fact, the anomalies with the largest amplitude occurred 22, 11, 9, and 1 d before and 3 d after the earthquakes.

Finally, based on evidence related to anomalies in radon, water elevation, heat sources, and ULF signals, we suggest that earthquakes in Taiwan and the surrounding region affect the ionosphere by means of the geochemical, acoustic, and electromagnetic channels that are described by the LAIC model. In particular, based on the earthquake precursor anomalies that were found in gas radon, water radon, and F⁻ data, it is possible that EQ1 and EQ2 influenced the ionosphere through the geochemical channel. At present, there is no evidence linking EQ1 and EQ2 with the observed anomalies through the acoustic and electromagnetic channels.