The correlation between ionospheric electron density variations derived from swarm satellite observations and seismic activity at the Australian–Pacific tectonic plate boundary

Remote Sens. 2023, 15, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing Article The Correlation between Ionospheric Electron Density Variations Derived from Swarm Satellite Observations and Seismic Activity at the Australian–Pacific Tectonic Plate Boundary Wojciech Jarmołowski 1,*, Paweł Wielgosz 1, Manuel Hernández-Pajares 2, Heng Yang 2, Beata Milanowska 1, Anna Krypiak-Gregorczyk 1, Enric Monte-Moreno 2, Alberto García-Rigo 2, Victoria Graffigna 2 and Roger Haagmans 3 1 Faculty of Geoengineering, University of Warmia and Mazury in Olsztyn, ul. Oczapowskiego 2, 10-719 Olsztyn, Poland; pawel.wielgo[email protected].pl (P.W.); beata.milanowsk[email protected].pl (B.M.); a.krypiak-gregor[email protected].pl (A.K.-G.) 2 Department of Mathematics, IonSAT, Universitat Politecnica de Catalunya, 08034 Barcelona, Spain; [email protected] (M.H.-P.); [email protected] (H.Y.); enric.mo[email protected] (E.M.-M.); [email protected] (A.G.-R.); [email protected] (V.G.) 3 ESTEC, European Space Agency, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands; [email protected] * Correspondence: [email protected] Abstract: Swarm electron density (Ne) observations from the Langmuir probe (LP) can detect ionospheric disturbances at the altitude of a satellite. Along-track satellite observations provide a large number of very short observations of different places in the ionosphere, where Ne is disturbed. Moreover, different perturbations occupy various Ne signal frequencies. Therefore, such short signals are more recognizable in two dimensions, where aside from their change in time, we can observe their diversity in the frequency domain. Spectral analysis is an essential tool applied here, as it enables signal decomposition and the recognition of composite patterns of Ne disturbances that occupy different frequencies. This study shows a high-resolution application of short-term Fourier transform (STFT) to Swarm Ne observations in the Papua New Guinea region in the vicinity of earthquakes, tsunamis, and related general seismic activity. The system of tectonic plate junctions, including the Pacific–Australian boundary, is located orthogonally to Swarm track footprints. The selected wavelengths of seismically induced ionospheric disturbances detected via Swarm are compared with the three sets of three-month records of seismic activity: in the winter solstice of 2016/2017, when seismic activity was highest, and in the summer solstice and vernal equinox of 2016, which were calmer. Moreover, more Swarm data records are analyzed at the same latitudes for validation purposes, in a place where there are no tectonic plate boundaries that are orthogonal to the Swarm orbital footprint. Additional validation is supplied through Swarm Ne observations from completely different latitudes, where the Swarm orbital footprint orthogonally crosses a different subducting plate boundary. Aside from the seismic energy, the solar radio flux (F10.7), equatorial plasma bubbles (EPBs), and geomagnetic ap and Dst indices are also reviewed here. Their influence on the ionospheric Ne is also found in Swarm observations. Finally, the Pearson correlation coefficient (PCC), applied to the pairs of 3-month time series created from Swarm Ne variations, seismic energy, ap, Dst, and F10.7, summarizes the graphical inspection of mutual correlations. It points to the predominant correlation of Swarm Ne disturbances with seismicity, especially during nighttime. We show that most of the Ne disturbances at a selected wavelength of 300 km correlate more with seismicity than with geomagnetic and solar indices. Therefore, Swarm LP can be assessed as being capable of observing the lithosphere–atmosphere–ionosphere coupling (LAIC) from the orbit. Citation: Jarmołowski, W.; Wielgosz, P.; Hernández-Pajares, M.; Yang, H.; Milanowska, B.; Krypiak-Gregorczyk, A.; Monte-Moreno, E.; García-Rigo, A.; Graffigna, V.; Haagmans, R. The Correlation Between Ionospheric Electron Density Variations Derived from Swarm Satellite Observations and Seismic Activity at the Australian–Pacific Tectonic Plate Boundary. Remote Sens. 2023, 15, x. https://doi.org/10.3390/xxxxx Academic Editors: Xuemin Zhang, Chieh-Hung Chen, Yongxin Gao and Katsumi Hattori Received: 12 October 2023 Revised: 11 November 2023 Accepted: 27 November 2023 Published: date Copyright: © 2023 by the authors. Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Remote Sens. 2023, 15, x FOR PEER REVIEW 2 of 37 Keywords: swarm satellites; electron density; seismically induced ionospheric disturbance; earthquake; power spectral density 1. Introduction Studies of seismically induced ionospheric disturbances (SIIDs) are typically limited to selected large earthquakes or to the list of events that have a significant magnitude and directly trigger local SIIDs. However, tectonic processes are very complex and related to various types of motions of the Earth’s lithosphere, as well as to the generation of electric charge and chemical reactions [1–3]. The accumulated energy is released from the lithosphere and generates acoustic waves and atmospheric gravity waves that may travel to the atmosphere through weak zones in the Earth’s crust that are located at tectonic plate boundaries [4]. These waves cause disturbances in the temperature, pressure, and electromagnetic field [5,6]. The generated atmospheric gravity waves propagate right up to the ionosphere. When these waves reach the ionosphere, they cause detectable changes in the electron density (Ne) and total electron content (TEC). A strong electric field can also be induced prior to an earthquake, and this is characterized by radioactive particle emanation (mainly radon) and increased conductivity, causing electric field anomalies [2,7]. There is also an electrostatic channel of seismo-ionospheric coupling described by Freund [8] and Hayakawa et al. [7], in which positive ionization holes are generated from the stressed ground. According to their origin, the ionospheric disturbances can be divided into two categories. The first one refers to the Sun’s activity, e.g., the solar flare eruption or coronal mass ejection, while the second one includes anthropogenic and natural events on the Earth. Both categories have different scales and duration times. The disturbances from Earth-based sources last for several hours or days, and their scale is about hundreds or thousands of kilometers [9,10]. Importantly, the ionospheric anomalies caused by earthquakes or volcano eruptions arise from a process called lithosphere–atmosphere–ionosphere coupling (LAIC) [11–15]. Due to several tens of satellites and thousands of stations, ground-based Global Navigation Satellite System (GNSS) data are the most promising for research on SIIDs. This fact comes from one main reason: ground-based GNSS observations provide the only reliable spatial view of SIIDs and enable recognition of their spatial shapes and movement velocity. A vast majority of undertaken research related to SIID analyses based on groundbased GNSS data was focused on the identification of different kinds of ionospheric waves and their propagation after large earthquake events [9,13,16–21]. Some of the authors apply statistical analyses of TEC responses to the earthquakes over a longer time span, which are typically a basis for studies on earthquake precursors and preseismic disturbances [22,23]. The application of ground-based GNSS data is focused on the locations where the stations provide spatially dense observations. There is, however, a limited number of dense GNSS networks. The ground GNSS stations have a heterogeneous distribution worldwide. Moreover, there are also many places along the subducting plate boundaries where GNSS stations are very sparse or unavailable. A low density of GNSS stations significantly impedes the identification of SIIDs in many interesting places and limits the study of them on a global scale. However, the observed increasing number of explorer satellite missions dedicated to global ionosphere observation can significantly extend the global approach. Low Earth orbiting (LEO) satellites can visit all the regions of the Earth, which is interesting for the global investigation of LAIC. An example of a completed satellite mission that was designed to monitor the variations in the atmosphere related to earthquakes is the French DEMETER (Detection of Electromagnetic Emissions Transmitted from Earthquake Regions) mission, active from 2004 to 2010. The DEMETER satellite had been orbiting at a Sun-synchronous circular orbit, having an approximate altitude of 710 km and inclination of 98.2°. The scientific Remote Sens. 2023, 15, x FOR PEER REVIEW 3 of 37 payload on the DEMETER satellite included a Langmuir probe (LP) measuring Ne and electron temperature. DEMETER provided observations of electromagnetic fluctuations related to earthquakes [24,25]. DEMETER data have also been examined in various ways to define the perturbations of the ionosphere caused by acoustic gravity waves from earthquakes [24–27]. The ionospheric disturbances associated with preseismic activity were also investigated in a study based on DEMETER data, conducted by Zhang et al. [28], and studies combining GNSS and DEMETER data conducted by Akhoondzadeh et al. [29] and Ibanga et al. [30]. The existing satellite missions that include LPs amongst the scientific payloads are the first European Space Agency (ESA) constellation mission for Earth observation called Swarm [31,32] and the Chinese–Italian space mission named the China Seismo-Electromagnetic Satellite (CSES) [33]. The Swarm constellation comprises three satellites, which compose the fifth Earth Explorer mission that was approved in ESA’s Living Planet Program, and they were successfully launched on 22 November 2013. The orbits of Swarm Alpha (Swarm A) and Charlie (Swarm C) have an initial altitude of 460 km and inclination of 87.4°, whereas Swarm Bravo (Swarm B) runs at an altitude of 530 km and at a 88° inclination angle. The research objective of the Swarm mission is to provide the best-ever survey of the geomagnetic field and electric field in the atmosphere using precise magnetometers and electric field instruments. Very frequently, Swarm data studies related to earthquakes incorporate magnetic data from all three satellites. There is published evidence of precursory signals in the magnetic data [34,35], as well as works referring to co-seismic signals in the magnetic field [36–38]. Zhu et al. [35] referred to the extended number of earthquakes, including ones of smaller magnitudes, and showed one example of spectral recognition of disturbed signals, which is also an objective of our current study. Ghamry et al. [39] analyzed the seismic precursory magnetic anomalies measured onboard Swarm together with different additional data, like thermal anomalies from other explorer satellites, like Terra and Aqua. Aside from the magnetic sensors, the Swarm payload also includes LPs and GNSS receivers for precise orbit determination, which can be used to observe in situ ionospheric Ne and topside TEC, respectively. Stanica et al. [40] jointly studied the variations in magnetic and electric field parameters as precursory to some moderate-magnitude earthquakes in Romania. A stronger event in Mexico was investigated by Marchetti and Akhoondzadeh [41] with the application of magnetic and Ne data, whereas Marchetti et al. [42] have studied seismic sequences in Italy using a variety of Swarm observations. In terms of earthquake number, an extended study was conducted by De Santis et al. [43], who prepared a statistical assessment of magnetic and Ne disturbances related to 12 earthquakes of various magnitudes. The satellite of the CSES mission that is dedicated to monitoring electromagnetic fields was launched on 2 February 2018. It is equipped with instruments that are similar to Swarm’s, capable of observing electric and magnetic fields from a Sun-synchronous orbit at an altitude of around 500 km, with an inclination of around 97.3°. There are already works based on CSES data related to selected earthquakes by Huang et al. [44] and groups of similar large earthquakes by Yan et al. [45]. Xiong et al. [46], as well as Yan et al. [45], analyzed CSES data focusing on precursory signals with respect to earthquakes. Marchetti et al. [42] have also investigated possible LAIC effects prior to the M = 7.5 earthquake using, among the others, Swarm and CSES datasets. De Santis et al. [47] have prepared a statistical spatiotemporal analysis of CSES-observed disturbed ionospheric signals due to an earthquake occurrence, which is quite analogous to the analyses made with Swarm data in [48]. He et al. [49] significantly extended the analyses of Swarm Ne disturbances in relation to seismic activity. They analyzed the global distribution of Swarm Ne disturbances and compared them to areas of higher and lower seismic activity rather than to specific earthquakes. Their analysis has a commonality with our current study, as we also try to focus on the continuous global seismic activity of various magnitudes. They analyzed global Ne disturbances statistically, with no spectral tools. In contrast, we try to analyze one important place locally and decompose disturbances in their frequency domain. The Remote Sens. 2023, 15, x FOR PEER REVIEW 4 of 37 advantage of LEO satellites is that they provide data in places that cannot be reached by ground GNSS observations. The orbital tracks of LEO can reach all seismically active zones of the world. A drawback of LEO is that the satellite data are dense only along the orbital tracks. The one-dimensional nature of LEO along-track Ne observations and sometimes the limited chance for spatiotemporal correlation due to the nonrepeating orbits leads to a strong requirement for spectral analysis for better recognition of the signals. The spectral analysis of residual signals can be found several times in the studies on seismically induced ionospheric disturbances in ground GNSS data [17,18,50]. It is possible that the processing of LEO data requires a spectral approach more than time series from a ground station do, as the disturbed signal is less pronounced in transient LEO observations. Fortunately, the frequency domain has been analyzed more frequently in LEO data in recent years [35,51–54]. These studies are encouraging us to take a further step and try a spectral classification of specific patterns that refer to different Ne and TEC signal disturbances. The proposed research aims at the investigation of LAIC via short-term Fourier transform (STFT) analysis of disturbed Ne, measured by Swarm over tectonic plate boundaries. This facilitates the analysis of seismically induced ionospheric disturbances with respect to larger earthquakes and their groups within the earthquake catalogue. The validation of disturbances detected by Swarm satellites in times of strong earthquake activity is based on the analysis of analogous data from the periods and places of lower seismic activity. The spectral decomposition of the disturbed signal can significantly help in the inspection of spectral behavior of ionosphere anomalies and the identification of the most frequent wavelengths of their occurrence. These determined spectral patterns, after more extended research in the future, can be classified as different types of disturbed signals, which was previously suggested in [17,18,51,52]. The extended spectral analyses of LEO electric field data can, in the future, help in the recognition of satellite data disturbances related to different seismic activities like mainshocks, aftershocks, preseismic activity, and also to moderate seismic activity in the vicinity of tectonic boundaries. The tectonic plate boundaries are weak places of the Earth’s lithosphere and can trigger preseismic, co-seismic, and post-seismic ionospheric disturbances. Swarm satellites visit these boundaries at different longitudes on the subsequent days and, due to their high speed, collect short samples of observations along their orbits. Therefore, the current study investigates two issues: firstly, a comparison of three-month-long records of Ne disturbances to the continuous seismic activity and not only to the largest earthquakes; secondly, a preliminary STFT recognition of Swarm Ne disturbances. In both cases, Swarmderived Ne are compared with continuous seismic records acquired from the United States Geological Survey (USGS) in the regions of interest. This study investigates nine three-month datasets, but primarily one example of a nearly latitudinally located boundary of the Australian and Pacific Plate that is rich in nearby minor plates and is very active. Swarm A/C and B tracks over this boundary are analyzed within a longitude range of around 30°. The dependence of the seismicity level on the season is very often discussed among researchers. Moreover, seasonal periodicities were already estimated in some works [55,56]. Therefore, at the beginning, we apply three-month-long series of nighttime Swarm data in Papua New Guinea from three different seasons of the year and compare them with the corresponding lists of earthquakes with a M ≥ 5.0. The fourth and fifth data series are from the Aleutian Islands at alternative latitudes and Africa at alternative longitudes. These series are analyzed to better validate Swarm’s sensitivity to different seismic activity levels and solar conditions. The remaining four datasets are derived from daytime Swarm passes in Papua New Guinea and analyzed with respect to daytime signal sensitivity. To observe the influence of various factors, Swarm Ne disturbances are analyzed in the time domain together with seismicity, solar radio flux (F10.7), geomagnetic ap and Dst indices, and the ionospheric bubble index (IBI), indicating equatorial plasma bubbles (EPBs) along the satellite tracks [57–59]. The primary graphical inspection of the correlation of Swarm Ne disturbances with seismicity or with ap, Dst, and F10.7 is summarized by the Pearson correlation coefficients (PCC), calculated between all pairs of Remote Sens. 2023, 15, x FOR PEER REVIEW 5 of 37 parameters and also between their 5, 7, and 10-day averages. PCC validated mutual correlations over time between Swarm Ne variations, seismicity, and solar and geomagnetic parameters. The most pronounced correlation is observed between Swarm Ne irregularities and the cumulative seismic energy at nighttime. This is promising for the future tracking of seismic activity over tectonic plate boundaries using spectrally analyzed LEO data but is at the same time challenging and needs more research employing global data. 2. Methodology This work is based on the calculation of time-variable power spectral density (PSD) of the high frequencies (short periods) of Swarm Ne 2 Hz data from LPs, starting from around 65 sec. We use wavelengths as frequency domain units, mainly for a better geometrical understanding of the disturbance size found at specific frequency. This unit is derived from the recalculation of Swarm data time units to distances along the orbit, using Swarm speed. It is therefore possible that two frequency domain units are applied alternately, i.e., period and wavelength. However, the term “period” is used also in reference to selected time epochs, which can be recognized from the context. The high-pass filtering of selected track parts of Swarm LP data is based on discrete Fourier transform (DFT), which is more universal than wavelets in the band-filtering processes, as it can decompose discrete signals in terms of sinusoids or complex exponentials. The signal of time function x(t) can be decomposed into frequency responses X(ω), sampled for the entire ω spectrum. These responses can also be transformed back exactly to x(t) via the inverse Fourier transform. The selected signal band can be subtracted from the data, which, contrary to simple methods of detrending, e.g., moving average or polynomial, removes exactly the part of the signal that we want to remove in terms of a frequency band. The nonspectral methods can remove arbitrary frequencies/wavelengths if the signal variances and spectral content differ in different places. The amplitude of signal responses in the frequency domain provides the PSD at individual wavelengths. The amplitude of spectral peaks usually varies with time, and this is the case with Swarm signals recorded along its orbital track. We have to carefully and accurately analyze the time-evolution of ionosphere disturbances measured by Swarm, because the speed of the satellite is fast with respect to sensed disturbances, and therefore, signal samples are composite but short. Time-variable PSDs plotted over time have the ability to map such changes in the frequency domain. The windowed STFT is a modification of the DFT method, which computes the spectrum of overlapping segments of the time series. The output of STFT is the short-term, time-localized evolutionary PSD. The data sequence to be transformed is multiplied using a window function, which is nonzero for only a short period of time [60,61]. The STFT of the signal is computed as the window is sliding along the time axis, resulting in a two-dimensional representation of the signal. Mathematically, this is written as 𝑆𝑇𝐹𝑇𝑥,𝑓,𝑡 =∑𝑥𝑛𝑤𝑛−𝑡𝑊𝑁−𝑓𝑛 𝑓=1,2,…,𝑁−1 𝑁−1 𝑛=0 (1) where 𝑊𝑁=𝑒𝑖2𝜋 𝑁 (2) and 𝑤𝑛 is a discretized window function, and t is its time index. The window applied in this study is the Tukey window. This window, also known as the cosine-tapered window, can be regarded as a cosine lobe of width 𝑟 2(𝑁−1) that is convolved with a rectangular window of width (1−𝑟 2)(𝑁−1). The equations defining Tukey window are Remote Sens. 2023, 15, x FOR PEER REVIEW 6 of 37 𝑤𝑇𝑈(𝑛)= { 1 2[1+𝑐𝑜𝑠(2𝜋 𝑟𝑛−1 𝑁−1−𝜋)] , 𝑛<𝑟 2(𝑁−1)+1 1 , 𝑟 2(𝑁−1)+1≤𝑛≤𝑁−𝑟 2(𝑁−1) 1 2[1+𝑐𝑜𝑠(2𝜋 𝑟−2𝜋 𝑟𝑛−1 𝑁−1−𝜋)] , 𝑁−𝑟 2(𝑁−1)<𝑛 } (3) where r is a value between 0 and 1. STFT of a data sequence 𝑥𝑛 describes its local spectral content near the moment t as a function of 𝑓 . Moving the center of the window 𝑤𝑛−𝑡 along the real line (of real units) allows for obtaining snapshots of the time-frequency behavior of 𝑥𝑛. STFT is used in the calculation of the amplitude of frequency spectrum for each sample. These spectra then form a 2-D shape that depends on the window size and overlay size. This corresponds to the computation of the STFT squared magnitude of the signal sequence 𝑥𝑛—that is, for a window width n: 𝑆𝑃𝑥𝑥,𝑓,𝑡 =|𝑆𝑇𝐹𝑇𝑥,𝑓,𝑡|2 (4) The size and shape of the analysis window can be modified depending on the purpose. A shorter window will produce more accurate results in terms of timing at the expense of the precision of frequency representation. A longer window will provide a more precise frequency representation at the expense of precision in time location. 3. Swarm and Ancillary Data The region of Papua New Guinea is located at the junction of four main tectonic plates: the Australian, Pacific, Eurasian, and Philippine Sea Plate. There are also several smaller plates, which makes this region especially seismically active. The main time period selected for the analysis due to the strong seismic activity is the winter period covering November and December of 2016 and January of 2017. The two additional threemonth periods for the validation are the spring period from February to April of 2016 and the summer period extending from May to July of 2016. The general M = 5.0 threshold for the selection of earthquakes in this work is assumed approximately on the basis of existing works on the sensitivity of ground GNSS data to seismic activity (see: [14,62]), as well as from the recent works on the sensitivity of Swarm magnetic data to earthquakes [37,42]. November of 2016 is relatively seismologically quiet and includes a week during which no earthquake of M ≥ 5.0 was recorded in the investigated region. December 2016 in Papua is very active seismologically, with two earthquake events having magnitudes of 7.8 and 7.9, both triggering tsunamis. January is not as quiet as November, including more earthquakes with M ≥ 5.0. It includes also one having a M = 7.9; however, this was deep, with no following aftershocks. It should be pointed out that based on the earthquake number, the tectonic plate junction seems to have been more relaxed in January 2017 than before December 2016. February–April of 2016 and May–July of 2016 are seismically lower active periods compared with to winter. However, in this especially active region, these periods also include a significant number of earthquakes exceeding 5.0 and 6.0 magnitudes. The earthquake records for all three seasons in the Papua New Guinea region are selected using the geographical window of 10°N–30°S and 130°E–170°E. Additionally, two other regions of Swarm data analysis are selected for validation and confutation purposes. These regions use geographical windows selecting Swarm Ne data at the equatorial latitudes first, but at a significantly different longitude, and second, at an absolutely different latitude, far from the equator. For the African region, we applied the window of 10°N–30°S and 0°E–40°E. The region of Aleutian Islands was selected using the window of 30°N– 70°N and 160°E–130°W. We used geodetic latitudes, as the difference with respect to geomagnetic ones is less significant here than the selection of the place that is free from a latitudinally aligned tectonic plate boundary, which is Africa. The data in both ancillary regions cover the primary period of Papua New Guinea data from November 2016 to January 2017. Remote Sens. 2023, 15, x FOR PEER REVIEW 7 of 37 The list of earthquakes is extracted from the Search Earthquake Catalog, which is available from the USGS website. There were almost 200 earthquakes in the Papua New Guinea region with a magnitude equal to or above M = 5.0 from the beginning of November 2016 to the end of January 2017 (Figure 1). A significant part of these earthquakes included aftershocks and occurred just after 8 and 17 December 2016. The orbits of Swarm satellites, which have an inclination between 87° and 88°, pass the boundary of two major subducting tectonic plates (Australian and Pacific) at an angle that is close to 80°. Figure 1. The earthquakes over a magnitude of 5.0 in the Papua New Guinea region from November 2016 to January 2017 and tectonic plate boundaries in Papua New Guinea Region. Three major earthquakes (7.8–7.9) in December and January are indicated with white labels. The analysis of the Swarm data in the Papua New Guinea region is based on the LP data from Swarm Level 1b Electric Field Instrument (EFIxLP). The longitude (λ) of the Swarm satellite pass does not repeat the next day at the same time, and moreover, it is not exactly the same. Therefore, we selected Swarm passes from consecutive days of the threemonth periods, having close longitudes and being recorded at a similar time of the day. The selection of Swarm A and C orbital tracks for the analysis of winter solstice starts from Swarm A and C trajectories recorded on 17 December 2016, which both took place 40 min. after the largest M = 7.9 earthquake occurred close to the eastern coast of New Ireland island in Papua New Guinea (10:51 UTC, which is 20:51 local time (LOC)). Then, for each satellite (A/C), an additional 90 Swarm tracks from the preceding and subsequent days are selected (1 November 2016–30 January 2017). The threshold difference in the 90 adjoined tracks in the time unit, in relation to that on 17.12, is set to a maximum of 4.4 h, whereas the allowed difference in longitude is a maximum of 20°. The range of latitudes Remote Sens. 2023, 15, x FOR PEER REVIEW 8 of 37 of the 91 selected tracks is approximately between 25°N and 45°S (Figure 2a). This selection produced 91 track sections, shown in Figure 2a, and the same rule with reference Swarm A/C tracks is applied in a more quiet region of the African Plate to select a similar number of tracks at very similar latitudes (Figure 2b). In order to keep Figure 2a transparent, the date labels are only shown for the eight selected days and are also presented in Figures 4–6, including the most significant disturbances within the selected period. In the third selected region, which is the junction of the Pacific Plate and the North American Plate located at completely different latitudes, the reference Swarm track is selected from the 17th day of the middle month (Figure 2c). The longitude differences and time differences of the adjoined tracks in relation to the reference track are the same as for the primary region. The Swarm orbital track selections in the Papua New Guinea region for the other seasons of 2016 are prepared in exactly the same way, starting with the track within very similar latitudes on the 17th day of the middle month, and adjoining remaining Swarm passes at close longitudes using the same longitude and time separation maxima. Additionally, the same rule is applied for the selection of the opposite daytime tracks that were investigated, aside from the above-mentioned nighttime or nearly nighttime passes. Neither the track selection in the summer solstice and vernal equinox nor the daytime tracks need to be plotted in the figures, as these tracks follow the same rule that can be seen in Figures 2a and 3a. They pass tectonic plate boundaries within very similar latitudes, similar longitudes, and similar times, dependent on the actual Swarm trajectory, but always refer to the selected half of the day. Figure 2. (a) Selection of descending (nighttime) Swarm A/C tracks over Papua New Guinea between 1 November 2016 and 30 January 2017. The tracks are selected to have the closest possible time of the day and longitude to the first selected track just after M = 7.9 mainshock. The red triangles denote earthquakes above M = 6.5. (b) Selection of Swarm A/C tracks at the same latitudes and local time in amore seismically quiet region. The red triangles denote earthquakes above M = 4.0. (c) Selection of Swarm A/C tracks at the high and middle latitudes. The red triangles denote earthquakes above M = 4.5. The selection of Swarm B tracks starts from the Swarm B trajectory recorded on 8 December 2016, approximately 30 min. before the largest M = 7.8 earthquake in the Solomon Islands region near Papua New Guinea, which occurred close to the western coast of San Cristobal Island (17:38 UTC—3:38 LOC). Then, an additional 90 Swarm tracks from the preceding and subsequent days (1 November 2016–30 January 2017) are selected. The range of latitudes of the 91 selected tracks is approximately the same as for the Swarm A/C pair. The allowed between-track differences in time and longitude with respect to a Remote Sens. 2023, 15, x FOR PEER REVIEW 9 of 37 trajectory of 8.12 are also the same. The tracks of Swarm B are grouped in a slightly different way than Swarm A/C, as the orbit repetition cycles of Swarm B (Figure 3a) differ in comparison to those for Swarm A/C. The selection of Swarm B passes in a more seismically quiet area is shown in Figure 3b. The selection of Swarm B passes at different latitudes is shown in Figure 3c. It should be pointed out that the reference track selection is unified with the rules of the Swarm A/C reference track selection, which is always selected on the 17th of the middle month. The selection of Swarm B tracks in the summer of 2016, spring of 2016, and during the daytime has followed the same rule of very similar latitudes, similar longitudes (±20°), and similar times, separated by 4.4 h from the reference track. Such an apparently long separation in time is necessary to find Swarm tracks with success on every day and to have the smallest number of gaps. Figure 3. (a) Selection of descending (nighttime) Swarm B tracks over Papua New Guinea between 1 November 2016 and 30 January 2017. The tracks are selected to have the closest possible time of the day and longitude to the first selected track just before M = 7.8 mainshock. The red triangles denote earthquakes above magnitude 6.5. (b) Selection of Swarm B tracks at the same latitudes and local time in a more seismically quiet region. The red triangles denote earthquakes above magnitude 4.0. (c) Selection of Swarm B tracks at the high and middle latitudes. The red triangles denote earthquakes above M = 4.5. Figures 2 and 3 present the rule of selection of Swarm A/C and B trajectories—one each day at a similar time of the day. The three-month sets of around 91 trajectories for different satellites, regions, and times of the day are then prepared for the comparison of ionospheric Ne perturbations to seismicity, solar activity, EPBs, and geomagnetic activity. Three seasons are analyzed in Papua New Guinea in the nighttime and daytime, which provides six three-month series. Two other series come from the winter solstice in Africa and Aleutian Islands. The other dataset includes daytime data in Papua New Guinea, with Ne irregularities analyzed at the same frequency as the nighttime data. Three other daytime analyses apply sampling at the higher frequency. 4. Spectral Patterns of Swarm Ne Disturbances during Nighttime of the 2016–1017 Winter From the winter set of the 91 trajectories within the three-month period, eight are selected for each Swarm spacecraft to illustrate the spectral characteristics of the most important moments in the spectrograms (see labels on Swarm passes in Figures 2 and 3). Besides that, the measures of the disturbance scale for the entire three-month comparisons with seismicity are also determined from the spectrograms as the maximum PSD of Ne disturbances at the selected wavelength (from here, simply called Swarm Ne PSD). The Remote Sens. 2023, 15, x FOR PEER REVIEW 16 of 37 The high-pass filtered (at 500 km wavelength ≈ 65 s period) in situ measured Ne in the 91 orbital passes of each of the three satellites are analyzed spectrally via STFT, and the spectrograms are prepared in the same way as in Figures 4–6. However, the presentation of 273 spectrograms in this research paper would be impossible. Therefore, the spectrograms have been sampled at a 300 km wavelength (~40 s), and the maxima from the vectors of the sampled PSDs have been collected as daily time series. These maxima are presented with the use of a light-blue area plot, and we also call these “Swarm Ne PSD” to maintain consistency along the entire study. Swarm Ne PSD values are applied in the assessment of Swarm’s sensitivity to seismic activity in the region over time. The selection of a 300 km wavelength is a preliminary choice, and it does not mean that 300 km is the most representative frequency domain wavelength for SIID detection. However, this wavelength looks to be approximately located in the middle of the disturbances found in Figures 4–6, and it is suspected to be sensitive to seismic activity. Therefore, Swarm Ne PSD sampled at 300 km from the spectrograms is used as an indicator of Swarm’s sensitivity to the seismic record, F10.7, EPB, and geomagnetic ap and Dst indices. Additionally, observing a high sensitivity of the LPs and quite dynamic changes in Swarm disturbances over the tectonic plate boundary, and assuming a complex and long-term evolution of seismic activity, the 5-, 7-, and 10-day moving averages are also calculated from Swarm Ne PSD to reveal the more long-term character of the Ne changes. Averaging in time can also mitigate two other disadvantages: strong short-term Ne variations affecting the correlated temporal Ne change and differences between the longitudes of consecutive orbital footprints. The 10-day averages are plotted as bold blue lines in Figures 7–15. The earthquakes are visualized by their depths and magnitudes. The magnitudes are multiplied by 10 to easily maintain the consistency of the axes in Figures 7–15. The ap index and the Dst index are multiplied by two due to the same reason. The 5-, 7-, and 10-day moving averages are also calculated for ap and Dst indices, and their 10-day averages are similarly plotted in Figures 7–15 using green and orange bold lines, respectively. The F10.7 parameter in the solar flux units (SFUs) originally have a compatible scale here, but the numbers of Swarm epochs with confirmed and unconfirmed plasma bubbles along the individual tracks have to be divided by 10. They are drawn in the bottom direction to be more readable amongst the other parameters. Furthermore, the Swarm Ne PSD maxima and their moving averages have arbitrary and nonrelated scales that are used to clearly visualize their most interesting relative variation. This does not affect the results, as the scaling has no influence on the correlation observation. The scales of Swarm Ne PSD and its moving average are equal for Figures 7, 8, and 10–12 (equatorial latitudes), increased six times in Figure 9 (high latitudes) in comparison to Figures 7 and 8, and increased a hundred times in Figures 13–15 (daytime, sampling at higher frequency), because the signals at different wavelengths have different amplitudes. Based on the earthquake magnitude and depth and knowing that the shallowest earthquakes are the ones that are most suspected of affecting the ionosphere [23], we calculated an equivalent of the so-called cumulative Benioff strain, which corresponds to the cumulative (square root sum) of the seismic energy [66–68]. The quantity calculated in this study is a kind of weighted cumulative sum, because the square roots of energy approximation are divided by the focal depths, i.e., 𝑐𝑠= ∑√10𝑀(𝑖)/𝐻(𝑖) 𝑛 𝑖=1 , where n is the number of earthquakes in the selected time window, and H(i) are their focal depths. This sum of approximate earthquake energy, calculated within selected 5-, 7-, and 10-day windows, will from here be called the cumulative sum of seismic energy. The quantity that is calculated this way reveals more long-term variation in seismicity, and similarly to the moving average referred to as Swarm Ne PSD, it is drawn using bold lines of the same color as its input quantity, i.e., a black bold line (Figures 7– 15). It is proportionally rescaled in Figures 7–15, except in Figure 9, where it is increased six times. Figures 7–15 in the following chapters include the labels with time windows of the Swarm track selections recalculated to both UTC and LOC. The first high peak of the Swarm Ne PSD of satellite A at the wavelength of 300 km in Figure 7a occurs before November 10, which is close to an earthquake with M = 5.9. Remote Sens. 2023, 15, x FOR PEER REVIEW 17 of 37 Another high peak comes after a relatively calm week on the first day with an earthquake with M = 5.0 (November 27, see also Figure 4c). Next, high peaks repeat twice before December 8. The other strong signals occur on the days of major earthquakes and one day after them, i.e., on December 8–9 and December 17–18, which is illustrated by the spectrograms in Figure 4d–g. The five particular days, i.e., the day of the first M = 5.0 (November 27), the days of major earthquakes, and the days after them, are denoted with blue dots in Figure 7a (but so are the other days presented in Figure 4). This more clearly reveals that Swarm A disturbances at a 300 km wavelength increase at that time. The same days that Swarm C passes are marked with large dots in Figure 7c, and the behavior of variations in the Swarm Ne PSD follows a similar trend as in the case of Swarm A. These observations suggest that Swarm Ne PSD variations in Swarm A and Swarm C in the analyzed time are noticeably correlated with M ≥ 5.0 seismic activity in the region of Papua New Guinea. The analogous Figure 7b for Swarm B also includes five blue dots indicating the first M = 5.0, two major earthquakes, and two days after the major earthquakes, as well as three dots denoting the remaining days from Figure 5. The Swarm Ne PSD of Swarm also increases here with respect to the neighboring small values, which confirms their correlation with seismicity. Nonetheless, not all Swarm B disturbances at large dots are the largest throughout the three-month period considered. The largest peaks for Swarm B, which are at the same time larger than all values observed here for the Swarm A/C pair, are observed on December 12 and January 8 (Figure 7b), the latter of which has the spectrogram in Figure 5h. This is an example spectrogram of a larger-scale ionospheric disturbance detected by Swarm B, and we should add that the spectrogram for December 12 is similar. The disturbances of this type exceed the area of the tectonic plate boundaries and dominate in Figure 7b over the signals suspected to be seismically triggered. The geomagnetic ap index reaches 4+ on January 8 but is small on December 14, and Dst also exhibits no variations. Thus, the geomagnetic origin of these large disturbances, like those in Figure 5h, is suspected but not confirmed. There are even no unconfirmed EPBs from the L2 IBI product in Figure 5h, and therefore, the analysis of such signals is open, with a strong indication of geomagnetic activity. The F10.7 in Figure 7 increases a few days before the first major earthquake on December 8, but also at the time for the third major earthquake on January 22. The first F10.7 increase occurring at the beginning of December coincides with repeating M ≥ 5 earthquakes before the major earthquake on December 8 and also with an increase in disturbed Swarm signals (Figure 7a–c). However, on January 22, the increase in F10.7 is co-located with the major earthquake, rather than with the Swarm A/C disturbed signals (Figure 7a,c). Nevertheless, some moderate disturbance of the Swarm B signal can be detected there (Figure 7b). The comparative analysis of Swarm Ne PSD with respect to the series of earthquakes in the Papua New Guinea region confirms the high sensitivity of the Swarm LP instrument and even the dynamic character of the high-frequency ionospheric Ne changes at around a 500 km altitude. In order to slightly average rapid the variations in residual high-frequency Swarm signals and look at their longer-term behavior, the moving average from 10 days is presented together in Figures 7–15. It is easy to guess from this averaged Swarm Ne PSD (bold blue lines in Figure 7) that a more long-term Swarm high-frequency signal increase starts several days before the time of the two major earthquakes in the case of Swarm A/C (Figure 7a,c). In the case of Swarm B, this long-term signal increase is not so rapid at the beginning, but it also starts before December 2016 (Figure 7b). These findings correspond to the findings that are related to a precursory character of Swarm signals in [43,48]. Remote Sens. 2023, 15, x FOR PEER REVIEW 18 of 37 Figure 7. Maxima of PSDs of residual Ne sampled at 300 km wavelength from individual tracks on different days (light-blue area plot) for (a) Swarm A, (b) Swarm B, (c) Swarm C in the nighttime in the winter of 2016–2017. The figures also present 10-day moving averages of Swarm Ne PSD (bold blue lines), earthquakes with M ≥ 5.0 (magnitude multiplied by 10—black stems with dots, depth— black stems with circles), cumulative sum of seismic energy within 10-day period described in Section 5 (bold black), 10-day average of ap index multiplied by 2 (bold green), 10-day average of Dst index multiplied by 2 (orange line), F10.7 (red line), confirmed EPB number divided by 10 (red bars), and unconfirmed EPB number divided by 10 (blue bars). Swarm passes denoted with blue dots have spectrograms in Figures 4–6. The label in the right-bottom corner describes the range of times of the day for the selected tracks. In order to validate Swarm LPs sensitivity over a seismically active region, an alternative region that is free from tectonic activity is also selected. It is intentionally selected to be far from the Pacific Plate and located over the African Plate, west of the Papua New Guinea region. This region does not include latitudinally arranged seismically active tectonic plate boundaries. Additionally, there were almost no earthquakes in the winter of 2016/2017 at the tectonic junction between the African, Somalia, and Lwandle plates. The tracks are selected there at a possibly similar period of the same days, i.e., passing the region at evening and night hours. The T selected tracks have a very similar range of geodetic latitudes to the Papua New Guinea region, which also leads to similar geomagnetic Remote Sens. 2023, 15, x FOR PEER REVIEW 19 of 37 latitudes in the case of these two regions. We have selected three months of tracks using the same rules of Swarm pass selection defined in Section 3, and these new tracks are shown in Figures 2b and 3b. The seismic relaxation, which was not checked before the track selection, turned out to be more perfect than expected. The scale of Swarm Ne PSD and its 10-day average applied here (Figure 8) is the same as in Figure 7. The main finding in Figure 8 is a very small level of disturbances at a 300 km wavelength in Ne signals from the Swarm A/C pair (Figure 8a,c). Figure 8b, aside from some surrounding data gaps in Swarm B observations, shows a potentially precursory disturbance to a very shallow earthquake on January 27. Summarizing Figure 8, the level of Swarm Ne disturbances resulting from almost no stronger seismic activity is, as expected, very low in the case of the African region that is located at exactly the same latitudes as Papua New Guinea. Moreover, the highest increase in Ne PSD in Figure 8b even unexpectedly correlates with the only shallow M ≥ 5.0 earthquake. The low Swarm Ne PSD at a 300 km wavelength in Figure 8 confirms that the Sun being at low latitudes and an EIA related to the equatorial latitudes cannot trigger large disturbances alone. The remaining question is if a latitudinally oriented active tectonic plate boundary at higher latitudes can contribute to Ne disturbances at Swarm altitudes. To analyze this issue, the Aleutian Islands are selected additionally, as this region includes a subduction boundary that is directed orthogonally to the Swarm orbit, similarly as in the Papua New Guinea region (Figures 2c and 3c). The winter of 2016/2017 is selected here, as was also the case for Papua New Guinea and Africa. The Swarm A/C and B nighttime tracks are selected and processed in the same way as previously selected data, and the Swarm Ne PSDs at 300 km are summarized in Figure 9. The magnitudes of Swarm Ne PSDs are lower here compared with the equatorial region, probably due to the generally lower electron concentration. The cumulative sum of seismic energy is also lower due to the lower seismic activity at that time. Therefore, both quantities are multiplied here by six in comparison to Figures 7–8. After multiplication by six, we can see at least three epochs of correlation of Swarm disturbances with the earthquake record at this tectonic plate boundary (Figure 9). The Swarm Ne PSD, as well as its moving average, is most often increased in Swarm A, B, and C signals between 10 and 20 of November, when the seismic activity is also the largest (Figure 9a–c). The second highest disturbance in Ne signal of all three Swarm spacecrafts occurs in the first week of December, together with the shallow earthquake on 1st December. The other interesting disturbances are recorded by Swarm A and C on January 18 and coincide exactly with the other shallow earthquake. It is worth noticing that the number of extreme 10-day average ap or Dst values does not correspond to the number of Swarm Ne PSD increases. Instead, the Swarm Ne PSD increases always approximately correspond to the groups of earthquakes or single shallow earthquakes, which is most evident in Figure 9b. Remote Sens. 2023, 15, x FOR PEER REVIEW 20 of 37 Figure 8. The same as in Figure 7, but for the Africa region of more quiet seismic conditions. Subfigure titles (a–c) refer directly to individual Swarm satellites. The earthquakes shown along the time axis have M ≥ 5.0. Remote Sens. 2023, 15, x FOR PEER REVIEW 21 of 37 Figure 9. The same as in Figure 7, but for the Aleutian Islands region, with tectonic plate boundary orthogonal to Swarm track at 55°N. Subfigure titles (a–c) refer directly to individual Swarm satellites. The earthquakes shown along the time axis have M ≥ 5.0. The scales of Swarm Ne PSD maxima and cumulative sum of seismic energy are multiplied by 6 in comparison to Figures 7 and 8. 6. Swarm Ne Disturbances during Nighttime of the 2016 Summer and Spring and M5+ Earthquakes In order to examine the sensitivity of LPs measuring onboard the Swarm constellation in different seasons of the year, two other seasons are selected, i.e., the summer solstice (May–July 2016) and the vernal equinox (February–April 2016). The Swarm data and earthquake data are selected using exactly the same rules as for the winter solstice season. Around 90 tracks from three months are selected—one for each nighttime. The region is the same as that in Figures 2a and 3a, and the same rules for latitudes, longitudes, and time have been applied. The starting reference track is always selected on the 17th day of the middle month, and then around 90 tracks are matched with the starting one, depending on the Swarm orbit availability in the selected spatiotemporal windows. It should be mentioned that limited ranges of longitude and time could sometimes cause a loss of single tracks from the set, and this is visible as gaps in the blue lines of the average Swarm NE PSD in Figures 7–15. However, the gaps are always small and have no significance for the entire study. As mentioned in Section 3, the track selection in the summer and in the Remote Sens. 2023, 15, x FOR PEER REVIEW 22 of 37 spring are not shown to save the space, and also because these sets have a similar spatial distribution to that from the winter. Figure 10 describes the summer case when a specific seismic activity level occurs. This level is lower in comparison to that from Figure 7, the number of earthquakes is smaller, and they only reach M = 6.4. Nevertheless, we can observe longer series around June 15–22 and shorter ones around July 25, when Dst does not indicate a storm, and the series occur instead in the phase of solar activity decrease. The second shorter series has indicated unconfirmed plasma bubbles at this time, but there are no plasma bubbles at the time of the first series. Both series of earthquakes coincide approximately with the largest Swarm A/C Ne disturbances at a 300 km wavelength (Figure 10a,c). There are smaller peaks of Swarm Ne PSDs that are precursory to the first earthquake series around June 9 in Swarm A/C, coinciding with increased geomagnetic activity, and the largest peak just is after the series (June 22), which happens at low ap/Dst indices. The other strong peaks of Swarm Ne PSDs are around July 16–19 (Swarm A and C), and the solar activity is also the strongest at that time. However, it should be noticed that at the time of the second earthquake group, which started around July 18–19, the solar activity increased in parallel to the seismic activity and Swarm disturbed signal occurrence. Swarm B also recorded disturbances at that time, and if we take a more detailed look at Figure 10b, we can easily find that all the small peaks from Swarm B coincide well with the seismic activity increase. Figure 11 describes the spring case, when, again, the seismic activity is lower than in the winter, and again includes specific increases in the form of earthquake groups. The maximum magnitude of an earthquake reaches 7.0 in this case. The most interesting epochs during these three spring months are the earthquake series that started around the beginning of April and lasted for almost three weeks and three very close-together and very shallow earthquakes on March 18. The longer earthquake series in April coincides with all other factors, i.e., the solar activity reaches the maximum at that time, ap 10-day average has a moderate increase here, as does Dst. Confirmed EPBs also occur on some days and coincide with the Swarm Ne PSD peaks. Therefore, in April 2016, it is hard to credibly assess the reason for the very high Swarm Ne PSD peaks, but at the same time, a seismic source is hard to exclude. Additionally, if we take a look at two slightly smaller, but still considered large, Swarm Ne PSD peaks in Figure 11a (March 19 and April 6), we can see that F10.7 is not high at that time, the ap and Dst are low or moderate, but groups of shallow earthquakes are very close to Swarm disturbances. Summarizing Figure 11, we can say that there are more cases of large Swarm Ne PSDs at 300 km that coincide with some characteristic seismic activity than those coinciding closely with the highest average ap, particularly anomalous Dst or solar radio flux. Remote Sens. 2023, 15, x FOR PEER REVIEW 23 of 37 Figure 10. The same as in Figure 7, but for the summer of 2016. Subfigure titles (a–c) refer directly to individual Swarm satellites. Remote Sens. 2023, 15, x FOR PEER REVIEW 24 of 37 Figure 11. The same as in Figure 7, but for the spring of 2016. Subfigure titles (a–c) refer directly to individual Swarm satellites. 7. Swarm Ne Disturbances during Daytime of the 2016–2017 Winter and Daytime of the 2016 Summer and Spring and M5+ Earthquakes The selection of Swarm trajectory parts in Sections 4–6 was focused on the evening or night hours, which was motivated by the severe observational conditions during the day. Indeed, it will be shown that the sensitivity of Swarm to Ne differs in the daytime. This section investigates daytime Swarm tracks, which, as in the former sections, can differ by some hours between the Swarm A/C pair and Swarm B due to the constellation geometry. However, each Swarm spacecraft visits the same region twice daily (one ascending and one descending track), and it is possible to find the passes closer to midnight, as well as those closer to noon. The Swarm Ne PSD sampled at a 300 km wavelength from the daytime spectrograms in the winter of 2016–2017 is presented in Figure 12. The Swarm Ne PSD at 300 km has a significantly smaller amplitude in the daytime, which is shown by low values of light-blue area plots and bold blue lines. It can, however, be concluded from the entire Figure 12 that Swarm LP signals at a 300 km wavelength have a small amplitude and poorly coincide with the earthquakes, as well as with the average ap, Dst, or F10.7 series. The only conclusion with regard to the Swarm A/C pair (Figure 12a,c) is that among the Swarm disturbances in the daytime, at most, some low precursory signals can be found, rather than co-seismic ones. This small sensitivity of Swarm LPs at a 300 km Remote Sens. 2023, 15, x FOR PEER REVIEW 25 of 37 wavelength during sunny conditions encouraged us to inspect the shorter wavelengths, which have, however, much lower amplitudes of the signal. The sampling of Swarm Ne PSDs at higher frequencies provided some interesting findings, which appeared worth presenting and are therefore shown in Figures 13–15. Swarm orbital track samples are selected for the same three analyzed seasons in the daytime, and the calculated STFT spectrograms are sampled at a very short wavelength of 80 km, corresponding to around 10 sec of Swarm flight. The signals have significantly smaller amplitudes at a 80 km wavelength, and therefore, the scale of the Swarm Ne PSD and its 10-day average presented in Figures 13–15 is increased a hundred times in comparison to the case of signals sampled at a 300 km wavelength in Figures 7, 8, and 10–12. Figure 13 presents some peaks of the Swarm Ne PSD, but their time location differs between the satellites. Figure 13a for Swarm A in the winter daytime shows one peak of Ne signal at 80 km, around November 25. This is the time of relaxation from earthquakes with M ≥ 5.0, and it is possible that a higher ap and larger negative Dst contribute to that disturbance. On the other hand, there is a second peak at the beginning of December, which coincides exactly with a shallow earthquake at the time of a low ap and Dst around zero. These two peaks raise the 10-day moving average of the Swarm Ne PSD in such a way that it looks like a precursor to the earthquake series. Figure 13b (Swarm B) shows the first peak in the time of the moderately frequent earthquake series, and the second one coincides well with the major earthquake on December 8, 2016, and has many smaller peaks during the neighboring days. Swarm C, analyzed in Figure 13c, observed the disturbances after the second major earthquake on December 17 in the time of aftershocks, but started detecting very small disturbed signals at the beginning of December. Figure 14 shows the summer daytime, where the largest disturbances reflected in Swarm Ne PSD sampled at a 80 km wavelength from Swarm A/C (Figure 14a,c) occur several days before the strongest earthquake series that is available here. Swarm B has a large Swarm Ne PSD peak at the beginning of May, which occurs together with several shallow earthquakes (Figure 14b). Figure 15 describes the daytime of Swarm LP measurements over Papua New Guinea in the spring of 2016. The Swarm A/C pair, flying close to each other, again have similar disturbed signals. The highest peaks for Swarm A/C are during the earthquake series that started at the beginning of April (Figure 15a,c). The solar radio flux starts to increase a few days later, but this is another example where a high F10.7 coincides well with the earthquake series, as well as with Swarm disturbances. The other peaks occur on March 6 and at the beginning of February, and both determine increases in the 10-day average that are approximately correlated with the cumulative sum of seismic energy. Figure 15b also shows the highest increase in the Swarm Ne PSD from satellite B around the beginning of March, which also coincides with increased seismic energy. Remote Sens. 2023, 15, x FOR PEER REVIEW 32 of 37 Figure 16. Pearson correlation coefficients of Swarm A/B/C Ne PSD with cumulative sum of seismic energy (black), ap (green), Dst (orange), and F10.7 (red), respectively (a) directly between daily Swarm Ne PSDs, cumulative sum for the short period of 5 days, ap, Dst, and F10.7, (b) using 5-day moving averages of all series except F10.7 (original daily F10.7 is used due to its apparent long-term character), (c) 7-day averages of all with original F10.7, and (d) 10-day averages of all with original F10.7. N and D denote nighttime and daytime tracks, respectively. The wavelength of Ne signal sampling is also given in km. 9. Conclusions The preliminary observations of the common spectral patterns of ionospheric disturbances are promising for the future spectral recognition of short, along-track disturbed signals. Therefore, this study suggests that continuous tracking of seismic activity from the satellite orbit and spectral analysis of LEO along-track signals are fundamental and can facilitate the division of disturbed signals of different kind, e.g., distinguish preseismic from co-seismic ones. This also suggests that the processing of ground GNSS observations for LAIC studies should also use spectral analysis and must be considered not only at the time and place of the strongest earthquakes. Remote Sens. 2023, 15, x FOR PEER REVIEW 33 of 37 The visual observation of the correlations between Swarm Ne disturbances, seismicity, ap index, Dst index, and solar radio flux in Figures 7–15 has been summarized by PCCs, calculated between Swarm Ne PSDs from three Swarm satellites and four remaining quantities: the cumulative sum of seismic energy, ap and Dst indices, and F10.7. The PCCs that were determined for the original data series (except seismicity, for which the shortest-term 5-day sum is used) are notable (≥0.3) only in three nighttime cases, but even these cases confirm the best correlation of the Swarm–seismicity pair. The 5-day averaging of all data series (except F10.7 due to its smooth long-term character in daily series) leads to a slight increase in correlations while still maintaining an advantage of Swarm–seismicity pairs. The increase of the averaging window to 7 days confirms the dominant correlation between Swarm Ne PSDs and weekly seismic energy, with the PCC reaching or exceeding 0.3–0.4 in all nighttime cases. At the same time, no significant increase in the PCC for the other data pairs is attained after 7-day averaging of the analyzed series. The averaging with a 10-day window increases all nighttime correlations of Swarm–seismicity pairs to 0.4–0.7, whereas the other pairs infrequently obtain an PCC equal to at most 0.3. These numbers strongly suggest that longer-term (several day) disturbances in Ne at a 300 km wavelength, recorded by Swarm 450 km over the tectonic plate boundary, are correlated predominantly with seismicity. Moreover, the correlation of seismicity with nighttime Ne variations is strong even without the determination of the respective time shift between these two phenomena. This can suggest that taking into account the time shift resulting from, e.g., the precursory character of Ne disturbances, can provide an even larger PCC. The first rough conclusion from the analysis of daytime Swarm LP observations is that the daytime LP Ne signal is less sensitive to ionospheric disturbances from the lithosphere at a 300 km wavelength. However, high PCC values at an 80 km wavelength in Papua New Guinea in the winter and spring support the suspicion that by analyzing very high frequencies, we can detect LAIC from the Swarm orbit, also in the daytime. The long-term disturbances approximated by 7and 10-day moving medians of the Swarm Ne PSD reveal a notably better correlation with persisting long-term seismic activity than individual Ne signal disturbances with individual large earthquakes. The persisting Swarm Ne disturbances often precede significant seismic activity, which can point to their precursory character. However, there are also many examples of disturbances occurring just after shallow earthquakes, at the end of an earthquake series, or in the middle of the series, which proves that this analysis over a longer time span is relevant, whereas a selective approach to seismicity studied with respect to LAIC can hide some important information. Author Contributions: Conceptualization, W.J., P.W., R.H., and M.H.-P.; Methodology, W.J., P.W., and M.H.-P.; Software, W.J.; Validation, H.Y., E.M.-M., and V.G.; Formal Analysis, B.M., A.G.-R., and A.K.-G.; Investigation, W.J., P.W., and A.K.-G.; Resources, A.K.-G. and A.G.-R.; Data Curation, W.J. and B.M.; Writing—Original Draft Preparation, W.J.; Writing—Review and Editing, B.M. and M.H.-P.; Visualization, W.J.; Supervision, P.W.; Project Administration, R.H.; Funding Acquisition, W.J. All authors have read and agreed to the published version of the manuscript. Funding: This study was partially funded from the COSTO project, which has received funding from the ESA contract 4000126730/19/NL/IA in the framework of the “Swarm + Coupling: HighLow Atmosphere Interactions” Call, and partially from the Research Grant no. DEC2021/41/B/ST10/03954, which has received funding from the National Science Centre (NCN) of Poland. Data Availability Statement: Swarm data are downloaded from ESA Data Access Service (https://swarm-diss.eo.esa.int, ftp://swarm-diss.eo.esa.int). The earthquake series is downloaded from the United States Geological Survey service (https://earthquake.usgs.gov/earthquakes/search/). The Solar Radio Flux data are downloaded from the website of the International Reference Ionosphere (IRI) of the Committee on Space Research (COSPAR) (http://irimodel.org/). The ap and Dst indices are obtained from the International Service of Geomagnetic Indices (http://isgi.unistra.fr/data_download.php). Tectonic plate boundaries are downloaded from a data Remote Sens. 2023, 15, x FOR PEER REVIEW 34 of 37 webpage (https://github.com/lcx366/PlateTectonic). The Ionospheric Bubble Indices are taken from the Swarm L2 IBI Product (Park et al. 2013). Acknowledgments: The authors wish to thank the other participants of the COSTO project for the joint cooperation on closely related studies that were important for this study. They are from the National Observatory of Athens (Anna Belehaki, Ioanna Tsagouri, and Evangelos Paouris) and from the Technical University of Munich (Michael Schmidt, Andreas Goss, and Eren Erdogan). The authors wish to thank the three anonymous reviewers for their in-depth and very constructive comments. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations DFT discrete Fourier transform EFI electric field instrument EIA equatorial ionization anomaly EPB equatorial plasma bubble F10.7 solar radio flux GNSS Global Navigation Satellite System IBI ionospheric bubble index LAIC lithosphere–atmosphere–ionosphere coupling LEO low Earth orbiting LOC local time LP Langmuir probe Ne electron density PCC Pearson correlation coefficient PNG Papua New Guinea PSD power spectral density SFU solar flux units SIID seismically induced ionospheric disturbance STFT short-term Fourier transform TEC total electron content USGS United States Geological Survey UTC Universal Time Coordinated λ geodetic longitude ϕ geodetic latitude References 1. Harrison, R.G.; Aplin, K.L.; Rycroft, M.J. Atmospheric electricity coupling between earthquake regions and the ionosphere. J. Atmos. Solar-Terr. Phys. 2010, 72, 376–381. 2. Pulinets, S.; Boyarchuk, K. Ionospheric Precursors of Earthquakes; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2005 3. Sorokin, V.M.; Yaschenko, A.K.; Hayakawa, M. A perturbation of DC electric field caused by light ion adhesion to aerosols during the growth in seismic-related atmospheric radioactivity. Nat. Hazards Earth Syst. Sci. 2007, 7, 155–163. 4. Prucha, J.J. Zone of weakness concept: A review and evaluation. In Basement Tectonics 8, Proceedings of the International Conferences on Basement Tectonics, Duluth, MI, USA, 1–11 August 1992; Bartholomew, M.J., Hyndman, D.W., Mogk, D.W., Mason, R., Eds.; Springer: Dordrecht, The Netherlands, 1992; Volume 2. https://doi.org/10.1007/978-94-011-1614-5_7. 5. Pulinets, S.A.; Boyarchuk, K.A.; Hegai, V.V.; Kim, V.P.; Lomonosov, A.M. Quasielectrostatic Model of atmosphere-thermosphere-ionosphere coupling. Adv. Space Res. 2000, 26, 1209–1218. https://doi.org/10.1016/S0273-1177(99)01223-5. 6. Pulinets, S.; Davidenko, D. Ionospheric precursors of earthquakes and Global Electric Circuit. Adv. Space Res. 2014, 53, 709–723. https://doi.org/10.1016/j.asr.2013.12.035. Remote Sens. 2023, 15, x FOR PEER REVIEW 35 of 37 7. Hayakawa, M.; Izutsu, J.; Schekotov, A.; Yang, S.-S.; Solovieva, M.; Budilova, E. Lithosphere–Atmosphere–Ionosphere Coupling Effects Based on Multiparameter Precursor Observations for February–March 2021 Earthquakes (M~7) in the Offshore of Tohoku Area of Japan. Geosciences 2021, 11, 481. https://doi.org/10.3390/geosciences11110481. 8. Freund, F. Pre-earthquake signals: Underlying physical processes. J. Asian Earth Sci. 2011, 41, 383–400. https://doi.org/10.1016/j.jseaes.2010.03.009. 9. Chou, M.-Y.; Cherniak, I.; Lin, C.C.H.; Pedatella, N.M. The persistent ionospheric responses over Japan after the impact of the 2011 Tohoku earthquake. Space Weather 2020, 18, e2019SW002302. https://doi.org/10.1029/2019SW002302. 10. Shinbori, A.; Otsuka, Y.; Sori, T.; Nishioka, M.; Perwitasari, S.; Tsuda, T.; Nishitani, N. Electromagnetic conjugacy of ionospheric disturbances after the 2022 Hunga Tonga-Hunga Ha’apai volcanic eruption as seen in GNSS-TEC and SuperDARN Hokkaido pair of radars observations. Earth Planets Space 2022, 74, 106 https://doi.org/10.1186/s40623-022-01665-8. 11. Calais, E.; Minster, J.B.; Hofton, M.; Hedlin, M.A.H. Ionospheric signature of surface mine blasts from Global Positioning System measurements. Geophys. J. Int. 1998, 132, 191–202. https://doi.org/10.1046/j.1365-246x.1998.00438.x. 12. Heki, K.; Otsuka, Y.; Choosakul, N.; Hemmakorn, N.; Komolmis, T.; Maruyama, T. Detection of ruptures of Andaman fault segments in the 2004 great Sumatra earthquake with coseismic ionospheric disturbances. J. Geophys. Res. Solid Earth 2006, 111, B09313. https://doi.org/10.1029/2005JB004202. 13. Jin, S.G.; Jin, R.; Li, J.H. Pattern and evolution of seismo-ionospheric disturbances following the 2011 Tohoku earthquakes from GPS observations. J. Geophys. Res. Space Phys. 2014, 119, 7914–7927. https://doi.org/10.1002/2014JA019825. 14. Jin, S.G.; Occhipinti, G.; Jin, R. GNSS ionospheric seismology: Recent observations evidences and characteristics. Earth-Sci. Rev. 2015, 147, 54–64. https://doi.org/10.1016/j.earscirev.2015.05.003. 15. Piscini, A.; Marchetti, D.; De Santis, A. Multi-Parametric Climatological Analysis Associated with Global Significant Volcanic Eruptions During 2002–2017. Pure Appl. Geophys. 2019, 176, 3629–3647. https://doi.org/10.1007/s00024-019-02147-x. 16. Dučić, V.; Artru, J.; Lognonné, P. Ionospheric remote sensing of the Denali Earthquake Rayleigh surface waves. Geophys. Res. Lett. 2003, 30, 1951. https://doi.org/10.1029/2003GL017812. 17. Occhipinti, G.; Rolland, L.; Lognonné, P.; Watada, S. From Sumatra 2004 to TohokuOki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes. J. Geophys. Res. Space Phys. 2013, 118, 3626–3636. https://doi.org/10.1002/jgra.50322. 18. Rolland, L.M.; Lognonné, P.; Astafyeva, E.; Kherani, E.A.; Kobayashi, N.; Mann, M.; Munekane, H. The resonant response of the ionosphere imaged after the 2011 off the Pacific coast of Tohoku Earthquake. Earth Planets Space 2011, 63, 853–857. https://doi.org/10.5047/eps.2011.06.020. 19. Su, Y.C.; Liu, J.Y.; Chen, S.P.; Tsai, H.F.; Chen, M.Q. Temporal and spatial precursors in ionospheric total electron content of the 16 October 1999 Mw7.1 Hector Mine earthquake. J. Geophys. Res. Space Phys. 2013, 118, 6511–6517. https://doi.org/10.1002/jgra.50586. 20. Yang, H.; Monte Moreno, E.; Hernández-Pajares, M. ADDTID: An alternative tool for studying earthquake/tsunami signatures in the ionosphere. Case of the 2011 Tohoku earthquake. Remote Sens. 2019, 11, 1894. https://doi.org/10.3390/rs11161894. 21. Zhao, B.; Yu, T.; Wang, M.; Wan, W.; Lei, J.; Liu, L.; Ning, B. Is an unusual large enhancement of ionospheric electron density linked with the 2008 great Wenchuan earthquake? J. Geophys. Res. 2008, 113, 304. https://doi.org/10.1029/2008JA013613. 22. Ke, F.; Wang, Y.; Wang, X.; Qian, H.; Shi, C. Statistical analysis of seismo-ionospheric anomalies related to Ms > 5.0 earthquakes in China by GPS TEC. J. Seismol. 2016, 20, 137–149. https://doi.org/10.1007/s10950-015-9516-x. 23. Le, H.; Liu, J.Y.; Liu, L. A statistical analysis of ionospheric anomalies before 736 M6.0+ earthquakes during 2002–2010. J. Geophys. Res. 2011, 116, A02303. https://doi.org/10.1029/2010ja015781. 24. Li, M.; Parrot, M. Statistical analysis of an ionospheric parameter as a base for earthquake prediction. J. Geophys. Res. Space Phys. 2013, 118, 3731–3739. https://doi.org/10.1002/jgra.50313. 25. Píša, D.; Parrot, M.; Santolík, O. Ionospheric density variations recorded before the 2010 Mw 8.8 earthquake in Chile. J. Geophys. Res. Space Phys. 2011, 116, A08309. https://doi.org/10.1029/2011JA016611. 26. Anagnostopoulos, G.C.; Vassiliadis, E.; Pulinets, S. Characteristics of flux time profiles, temporal evolution, and spatial distribution of radiation belt electron precipitation bursts in the upper ionosphere before great and giant earthquakes. Ann. Geophys. 2012, 55, 21–36. https://doi.org/10.4401/ag-5365. 27. Ryu, K.; Lee, E.; Chae, J.S.; Parrot, M.; Pulinets, S. Seismo-ionospheric coupling appearing as equatorial electron density enhancements observed via DEMETER electron density measurements. J. Geophys. Res. Space Phys. 2014, 119, 8524–8542. https://doi.org/10.1002/2014JA020284. 28. Zhang, X.; Shen, X.; Parrot, M.; Zeren, Z.; Ouyang, X.; Liu, J.; Qian, J.; Zhao, S.; Miao, Y. Phenomena of electrostatic perturbations before strong earthquakes (2005–2010) observed on DEMETER. Nat. Hazards Earth Syst. Sci. 2012, 12, 75–83. https://doi.org/10.5194/nhess-12-75-2012. 29. Akhoondzadeh, M.; Parrot, M.; Saradjian, M.R. Electron and ion density variations before strong earthquakes (M > 6.0) using DEMETER and GPS data. Nat. Hazards Earth Syst. Sci. 2010, 10, 7–18. https://doi.org/10.5194/nhess-10-7-2010. 30. Ibanga, J.I.; Akpan, A.E.; George, N.J.; Ekanem, A.M.; George, A. Unusual ionospheric variations before the strong Auckland Islands, New Zealand earthquake of 30th September, 2007. NRIAG J. Astron. Geophys. 2018, 7, 149–154. https://doi.org/10.1016/j.nrjag.2017.12.007. 31. Diego, P.; Coco, I.; Bertello, I.; Candidi, M.; Ubertini, P. Ionospheric Plasma Density Measurements by Swarm Langmuir Probes: Limitations and possible Corrections. Ann. Geophys. Discuss. 2019, https://doi.org/10.5194/angeo-2019-136. Remote Sens. 2023, 15, x FOR PEER REVIEW 36 of 37 32. Luo, X.; Xiong, C.; Gu, S.; Lou, Y.; Stolle, C.; Wan, X.; Liu, K.; Song, W. Geomagnetically conjugate observations of equatorial plasma irregularities from Swarm constellation and ground-based GPS stations. J. Geophys. Res. Space Phys. 2019, 124, 3650–3665. https://doi.org/10.1029/2019JA026515. 33. Liu, C.; Guan, Y.; Zheng, X.; Zhang, A.; Diego, P.; Sun, Y. The technology of space plasma in-situ measurement on the China Seismo-Electromagnetic Satellite. Sci. China Technol. Sci. 2019, 62, 829–838. https://doi.org/10.1007/s11431-018-9345-8. 34. Zhu, K.; Li, K.; Fan, M.; Chi, C.; Yu, Z. Precursor analysis associated with the Ecuador earthquake using Swarm A and C satellite magnetic data based on PCA. IEEE Access 2019, 7, 93927–93936. https://doi.org/10.1109/ACCESS.2019.2928015. 35. Zhu, K.; Fan, M.; He, X.; Marchetti, D.; Li, K.; Yu, Z.; Chi, C.; Sun, H.; Cheng, Y. Analysis of Swarm Satellite Magnetic Field Data Before the 2016 Ecuador (Mw = 7.8) Earthquake Based on Non-negative Matrix Factorization. Front. Earth Sci. 2021, 9, 621976. https://doi.org/10.3389/feart.2021.621976. 36. Akhoondzadeh, M. Seismo-Magnetic Field Anomalies Detection Using Swarm Satellites (Alpha, Bravo and Charlie). Int. Arch. Photogramm. Remote Sens, Spat. Inf. Sci. 2019, XLII-4/W18, 45–49. https://doi.org/10.5194/isprs-archives-XLII-4-W18-45-2019. 37. Marchetti, D.; De Santis, A.; Jin, S.; Campuzano, S.A.; Cianchini, G.; Piscini, A. Co-Seismic Magnetic Field Perturbations Detected by Swarm Three-Satellite Constellation. Remote Sens. 2020, 12, 1166. https://doi.org/10.3390/rs12071166. 38. Marchetti, D.; De Santis, A.; D’Arcangelo, S.; Poggio, F.; Jin, S.; Piscini, A.; Campuzano, S.A. Magnetic Field and Electron Density Anomalies from Swarm Satellites Preceding the Major Earthquakes of the 2016–2017 Amatrice-Norcia (Central Italy) Seismic Sequence. Pure Appl. Geophys 2020, 177, 305–319. https://doi.org/10.1007/s00024-019-02138-y. 39. Ghamry, E.; Mohamed, E.K.; Abdalzaher, M.S.; Elwekeil, M.; Marchetti, D.; De Santis, A.; Hegy, M.; Yoshikawa, A.; Fathy, A. Integrating Pre-Earthquake Signatures from Different Precursor Tools. IEEE Access 2021, 9, 33268–33283. https://doi.org/10.1109/ACCESS.2021.3060348. 40. Stanica, D.A.; Stanica, D.; Błęcki, J.; Ernst, T.; Jozwiak, W.; Słominski, J. Pre-seismic geomagnetic and ionosphere signatures related to the Mw5.7 earthquake occurred in Vrancea zone on September 24, 2016. Acta Geophys. 2018, 66, 167–177. https://doi.org/10.1007/s11600-018-0115-4. 41. Marchetti, D.; Akhoondzadeh, M. Analysis of Swarm satellites data showing seismo-ionospheric anomalies around the time of the strong Mexico (Mw = 8.2) earthquake of 08 September 2017. Adv. Space Res. 2018, 62, 614–623. https://doi.org/10.1016/j.asr.2018.04.043. 42. Marchetti, D.; De Santis, A.; Shen, X.; Campuzano, S.; Perrone, L.; Piscini, A.; Giovambattista, R.D.; Jin, S.; Ippolito, A.; Cianchini, G.; et al. Possible Lithosphere-Atmosphere-Ionosphere Coupling effects prior to the 2018 Mw = 7.5 Indonesia earthquake from seismic, atmospheric and ionospheric data. J. Asian Earth Sci. 2020, 188, 104097. https://doi.org/10.1016/j.jseaes.2019.104097. 43. De Santis, A.; Marchetti, D.; Spogli, L.; Cianchini, G.; Pavón-Carrasco, F.J.; Franceschi, G.D.; Di Giovambattista, R.; Perrone, L.; Qamili, E.; Cesaroni, C.; et al. Magnetic Field and Electron Density Data Analysis from Swarm Satellites Searching for Ionospheric Effects by Great Earthquakes: 12 Case Studies from 2014 to 2016. Atmosphere 2019, 10, 371. https://doi.org/10.3390/atmos10070371. 44. Huang, H.; Yan, R.; Liu, D.; Xu, S.; Lin, J.; Guo, F.; Huang, J.; Zeren, Z.; Shen, X. The variations of plasma density recorded by CSES-1 satellite possibly related to Mexico Ms 7.1 earthquake on 8th September 2021. Nat. Hazard Res. 2022, in press. https://doi.org/10.1016/j.nhres.2021.12.002. 45. Yan, R.; Shen, X.H.; Huang, J.P.; Wang, Q.; Chu, W.; Liu, D.P.; Yang, Y.Y.; Lu, H.X.; Xu, S. Examples of unusual ionospheric observations by the CSES prior to earthquakes. Earth Planet. Phys. 2018, 2, 515–526. http://doi.org/10.26464/epp2018050. 46. Xiong, P.; Long, C.; Zhou, H.; Battiston, R.; De Santis, A.; Ouzounov, D.; Zhang, X.; Shen, X. Pre-Earthquake Ionospheric Perturbation Identification Using CSES Data via Transfer Learning. Front. Environ. Sci. 2021, 9, 779255. https://doi.org/10.3389/fenvs.2021.779255. 47. De Santis, A.; Marchetti, D.; Perrone, L.; Campuzano, S.A.; Cianchini, G.; Cesaroni, C.; Di Mauro, D.; Orlando, M.; Piscini, A.; Sabbagh, D.; et al. Statistical correlation analysis of strong earthquakes and ionospheric electron density anomalies as observed by CSES-01. Il Nuovo C. 2021, 44, 119. https://doi.org/10.1393/ncc/i2021-21119-1. 48. De Santis, A.; Marchetti, D.; Pavón-Carrasco, F.J.; Cianchini, G.; Perrone, L.; Abbattista, C.; Alfonsi, L.; Amoruso, L.; Campuzano, S.A.; Carbone, M.; et al. Precursory worldwide signatures of earthquake occurrences on Swarm satellite data. Sci. Rep. 2019, 9, 20287. https://doi.org/10.1038/s41598-019-56599-1. 49. He, Y.; Zhao, X.; Yang, D.; Wu, Y.; Li, Q. A study to investigate the relationship between ionospheric disturbance and seismic activity based on Swarm satellite data. Phys. Earth Planet. Inter. 2021, 323, 106826. https://doi.org/10.1016/j.pepi.2021.106826. 50. Sunil, A.S.; Bagiya, M.S.; Reddy, C.D.; Kumar, M.; Sai Ramesh, D. Post-seismic ionospheric response to the 11 April 2012 East Indian Ocean doublet earthquake. Earth Planets Space 2015, 67, 37. https://doi.org/10.1186/s40623-015-0200-8. 51. Balasis, G.; Daglis, I.A.; Georgiou, M.; Papadimitriou, C.; Haagmans, R. Magnetospheric ULF wave studies in the frame of Swarm mission: A time-frequency analysis tool for automated detection of pulsations in magnetic and electric field observations. Earth Planet. Space 2013, 65, 18. https://doi.org/10.5047/eps.2013.10.003. 52. Balasis, G.; Papadimitriou, C.; Daglis, I.A.; Pilipenko, V. ULF wave power features in the topside ionosphere revealed by Swarm observations. Geophys. Res. Lett. 2015, 42, 6922–6930. https://doi.org/10.1002/2015GL065424. 53. Balasis, G.; Papadimitriou, C.; Boutsi, A.Z. Ionospheric response to solar and interplanetary disturbances: A Swarm perspective. Philos. Trans. R. Soc. A 2019, 377, 20180098. https://doi.org/10.1098/rsta.2018.0098. 54. Jarmołowski, W.; Belehaki, A.; Hernández Pajares, M.; Schmidt, M.; Goss, A.; Wielgosz, P.; Yang, H.; Krypiak-Gregorczyk, A.; Tsagouri, I.; Paouris, E.; et al. Combining Swarm Langmuir probe observations, LEO-POD-based and ground-based GNSS Remote Sens. 2023, 15, x FOR PEER REVIEW 37 of 37 receivers and ionosondes for prompt detection of ionospheric earthquake and tsunami signatures: Case study of 2015 ChileIllapel event. J. Space Weather Space Clim. 2021, 11, 58. https://doi.org/10.1051/swsc/2021042. 55. Dutilleul, P.; Johnson, C.W.; Bürgmann, R.; Wan, Y.; Shen, Z.-K. Multifrequential periodogram analysis of earthquake occurrence: An alternative approach to the Schuster spectrum, with two examples in central California. J. Geophys. Res. Solid. Earth 2015, 120, 8494–8515. https://doi.org/10.1002/2015JB012467. 56. Zhan, Z.; Shearer, P.M. Possible seasonality in large deep-focus earthquakes. Geophys. Res. Lett. 2015, 42, 7366–7373. https://doi.org/10.1002/2015GL065088. 57. Jin, Y.; Xiong, C.; Clausen, L.; Spicher, A.; Kotova, D.; Brask, S.; Kervalishvili, G.; Stolle, C.; Miloch, W. Ionospheric plasma irregularities based on in situ measurements from the Swarm satellites. J. Geophys. Res. Space Phys. 2020, 124, e2020JA028103. https://doi.org/10.1029/2020JA028103. 58. Kepkar, A.; Arras, C.; Wickert, J.; Schuh, H.; Alizadeh, M.; Tsai, L.-C. Occurrence climatology of equatorial plasma bubbles derived using FormoSat-3/COSMIC GPS radio occultation data. Ann. Geophys. 2020, 38, 611–623. https://doi.org/10.5194/angeo38-611-2020. 59. Zakharenkova, I.; Astafyeva, E.; Cherniak, I. GPS and in situ Swarm observations of the equatorial plasma density irregularities in the topside ionosphere. Earth Planet. Space 2016, 68, 120. https://doi.org/10.1186/s40623-016-0490-5. 60. Harris, F.J. On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform. Proc. IEEE 1978, 66, 66–67. 61. Stankovic, L. 2015. Digital Signal Processing with Selected Topics; An Amazon.com Company, North Charleston, South Carolina, USA: 2015; p. 820. 62. Perevalova, N.P.; Sankov, V.A.; Astafyeva, E.I.; Zhupityaeva, А.S. Threshold magnitude for ionospheric TEC response to earthquakes. J. Atmos. Sol. Terr. Phys. 2014, 108, 77–90. https://doi.org/10.1016/j.jastp.2013.12.014. 63. Astafyeva, E.; Shults, K. Ionospheric GNSS imagery of seismic source: Possibilities, difficulties, and challenges. J. Geophys. Res. Space Phys. 2019, 124, 534–543. https://doi.org/10.1029/2018JA026107. 64. Goto, S.I.; Uchida, R.; Igarashi, K.; Chen, C.-H.; Kao, M.; Umeno, K. Preseismic ionospheric anomalies detected before the 2016 Taiwan earthquake. J. Geophys. Res. Space Phys. 2019, 124, 9239–9252. https://doi.org/10.1029/2019JA026640. 65. Rajesh, R.; Catherine, J.K. The first evidence for the detection of CIDs masked by equatorial plasma bubbles from GPS-TEC data. J. Geophys. Res. Space Phys. 2022, 127, e2021JA029798. https://doi.org/10.1029/2021JA029798. 66. Kawada, Y.; Nagahama, H. Cumulative Benioff strain-release, modified Omori’s law and transient behaviour of rocks. Tectonophysics 2006, 424, 157–166. https://doi.org/10.1016/j.tecto.2006.03.032. 67. Rundle, J.B.; Klein, W.; Turcotte, D.L.; Malamud, B.D. Precursory seismic activation and critical-point phenomena. Pure Appl. Geophys. 2000, 157, 2165–2182. https://doi.org/10.1007/PL00001079. 68. Varnes, D.J. Predicting earthquakes by analyzing accelerating precursory seismic activity. Pure Appl. Geophys. 1989, 130, 661– 686. https://doi.org/10.1007/BF00881603. 69. Park, J.; Noja, M.; Stolle, C.; Lühr, H. The Ionospheric Bubble Index deduced from magnetic field and plasma observations onboard Swarm. Earth Planet. Space 2013, 65, 13. https://doi.org/10.5047/eps.2013.08.005. 70. Kuo, C.L.; Huba, J.D.; Joyce, G.; Lee, L.C. Ionosphere plasma bubbles and density variations induced by pre-earthquake rock currents and associated surface charges. J. Geophys. Res. 2011, 116, A10317. https://doi.org/10.1016/j.dib.2019.104513. Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.