© W. Jarmołowski et al., Published by EDP Sciences 2021

1 Introduction – state of the art in detection of seismic ionospheric disturbances

1.1 Ground GNSS and ionosonde data in the detection of ionospheric disturbances induced by earthquakes and tsunamis

The change of the ionosphere related with a seismic event can be observed in the parameters such as: total electron content (TEC), electron density (ED), temperature, ion density, and ion temperature. The atmospheric electric field generated on or near the ground surface during the seismic period is assumed to cause ionospheric anomalies (Pulinets et al. 2000, Pulinets & Davidenko, 2014). The vertical displacement in the atmosphere, induced by Rayleigh waves, produces an acoustic wave that propagates vertically to the ionosphere by dynamic coupling. During the time of upward propagation, the effect of the conservation of kinetic energy and decreasing of the atmospheric density strongly amplify this acoustic wave. This wave reaches the ionosphere and creates strong variations in the plasma density and plasma velocity. The propagation properties of Rayleigh waves can be estimated by observing the ionosphere, which is shown by Artru et al. (2004). Recently, an increasing number of observations on seismic ionospheric anomalies were reported, and these observations are based on different instrumentations (Jin et al., 2010, 2014; Occhipinti et al., 2013). GNSS is an efficient tool for observing traveling ionospheric disturbances (TID) and determining the 3D distribution of the ionospheric electron density (Liu (Tiger) et al., 2018). The earthquake (EQ) is associated with various electromagnetic signals, which might be detected by GNSS, also as EQ precursors (Ouzounov et al., 2018). The neutral atmospheric parameters and ionospheric TEC can be precisely estimated from ground-based and spaceborne GNSS, as well as by spaceborne GNSS radio occultation observations. The tropospheric and ionospheric delays derived from ground-based GNSS networks, precise orbit determination (POD) GNSS receivers onboard LEO satellites, and spaceborne GNSS Radio Occultation provide a chance to monitor and understand seismic ionospheric disturbances, which may give us a new perspective to observe and forecast the EQs in the future. Komjathy et al. (2012) show that the real-time global TEC monitoring network observing the acoustic and gravity waves generated by the EQ and/or tsunami could give potential warnings. In turn, Garcia et al. (2005) have even proposed a three-dimensional ionospheric tomography for the purpose of TID detection. Applying a complementary method, Maruyama et al. (2011) report tsunami-induced disturbances of the ionosphere by means of the ionosonde stations. Another example of ionosonde application to the analysis of seismically triggered ionosphere disturbances can be found in Liu & Sun (2011). Both above works refer to the largest M_w = 9.1 Tohoku event, however, Chum et al. (2016) uses ionosonde data as ancillary in the investigation of M_w = 7.8 EQ in Nepal 2015.

Co-seismic ionospheric disturbances have been investigated in numerous studies related to the EQs. The example is given in Dučić et al. (2003), who detected TEC perturbations related to the acoustic waves associated with the Rayleigh waves. These disturbances have been measured by GNSS receivers of dense Californian GNSS network following the Denali EQ in 2002. Lognonné et al. (2006) investigated the detection of seismic waves from ground GPS signal and referred to this and other EQs. Recently, more comprehensive studies of seismic ionospheric effects for the 2011 Tohoku-Oki EQ became possible using GNSS observations collected by around 1200 GEONET stations (Liu et al., 2011; Rolland et al., 2011; Tsugawa et al., 2011; Occhipinti et al., 2013). Astafyeva & Shults (2019) have even applied GEONET to the analysis of EQ source location and have presented possibilities and difficulties of such an approach. The extremely dense ground TEC observations performed by the GEONET network enabled clear distinguishing of acoustic-gravity waves generated at the epicenter, acoustic waves coupled with the Rayleigh wave, and gravity waves coupled with the following tsunami (Occhipinti et al., 2013; Jin et al., 2014; Chou et al., 2020). The southeastward ionospheric disturbances, possibly resulting from the coupling between tsunami and atmosphere, have also been detected (Galvan et al., 2012; Occhipinti et al., 2013). Different propagation velocities, which are related individually to the seismic Rayleigh waves, acoustic waves and tsunami-generated gravity waves were also discussed in several other works (Heki et al., 2006; Jin et al., 2014, 2015; Yang et al., 2019). Thus, we know that TEC disturbances spread out to almost all directions from the epicenter, with the speed of several hundreds of meters to several kilometers per second. Pre-seismic ionospheric anomalies appearing several hours to several days before the main shock have also been reported by lots of researchers as observed in TEC or other ionospheric parameters (Liu et al., 2004; Heki & Enomoto, 2013, De Santis et al., 2019; Goto et al., 2019, 2006a). Pre-seismic TIDs are often reported in time of the probability of a new EQ precursor. For instance, Liu et al. (2006a) showed changes of GPS TEC anomalies five days before the EQs having a magnitude larger than 5.0. It should be pointed that seismic TIDs have been statistically proved to be dependent on the magnitude and on the focal depth of EQs (Le et al., 2011). EQs with a focal depth of less than 60 km and magnitude greater than 6.0 provide presumably the major disturbances in ionospheric parameters because such EQ zones at shallower depth have a gigantic amount of energy. Su et al. (2013) reported that for shallow EQs (depth <40 km), the variation of TEC exhibited a significant increase before the EQs.

The relationship between tsunamis and ionospheric disturbances was firstly noticed by Hines (1972), who suggested that the ionospheric signature could be employed as an early warning of the occurrence of a tsunami. The reason for ionospheric signatures of a tsunami is a displacement of the atmosphere on the surface of the sea, which induces atmospheric waves that propagate up to the thermosphere, and subsequently generate amplified perturbations in the ionosphere. Liu et al. (2006b) reported significant theoretical and observational results of tsunami-driven ionospheric signatures caused by the perturbations of TEC after the deadliest Sumatra M_w = 9.2 EQ of 26 December 2004. Artru et al. (2005) used ground GNSS network data and detected small-scale ionospheric disturbances of gravity waves generated by the tsunami. This tsunami was triggered by Peru M_w = 8.2 EQ in 2001 and reached Japan several hours later. These ionospheric disturbances observed above Japan indicated that their time delay after the EQ coincides with the propagation time of the tsunami. Savastano et al. (2017) used a number of local GNSS receivers in Hawaii to detect tsunami-induced ionospheric anomalies after M_w = 7.7 Haida-Gwaii EQ in 2012. They had found typical variations in TEC when tsunami-induced ionospheric waves arrived on the islands. Rolland et al. (2010) analyzed perturbed ionospheric signal offshore Hawaii in the same way, using a limited number of receivers, but after three different EQ/tsunami events. They performed a preliminary spectral analysis of the detected disturbances. A significantly larger number and better spatial distribution of ground GNSS data have been applied by Azeem et al. (2017). They accurately presented the spatio-temporal evolution of TEC perturbations, which came to the USA coast and reached CORS GNSS network 12 h after Tohoku EQ.

The magnitude of EQ and time of perturbances, as well as their relation with ionospheric response size have been analyzed by Astafyeva et al. (2013). Perevalova et al. (2014) examined the EQs having various magnitudes in different regions of the Earth, and found that there was a threshold magnitude (around M_w = 6.5) below which there were no observed seismic TEC disturbances in ground GNSS data. Ke et al. (2016) performed another statistical analysis of seismo-ionospheric anomalies related to the EQs in China whose magnitude exceeds some selected level. The geomagnetic conditions also significantly influence the amplitude of TID (Heki & Ping, 2005; Afraimovich et al., 2010). The ionospheric disturbance amplitude is highly relative in disturbed geomagnetic days, and the filtered TEC series is affected by the geomagnetic variations. The coincidence of the tsunami TID source location and the tsunami location shows that the ionospheric TEC recorded by local ground-based GNSS receivers can be used to confirm the tsunami occurrence, find the tsunami location, and support the tsunami early warning (Liu et al., 2020).

1.2 Low-Earth-Orbit (LEO) satellite data in the detection of ionospheric disturbances induced by earthquakes and tsunamis

The example of a completed mission designed to monitor the variations in the atmosphere related to the EQs was the French DEMETER (Detection of Electromagnetic Emissions Transmitted from EQ Regions) mission (2004–2009). The scientific payloads on the DEMETER satellite were the Langmuir probe (LP), measuring ED and electron temperature, and the ion analyzer (IAP). Among other payloads, the Electric Field Instrument (fr. ICE) and magnetic search coil instrument (fr. IMSC) were responsible for the determination of the electric and magnetic field variations related to the EQs. DEMETER provided observations of electromagnetic fluctuations related to the EQs (Venkatraman & Heelis, 2000; Píša et al., 2011; Li & Parrot, 2013). DEMETER data have also been examined in various ways to define the perturbations of the ionosphere caused by acoustic gravity waves from the EQs (Anagnostopoulos et al., 2012; Píša et al. 2011; Li & Parrot 2013; Ryu et al. 2014). An extensive, statistical approach to DEMETER data by Parrot and Li (2018) has shown that significant number of satellite data can be useful in analyzing EQ precursors. It was concluded that a combination of satellite and ground data could be particularly powerful in detections and warnings. The precursory character of TIDs was also reported in the study combining GNSS and DEMETER data by Akhoondzadeh et al. (2010), Zhang et al. (2012), and Ibanga et al. (2018).

The existing missions, including LPs among the scientific payload, are the Swarm constellation (Diego et al., 2019; Luo et al., 2019) and a very recent China Seismo-Electromagnetic Satellite (CSES) (Liu et al., 2019). Swarm is the fifth Earth Explorer mission approved in ESA’s Living Planet Program and was successfully launched on November 22, 2013. The research objectives of the Swarm mission are to provide the best-ever survey of the geomagnetic field and electric field in the atmosphere using a constellation of three identical satellites carrying precise magnetometers and electric field instruments. It is very common for EQ-related Swarm studies to incorporate magnetic data obtained from these three satellites. There are published studies of precursory signals in magnetic data (Zhu et al., 2019) and as papers referring to co-seismic disturbances of the magnetic field (Akhoondzadeh, 2019; Marchetti et al., 2020a). However, Swarm’s payloads also include LP sensors and POD GPS receivers, and therefore seismic signals in the electric field were also investigated using these data. Stanica et al. (2018) studied jointly the variations of magnetic and electric parameters as precursory to some quite low magnitude EQs. A stronger event was investigated by Marchetti & Akhoondzadeh (2018), and pre-seismic and co-seismic characteristics in magnetic and ED data were shown. De Santis et al. (2019) prepared a statistical assessment of magnetic and ED disturbances related to 12 EQs with various magnitudes. Pulinets & Ouzounov (2018) have applied Swarm B ED data for the analysis of Chile-Illapel EQ and have provided preliminary observations of pre-seismic disturbances. Finally, Marchetti et al. (2020b) investigated possible lithosphere-atmosphere-ionosphere coupling effects prior to the M_w = 7.5 EQ using, among others, Swarm and CSES datasets.

The constellation of FORMOSAT-3/COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) satellites has been applied for TID detection in the ionosphere. Coisson et al. (2015) presented the first experimental radio occultation observation of a gravity wave excited by a tsunami during the 2011 Tohoku event. Tulasi Ram et al. (2017) gathered simultaneous observations of ionospheric perturbations from ground-based GPS vertical TEC (VTEC) and COSMIC radio occultation measurements. They also estimated Rayleigh wave velocity by the application of spectral methods to ground GNSS data. Sun (2019), in turn, analyzed ground-based GNSS TEC and ED from FORMOSAT-3/COSMIC and prepared spectral analysis also for the perturbations in COSMIC data. There are also some evidence of the usability of other LEO signals in detecting seismically triggered waves in the ionosphere, mainly in reference to Tohoku 2011, which due to its magnitude, still stands as the most effective case study. Garcia et al. (2014) have presented high-frequency atmospheric density perturbations of GOCE accelerometer and thruster, at the height of 270 km, induced by gravity waves created by the Tohoku-Oki tsunami propagating in the Pacific Ocean. Yang et al. (2014) showed that significant TEC fluctuations could be observed by the K/Ka-band ranging (KBR) signals approximately 20 min after the mainshock of the Tohoku-Oki EQ. At that time, the GRACE satellites were orbiting the Earth at an altitude of 450 km.

2 Chile-Illapel earthquake and tsunami of September 16, 2015: the largest event in time of Swarm, but little explored in the ionosphere from satellite and ground data

2.1 The characteristics of Chile-Illapel region, earthquake and tsunami

The largest EQ/tsunami event during the Swarm mission was recorded in Chile. The EQ on September 16 in 2015, occurred at 22:54 UTC. The EQ was located 46 km offshore from Illapel in Chile, which is also 169 km northerly from Valparaiso in Chile (31.570°S, 71.654°W, 24.9 km depth – USGS location). This EQ was preceded by several significant EQs, which have indicated increased seismic activity in the Chilean region before the mainshock (Fig. 1). There were also a significant number of aftershocks that occurred after M_w = 8.3 EQ, and the strong seismic activity persisted there for five days. In general, Figure 1b presents a number of EQs close to the time of mainshock, which has magnitudes starting from 5.0 and is placed no more than 1000 km away from Chile-Illapel. All these EQs are, in fact, located along the same subducting edge of the South American and Nazca tectonic plates (see Figs. 1a and 2).

Fig. 1 Seismic activity from magnitude 5.0 in region of South America (close to Chile-Illapel) during September 2015. (a) Map of EQs in September 2015 in the neighborhood of Chile-Illapel. (b) Plot of EQs with respect to mainshock. Magenta stems denote EQ magnitude multiplied by 10, brown stems denote hypocenter depth. Source of subfigure (a) and data: USGS.

Fig. 2 Residuals from TG and DART readings of Chile Illapel 2015 tsunami (September 16th, 2015) recorded 15.5 h after source EQ.

The tsunami had started from the mainshock (M_w = 8.3) at 22:54 UTC on September 16, 2015. Its wave is the maximum found among all recorded at tide-gauges by the NOAA Tsunami Warning Center, and it reaches 467 cm. Chile-Illapel tsunami event has been recorded at 156 tide-gauges and buoys. Figure 2 presents a snapshot moment of tsunami wave recorded by many tide-gauges in the afternoon of September 17, 2015. The inducing EQ was late in the local evening of the previous day, so the tsunami lasted all day of September 17. The times of arrival of this tsunami can be observed in Figure 3, where we can find that the whole next day (September 17, 2015) was necessary for its travel across the Pacific Ocean.

Fig. 3 Travel Time Map of the Chilean 2015 tsunami (September 16th, 2015). Source: NOAA, U.S.TWS: https://www.tsunami.gov/previous.events/?p=09-16-2015-Chile(M8.3).

The analyzed Chile region is located at the western border of the South-American tectonic plate, which is rich in a number of smaller, minor plates like the North Andes Plate and Altiplano Plate (Fig. 2). The western edge of the South-American plate is a convergent boundary with the subducting Nazca Plate. The eastward-moving and more dense Nazca Plate is subducting under the South-American Plate, along the Pacific coast, at a rate of 7.7 cm per year. Altiplano Plate is a small tectonic plate, which was uplifted primarily because of crustal thickening produced by horizontal shortening of a thermally softened lithosphere. Altiplano Plate is located between the Nazca Plate in the southwest (the border is the Atacama Trench) and the South-American Plate that surrounds it from the other sides (Fig. 2).

2.2 Study objective – the ionosphere above Chile-Illapel EQ/tsunami revisited by Swarm, GRACE, ground GNSS, and ionosondes

Subsequent to the Tohoku 2011 EQ/tsunami, the ionosphere in the vicinity of the epicenter has been intensively investigated, whereas only single preliminary works focused on the ionospheric responses to Chile-Illapel EQ have been shown. The figures of concentric TIDs after Tohoku appear in many studies (Jin et al., 2014; Yang et al., 2019; Chou et al., 2020). The Chile EQ/tsunami event in 2015 has been investigated notably less frequently than the Tohoku event, also due to the smaller impact on the ionosphere. However, it can be suspected that more sparse GNSS networks in Chile and Argentina can also contribute to this fact. It was quite transparently indicated in Section 1.1 that ground GNSS data, due to several tens of satellites and thousands of stations, are the most encouraging to research on seismic TIDs. This fact comes from one main reason, i.e., ground GNSS networks provide the only reliable spatial view on TIDs and enable to inspect their spatial shapes, spatial relations, and speeds. However, ground GNSS stations have a heterogeneous distribution worldwide, and only a few seismically interesting active zones are rich in dense ground GNSS networks.

Sparse distribution of ground GNSS stations significantly impedes the identification of TIDs in many interesting areas and limits studies on their correlation in a wider range. By means of tens of receivers from the Centro Seismológico Nacional (CSN) and the International GNSS Service (IGS) network, Reddy et al. (2017) report the characteristics of ionospheric perturbations induced by shock acoustic waves of the EQ/tsunami, proving the advantages of GPS observation networks in imaging tsunami-driven ionospheric perturbations. Nevertheless, they strongly limit the area, a number of stations, and GNSS satellites in the analysis, and from their study, it was difficult to discern different velocities of the ionospheric TEC perturbations. Also, the timeframe of their study was limited to the mainshock, which is typical in the studies of seismic ionospheric disturbances. Grawe & Makela (2017) show ionospheric perturbations appearing in filtered GPS-derived TEC determined at Hawaii stations after Chile-Illapel EQ, however, these disturbances are found with respect to only one GPS satellite. Lately, Ravanelli et al. (2021) have studied Chile-Illapel event with the use of ground GNSS data. They found the velocities of the gravity waves triggered by the mainshock using the keograms based on the data from selected GNSS PRNs. They have also shown some preliminary scatter representation of these ionospheric disturbances for the first time, but it is suspected that the use of many epochs for displaying by the scatter makes it not very accurate. The current study analyzes more moments related to the Chile event, does not select individual PRNs for the keograms, and shows more evident seismically triggered ionospheric disturbances in the scatters and keograms based on ground GNSS data presents also seismic TIDs along the LEO tracks. Pulinets & Ouzounov (2018) have analyzed Swarm B along-track ED data related to Chile-Illapel EQ, and they have found the largest disturbance two days before the mainshock. The analyses of Swarm, other LEO data and ground data altogether, can provide interesting observations of ionospheric disturbances not only directly related to the mainshock and tsunami, but also occurring during entire periods of enhanced seismic activity and at larger distances from EQ epicenter.

The objective of this study is to explore more the Chile-Illapel EQ and tsunami on September 16, 2015 and extend the works that can be found e.g., in Reddy et al. (2017) or He & Heki (2016). The current study analyzes a much longer time period referring to the time before and after the M_w = 8.3 EQ, which can give more chances for observing potential pre-seismic and post-seismic signals and provide additional observations useful for the potential detection and warning system (Komjathy et al., 2012). A wider area with respect to mainshock epicenter location is analyzed, including the entire area of South America, southern North America, and Pacific Ocean. A significant extension of this Chile-Illapel 2015 analysis is based on several satellite and ground data types, with a special focus on Swarm and ground GNSS, including GRACE, ionosondes, and seismic records. The tomographic voxel model is also applied in the research, and 4D ED gradients are extracted to reveal the post-seismic TIDs gathered in the model. The satellite data used in the research are Swarm LP in-situ ED, topside TEC from Swarm POD GNSS, and TEC difference based on dual-frequency precise ranging between GRACE satellites. The ground data applied in this study are ground GNSS data from around 200 stations belonging to the Chilean CSN network and Argentina Instituto Geográfico Nacional (IGN) network, as well as two Digisonde stations located in nearby locations with respect to Chile-Illapel event epicenter and evolution, as well as with respect to interesting Swarm passes.

The presented study includes a number of mathematical techniques, including also spectral analyses by short-term Fourier transform (STFT) and continuous wavelet transform (CWT), which are here applied for the first time to Swarm (LP and POD TEC) and GRACE KBR along-track data. Also, combined lately developed techniques (Yang et al., 2019) have been applied in the detection of TIDs. The processing of extended scope of the data with advanced or composed techniques enabled the detection of many interesting pre-seismic and post-seismic ionospheric disturbances of different sizes, at different distances from the event and within different timeframes. Summarizing, this work includes a review of detectable TIDs for an extended timeframe, in a close neighborhood of EQ and more distant places, and proves that wide area, different altitudes, and extended timeframe are interesting in the analysis of EQ/tsunami event in the context of future prompt detection of the tsunami.

The Chile EQ/tsunami event in 2015 occurred in a period when the magnetosphere was recovered from three consecutive geomagnetic storms on 7, 9, and 12 September 2015, and the values of the Dst-index were still at the level of −20 nT to −30 nT. Such conditions might trigger perturbation in the ionosphere. To examine the state of the ionosphere, the ionospheric characteristics observed from the Cachoeira Paulista Digisonde station are presented in Figure 4. Cachoeira Paulista Digisonde is in the southern hemisphere, at low latitudes (−22.70 S, 315.0 E), at a relatively short distance from the EQ epicenter, and it is expected that large scale ionospheric conditions recorded by this Digisonde station, are representative of the conditions in the wider region of the low latitude zone in the American continent. The data are retrieved by the GIRO (Galkin et al., 2006) and are manually scaled for higher accuracy. Figure 4 presents the relative deviation together with the relative standard deviation of the foF2 characteristic for reference. This time plot verifies that ionospheric conditions were quiet during the days around the EQ. The relative deviation and the relative standard deviation of TEC from GNSS receivers calculated for Cachoeira Paulista location with the UQRG model (Hernández-Pajares et al., 1999) are also shown in Figure 4. For a few days before and after the EQ, no significant perturbation in TEC and foF2 is observed. A systematic increase in the ionization, seen as a positive deviation in foF2 and TEC, from 24 September until the end of the month is the result of the fast-growing value of the F10.7 index. However, this occurs 9 days after the EQ, and it does not affect the results presented in this paper.

Fig. 4 From top to bottom: the Dst index, and the relative deviation of TEC and foF2 values with respect to monthly medians (black lines) over Cachoeira Paulista (Brazil) during the whole month of September 2015. The color lines represent the relative standard deviation (+/−) of the values taken into account for the estimation of the medians in each case. The latter is included as an indicator of the ionospheric variability in normal conditions that are used as a threshold for identifying significant ionospheric disturbances in a large scale perspective.

3 Numerical study by applying data combinations

3.1 Close to the epicenter

3.1.1 Detection of TEC disturbances from the network of ground GNSS stations

The composed process for characterizing the TIDs associated with EQ and tsunami is done in two stages (see details in Yang et al., 2019). The first stage consists of the lines of sight computation and GNSS data preprocessing, i.e., obtaining the ionospheric combination from the GNSS carrier phase data, filtering the trend and bias, and mapping the detrended VTEC data, see details in Hernández-Pajares et al. (2006). The technique of detrending employs double-time difference to eliminate the trends in the ionospheric combinations. The selection of the double-time difference delays is justified on the basis of Rolland et al. (2011) and Galvan et al. 2012. Therefore, in this work, we selected double-difference periods of 60 s. and 300 s. The 300-s interval allows for a high sensibility to the tsunami-induced TIDs with the dominant periods in the range of 400–1200 s. The 60-s interval is used to filter the trend and preserve the TIDs with the periods centered around the range of 80–240 s. The TIDs in this range are mostly triggered by the Rayleigh waves and the corresponding acoustic resonance after the EQ. In order to better characterize the gravity waves, we assumed that the most frequent dominant TIDs produced by the interaction between the neutral and ion particles occur below hmF2 (see Hernández-Pajares et al., 2006), and thus we took 200 km as the mean effective height of the TIDs, which is lower than the average hmF2.

Chile-Illapel EQ and tsunami induced the ionospheric perturbations recorded by tens of permanent GNSS receivers, i.e., the Chilean CSN network and Argentinian IGN network. By combining these two GNSS networks in the region surrounding the event, we have almost 200 stations at disposal, but the exact number depends on the investigated area. This is much less in comparison to the Tohoku 2011 case, where the density of IPPs obtained from 1200 GEONET stations distributed over Japan is large. Nevertheless, this study aims to prove the detectability of clear and informative observations also from a smaller number and density of stations. In this manuscript, the detailed propagation characteristics of the EQ/tsunami-driven TIDs are studied, and we can show spatial shapes and temporal evolution of ionospheric disturbances. The quantity providing elementary observations of ionospheric disturbances is L_G-F = L₁ − L₂, the geometry-free linear combination of both frequencies carrier phase measurements (also known as ionospheric combination L_I, see Hernández-Pajares et al., 2006, 2011; Krypiak-Gregorczyk & Wielgosz, 2018). Figure 5 presents the region close to the epicenter, divided into four quadrants, according to ADDTID technique assumptions. The division is based on the assumption of different propagation directions and expected characteristics of detected TIDs in different types of area i.e., land and ocean.

Fig. 5 GNSS receiver distribution in Chile and Argentina. Cyan crosses denote the location of the GNSS receivers, and the red star denotes the epicenter of the Mw 8.3 EQ, on September 16, 2015. The grid cells (a)–(d) denote four subnetworks with a geographical grid of 10° × 10°. The ocean depth is shown by the cold color map, originating from the Earth observatory of the US NASA blue marble with topography and bathymetry.

The IPP trajectories at 200 km altitude, recorded from 9 min before the EQ to 21 min after the EQ, are shown in Figure 6a. Each color represents the IPPs set from individual GPS satellites. The location of the epicenter is indicated as a red star. The IPPs in green, denoting measurements to the GPS satellite of PRN 12, have recorded the ionospheric variation around the epicenter. For this satellite, the observational data from GNSS receivers whose IPP locations are at a distance of less than 250 km from the epicenter (marked by the black circle in Figure 6a), are selected in order to describe the relationship between the ionospheric perturbation and EQ and tsunami. The receiver observations were also filtered out under the losses of lock on GPS. Figure 6b depicts the slant TEC from these receivers to PRN 12 by means of L_G-F in length units, which are arbitrarily aligned in order to clarify the ionospheric variations. Some of arbitrarily aligned STECs around the epicenter show a large “N-shaped” TEC depletion of about 2 TECU appeared about 9 min after the EQ, in agreement with the 10 min of the arrival time of the acoustic waves associated with the seismic signal up to 221 km height reported in Artru et al. (2004). These anomalies are similar to a seismic signature and the ionospheric disturbances following the depletion. Yang et al. (2019) also reported similar observations of the sudden decrease of TEC as the response of the ionosphere to the 2011 Japan Tohoku EQ/tsunami.

Fig. 6 (a) IPP locations at the height of 200 km, corresponding to the GNSS receivers, from 22:45 to 23:15 UTC, the color codes of cycles for the IPP locations of different GPS satellites, the red star for the location of epicenter, (b) time evolution of the arbitrarily aligned slant TECs, for the GPS satellite of PRN 12, at the time of 22:30–24:00 UTC, the color codes of lines correspond to different GNSS receivers, the black dash line indicates the time of the EQ at 22:54 UTC.

A typical way to describe TIDs spatiotemporally and estimate their speeds is the keogram of the disturbances, presented in Section 3.1.3, together with the signatures of ionospheric disturbances detected by Swarm LP ED data. The spatiotemporal comparison of disturbances in LEO along-track and ground GNSS data enables to prove sensitivity of LP observations and validate detected ionospheric anomalies onboard Swarm satellites.

3.1.2 Detection and characterization of TIDs during 2015 Chile EQ/Tsunami by ADDTID

In the second stage, we detect and characterize the parameters of TIDs by means of the ADDTID algorithm. ADDTID (see details in Yang et al., 2017, 2019) uses data from a dense GNSS receiver network and can detect and characterize TIDs from the IPPs. The algorithm creates a dictionary of possible waves that encompass the geographical area under study. The estimation is done using only the directly measured detrended VTEC at the expected spectral range and IPP locations and solving a convex optimization problem with a regularization term. This algorithm determines one or few TIDs as the elements of the dictionary that best approximate the observations. Despite a small number of stations, the TIDs during the 2015 Chile EQ/Tsunami can be studied by the ADDTID method (see also its application to dense GNSS network in Yang et al., 2017, 2019). In order to analyze the local behavior of TIDs, the network has been divided into 4 geographical grids of 10°× 10° (Fig. 5). The propagation parameters of the local TIDs from one hour before the EQ to one hour after the EQ, are described in the polar plots in Figure 7. Figure 8 describes the time evolution of the estimated TID velocity distribution from 22:00 to 24:00 UTC for the subnetworks (a), (b), (c), and (d), respectively. These plots show MSTIDs moving at the highest velocities (above 300 m/s) towards the deep sea and further away from the epicenter, fully in agreement with the expected tsunami-related ionospheric signatures. It should be noted that the intensity of the TIDs is taken as the weights for the heights of histograms, which then can be easily read. The results estimated by ADDTID clearly show a set of strong TIDs appearing a few min after the EQ and later (see dark blue and light blue points, in Fig. 7). Also, the occurrences of the strong local TIDs in the grids of the network (subplots (a) and (c) pointing mainly to the sea) can give a clue to indicate the current location of the tsunami.

Fig. 7 Detected MSTIDs of each subnetwork. Parameter evolution of detected TIDs before the EQ during 22:00–24:00 UTC, representing non-uniformly scaled polar plots for TIDs of the corresponding grid networks (a)–(d), the location representing for the azimuth vs. velocity of 0–360 degree and 0–900 m/s, the circle size for wavelength of 50–600 km, and the color code for different time intervals.

Fig. 8 Time evolution of TID velocity distribution during 22:00–24:00 UTC, for the subnetwork (a), (b), (c), and (d), respectively. The scale of the velocity axis ranges from 0 m/s to 600 m/s.

3.1.3 Comparison of detected pre-seismic and co-seismic ionosphere disturbances by Swarm LP ED, Swarm POD TEC, and ground GNSS data

This section analyzes the ionosphere disturbances detected mainly by Swarm LP ED and POD TEC in three time periods, selecting the time when interesting Swarm passes over the Chile-Illapel vicinity are available. The EQ of magnitude 8.3 took place on September 16, 2015 at 22:54, however, the orbits of Swarm B just after this EQ were more than 2000 km to the west of this event, and Swarm A/C orbits were distant even more than 3000 km. This is too far for the observation of starting EQ/tsunami event. Swarm satellites A/C passed this region on September 16, 2015 around 20:00 UTC, whereas Swarm B passed the region around 22:00 UTC. This latter track of Swarm B is investigated later in Section 3.1.3 in detecting close pre-seismic TIDs. The closest available Swarm track after the mainshock in terms of close time and location is recorded on September 17, 2015, 7:25 AM (Fig. 19), and this time also corresponds to tsunami propagation time. This pass of Swarm A took place around 8.5 h after M_w = 8.3 EQ when there were many aftershocks, and the tsunami was in the middle phase of the propagation in the ocean.

LEO along-track data are compound 1D signals that can be analyzed by 1D STFT, an evolutionary form of discrete Fourier transform (DFT) (Harris, 1978). In this study, Swarm LP ED and POD TEC signals are high-pass filtered by DFT. The classification of the spectral characteristics of disturbing along-track signals by STFT is supported by their simultaneous search in ground GNSS observations, which gives an opportunity for the validation of spectral recognition. STFT spectral approach to along-track Swarm data gives an opportunity for distinguishing the signals of different origins. Firstly, we analyze Swarm B pass in UTC evening on September 10, 2015, which is six days before the mainshock. Secondly, we study Swarm B at the moment, just before the Chile-Illapel mainshock and analyze ground GNSS data also up to an hour after this EQ on September 16, 2015. Thirdly, we inspect the closest possible Swarm A/C pass after the mainshock in the morning of September 17, 2015. At these three moments, Swarm detection of a possible TID and its validation by ground GNSS data is made in the following three steps. Firstly, STFT spectrograms are made for the selected LP ED high-pass filtered observations, using 70 s. Tukey window. The ED measured by LP and POD TEC is high-pass filtered by DFT decomposition, and the signals that have a period larger than 70 s along the track at its altitude are removed. This time interval corresponds to around 500 km or 4.5° along the Swarm orbit, which is assumed as the upper limit for disturbances interesting from the seismic point of view (see e.g., Rolland et al., 2011; Chou et al., 2020). Additionally, we can see in Pulinets & Ouzounov (2018, p. 4–25) that this period seems reasonable and large enough with respect to the seismic disturbance observed by them on September 14, 2015, and small enough with respect to Equatorial Ionospheric Anomaly (EIA) crests and Sun-induced features of ED along-track record. It can be noticed at first look that detected potential anomalies on September 14, 2015 are of higher frequency than different characteristic components of EIA. It should be stressed that Fourier-based filtering avoids the remaining long-wavelength features of the signal with respect to gradient-based methods, which are approximate filters. The PSD from the spectrogram is sampled at a frequency of 50 s in order to visualize large PSD values along the Swarm track together with a scatter plot of ground GNSS data. Secondly, the keogram based on band-pass filtered L_G-F combination of phase data (Hernández-Pajares et al., 2006, 2011; Krypiak-Gregorczyk & Wielgosz, 2018), and distances calculated from the epicenter is prepared and overlayed with the selected Swarm LP ED spectrogram section in order to identify the potential spatiotemporal correlation of TID in these two data sets. The L_G-F is band-pass filtered by DFT between the period from 1800 s to 240 s, which is a slightly wider range than in Section 3.1.1. However, a larger interval is intentional for the overview keograms of larger areas and timeframes, where we expect a wider spectrum of disturbances to recognize. The primary assessment of character and origin of disturbances can be found in Section 3.1.2 and continued here with the help of keograms and scatters. In the third step, the scatter plots describe several min. of L_G-F from ground GNSS data and respective several min. of ED spectrogram section along the Swarm orbit are directly spatiotemporally compared. The disturbances observed by Swarm and reflected in large PSD along the orbital track can be directly visually compared to TIDs present in ground GNSS observations. TIDs in-ground GNSS data are visible as colored strips (e.g., Fig. 9) with positive slope, which travels outward from the seismic sources, with all velocities around 0.3 km/s., which makes them most similar to gravity waves (Galvan et al., 2012; Jin et al., 2014; Occhipinti et al., 2013). The disturbances in Swarm ED can be found at the selected frequencies of spectral patterns. It should be stressed that IPPs, which indicate the location of L_G-F phase combination data are referred in this case to the altitude of 450 km, which gives coincidence with in-situ LP data from Swarm A, and approximate coincidence with slightly higher orbit of Swarm B.

Fig. 9 Keogram of bandpass filtered LG-F phase combination from ground GNSS stations of Chilean and Argentinian networks, recalculated to TEC units, in the UTC evening on September 10, 2015. All distances of IPP tracks are calculated with respect to the position of the mainshock. Large dots denote Swarm A and B passes with distances respective to the mainshock location, and their color indicate size of power PSD sampled at 50 s. from the spectrogram. The blue-red strips denote different TIDs.

The sample day September 10, 2015, where evidence of medium-scale ionospheric irregularities was found, precedes the mainshock by six days. As investigated in Section 1.1, different authors select different periods before the main EQs to observe potential pre-seismic ionospheric anomalies (see e.g., De Santis et al., 2019; Liu et al., 2006b), and the discussion in this field is open. It is also possible that these observations still correspond to the local winter, which shows typical MSTIDs climatology during daytime and equatorward propagation (see for instance, Hernández-Pajares et al., 2012). The selection of September 10, 2015 is based on the preliminary observation of simultaneous Swarm LP and ground GNSS anomalies in the keogram in Figure 9, which presents blue-red strips at different distances from the mainshock location. Some ionospheric disturbances in Figure 9 coincide with a larger Swarm PSD at period 35–40 s, indicating in-situ ED disturbance. The spatial directions of the GNSS satellite-station phase measurements are often do not coinciding with TIDs. Therefore, the larger number of ground GNSS TEC determination points (IPPs), the better are TIDs observations. There are two passes of two Swarm satellites accessible in the time epoch shown in Figure 9. The disturbances in-ground GNSS data at the time of Swarm A pass slightly disappear for some time and come back after this pass. However, Swarm B crosses exactly the TIDs just after 22:00, when the preceding storm from 9.09 calms down, and the Dst index is around −20 nT (Fig. 4, top). The relative TEC changes are also small at this time (Fig. 4, middle). After signal high-pass filtering to 70 s, we can find a disturbance in Swarm LP ED, which can be analyzed by the spectrogram (Fig. 10c). This disturbance has a frequency starting from around 40 s, which then simultaneously decreases and increases, so its frequency splits into different frequencies. We can see in Figure 10b that TEC anomaly over the Altiplano minor plate is separated from EIA and that ED anomalies detected by Swarm are not necessarily co-located with the highest TEC anomalies.

Fig. 10 (a) Swarm B LP high-pass filtered ED data, (b) STFT spectrogram of residual ED, and (c) the closest UQRG GIM epoch with overlayed Swarm B along-track PSD section of ED at frequency 50 s. South America region on September 10, 2015. Black line (top-left and bottom) denotes the time of passing the latitude of mainshock.

A preliminary coincidence of Swarm B disturbance and TIDs in-ground GNSS data (Fig. 9) must be inspected by the scatter plot because the keogram indicates only a common distance of the observations from the epicenter and not necessarily a common location. The PSD from the spectrogram sampled at 50 s frequency is visualized along several min. of the Swarm B track, together with the same amount of ground GNSS L_G-F data in Figure 11. We can see clear TIDs located easterly from the Altiplano plate edge and exactly co-located with Swarm pass and the disturbance of changing frequency presented in the spectrogram in Figure 10c. The location of these TIDs parallel to the Altiplano edge and their speed similar to gravity waves makes their seismic origin fairly probable.

Fig. 11 Scatter map of around 16 min of band-pass filtered phase LG-F combination from ground GNSS stations of Chilean and Argentinian networks (black triangles), recalculated to TEC units, on September 10, 2015. Larger dots indicate Swarm B along-track power spectrum of LP ED disturbances sampled from the spectrogram at 50 s. Magenta square is the last Swarm position indicating the direction of flight. Blue lines are tectonic plate boundaries. Black rectangles indicate places of the most evident TID. Yellow triangle is the epicenter of Chile-Illapel EQ.

The closest Swarm track preceding the mainshock is the Swarm B track occurring on September 16, 2015 around 21:50 UTC, which is 1 h before the mainshock. The TIDs preceding the mainshock have again the speeds similar to gravity waves (Galvan et al., 2012; Jin et al., 2014; Occhipinti et al., 2013). Some of the TIDs, which are noticeable in the keogram before the mainshock (Fig. 12), coincide in time and distance from the epicenter with Swarm B passing across this region. The spectrogram of Swarm B in-situ measured ED (Fig. 13c) and accompanying UQRG TEC map (Fig. 13b) show ionospheric disturbance and TEC enhancement above the Altiplano minor plate, respectively. The period of the anomaly centered rather at 50 s, decreases at the end but also has some higher frequencies of lower power accompanying the most evident one. The scatter plot of L_G-F in Figure 14 presents some correlated disturbances slightly westerly and easterly from the Swarm B track. However, it should be pointed that the density of IPPs is not satisfactory in the area of Altiplano plate and these disturbances probably could also cross the Swarm footprint in case of more dense data. It is probable that two places of noticeable TIDs in Figure 14 surrounded by black rectangles are continuous concentric rings of perturbations, and their ionospheric response is here broken due to the small number and adverse directions of STEC measurements. Quite good and long continuity of these disturbances can, however, be found in the keogram in Figure 12.

Fig. 12 The same as in Figure 9, in the UTC evening on September 16, 2015. The vertical lines with label denote Chile-Illapel Mw = 8.3 mainshock and aftershocks. The blue-red strips in the right-bottom corner are TID directly related to the mainshock. The other blue-red strips denote other pre-seismic ionosphere disturbances.

Fig. 13 The same as in Figure 10. Swarm B, South America region on September 16, 2015, less than 1 h before the mainshock.

Fig. 14 The same as in Figure 11, on September 16, 2015, 1.1 h before the mainshock.

One hour later, when no Swarms are available, the M_w = 8.3 EQ has generated TIDs easily detectable in the band-pass filtered ground GNSS data. TIDs are mostly composed of circular waves centered at the EQ epicenter, and with velocity around 0.5 km/s are most probably due to acoustic waves induced by the Rayleigh waves, by dynamic coupling, and propagating upward in the atmosphere strongly amplified by the combined effects of the decrease of atmospheric density and the conservation of kinetic energy density (Occhipinti et al., 2018). The TIDs triggered by the mainshock are easily noticeable in the right-bottom corner of Figure 12. These disturbances, which have higher speeds than TIDs before the mainshock, can also be presented in the scatter plot as concentric traveling rings (Fig. 15). Several close aftershocks occurred just after the mainshock, indicated by red vertical lines in Figure 12.

Fig. 15 Scatter map of two close epochs of band-pass filtered phase LG-F combination from ground GNSS stations of Chilean and Argentinian networks, recalculated to TEC units, on September 16, 2015, 0.7 h after the mainshock. Blue lines are tectonic plate boundaries. Black outline marks places of the most evident TID.

The closest possible Swarm passes after the Chile-Illapel EQ of M_w = 8.3 are available in the morning of September 17, 2015, 8 h after the mainshock. The first pass of Swarm A over this region was at 7:22 AM UTC, which is 3:22 AM of Chilean time. The Sun was on the opposite side, and the TEC values were small at that time (Fig. 17b). Much lower Sun activity directly reflects in much weaker red-blue strips in the keogram (Fig. 16). TIDs can be found for the last time in Figure 16 around 7:00–8:00 AM UTC and are almost undetectable in the next hours. In Figure 16 we see again that not all station-PRN sounding directions can cross TIDs, and red-blue strips are often separated by undisturbed, green IPPs. Therefore, many authors select GNSS satellites in the keogram creation to have an effect of continuous TIDs, which is avoided in this work to better show the ionosphere behavior. The ED spectrogram and scatter plot are prepared only for the time of the first pass of Swarm A around 7:22 UTC (Fig. 16, leftmost track), which as the only one reaches some blue-red strips moderately visible in the keogram. The spectrogram of this pass (Fig. 17c) shows clear disturbances over the Altiplano plate and over the continent (Fig. 17b), which are of a smaller amplitude in relation to that observed during afternoon hours in the previous figures (e.g. Fig. 10).

Fig. 16 The same as in Figure 9, in the UTC morning on September 17, 2015. The blue-red strips denote different TIDs, which are less detectable in the morning.

Fig. 17 The same as in Figure 10. Swarm A, South America region on September 17, 2015.

POD TEC signals are expected to provide similar disturbances in relation to LP ED data, with some exceptions. The STEC measurements between Swarm A spacecraft and GNSS satellites point in spatially different directions. Therefore, in-situ LP detection of ionosphere perturbation can be absolutely not reflected in POD STEC measurement in case when its direction does not coincide with disturbing wave direction. However, there is always some chance, and it is worth analyzing Level 2 POD TEC product. STEC is recalculated to VTEC at IPP altitude 50 km above the Swarm A. In STFT analysis, we use VTEC with exactly the same parameters as for LP data analysis. The difference is only in the power spectrum maximum range in the plots, which is different due to different data types. The position of GPS IPPs with respect to Swarm A orbit can be found in Figure 18b, as respective locations of colored PSD sections and magenta lines. The POD TEC measurement from Swarm A corresponding best to LP ED signal is taken to PRN 16, which flies west of Swarm A (Fig. 18b).

Fig. 18 The same track as in Figure 17, (a) Swarm A – GPS 16 POD TEC high-pass filtered data, (b) the closest UQRG GIM epoch with overlayed Swarm A along-track PSD section of TEC at frequency 50 s, and (c) STFT spectrogram of residual POD TEC. Black line (a) and (c) denotes the time of passing the latitude of mainshock.

The scatter plot of L_G-F in the morning of September 17, 2015 finds TIDs directed to the ocean, which does not coincide exactly with the disturbances along with the Swarm A track (Fig. 19). These presumable TIDs are located quite close to and westerly from the epicenter of Chile-Illapel, over the ocean, which seems not to be related with the expected seasonal MSTIDs in the winter-spring transition due to its occurrence in the local night (Hernández-Pajares et al., 2012), but it does not exclude their relation with the tsunami. However, the gaps in IPP data significantly impede the exact recognition of their shapes. These disturbances, however, are away from the strong Sun activity, which additionally proves their seismic origin. The disturbance along the track of Swarm A in Figure 17, which can be recognized clearly with the use of STFT, can suggest that in-situ LP observation can be more sensitive to such anomalies than TEC calculated along the entire way between ground stations and GNSS satellites. The POD TEC measured along the direction LEO-GNSS (Fig. 18) also appears as more fortunate than ground measurements from Figure 19, however, more factors can play a role there, like e.g., reference altitude of measurements or data filtering band.

Fig. 19 The same as in Figure 11, in the morning of September 17, 2015.

3.1.4 Potential Chile 2015 tsunami signature in the directly estimated 4D electron density gradient

The combination of co-located ground-based and LEO-based GNSS measurements under the ionospheric tomographic model has been a powerful approach for studying the 4D ionospheric footprint of Space Weather at global scale (Hernández-Pajares et al., 1998). Indeed, the Tomographic Ionosphere model (TOMION) (Hernández-Pajares et al., 1998, 1999) is the baseline of the UPC VTEC GIMs produced for the International GNSS Service (IGS) since June 1st, 1998, which have been improved in comparison to their first release (Roma-Dollase et al., 2018). TOMION is directly based on the available measurements expressed in terms of voxel basis functions. The model can capture localized ionospheric electron content features, modeling properly the 4D components of the electron content gradients with the ionospheric-mapping-function. Mean ED of the voxels is obtained from the ionospheric geometry-free combination L_G-F of carrier phases only (see Hernández-Pajares et al., 2011). The solution is obtained thanks to the Kalman filter integrated into TOMION, treating the sun-fixed ED as a random walk and estimating the carrier phase bias (or ambiguity) as a random variable.

Very recently, an extension of the applicability of the TOMION model (multiTOMION) has been done by combining GNSS-POD LEOs-, GNSS ground- and vessel-, and DORIS-based measurements with the corresponding improvement of VTEC Global Ionospheric Maps (Hernández-Pajares et al., 2020). This extension can provide a direct estimation of 4D ED gradients as a good tool for future tsunami detection and characterization, augmented with the use of LEO constellations much larger than the Swarm one, such as e.g., Spire CubeSats. The larger constellations will assure higher co-location with potential future tsunamis. TOMION model is applied in this study to the calculation of 4D ED gradient directly from the voxel-based tomography, in the analysis of Chile-Illapel EQ and the following tsunami. We have run TOMION in an extremely sensitive mode: with high horizontal resolution (1° × 1° in right ascension and latitude), high temporal resolution (30 s), increased random walk noise, and the standard two-layer models suitable in ground-based tomographic modeling. Due to the lack of co-location of the Swarm constellation and other LEOs, we supported the model with the available multi-GNSS measurements from CSN, IGN, IGS, and other receivers available in this region. By applying this approach, we are able to find significant horizontal bottom ED gradients on September 16, 2015 after the EQ, following a radial pattern, which is compatible with the tsunami signature (Fig. 20b). Figure 20 a presents quite randomly distributed horizontal ED components during a similar time period of the previous day.

Fig. 20 Horizontal component of ED gradient corresponding to the bottom layer of voxels, the height of which ranges between 110 and 790 km and representing pierce points at mid-height of 450 km, (b) during day 2015-09-16: 16.5 min after EQ main shock, and (a) 1 day before the EQ at a similar time (2015-09-15, 23:11:00), computed with TOMION from the CSN and IGN networks of ground GNSS receivers.

3.2 Far away from the epicenter

3.2.1 Validation of Swarm co-seismic ionosphere disturbances (CID) detection with ionosonde data

The purpose of this study part is to examine Digisonde and Swarm coincident observations and consider the possible connection to EQ-triggered perturbations. The investigation of simultaneous observations in the bottom side and topside ionosphere can lead to results regarding the connection between medium-scale TIDs, spread F, and topside ED gradients with the effects triggered by the EQ at various ionospheric heights. In order to achieve this goal, the first step is to ingest the Digisonde-bottom side and the Swarm in situ-topside observation in a continuous ED distribution from the height of 90 km up to the altitude of the Swarm satellite. To that effect, Swarm-Digisonde coincidences were identified in order to have simultaneous and co-located observations. A Swarm-Digisonde coincidence is defined as the satellite pass above an area around the Digisonde with a radius of 2.5 degrees (≈250 km). Each Swarm pass above the area of coincidence has a maximum duration of 90 s. For each coincidence, the ED profile is calculated from the bottom side ionosphere up to Swarm A altitude by applying the α-C-SD model (Belehaki et al., 2021). The model assumes that above the hmF2 height and up to the Swarm altitude, the distribution of the ionospheric ED can be approximated with an α-Chapman function with variable scale height. The maximum ED at hmF2 measured by the Digisonde and the ED provided by Swarm LP are the two anchor points that define the shape of the ED distribution. The scale height at the hmF2 is calculated from the Digisonde bottom side profile (Huang & Reinisch, 1996), and it is provided as an ionospheric characteristic in the SAO files (Gamache et al., 1996). The scale height at the Swarm height is calculated with the α-Chapman function so that the boundary condition defined by the LP ED measurement at Swarm altitude is satisfied. The results are validated by comparing the reconstructed TEC with the GNSS TEC. In the vast majority of cases, the reconstructed total TEC matches with extremely high accuracy the GNSS TEC, providing strong evidence that the assumption for the α-Chapman ED approximation in the topside part between the hmF2 and the Swarm altitude is a valid one (Belehaki et al., 2021).

In this study, the α-C-SD model is applied to reconstruct the ED profile during Swarm A passes over Digisondes in a wide region around the EQ epicenter. The variability of the ionospheric characteristics in the bottom side and topside ionosphere is considered as an indicator for lower and upper atmosphere interactions. Based on the abovementioned criteria, there are no coincident Swarm data with Cachoeira Paulista station. Therefore, the results are presented for two Swarm A passes over the closest Boa Vista Digisonde, several hours after the main EQ shock. Swarm A passes over the Digisonde occurred 20 and 32 h after the mainshock. The location of the two Digisondes in respect to the EQ epicenter is given in Figure 21.

Fig. 21 Location of Boa Vista (BVJ03) and Cachoeira Paulista (CAJ2M) Digisondes in respect to the EQ epicenter.

Figure 22 presents Boa Vista Digisonde recordings from September 16, 2015 to September 18, 2015 as isoionic contour plots and as time-plots of the following ionospheric characteristics:

• The F2 layer critical frequency, foF2

• The peak height of the F2 layer, hmF2

• The minimum virtual height of the F layer, h′F

• The range spread of F-layer, QF

• The frequency spread between fxF2 and fxI, FF

Fig. 22 Data from Boa Vista Digisonde ionograms: (a) The isoionic contour plots, (b) the critical characteristics foF2 (blue line), h′F (black line), hmF2 (red line), and (c) the spread F parameters FF (blue line) and QF (red line).

All the ionograms used to extract the ionospheric characteristics were manually scaled for higher accuracy. The cadence of the vertical incidence soundings in Boa Vista Digisonde is 5 min. The vertical solid line in Figure 22 corresponds to the occurrence time of the main EQ shock, while the vertical dashed lines correspond to the Swarm A passes over the Boa Vista Digisonde coincidence area.

Shortly after the EQ, the bottom side ionosphere shows signatures of irregularities since spread F (FF and QF) is recorded in the ionograms consequently for at least 8 h, from 23:00 UTC on September 16, 2015 to 09:00 UTC on September 17, 2015. During this interval, signatures of medium-scaled TIDs (MSTIDs) are observed in the isoionic contours (Fig. 22a). The next intensification of irregularities appears from around 20:00 UTC on September 17, 2015 to around 06:00 UTC on September 18, 2015, manifested as high-frequency fluctuations in the isoionic contours, as strong fluctuations in h′F and hmF2 and high values of QF and FF characteristics.

The data collected during the first Swarm A pass over Boa Vista Digisonde on September 17, 2015 from 19:11:20 UTC to 19:12:37 UTC are shown in Figure 23a with the following order from top to bottom: the scale height at Swarm A altitude calculated with the α-C-SD model, the in-situ ED from Swarm LP experiment and the Swarm A irregularity indices of gradient ED at 20 km, 50 km, and 100 km (Jin et al., 2019). The Digisonde ED profile (SAO profile) and the topside ED profile reconstructed with the α-C-SD model in a red color trace (Swarm profile) are shown in Figure 23b together with the in-situ Swarm ED data in green color, observed during the Swarm A pass over the Digisonde. It should be noted here that the Digisonde ED profile is calculated from ionogram inversion only up to the hmF2 height. Its topside part is the result of the Huang–Reinisch model (Huang & Reinisch, 2001), which assumes a constant scale height. However, it should be taken into account that this is only an approximation and it is shown here as a reference. During this short duration of Swarm A pass over the Digisonde, the difference between the SAO profile and the Swarm reconstructed profile is very pronounced. This is the result of the different scale heights. The SAO profile reconstruction assumes that the scale height in the topside is constant and equal to the scale height at hmF2. The Swarm profile reconstruction assumes that the scale height varies so that the electron density calculated with the α-Chapman approximation becomes equal to the measured electron density from Swarm LP. The SAO scale height at hmF2 is 104.20 km, whereas the scale heights calculated with the Swarm reconstruction α-C-SD model get values in the range of (25.18 ± 8.52) km during the 90 s of Swarm pass over the Digisonde. This large decrease of the topside scale height indicates that, at least locally, the electron density above the hmF2 and up to ~450 km decreases fast with altitude. It must be noted that this local decrease in the topside ED is observed simultaneously with fluctuations in the virtual scale height of the F layer and in the isoionic contour plots, which might be indicative of MSTIDs in the bottom side ionosphere.

Fig. 23 Swarm and Digisonde observations on September 17, 2015 (a) the topside scale-height at Swarm altitude calculated with the α-C-SD model, the Swarm electron density Ne and the Swarm irregularity indices gradient ED at 20 km, 50 km and 100 km during the 90 s of the satellite pass over the Digisonde from 19:11:19 UTC to 19:12:37 UTC. (b) The Digisonde electron density profile and the reconstructed topside electron density with the α-C-SD model in a red color trace; the green color trace represents the insitu Swarm ED observed during the Swarm pass over the Digisonde.

During the 90 s of coincidence, the topside Grad Ne index at the range of 20 km shows systematic fluctuations, which faint at larger ranges, as seen in the topside Grad Ne index at 100 km. These observations indicate upward propagation of plasma irregularities which break into small structures in the topside ionosphere (Kuo et al., 1998; Eltrass et al., 2014, 2016). These waves may be initiated by the numerous aftershocks that occurred for the whole day after the mainshock (as presented in Fig. 1), produced by the conversion of the crust oscillations into acoustic waves and then in variations in the ionospheric electron density.

The second Swarm A pass over Boa Vista Digisonde occurred 32 h after the main EQ shock, from 0655:11 UTC to 0656:28 UTC on September 18, 2015. As concluded from the inspection of Figure 22, this second pass occurred when the ionosphere was relatively unperturbed since the QF and FF indices were close to zero and the ionospheric characteristics h′F and hmF2 did not show any fluctuations. The top panel (Fig. 24a) provides the time plots of the characteristics and calculated parameters during the 90 s. of Swarm A pass over the Digisonde. The lower panel (Fig. 24b) presents the ED profile from the Digisonde (SAO profile) and from the α-C-SD model (Swarm profile) together with the insitu ED data from Swarm LP, as in Figure 23. Some very weak disturbances are recorded during the Swarm A pass over the Digisonde, which are also reflected in the scale height.

Fig. 24 Same as Figure 23 for the satellite pass over Bua Vista Digisonde which occurred from 06:55:11 UTC to 06:56:28 UTC on September 18, 2015.

For this second Swarm A pass over the Boa Vista Digisonde, the reconstructed topside ED profile agrees relatively well with the Digisonde reconstructed ED profile. The SAO scale height at hmF2 is 38.0 km, whereas the scale heights calculated with the Swarm reconstruction model get values in the range of (41.38 ± 7.80) km during the 90 s of Swarm pass over the Digisonde. In other words, the SAO scale height is within the range of values of the topside scale height. This agreement indicates that the ionospheric density in the layer above hmF2 and up to ~450 km is relatively stable. The ionospheric conditions during this second pass are very different in the first Swarm pass over Boa Vista, where the reconstructed topside ED profile differs significantly from the Digisonde reconstructed ED profile, providing an indication that the proposed method is efficient for the detection of irregularities in the topside ionosphere.

3.2.2 TIDs detection on board of GRACE and Swarm in the ocean area

The GRACE mission is equipped with Ka-Band and K-Band ranging systems (KBR). The system of KBR tracking measures the distance between the two GRACE satellites (satellite-to-satellite tracking – SST) with very high precision. Each of the two satellites transmits a carrier signal and measures the phase of the carrier generated by the other satellite relative to the signal it is transmitting. The sum of the phases generated is proportional to the range change between the satellites, while the phase variation due to long-term instability in each clock cancels out. The KBR system can measure the range (with a bias) to the μm level. However, the observations suffer from the ionospheric delay caused by the electrons between the transmitting satellite T and the receiving satellite R. TIDs which are induced by tsunamis cause changes in ED which can be observed by the GRACE mission, in case when the satellite passes these TIDs. The along-track measurements of selected observation techniques detecting the gravity wave induced by EQ or tsunami and, thus, variations due to TIDs, can be analyzed by means of continuous wavelet transform (DWT) (Schmidt, 2001). In order to determine whether these signal variations are present in the observations from Swarm and GRACE at the time when they pass the tsunami region far from the epicenter, a two-step procedure is applied:

1. In a first step, a tsunami propagation model is used to determine the approximate time when the satellites crossed the tsunami in the Pacific Ocean.

2. In a second step, the LP observations measured during the tsunami crossing are analyzed for signal variations with periods of 15–30 s. For this purpose, a wavelet transform is used.

There are co-locations of the Swarm and GRACE satellites over the region where the tsunami was spreading 17 h after the mainshock, in the time between 15:15 and 15:50 UTC.

Using the wavelet transform the ED observations of the LP sensors on Swarm A, and C satellites, as well as the K/Ka-Band observations from GRACE are analyzed for the time periods shown as co-located in Figure 25. As the aim of this study is to detect the TID’s in the observations of ED, the change of the signal as the centered differences

$$\mathrm{f}\left(\mathrm{t}\right)=\frac{1}{2\mathrm{h}}\left(\mathrm{y}\left(\mathrm{t}-∆\mathrm{t}\right)\mathrm{ }-\mathrm{y}\left(\mathrm{t}+∆\mathrm{t}\right)\right)$$(1)

of the observations are computed at first in order to eliminate any biases from the signal. In equation (1) y(t) is the observed ED at the time t and ∆t is the temporal sampling of the observed data; in case of the Swarm satellites A, B and C the sampling interval is fixed to ∆t = 1 s for the LP observations. The KBR observations of GRACE are initially sampled with a temporal distance of ∆t = 5 s. Figures 26a–26c) shows the centered differences of Swarm A and Swarm C LP data and the GRACE KBR instrument. The scalograms calculated for LP ED observations from the satellites Swarm A and C co-located with the tsunami between 15:15 UTC and 15:35 UTC on September 17, 2015 are presented in Figures 26d–26e. As it can be seen from the scalograms, signal variations with significant periods around 15 s are detectable, which may be related to the TIDs induced by the tsunami (Yang et al., 2019). Since the Swarm satellites A and C follow approximately one after the other, with a small shift in the ground tracks (cf. Fig. 25), the two scalograms show similar characteristics. Differences in the scalograms may result from their distance in the orbit and thus, from a shift of the phase of the TID, as well as from slightly different angles at which the satellite crosses the TID. The GRACE satellites cross the tsunami wave some min after the two Swarm satellites, i.e., between 15:30 UTC and 15:50 UTC. Figure 26f shows the scalogram of the differenced GRACE K/Ka-band observations, where signal variations with periods of 10 s and 30 s show significant energy distribution, similar to those of the two Swarm satellites. The periods in the range of 15 sec appear in all scalograms at minute 33, as well as the periods in the range of 25 s, which appear in Figures 26d–26e at the same time and in 26f with a time offset, should be emphasized at this point. The differences visible in the scalograms of the different satellites may occur due to the different altitudes and may be caused by different times and directions of propagation.

Fig. 25 Groundtrack of the Swarm A and C satellites plotted with blue and red color, respectively, for the time between 15:15 UT and 15:35 UT on September 17, 2015 and the groundtrack of the GRACE satellites plotted in green color for the time span between 15:30 UT and 15:50 UT on September 17, 2015.

Fig. 26 Scalograms of the differenced LP observations of Swarm A and C and the differenced Ka-band observations of GRACE for a selected time interval of Tsunami crossings.

4 Summary and conclusions

The objective of this work was to demonstrate the synergistic effect of a wider scope of data sources and techniques worth considering together, as well as a wider spatiotemporal window analyzing seismically triggered changes in the ionosphere and changes correlated with tsunami. One general conclusion referring to the data is that LEO satellite data have a notably wider spatial range with respect to ground GNSS, and this advantage should be exploited. The second observation related to the data is that the more data, the lower uncertainty and wider the scope of different complementary observations. The general conclusion referring to the techniques is that filtering and spectral analysis play crucial role with respect to ground and especially LEO data, as the ionospheric signals are composed of many externally (Sun) and internally induced (Earth) anomalies. The general conclusion referring to spatio-temporal issues is that focusing only on the local dense GNSS network and close vicinity of mainshock significantly limits our knowledge about composed coupling processes between the lithosphere and ionosphere. The very detailed analysis of shape and velocities of TIDs after Tohoku EQ/tsunami in 2011, based on 1200 Japanese GEONET stations, cannot be repeated with the same precision here, as the number of used Chilean and Argentinian stations is around 200 and the analyzed region is larger. Nevertheless, the extension of the observation time and data scope gives us a very interesting preliminary inspection of the most promising timeframes and directions of future research.

Ground GNSS data have shown many different responses of seismic activity to the ionosphere, including ionospheric anomalies before the mainshock that have speeds around 0.3 km/s, clear and faster (around 0.5–0.6 km/s) response to the mainshock (M_w = 8.3), as well as post-seismic disturbances, presumably also related with the tsunami. The traveling character of these disturbances and their relation with seismicity are evident in the keograms in Section 3.1.3. Swarm along-track data also reveals ionospheric anomalies of a small period (<500 km along the track), and their sizes and specific spectral patterns can be found by spectral analysis. The inspection of ground GNSS data supplies a very reliable validation for the disturbed Swarm signals, as selected Swarm disturbances are spatiotemporally consistent with TIDs clearly seen in ground GNSS data (e.g. Fig. 11). Therefore the spectral patterns of Swarm disturbances in ED and POD TEC are especially worth further investigation and comparison with ground GNSS for different EQ/tsunami events.

The combined analysis of coincident observations recorded over the region affected by Chile-Illapel EQ, from ground-based ionosonde sounders and the Swarm LP instrument was carried out in Section 3.2.1 in order to understand the lower and upper atmosphere coupling process. Disturbances detected with ground-based ionograms were classified as medium scale TIDs, frequency spread-F, and range spread F. Swarm A detects coherent disturbances as highly variable gradients in the ED and highly variable scale height. This observation provides evidence of ED irregularities extended from the bottom side to the lower altitudes of the topside ionosphere. Swarm A detects coherent disturbances as variable gradients in the ED and variable scale height. This observation indicates upward propagation of plasma irregularities which break into small structures in the topside ionosphere, a mechanism for the cascading of ionospheric irregularities proposed by Kuo et al. (1998), Eltrass et al. (2014, 2016). The results indicate that although Swarm LP can detect in situ plasma perturbations, the simultaneous analysis of Digisonde and Swarm coincident data provide an improved nowcasting information for plasma disturbances triggered from the bottom side ionosphere and consequently an improved characterization of the whole ionosphere – lower atmosphere – lithosphere dynamics.

Simultaneous observation of the tsunami-induced TIDs on different satellites and the study of the horizontal coupling processes can determine the direction of propagation of the tsunami and its impact on the coastline. Section 3.2.2 describes the application of the wavelet transform, and similar signatures with short period lengths detected in the K/Ka-band observations of GRACE and in Swarm LP observations (altitude 450 km), in time of Chile-Illapel tsunami propagation. The satellites of GRACE and Swarm missions flew successively over the tsunami region. The analysis of in-situ observations from the satellites can be used to determine the position of the tsunami as it propagates across the ocean. The roadmap studies provided the first opportunity to observe the tsunami-induced gravity waves in the ionosphere, however, their periods need more further work to be better confirmed.

The velocities of the detected seismically triggered TIDs have been analyzed in Section 3.1.2 by the ADDTID method. The multi-TID blind analysis by ADDTID is done from around 1 h before the mainshock to 1 h after and confirms the predominant concentration of seismic TIDs in the north-east direction from the epicenter. This is consistent with keograms and scatters representations of TIDs in Section 3.1.3. Nevertheless, there is also an evident increase of TID occurrence westerly from the epicenter, as detected by ADDTID and in the scatter plots of L_G-F, and some of them can be related to the tsunami, due not only to counter direction regarding the epicenter but also to the correlation of the increased velocity with increased sea depth. EQ/tsunami compatible ionospheric disturbances located westerly from the epicenter, over the ocean, have also been detected in the 4D ED gradient, directly based on the voxel-based tomography model, as summarized in Section 3.1.4.

Dense ground GNSS network availability is limited to a few regions on the globe. Swarm and other LEO satellite data, in turn, are globally distributed. Additionally, due to the recent developments, the usage of GNSS receivers will be possible very soon, on board tens or even hundreds of CubeSats developed for ionospheric sounding. Therefore, the detection of TIDs from the POD GNSS receiver onboard the LEOs and CubeSats, or from the other instruments (LP, KBR) can be especially useful for the detection of ionospheric responses to EQs and tsunami. The spectral decomposition presented in Sections 3.1.3 and 3.2.2 can help in the analysis of short along-track signals. The potential warning obtained from the LEOs can then be further analyzed in a second prompt round of EQ/tsunami confirmation and characterization. This can be done using a combination of ground GNSS data, ionosondes, and 4D ED gradient techniques. Therefore, current preliminary studies can be helpful in the creation of the system of wide-area detection of seismic ionospheric anomalies for prompt detection of tsunamis.

Acknowledgments

The study was funded from COSTO project, which has received funding from the ESA contract 4000126730/19/NL/IA, in the framework of the “Swarm+ Coupling: High-Low Atmosphere Interactions” Call. Anna Belehaki acknowledges financial support provided by the AFRL grant award FA9550-19-1-7019. The editor thanks Dimitar Ouzounov and Sergey Pulinets for their assistance in evaluating this paper.

References

All Figures

	Fig. 1 Seismic activity from magnitude 5.0 in region of South America (close to Chile-Illapel) during September 2015. (a) Map of EQs in September 2015 in the neighborhood of Chile-Illapel. (b) Plot of EQs with respect to mainshock. Magenta stems denote EQ magnitude multiplied by 10, brown stems denote hypocenter depth. Source of subfigure (a) and data: USGS.
In the text

	Fig. 2 Residuals from TG and DART readings of Chile Illapel 2015 tsunami (September 16th, 2015) recorded 15.5 h after source EQ.
In the text

	Fig. 3 Travel Time Map of the Chilean 2015 tsunami (September 16th, 2015). Source: NOAA, U.S.TWS: https://www.tsunami.gov/previous.events/?p=09-16-2015-Chile(M8.3).
In the text

Fig. 4 From top to bottom: the Dst index, and the relative deviation of TEC and foF2 values with respect to monthly medians (black lines) over Cachoeira Paulista (Brazil) during the whole month of September 2015. The color lines represent the relative standard deviation (+/−) of the values taken into account for the estimation of the medians in each case. The latter is included as an indicator of the ionospheric variability in normal conditions that are used as a threshold for identifying significant ionospheric disturbances in a large scale perspective.

In the text

Fig. 5 GNSS receiver distribution in Chile and Argentina. Cyan crosses denote the location of the GNSS receivers, and the red star denotes the epicenter of the Mw 8.3 EQ, on September 16, 2015. The grid cells (a)–(d) denote four subnetworks with a geographical grid of 10° × 10°. The ocean depth is shown by the cold color map, originating from the Earth observatory of the US NASA blue marble with topography and bathymetry.

In the text

Fig. 6 (a) IPP locations at the height of 200 km, corresponding to the GNSS receivers, from 22:45 to 23:15 UTC, the color codes of cycles for the IPP locations of different GPS satellites, the red star for the location of epicenter, (b) time evolution of the arbitrarily aligned slant TECs, for the GPS satellite of PRN 12, at the time of 22:30–24:00 UTC, the color codes of lines correspond to different GNSS receivers, the black dash line indicates the time of the EQ at 22:54 UTC.

In the text

	Fig. 7 Detected MSTIDs of each subnetwork. Parameter evolution of detected TIDs before the EQ during 22:00–24:00 UTC, representing non-uniformly scaled polar plots for TIDs of the corresponding grid networks (a)–(d), the location representing for the azimuth vs. velocity of 0–360 degree and 0–900 m/s, the circle size for wavelength of 50–600 km, and the color code for different time intervals.
In the text

	Fig. 8 Time evolution of TID velocity distribution during 22:00–24:00 UTC, for the subnetwork (a), (b), (c), and (d), respectively. The scale of the velocity axis ranges from 0 m/s to 600 m/s.
In the text

Fig. 9 Keogram of bandpass filtered LG-F phase combination from ground GNSS stations of Chilean and Argentinian networks, recalculated to TEC units, in the UTC evening on September 10, 2015. All distances of IPP tracks are calculated with respect to the position of the mainshock. Large dots denote Swarm A and B passes with distances respective to the mainshock location, and their color indicate size of power PSD sampled at 50 s. from the spectrogram. The blue-red strips denote different TIDs.

In the text

	Fig. 10 (a) Swarm B LP high-pass filtered ED data, (b) STFT spectrogram of residual ED, and (c) the closest UQRG GIM epoch with overlayed Swarm B along-track PSD section of ED at frequency 50 s. South America region on September 10, 2015. Black line (top-left and bottom) denotes the time of passing the latitude of mainshock.
In the text

Fig. 11 Scatter map of around 16 min of band-pass filtered phase LG-F combination from ground GNSS stations of Chilean and Argentinian networks (black triangles), recalculated to TEC units, on September 10, 2015. Larger dots indicate Swarm B along-track power spectrum of LP ED disturbances sampled from the spectrogram at 50 s. Magenta square is the last Swarm position indicating the direction of flight. Blue lines are tectonic plate boundaries. Black rectangles indicate places of the most evident TID. Yellow triangle is the epicenter of Chile-Illapel EQ.

In the text

	Fig. 12 The same as in Figure 9, in the UTC evening on September 16, 2015. The vertical lines with label denote Chile-Illapel Mw = 8.3 mainshock and aftershocks. The blue-red strips in the right-bottom corner are TID directly related to the mainshock. The other blue-red strips denote other pre-seismic ionosphere disturbances.
In the text

	Fig. 13 The same as in Figure 10. Swarm B, South America region on September 16, 2015, less than 1 h before the mainshock.
In the text

	Fig. 14 The same as in Figure 11, on September 16, 2015, 1.1 h before the mainshock.
In the text

	Fig. 15 Scatter map of two close epochs of band-pass filtered phase LG-F combination from ground GNSS stations of Chilean and Argentinian networks, recalculated to TEC units, on September 16, 2015, 0.7 h after the mainshock. Blue lines are tectonic plate boundaries. Black outline marks places of the most evident TID.
In the text

	Fig. 16 The same as in Figure 9, in the UTC morning on September 17, 2015. The blue-red strips denote different TIDs, which are less detectable in the morning.
In the text

	Fig. 17 The same as in Figure 10. Swarm A, South America region on September 17, 2015.
In the text

	Fig. 18 The same track as in Figure 17, (a) Swarm A – GPS 16 POD TEC high-pass filtered data, (b) the closest UQRG GIM epoch with overlayed Swarm A along-track PSD section of TEC at frequency 50 s, and (c) STFT spectrogram of residual POD TEC. Black line (a) and (c) denotes the time of passing the latitude of mainshock.
In the text

	Fig. 19 The same as in Figure 11, in the morning of September 17, 2015.
In the text

	Fig. 20 Horizontal component of ED gradient corresponding to the bottom layer of voxels, the height of which ranges between 110 and 790 km and representing pierce points at mid-height of 450 km, (b) during day 2015-09-16: 16.5 min after EQ main shock, and (a) 1 day before the EQ at a similar time (2015-09-15, 23:11:00), computed with TOMION from the CSN and IGN networks of ground GNSS receivers.
In the text

	Fig. 21 Location of Boa Vista (BVJ03) and Cachoeira Paulista (CAJ2M) Digisondes in respect to the EQ epicenter.
In the text

	Fig. 22 Data from Boa Vista Digisonde ionograms: (a) The isoionic contour plots, (b) the critical characteristics foF2 (blue line), h′F (black line), hmF2 (red line), and (c) the spread F parameters FF (blue line) and QF (red line).
In the text

Fig. 23 Swarm and Digisonde observations on September 17, 2015 (a) the topside scale-height at Swarm altitude calculated with the α-C-SD model, the Swarm electron density Ne and the Swarm irregularity indices gradient ED at 20 km, 50 km and 100 km during the 90 s of the satellite pass over the Digisonde from 19:11:19 UTC to 19:12:37 UTC. (b) The Digisonde electron density profile and the reconstructed topside electron density with the α-C-SD model in a red color trace; the green color trace represents the insitu Swarm ED observed during the Swarm pass over the Digisonde.

In the text

	Fig. 24 Same as Figure 23 for the satellite pass over Bua Vista Digisonde which occurred from 06:55:11 UTC to 06:56:28 UTC on September 18, 2015.
In the text

	Fig. 25 Groundtrack of the Swarm A and C satellites plotted with blue and red color, respectively, for the time between 15:15 UT and 15:35 UT on September 17, 2015 and the groundtrack of the GRACE satellites plotted in green color for the time span between 15:30 UT and 15:50 UT on September 17, 2015.
In the text

	Fig. 26 Scalograms of the differenced LP observations of Swarm A and C and the differenced Ka-band observations of GRACE for a selected time interval of Tsunami crossings.
In the text