Ionospheric response to the 2020 Samos earthquake and tsunami

On 30 October 2020 at 11:51 UT, a magnitude 7.0 earthquake occurred in the Dodecanese sea (37.84°N, 26.81°E, 10 km depth) and generated a tsunami with an observed run‑up of more than 1 m on the Turkish coasts. Both the earthquake and the tsunami produced acoustic and gravity waves that propagated upward, triggering co‑seismic and co‑tsunamic ionospheric disturbances. This paper presents a multi‑instrumental study of the ionospheric impact of the earthquake and related tsunami based on ionosonde data, ground‑based Global Navigation Satellite Sys‑ tems (GNSS) data and data from DORIS beacons received by Jason3 in the Mediterranean region. Our study focuses on the Total Electron Content to describe the propagation of co‑seismic and co‑tsunami ionospheric disturbances (CSID, CTID), possibly related to gravity waves triggered by the earthquake and tsunami. We use simultaneous vertical ionosonde soundings to study the interactions between the upper and lower atmosphere, highlighting the detection of acoustic waves generated by the seismic Rayleigh waves reaching the ionosonde locations and propagating verti‑ cally up to the ionosphere. The results of this study provide a detailed picture of the Lithosphere‑Atmosphere–Iono‑ sphere coupling in the scarcely investigated Mediterranean region and for a relatively weak earthquake.

Ionospheric perturbations generated by earthquakes are called co-seismic ionospheric disturbances (CSID).They can be detected in the vicinity of an earthquake's epicenter, but also up to 6000-9000 km away from it, when exceptional events occur like the 2011 Tohoku earthquake (Chum et al. 2012).The near-field (i.e., within ~ 1000 km from the epicenter) CSID usually exhibit an interference of ionospheric disturbances generated by acoustic waves due to co-seismic crustal vertical displacements and Rayleigh surface waves (e.g., Davies and Baker 1965;Row 1966;Calais and Minster 1995;Afraimovich et al. 2001;2010;Artru et al. 2004;Heki and Ping 2005;Liu et al. 2006aLiu et al. , 2011;;Rolland et al. 2011aRolland et al. , 2013;;Perevalova et al. 2014;Astafyeva and Shults 2019).As the acoustic waves propagate at the speed of sound, i.e., ~ 330 m/s at the surface and 800-1000 m/s at the ionospheric altitude of 250-300 km, the CSID are typically registered ~ 7-9 min after an earthquake and are often N-shaped, partly repeating the waveform of the acoustic waves with compression and rarefaction phases.In some rare cases, free gravity waves associated with ground displacements can also be observed after major earthquakes (e.g., Rolland et al. 2011b).
Additionally, tsunamis generate internal gravity waves that, unlike acoustic waves, propagate obliquely and have lower vertical velocities because of the gravity action.Therefore, these disturbances reach the ionosphere altitudes of 250-400 km ~ 45-60 min after they are generated on the water surface.The co-tsunamic ionospheric disturbances (CTID) mainly match the period, velocity, and propagation direction of the tsunamis causing them, and are usually characterized as quasi-periodic structures with typical periods between 10 and 30 min (e.g., Artru et al. 2005;Liu et al. 2006b;Rolland et al. 2010;Makela 2011;Galvan et al. 2012;Occhipinti et al. 2013;Komjathy et al. 2016).CTID can be detected in the vicinity of earthquakes' epicenters, but also in the open sea/ocean and/or upon their arrival on shores.
CSID and CTID can be detected by using groundbased and space-borne data.Typically, they appear as wave-like perturbations of the plasma visible in the Total Electron Content (TEC) and in electron density variations, measured by High Frequency (HF), Very High Frequency (VHF) and L-band signals.
The earthquake triggered a tsunami that caused significant damages on the north coast of Samos and in Sigacik Bay where one fatality was reported (Yalciner et al. 2020), but also affected the wider coastal area of Izmir Province of Turkey (Fig. 1).This was the largest tsunami to hit the Eastern Mediterranean Sea since 1956, despite the high tsunami hazard of the region, with around three hundred historical tsunamis recorded since ancient times (Ambraseys 1962).Early field investigations allowed the collection of numerous testimonies, along with runup and run-in field evidence (Triantafyllou et al. 2020;Aksoy 2021;Dogan et al. 2021;Kalligeris et al. 2021).The waves were also recorded by various video sources, providing strong constraints, especially temporal, for the tsunami modelling (e.g., https:// www.in.gr/ 2020/ 10/ 30/ world/ seism os-thala ssa-vgike-sti-steria-sti-smyrni-ormit ika-nera-mesa-stin-poli/). Regrettably, there was no tide gauge data available in the near field of the earthquake.The run-up height was over one meter in several locations, and oscillations with periods of around ten minutes were observed, resulting in up to 61 run-ups with significant land flooding (https:// www.ngdc.noaa.gov/ hazard/ tsu_ db.shtml).
In this work, we provide a full model of the 2020 Samos tsunami, and subsequently focus on the ionospheric effects observed in regions up to 1300 km apart from the epicenter thanks to multiple observing platforms (GNSS receivers, ionosondes and Digisonde stations, beacon receiver on board of Swarm and Jason3 satellites).Our study aims to fill the knowledge gap of the scarcelystudied ionospheric response to smaller earthquakes (M w < 8.0), with the peculiarity of analyzing a case study that occurred in a narrow region, the Mediterranean area.

Tsunami modelling
We modelled the tsunami using the TAITOKO code, an in-house numerical code developed in the late 1990s for tsunami early warning and widely used for tsunami hazard prevention (Heinrich et al. 1998;Hébert et al. 2001).To simulate the ground deformation caused by the earthquake, we assumed it to be the static surface displacement of a homogeneous elastic half-space generated by an instantaneous rectangular dislocation.The code uses Okada's equations (Okada 1985) to simulate this deformation according to the seismic parameters of the estimated source.It then transforms this deformation through the water column, as described in Glimsdal et al. (2013), and propagates it using either Boussinesq or shallow-water equations (Jamelot et al. 2019;Heinrich et al. 2021).
For this study, we used the fault model calculated by Ganas et al. (2021).The fault plane is well-constrained by inverting geodetic data, i.e.Synthetic Aperture Radar Interferometry (InSAR) and GNSS data, constrained by seismological, and geological data, including hypocenter localization, focal mechanism and aftershocks, faults maps, and field observations.The modelled fault plane strikes N276°E, i.e., roughly east-west, dipping north on a 37°-dipping plane and produces a rupture of pure normal rake (− 90°) with a mean slip of 1.7 m on a rectangular patch measuring 37 km in length and 15 km in width.The center of the patch is located at 26.678° E, 37.829° N, and has a centroid depth of 5.1 km (Fig. 1).
To better model the local effects of bathymetry, the TAITOKO code allows the use of multiple nested grids, with coarse grids over deep-water regions and fine grids over coastal regions.For the large-scale propagation, we used the 2020 European Marine Observation and Data Network (EMODNet, https:// emodn et.eu) bathymetry, with a resolution of 120 m.Given that the water depth rarely exceeds 1 km in the Aegean Sea, which is much less than the typical tsunami wavelength of ~ 10 km for this event, despite the limited fault rupture and small basin, we used the shallow water approximation for the numerical modelling.We produced local high-resolution (5 × 5 m) bathymetric grids by interpolating depths retrieved from nautical soundings in the area in four harbors where reliable observations of the tsunami waves were available (Triantafyllou et al. 2020;Aksoy 2021;Dogan et al. 2021;Kalligeris et al. 2021): Vathy Bay, where strong amplification was observed, Karlovasi Harbor, where early observations were made, Kusadasi Harbor, Tiparcik Harbor, and Siğacik Harbor, where video records show strong constructive interference.In these areas, we assume the accuracy to be metric, as required for nautical charts.We also incorporated an onshore 5 × 5 m Digital Elevation Model to retrieve the locations of streets and city blocks, enabling us to precisely match run-up simulations with observations in urban areas.To validate the accuracy of our models, we compared the simulated run-in and runup with field evidence, witness testimonies, pictures, videos, and digital recordings of existing far-field tide gauges.

Ionospheric response
To investigate the ionospheric effects observed in the region of interest, we exploited data acquired by different observing platforms, both ground-based (GNSS receivers, ionosondes and Digisonde stations) (Fig. 2) and space-based (DORIS/ Jason3).
During the earthquake, eight GNSS (GPS and GLO-NASS) satellites were visible by ground-based GNSS receivers covering the region of interest (Fig. 2).To search for the ionospheric TEC response to the 2020 Samos earthquake and the following tsunami, we applied two methods: a 6-28 min window as band pass running mean filter (Astafyeva et al. 2009) and the VARION (Variometric Approach for Real-Time Ionosphere Observation), algorithm (Savastano et al. 2017).VARION is an open source, entirely Python-based software (https:// github.com/ giorg iosav astano/ VARION), described by Savastano et al. (2017) andRavanelli et al. (2021).It is based on the computation of the integral over a certain interval of the time differences of geometry-free combinations of carrier-phase measurements from a stand-alone GPS station, that reads: where t i represents the initial time of the considered period, L S 4R is the geometry free combination of the carrier phase measurements calculated considering the receiver R and the satellite S and f 1 , f 2 are the GPS L1 and L2 signal frequencies, respectively.
VARION makes use of the standard orbit and clock products, and it is based on a thin layer approximation of the ionosphere.As the time difference of the carrier phase measurements, the effect of the inter-frequency biases on TEC evaluation can be ignored since they can be considered constant along every single arc if no cycles slips occur.The long-term trend from TEC VAR time- series is removed, computing the residuals with respect to a 10th order polynomial fit.
Time series of TEC VAR obtained from a set of GNSS receivers located close to the epicenter are plotted with respect to the minutes after the earthquake main shock to look for the ionospheric response to earthquake and tsunami.
From the network of GNSS receivers in the considered region (described in Fig. 2), we also applied the TOMION technique (TOmographic Model of the (1) IONosphere; Hernández-Pajares et al. 2020a) to estimate the electron density from GNSS data.TOMION has been run in a high resolution mode from GNSS data, assuming the electron density (Ne) distributed in voxel of 1° × 680 km (2 layers) × 30 s in Right Ascension (α) Latitude (φ) × Height (h) × Time (t).And, in order to better adapt it to the expected local variations of electron density due to the tsunami, the random walk process noise has been increased one order of magnitude regarding the value used for computing the Global Ionosphere Maps (GIMs) for International GNSS Service, where the Universitat Politecnica de Catalunya Quarter-of-an-hour time resolution Rapid GIM (UPC GIM UQRG) is the slightly best behaving one (Fig. 3).
This technique allows us to estimate the Ne time derivative, latitudinal and altitudinal gradients at the two layers during the earthquake day and the previous ones.The temporal and spatial gradients are given as maps and videos to show up the ionospheric response.
Beyond the ground-based dual-frequency GNSS measurements, input for the tomographic estimation with TOMION, we have also analyzed detrended DORIS (Doppler Orbitography by Radiopositioning Integrated on Satellite) dual-frequency measurements.The corresponding doppler integrated measurements were gathered from the Jason-3 LEO satellite from the groundbased transmitters, considering an effective height of 250 km and a detrending time of 20 s (Yang et al. 2022).
INGV operates an AIS-INGV ionosonde in Gibilmanna, Sicily (37.9°N, 14.0° E) (www.eswua.ingv.it).Our analysis is based on the ionograms inspection to look for unusual deformation of the traces, possibly due to the earthquake.In particular, we looked for Multiple-Cusp Signatures (MCS) caused by traveling co-seismic ionospheric disturbances (Maruyama et al. 2012(Maruyama et al. , 2016) ) generated by the propagating surface Raylegh waves.

Tsunami modelling
The calculated runup height at Seferihisar, Turkey, and Vathy, Greece, was found to be around 3 m, which is significantly higher than the surface motion generated by the earthquake.This emphasizes the significance of local effects on the runup and the importance of simulations in anticipating potential damage and secondary effects such as ionospheric perturbations.Bays have the ability to amplify waves, while rugged coastlines like Tiparcιk can cause complex and sometimes constructive interferences, leading to significant amplification of the waves which then propagate into the Aegean Sea.
The model accurately reproduces the series of waves recorded in the inundation videos.The oscillations in the model have a period of around 15 min at Vathy and remain stable over time.This consistency is observed across multiple simulations with varying sources, indicating that the coastline plays a significant role in determining the wave parameters.The surface propagation of the tsunami wave aligns with the evolution of the electron density vertical gradient component, shifted by a vertical propagation delay of 41 min and 30 s empirically deduced from the starting time of the perturbed vertical gradient component (see corresponding animation in the supplementary material).Given a boundary altitude of 790 km between the two ionospheric layers, this delay suggests a vertical speed slightly exceeding 300 m/s.This speed aligns with the characteristics of acoustic waves at ground, while the acoustic speed in the ionosphere is significantly higher.However, the perturbed vertical gradient in electron density lasting over 20 min aligns more closely with the characteristics of gravity waves that follow.

Ionospheric response from GNSS and satellite in-situ observations
Even if we have not found a clear co-seismic signature in the TEC analysis, commonly reported to be N-wavelike disturbance occurring ~ 8-9 in after the earthquake, we notice a perturbation in the satellite G07 data.This might be related to the earthquake or the tsunami, as it seems to correlate with the epicenter on a Travel-time diagram (TTD, Fig. 4a).This TTD is done for an ionospheric height, Hion, equal to 250 km, and the US Geological Survey epicenter is taken as the source.From this diagram, we estimate the apparent horizontal velocity to be around 120 m/s.CTIDs are known to have the same or very similar speed as the tsunami that generates them.A velocity of 120 m/s would be consistent with the expected propagation velocity of a tsunami in the shallow waters of the Aegean Sea, given by √(gh), where h represents the water depth and g the gravitational acceleration, giving an expected value of around 100 m/s for a 1-km-deep basin.However, the CTID are driven by gravity waves and are expected about 40-60 min after the tsunami generation, while, on Fig. 4a one can see the perturbation occurs ~ 10-12 min after the earthquake.The TEC VAR derived from the same satellite confirms the perturbation detection at a similar time with a comparable estimated velocity (Fig. 4c).It should be noted that changing the altitude of detection changes (slightly) the location of sub-ionospheric points, and, therefore, the TTD figure.We chose H ion = 250 km because we obtain the best alignment in the perturbations detected by different stations.
Another satellite (G30) had much lower elevation angles, it only started to be visible for the GNSS receivers at the time of the earthquake.For these observations, we take H ion = 180 km.However, the signal is less clear in the data of this satellite, and the correlation with the epicenter or a tsunami source is more difficult to prove.Anyhow, the analysis shows a TEC perturbation occurring after about 30 min, with a propagation velocity inferior to 100 m/s (Fig. 4b).Also in this case the TEC VAR detects the signature of ionospheric response moving with a similar velocity (Fig. 4d).
We compute the gradients of Ne by TOMION with respect to time and to longitude, latitude and height for the earthquake day and for the previous day (for comparison), without any background model and from the dual-frequency GNSS carrier phase measurements provided by permanent receivers in the earthquake region (see Fig. 2).Two ionospheric shells have been considered, to guarantee sensitivity to the vertical drift under the available ground-based observation geometry, taking into account both ionosphere and plasmasphere (between 110 and 790 km, and between 790 and 1470 km).
The ionospheric response to the earthquake is evident in the vertical component of the Ne gradient estimated from the voxels-based tomographic solution among the two layers, with central heights at 450 km and 1130 km, respectively (Fig. 5), while no clear signature has been detected in the latitudinal, longitudinal and temporal gradients (not shown).In this way, we can coarsely assign the height to such a vertical gradient component in between both layers, at 790 km height.The perturbation in Fig. 5 is observed starting from about 45 min after the main shock of the earthquake and the detection at 790 km suggests it might be the signature of the gravity waves generated by the tsunami waves, causing detectable effects up to the thermosphere.The works by Vadas 2007 andVadas et al. (2015) report electron density perturbations due to dissipation of the gravity waves at thermosphere altitudes which, by releasing momentum and energy into the ionosphere, causes an acceleration of the plasma in the propagation direction.However, since the ion-neutral collisional frequency at 790 km height is too small for their direct interaction, we consider also the possibility that the dissipation of gravity wave energy and momentum into the ionosphere occurs in the lower altitudes and some other processes, such as electrodynamical forces and momentum transfer along the magnetic field lines within the ionosphere, propagate the deposited energy and momentum to the ionosphere at higher altitudes.
Our result, obtained with the GNSS-based ionospheric tomographic provided by TOMION previously assessed in similar conditions vs different techniques and scenarios (Kotov et al. 2018(Kotov et al. , 2019;;Jerez et al. 2020;Hernández-Pajares et al. 2020a, 2020b, 2022;Wielgosz et al. 2021;Porayko et al. 2023), would expand the observing altitude of this kind of perturbation, beyond the so-far assumed heights of 400-500 km.This result, which should be confirmed in future events, would be compatible with the detection of potential tsunami ionospheric signature in Yang et al. (2022) in the Total Electron Content measurements above Swarm LEOs, flying at 460-510 km.Indeed, in such a work the tsunami perturbation was detected at both in-situ (Langmuir probe) electron density and in the line-of-sight path integral of the electron density (i.e.TEC) from the Swarm POD GNSS receivers.The end of the perturbation at 500 km, i.e. at a very short length above LEO, should have not produced any signal at TEC, in front of the typical topside distribution of electron density up to the end of the ionosphere and beyond, the plasmasphere.But this was not the case.This result can be compatible with the TEC enhancement of over 2 TECu magnitude observed above 530 km near the Tonga volcanic eruption site by Liu et al. (2023).
Looking at the sequence of the maps of the vertical component of the Ne gradient detected with TOMION after the main shock (Fig. 6 and full sequence in the supplementary material), we identify the Ne gradients maxima (starting about 45 min after the earthquake) appearing on the south from the epicenter, as confirmed also by the ΔTEC analysis shown in Fig. 4e,  f.Such a result agrees with what is expected in the northern hemisphere due to the attenuation of the ionospheric disturbances propagating northward caused by the geomagnetic field configuration (Astafyeva 2019).
Also Fig. 7 confirms that the southward propagation of the ionospheric disturbances is consistent with its continuation later on as observed on the DORIS track at 13:19 UT.The onset of the ionospheric response is observed in Athens (260 km from the epicenter) in the ionogram of 12:00 UT, showing MCS at F2 in the virtual height range between 300 and 400 km and in the frequency spanning from 6 to 7 MHz (Fig. 8a).The seismogram acquired by the seismic station ATHU, located 274 km from the epicenter, allows picking the Rayleigh waves arrival at 11:53:03 UT.
The ionosonde in Nicosia (640 km from the epicenter) recorded MCS in the F1 and F2 traces at 12:02 in the virtual height range between 150 and 250 km and in the frequency spanning from 3 to 4 MHz (Fig. 8b).The seismogram acquired by the station LFK, located at 659 km from the epicenter, allows picking the Rayleigh waves arrival at 11:55:50 UT.
The ionograms recorded in Gibilmanna (1100 km from the epicenter) show MCS on F1 trace at 12:15 UT at virtual altitudes between about 100 km and 150 km and frequency ranging between 3 MHZ and 4 MHz (Fig. 8d).The seismogram acquired by the station GIB, located at 1127 km from the epicenter, allows picking the Rayleigh waves arrival at 11:58:31 UT.
The ionograms recorded by the ionosonde operating in Rome (1300 km from the epicenter) report the The times of MCS detection identified by all the considered ionosondes are compatible with the arrival of acoustic waves in the ionosphere launched by the Rayleigh surface waves.Table 1 summarizes the ionospheric response observed by the ionosonde and the digisondes considered in our analysis.
Similarly to what was done by Maruyama et al. (2012), we estimated the propagation velocity of the Rayleigh waves from the appearance of MCS on the ionograms, considering the sounding cadence as error bars associated with the arrival time of the waves in the ionosphere.

Summary and conclusions
It is now well known that nearly all earthquakes with M w > 6.8 can generate disturbances in the ionosphere (Perevalova et al. 2014).For instance, the 2020 Samos earthquake, that occurred in the Dodecanese sea, and the associated tsunami, caused perturbations in the ionosphere observed by a network of ground-based GNSS receivers and ionosondes located in the Mediterranean area and by Jason 3 flying over the same region.
Our results show: The GNSS TEC variations southward of the epicenter, appearing from 10 to 30 min after the main shock and traveling with a speed ranging from 100 to 120 m/s.The time delay from the earthquake is too long to be attributed to CSID caused by acoustic waves (8-9 min), but too short to be related to gravity waves (~ 45 min).The velocity is compatible with a tsunami propagating in the sea with shallow bathymetry, as in the Mediterranean case, so the TEC variation could be a CTID in the form of ionospheric response to acoustic-gravity waves passage.
The tomographic reconstruction of the electron density from the GNSS TEC identifies an enhancement of one order of magnitude (with respect to the previous day) in the vertical component of Ne gradient about 45 min after the earthquake and at an altitude coarsely assigned to 790 km.We interpret the plasma uplift as the signature of the propagation up to the thermosphere of the gravity waves generated by the Tsunami waves.Indeed, according to previous computations the dissipation of the gravity waves at thermosphere altitudes, releasing momentum and energy to the ionosphere, causes an acceleration of the plasma in the propagation direction that could support the CTID signature we observed in the topside ionosphere.Moreover, recent findings from  LEO satellite observation confirm the possibility to detect CTID above 500 km.
The tomographic map sequence of the vertical component of Ne gradient indicates that from about 45 min after the earthquake onwards, the electron density maxima appear south of the epicenter, in approximately the same location as already identified by the ΔTEC analysis.The DORIS/Jason 3 observations also confirm a southward propagation.Such a result is in accordance with the expected attenuation of the CSID in the Northern Hemisphere due to the geomagnetic field configuration (Astafyeva 2019).
From the Multiple Cusp Signatures observed on the ionograms of the considered network, we clearly observed the propagation of the acoustic waves triggered by the Rayleigh surface waves detected by nearby seismic stations.From the times of the MCS appearance with respect to the distance from the epicenter, we derived a velocity of propagation of the acoustic waves in the ionosphere of 2225 ± 829 m/s, in agreement with the Rayleigh waves velocity estimated by the seismograms, of 2940 ± 390 m/s.
In conclusion, our study reports an original description of the ionospheric response to the earthquake, the related tsunami and the Rayleigh wave propagation due to the Samos earthquake on 30 October 2020 in the Mediterranean area.The novel aspects of our multiinstrumental analysis are: (a) the evidence of topside ionospheric disturbances possibly associated with effects up to the thermosphere due to the dissipation of the gravity waves induced by the tsunami; (b) the signature of TEC variation by multi-technique analysis likely related to acoustic-gravity waves induced by the tsunami; (c) the detection of the acoustic waves propagation in the ionosphere through the analysis of MCS in the ionograms recorded by the selected ionosondes.
Our work aims to contribute to the general understanding of the ionospheric sensitivity to earthquakes and tsunamis, with a special focus on the Mediterranean region, also in view of the possible development of tsunami early warning tools.

Fig. 1
Fig.1Modeled vertical permanent deformation (negative values corresponds to subsidence) due to the Samos Earthquake (colored and with contour lines), using fault parameters ofGanas et al. (2021)

Fig. 2
Fig. 2 Location of the instruments considered in the study: ionosonde and digisondes (red stars), GNSS receivers (green dots).The epicentre is identified as a blue diamond

Fig. 5
Fig. 5 Example of vertical component of the Ne gradient estimated from TOMION among the two-layers, with central heights at 450 km and 1130 km during the earthquake day (top) and during the previous day (bottom) in the bin 26°-28°N, 36°-38° E

Fig. 6
Fig. 6 Sequence of maps of the vertical component of the Ne gradient by TOMION after the main shock Figure 9 reports for each considered ionosonde the time of the last quiet ionogram (empty circle) and the time of the first MCS ionogram (filled circle) with respect to the distance of the ionosondes from the epicenter.The three blue lines are obtained by linear fitting of the timing of last quiet ionogram (dashed line on the left), the timing of the first disturbed ionogram (dashed line on the right) and the mid-point between the two timings (continuous line).The resultant slopes indicate a velocity of 2225 ± 829 m/s, in accordance with the Rayleigh waves velocity estimated by the seismograms (red line) equal to 2940 ± 390 m/s.

Fig. 9
Fig. 9 Travel-time diagram of CSID causing MCSs (blue) and of Rayleigh waves speed (red).The times of the last quiet ionogram and the first MCS ionogram are represented by an empty and a filled circle, respectively.The bars refer to the ionograms cadence.The red star indicates the earthquake time.The three blue lines are obtained by linear fitting of the timing of last quiet ionogram (dashed line on the left), the timing of the first disturbed ionogram (dashed line on the right) and the mid-point between the two timings (continuous line)

Table 1
Ionospheric and seismic parameters related to MCS ionograms