Ionospheric Response to the 6 February 2023 Turkey–Syria Earthquake

Abstract

Two strong earthquakes occurred in Turkey on 6 February 2023, at 01:17:34 (nighttime, Mw = 7.8) and at 10:24:50 UT (daytime, Mw = 7.5). The seismo-ionospheric impact is an important part of the near-Earth environment state. This paper provides the first results on the ionospheric effects associated with the aforementioned earthquakes. We used data from global navigation satellite system (GNSS) receivers and ionosondes. We found that both earthquakes generated circle disturbance in the ionosphere, detected by GNSS data. The amplitude of the ionospheric response caused by daytime M7.5 earthquake exceeded by five times that caused by nighttime M7.8 earthquake: 0.5 TECU/min and 0.1 TECU/min, respectively, according to the ROTI data. The velocities of the earthquake-related ionospheric waves were ~2000 m/s, as measured by ROTI, for the M7.5 earthquake. TEC variations with 2–10 min periods showed velocities from 1500 to 900 m/s as disturbances evolved. Ionospheric disturbances occurred around epicenters and propagated to the south by means of 2–10 min TEC variations. ROTI data showed a more symmetric distribution with irregularities observed both to the South and to the North from 10:24:50 UT epicenter. The ionospheric effects were recorded over 750 km from the epicenters. Ionosonde located 420/490 km from the epicenters did not catch ionospheric effects. The results show significant asymmetry in the propagation of coseismic ionospheric disturbances. We observed coseismic ionospheric disturbances associated with Rayleigh mode and acoustic modes, but we did not observe disturbances associated with acoustic gravity mode.

1. Introduction

The 6 February 2023 earthquakes in Turkey and Syria were among of the largest disasters affected the region for many years. According to the United States Geological Survey (USGS, https://earthquake.usgs.gov/ (accessed on 20 April 2023)), the first major event, with magnitude Mw = 7.8 (also known as Pazarcik earthquake), occurred at 01:17:34 UTC in southern Turkey, near the northern border of Syria. It was the strongest earthquake in Turkey since 1939. The epicenter was located at 37.220°N 37.019°E. The focal mechanism corresponded to strike–slip faulting at a depth of about 10 km.

Approximately nine hours later, the M7.8 event was followed by a M7.5 shock (known as the Elbistan earthquake). The M7.5 earthquake occurred at 10:24:50 UTC and was located at 38.016°N 37.206°E (95 km north-northeast from the first one); the hypocenter was at a depth of 15.0 km. The M7.8 and M7.5 earthquakes triggered aftershock sequences. More than 100 aftershocks with magnitude M4.5 or greater were recorded within 24 h of the M7.8 earthquake (https://earthquake.usgs.gov/ (accessed on 20 April 2023)).

Lithosphere–atmosphere–ionosphere coupling is a complex field which includes several channels for energy transfer as well as mechanisms of chemical composition change. Among energy transfer channels, acoustic gravity wave (AGW) channels cause a wavelike structure in the electron density at the ionosphere peak height. The Earth’s surface displacement during earthquakes causes atmospheric perturbations (including AGW) which then propagate upward.

Acoustic gravity waves are known to be a source of ionospheric electron density disturbances. Liu et al. [1] made an estimate of how ground level displacement should be amplified at ionospheric height and found that the amplification factor is between 4 and 5 orders. While the displacement of the Earth’s surface is rather small, such an amplification can cause noticeable oscillations of the neutral gas in the upper atmosphere. When these oscillations are coupled with ionospheric particles, one can see their manifestation by means of ionospheric sensors.

Earthquakes can produce significant ionospheric effects [2,3,4], including effects due to secondary events like tsunami [5]. Among ionospheric effects, we consider variations in total electron content that reflect a propagating disturbance. The effect scale is connected with the earthquake magnitude and focal mechanism [6]; additionally, a threshold may exist when an ionospheric response appears [7]. The typical traveling ionospheric disturbance speed and periods are associated with the particular mechanism.

There are different kinds of ionospheric irregularities associated with earthquakes. Scientists distinguish among them using the velocity at which the irregularities propagate from the epicenter. The fastest is the irregularity associated with Rayleigh waves (~2000–3000 m/s). Rayleigh waves propagate from the epicenter in the solid Earth, and every piece of displaced surface serves as a secondary source of atmospheric fluctuations. The fluctuations propagate upward, causing ionospheric disturbances. The AGW propagating from the epicenter takes a longer time to reach the same region in the ionosphere, since propagation occurs entirely in the atmosphere. Acoustic (~700–1200 m/s) and gravity (~150–300 m/s) modes are distinguished by means of velocity.

Yasyukevich et al. [8] showed that the ionospheric effects associated with Rayleigh waves had north–south asymmetry for the Tohoku earthquake. Kherani et al. [5] simulated the coseismic TEC disturbances from the Tohoku earthquake and showed that AGW were observed with higher amplitudes in west–east directions and with smaller amplitudes in south–north, even though observations showed a more symmetric AGW distribution. The energy distribution and asymmetry of different ionospheric modes provide important knowledge on seismo-ionospheric interaction and additional information on wave propagation.

The current paper provides the first information on the ionospheric response caused by the 6 February 2023 Turkey–Syria Earthquake. We report an analysis of the ionospheric effects observed with ionospheric sensors, such as dual-frequency GNSS receivers and ionosondes.

2. Materials and Methods

There are different tools to study the ionospheric response to different impacts. Global navigation satellite systems (GNSS) provide sensitive measurements of the ionospheric total electron content (TEC) [9], which is calculated from dual-frequency GNSS phase measurements Iφ.

$$I\phi =\frac{1}{40.308}\frac{f_1^2{f}_2^2}{f_1^2-f_2^2}\left({\lambda }_1{L}_1-{\lambda }_2{L}_2+K+\sigma \phi \right)$$

(1)

where f₁, f₂ are the GNSS operating frequencies, λ₁L₁ and λ₂L₂ are the phase ranges at L₁ and L₂, K is ambiguity, and σφ is TEC noise due to L₁/L₂ phase noises.

Since a pioneering paper by Calais and Minster [4] was published, dual-frequency phase TEC has been widely used to study earthquake-related ionospheric waves [10,11]. Phase TEC series contain a trend caused by GNSS satellite motion, navigation signals travel different paths in the ionosphere as a satellite moves; the longer the path, the larger the TEC that is observed. To exclude this trend and select TEC variations within different bands (2–10 min, 10–20 min, 20–60 min), we used smoothing splines and running average filters, correspondingly, because these procedures showed the best performance to distinguish the effects of traveling ionospheric disturbances from raw TEC series [12].

Besides TEC variations, we used rate-of-TEC index (ROTI) [13,14]. ROTI is a root mean square of the TEC rate within a 5 min interval.

$$ROTI=\sqrt{{\langle \left(\frac{dI\phi }{dt}\right)}^2\rangle -{\langle \frac{dI\phi }{dt}\rangle }^2},$$

(2)

GNSS data could be used to derive coordinates based on precise point positioning solutions [15], which are used to estimate shifts in the Earth’s crust [16], and hence, the magnitude of AGW induced by such motion.

In addition to GNSS, we used vertical sounding ionosondes that probe the ionosphere with a set of frequencies and collect propagated signals that are compiled in an ionogram. For vertical sounding, ionograms represent amplitude–height–frequency characteristics. To obtain the electron density profile, the inverse problem is solved; hence, the ionospheric critical frequency, foF2, is usually more accurately defined than the peak height, hmF2. We used automatic ionogram scaling based on ARTIST-5 software.

The GNSS gives good spatial and time resolutions (the data sample rate is 30-s), while the ionosondes gives better understanding of how ionospheric absolute parameters evolve during earthquakes.

We used the SIMuRG on-line tool (https://simurg.space/ (accessed on 20 April 2023)) for data collection, processing, and visualization [17]. SIMuRG provides TEC-based data products. We involved the International GNSS Service receiver network [18] and Turkish National Permanent GPS Network (TNPGN) [19] for GNSS data and Global Ionospheric Radio Observatory [20] for ionosonde data. Experimental facilities in the region included 296 GNSS receivers (148 of global GNSS network and 148 of TNPGN) and 1 ionosonde (Nicosia) (35.0°N, 33.2°E) (see Figure 1). We also used distant ionosondes Eglin (30.5°N, 86.5°W) and El Arenosillo (37.1°N, 6.7°W) as a reference. The ionogram sample rate was 5 min for the Nicosia site and 15 min for the other sites.

Figure 1 also shows the focal mechanisms for the studied earthquakes according to the United States Geological Survey (for information about the 01:17 UT and 10:24 UT earthquakes, see https://earthquake.usgs.gov/earthquakes/eventpage/us6000jlqa/executive (accessed on 20 April 2023) and https://earthquake.usgs.gov/earthquakes/eventpage/us6000jllz/executive (accessed on 20 April 2023), respectively).

To analyze the velocities of waves, we used distance–time (or travel-time) diagrams [21]. The distance–time diagrams are two-dimensional maps of TEC variations or ROTI vs. earthquake shock time and the distance between an observation (ionospheric pierce point) and the epicenter. When we observed a wave/disturbance with a circular front, we recorded a line on the distance–time diagram. The slope corresponds to the wave/disturbance velocity. We sorted data in such a way that points with higher amplitudes are on top so that we can see the entire effect.

While the distance–time diagram method is well known, we suggest using another technique that gives statistical estimates for propagation velocity. First, we estimated the propagation time and calculated a velocity for each line of sight. Second, we assumed that the observed velocity should fit a normal distribution and we found its parameters.

In order to obtain the data for normal distribution fit, we found the time when individual line-of-sight registered an irregularity. We used the minimum rate-of-TEC observed at individual line-of-sight as a criterion for irregularity presence. Only minima below 0.25 TECu were taken into account. Since we used TEC difference at two successive time samples with a 30-s timestamp, 0.25 TECu corresponded to a rate of 0.5 TECu/min.

We used four satellites that show earthquake effects in the rate-of-TEC: G17, G14, G24, and E08. The data from these four satellites show the effect that was reproduced at different GNSS receivers.

We subtracted 10:35:00 UT from the time when the minimum rate-of-TEC (below 0.25 TECu) was detected and used it as disturbance propagation time from the ionospheric region above the epicenter. We obtained the ionospheric pierce point location for the rate-of-TEC minimum at each line of sight. Then, we calculated the distance between these locations and the ionospheric pierce point above the epicenter at a height of 300 km. The velocity was obtained based on the distance and the propagation time. Then, we fit the velocity distribution by a normal distribution fit.

Geomagnetic activity was calm: Kp was mostly around 1.0 for February 5 and around 3.0 for February 6. To avoid the influence of geomagnetic and solar activities, we analyzed local variability around the times of the shocks, rather than performing a day-to-day comparison.

3. Results

Figure 2 shows snapshots for GNSS data products for several selected times for the 6 February 2023 10:24 UT earthquake. We used 300 km to plot ionospheric pierce points on the maps. The panels show (from top to bottom) ROTI, 2–10 min TEC variations, 10–20 min TEC variations, and 20–60 min TEC variations. Figure 2 shows a well pronounced effect by means of ROTI. The earthquake might have been accompanied with TEC variations with periods of 2–10 min (acoustic mode) and 10–20 min (gravitational mode). The coseismic ionospheric effects were less pronounced in the TEC variations due to “background” TEC variations with similar amplitudes existing at the time. A further analysis was done to check whether an earthquake caused effects in the 2–10 min and 10–20 min bands. It was found that the 20–60 min TEC variations did not produce any noticeable effects, i.e., there were no TEC variations associated with earthquakes above the background or variations with the specific distribution.

We observed the first ionospheric effects 10 min after the shocks (see Supplementary Material File S1). According to the ionosonde data, Nicosia station hmF2 was ~300 km from the epicenter, so we can estimate a vertical propagation velocity of 500 m/s given 10 min (600 s) propagation time. This value corresponds to the average speed of sound Cs = 550 m/s between the Earth’s surface and the ionosphere, i.e., the average of the speed near the Earth’s surface, Cs~300 m/s, and that at a height of 300 km, Cs~800 m/s [22]. It might be that the gravitational branch of AGW is responsible for this.

The maximal amplitudes for ROTI were recorded to the north from the epicenter. The 2–10 min TEC variation data show that ionospheric disturbances were seen farther to the south than to the north. We will show this further using the distance–time method. The amplitude of southward propagated disturbances was smaller than the amplitude of the disturbances recorded to the north from the epicenter.

We discovered that the 10:24 UT earthquake had a shape similar to the ionospheric effects: elongated in longitudinal direction (https://earthquake.usgs.gov/earthquakes/eventpage/us6000jlqa/map (accessed on 20 April 2023)). The disturbed region was split, and at 10:45 UT, we observe two different regions: one the north and another to the south from the epicenter.

Figure 3 compares the amplitudes of ionospheric effects for the earthquakes that occurred at 01:17 UT (upper panels) and 10:24 UT (bottom panels). The panels correspond to the times when the maximal amplitudes of the ionospheric disturbance were recorded. We observe a southward co-seismic disturbance (using 2–10 min TEC variations) propagation for the 01:17 UT earthquake, similar to that of the 10:24 UT earthquake.

The 10:24 UT earthquake was accompanied by much stronger ionospheric effects. For the 10:24 UT shock, amplitudes of 2–10 TEC variations and ROTI reached up to 0.5 TECU and more than 0.5 TECU/min, correspondingly, while values of 0.1 TECU and 0.15 TECU/min were reached for the 01:17 UT shock. Ionospheric effects delayed the shocks by 11 min (observed at 01:28 UT and 10:35 UT).

The effects lasted around 10 min and were observed in the region of (30–40°N, 30–40°E). The direction of the ionospheric disturbance propagation was preferably southward according to 2–10 TEC variations data and northward according to ROTI data. Ionospheric disturbances were observed as far as 750 km from the epicenters. The lack of GNSS receivers in the region compelled us to estimate the horizontal pattern of irregularities. Data to the north show enhanced ROTI values beyond the black sea coastline (see the middle plot for ROTI). Using only the global GNSS receiver network, we observed enhanced ROTI only over land; meanwhile, TNPGN revealed the true horizontal distribution of the irregularity. Enhanced ROTI values were observed in the South and in the North; however, the center of the enhanced ROTI region did not coincide with the earthquake epicenter (see left plot for ROTI). TEC variations with 2–10 min periods showed almost South propagation (see right panels in the Figure 3).

In order to estimate the velocities of ionospheric disturbances, we calculated distance–time diagrams (Figure 4) for the data used in Figure 2. Each snapshot (like shown in Figure 2) provides a single column of data in Figure 4.

Figure 4 also shows that the first manifestation of the ionospheric effects appeared at 10:35 UT (~10 min after the 10:24:50 UT earthquake) and lasted about 10 min. Increased ROTI values were seen as far as 750 km from the earthquake location. The earthquake effect was isolated from other regions with higher ROTI. Horizontal lines of increased ROTI values were due to geostationary satellite data that appeared at the same location for any epoch. We plotted slant lines to estimate the velocity and obtained values of 2000 m/s for ROTI data and 900–1500 m/s for 2–10 min TEC variations. There was no noticeable pattern in the 10–20 min TEC variations that allows velocity estimation. This information allowed us to conclude that there were no co-seismic disturbances in the 10–20 min band.

We estimated the velocity of vertical propagation as 500 m/s, given that the ionospheric maximum was at 300 km. The middle panel of Figure 4 shows the horizontal propagation of ionospheric disturbance using 2–10 min TEC variations. We attribute this disturbance to the acoustic branch of the AGW, because its period was estimated as 5–6 min (period is estimated using distance between oblique lines) and its velocities as 900 m/s and higher. This velocity corresponds to the speed of sound at ionospheric maximum heights. Acoustic waves have frequencies higher than the acoustic cutoff frequency ωa or periods less than the acoustic cutoff period, Ta, that is, less than ~10 min in the ionosphere.

The distance–time diagram for 2–10 min TEC variations shows that the 01:17 UT (Figure 5) and the 10:24 UT (Figure 4) earthquakes were accompanied by ionospheric disturbance with velocities of 900–1500 m/s. We see less pronounced effects compared to the 10:24 UT quake (see middle panel in Figure 4), which might be connected to the lower electron density during nighttime, as well as shorter delay between tremors (~2 min for 01:17 UT earthquake and ~10 min for 10:24 UT according to ANTO seismometer). We changed the color scale to make the effect visible. The ROTI data for the 10:24 UT earthquake show a velocity of 2000 m/s. The ROTI data for the 01:17 earthquake did not show any pattern that allowed us to estimate the velocity, and a plot of these data is not provided in the current paper. Ionospheric disturbance by means of 2–10 min TEC variations showed that the velocity evolved from 1500 m/s to 900 m/s after one wavelength had passed (see dashed, solid, and dotted slant lines in the middle panel in Figure 4 and the bottom panel in Figure 5).

In Figure 6, we present the result of the alternative velocity estimation method for co-seismic disturbance, as discussed in Section 2. The velocity values as bins and the black curve show the normal distribution with a standard deviation of 420 m/s and mean values of 2027 m/s. Our statistical estimate agrees with the velocity estimation based on the distance–time diagram.

The velocities correspond to Rayleigh mode [2,23] by means of ROTI and acoustic mode by means of 2–10 min TEC variations. This indicates that the ionospheric disturbances may have been caused by different sources: (a) by an acoustic wave generated either from the propagating Rayleigh waves or (b) from the ground displacement at the epicenter. Velocities of 1000–1500 m/s of co-seismic TEC disturbances are very typical for earthquakes in Turkey [2,24]. These waves are assumed to relate to direct acoustic waves from the focal area. Rolland et al. [24] confirmed this assumption with 3D modeling. Thus, the 2–10 min TEC variations were caused by an acoustic wave generated at the epicenter. ROTI revealed disturbances caused by the Rayleigh wave.

We also check the north–south asymmetry of the 2–10 min TEC variations using distance–time diagrams, as shown in Figure 5. We plotted a distance–time diagram for data from the north of the epicenter (top panel) and another diagram for the data from the south (bottom panel). We see that the 2–10 min TEC variations were mostly observed in the South. We see the same behavior for the 10:24 UT earthquake: the propagation direction was mainly southward.

The ROTI data showed more symmetric response for the 10:24 UT earthquake, the effects of which were seen to the north and to the south of the epicenter.

From Figure 5, we can estimate the period for ionospheric irregularities as 5–6 min for the 01:17 UT earthquake. The same period was observed for the 10:24 UT earthquake (see middle panel in Figure 4). The estimation was done based on the distance between the dashed and dotted oblique curves. The same period estimate holds for the 10:24 UT earthquake.

ROTI behaved differently, and the hypothesis is that it was driven by irregularities of smaller scale. Smaller scale irregularities manifest themselves as fast TEC changes, so the filtration smooths them out. To examine this, we considered individual receiver-satellite lines-of-sight. Each line-of-sight provides a TEC time series. To reduce the influence of different sounding geometries, we selected a single satellite, E08, whose data showed the maximal number of disturbed series and corresponded to the area to the north of the 10:24:50 epicenter. Data for E08 are shown in Figure 7.

We calculated the dTEC values by subtracting the TEC from its previous value at every epoch. To remove the influence of linear trend of the TEC, we subtracted the average value from the corresponding dTEC series: a linear term in TEC corresponds to a constant term in dTEC. We used a 0.5 TECu threshold to limit the disturbed series number. The series were ordered by the time of first maximum occurrence. We also visualized TEC series with the linear trend removed to see the signature of the disturbance in the TEC series. The ROTI data are given for reference to estimate how long a given particular feature in the TEC series can persist in ROTI. The data are presented as three panels in Figure 7. We see that typical dTEC had enhanced values with a time range less or equal 2 min. This explains the different effects shown by ROTI and the 2–10 min variation data. ROTI can catch up to faster changes of TEC than 2–10 min variations. Restored TEC showed the form of the ionospheric response. The quadratic trend remained in the data, but it did not prevent us from identifying a response shape.

Ionosonde data did not show any significant foF2 or hmF2 changes or variations for the earthquake day (Figure 8). We did not observe the significant wave-like variations that we would expect in foF2 or hmF2. On the ionograms, we also did not notice any features that would be typical for traveling ionospheric disturbances.

Figure 8 demonstrates a significant increase in the daytime foF2 for the earthquake day compared to previous days. However, this was caused by increased solar flux. The F10.7 index from the OMNIWeb service showed an increase of F10.7 flux (50 units for 3 days) starting on February 6.

We applied the same trend extraction procedure to foF2 and hmF2 data as we did for the TEC series. The same band filtration procedure was not available due to the different temporal resolutions and lower sampling rate. Detrended data still showed no visible effects according to foF2 and hmF2 data.

We also studied the GNSS receiver coordinates based on precise point positioning solutions. We calculated positioning errors in the X (east), Y (north), and Z (up) directions for every time step. The data were subjected to the distance–time method described above. The results are presented in Figure 9. Each receiver contributed a single row to each panel. We did not observe an increase in positioning errors linked to earthquake times. There were a couple of receivers close to the epicenters (within 100 km) that stopped operating soon after the earthquakes.

4. Discussion

Based on GNSS and ionosonde data, we studied the ionospheric response to the earthquakes that occurred in Turkey on 6 February 2023. There were two major shocks: the first strong earthquake (with magnitude M7.8) occurred at 01:17:34 UT and the second one (with magnitude M7.5) occurred at 10:24:50 UT. The magnitude of both earthquakes exceeded the known threshold of M6.5, below which there were no pronounced earthquake-induced GNSS-TEC disturbances [7]. On the other hand, both shocks had a focal mechanism corresponding to strike–slip faulting, i.e., there were no large vertical displacements of the ground during the earthquakes. Such earthquakes may not generate ionospheric disturbances. However, Astafyeva et al. [25] showed that strong strike–slip earthquakes can also produce TEC disturbances with sufficiently big amplitudes. Thus, we would expect that both major earthquakes would cause noticeable disturbances in the ionospheric GNSS-TEC.

GNSS data analysis indeed revealed noticeable effects in ROTI and 2–10 TEC variations after the first strong (01:17:34 UT, M7.8) shock and after the second one (10:24:50 UT M7.5). The average period of TEC disturbances was ~ 5 min. This means that the TEC disturbances could have been caused by an acoustic wave generated in the epicenter. According to the authors of [26], acoustic waves have frequencies ω > ωa and periods T < Ta = 2π/ωa, where ωa= g·γ/2Cs is the acoustic cutoff frequency, g ≈ 9.8 m/s² is the acceleration of gravity, Cs is the speed of sound, and γ is the ratio of specific heats. In the ionosphere at altitudes higher than 200 km [22], Cs ≈ 800 m/s and γ ≈ 1.67, so Ta ≈ 10.2 min. Thus, the TEC disturbances with periods of ~ 5 min corresponded to acoustic waves.

The amplitude of 2–10 min TEC disturbances for the 10:24 UT earthquake varied from 0.1 to 0.5 TECU. These values were in good agreement with the results of other researchers [2,24,25]. In particular, Astafyeva et al. [25] showed that strike–slip earthquakes with magnitudes from M7.5 to M8.0 can generate TEC disturbances with average amplitudes from 0.2 to 1.2 TECU.

We observed that the 10:24 UT earthquake was accompanied by much stronger ionospheric effects compared with the 01:17 UT earthquake. This might be explained by the fact that the 01:17 UT earthquake corresponded to a nighttime ionosphere with a background TEC lower than that at 10:24 UT. Another possible cause may be neutral wind, the direction of which during the day differs from the direction during the night. It is possible that at night, the wind weakens the AGWs and prevents their propagation up to ionospheric heights. Yet another reason might be connected with the earthquake. According to ANTO seismometer data, the nighttime earthquake had a shorter delay between tremors (~2 min) than the 10:24 UT earthquake (~10 min).

The observed 10-min time delay between the shocks and the ionospheric response agree with the results reported previously. Astafyeva et al. [3] showed a ~500 s time delay between earthquake time and TEC variation. The TEC variations considered by Astafyeva et al. [3] were obtained by removing trends with a 20 min running average. Calais and Minster [4] considered 3–10 min TEC variations (“band pass filtered between periods of 10 and 3 min”) and found that the amplitude of earthquake-related variations increases 10–30 min after the shock. They considered several lines-of-sight, since the receiver network was not dense. Some of the lines-of-sight were farther from the epicenter than others. A slower arrival time corresponds to more distant satellite-receiver lines-of-sight. However, the lower limit, 10 min, corresponded to our results.

Liu et al. [1] studied the M8.0 Wenchuan Earthquake and used TEC variations with periods greater than 10 min. Meng et al. [27] performed a simulation to study electron density perturbations caused by earthquake-induced AGW. Their simulation for the Tohoku earthquake showed that 9 min are enough for a disturbance to reach 200 km and produce 20% electron density variations that would be seen in the TEC variations.

Our results showed that ionospheric disturbances propagated more to the south by means of 2–10 min TEC variations. Based on ROTI data, the disturbance was observed both to the south and to the north of the 10:24 UT earthquake epicenter. We attributed the disturbance shown by ROTI to the disturbance generated by Rayleigh waves. There was no other evidence, however, that the velocity of 2000 m/s significantly exceeded the speed of sound in the ionosphere, i.e., at an altitude of 300 km, Cs ~800 m/s [22]. Therefore, it is assumed that ionospheric disturbances with velocities of more than 2000 m/s were caused by surface Rayleigh waves propagating from the epicenter [28,29,30]. This assumption was confirmed by numerical modeling [30]. So, we can attribute the 2000 m/s disturbance to Rayleigh waves.

Yasyukevich et al. [8] showed southward propagation of Rayleigh-mode-associated TEC disturbance for the Tohoku earthquake. The asymmetry could be connected with the particular pattern of vertical displacement of the Earth’s surface, i.e., larger ground vibrations cause larger ionospheric effects. The southward propagation of co-seismic ionospheric disturbances was previously detected for other earthquakes in Turkey [2,24].

For atmospheric waves, the asymmetry can be explained by the influence of the geomagnetic field. AGWs generate ionospheric disturbances due to collisions of neutral and charged particles. The charged particle motion across the geomagnetic field line is restricted. In the direction perpendicular to the geomagnetic field, charged particles cannot move together with neutral ones, and ionospheric disturbances will be suppressed. According to the IGRF-13 model [31] (https://www.ngdc.noaa.gov/geomag/magfield.shtml (accessed on 20 April 2023)), geomagnetic declination and inclination for the area of Turkish earthquakes are 5.8° and 55.2°, respectively. This means that southward AGWs propagate parallel to the geomagnetic field, while northward, eastward, and westward AGWs propagate perpendicular to the geomagnetic field. Thus, the geomagnetic field permits the southward propagation of ionospheric disturbances. Based on 3D modeling, Rolland et al. [24] reproduced such asymmetry of ionospheric disturbances during the 2011 Van earthquake, Turkey. Furthermore, they showed that the spatial distribution of ionospheric disturbance amplitude is also controlled by the geomagnetic field. For the ionospheric disturbances associated with the Rayleigh waves, the most probable explanation is in the peculiarities of Rayleigh waves on the Earth’s surface [8].

The Nicosia ionosonde did not show wave-like peculiarities in the foF2 or hmF2 data. The critical frequency, foF2, is a local parameter rather than an integral (like TEC), and we would expect a more profound effect for ionosondes, because we avoided smoothing by integrating data from different ionospheric regions. However, we had lower time resolution for the ionosonde data (5 min) and most probably could not catch the fast processes that were recorded by the GNSS network.

We found that on the earthquake day and the days after, foF2 was higher for daytime than for previous days. Sometimes, changes in general ionosphere dynamics are considered to be earthquake related [32,33]. We analyzed the foF2 time series retrieved from two distant ionosondes Eglin and El Arenosillo stations. While there was a significant difference in the foF2 dynamics and values for different ionosondes, the same peculiarities (an increase in daytime values) were observed at distant ionosondes, so the effects seemed to be a global one (connected to increased solar flux) and not linked to the earthquake. F10.7 increased from 140 s.f.u. on February 5 and 145 s.f.u. on February 6 to 200 s.f.u. on February 9.

5. Conclusions

Two strong earthquakes occurred in Turkey and Syria on 6 February 2023, with magnitudes of M7.8 and M7.5—one at nighttime and the other in the daytime. Both generated disturbances in the ionosphere, as detected by GNSS data. The amplitude of the ionospheric response for the daytime, M7.5 earthquake exceeded those of the nighttime M7.8 earthquake by five times, i.e., 0.5 TECU/min and 0.15 TECU/min according to the ROTI data. The velocities of the earthquake-related ionospheric waves were ~2000 m/s according to the ROTI data. We also employed a method that allowed us to estimate the velocity of the ionospheric disturbances using different approach than the distance–time method. For now, this method uses the time when the first effect is observed in the ionosphere by means of rate-of-TEC data to select the effect. Further work could be done to avoid these constraints. The results of our method agree with our velocity estimation using the distance–time method. A velocity of 2000 m/s corresponds to the ionospheric disturbances caused by Rayleigh mode.

TEC variations with 2–10 min periods showed velocities from 1500 to 900 m/s. This corresponds to the acoustic mode. We can conclude that the disturbances caused by different modes had different time scales. The Rayleigh mode TEC response is shorter, i.e., around 2 min. Acoustic mode is accompanied with TEC variations with 5–6 min periods. Ionospheric disturbances propagated farther to the south.

The ionospheric effects were recorded as far as 750 km from the epicenter. Ionosonde data located 420/490 km from the epicenters did not catch ionospheric effects. We did not notice any acoustic gravity mode of atmospheric waves.

We observed that the irregularities developed around the epicenter of the 10:24 UT earthquake and split into two different regions around 10:42 UT. At 10:45 UT, two different regions with enhanced ROTI could be well distinguished. The geometry of the earthquake affects the geometry of the ionospheric effect: longitude elongation of ionospheric effect occurred, region of enhanced ROTI values splitting into two regions to the north and to the south of the epicenter.

We observed a symmetric ROTI response to the 10:24 earthquake and an asymmetric response by means of 2–10 TEC variations. This agrees with the simulations performed by the authors of [24]. The 2–10 TEC variations were caused by the acoustic mode of AGW, and the geomagnetic field allowed only southward propagation to occur. The enhanced ROTI was caused by the upward propagation of AGW caused by the Rayleigh wave that travelled into the Earth’s crust. Since the geometry of the geomagnetic field was nearly the same around the epicenter, the ionospheric response was expected to be symmetric; the ROTI data confirmed this.