Anisotropy of the Near‐Field Coseismic Ionospheric Perturbation Amplitudes Reflecting the Source Process: The 2023 February Turkey Earthquakes

The East Anatolian Fault in southern Turkey ruptured on 6 February 2023, causing a Mw 7.8 earthquake (EQ1), one of the largest strike‐slip events recorded on land. ∼9 hr later, earthquake of Mw 7.7 (EQ2) occurred to the north of EQ1. We investigate here near‐field coseismic ionospheric perturbations (CIP) caused by acoustic waves (AWs) excited by coseismic vertical crustal movements. We find that observed CIP periods were somewhat longer for EQ1 than EQ2. EQ1 also showed azimuthal dependence in CIP amplitudes that cannot be explained by known factors such as geomagnetism and line‐of‐sight geometry. Numerical experiments revealed that CIP by EQ1 can be well reproduced by assuming a suite of sources along the fault that successively ruptured. Small but significant dependence of amplitudes and periods on azimuths were caused by interference of AWs from multiple sources. We also found that CIP amplitudes of strike‐slip earthquakes tend to be lower than dip‐slip earthquakes.

Plain Language Summary Fault dislocations in earthquakes cause vertical movement of the surface and excite acoustic waves (AWs) in the overlying atmosphere which propagate upward with increasing amplitudes.Such waves reach the ionosphere and disturb electron density distribution there, causing disturbances in numbers of electrons along the line-of-sights connecting ground GNSS receivers and satellites.So far, such near-field co-seismic ionospheric perturbations are modeled by assuming single acoustic pulse from the surface, although large earthquakes often involve ruptures of multiple fault segments spanning hundreds of kilometers.Here we demonstrate that interference of AWs from these multiple sources makes differences in the perturbations amplitudes and periods at GNSS stations in different azimuths from the epicenter.
• Mw 7.8 and Mw 7.7 February 2023 Turkey earthquakes produced typical ionospheric perturbations • For the first earthquake, we found small but significant difference in the perturbation amplitudes that cannot be explained by known factors • By assuming multiple sources along the fault and interference of acoustic waves from them, we can explain the observed azimuthal asymmetry

Supporting Information:
Supporting Information may be found in the online version of this article.
and demonstrated that the sum of CIP from those sources explains complicated waveforms of the observed CIP.Based on the azimuthal distribution of near field CIP, Bagiya et al. (2020) identified discrete seismic sources along the Tohoku 2011 rupture during the initial 60s of the event.
Here we analyze near-field CIP of EQ1.We demonstrate, for the first time, that CIP generated by relatively small (<Mw 8) earthquakes could also contain contributions from a suite of sources along the fault, and that CIP shows variety of amplitudes and periods for different satellite-station pairs due to combinations of sub-CIPs from multiple sources with different time lags.Further, we show that the CIP of EQ2 has much larger amplitude and slightly shorter period than EQ1 and explain these differences assuming a single source and higher background ionospheric electron densities.

Seismic Sources
A lithospheric block bounded by EAF to the south and the North Anatolian Fault to the north, referred to as the Anatolian block, is extruded westward causing left-lateral movement on these faults.EQ1 and EQ2 occurred as strike-slip earthquakes on the southwestern segment of EAF. Figure 1a shows their epicenters together with the three sub-segments that ruptured during the two events.The source time functions by US Geological Survey (2023) (Figure 1a) revealed that moment release occurred over periods of ∼130 and ∼70 s for EQ1 and EQ2, respectively.Figure 1b shows the slip distributions of segments 2 and 3 of EQ1 and segment 1 of EQ2 from US Geological Survey (2023).
The slip distribution maps also provide contours indicating times of dislocation (US Geological Survey, 2023).
The source time functions in Figure 1a suggest that the energy release took place in at least two episodes during EQ1 while it was mainly released in a single episode during EQ2.Colored stars in Figure 1b are the sources supposed to have contributed to generate AWs and consequent ionospheric perturbations.Actually, such excitation occurred continuously in space rather than at three distinct locations, but we simplified the situation by representing the whole process with just three sources in performing numerical experiment explained in Section 3.

GNSS-TEC Data
The CIP generated in the near-field region (500-600 km of the seismic source) are well captured using GNSS as changes in total electron content (TEC) (Bagiya et al., 2019 and references therein).In the present study, we used data from GNSS stations in Turkey, Israel, Cyprus, and Iraq (Figure S1 in Supporting Information S1).Technical details on derivation of TEC along the line-of-sight (LOS) from a satellite to a GNSS receiver can be found in recent review articles, for example, Astafyeva (2019) and Heki (2021).Here we use TEC derived from the phase measurements.In order to highlight CIP, TEC time series in Figure 2 were detrended using 2-8 mHz bandpass filter.
The black dotted circle in Figure S1 in Supporting Information S1 shows the region within ∼600 km from the EQ1 epicenter.TEC represents electron density integrated along a LOS.In GNSS-TEC studies, we often calculate its ionospheric pierce point (IPP), assuming a thin layer at altitude 300 km, and let them represent the measurement points (Figures 2c and 2d).
Space weather affects the global ionosphere on various spatial and temporal scales.The geomagnetic storm condition during EQ1 and EQ2 was examined using two indices.The SYM-H index shows the measure of the storm time ring current.The SYM-H index values were −26 and −32 nT respectively at the occurrence times of EQ1 and EQ2.Also 3-hr Kp-index, reflecting disturbances in mid-latitude regions, were less than 5 then.These indices suggest calm geomagnetic activities (https://www.swpc.noaa.gov/products/planetary-k-index).

Observed CIP
Clear CIP for both EQ1 and EQ2 are recorded at GNSS stations, mainly located to the south of the epicenters.They are shown as band-pass filtered time series in Figure 2. Alphabets before satellite numbers indicate different GNSS, that is, G, R, E, and C for GPS, GLONASS, Galileo, and BDS, respectively.Sampling intervals are 30 s, except for those shown in Figure 3 (1 s).GPS satellites G31 and G28 recorded significant CIP during EQ1.

Synthetic CIP
Slant TEC (STEC) change time series for a given satellite-station pair by an AW from a point source on the surface can be synthesized by the following three steps (Heki & Fujimoto, 2022).(a) Ray tracing of AW propagating upward in various launch angles over the source, (b) calculating ionospheric electron density anomalies in three-dimensional (3D) space filled with blocks, (c) integrating electron density along the real GNSS satellitestation LOS penetrating the 3D blocks.By repeating (b) and (c) for various epochs, we could synthesize STEC curves.
In step (a), we used the AW speed from US standard atmosphere 1976.Effects of geomagnetic field and the LOS geometries are considered in steps (b) and (c), respectively.The inclination and declination values at IPP altitude of 300 km are obtained from IGRF-13 model (Alken et al., 2021).An example of the steps (a) and (b) is given for three epochs in Figure S5 in Supporting Information S1.There, we assumed an N-shaped disturbance with a period of 2 minutes.We also assumed the Chapman distribution of ionospheric electron density profile at the peak height of 300 km.Because of the northward/downward geomagnetic fields, electron density anomalies are much suppressed at the northern side of the epicenter, causing strong north-south asymmetry.In (c), we calculated time-variable geometry of LOS of specific satellite-station pairs and integrated electron density anomalies along LOS every 15 s to draw synthetic STEC time series.Examples for EQ1 and EQ2 are shown in Figure 3. Bagiya et al. (2019) summarized non-seismic forcing mechanisms affecting CIP amplitudes into three, that is, geomagnetic field-neutral wave orientation factor (GNF), satellite geometry factor (SGF), and electron density factor (EDF).In the numerical synthesis presented here, GNF and SGF are considered in steps (b) and (c), respectively.Figures S6a and S7 in Supporting Information S1 show these factors for the present case.In our numerical simulation, we gave an arbitrary scaling factor to adjust synthetic curve amplitudes to those observed.
In Figure 3 we consider the difference in EDF using a factor of 4 when we compare simulation results of EQ1 and EQ2.EQ2 occurred in the local noon time while EQ1 occurred in the early morning time.Global ionospheric maps indicate that vertical TEC for EQ2 is ∼4 times as large for EQ1 (Figure S6b and S6c in Supporting Information S1).

Azimuth Dependence of CIP Amplitudes
CIP amplitudes depend on the magnitude and mechanism of the earthquake, that is, a larger dip-slip earthquake makes larger CIP (Cahyadi & Heki, 2015).They are also influenced by the three non-seismic factors, geomagnetic field, LOS geometry and background electron density, as discussed in the previous section.However, here we point out that the variations in CIP amplitudes by EQ1 cannot be explained by such known factors.Figure S8a in Supporting Information S1 compares observed and synthesized STEC changes for three different satellite-station pairs.There, the IPPs of the three cases are at similar distances from the epicenter but in three different azimuths, southwest (nico-E09), south (bshm-R04), and southeast (bshm-C16).The latter observations show ∼1.7 times as large CIP amplitudes as the former.However, synthesized TEC changes do not reproduce such large differences.This is also supported by spatial distributions of the factors for geomagnetism (GNF) and LOS geometries (SGF).Figures S6a and S7 in Supporting Information S1 indicate that these two factors did not show any substantial differences for these three satellites (Figures S6a and S7 in Supporting Information S1), that is, these non-seismic factors are not responsible for the observed amplitude variations.These arguments suggest that some unknown factor(s) controls CIP amplitudes.We speculate it is the interference of AWs from multiple sources and verify this using numerical simulations.
Figure 1 shows that the rupture started in segment 1 (S1 in the figure) and propagated to segments 2 and 3 by ∼3.2 km/s (Melgar et al., 2023).AWs might have been excited continuously from sources moving together with rupture along the fault.Although the entire process is continuous, we represent it with three distinct points shown as AW-sources1, 2, and 3 which are the centroids of large slip zones.Slips in segment 1 are small, and we neglect them.AW-source1 corresponds to the largest dislocation patch of segment 2 with the largest slip occurring ∼20 s after the rupture onset.Sources2 and 3 are assumed to have generated AWs simultaneously ∼40 s after the rupture onset.These two sources correspond to the later peak in the source time function of EQ1 (Figure 1a).Here we assume that the amplitudes of AWs from sources 2 and 3 are 2/3 of that from AW-source1 considering the ratios of total seismic moments represented by the three sources.
In Figure 3a, we assumed that the nico-E09 pair recorded CIP as the sum of the sub-CIPs from these three sources, shown in three different colors with time lags.Considering the small CIP amplitudes (<1% of the background TEC), we assumed the linearity.The time lags reflect the rupture time lag and the propagation time lag.In this case, the former is ∼20 s, and the latter may vary from 0 to ∼90 s depending on the azimuth of IPP (Figure S9 in Supporting Information S1).In Figure 3a case, the IPP is located to the southwest of the epicenter, that is, the western extension of the ruptured fault, and this maximized their propagation time lags.Here, sub-CIP by AW-source3 is generated 20 s later than AW-source1, it arrives at the IPP earlier by ∼70 s (90-20 s).On the other hand, sub-CIP from AW-source2 arrives much later than the other two and does not contribute to increase the total amplitudes.Such time lags contribute to make the total CIP period longer without increasing the total amplitude relative to individual sub-CIPs.
Figure 4 shows two examples, bshm-R04 and bshm-C16.Their IPPs are located south and southeast of the epicenter, respectively, and CIPs from sources 1-3 arrive there with less time lags than Figure 3a case, increasing the amplitude of the total CIP.Although single-source model cannot explain the larger CIP amplitudes for these two cases (Figure S8 in Supporting Information S1), the three-source model can well reproduce the factor ∼2 difference in the amplitude (Figure 4).
Regarding EQ2, its CIP for nico-E02 (Figure 3c) can be well simulated using a single source excited 10 s after the rupture onset along segment1, reflecting its compact source time function (Figure 1a).We neglect small peaks in the source time function in our simulation.Because the spatial coverage of the EQ2 slips is smaller, assuming an additional source for the secondary peak would not make a delay in AW arrival large enough to influence the waveform and the amplitude.Seismic moment by EQ2 is comparable to that represented by AW-source1 of EQ1.However, EQ2 occurred near the local noon, and background vertical TEC showed ∼4 times as high background TEC as EQ1 (Figures S6b and S6c in Supporting Information S1).Hence, we assumed 4 times larger AW amplitude for EQ2 than the AW-source1 of EQ1.This well explains the observed CIP amplitude in Figure 3c.Period of the EQ2 CIP (Figure 3a) is ∼2.7 min, shorter than ∼3.4 min for EQ1 (Figure 3c).This relatively short period of CIP can also be understood as one of the single source CIP characteristics.

Concluding Remarks
The goal of this study is to emphasize that the earthquake source process also influences the amplitudes and periods of CIP for relatively small earthquakes (<Mw 8), together with known factors such as geomagnetism and LOS geometry.On the other hand, we do not attempt to reproduce exactly the same CIP amplitudes and periods for various satellite-station pairs with numerical simulations.Hence, our assumptions are gross, and our simulation results contain uncertainties as explained below.
At first, we modeled the continuous AW excitation due to long rupture with only three sources.This might be an oversimplification, but two prominent peaks in the EQ1 moment release curve (Figure 1a) may justify it to some extent.The second point would be the difference of locations of surface uplift and fault slips.We assumed the positions of AW sources at points of large slips.However, in strike-slip earthquakes, a pair of uplift and subsidence occurs where the slip vectors change their lengths along the fault, for example, at the termination points of the rupture.This may result in errors in the AW sources by up to a few tens of kilometers.Such errors may be responsible for small differences of arrival times between observed and synthesized CIP seen in Figures 3 and 4.
As the third point, our numerical simulation does not include estimation of absolute electron density anomalies.Hence, disturbance amplitudes of the synthesized TEC curves are meaningful only in relative sense, for example, between EQ1 and EQ2, and between different satellite-station pairs.More realistic physical approaches, including the quantification of atmospheric density perturbation caused by surface uplift, would let us synthesize CIP having absolute amplitudes (Carbone et al., 2021;Gao et al., 2023).
Though the initial strong positive anomalies of CIP for the bshm-R04 and bshm-C16 pairs were effectively reproduced by introducing three wave sources (Figure 4), the subsequent decreases were not clear in the observed CIP.In the N-shaped disturbance we assumed from individual AW-sources, the first half (positive part) represents the compression of neutral atmosphere by crustal uplift.However, the second half (negative part) are weaker, and we would need more precise modeling to well reproduce them with simulations.
At last, we compare the CIP amplitudes of the strike-slip earthquakes, that is, the two Turkish earthquakes studied here and the north Sumatra 2012 doublet events, with those in dip-slip earthquakes.Figure S10 in Supporting Information S1 suggests that amplitudes of the disturbances by strike-slip events are somewhat lower than dip-slip earthquakes by a factor of ∼1/2.This negative shift might reflect smaller vertical surface movements by strike-slip faulting than dip-slip faulting.However, it is difficult to extract a statistically significant conclusion with the current data set, that is, we need further cases to detect a significant difference between the two groups.
are shown in Figure2cand FigureS2in Supporting Information S1, respectively.Because the earthquake occurred in the northern hemisphere, CIP amplitudes are larger on the southern side of the epicenter.GNSS satellites passing just above the epicenter or to the north do not capture clear CIP signatures.They are influences of non-seismic forcing mechanisms caused by geomagnetism and satellite geometry and are discussed in the next section.For EQ2, many GNSS satellites (i.e., G17, G24, G19, E02) captured CIP.Their IPP tracks are shown in Figure2dand partly in FigureS3in Supporting Information S1.Again, we find stronger CIP signals from the south.Clear CIP do not show up in regions farther than ∼600 km from the epicenter.Their propagation velocities, shown in the hodochrone plots in FigureS4in Supporting Information S1, confirm their AW origin.

Figure 1 .
Figure 1.(a) Map showing the East Anatolian Fault, North Anatolian Fault and locations of EQ1 and EQ2 epicenters and their focal mechanisms.Both events ruptured three segments, shown as S1, S2, and S3 (segment 1, 2, and 3).Moment rate functions for both events (US Geological Survey, 2023) are also shown.(b) Finite fault models from US Geological Survey (2023).The contours in each plot show the rupture propagation in second.Significant slip occurred on S2 and S3 in EQ1 and on S1 in EQ2.Colored stars in panels are the acoustic wave sources assumed in this study to generate synthetic ionospheric perturbations.

Figure 2 .
Figure 2. Near field ionospheric responses to (a) Mw 7.8 EQ1 and (b) Mw 7.7 EQ2.Red lines show their rupture onset times.(c, d) Show ionospheric pierce point tracks for time windows shown in panels (a, b), respectively.Time series and tracks are labeled with GNSS station-satellite pairs.Each track is marked with star showing rupture start times (red-EQ1, yellow-EQ2).Dotted circle shows the ∼600 km distance from the EQ1 epicenter.

Figure 3 .
Figure 3. Synthetic coseismic ionospheric perturbations (CIP) for the nico-E09 and nico-E02 pairs for the (a) EQ1 and the (c) EQ2 compared with observed Slant TEC time series.For EQ1, we approximated the entire rupture by three discrete sources of acoustic waves (AWs), AW-source1, AW-source2, and AW-source3, whose positions are shown in panel (b).They were assumed to have ruptured 20 s (AW-source1) and 40 s (AW-source2 and AW-source3) after the rupture onset.For more information, please refer finite fault models shown in Figure1b.The AW-source2 and 3 are assumed to have an amplitude 2/3 of AW-source1.The sum of the synthesized disturbances explains the longer period of the total CIP than the individual sub-CIPs.Ionospheric pierce point track of E09 is for time window of 0.5-2.0UT, same as (a).Black disks along the track show hourly time marks, that is, 1.0 UT and 2.0 UT.For EQ2, considering just a single peak in its moment release (Figure1a), we assumed just one source (d) that ruptured 10 s after the rupture onset.

Figure 4 .
Figure 4. Synthetic coseismic ionospheric perturbations (CIP) for (a) bshm-R04 and (b) bshm-C16 pairs for EQ1.Similar to Figure 3a, we approximated the entire rupture using three discrete sources of acoustic waves (AWs), AW-source1, AW-source2, and AW-source3 in panel (c).Unlike the nico-E09, shown as black curves, sub-CIPs from the three sources have less time lags, and the total CIP have significantly larger amplitudes.Ionospheric pierce point tracks are between 0.5 UT-2.0 UT.Black disks along the tracks show hourly time marks, that is, 1.0 UT and 2.0 UT.