Co-Seismic Magnetic Field Perturbations Detected by Swarm Three-Satellite Constellation

The first 5.3 years of magnetic data from three Swarm satellites have been systematically analyzed, and possible co-seismic magnetic disturbances in the ionosphere were investigated just a few minutes after the occurrence of large earthquakes. We preferred to limit the investigation to a subset of earthquakes selected in function of depth and magnitude. After a systematic inspection of the available data around (in time and space) the seismic events, we found 12 Swarm satellite tracks with co-seismic disturbances possibly produced by ten earthquakes from Mw5.6 to Mw6.9. The distance of the satellite to the earthquake epicenter corresponds to the measured distance-time arrival of the disturbance from the surface to the ionosphere, confirming that the identified disturbances are most likely produced by the seismic events. Secondly, we found a good agreement with a model that combined a propagation of the disturbance to the F2 ionospheric layer with an acoustic gravity wave at a velocity of about (2.2 ± 0.3) km/s and a second faster phenomenon that transmits the disturbance from F2 layer to the Swarm satellite with a velocity of about (16 ± 3) km/s as an electromagnetic scattering propagation.


Introduction
This paper aimed to search for perturbations of the geomagnetic field in the ionosphere possibly produced by the occurrence of an earthquake in the lithosphere. The possibility that the ionosphere could be perturbed just a few minutes after the occurrence of a medium (Mw6+) and large (Mw7+) earthquakes is well known. For example, Yeh and Liu, 1974 [1] presented a model for the generation, propagation and response of the ionosphere to acoustic gravity waves. They also provided experimental data by supporting their model. Yamamoto, 1982 [2] theorized a model to describe the acoustic gravity waves induced by submarine earthquakes. Submarine earthquakes deserve special attention due to their potential to generate tsunamis. Lognonne et al., 2006 [3] observed the impact of the acoustic waves by the Global Positioning System (GPS) network 660-670 s (~11 min) after the arrival time of Rayleigh waves at the ground for the Mw7.9 November 3rd, 2002, Alaska earthquake. Hao et al., 2012 [4] reported not only ionospheric Doppler caused by the Mw9.0 Tohoku 2011 earthquake, but deepest events that could produce a disturbance in the ionosphere because of their significant released energy. The number of the selected events was 1830 out of 2437 downloaded earthquakes, represented by a red star/dot in Figure 1. In particular, in Figure 1b it is possible to note the exclusion of the deepest events. The selected earthquakes are a sub-set in which it is more reliable to search for effects induced in the ionosphere, but future works can also investigate other events or subsets.

Swarm Satellite Magnetic Data
In this paper, the data from the ESA Swarm mission have been investigated. The main goal of the mission was to measure and monitor the Earth's magnetic field with global coverage and better understand the electromagnetic phenomena on our planet [17]. Swarm is a three-identical satellite constellation successfully launched on 22 November 2013. The mission is still in orbit, and the present objective is to extend the mission for an entire solar cycle of 11 years. Most of the instruments continue to work in nominal conditions. Currently, more than five years of data are available.
The three satellites follow quasi-polar orbits: two of them (called Alpha and Charlie) fly close to each other at a lower orbit (around 450 km) and the other one, called Bravo, currently has a separation of about 90° longitude, and it is placed at a higher orbit, around 510 km. In the first months of the mission, all the satellites orbited close to each other.
Each satellite is equipped with the following payloads: Vector Field Magnetometer (VFM), Absolute Scalar Magnetometer (ASM), electric field instrument (EFI), accelerometers, GNSS receivers and laser retroreflectors. The VFM and ASM measure the orientation and intensity of the magnetic field. These sensors are placed on a 4-m length boom to avoid possible interferences with the satellite platform. The VFM sensor, placed half-way along the boom, is a fluxgate magnetometer with threeaxis compact spherical coil with a three-axis compact detector coil inside. The magnetic field intensity calibration of the VFM is based on the data acquired by the ASM instrument placed at the end of the boom. For the Charlie satellite, the calibration is set on the data of Alpha ASM, as its ASM has not functioned since November 2013. The orientation of the VFM was precisely monitored with three star-tracker cameras placed on the same optical bank of the magnetometer.
The Swarm mission has already provided vast results; among them, a map of the lithospheric magnetic field with a resolution up to 250 km [18] and the oceanic tidal effect on geomagnetic field and electromagnetic induction generated by oceanic circulation (e.g., [19][20][21]). Furthermore, some of the authors of this paper found a series of anomalies in Swarm magnetic field and electron density data that could be related to the preparation phases of some earthquakes around the world (e.g., Mw7.8 Nepal 2015 [22], Mw8.2 Mexico 2017 [23], 2018, Mw6.5 Italy 2016 [24], Mw7.3 Iran 2017 [25], Mw7.5 Indonesia 2018 [26]). A recent paper by De Santis et al., 2019 [27] provided for the first time a

Swarm Satellite Magnetic Data
In this paper, the data from the ESA Swarm mission have been investigated. The main goal of the mission was to measure and monitor the Earth's magnetic field with global coverage and better understand the electromagnetic phenomena on our planet [17]. Swarm is a three-identical satellite constellation successfully launched on 22 November 2013. The mission is still in orbit, and the present objective is to extend the mission for an entire solar cycle of 11 years. Most of the instruments continue to work in nominal conditions. Currently, more than five years of data are available.
The three satellites follow quasi-polar orbits: two of them (called Alpha and Charlie) fly close to each other at a lower orbit (around 450 km) and the other one, called Bravo, currently has a separation of about 90 • longitude, and it is placed at a higher orbit, around 510 km. In the first months of the mission, all the satellites orbited close to each other.
Each satellite is equipped with the following payloads: Vector Field Magnetometer (VFM), Absolute Scalar Magnetometer (ASM), electric field instrument (EFI), accelerometers, GNSS receivers and laser retroreflectors. The VFM and ASM measure the orientation and intensity of the magnetic field. These sensors are placed on a 4-m length boom to avoid possible interferences with the satellite platform. The VFM sensor, placed half-way along the boom, is a fluxgate magnetometer with three-axis compact spherical coil with a three-axis compact detector coil inside. The magnetic field intensity calibration of the VFM is based on the data acquired by the ASM instrument placed at the end of the boom. For the Charlie satellite, the calibration is set on the data of Alpha ASM, as its ASM has not functioned since November 2013. The orientation of the VFM was precisely monitored with three star-tracker cameras placed on the same optical bank of the magnetometer.
The Swarm mission has already provided vast results; among them, a map of the lithospheric magnetic field with a resolution up to 250 km [18] and the oceanic tidal effect on geomagnetic field and electromagnetic induction generated by oceanic circulation (e.g., [19][20][21]). Furthermore, some of the authors of this paper found a series of anomalies in Swarm magnetic field and electron density data that could be related to the preparation phases of some earthquakes around the world (e.g., Mw7.8 Nepal 2015 [22] [27] provided for the first time a worldwide statistical correlation of Swarm electron density and magnetic anomalies with earthquake occurrences, finding Remote Sens. 2020, 12, 1166 4 of 21 some distinct concentrations of anomalies that preceded the earthquakes from some days to several months; additionally, a linear dependency was established between the logarithm of the anticipation time of the anomaly and the earthquake magnitude, confirming for satellite data the empirical Rikitake law [28], which was already introduced for ground precursors. The Swarm magnetic data have been downloaded from the free accessible FTP repository (swarm-diss.eo.esa.int) in Level 1B as CDF files with the filename SW_OPER_MAGx_LR_1B_[ . . . ], where x stands for A, B or C satellite and the filename is completed with the starting and ending timestamp of the file and the data version (408 for most of this work and 505 for the most recent data). MAG is the label that identifies the magnetic product and LR stands for Low Resolution (1 Hz) that is a product re-sampled from the original one (HR, i.e., High Resolution: 50 Hz) at the GPS seconds.

Methods to Analyze the Magnetic Satellite Data
A series of MATLAB ® routines have been developed to analyze the Swarm Magnetic data. The first script reads the earthquake catalogue and selects the events in agreement with the magnitude-depth requirements described in Section 2.1 (all earthquakes with Mw ≥ 5.5 and depth according to a linearly varying threshold from 20 km depth for Mw = 5.5 and Mw = 7.0 at 150 km depth). For each event, the algorithm selects the magnetic satellite data on the same date. The polar regions above a 70 • geographical latitude are cut, as they are not suitable for this study due to the typical perturbation of the auroral activity. The goal of this work was to analyze the anomalous magnetic activity associated with earthquakes; thus, other magnetic perturbations, such as the geomagnetic storms, must first be rejected. After that, the script calculates the spatial distance of each sample from the epicentre and the time difference between the acquisition time and the earthquake origin time. If there are some samples close in space-time to the seismic event (distance ≤ 150 km and time difference from 10 min before until 20 min after the earthquake), the script selects and runs an automatic analysis of the close track.
Two types of automatic analyses have been performed: the first one is called Method 1, and it is the same proposed in previous works on Swarm data related to earthquakes (e.g., [22,24]. This technique consists of computing the first derivative of the data along the track approximated as the first difference divided for the time interval between two consecutive samples. After this step, a cubic spline is removed and the residuals of X, Y, Z (geomagnetic components) and F (intensity) are plotted with a fixed maximum range of ± 1 nT/s. The second method (so-called Method 2) has been developed to avoid a possible over-fitting problem of the first method, i.e., that the spline could remove the eventual perturbation induced by the supposed gravity wave. It consists of the subtraction of the International Geomagnetic Reference Field (IGRF-12) model [29] from the original data and the plotting of the residuals in X, Y, Z and F. The IGRF-12 is a global Earth's magnetic field model, which is expected to be sufficient to highlight the studied disturbances. A five-degree polynomial has been further applied to remove the residual trend of each track.
The output figures of both methods are presented in the Supplementary Materials. Each figure presents a map of the geographical region crossed by the satellite, the projection of the satellite track on the ground (brown for normal conditions and red for some abnormal conditions on the satellite) and the earthquake epicenter (green star). The title of the figure provides the satellites, the time in UTC and Local Time and the geomagnetic indices Dst and ap during the acquisition time. The residuals of X, Y, Z and F are reported together with the Fast Fourier Transform analysis of the signal that helps to select the best cut-off frequency for the following signal processing (for the second method, the mean of the signal is subtracted before calculating the FFT). The figures are grouped in tables that refer to each seismic event.
All the satellite data of the selected track are stored in a text file for further analysis. The produced material is visually investigated and the tracks that present some possible features linkable to the seismic events are then re-analyzed and plotted again. For this final analysis (presented in Figures 3-14), only the Y component is considered: after the subtraction of IGRF-12 model, a high pass filter at 0.25 Hz Remote Sens. 2020, 12, 1166 5 of 21 was applied and both quantities were then plotted. The frequency cut-off was chosen after multiple tests as the best value to extract the anomalous signal possibly related to the seismic event. Taking into account that the sampling rate of the data was 1 Hz, the available band was 0.25 Hz-0.50 Hz (following Nyquist theorem).

Results of the Analysis
From the launch of the mission on November 22 nd , 2013 until April 2019, the Swarm satellite data have been analyzed to search when one or more satellites passed above the epicentral area during an earthquake occurrence. To be conservative, a simultaneous passage was defined by a circular area of 150 km around the epicentre and a time window that spanned 10 min before the event until 20 min after. With such criteria, 37 satellite passages (tracks) were selected. For each case, method 1 and 2 were applied, searching for possible "interesting" features in the residual of at least one method. We finally noticed that even if the algorithm selected the space-time tracks closer to the epicenter, in some cases, the satellites passed too early or too late above the epicentral area, losing the chance to detect a co-seismic signal. After visually inspecting the residual plots, 12 tracks were considered the most likely candidate containing a co-seismic effect. Figure 2 reports the distribution map of the earthquakes and the ionospheric disturbances possibly associated with the seismic events.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 20 pass filter at 0.25 Hz was applied and both quantities were then plotted. The frequency cut-off was chosen after multiple tests as the best value to extract the anomalous signal possibly related to the seismic event. Taking into account that the sampling rate of the data was 1 Hz, the available band was 0.25 Hz-0.50 Hz (following Nyquist theorem).

Results of the Analysis
From the launch of the mission on November 22 nd , 2013 until April 2019, the Swarm satellite data have been analyzed to search when one or more satellites passed above the epicentral area during an earthquake occurrence. To be conservative, a simultaneous passage was defined by a circular area of 150 km around the epicentre and a time window that spanned 10 min before the event until 20 min after. With such criteria, 37 satellite passages (tracks) were selected. For each case, method 1 and 2 were applied, searching for possible "interesting" features in the residual of at least one method. We finally noticed that even if the algorithm selected the space-time tracks closer to the epicenter, in some cases, the satellites passed too early or too late above the epicentral area, losing the chance to detect a co-seismic signal. After visually inspecting the residual plots, 12 tracks were considered the most likely candidate containing a co-seismic effect. Figure 2 reports the distribution map of the earthquakes and the ionospheric disturbances possibly associated with the seismic events.  Table 1 reports the list of the earthquakes that present a possible magnetic disturbance in the ionosphere with the Swarm satellite characteristics in the moment of disturbance detection. The measured time is that after the origin time of the earthquake in which the anomaly starts to be detected by the Swarm satellite. The table also reports the height in km of the ionospheric F2 electron density peak (hmF2) estimated by the IRI-2016 model [30] and the air temperature at 2 m retrieved from the European Centre for Medium-range Weather Forecast (ECMWF) ERA-Interim archive, which will be used for the discussion and interpretation of the results. Table 1. List of the earthquakes that created a potential magnetic disturbance in the ionosphere as detected from the Swarm satellite. The orbital characteristics are reported (satellite, direction of the semi-orbit: D for downward orbit, U for upward orbit, altitude, latitude and longitude) in the moment of disturbance detection. The height of the F2 ionospheric layer (hmF2) and the surface air temperature are also provided.  Table 1 reports the list of the earthquakes that present a possible magnetic disturbance in the ionosphere with the Swarm satellite characteristics in the moment of disturbance detection. The measured time is that after the origin time of the earthquake in which the anomaly starts to be detected by the Swarm satellite. The table also reports the height in km of the ionospheric F2 electron density peak (hmF2) estimated by the IRI-2016 model [30] and the air temperature at 2 m retrieved from the European Centre for Medium-range Weather Forecast (ECMWF) ERA-Interim archive, which will be used for the discussion and interpretation of the results. Table 1. List of the earthquakes that created a potential magnetic disturbance in the ionosphere as detected from the Swarm satellite. The orbital characteristics are reported (satellite, direction of the semi-orbit: D for downward orbit, U for upward orbit, altitude, latitude and longitude) in the moment of disturbance detection. The height of the F2 ionospheric layer (hmF2) and the surface air temperature are also provided.

Earthquake Characteristics
Satellite Characteristics The Swarm magnetic data during the earthquake occurrence are represented following the same order of Table 1, i.e., a chronological one, from Figures 3-13. Each figure is horizontally aligned with the origin time of the earthquake. In the upper part is shown in red the residual of the Y component of the magnetic field with respect to the IGRF-12 model. The Y component was selected as it is the most sensible for the internal perturbations ( [22,31]). In the lower graph, a 0.25 Hz high-pass filter was applied, and the filtered data are shown in blue. Here, we did not subtract any polynomial, as the filter was sufficient to remove the residual trend. A black vertical line underlines the position where the satellite crosses the earthquake latitude if it is after the earthquake occurrence. An orange line represents the latitude of the satellite. Note that only the data inside −70 • and +70 • geographical latitude are shown. More details are provided in the Supplementary Materials, where Figures S1-S10 show the residuals of the other components as well as the total intensity of the magnetic field, as deduced using method 1 and method 2. Figure S11 also shows an example of a track without any evidence of possible co-seismic magnetic disturbance. The reasons for the lack of disturbance in the latter case could be several: the earthquake happened in the Atlantic ocean ridge and the ocean water depth at the epicentre is of several kilometres (about 3.4 km), thus possibly attenuating the mechanical effect (also considering the low magnitude of the event: Mw 5.5); another reason could be an intrinsic disturbance in the magnetic signal all over the track that could cover the detection of any small co-seismic effect. On the other hand, in the Figures S1-S10, the magnetic disturbance cannot be explainable merely by geomagnetic activity, which was low at the acquisition time; nevertheless, it could be induced by some other geophysical source (such as electric and magnetic fields produced by ocean currents-Grayver et al., 2016-or tectonic stress along the ocean drift). mission that we analyzed.
An earthquake happened offshore Chile on 2015-11-11 at a depth of 10 km between the Nazca and South-American plates. In this case, a perturbation was visible just north of the epicentre and we cannot exclude that this could be an electromagnetic phenomenon coming directly from the fault. There was a certain symmetry of the signal with respect to the epicentre latitude. The ionospheric disturbance was possibly associated with the high magnitude earthquake (Mw 6.9) considered in this case ( Figure 6).   Figure S1). mission that we analyzed.
An earthquake happened offshore Chile on 2015-11-11 at a depth of 10 km between the Nazca and South-American plates. In this case, a perturbation was visible just north of the epicentre and we cannot exclude that this could be an electromagnetic phenomenon coming directly from the fault. There was a certain symmetry of the signal with respect to the epicentre latitude. The ionospheric disturbance was possibly associated with the high magnitude earthquake (Mw 6.9) considered in this case ( Figure 6).  The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. This anomaly has been selected as it is the wider after the earthquake occurrence and it is the corresponding anomaly of the previous track (acquired after few seconds). In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S1). The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. This anomaly has been selected as it is the wider after the earthquake occurrence and it is the corresponding anomaly of the previous track (acquired after few seconds). In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S1). The vertical green line is the candidate point where the earthquakeinduced phenomena are supposed to start. We selected this anomaly because it is the first after the earthquake occurred, while those occurred after 7 min are localized at high latitudes, and thus more likely due to magnetic auroral activity. The data after +70° of geographical latitude are cut. The satellite passed just before the earthquake over the latitude of the epicentre. After little more than 10 min, the satellite overpassed the 70° north geographical latitude, and thus, the data are not shown. In the Supplementary Materials, more details and other available Swarm magnetic data are represented ( Figure S2).

Figure 6.
Residual of the Y component of the Swarm Alpha magnetic field after the occurrence of the Mw 6.9 2015-11-11 earthquake. The vertical black line in the lower plot represents the point where the satellite crosses the earthquake latitude. The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. We selected this anomaly because it is the first on the track after the earthquake. Moreover, we noticed that its shape is different from the other anomalies, as the amplitude is almost constant (peak-to-peak about 1nT) during its duration (little The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. We selected this anomaly because it is the first after the earthquake occurred, while those occurred after 7 min are localized at high latitudes, and thus more likely due to magnetic auroral activity. The data after +70 • of geographical latitude are cut. The satellite passed just before the earthquake over the latitude of the epicentre. After little more than 10 min, the satellite overpassed the 70 • north geographical latitude, and thus, the data are not shown. In the Supplementary Materials, more details and other available Swarm magnetic data are represented ( Figure S2). Mw 5.9 2014-02-18 earthquake. The vertical black line in the lower plot represents the point where the satellite crosses the earthquake latitude. The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. This anomaly has been selected as it is the wider after the earthquake occurrence and it is the corresponding anomaly of the previous track (acquired after few seconds). In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S1). The vertical green line is the candidate point where the earthquakeinduced phenomena are supposed to start. We selected this anomaly because it is the first after the earthquake occurred, while those occurred after 7 min are localized at high latitudes, and thus more likely due to magnetic auroral activity. The data after +70° of geographical latitude are cut. The satellite passed just before the earthquake over the latitude of the epicentre. After little more than 10 min, the satellite overpassed the 70° north geographical latitude, and thus, the data are not shown. In the Supplementary Materials, more details and other available Swarm magnetic data are represented ( Figure S2).

Figure 6.
Residual of the Y component of the Swarm Alpha magnetic field after the occurrence of the Mw 6.9 2015-11-11 earthquake. The vertical black line in the lower plot represents the point where the satellite crosses the earthquake latitude. The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. We selected this anomaly because it is the first on the track after the earthquake. Moreover, we noticed that its shape is different from the other anomalies, as the amplitude is almost constant (peak-to-peak about 1nT) during its duration (little Figure 6. Residual of the Y component of the Swarm Alpha magnetic field after the occurrence of the Mw 6.9 2015-11-11 earthquake. The vertical black line in the lower plot represents the point where the satellite crosses the earthquake latitude. The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. We selected this anomaly because it is the first on the track after the earthquake. Moreover, we noticed that its shape is different from the other anomalies, as the amplitude is almost constant (peak-to-peak about 1nT) during its duration (little more than 1 min), which is different from the other anomalies that increase, reach a maximum and then decrease. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S3).
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 20 more than 1 min), which is different from the other anomalies that increase, reach a maximum and then decrease. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S3). In the selection, we take into account that the disturbance between 0 and 1 min appears too early to be an acoustic gravity wave with respect to the earthquake origin time, while the selected one is more reliable. The disturbance between 1 and 6 min seems to be characteristic of this orbit. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S4).  Additionally, a zoom of the track around the selected point is provided. In the selection, we take into account that the disturbance between 0 and 1 min appears too early to be an acoustic gravity wave with respect to the earthquake origin time, while the selected one is more reliable. The disturbance between 1 and 6 min seems to be characteristic of this orbit. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S4).
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 20 more than 1 min), which is different from the other anomalies that increase, reach a maximum and then decrease. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S3). In the selection, we take into account that the disturbance between 0 and 1 min appears too early to be an acoustic gravity wave with respect to the earthquake origin time, while the selected one is more reliable. The disturbance between 1 and 6 min seems to be characteristic of this orbit. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S4).  Figure 6 but with lower amplitude. The anomalies after 14 min are probably due to the south auroral magnetic activity. About 18 min after the earthquake, the satellite went below -70° geographical latitude, thus, the data are not shown. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S6).  In this case, we excluded the anomalies closer to the time of the earthquake occurrence (as in previous cases) and we selected the first anomaly after earthquake. Other anomalies are present in the track, while the selected one seems to be symmetric with respect to the epicenter. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S8). shown. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S6).  In this case, we excluded the anomalies closer to the time of the earthquake occurrence (as in previous cases) and we selected the first anomaly after earthquake. Other anomalies are present in the track, while the selected one seems to be symmetric with respect to the epicenter. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S8). In this case, we excluded the anomalies closer to the time of the earthquake occurrence (as in previous cases) and we selected the first anomaly after earthquake. Other anomalies are present in the track, while the selected one seems to be symmetric with respect to the epicenter. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S8).
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 20 Figure 11. Residual of the Y component of Swarm Bravo magnetic field after the occurrence of the Mw 6.4 2019-04-05 earthquake. The satellite passed above the epicenter some minutes before the earthquake occurrence; thus, no vertical black line is represented here (see previous figures where this line was plotted for more details). The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. We selected this point, which is the start of a long "train" of anomalies in the track that persist for more than 9 min; other anomalies start after about 15 min, and thus cannot be generated by acoustic gravity waves. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S9).  . The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. We selected this point, which is the start of a long "train" of anomalies in the track that persist for more than 9 min; other anomalies start after about 15 min, and thus cannot be generated by acoustic gravity waves. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S9).
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 20 Figure 11. Residual of the Y component of Swarm Bravo magnetic field after the occurrence of the Mw 6.4 2019-04-05 earthquake. The satellite passed above the epicenter some minutes before the earthquake occurrence; thus, no vertical black line is represented here (see previous figures where this line was plotted for more details). The vertical green line is the candidate point where the earthquake-induced phenomena are supposed to start. We selected this point, which is the start of a long "train" of anomalies in the track that persist for more than 9 min; other anomalies start after about 15 min, and thus cannot be generated by acoustic gravity waves. In the Supplementary Materials, more details and the other available Swarm magnetic data are represented ( Figure S9).

Discussion and Interpretation
The 12 selected tracks associated to 10 seismic events have been analyzed to search for a possible linear regression between the time interval of the arrival of a supposed co-seismic gravity wave in the ionosphere and the distance between the epicentre and the Swarm satellite, as represented in Figure 14. The intercept has been fixed to the origin of the axis as the time to cover a null distance needs to be zero. The obtained slope is (0.00163 ± 0.00014) min/km that corresponds to a propagation velocity of (10.04 ± 0.85) km/s. The R 2 = 0.92 denotes a good agreement with this linear regression, confirming that the two quantities plotted in the graph seem to be linearly correlated. Although this correlation is pretty good, the propagation speed for this phenomenon seems to be higher than we might expect for this kind of waves. We propose that actually this could be caused by two different involved phenomena: one from the Earth's surface that reaches the maximum electron density altitude (hmF2 or altitude of ionospheric layer F2), and another physical phenomenon that propagates from this altitude to the satellite position. The acoustic gravity waves induced by seismic activities are only a type of a great variety of phenomena that could affect the ionosphere. For this reason, several disturbances affect these figures. Therefore, to measure the time of the arrival of the possible disturbance induced by the seismic event, a particular change in the waveform, concentrating in the part of the track that started after some minutes of the seismic event, was selected. We excluded the first minute of each track, because an acoustic gravity wave requires about 3.5 min to reach Swarm altitude (~450km or more). The latitude of the satellite was also taken into account to understand if the disturbance could be due to the standard geomagnetic activity of the Earth's polar regions. We admit that this interpretation is not totally objective and unique, although sufficiently reasonable. For a quick comparison and validation, we also present in the Appendix A an automatic approach that is totally objective and practically confirms our results.
For the earthquake that occurred on 2014-02-18, only the tracks of Bravo and Charlie (Figures 3  and 4) are provided as Swarm Alpha was not in nominal conditions during the passage over the epicentral area. Most of the data present a discrepancy between ASM and VFM measurements (Flag_F = 16), and the attitude information was interpolated with the nearest points for lack of data from two or all three star-cameras. These conditions were probably due to the early stage of the mission that we analyzed.
An earthquake happened offshore Chile on 2015-11-11 at a depth of 10 km between the Nazca and South-American plates. In this case, a perturbation was visible just north of the epicentre and we cannot exclude that this could be an electromagnetic phenomenon coming directly from the fault. There was a certain symmetry of the signal with respect to the epicentre latitude. The ionospheric disturbance was possibly associated with the high magnitude earthquake (Mw 6.9) considered in this case ( Figure 6).

Discussion and Interpretation
The 12 selected tracks associated to 10 seismic events have been analyzed to search for a possible linear regression between the time interval of the arrival of a supposed co-seismic gravity wave in the ionosphere and the distance between the epicentre and the Swarm satellite, as represented in Figure 14. The intercept has been fixed to the origin of the axis as the time to cover a null distance needs to be zero. The obtained slope is (0.00163 ± 0.00014) min/km that corresponds to a propagation velocity of (10.04 ± 0.85) km/s. The R 2 = 0.92 denotes a good agreement with this linear regression, confirming that the two quantities plotted in the graph seem to be linearly correlated. Although this correlation is pretty good, the propagation speed for this phenomenon seems to be higher than we might expect for this kind of waves. We propose that actually this could be caused by two different involved phenomena: one from the Earth's surface that reaches the maximum electron density altitude (hmF2 or altitude of ionospheric layer F2), and another physical phenomenon that propagates from this altitude to the satellite position. A sketch of the proposed propagation mechanism is shown in Figure 15. A simple model was implemented with a vertical propagation with velocity v1 until hmF2 and a second velocity v2 from the projection of the epicentre at hmF2 until the satellite position. The altitude of the F2 layer was determined by the IRI-2016 model and listed in Table 1. We supposed that the velocity v1 and v2 were equal for all the events, and thus, we determined which combination of v1 and v2 can reduce the difference with the measured times. "v1" and "v2" were selected in the range between 0.10 km/s and 50.00 km/s with a step of 0.01 km/s. We also take into account that the different air temperature T affects the mechanical velocity of propagation by a factor of √ [32]. Thus, the surface (2 m height) temperature has been retrieved from the climatological archive ECMWF ERA-Interim for the 10 earthquakes, and the velocity v1 was multiplied by /298 . The division by 298K assures to obtain the result referred to the standard temperature of 25°C. The best combination that fits our data (i.e., the one with the minimum residuals) was v1 = 2.20 km/s (at temperature of 25°C) and v2 = 16.07 km/s. From the dispersion of the measurements versus the model, we estimated an uncertainty on the velocities of about Δv1 = 0.3 km/s and Δv2 = 3 km/s. We note that the velocity v1 is in agreement with the known velocities in the literature for acoustic gravity waves (e.g., [9]). A comparison between the calculated and the measured times needed to propagate from epicentre to the satellite is represented in Figure 16. The red line represents the ideal case, i.e., when the model was identical to the measurements. The sum of the residual between the time estimated from the model and the measured one from satellite data was reported. The calculated values for the model with two velocities used the best-obtained velocity estimations. We note that in the second model (b), two cases were correctly described by the model as they are on the ideal trend, while the others were under-or overestimated. As also underlined by the lower value of the sum of the residuals, the model with two velocities estimates better the time of propagation of the disturbance from the epicenter to the satellite, with a lower dispersion of data in the plot and closer to the ideal case (red line). This statistical evidence deposes in favour of our assumption of a gravity wave that mechanically propagates until F2 layer making oscillating the high number of free electrons in this A sketch of the proposed propagation mechanism is shown in Figure 15. A simple model was implemented with a vertical propagation with velocity v 1 until hmF2 and a second velocity v 2 from the projection of the epicentre at hmF2 until the satellite position. The altitude of the F2 layer was determined by the IRI-2016 model and listed in Table 1. We supposed that the velocity v 1 and v 2 were equal for all the events, and thus, we determined which combination of v 1 and v 2 can reduce the difference with the measured times. "v 1 " and "v 2 " were selected in the range between 0.10 km/s and 50.00 km/s with a step of 0.01 km/s. We also take into account that the different air temperature T affects the mechanical velocity of propagation by a factor of √ T [32]. Thus, the surface (2 m height) temperature has been retrieved from the climatological archive ECMWF ERA-Interim for the 10 earthquakes, and the velocity v 1 was multiplied by √ T/298K . The division by 298K assures to obtain the result referred to the standard temperature of 25 • C. The best combination that fits our data (i.e., the one with the minimum residuals) was v 1 = 2.20 km/s (at temperature of 25 • C) and v 2 = 16.07 km/s. From the dispersion of the measurements versus the model, we estimated an uncertainty on the velocities of about ∆v 1 = 0.3 km/s and ∆v 2 = 3 km/s. We note that the velocity v 1 is in agreement with the known velocities in the literature for acoustic gravity waves (e.g., [9]).

Figure 15.
Possible mechanism of propagation of the disturbance from the epicentre to the Swarm satellite. From the lithosphere to the F2 ionospheric layer, an acoustic gravity wave is probably the mechanism, and from F2 to the Swarm position, some electromagnetic scattering could be involved (see the text for further discussion).   Figure 15. Possible mechanism of propagation of the disturbance from the epicentre to the Swarm satellite. From the lithosphere to the F2 ionospheric layer, an acoustic gravity wave is probably the mechanism, and from F2 to the Swarm position, some electromagnetic scattering could be involved (see the text for further discussion).
A comparison between the calculated and the measured times needed to propagate from epicentre to the satellite is represented in Figure 16. The red line represents the ideal case, i.e., when the model was identical to the measurements. The sum of the residual between the time estimated from the model and the measured one from satellite data was reported. The calculated values for the model with two velocities used the best-obtained velocity estimations. We note that in the second model (b), two cases were correctly described by the model as they are on the ideal trend, while the others were under-or overestimated. As also underlined by the lower value of the sum of the residuals, the model with two velocities estimates better the time of propagation of the disturbance from the epicenter to the satellite, with a lower dispersion of data in the plot and closer to the ideal case (red line). This statistical evidence deposes in favour of our assumption of a gravity wave that mechanically propagates until F2 layer making oscillating the high number of free electrons in this region and starting a scattering that propagates at least to the Swarm satellite altitudes (about 440-500 km above sea level). Considering also some literature [33,34], we tested the possibility that the second propagation follows the magnetic field lines strictly. We found that the results were worse than the simple straight propagation, thus, even if the perturbation could affect the magnetic field, the disturbance probably did not strictly follow the magnetic field lines. Figure 15. Possible mechanism of propagation of the disturbance from the epicentre to the Swarm satellite. From the lithosphere to the F2 ionospheric layer, an acoustic gravity wave is probably the mechanism, and from F2 to the Swarm position, some electromagnetic scattering could be involved (see the text for further discussion). To finally compare the two approaches describing the data (one-velocity or two-velocity model), we used the Akaike Information Criterion (AIC) developed by Akaike, 1973 [35]. In case of few parameters a "corrected" index AICc can be used. This criterion was used to choose the best model To finally compare the two approaches describing the data (one-velocity or two-velocity model), we used the Akaike Information Criterion (AIC) developed by Akaike, 1973 [35]. In case of few parameters a "corrected" index AICc can be used. This criterion was used to choose the best model (the one that presents the minimum AIC, or AICc, later described) and select which variables were really fundamental to the system. The AIC is defined as: where K is the number of estimated parameters, i.e., the degree of freedom and L is the maximum of the likelihood function. For this paper, we used the following equation to estimate the log(L) : where n is the number of samples and RSS is the Residual Sum of Squares. If the dataset on which the model is evaluated contains a small number of elements, as in our case, a correction factor must be introduced, estimating a second coefficient, called AICc: To calculate the values of AIC and AICc for our two models, we followed the approach of Snipes and Taylor, 2014 [36], even if they applied this criterion to a vastly different case study. Table 2 reports the parameters and the AIC and AICc coefficients obtained for the two models. The number of parameters has been increased by one to consider an error coefficient in the model. Despite the height of the F2 layer (hmF2) and the temperature are two parameters used in the second model, they do not add any degree of freedom to the same model, since they are fixed. The AICc (and even AIC) suggests that the model with two velocities better describes the magnetic disturbance identified in the 12-track dataset. From a physical point of view, this is promising, since the model with one velocity gives a speed that is not compatible with known phenomena of the interaction of lithosphere, atmosphere and ionosphere; otherwise, the first velocity of the second model is compatible with an acoustic gravity wave.

Conclusions
After a systematic analysis of the first 5.3 years of Swarm magnetic data around the Mw5.5+ earthquakes, we found that 12 Swarm satellite tracks recorded disturbances following earthquakes of magnitude Mw5.6 to 6.9. We underline that we limited our search for possible co-seismic effects to a subset of earthquakes selected in function of their depth and magnitude. Even if this subset could exclude some possibility to record some additional effect, we were driven by the necessity to exclude the deeper and lower magnitude events that are unfeasible to produce any seismic effect in the ionosphere. In addition, we do not expect particular difference with other selection criteria.
However, more energetic seismic events were recorded during this period (included three Mw8+ earthquakes: two in Chile and one in Mexico), but we did not find possible co-seismic ionospheric disturbances produced by these large events. This is due to the fact that the Swarm satellites did not fly just over the epicentral area at the time of the earthquake occurrence. This explains why we found geomagnetic ionospheric disturbances only for almost medium magnitude earthquakes.
The linear regression of the time versus the distance between the earthquake occurrence and the record of the disturbance by Swarm satellites seems to confirm that the seismic events possibly induce the disturbances.
To explain these disturbances, a model that combines two different propagation ways has been proposed with promising result. The model considers an acoustic gravity wave as the way to propagate the earthquake disturbance up to the F2 ionospheric layer, i.e., the layer with the maximum ionisation. The impact of the acoustic gravity waves on this layer could induce some electromagnetic scattering that propagates in the upper ionospheric layer, reaching the Swarm satellites. The propagation velocities, according to our data, are (2.2 ± 0.3) km/s with a 25 • C of surface air temperature for the acoustic gravity waves and (16 ± 3) km/s for the electromagnetic scattering propagation. We note that the approximation to the same velocity for all the seismic events is a limitation and, in particular, it does not take into account, e.g., the different local time: future improvements could introduce the importance of the local time together with some other parameters that would take into account the different chemical and physical conditions of the atmosphere and ionosphere where the phenomena propagate. This double-velocity model has better statistical reliability, providing a lower AIC than the one with only one velocity.
In order to validate the mostly subjective criteria chosen in this work to select the co-seismic anomalies, we developed an automatic algorithm that detects the co-seismic perturbations based on objective criteria (see Appendix A). The principal results of this paper were also confirmed by this automatic approach to detect the anomalies. Nevertheless, the automatic system was not yet able to achieve the quality with which a researcher can evaluate the signal of the satellite track to discriminate and choose the effects and causes of the detected disturbances. A future improvement of the automatic procedure can include the analysis of the signal frequency for not only searching a disturbance along the satellite track but also some particular frequencies without usually presenting in other parts of the track.
A larger constellation of Low Earth Orbit satellites equipped with magnetometers could improve the coverage, i.e., the capability to record and deeper understand this type of signal. The recent launch of the Chinese Seismo-electromagnetic satellite in 2018 could contribute to this endeavour greatly.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/7/1166/s1. Figures S1-S11: The output tracks where some possible co-seismic effects on magnetic field have been detected and one example where we do not found any evidence of possible earthquake-source disturbance.  The distance from the epicentre and satellite versus the time for the 31 tracks is represented in Figure A1. The time to reach the epicentre from the hypocenter has been estimated by a velocity of 4 km/s and subtracted from the data. The linear regression provides a slope of (2.60 ± 0.41) 10 −3 min/km that corresponds to a mean propagation velocity of (6.4 ± 1.0) km/s. The obtained velocity is quite similar to the one calculated in the main text, but the dispersion of the data is higher (in fact, the correlation coefficient is much lower). track, i.e., rms/RMS > kt. This threshold corresponds to about two standard deviations commonly used in similar approaches. Thresholds of 1.5, 1.75, 2.25 and 2.5 have been tested too, but they provide worst results (i.e., lower correlation coefficient of the results) with respect to kt = 2.0.
The automatic procedure found anomalies in 31 tracks out of 37 investigated tracks, i.e., 83.8%. This number is much higher than the selected anomalies in the main text (12). Despite this, we need to note that several anomalies are not reliable; for example, they occur too early in relation to the distance from the earthquake epicenter or too late and close with respect to the earthquake epicenter. The latter case could be a continuation of the co-seismic effect in ionosphere after its arrival.
The distance from the epicentre and satellite versus the time for the 31 tracks is represented in Figure 17. The time to reach the epicentre from the hypocenter has been estimated by a velocity of 4 km/s and subtracted from the data. The linear regression provides a slope of (2.60 ± 0.41) 10 −3 min/km that corresponds to a mean propagation velocity of (6.4 ± 1.0) km/s. The obtained velocity is quite similar to the one calculated in the main text, but the dispersion of the data is higher (in fact, the correlation coefficient is much lower). In analogy with the analyses of the main text, the model with two velocities has been tested with the same approach, using all the 31 anomalies and taking into account the air temperature and the height of hmF2 in addition to the simple model with one velocity. Figure 18 shows the results obtained by the models with one and two velocities for the automatically detected anomalies. Additionally, in this case, the model with two velocities presents lower residuals with respect to the first one. We note that the residuals are about 10 times more than the corresponding calculus in the Figure A1. Time after the earthquake occurrence versus the distance between the epicentre and the satellite for the 31 anomalies automatically detected. The time of propagation of the seismic wave from the hypocenter to the epicentre has been subtracted.
In analogy with the analyses of the main text, the model with two velocities has been tested with the same approach, using all the 31 anomalies and taking into account the air temperature and the height of hmF2 in addition to the simple model with one velocity. Figure A2 shows the results obtained by the models with one and two velocities for the automatically detected anomalies. Additionally, in this case, the model with two velocities presents lower residuals with respect to the first one. We note that the residuals are about 10 times more than the corresponding calculus in the main text ( Figure 16). This is partly due to the higher number of cases (31 instead of 12), the longer analyzed time (until about 30 m instead of 20 m) and, finally, to the lower quality of the co-seismic anomaly detection of the automatic approach. The best obtained velocities are v 1 = 0.7 km/s (at temperature of 25 • C) and v 2 = 14.8 km/s. The first velocity is very low with respect to the one obtained with the supervised track selection (2.2 km/s) and also too low for a gravity wave. The second velocity is closer to the result in the main text (16.1 km/s). Furthermore, this automatic approach confirms that v 2 is greater than v 1 and, finally, confirms the order of magnitude of the results.
analyzed time (until about 30 m instead of 20 m) and, finally, to the lower quality of the co-seismic anomaly detection of the automatic approach. The best obtained velocities are v1 = 0.7 km/s (at temperature of 25 °C) and v2 = 14.8 km/s. The first velocity is very low with respect to the one obtained with the supervised track selection (2.2 km/s) and also too low for a gravity wave. The second velocity is closer to the result in the main text (16.1 km/s). Furthermore, this automatic approach confirms that v2 is greater than v1 and, finally, confirms the order of magnitude of the results.