Multi-Dimension and Multi-Channel Seismic-Ionospheric Coupling: Case Study of Mw 8.8 Concepcion Quake on 27 2010

: GPS radio occultation (RO) technology can fully describe the subtle structure of the ionosphere. This paper discusses the dynamic abnormity observed by the RO data from the Constellation Observing System for Meteorology Ionosphere and Climate (FORMOSAT-3/COSMIC) before the great earthquake case in Concepcion, Chile (27 February 2010, Mw 8.8). Traditional ground-based GPS monitoring was considered as the external conditions and references to the excitation response. Using kriging interpolation, the global Nmf2 map (GNM) was ﬁrst constructed to study the ionosphere deviation from the normal state. Successively, the ionosphere abnormality in the F2 region (Nmf2), vertical structure (RO proﬁles), and multiple heights (electron density) of traveling are unfolded. The Nmf2 disturbances in the possibility of seismic inﬂuences were excluded from non-seismic noise factors, including the external input (e.g., space weather activity, 15 February) and meteorological events (e.g., lower atmospheric forcing in quiet periods). However, the results show that there were apparent local Nmf2 perturbations for up to 5 h in the epicenter area on 21 and 25 February. The disturbances of the RO proﬁles and the interaction of other layers of the ionosphere implied the ﬂuctuation signals of prominent long-wavelength ﬂuctuations >50 km in the F layer. The ionospheric ﬂuctuates wildly, and these wave signals considered as the trace of gravity wave propagating upward are mainly distributed at the elevation of 200–300 km. The simultaneous reaction of GNSS TEC further evidenced the potential possibility of acoustic gravity by the COSMIC RO proﬁles, reﬂecting the compounding couplings of seismo-ionosphere effects. In terms of the presentation of VLF radiation noise and the aerosol ion clusters, the electromagnetic and chemical channels have been previously completed by DEMETER and Terra/Aqua satellites. These ﬁndings implied the great potential of the FORMOSAT-7/COSMIC-2 system (now in the testing phase), with ~5000 soundings to investigate the subtle atmospheric stratiﬁcation.


Introduction
When an earthquake occurs, the earth's crust releases a large amount of energy, which acts on the earth's surface and affects the ionosphere [1]. The progress in statistical analyses suggests an undisputed connection between the quakes and the existence of ionospheric anomalies in the ionosphere. However, there is still an uncertainty in earthquake anomaly reports due to a lack of unified, indisputable, and explicit definitions, how they are generated, and why they spread before and after the quake over the seismogenic zones. the number of occultation profiles has decreased below 1000. The US Air Force Space Test Program successfully launched six FORMOSAT-7/COSMIC-2 satellites into a 24 deg inclination earth orbit on 25 June 2019. Because of the recent completion of data calibration and validation, the COSMIC generation-2 system now has massive data-more than the initial generation of products. The enhanced performance has potential for the ionospheric coupling mechanism and the detection of abnormal effects [20].
At 6:34 GMT, 27 February 2010, a shallow devastating quake with a magnitude of 8.8 on the Richter scale hit near Concepcion, Chile. The earthquake's focal depth was only 55 km, and the epicenter was in the waters near Maule, 320 km southwest of Santiago, Chile. This event is a typical thrust plate event. Nazca Plate, located in the west of the earthquake area, subducted rapidly to the South American plate at a compression rate of 65 mm/yr. The action resulted in the gradual accumulation of interplate stress until surface fracture [21]. After the earthquakes, the city moved 3 m westward. In 2010, the COSMIC system could stably provide~1500 detection profiles every day. Rich profile data provides a possibility for the vertical structure detection of the seismo-ionosphere effect [22].
Employing a case study of the Mw 8.8 Concepcion earthquake on 27 February 2010, the COSMIC occultation observation system's ability to monitor the ionosphere disturbance was studied. Based on statistical analysis, the background variation's effectiveness is verified, and a strict anomaly condition is given. After that, two aspects of the ionospheric anomalies in the F2 layer and other heights were investigated successively. On the one hand, by comparing the observed value with the constructed global Nmf2 map (GNM), the abnormal Nmf2 point fields over the detected area (10 • E-130 • E, 10 • N-60 • S) were detected. On the other hand, by observing the electron density profiles with a large deviation, the ionosphere fluctuation and anomalies in different heights were analyzed. This research, we believe, has practical and referential significance for the measurement model of the ionosphere and climate variations when FORMOSAT-7/COSMIC-2 is fully operational in the future.

Solar and Geomagnetic Data
The ionosphere is easily disturbed by solar and geomagnetic activities. Solar and geomagnetic activities should be analyzed prior to ionosphere anomaly detection before earthquakes [23]. The solar activity is reflected by the solar flares and radio flux F10.7. The geomagnetic activity is depicted by the disturbance storm time (Dst), auroral electrojet (AE), and Kp indices to present the equatorial to mid-high latitudes geomagnetic activity. The F10.7 and solar flares data are provided by the space center of the Chinese Academy of Sciences (http://www.sepc.ac.cn, accessed on 1 June 2021), with a time resolution of 1 d. A quiet solar state is less than 100 SFU and an M level for F10.7 and solar flares, respectively [24]. The Dst, AE, and Kp index are provided by the Kyoto geomagnetic data center, Japan (http://wdc.kugi.kyoto-u.ac.jp/dstae/index.html, accessed on 1 June 2021), and the Goddard Space Flight Center (https://omniweb.gsfc.nasa.gov/form/dx1.html, accessed on 1 June 2021), respectively. The time resolution of these geomagnetic indices is 1 h, with −30-30 nT for Dst and 0-4 for Kp under the stable condition of terrestrial magnetism [6,7,25].

Ionosphere Data
Electron density profiles, provided by the FORMOSAT-3/COSMIC, were used to investigate the F2 peak electron density (Nmf2), vertical structure oscillation (RO profiles), and ionosphere variations in multiple heights before the Concepcion shock. Over 1500 ionospheric occultation profiles available were acquired daily (https://www.cosmic. ucar.edu, accessed on 1 June 2021). In previous studies, comparisons of COSMIC radio occultation with ground-based ionospheric observations have been made. The verification results showed that the electron density profiles were consistent with the measurements from ionosondes and incoherent scatter radars (ISR) [26,27].
In terms of ground-based systems, the TEC from signal station receivers and global ionospheric maps (GIMs) are commonly used ionospheric disturbance analysis tools. The Remote Sens. 2021, 13, 2724 4 of 22 former is dominated by GNSS or ionosonde networks (e.g., CORS observations and the madrigal database), covering the continental regions of the earth. By comparison, the GIMs from International GNSS Service (IGS) were obtained by weighing the products for each Ionosphere Analysis Centers with different solution strategies. The GIM model takes data from more than 400 stations and then interpolates and smooths them to produce ionospheric TEC maps on one-or two-hour increments (ftp://cddis.gsfc.nasa.gov/gnss/ products/ionex/, accessed on 1 June 2021). Covering ±180 degrees of longitude and ±87.5 degrees of latitude with spatial resolutions of 5 and 2.5 degrees, respectively, its performance in marine areas can also be guarded [28]. Due to the Concepcion quake located in the coastal area, both the COSMIC RO ionospheric responses over the sea and land all need to be verified. Thus, the GIM TEC was further used to make a feasibility analysis of multi-dimensional and -channel seismo-ionosphere coupling by COSMIC RO profiles.

Data Preprocessing and Background Field Construction
The original profile database involved no quality control. Due to meteorological, instrument quality, and other factors, some irregular RO profiles were inevitable and controlled. First, the unqualified data with strong fluctuation were removed; for the detailed preprocessing, see Appendix A. For the final-presented analysis, around 20,000 observations were classified for the next procedure of constructing the Nmf2 model.
The common calculation method of the ionosphere's spatial distribution reference value includes a COSMIC empirical model value but usually needs to use long-term data [29]. Another method is to establish the global or regional model of electronic density quickly before the earthquake, after preprocessing RO profiles from COSMIC occultations. A kriging interpolation method was used here to obtain the fast Global Nmf2 maps. This technique can yield valuable information about the ionosphere's behavior and its evolution, with particular attention paid to the evolution of spatial gradients [30,31]. The ionosphere parameter (e.g., TEC and Ne) had been previously kriging interpolated and drawn mainly with the combination of ground-based techniques. Stanislawska et al. [32] firstly introduced the data of only 14 GNSS stations to elaborate TEC maps over Europe with a time rate of 1 h. Furthermore, Minkwitz et al. [33] proved the possibility of constructing a 3-D electron density region by a simple kriging approach. Regardless of the reconstructed TEC or electron densities, promising improvements have been achieved through the validation of the ionosondes Ne and GNSS TEC [34]. The method can not only solve the problem of overly sparse sampling, especially in situations of data scarcity but also give the estimation deviation of each point. To detect the seismic-ionosphere disturbed signals, we constructed the Nmf2 models in a short period before an earthquake by the kriging interpolation algorithm. For more specific operations, please see Appendix B.
Following the usual representation of the global ionospheric map such as IGS GIMs and IRI GNMs, these file products are provided with the Geographic Frame Coordinate System. However, inevitably, changes in the TEC or Nmf2 maps are influenced by the factors of the earth's rotation, the earth's pole moving, etc. If their strategy is adopted, the data sparsity of the ionosphere retrieved from the COSMIC constellation will greatly affect the result's reliability (see Supplementary Figure S1). In this paper, we introduce an Inertial Coordinate System of the Geomagnetic Sun-fixed Frame Coordinate System (local time and geomagnetic latitude). In the sun-fixed coordinate system, where the impact factors of the earth's rotation and the earth's pole moving are ruled out, the diurnal variation of the ionosphere tends to be stable (see Supplementary Figure S2). Therefore, the data of the ionosphere NmF2 are easily modeled [12]. We first extracted the COSMIC NmF2 data within 15 days before the global earthquake and divided the data in a day into 12 periods, with 2 h as the time interval. The data were then classified in the same period into the geomagnetic sun-fixed coordinate system, forming a total of 12 ionosphere GNMs based on kriging interpolation.
After constructing the Nmf2 models, the observation value was taken as the inspection value to compare its variations with the ionosphere's normal state. The deviation degree between NmF2 observations and model values is expressed as follows: where NmF2 obs and NmF2 model represent the observation value and GNM value, respectively. According to the deviation degree with a large δ value, some irregular electron density profiles may exhibit a damaged structure. Such profiles may also occur inside ionospheric storms or be caused by medium/small-scale irregularities [34]. The specific flow chart is shown in Figure 1. After data preprocessing, the global NmF2 models were constructed as background references, and their validity was tested by comparing global ionosphere maps (GIMs). The local F-layer anomaly in seismic effects was explored with corresponding research on RO profile fluctuations. Later, the ionosphere anomalies in different elevation and the simultaneous reaction of GIM TECs was examined successively.

Feasibility Study on GNM Nmf2
After processing the electronic density profile of COSMIC based on the kriging interpolation algorithm, the GNM models and the fitting error (root mean square error) were obtained. Figure 2 shows the Nmf2 spatial variation for a period of 0-2 h. In accordance with the spatial distribution of TEC, there was a high electron density on both sides of the

Feasibility Study on GNM Nmf2
After processing the electronic density profile of COSMIC based on the kriging interpolation algorithm, the GNM models and the fitting error (root mean square error) were obtained. Figure 2 shows the Nmf2 spatial variation for a period of 0-2 h. In accordance with the spatial distribution of TEC, there was a high electron density on both sides of the geomagnetic equator in the range of 60 • W-120 • E, 30 • N-30 • S. This "double hump" structure was influenced by the effect of the magnetic equatorial fountain [35].
Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 24 geomagnetic equator in the range of 60° W-120° E, 30° N-30° S. This "double hump" structure was influenced by the effect of the magnetic equatorial fountain [35]. To prove the reliability of Nmf2 construction, the applicability Nmf2 model in different latitudes was verified using the GIM data and the Nmf2 observations. To quantitatively assess the Abel retrieval and kriging interpolation, we set the cell of the GNM with the same size as the GIM (Global Ionosphere Map). All the cells were compared with the corresponding GIM TEC, as shown in Figure 3. The main panel represents the case when the Nmf2 and TEC scatter on the global level (87.5° N-87.5° S, 180° W-180° E) is statistical. The attached top/bottom panel shows the case when both cell values in the mid-low latitude (60° N-60° S) and the high latitude (90°-60° N or 90°-60° S) are accounted for, respectively. The corresponding correlation coefficients (R) are presented in respective panels. The GNMs (attached top panel) in the mid-low latitude showed good modeling accuracy, with the correlation coefficient increased from 0.9165 (global in the main panel) to 0.9654. In the bottom panel, their correlation became weaker with R = 0.6779. This is because the zones with weak correlations were all adjacent to high-latitude polar regions. Note that, along the descent latitude, the spatial sparsity and error amplification of Nmf2 observations led to Nmf2 models in high-latitude cells (especially in the polar zone) equal to zero or even less than zero, and this is nonexistent in reality. The results presented that the GNM Nmf2 in the mid-low latitude with a superior correlation for GIM TEC had a better modeling accuracy than the other zones. Naturally, this is suitable for pre-earthquake ionosphere analysis.  To prove the reliability of Nmf2 construction, the applicability Nmf2 model in different latitudes was verified using the GIM data and the Nmf2 observations. To quantitatively assess the Abel retrieval and kriging interpolation, we set the cell of the GNM with the same size as the GIM (Global Ionosphere Map). All the cells were compared with the corresponding GIM TEC, as shown in to 0.9654. In the bottom panel, their correlation became weaker with R = 0.6779. This is because the zones with weak correlations were all adjacent to high-latitude polar regions. Note that, along the descent latitude, the spatial sparsity and error amplification of Nmf2 observations led to Nmf2 models in high-latitude cells (especially in the polar zone) equal to zero or even less than zero, and this is nonexistent in reality. The results presented that the GNM Nmf2 in the mid-low latitude with a superior correlation for GIM TEC had a better modeling accuracy than the other zones. Naturally, this is suitable for pre-earthquake ionosphere analysis. The initial GNM model was output under the inertial sun-fixed coordinate system (local time and geomagnetic latitude), and the above operations for analyzing external conditions were completed in this system. In the following pre-earthquake ionosphere analysis, the background variations are analyzed based on an earth-fixed assumption to capture the daily abnormal signals of daily periods. cause the zones with weak correlations were all adjacent to high-latitude polar regions. Note that, along the descent latitude, the spatial sparsity and error amplification of Nmf2 observations led to Nmf2 models in high-latitude cells (especially in the polar zone) equal to zero or even less than zero, and this is nonexistent in reality. The results presented that the GNM Nmf2 in the mid-low latitude with a superior correlation for GIM TEC had a better modeling accuracy than the other zones. Naturally, this is suitable for pre-earthquake ionosphere analysis.

Space Environment Analysis
The active solar and geomagnetic state will cause a great disturbance to the ionosphere [36,37]. The variations of these two factors should be analyzed before analyzing ionosphere disturbance. Figure 4 shows the changes in the Dst, Auroral electrojet (AE), and Kp indices 15 days before the Concepcion earthquake. It can be seen that the geomagnetic activity is relatively active 12 days before the earthquake (February 15) and 11 days before the earthquake (February 16), and the Dst index drops below −30 nT. Simultaneously, the Kp index was 4, with the AE index (expansion phase) dramatically expanding to 720 nT. In these two periods, 12 h on February 15 and 6 h on February 16 reach the level of medium and small geomagnetic storms. The initial GNM model was output under the inertial sun-fixed coordinate system (local time and geomagnetic latitude), and the above operations for analyzing external conditions were completed in this system. In the following pre-earthquake ionosphere analysis, the background variations are analyzed based on an earth-fixed assumption to capture the daily abnormal signals of daily periods.

Space Environment Analysis
The active solar and geomagnetic state will cause a great disturbance to the ionosphere [36,37]. The variations of these two factors should be analyzed before analyzing ionosphere disturbance. Figure 4 shows the changes in the Dst, Auroral electrojet (AE), and Kp indices 15 days before the Concepcion earthquake. It can be seen that the geomagnetic activity is relatively active 12 days before the earthquake (February 15) and 11 days before the earthquake (February 16), and the Dst index drops below −30 nT. Simultaneously, the Kp index was 4, with the AE index (expansion phase) dramatically expanding to 720 nT. In these two periods, 12 h on February 15 and 6 h on February 16 reach the level of medium and small geomagnetic storms.
According to the solar flare data, the M flares occurred 21-19 and 15 days before the Concepcion earthquake (February 6-8 and February 12) but without the proton event. Besides these periods, the solar activity before the shocks were relatively quiet (F10.7 < 100 SFU), with weak force fluctuations on the ionosphere.

Ionospheric Structure Anomaly Before Earthquake
The changes in the observed values of Nmf2 before the earthquake relative to the modeled values in the detection area (60° S-20° N, 90° W-50° W) were first explored. Due to the complex and dynamic interactions, the δ value varied when time and region differed. After constructing GNM models, the statistics of δ values from Nmf2 observations are necessary. As shown in Figure 5, the δ values in the detection zone covering the Con- According to the solar flare data, the M flares occurred 21-19 and 15 days before the Concepcion earthquake (February 6-8 and February 12) but without the proton event. Besides these periods, the solar activity before the shocks were relatively quiet (F10.7 < 100 SFU), with weak force fluctuations on the ionosphere.

Ionospheric Structure Anomaly before Earthquake
The changes in the observed values of Nmf2 before the earthquake relative to the modeled values in the detection area (60 • S-20 • N, 90 • W-50 • W) were first explored. Due to the complex and dynamic interactions, the δ value varied when time and region differed. After constructing GNM models, the statistics of δ values from Nmf2 observations are necessary. As shown in Figure 5, the δ values in the detection zone covering the Concepcion epicenter were extracted from 13 to 26 February. After data preprocessing of the COSMIC profiles, the filtered RO events were collected to compute the δ values, and their two-fold average value (0.3433) was made as to the Nmf2 anomaly detection criterion. Here, the threshold results of δ threshold = 0.3433 were used as an ionosphere boundary limitation for systematic error and abnormal disturbance. If the δ value was lower than the δ threshold , the ionosphere was regarded as in the normal state; otherwise, it was considered to be affected by large energy fluctuations. It can be seen that, while most derivation was within bounds, Nmf2 anomalies could be detected at certain times, even those far exceeding the threshold. These responses indicated strong interactions with external factors.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 24 two-fold average value (0.3433) was made as to the Nmf2 anomaly detection criterion.
Here, the threshold results of δthreshold = 0.3433 were used as an ionosphere boundary limitation for systematic error and abnormal disturbance. If the δ value was lower than the δthreshold, the ionosphere was regarded as in the normal state; otherwise, it was considered to be affected by large energy fluctuations. It can be seen that, while most derivation was within bounds, Nmf2 anomalies could be detected at certain times, even those far exceeding the threshold. These responses indicated strong interactions with external factors. Four daily time phases (0-6 UT, 6-12 UT, 12-18 UT, and 18-24 UT) in this section are defined. Within the detection zone (60° S-20° N, 90° W-50° W), the daily latitude-time variations impending the earthquake for determining anomalous Nmf2 are analyzed(see Figure 6). The δ value of Nmf2 disturbance detected at most observation points is less than 0.2, while excess outliers appear in three phases: 6-18 h on February 15, 12-18 h on February 21, and 6-12 h on February 25. To separate the influence of internal/external factors from the seismo-ionosphere effects, the uniqueness of regional enhancement near the epicenter needs to be further explained. Responding to RO Nmf2 δ > 0.3433, all the disturbing events within 15 days prior to the shock were extracted. Three episodes were delineated: the magnetic-influence episode (February 15), the local effect episode (February 21 and 25), and the other quiet episodes. Based on the spatial distribution in RO in Figure 7, strongly disturbing events and local effects were studied in each episode.   Figure 6). The δ value of Nmf2 disturbance detected at most observation points is less than 0.2, while excess outliers appear in three phases: 6-18 h on February 15, 12-18 h on February 21, and 6-12 h on February 25. To separate the influence of internal/external factors from the seismo-ionosphere effects, the uniqueness of regional enhancement near the epicenter needs to be further explained. Responding to RO Nmf2 δ > 0.3433, all the disturbing events within 15 days prior to the shock were extracted. Three episodes were delineated: the magnetic-influence episode (February 15), the local effect episode (February 21 and 25), and the other quiet episodes. Based on the spatial distribution in RO in Figure 7, strongly disturbing events and local effects were studied in each episode. In recent years, among the Mw 8+ earthquakes in the South American plate, the ionosphere anomaly usually occurs within 8 days before the earthquake [7]. There is an obvious correlation between different ionosphere parameters [39]. For instance, several ionosphere parameters (GPS-TEC, ion+, electronic temperature K, etc.,) anomalies in the  Figure 7A shows the global ionospheric disturbances on February 15. The results show that a wide range of ionospheric anomalies occurred worldwide. With the active geomagnetic activity and occurrence well before the shock (12 days prior to the shock), the geomagnetic activity is the main factor of the Nmf2 anomalies. By contrast, in Figure 7B, with the sporadic points outside the seismogenic area, there was a noticeable disturbance in the F2 layer near the epicenter on February 21 and February 25. A pronounced local effect on both days is clearly visible.
Lastly, the F2 layer reaction in the quiet episode was inspected (see Figure 7C). There was a significant ionosphere variation in the South American sector and North America high latitude within the 15 days before the shock. In the Atlantic sector (near 0 • E), the disturbed Nmf2 could be captured as well. These Nmf2 anomalies have a rather similar spatial distribution of ionosphere scintillation [38]. Therefore, the 15-day investigation showed a unique regional enhancement on February 21 and February 25 near the Concepcion epicenter. The large-scale Nmf2 abnormal and sporadic ionosphere scintillation in other periods, however, had no such mesoscale local effect as an important distinction of "source" type perturbations from other anomalies. The next section discusses the disturbance of electron density in the F2 layer and other ionosphere heights near the earthquake area on February 21 and February 25.
In recent years, among the Mw 8+ earthquakes in the South American plate, the ionosphere anomaly usually occurs within 8 days before the earthquake [7]. There is an obvious correlation between different ionosphere parameters [39]. For instance, several ionosphere parameters (GPS-TEC, ion+, electronic temperature K, etc.,) anomalies in the south of the epicenter six days before the events were detected. Near the same geological fault, the ionosphere anomaly drifts from the south-central area to the epicenter north before the Chilean earthquake (April 1, 2014, Mw 8.2) have been previously recorded [7].
At 12-18 UT, February 21, when seismo-ionosphere effects occurred, both disturbed and quiet intervals were divided. Figure 8 shows that nine COSMIC occultation events occurred in the disturbance interval (14)(15)(16)(17)(18). Except for the northeast δ of <δ threshold , the occultation points in other directions all had obvious disturbance characteristics. In this interval, the strongest Nmf2 anomaly area of the F2 layer appeared in the southeast of the epicenter (δ > 0.5). Unfortunately, in the quiet interval (12-14 UT), only five COSMIC occultation events appeared in the south of the epicenter. Our deep exploration of the northern ionosphere disturbance characteristics was hindered. The importance of the spatial distribution and sparsity of the COSMIC occultation events is self-evident. During the disturbance interval, the ionosphere's abnormal state in other directions as well as the central and southern parts of the epicenter needs to be further confirmed. Figure 9 shows 10 COSMIC occultation events near the epicenter from 06:00 to 12:00 on February 25. According to the same definitions in Figure 8, six and four COSMIC occultation observation points were captured, respectively, in the disturbance (07-12 UT) and the quiet interval (06-07 UT). During the disturbance interval, the COSMIC occultation event distribution presented large-scale atmospheric oscillations in the F2 layer. A sharp energy exchange of lithosphere-atmosphere may have occurred for up to 5 h over the ionosphere. With the intense disturbance (δ > 70%) at 10:32 UT and 11:54 UT, the ionosphere's reaction reached its peak when it diffused around the atmosphere.        By COSMIC RO profiles, the change in electron density before some shocks (e.g., the Tōhoku event [40]) could be recorded in a response to violent fluctuations. Because the ionosphere may still vibrate even during quiet conditions, a reasonable strategy was adopted to estimate the limitation of "normal" fluctuations. Figure 10a shows that a slight fluctuation in the ionospheric profiles was a common phenomenon (−2.75-2.93 × 10 3 el·cm −3 ). Most of them were within a certain range, and other outliers were in an irregular fluctuation.
Using the singular spectrum analysis filter (see [5] for details), the positive/negative highfrequency signals of the profiles remained. We calculated the double averaging value of the positive high-frequency data as the upper limitations. With the same definition, the double averaging value of the negative high-frequency signals was regarded as the lower limit. Liu et al. [40] revealed that the acoustic gravity wave excited by the Rayleigh wave can cause prominent long wave fluctuations of 50-105 km in the vertical electron density profile. In panels (b) and (c), three profiles, #1~#3 (the strongest Nmf2 disturbance on February 21, see Figure 8), were analyzed with a similar feature. The periodic or quasiperiodic oscillation was quite prominent in Profiles #1 and #3, with long-wave atmospheric oscillations, even to a height >400 km. On the contrary, the fluctuation in Profile #2 failed to be regarded as an abnormal state. These fluctuations leave traces that are triggered by the bottom atmosphere, propagate upward to the F layer, and ultimately disappear in the upper ionosphere. These features imply potential evidence of acoustic transmission in the seismogenic area. Though restricted to sparse occultation events, the acoustic gravity wave's specific propagation mode was hard to determine. We are optimistic that, with dense RO events in the following 2-generation COSMIC, an enhanced performance is conducive to collecting more interesting signals.
By COSMIC RO profiles, the change in electron density before some shocks (e.g., the Tōhoku event [40]) could be recorded in a response to violent fluctuations. Because the ionosphere may still vibrate even during quiet conditions, a reasonable strategy was adopted to estimate the limitation of "normal" fluctuations. Figure 10a shows that a slight fluctuation in the ionospheric profiles was a common phenomenon (−2.75-2.93 × 10 3 el·cm −3 ). Most of them were within a certain range, and other outliers were in an irregular fluctuation. Using the singular spectrum analysis filter (see [5] for details), the positive/negative high-frequency signals of the profiles remained. We calculated the double averaging value of the positive high-frequency data as the upper limitations. With the same definition, the double averaging value of the negative high-frequency signals was regarded as the lower limit. Liu et al. [40] revealed that the acoustic gravity wave excited by the Rayleigh wave can cause prominent long wave fluctuations of 50-105 km in the vertical electron density profile. In panels (b) and (c), three profiles, #1~#3 (the strongest Nmf2 disturbance on February 21, see Figure 8), were analyzed with a similar feature. The periodic or quasi-periodic oscillation was quite prominent in Profiles #1 and #3, with longwave atmospheric oscillations, even to a height >400 km. On the contrary, the fluctuation in Profile #2 failed to be regarded as an abnormal state. These fluctuations leave traces that are triggered by the bottom atmosphere, propagate upward to the F layer, and ultimately disappear in the upper ionosphere. These features imply potential evidence of acoustic transmission in the seismogenic area. Though restricted to sparse occultation events, the acoustic gravity wave's specific propagation mode was hard to determine. We are optimistic that, with dense RO events in the following 2-generation COSMIC, an enhanced performance is conducive to collecting more interesting signals. Two groups of COSMIC occultation events by the ionosphere maximum (δ > 0.7, Prof#4/5/6) and minimum disturbances (δ < 0.1, Prof#7/8/9) in Figure 9 are divided. The Two groups of COSMIC occultation events by the ionosphere maximum (δ > 0.7, Prof#4/5/6) and minimum disturbances (δ < 0.1, Prof#7/8/9) in Figure 9 are divided. The corresponding electron density profiles were extracted to judge the vertical structure change in the ionosphere on the second day before the case. During the nighttime (local time), the "normal" ionosphere structure information was captured by the Prof#7-9 (around 6 a.m.) COSMIC event group. A typical F layer and E layer were visible in this group. Yet such a familiar ionosphere structure was collapsed in the matinal time. At the Prof#4-6 (collected at 10-12 a.m.) group, the vague stratification of the F layer (F1/F2 layer) emerged.
In the right panel of Figure 11, under the control of day-night transmission and the particle precipitation, the F2 layer peak electron density at Prof#4-6 recorded the largest value. However, in the adjacent space and time, the Nmf2 at Prof#7-9 is significantly lower than the former group, consistent with the ionospheric abnormal enhancement at two points. In [41], it was found that the ionosphere peak layer height (hmF2) enhanced abnormally before the earthquake. Based on the hmF2 the of both event groups in the right panel, the extraordinary feature before the Concepcion earthquake was not prominent. corresponding electron density profiles were extracted to judge the vertical structure change in the ionosphere on the second day before the case. During the nighttime (local time), the "normal" ionosphere structure information was captured by the Prof#7-9 (around 6 a.m.) COSMIC event group. A typical F layer and E layer were visible in this group. Yet such a familiar ionosphere structure was collapsed in the matinal time. At the Prof#4-6 (collected at 10-12 a.m.) group, the vague stratification of the F layer (F1/F2 layer) emerged. In the right panel of Figure 11, under the control of day-night transmission and the particle precipitation, the F2 layer peak electron density at Prof#4-6 recorded the largest value. However, in the adjacent space and time, the Nmf2 at Prof#7-9 is significantly lower than the former group, consistent with the ionospheric abnormal enhancement at two points. In [41], it was found that the ionosphere peak layer height (hmF2) enhanced abnormally before the earthquake. Based on the hmF2the of both event groups in the right panel, the extraordinary feature before the Concepcion earthquake was not prominent. The study of the variations of the ionosphere's hierarchical structure was continually carried out. For the 06-12 UT period on February 25 (two days before the earthquake), analysis was done on the electron density data for five heights from 100 to 500 km (interval: 100 km). Following the GNM construction, the same strategy was used to construct each height's electronic density maps to study the ionosphere disturbance at a different elevation. The results are shown in Figure 12. At 100 km, the distribution of RO events with drastic variations has no special rule. This may be related to the sporadic E-layer signatures or frequency collisions of electrons with a high-density neutral atmosphere. The imitation of the Abel inversion is another potential cause. The error accumulates as the profile goes to lower altitudes. Figure 12a shows what may be an effect of the error in the inversion process [42]; at the height of 200 km and 300 km, the electron concentration anomalies are consistent with the abnormal changes in the F2 layer. The vertical ionospheric structure is affected by the change in longitude and latitude. Especially for events with an extensive variation range, they are concentrated in the south area near the epicenter, slightly different from those in Figure 9. Generally, by the interaction of each ionization, the prominent long-wavelength fluctuation when traveling induced the amplification of the AGW effect in the F2 layer. The observation of the obvious local oscillation between 200 and 300 km in height here further confirms the idea. With the increase in height, the energy of the acoustic gravity wave gradually attenuated and the upper ionosphere gradually returned to a normal state. The study of the variations of the ionosphere's hierarchical structure was continually carried out. For the 06-12 UT period on February 25 (two days before the earthquake), analysis was done on the electron density data for five heights from 100 to 500 km (interval: 100 km). Following the GNM construction, the same strategy was used to construct each height's electronic density maps to study the ionosphere disturbance at a different elevation. The results are shown in Figure 12. At 100 km, the distribution of RO events with drastic variations has no special rule. This may be related to the sporadic E-layer signatures or frequency collisions of electrons with a high-density neutral atmosphere. The imitation of the Abel inversion is another potential cause. The error accumulates as the profile goes to lower altitudes. Figure 12a shows what may be an effect of the error in the inversion process [42]; at the height of 200 km and 300 km, the electron concentration anomalies are consistent with the abnormal changes in the F2 layer. The vertical ionospheric structure is affected by the change in longitude and latitude. Especially for events with an extensive variation range, they are concentrated in the south area near the epicenter, slightly different from those in Figure 9. Generally, by the interaction of each ionization, the prominent long-wavelength fluctuation when traveling induced the amplification of the AGW effect in the F2 layer. The observation of the obvious local oscillation between 200 and 300 km in height here further confirms the idea. With the increase in height, the energy of the acoustic gravity wave gradually attenuated and the upper ionosphere gradually returned to a normal state. In a specific coupling way, the satellite records the complex changes of the ionosphere over the seismogenic area. In [43], the "Lithosphere-atmosphere-ionosphere coupling" system was proposed to find convincing evidence of the ionosphere effect. At present, there are three main pathways of Lithosphere-Atmosphere-Ionosphere Coupling [44][45][46][47][48][49][50]: chemical, electromagnetic, and acoustic. The chemical pathway mainly refers to the influence of changes in geochemical parameters (gas or radon diffusion, change of water level, etc.,) on the ionosphere. Such abnormal phenomena are caused by gas release before the earthquake. Under the action of the crustal stress accumulation, the gas released from the fractured anode pore scatters to the surface. The ionization process is then activated to form the hydrated ion clusters at aerosol size (1-3 um) during ion-induced nucleation [1]. Based on the aerosol optical depth (AOD) database, the abnormal aerosol changes of the event were detected [41]. The radioactivity of the aerosol ion clusters causes neutral gas ionization. Consequently, the physical properties and electrical balance In a specific coupling way, the satellite records the complex changes of the ionosphere over the seismogenic area. In [43], the "Lithosphere-atmosphere-ionosphere coupling" system was proposed to find convincing evidence of the ionosphere effect. At present, there are three main pathways of Lithosphere-Atmosphere-Ionosphere Coupling [44][45][46][47][48][49][50]: chemical, electromagnetic, and acoustic. The chemical pathway mainly refers to the influence of changes in geochemical parameters (gas or radon diffusion, change of water level, etc.,) on the ionosphere. Such abnormal phenomena are caused by gas release before the earthquake. Under the action of the crustal stress accumulation, the gas released from the fractured anode pore scatters to the surface. The ionization process is then activated to form the hydrated ion clusters at aerosol size (1-3 um) during ion-induced nucleation [1]. Based on the aerosol optical depth (AOD) database, the abnormal aerosol changes of the event were detected [41]. The radioactivity of the aerosol ion clusters causes neutral gas ionization. Consequently, the physical properties and electrical balance change by ion introduction in the local atmosphere. The reaction leads to increasing or decreasing conductivity [51]. The electromagnetic path mainly refers to the low-frequency electromagnetic radiation caused by the microstructure of rocks or by electromagnetic noise in the seismogenic area. The inducement is produced by the release of electric charge onto the fracture surface of the seismic source and affects the ionosphere [46]. These electromagnetic signals propagate along with the atmospheric boundary layer and penetrate the orbit height of the satellite, which is directly recorded by the satellite. The recording of the electromagnetic signal before the Concepcion earthquake was completed by the Demeter satellite [52,53]. In the near-earthquake stage, there was a noise band in the electric field radiation of the VLF near the epicenter, and the radiation appeared in the cross of magnetic conjugate region along the magnetic field line. The strong electric field radiation affects the plasma distribution in the ionosphere [54].
These two coupling paths, shown in Figure 13, were detected early by the remote sensing satellite (Terra/Aqua) and the electromagnetic satellite (DEMETER), respectively [39,52]. The acoustic coupling path was further observed based on the data of the COSMIC occultation. Because of the NOAA/AVHRR satellite thermal images in operation, the infrared excitation from the earth's crust in Northeast China with the largest linear structures and fault system was captured. The double variations of CO 2 concentration and the increase in temperature of about 3 • were subsequently recorded. The greenhouse gases excited the acoustic-gravity waves, and the changes in the intensity of the optical ionospheric glow lasted dozens of minutes. Coincidentally, the ionospheric disturbances appeared with spatial and temporal scales, corresponding to the spectrum of AGW excitation. In the RO observations, the seismicity influence on the ionospheric anomaly in the F region by AGWs is expected to be confirmed. If the development of the AGW effect is analyzed in time, we may find the outlier fluctuation of the electron density over the seismically active area in the transmission process of the upward AGW waves. A similar response can be evidenced by TEC variations from GIMs prior to the Concepcion earthquake in Figure 14.
The TEC anomalous variations before the Mw 8.8 Concepcion earthquake on February 27, 2010, were analyzed. As shown in Figure 14, the highest daily total electron content usually occurs from universal time (UT) 16:00 to UT 20:00 (about local time (LT) 11:00- According to the 3D numerical model of the excitation of AGWs [10], the lithospheric gas source and propagation into the ionosphere were investigated. If a bell-shaped lithospheric gas source is characterized with the dimensions of 100 × 100 km 2 , the AGW excitation can cause a maximum change in the F region over 70 min. The peaks of the AGW wavelength and relative variations to the electron density (Ne) then reach dozens of km and percentage points, respectively. In Figures 10 and 11, the vertical ionosphere structure has a large vibration, covering from the bottom to the upper ionosphere. The oscillation amplitude increases with the uplift of height in a similar effect as the acoustic gravity wave. A clear feature set of the RO fluctuation, namely, a duration of around 1 h (14-15 UTC), a wavelength length of >50 km, a location at the F layer, and a relative variation (δ > 0.5), is consistent with these AGW descriptions. On this basis, we infer that the dissemination of the lithosphere's rupture may induce the acoustic gravity wave upward to the atmosphere. The reactions alter the temperature and density of the plasma layer and the subsequent structural damage of the ionosphere. Therefore, the ionosphere anomalies can also be characterized for AGWs in observed profiles prior to quakes, in addition to after the event. Moreover, AGWs can even be responsible for the change in electromagnetic wave characteristics in the "Earth-Ionosphere" waveguide. Today, plasma instability caused by AGWs excited by lithospheric gas is still worth studying. To satisfy the unified physical phenomena observed before earthquakes, different channels are under development in statistical cases and specific quakes. If RO profiles are amassed (the COSMIC-2 system in full operation), the time series at each profile location will further confirm the AGW coupling.
If the development of the AGW effect is analyzed in time, we may find the outlier fluctuation of the electron density over the seismically active area in the transmission process of the upward AGW waves. A similar response can be evidenced by TEC variations from GIMs prior to the Concepcion earthquake in Figure 14. ground-based GNSS monitoring. Due to the sparsity of COSMIC products, the detailed mechanism by which AGWs trigger PEIAs in dissemination is yet to be determined. With the arrival of advanced COSMIC-2 products, its capability of producing ~5000 high-accuracy soundings per day will promote the in-depth excavation of Lithosphere-Atmosphere-Ionosphere Coupling. Relevant investigations on the specific propagation mode of AGWs, the ionospheric traveling disturbance, and the damage degree of the ionosphere's structure are expected in the future.

Conclusions
In this paper, the global RO profiles from the FORMOSAT-3/COSMIC system were used to construct global maps of the ionosphere electron density 15 days before Chile's The TEC anomalous variations before the Mw 8.8 Concepcion earthquake on 27 February 2010, were analyzed. As shown in Figure 14, the highest daily total electron content usually occurs from universal time (UT) 16:00 to UT 20:00 (about local time (LT) 11:00-15:00). A magnitude increases to more than 25 TECU during active conditions. During quieter conditions, the average TEC usually ranges from 15 to 20 TECU (Figure 14a). In addition, we found that the ionospheric TECs in the area of the epicenter fluctuate dramatically in the approximate period with Nmf2 disturbances. Obvious TEC anomalies were observed before the shocks, on February 15,20,21,and 25, and all the ionospheric TEC anomalies were positive ( Figure 14b). As shown in Figure 14c,d, on 21 February 2010 (six days before the shock), an apparent TEC increase in the southeast region of Chile and the Weddell Sea from UT 10:00 to 12:00 was captured. The anomalous amplitude was about 5 TECU. Four days later, on 25 February 2010, a strong, positive, TEC anomaly (>6 TECU) was observed over the region of the seismogenic zone from UT 08:00 to 10:00. The ionospheric variable characteristics before the Concepcion earthquake conform to the LAIC model with an identical spatial distribution of TEC and Nmf2 disturbance peaks.
There are differences in the pre-earthquake ionosphere anomalies (PEIA) recorded in GIM and COSMIC. Compared with COSMIC occultation observations, the spatio-temporal continuity of GIMs (in the presentation of the planar ionosphere structure) is obviously better. The size and intensity of the TEC anomalies, however, are uncertain, especially in the ocean [55]. In addition to the difficulty in determining the TEC boundaries, the lack of GNSS sites in marine areas gives rise to another major challenge. By contrast, the COSMIC profiles can effectively reflect the ionospheric structure of Nmf2 and the electron density at different heights, including the oscillation mode and the uplift or descent in the F2 layer. These RO events are conducive to capturing weak traces of PEIA atmospheric dynamics but are not conducive to collecting commotion signals by traditional ground-based GNSS monitoring. Due to the sparsity of COSMIC products, the detailed mechanism by which AGWs trigger PEIAs in dissemination is yet to be determined. With the arrival of advanced COSMIC-2 products, its capability of producing~5000 high-accuracy soundings per day will promote the in-depth excavation of Lithosphere-Atmosphere-Ionosphere Coupling. Relevant investigations on the specific propagation mode of AGWs, the ionospheric traveling disturbance, and the damage degree of the ionosphere's structure are expected in the future.

Conclusions
In this paper, the global RO profiles from the FORMOSAT-3/COSMIC system were used to construct global maps of the ionosphere electron density 15 days before Chile's Concepcion earthquake (27 February 2010). By kriging interpolation, the disturbances of the Nmf2 parameter and the ionosphere's structure at different heights in the 15 days before the earthquake were compared and analyzed. The results show a noticeable disturbance in the F2 layer near the epicenter on February 21 and 25, and the abnormal area had an apparent local effect. Compared with other periods, the maximum disturbing momentum of Nmf2 over the epicenter of the two days was more than 50%. The ionosphere anomaly had an absolute directional difference. The disturbance peak point appeared in the north of the epicenter on February 25, with the disturbance point occurring in the southeast of the epicenter on February 21. Simultaneously, the electronic density profiles at these COSMIC occultation points were detected to have irregular oscillations below the F layer. The difference in the vertical structure of the ionosphere resulted in different responses to earthquakes at separate heights. The local effect was the most obvious between 200 and 300 km. The ionosphere disturbance at a height over 400 km gradually weakened. Influenced by the neutral atmosphere, the ionosphere in the bottom layer presented an irregular disturbance. Jointed by the similar response of GIM TEC, the traveling ionospheric disturbance was confirmed to appear across the bottom to the upper layer. These phenomena proved the possibility of the multi-channel coupling effect of the seismo-ionosphere.
Author Contributions: Conceptualization, writing and funding acquisition, K.S. and J.G.; investigation, validation, Y.Z. and W.L.; methodology and resources, Q.K.; software and formal analysis, T.Y. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare that they have no competing interests.

Appendix A
The following content in Appendix A is the supplemental descriptions in the Section 2 on how to control the quality of RO profiles.
1. The 9-point moving average proceeds to an electronic density (Ne) series from the bottom to the top height of the designated RO profile. After that, the average deviation between the Ne of each point and its initial value is where Ne i and Ne i are the observation and 9-point moving average of the electron density at point i above 110 km, and N is the sampling number per profile. For AD > 0.25, the unqualified profile is considered to be at a high-noise level [56]. 2. The hmF2 exists at a 150-450 km height: 150 km < hmF2 < 450 km (A2) 3. In the RO dataset, there are also profiles without obvious Nmf2 and hmf2, which is unreasonable and needs to be removed. For the quality control, two slopes (slope1 and slope2) at a height of 150-450 km are defined and should be satisfied as where h 1 (~450 km usually) and h 2 (~150 km usually) represent the top and the lowest height at each profile, respectively, and the Ne (h 1 ) and Ne (h 2 ) are the electron densities at the height of h 1 and h 2 , respectively. Here, the maximum electron density at the height 150-450 km defaults to the Nmf2 with a corresponding hmF2 as a maximum height in the F2 layer. 4. Except for meeting the conditions of (A1)-(A3), the difference between longitudes and latitudes at the heights of 100 and 500 km should be less than 5 • , and if not, it must be filtered. The condition (A4) is expressed as |Lon 100 km − Lon 500 km | < 5 • | Lat 100 km − Lat 500 km |< 5 • (A4)

Appendix B
The following content in Appendix B is the detailed descriptions in the Section 2 about the thorough operation of kriging interpolation. Based on the variability and correlation of variables, this method estimates the value of regionalized variables in a limited area in an unbiased optimal way; from the perspective of interpolation, it aims to find the linear optimal and unbiased internal estimation for scattered and sparse data. In the appendix, the technique will be described in the combination of variograms γ and prior sampling value Z(x i ), to determine the unknown values at different locations.

The ordinary kriging interpolation
This procedure uses known sample values Z(x i ), and variograms γ, to determine the unknown values at different locations as a linear combination of these known values. This technique can be written by the expression as: where Z(x 0 ) is the estimated value at the wanted position x 0 , Z(x i ) is the true value at the xi position, and λ i represents the weight function that only depends on the distance of the estimation by kriging equations and was restricted by: where µ is a Lagrange multiplier; γ ij and γ i0 represent the variogram values corresponding to the vectors between the observed positions (x i to x j ) and estimated positions (x i to x 0 ), respectively. These quantities can be calculated by fitting the experimental variogram with the theoretical model. Under the constraint condition of Equation (A6), λ j can be solved and then substituted into Equation (A5), in case the wanted values at the x 0 location can be interpolated from the prior vales of Z(x i ).

The calculation of variograms
The method for determining the variograms is based on minimizing the variance of the variable as: 2γ(|h|) = var(Z(x + h) − Z(x)) (A7) where var represents the variance function, γ is the variogram, and x and x + h are two close positions. By modeling its variogram under the assumption that Z is locally stationary, E(Z(x + h) − Z(x)) = 0 (A8) where E the expectation operator. In practice, the number of samples is always limited, so the variogram can be regarded as the experimental variogram composed of finite measurements, as follows: where γ*(h) stands for the estimation value of theoretical variogram γ*(h), and N is the pair number of experimental measurements with a certain distance (h) apart, minimizing the variance estimate with respect to the observations of Z(x i ). According to Equation (A5), the experimental variogram γ*(h) with the space distance (h) is calculated. In the following section, after the selection and fitting of the appropriate theoretical model is complete, the regional kriging interpolation is finalized toward the unknown Z(x 0 ).

The variograms with the Gaussian psill
According to whether there is a psill value, the variogram of the theoretical models can be divided into two categories. It is generally believed that the regional variables of the psill model meet the assumption of the second-order stationary assumption, while the other is only considered to satisfy the intrinsic hypothesis. Here, we adopt the Gaussian psill model with good behavior in spherical fitting and can be written as: where a is the range (set in advance), recorded by horizontal variogram coordinates beyond which a stability state is reached. Based on the γ*(h) (approximate to γ(|h|)) obtained in Equation (A5), the parameters C 0 and C are calculated using the least square fitting method, and the theoretical variograms model can be determined.
With respect to the observed value Z(x i ) with the restriction on the λ function (Equation (A1)), the system of kriging interpolation with the Gaussian psill is obtained.