# Efficiency of Multi-Terminators Method to Reveal Seismic Precursors in Sub-Ionospheric VLF Transmitter Signals: Case Study of Turkey–Syria Earthquakes Mw7.8 of 6 February 2023

Mohammed Y. Boudjada, Patrick H. M. Galopeau, Sami Sawas, Giovanni Nico, Hans U. Eichelberger, Pier F. Biagi, Michael Contadakis, Werner Magnes, Helmut Lammer, Wolfgang Voller

► To cite this version:

Mohammed Y. Boudjada, Patrick H. M. Galopeau, Sami Sawas, Giovanni Nico, Hans U. Eichelberger, et al.. Efficiency of Multi-Terminators Method to Reveal Seismic Precursors in Sub-Ionospheric VLF Transmitter Signals: Case Study of Turkey–Syria Earthquakes Mw7.8 of 6 February 2023. Geosciences, 2025, 15 (7), pp.245. 10.3390/geosciences15070245. hal-05297673

HAL Id: hal-05297673

https://hal.science/hal-05297673v1

Submitted on 5 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Efficiency of Multi-Terminators Method to Reveal Seismic Precursors in Sub-Ionospheric VLF Transmitter Signals: Case Study of Turkey–Syria Earthquakes Mw7.8 of 6 February 2023

Mohammed Y. Boudjada 1,* , Patrick H. M. Galopeau 2, Sami Sawas 3, Giovanni Nico 4, , Hans U. Eichelberger 1, Pier F. Biagi 5, , Michael Contadakis 6, Werner Magnes 1, , Helmut Lammer 1 and Wolfgang Voller 1

<figure><img src="image_2.png" /></figure>

check for
updates

Academic Editor: Dimitrios

Nikolopoulos

Received: 24 March 2025

Revised: 19 June 2025

Accepted: 21 June 2025

Published: 1 July 2025

**Citation:** Boudjada, M.Y.; Galopeau, P.H.M.; Sawas, S.; Nico, G.; Eichelberger, H.U.; Biagi, P.F.; Contadakis, M.; Magnes, W.; Lammer, H.; Voller, W. Efficiency of Multi-Terminators Method to Reveal Seismic Precursors in Sub-Ionospheric VLF Transmitter Signals: Case Study of Turkey–Syria Earthquakes Mw7.8 of 6 February 2023. *Geosciences* **2025**, *15*, 245. https://doi.org/10.3390/geosciences15070245

**Copyright:** © 2025 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license
(https://creativecommons.org/
licenses/by/4.0/).

## Abstract

This work presents an analysis of the sub-ionospheric VLF transmitter signal disturbances which were detected more than one week before the Turkey–Syria EQ occurrence. We have applied the multi-terminator method when considering amplitude and phase variations of the TBB transmitter signal (Turkey), selected because of a good signal to noise ratio for the amplitude, a stable phase variation, and a ray-path propagation crossing the pre-seismic sensitive region, estimated from the combination of the Dobrovolsky area and the Fresnel zone. New spectral features, i.e., inflexions and jumps, are considered in this study, besides the minima and maxima investigated in. The spectral occurrence probabilities are derived at three specific locations: Graz facility, TBB station and EQ epicenter. We show that two main precursors occurred from 27 to 30 January, and from 31 January to 3 February. More important are the prior precursors detected from 23 January to 25/26 January, where anomaly fluctuations were found to be similar to those at the EQ epicenter area, approximately. A forecasting model is proposed, in which the main steps can provide, in the presence of spectral anomalies, first hints regarding the longitudinal locations of the seismic preparation zone.

**Keywords:** radio earthquakes precursors; sub-ionospheric VLF propagation; Turkey–Syria Earthquake 6 February 2023

## 1. Introduction

The pioneer investigations of Gokhberg et al. [1] and Warwick et al. [2] motivated and boosted interest in the study of the origins of seismo-electromagnetic phenomena. Main reviews on this subject have been reported in [3–6]. Four aspects are introduced hereafter: the first reports on seismo-ionospheric precursors; the second concerns the terminator method and its use for other EQs by research groups; the third insists on the physical

<footnote n="1">1 Space Research Institute, Austrian Academy of Sciences, 8042 Graz, Austria; hans.ulrich.eichelberger@oeaw.ac.at (H.U.E.); werner.magnes@oeaw.ac.at (W.M.); helmut.lammer@oeaw.ac.at (H.L.); wolfgang.voller@oeaw.ac.at (W.V.)</footnote>

<footnote n="2">2 Laboratoire Atmosphère, Milieux, Observations Spatiales-Centre National de la Recherche Scientifique, UVSQ Université Paris-Saclay, 78280 Guyancourt, France; patrick.galopeau@latmos.ipsl.fr</footnote>

<footnote n="3">3 Institute of Communication Networks and Satellite Communications, Graz University of Technology, 8010 Graz, Austria; s.sawas@zoho.com</footnote>

<footnote n="4">4 Institute for Applied Mathematics (IAC), National Research Council of Italy (CNR), 70126 Bari, Italy; g.nico@ba.iac.cnr.it</footnote>

<footnote n="5">5 Department of Physics, University of Bari, 70126 Bari, Italy; pf.biagi@gmail.com</footnote>

<footnote n="6">6 Department of Geodesy and Surveying, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece; mcont@topo.auth.gr</footnote>

<footnote>* Correspondence: mohammed.boudjada@oeaw.ac.at</footnote>
mechanisms at the origin of sub-ionospheric VLF anomalies associated with EQs; and the fourth summarizes the recent investigations into the Turkey–Syria seismic precursors.

## 1.1. Acoustic and Chemical Seismo-Ionospheric Precursor Channels

The interest in very low frequencies (VLFs, from 3 kHz to 30 kHz) is principally due to their propagation in the Earth's wave guide, defined by the lithosphere and the ionosphere, with a low signal attenuation over long distances up to 20,000 km [7]. Navigation and communication facilities established powerful VLF ground-based transmitter stations, mainly to connect and network with military submarines [3,8]. EQ anomalies linked to VLF sub-ionospheric propagation were first reported by Gokhberg et al. [9], Hayakawa et al. [10] and Molchanov and Hayakawa [11]. These ground-based outcomes allowed later the provision of basic physical processes which occur in the lithosphere, the atmosphere and the ionosphere, the so-called 'LAI coupling' [12]. Hence, micro-fracturing, liquid diffusions and pressure variations in the seismic preparation zone engender free electrons and electric charges at the Earth's surface. These lithospheric pre-seismic disturbances affect the conductivity of the atmosphere, where weak amplitude of the acoustic gravity wave (AGW) appears [13–15]. The decrease in atmospheric density leads to an intensification of the AGW amplitude when the altitude increases. Turbulence related to AGW amplification generates, in the upper part of the atmosphere, enhancements of gravity waves (GWs) and planetary waves (PWs). These pre-seismic ionospheric disturbances are linked to GW and PW effects [16,17]. These phenomena at the origin of pre-seismic VLF sub-ionospheric disturbances belong to the so-called acoustic channel [3].

It is important to emphasize the chemical channel, another physical process which takes place in the atmosphere, besides the AGW turbulence earlier reported. This concerns the radon ionization of the atmosphere, where the radon is considered as a radioactive element broadly exploited for geophysical studies. In the case of the Kobe EQ, a radon concentration enhancement was recorded two months before the seismic event occurrence [18], and the investigation of the radon daily residual values for Kobe EQ suggested a log-periodic oscillation model [19], where infrequent and occasional radon variations were linked to the local stresses instead of the Kobe main fault. The significant ionization of the atmospheric boundary layer furnishes punch power to the synergic development of coupling processes, where the radon is the first link in the chain, followed by the ionosphere, and the magnetosphere [20,21]. Recent publications insisted on the role of the radon physical process linked to seismic precursors for large earthquakes [22], short-term forecasting for earthquakes [23], and meteorological anomalies during EQ preparation [24].

## 1.2. Study of Earthquake Precursors Using the Terminator Method

The Kobe earthquake in Japan occurred on 17 January 1995 with a magnitude of Mw = 7.2. Anomalous electromagnetic occurrences were reported at frequency ranges from ELF to VHF [3]. The pioneer work of Hayakawa et al. [10] made evident particular spectral features linked to the amplitude and the phase of very low frequency (VLF) sub-ionospheric transmitter signals in the case of Kobe EQ. On this occasion, clear and sharp time shifts in sunrise and sunset terminators occurred a few days before the EQ. Such terminator time shifts are missing in the absence of an earthquake.

The basic concept of the terminator method involves a radio technique, where the VLF signal is emitted by a transmitter station, and detected after propagation by a VLF reception facility. The transmitter signal propagates in the Earth's waveguide before it attains the reception facility. The detected radio signal is modulated by the Earth's rotation, where the ionospheric E-layer and D-layer have crucial effects on the signal propagation, respectively, during the night and during the day [7]. At terminator occasions, the transmitter amplitude
and phase signals exhibit two minima, one at the sunrise terminator and the other one at the sunset terminator [3]. The terminator method is linked to the time shifts of these minima with regard to the terminator, i.e., sunrise or/and sunset, considered as a seismic precursor anomaly by Hayakawa et al. [10]. Later, several research groups showed the presence of such anomalies using the terminator method in the case of EQs which occurred in: Pakistan Mw = 7.2 [25], India Mw = 9.0 [26] and Mw = 6.0 [27], Nepal Mw = 7.3 [28], and Japan Mw = 9.0 [28,29]. These studies allowed the introduction of new parameters related to ionospheric D-layer behavior, i.e., formation and fading times at terminators of the D-layer [27], along with the use of numerical simulation based on an ionospheric model [30] to reproduce the variation in the VLF signal during its propagation along the ray path, by considering the effective reflection height ($h'$) and effective steepness ($\beta$) parameters [28]. Recently, Politis et al. [31] investigated six-years' (2014–2020) of VLF terminator time shifts observed during seismic events in Japan. They introduced the effective EQ magnitude ($M_{eff}$), defined as the total seismic energy affecting the midpoint of a given daily path. The applied statistical methodology [32] yielded the selection of terminator time anomalies, and the processing of a daily $M_{eff}$ and surface wave magnitude [33]. The investigation of the cross-correlation between effective EQ magnitudes and terminator time anomalies allowed the exhibition of maximum values prior to seismic occurrences, with a lag between –7 days and –1 day.

## 1.3. Physical Mechanisms at the Origin of Subionospheric VLF Disturbances Linked to EQs

Molchanov and Hayakawa [11] investigated about ten EQs with magnitude higher than Mw = 6.0. They reported similar terminator time shifts and suggested gravity waves (GWs) as physical processes to explain such seismic anomalies. Previously, Russian scientists [34,35] suggested GWs and acoustic gravity waves (AGWs) as the media connecting lithospheric seismic gravitational vibrations affecting primarily the atmosphere and, later, the Earth's ionosphere. All these disturbances occurred before the seismic event occurrences. These basic ideas regarding GWs and AGWs allow for further investigations of the lithospheric atmosphere-ionosphere coupling, as detailed in Pulinets and Boyarchuk [5]. Terminator time shift anomalies are the pointers and manifestations of the VLF wave propagation disturbances linked to LAI coupling impacts. Yoshida et al. [36] provided a terminator time generation mechanism to explain the ionospheric effect on VLF wave propagation for a moderately short path of about 2000 km. Sky and ground waves are considered in this mechanism [36] as the source of the terminator times resulting from the destructive interference between both waves. In this context, the authors applied a wave-hop method developed from ray theory [37] and derived the amplitude and phase variations at terminator times. They found that the change in VLF wave signal is linked to a decrease in reflection height by 3 km, resulting in the shift to earlier hours for the sunrise terminator and to later hours for the sunset terminator. More recently, Rapaport et al. [38] characterized seismic disturbances using the spectral method to study periodicities in VLF signal. They showed the presence of quasi-wave oscillations (i.e., periods of 4–10 min and 20–25 min) associated with atmospheric gravity waves. The authors emphasized the use of a synergetic approach [39] allowing the determination of entropy and information [40,41] for the characterization of seismic precursors.

## 1.4. Anomalous Seismic Precursor Reported for Turkey–Syria EQ of 6 February 2023

A series of damaging and injurious EQs struck the southern part of Turkey on 6 February 2023, close to the Turkish–Syrian border. It is considered as the greatest seismic event since 1939 in this region, with a magnitude of Mw = 7.8 [42]. Several studies reported earthquake precursors using space and ground-based observations. Contadakis et al. [43]
analyzed TEC values, recorded by the EUREF Network (http://www.epncb.oma.be; accessed on 12 February 2025), allowing a large coverage of the Mediterranean area, both before and after the mainshock. The authors found that the high-frequency limit of the ionospheric turbulence content progressively increases before the EQ. Similar results have been reported by Contadakis et al. [44] in the case of the Crete (Greece) EQ Mw = 6.5 of 12 October 2013. TEC derived from GNSS receivers and ionosondes have been investigated by Vesnin et al. [45], showing the existence of co-seismic ionospheric disturbances presenting significant asymmetry in propagation. Principally using space satellite observations, like CSES1 and Swarm, combined with Earth's atmospheric data derived from Google Earth Engine (GEE) and the Giovanni NASA data service, Akhoondzadeh and Marchetti [46] showed the presence of precursors at the end of January, and at the beginning of February 2023. Liu et al. [47] analysed ionospheric plasma parameters recorded onboard CSES and Swarm satellites, from 11 January to 14 February 2023. The authors reported irregularities in electron and $O^+$ densities, principally when satellite trajectories were above the epicenter regions. In both cases (i.e., CSES and Swarm satellites), anomalies appear a few days before the EQ occurrence, i.e., 30 January 2023. Similar results were also reported by Boudjada et al. [48], where only ground-based facilities were used, i.e., the Graz reception station (Austria) and the VLF transmitter (Turkey), in the vicinity of the EQ epicenter. Previous investigations confirmed the presence of common precursor periods despite numerous different reception systems from ground and space, and the different methods applied.

In this paper, we begin in the next Section to give more details of the multi-terminators method based on the main outcomes of Boudjada et al. [49] and Boudjada et al. [48]. This leads to emphasis in Section 3 on the new features derived from the amplitude and the phase of the TBB transmitter signal, as detected by the Graz VLF facility, and quantification of the occurrence probability of spectral anomaly features. In Section 4, the results and analysis of VLF sub-ionospheric spectral anomalies are discussed. The main outcomes and future perspectives are outlined in Section 5.

## 2. Multi-Terminators Method

In addition to the conventional paper of Hayakawa et al. [10], new attributes have been advanced in the context of the multi-terminator method applied to study EQs which occurred in the eastern-southern part of the area covered by the Graz facility (Austria), two in Croatia (22 March 2020 Mw = 5.4 and 29 December 2020 Mw = 6.5) and one in Turkey (6 February 2023 Mw = 7.8). This method has been introduced in [49], and further detailed in [48], as hereafter summarized. The originality of the multi-terminator method is its capability to provide indications and hints about the geographical locations of the precursor preparation zones. Such precursor areas are derived from anomalies in transmitter signals during their propagations in the seismically disturbed ionosphere. A close connection exists in this context between space ionosphere variability (i.e., due to seismic perturbations) and the ground detection of transmitter signals. As reported in Pulinets and Boyarchuk [5], the disturbed ionospheric area is comparable to the Dobrovolsky preparation zone.

### 2.1. *Case Study of 22 March 2020 Mw = 5.4 and 29 December 2020 Mw = 6.5 Croatia EQs*

In the first paper, Boudjada et al. [49] reported on EQs which occurred close to Zagreb (Croatia) at a distance of about 200 km from the Graz VLF facility (Austria). They used amplitude signals of four transmitters (i.e., ICV, ITS, TBB and RRO) as recorded at the Graz facility (15.43° E, 47.06° N, Austria) by the UltraMSK-1 reception system and by the INFREP network. Precursor anomalies are linked to the presence of time shifts at the Graz terminators. Three new features were investigated: (a) the time shift anomalies were compared to the local terminators and to the sunrise and sunset terminators of the transmitter
stations (see Figure 3 and Figure 4 in Boudjada et al. [49]); (b) a geometrical approach was initiated, in which VLF transmitters are supposed to isotropically and circularly emit radiation, permitting to a delineation of the so-called 'presumed pre-seismic areas' (see Figure 5 and Figure 6 in Boudjada et al. [49]); (c) the estimation of the anomaly correlation degree by combining two transmitter signals, leading to quantification of the maximum correlation proportion, and the corresponding lag and date (see Table 2, Boudjada et al. [49]). The multi-terminators method in the first paper allowed the conclusion that the presumed pre-seismic areas were confined to the Dobrovolsky preparation zone [50], and that the anomalies appear on 10 March 2020 and 18 December 2020, respectively, 12 days and 11 days before the EQs.

## 2.2. Case Study of 6 February 2023 Mw = 7.8 Turkey Syria EQ

The second paper concerns the EQ which occurred on 6 February 2023 with a magnitude of Mw = 7.8. The TBB transmitter signal was recorded by the UltraMSK-2 reception system [51] to analyze the ionospheric disturbances before the EQ. The VLF emitter (Bafa, Turkey; 27.31° E, 37.40° N) and the EQ epicenter (Gaziantep, Turkey; 37.17° N, 37.08° S) belong to the same region, i.e., the Anatolian region, and are separated by a distance of about 865 km. The multi-terminators method has also been applied in the second paper with new features: (a) further spectral features in the TBB signal, like maxima and jumps (see Figure 4, Boudjada et al. [48]) were studied; (b) each anomaly was linked to the Graz facility terminator or transmitter station terminator (see Figure 5 and Figure 6, and Table 4 in Boudjada et al. [48]); (c) the pre-seismic geographical coordinates were derived from minima and maxima terminators (see Table 5 in Boudjada et al. [48]). These newly introduced spectral features showed the relationship between precursor anomalies and seismic geographical regions.

# 3. Occurrence Probability of Spectral Anomaly Features

We analyze in this Section the amplitude and phase variations of TBB transmitters as detected by the Graz facility. The distances and the sensitive regions, as derived from the Dobrovolsky preparation area and Fresnel zone, are given in Section 3.1. An overview of the TBB amplitude and phase variations is given in Section 3.2, followed by analysis of minima and maxima features in Section 3.3, and jumps and inflexions, i.e., new anomalies, in Section 3.4. Both the methods of Hayakawa et al. [10] and Boudjada et al. [48] are considered in Section 3.5. The occurrence probability of the observed spectral anomalies is detailed in Section 3.6, followed by the estimation of the spectral anomaly time shifts in Section 3.7. A period of 17 days is analysed in this work, 14 days from 23 January to 5 February, the day of the EQ on 6 February, and 2 days from 7 to 8 February 2023.

## 3.1. Geographical Location of Turkey–Syria Earthquake

Figure 1 displays the geographical locations of the Graz facility (GRZ, Austria), TBB transmitter station (TBB, Turkey), and earthquake epicenter at Gaziantep (EQ, Turkey). It can be observed that the GRZ facility is localized in the northern part at a latitude of about 47° N, and the TBB station and EQ epicenter in the southern part at a latitude in the order of 37° N. The distances GRZ–EQ, GRZ–TBB and TBB–EQ are, respectively, in the order of 2080 km, 1450 km, and 860 km. We define a pre-seismic sensitive region using two well-known relationships. The Dobrovolsky preparation zone is estimated from $ρ = 10^{0.43M}$, where M is the EQ magnitude and ρ is the Dobrovolsky radius shown in Figure 1 with red line arrows. We find ρ equal to 2260 km when considering that Mw = 7.8 in the investigated seismic event. The second relationship is the elliptic Fresnel zone, where the minor semi-axis *b* is derived from $b \approx [\lambda. D/2]^{0.5}$ *D* is the distance between the TBB
emitter and GRZ receiver equal to 1449 km, and $\lambda$ is the transmitter wavelength of 77 km. In this relationship, we neglect the ratio $(h/D)^2$, which is found to be much smaller than 1, where $h$ is the average reflection height equal to 80 km. The orange elliptical curve in Figure 1 defines the fifth Fresnel zone, with $b$ equal to about 1180 km. The pre-seismic sensitive region is found to be further extended when one takes into consideration the Dobrovolsky preparation zone and limited with regard to the Fresnel elliptic area. Despite the difference between both areas, the ray paths of the TBB transmitter signal are mainly inside the sensitive region, enabling high detection probabilities of precursor anomalies.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Geographical locations of EQ epicenter (i.e., EQ), TBB transmitter (i.e., TBB) and Graz facility (i.e., GRZ). All distances between these three geographical points are indicated by black dashed lines. The pre-seismic sensitive regions are derived from the Dobrovolsky preparation zone and the Fresnel area, defined by red arrow lines and orange elliptic curve, respectively. The Dobrovolsky radius, $\rho$, is calculated from the relationship $\rho = 10^{0.43M}$, where M is the EQ magnitude, and the semi-axis $b$ of the elliptic Fresnel is estimated from $b \approx [\lambda \cdot D/2]^{0.5}$, with D the distance TBB-GRZ stations, equal to 1449 km, and $\lambda$ the transmitter wavelength, equal to 77 km.</figcaption></figure>

## 3.2. Overview of Amplitude and Phase Variations of TBB Transmitter Signal

Hereafter, we give an overview of the TBB transmitter signal variation. It is important to report that a new reception system (i.e., UltraMSK-2) was installed at the end of 2021 at the Graz facility, tested in the first semester of 2022, and has been in full operation since September 2022. More details about UltraMSK-2 software can be found on https://www.ultramsk.com (accessed on 18 June 2025). The absence of intense seismic events, since September 2022, did not allow us to consider other EQs. The occurrence of the Turkey-Syria EQ provides us with the first good opportunity to use UltraMSK-2 observations. Only the TBB transmitter signal has been investigated in this work because its ray paths occurred in the Dobrovolsky preparation zone [50] and its location is not far from the EQ epicenter location. Other transmitters, particularly the ICV and ITS VLF transmitters, have not been included, because their ray paths were outside the Dobrovolsky preparation zone, and
their phase variations exhibit several rotations per day, contrary to the TBB phase, which is almost constant during the day and nighttime intervals.

Figure 2 displays the daily recorded TBB signal fluctuations in the period from 29 January (29 DOY (day of the year 2023)) to 7 February (39 DOY) 2023. The time of the earthquake is marked by red vertical dashed lines for the amplitude and the phase signals. The Graz facility and TBB station terminators are indicated, respectively, by magenta and green vertical dashed lines. The two upper panels (i.e., with 10 sub-panels) of Figure 2 illustrate the amplitude variation, where the beginning and the end of the night fluctuations are marked by the sunset and sunrise time terminators, and vice versa for the day variations. At terminator times, the level of the amplitude is, on average, in the order of –75 dB, very close to the background noise. During the night and day periods, the amplitude signals are, respectively, intense at about –45 dB and decrease to less than –55 dB. The two lower panels of Figure 2 (i.e., with 10 sub-panels) display the phase variations of the TBB transmitter signal. As in the case of the amplitude fluctuations, magenta, green and red vertical dashed lines indicate, respectively, the TBB terminators, Graz terminator and time of EQ occurrence. Sunrise and sunset terminators separate day and night observations with the particular presence of abrupt variations of phase angle, as detailed in Section 3.4. During the investigated period, from 23 January 2023 to 8 February 2023, solar and geomagnetic activities were generally low. The value of the mean 3 h-Kp-index is principally found to be lower than 5, corresponding to low geomagnetic activity. The solar effect is mainly associated with solar flares, where the amplitude and phase are intensified during several minutes, mostly during the day, and not during night observations. Only on two occasions were solar flares recorded on four days, at the end of January (i.e., 30 DOY and 31 DOY), and close to the EQ occurrence (i.e., 36 DOY and 37 DOY). It can be seen in Figure 2, for the previously cited days, the enhancements of the amplitude and the phase at around mid-day.

<figure><img src="image_4.png" /></figure>

**Figure 2.** Daily amplitude (first and second panels) and phase (third and fourth panels) variations of the TBB transmitter signal from 29 January (29 DOY) to 7 February (38 DOY) 2023. The horizontal axis indicates the time (expressed in UT), and the vertical axis designates the amplitude (expressed in dB) and phase (expressed in degree). Vertical magenta and green dashed lines specify, respectively, the TBB station and Graz facility sunrise and sunset terminators. The occurrence time of the earthquake on 6 February 2023 is shown by a red vertical dashed line.
### 3.3. Minima and Maxima Anomalies at Sunrise and Sunset Terminators

Figure 3 shows the variation of extrema by day of the year (DOY) 2023, as investigated by Boudjada et al. [48]. The colored circles indicate minima (red) and maxima (blue) for the amplitude, and minima (black) and maxima (green) for the phase. The day of the earthquake (i.e., 6 February 2023—37 DOY) is indicated by a red vertical dashed line. The terminator variations of Graz VLF facility (47.03° N, 15.46° E, Graz, Austria), TBB transmitter station (37.40° N, 27.31° E, Baffa, Turkey) and EQ epicenter (37.17° N, 37.08° E, Gaziantep, Turkey) are signaled, respectively, by green, magenta and black lines. From the 284 extrema recorded in 17 days (i.e., from 23 January to 8 February 2023), 167 extrema were observed at sunrise (upper panel of Figure 3) and 117 at sunset (lower panel of Figure 3).

<figure><img src="image_5.png" /></figure>

Figure 3. Minima and maxima (expressed in hours) derived from the TBB transmitter signal for days of the year 2023, as reported by Boudjada et al. [48]. The horizontal and vertical axes correspond,
respectively, to the DOY 2023, and the terminator times expressed in hours. The extrema observed at sunrise and sunset terminators are displayed, respectively, in the upper and the lower panels. The minima and maxima of the amplitude are designated, respectively, by red and blue colored circles. Minima and maxima of the phase are indicated by black and green colored circles. The earthquake day (37 DOY) is designated by the vertical red dashed line. The investigated period starts on 23 January 2023 (23 DOY) and ends on 8 February 2023 (39 DOY). The green, magenta, and black colored lines indicate the sunrise or sunset terminator variations, respectively, at the Graz VLF facility (47.03° N, 15.46° E), TBB transmitter station (37.40° N, 27.31° E), and the earthquake epicenter (37.17° N, 37.08° E).

Table 1 lists the number of minima and maxima recorded at sunrise and sunset terminators as derived from the amplitude and the phase of the TBB transmitter signal. These extrema are displayed in Figure 5 and Figure 6 in Boudjada et al. [48], respectively, for amplitude and phase.

<table><caption>**Table 1.** List of extrema recorded from 23 January to 8 February 2023 as reported by Boudjada et al. [48].</caption><thead><tr><th rowspan="2"></th><th colspan="2">Sunrise</th><th colspan="2">Sunset</th></tr><tr><th>Minima</th><th>Maxima</th><th>Minima</th><th>Maxima</th></tr></thead><tbody><tr><td>Amplitude</td><td>44</td><td>34</td><td>29</td><td>21</td></tr><tr><td>Phase</td><td>50</td><td>39</td><td>27</td><td>40</td></tr></tbody></table>

### 3.4. Jumps and Inflexions as Spectral Features

In Boudjada et al. [48] only extrema (minima and maxima) have been exploited to study the terminator anomalies as indicated in Section 3.3 and displayed in Figure 3. Hereafter, we show zooms of the anomalies recorded on 1 February 2023 around sunrise and sunset terminators to clarify how spectral features are derived from daily amplitude and phase variations. This leads us to introduce two new features, both observed around terminators: the first is of the jump spectral type, which occurs in the TBB phase signal where the phase springs from -180° to +180° in about one minute, and the second is the spectral inflexion, which occurs in the TBB amplitude signal, where one can note a brutal change, in a short time, of the amplitude slope signs, i.e., from positive to negative and vice versa.

Figures 4 and 5 display amplitude and phase variations, respectively, as recorded on 1 February 2023 (32 DOY). UltraMSK-2 reception leads to the recording of three measurements per minute, one each 20 s, for the amplitude and the phase transmitter signal. The observations in Figures 4 and 5 are derived from the 8640 daily measurements. All spectral features occur around TBB and Graz sunrise and sunset terminators. Minima and maxima are indicated by red and blue colored circles and jumps and inflexions, both new anomalies, are designed by orange-colored circles. The jumps appear only in the TBB phase signal, as shown in the zoomed part of Figure 5, and the inflexions can be seen in the zoomed part of Figure 4.

Table 2 lists the inflexion spectral features, where the sunrise or sunset terminator (first column), the observation date (second column), the time in UT (third column), the amplitude in dB (fourth column), the phase in degree (fifth column), and finally the drift rate of the inflexion (sixth column) are given. For the investigated period from 23 January to 8 February 2023, the daily inflexion type has been observed 13 times at sunrise and 18 at sunset. Such spectral types are absent on two occasions, i.e., on 26 January and 28 January.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Amplitude of TBB transmitter signal recorded on 1 February 2023, five days before EQ occurrence. Horizontal and vertical axes indicate, respectively, the observation time (expressed in DOY) and the amplitude level (expressed in dB). Sunrise and sunset terminators for TBB transmitter station and Graz facility are designated, respectively, by magenta and green vertical dashed lines. Zoomed parts are shown for sunrise (left small-panel) and sunset (right small-panel), where minima, maxima and inflections, i.e., spectral anomalies, are specified by red, blue and orange circles, respectively.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

**Figure 5.** The first panel shows the phase of the TBB transmitter signal, recorded on 1 February 2023, where horizontal and vertical axes indicate, respectively, the observation time (expressed in DOY) and the phase (expressed in degree). As in Figure 4, magenta and green vertical dashed lines correspond to sunrise and sunset terminators for TBB transmitter station and Graz facility. Zoomed parts are shown for sunrise (left small-panel) and sunset (right small-panel), where minima, maxima and jumps, i.e., spectral anomalies, are specified by red, blue and orange-colored circles, respectively.
<table><caption>**Table 2.** List of inflexion spectral type recorded in the period from 23 January 2023 to 08 February 2023.</caption><thead><tr><th>Terminator</th><th>Date</th><th>Time (UT)</th><th>Amplitude (dB)</th><th>Phase (Degree)</th><th>Drift Rate (dB/min)</th></tr></thead><tbody><tr><td>Sunset</td><td>23 January</td><td>16:14</td><td>-57.05</td><td>+157.66</td><td>+0.240</td></tr><tr><td>Sunset</td><td>24 January</td><td>16:34</td><td>-53.68</td><td>+148.75</td><td>+0.100</td></tr><tr><td>Sunrise</td><td>25 January</td><td>04:42</td><td>-51.05</td><td>+137.09</td><td>-0.180</td></tr><tr><td>Sunset</td><td>25 January</td><td>16:15</td><td>-57.47</td><td>+157.47</td><td>+0.345</td></tr><tr><td>Sunset</td><td>26 January</td><td>16:18</td><td>-53.54</td><td>+170.54</td><td>+0.360</td></tr><tr><td>Sunset</td><td>28 January</td><td>14:41</td><td>-63.03</td><td>+095.44</td><td>-0.510</td></tr><tr><td>Sunset</td><td>29 January</td><td>14:18</td><td>-62.00</td><td>+028.89</td><td>-0.420</td></tr><tr><td>Sunset</td><td>29 January</td><td>15:14</td><td>-62.51</td><td>-034.17</td><td>+0.405</td></tr><tr><td>Sunrise</td><td>30 January</td><td>06:48</td><td>-62.51</td><td>+088.32</td><td>-0.150</td></tr><tr><td>Sunset</td><td>30 January</td><td>16:17</td><td>-52.01</td><td>-109.93</td><td>+0.120</td></tr><tr><td>Sunrise</td><td>1 February</td><td>06:40</td><td>-60.50</td><td>+093.21</td><td>-0.120</td></tr><tr><td>Sunset</td><td>1 February</td><td>14:32</td><td>-64.78</td><td>+070.22</td><td>-0.480</td></tr><tr><td>Sunset</td><td>1 February</td><td>16:22</td><td>-52.01</td><td>-126.67</td><td>+0.225</td></tr><tr><td>Sunrise</td><td>2 February</td><td>04:49</td><td>-48.01</td><td>-143.11</td><td>-0.150</td></tr><tr><td>Sunset</td><td>2 February</td><td>14:09</td><td>-62.04</td><td>+029.47</td><td>+0.120</td></tr><tr><td>Sunset</td><td>2 February</td><td>16:21</td><td>-54.03</td><td>-126.33</td><td>+0.360</td></tr><tr><td>Sunrise</td><td>3 February</td><td>04:40</td><td>-46.51</td><td>-149.41</td><td>-0.240</td></tr><tr><td>Sunset</td><td>3 February</td><td>16:29</td><td>-53.50</td><td>-126.56</td><td>+0.330</td></tr><tr><td>Sunrise</td><td>4 February</td><td>04:59</td><td>-50.53</td><td>-127.26</td><td>-0.450</td></tr><tr><td>Sunrise</td><td>4 February</td><td>07:11</td><td>-66.52</td><td>+047.48</td><td>+0.630</td></tr><tr><td>Sunrise</td><td>5 February</td><td>04:54</td><td>-48.97</td><td>-130.80</td><td>-0.330</td></tr><tr><td>Sunrise</td><td>5 February</td><td>07:18</td><td>-67.25</td><td>+040.03</td><td>-1.860</td></tr><tr><td>Sunset</td><td>5 February</td><td>14:31</td><td>-65.07</td><td>+028.18</td><td>+0.810</td></tr><tr><td>Sunset</td><td>6 February</td><td>15:18</td><td>-66.26</td><td>-059.87</td><td>+1.230</td></tr><tr><td>Sunrise</td><td>7 February</td><td>04:58</td><td>-51.40</td><td>-116.96</td><td>-0.585</td></tr><tr><td>Sunrise</td><td>7 February</td><td>06:54</td><td>-65.94</td><td>+079.02</td><td>+2.460</td></tr><tr><td>Sunset</td><td>7 February</td><td>16:28</td><td>-55.51</td><td>-121.42</td><td>0.210</td></tr><tr><td>Sunrise</td><td>8 February</td><td>04:52</td><td>-48.96</td><td>-124.88</td><td>-0.180</td></tr><tr><td>Sunrise</td><td>8 February</td><td>07:01</td><td>-67.02</td><td>+056.95</td><td>+0.630</td></tr><tr><td>Sunset</td><td>8 February</td><td>14:00</td><td>-58.51</td><td>+036.26</td><td>-0.210</td></tr><tr><td>Sunset</td><td>8 February</td><td>16:31</td><td>-50.50</td><td>-114.36</td><td>+0.270</td></tr></tbody></table>

Table 3, similarly to Table 2, indicates the jump spectral features, where we find that five jumps occur at sunset and twelve at sunrise. As in the case of inflexion type, jumps are missing on two days, 27 and 31 January 2023. It is important to note that the jump type shows mainly a positive slope, i.e., when the time increases, the phase springs from negative to positive values, as shown in the sixth column of Table 3. Only at three occasions have the slopes been negative: (a) on 29 January at sunrise, with a decline phase from +179.54° to -179.37°; (b) on 31 January at sunset, with a decrease from +176.44° to -177.57°; and (c) one day before the EQ occurrence at sunrise, with a drop from +179.81° to -178.79°.

It is interesting to note that the reported inflexion slopes (see sixth column in Table 2) exhibit positive and negative maxima at three occasions: (a) 28 January at sunset with a value of -0.510 dB/min; (b) 1 February at sunrise equal to -0.480 dB/min; and (c) 4 February at sunrise with a value of +0.630 dB/min. One day before (5 February) and after (7 February) the EQ, the slopes are equal to -1.86 dB/min and +2.46 dB/min, respectively. Similarly, the drift rates of the jumps (see sixth column in Table 3) also displays positive and negative maxima: (a) 27 January at sunset with an amount of +1078.44 degree/min; (b) 30 January at sunrise with a value of +1071.75 degree/min' and (c) 2 February at sunrise with a drift of about +1071.39 degree/min. As for the inflexion slopes, on 5 and 7 February, values equal to -1075.80 degree/min and +1073.04 degree/min were recorded.
<table><caption>**Table 3.** List of jump spectral type recorded in the period from 23 January 2023 to 8 February 2023.</caption><thead><tr><th>Terminator</th><th>Date</th><th>Time (UT)</th><th>Amplitude (dB)</th><th>Phase (Degree)</th><th>Drift Rate (Degree/min)</th></tr></thead><tbody><tr><td rowspan="2">Sunset</td><td>23 January</td><td>15:52:40</td><td>-67.05</td><td>-179.96</td><td></td></tr><tr><td>23 January</td><td>15:53:00</td><td>-66.63</td><td>+179.76</td><td>+0539.58</td></tr><tr><td rowspan="2">Sunset</td><td>24 January</td><td>16:01:40</td><td>-62.29</td><td>-178.73</td><td></td></tr><tr><td>24 January</td><td>16:02:00</td><td>-62.05</td><td>+179.80</td><td>+0537.80</td></tr><tr><td rowspan="2">Sunset</td><td>25 January</td><td>16:00:40</td><td>-64.73</td><td>-178.41</td><td></td></tr><tr><td>25 January</td><td>16:01:00</td><td>-54.30</td><td>+179.90</td><td>+0537.51</td></tr><tr><td rowspan="2">Sunset</td><td>27 January</td><td>16:04:00</td><td>-54.25</td><td>-179.82</td><td></td></tr><tr><td>27 January</td><td>16:04:20</td><td>-66.63</td><td>+179.66</td><td>+1078.44</td></tr><tr><td rowspan="2">Sunrise</td><td>29 January</td><td>06:02:40</td><td>-71.35</td><td>+179.54</td><td></td></tr><tr><td>29 January</td><td>06:03:00</td><td>-71.12</td><td>-179.37</td><td>-0538.37</td></tr><tr><td rowspan="2">Sunrise</td><td>30 January</td><td>05:39:20</td><td>-67.80</td><td>-179.98</td><td></td></tr><tr><td>30 January</td><td>05:39:40</td><td>-68.08</td><td>+177.27</td><td>+1071.75</td></tr><tr><td rowspan="2">Sunrise</td><td>31 January</td><td>05:33:20</td><td>-71.60</td><td>-174.09</td><td></td></tr><tr><td>31 January</td><td>05:33:40</td><td>-71.77</td><td>+179.18</td><td>+1059.81</td></tr><tr><td rowspan="2">Sunset</td><td>31 January</td><td>15:51:40</td><td>-71.52</td><td>+176.44</td><td></td></tr><tr><td>31 January</td><td>15:52:00</td><td>-71.27</td><td>-177.57</td><td>-0531.02</td></tr><tr><td rowspan="2">Sunrise</td><td>1 February</td><td>05:26:40</td><td>-71.49</td><td>-178.60</td><td></td></tr><tr><td>1 February</td><td>05:27:00</td><td>-72.41</td><td>+173.16</td><td>+0527.64</td></tr><tr><td rowspan="2">Sunrise</td><td>2 February</td><td>05:22:00</td><td>-63.00</td><td>-179.98</td><td></td></tr><tr><td>2 February</td><td>05:22:20</td><td>-63.32</td><td>+177.15</td><td>+1071.39</td></tr><tr><td rowspan="2">Sunrise</td><td>3 February</td><td>05:17:20</td><td>-61.27</td><td>-179.68</td><td></td></tr><tr><td>3 February</td><td>05:17:40</td><td>-61.46</td><td>+177.27</td><td>+1070.85</td></tr><tr><td rowspan="2">Sunrise</td><td>4 February</td><td>05:21:40</td><td>-72.21</td><td>-179.56</td><td></td></tr><tr><td>4 February</td><td>05:22:00</td><td>-73.23</td><td>+176.77</td><td>+0534.50</td></tr><tr><td rowspan="2">Sunrise</td><td>5 February</td><td>05:20:40</td><td>-70.08</td><td>-177.59</td><td></td></tr><tr><td>5 February</td><td>05:21:00</td><td>-70.72</td><td>+179.45</td><td>+0535.56</td></tr><tr><td rowspan="2">Sunrise</td><td>5 February</td><td>06:08:00</td><td>-58.33</td><td>+179.81</td><td></td></tr><tr><td>5 February</td><td>06:08:20</td><td>-58.45</td><td>-178.79</td><td>-1075.80</td></tr><tr><td rowspan="2">Sunrise</td><td>6 February</td><td>05:23:40</td><td>-76.70</td><td>-175.16</td><td></td></tr><tr><td>6 February</td><td>05:24:00</td><td>-77.53</td><td>+167.10</td><td>+0513.39</td></tr><tr><td rowspan="2">Sunrise</td><td>7 February</td><td>05:21:00</td><td>-65.21</td><td>-179.00</td><td></td></tr><tr><td>7 February</td><td>05:21:20</td><td>-65.51</td><td>+178.68</td><td>+1073.04</td></tr><tr><td rowspan="2">Sunrise</td><td>8 February</td><td>05:21:20</td><td>-63.89</td><td>-175.93</td><td></td></tr><tr><td>8 February</td><td>05:21:40</td><td>-63.74</td><td>+179.90</td><td>+1067.49</td></tr></tbody></table>

### 3.5. Terminator Method (TM) Versus Multi-Terminator Method (MTM)

The fundamental paper of Hayakawa et al. [10] provides the method to identify the presence of precursor anomalies in the sub-ionospheric transmitter signal. It is important to note that, for the TM method, the terminator time shifts only concern minima around the local reception station. This method has been improved and advanced by Boudjada et al. [48] and Boudjada et al. [49] by including other spectral anomalies at the local reception facility, at the emitter transmitter, and close to the EQ epicenter.

The difference between both methods can be clearly established by considering the spectral anomalies in Figures 4 and 5 for the TBB signal observed on 1 February 2023. In Figure 4, we find that two minima occurred at sunrise and one minimum at sunset, a total of three minima, when applying the TM method and, from Figure 5, a total of six minima (i.e., three at sunrise and three at sunset). With the MTM method, one needs to add to the number of minima the spectral anomalies, two maxima and three inflections (see lower panels of Figure 4) and four maxima and two jumps (see lower panels of Figure 5). Therefore, 9 minima spectral anomalies are derived from the TM method, and 19 anomalies (i.e., 9 minima, 6 maxima, 2 inflections and 2 jumps) from the MTM method.
In this study, we have analyzed 138,240 spectral observations from 23 January to 8 February 2023, where only 349 anomalies have been found, representing 0.25% of the total measurements. The percentage of spectral anomalies are 43% minima, 38.4% maxima, 9.7% jumps, and 8.9% inflections when considering the 349 anomalies.

### 3.6. Occurrence Probability of Spectral Features Observed at Sunrise and Sunset Terminators

In this sub-section, we intend to compute the anomaly time occurrence times by separating the sunrise and sunset spectral anomalies. A specific time interval, called hereafter bin size, is used to find the anomaly occurrence probability every 12 min in the interval 04.40–07.60 UT for sunrise, and 13.20–16.80 UT for sunset. The bin size of 12 min allows study of the anomaly occurrence 5 times per hour (5 × 12 min = 60 min) corresponding to a ‘longitude bin size’ of about 3° (i.e., for 15°/5). For this, we use the following histogram function $H_v$ to compute the occurrence probability $P$ [52]:

$$H_v = \sum_{i=0}^{n-1} P(F_i, v), \qquad v = 0, 1, 2, \dots, \left[ \frac{Max - Min}{Binsize} \right] \quad (1)$$

where $F_i$ is the value of the element of $i$ index with $0 \le i < n$. The number of occurrences is given by $v$, as indicated in Equation (1), which depends on the maximum (Max), the minimum (Min) and the histogram bin size (Binsize). The occurrence probability $P$ is given by the following:

$$P(F_i, v) = \begin{cases} 1 & v \le (F_i - Min)/Binsize < v + 1 \\ 0 & \text{Otherwise} \end{cases} \quad (2)$$

Figure 6 displays the probability of anomaly occurrences derived for the spectral features observed at sunrise and sunset terminators. A total of 349 anomalies have been examined corresponding to 150 minima, 134 maxima, 31 inflections and 34 jumps. The Max–Min interval (expressed in UT hours) is, for the sunrise terminators, [4.38–7.80], corresponding to [04:23–07:48] (left panel in Figure 6), and, for the sunset terminators [13.21–16.98], corresponding to [13:12–16:59] (right panel in Figure 6). The bin size has been fixed to 12 min (0.20 in hours) for both histograms.

The three colored rectangles shown in Figure 6 correspond to the beginning and the end of the terminators at the Graz VLF facility (blue dashed rectangle), TBB transmitter station (green dashed rectangle) and EQ epicenter (red dashed rectangle). For the investigated period from 23 January to 8 February 2023, the terminator time intervals are found as follows:

* Graz at sunrise [06.53–06.18] or [06:32–06:11] UT, and at sunset [15.78–16.21] or [15:47–16:13] UT,

* TBB transmitter station at sunrise [05.31–05.08] or [05:19–05:05] UT, and at sunset [15.40–15.71] or [15:24–15:43] UT,

* EQ epicenter at sunrise [04.66–04.43] or [04:40–04:26] UT, and at sunset [14.76–15.08] or [14:46–15:05] UT.

The probability exhibits three occurrence maxima for the sunrise terminators (left panel in Figure 6), the first close to the EQ epicenter terminator between [04:45 and 04:57] UT with a value of 24, the second closer to the TBB transmitter terminator [05:20 and 05:31] UT with a value of 28, and the third overlapping the Graz facility terminator [06:17 and 06:28] UT with a value of 15. The second or the third occurrence peak can be assumed to be linked, respectively, to the transmitter station and the Graz facility terminators. However, the first peak with a value of 24 appears earlier when compared to the TBB terminators
(green rectangle in the left panel of Figure 6) and seems to be generated in the preparation
EQ zone in the eastern part of the area covered by the transmitter station.

<figure><img src="image_8.png" /></figure>

**Figure 6.** Occurrence probabilities of spectral anomalies recorded at sunrise (left panel) and sunset (right panel) as derived from TBB signal variations. The horizontal axis indicates the terminators (expressed in hours) and the vertical axis the number of occurrences per bin size. The red, green and blue rectangles designate the beginning and the end of the terminators, respectively, for EQ epicenter, TBB transmitter station and Graz VLF facility.

The occurrence in the case of the sunset terminators (right panel of Figure 6) displays
peaks, but with maximum values less than 17. Similarly, to sunrise occurrence, sunset
also has three occurrence maxima: (a) between [14:12 and 14:24] UT with a value of 6;
(b) half-overlapping the EQ epicenter terminators [14:48–15:00] UT with a value of 14; and
(c) partially covering the TBB transmitter and Graz terminators [15:47–15:59] UT with a
value of 16. It is interesting to note that the first peak, with a value of 14, is totally related to
the EQ preparation zone which seems to be located in the eastern part of the EQ epicenter.
Despite the small value of the first peak (i.e., equal to 5), the corresponding terminators
are situated in the eastern part of the EQ epicenter. This gives a first idea about the size of
this preparation zone, found in the eastern part of the area covered by the TBB station, and
which extended from 14:12 UT to 15:15 UT, i.e., to about one hour, corresponding to 16° in
geographical longitude.

The occurrence of the spectral feature anomalies, as shown in Figure 6, leads us to
confirm the presence of peaks which can be assumed as precursors, since they arrive
approximately at the EQ epicenter terminators. The sunset spectral features show evident
precursors, one minor at round 14:18 and another major at around 14:54. For the sunrise
occurrence, it is not simple to recognize and distinguish the precursors, since the peak,
around 04:51, is found between the TBB station and EQ epicenter sunrise terminators.
### 3.7. Estimation of Spectral Anomaly Time Shifts

It is clear from Figure 6 that spectral anomalies derived from the multi-terminators method lead us to reveal time intervals linked to the EQ preparation zone. This can be considered as a first stage towards a better knowledge of the efficiency of this method. A second stage is proposed in this sub-section and which consists in estimating the spectral anomaly time shifts, defined as the time difference (shown hereafter as $\Delta$) between the Graz sunrise/sunset terminators and the occurrence time of the spectral anomalies. This introduced quantity leads to daily quantification of when the spectral anomalies appear when compared to the daily Graz terminator. Hence, if $\Delta$ time difference is positive or negative, then the spectral features are, respectively, localized in the eastern or western part of the Graz facility. This quantity $\Delta$ is calculated for seven intervals, defined as follows:

* Group 1 are anomalies which occur after 06:32 UT at sunrise and after 16:13 at sunset.

* Group 2 are close to the Graz terminators, i.e., [06:11–06:32] UT at sunrise and [15:47–16:13] UT at sunset.

* Group 3 are between Graz and TBB terminators, i.e., [05:19–06:11] UT at sunrise and [15:43–15:47] UT at sunset.

* Group 4 are nearby TBB terminators, i.e., [05:05–05:19] UT at sunrise and [15:24–15:43] UT at sunset.

* Group 5 are between TBB and EQ terminators, i.e., [04:40–05:05] UT at sunrise and [15:05–15:24] UT at sunset.

* Group 6 are close to EQ epicenter terminators, i.e., [04:26–04:40] UT at sunrise and [14:46–15:05] UT at sunset.

* Group 7 are before 04:26 UT at sunrise and before 14:46 UT at sunset.

Figures 7 and 8 display $\Delta$ time difference (expressed in minutes) versus the day of the year (DOY) from 23 January to 8 February 2023. The amplitude spectral anomalies for sunrise and sunset are shown, respectively, in the upper and lower panels of Figure 7. Similarly, the sunrise and sunset phase spectral anomalies are illustrated in the upper and lower panels of Figure 8.

The red and blue colored symbols indicate in both figures the minima and maxima spectral anomalies, respectively. The black colored symbols design the inflexions and jumps spectral anomalies, respectively, in Figures 7 and 8. It can be seen that $\Delta$ time difference for the Graz facility is equal to zero (shown by green horizontal dashed line), for the TBB transmitter station $\Delta$ is close to +50 min (designed by magenta horizontal dashed line), and for the EQ epicenter $\Delta$ is about +85 min (indicated by a red horizontal dashed line).

#### 3.7.1. Amplitude Sunrise and Sunset Spectral Anomaly Time Shifts

In the upper panel of Figure 7, one can see minima anomalies (i.e., red colored symbols) around –50 min, followed by maxima (i.e., blue colored symbols) at about +25 min, and again minima around +50 min, maxima about +70 min, followed by minima around +90 min. These minima and maxima are directly linked to the TBB signal propagation in the waveguide due to the combined effect of the sky and ground waves.

The quasi-parallel lines exhibit significant irregularities, between –50 min and –20 min, starting from 27 January to 31. This is a clear precursor signature which affects sunrise terminators at Graz facility. In the upper panel of Figure 7 is indicated, in a blue rectangle, the first main precursor MP1. Simultaneously, there is a discontinuity in the maxima line on 27 January (i.e., absence of one maximum) at about $\Delta$~+50 min close to the TBB transmitter station, and also in the minima line on 28 and 29 January when $\Delta$ is about +25 min. However, maximum anomalies display a regular variation at the EQ epicenter terminators at around +90 min. Lower than $\Delta$~+10 min, one finds that three groups, i.e., group 1 ('+' symbol), group 2 ('*' symbol), and inflexion anomalies ('●' symbol),
have similar fluctuations, with a minimum on 25 January, followed by a maximum on 29 January, and again a minimum on 2 February 2023. In the Δ interval between [+10 min, +60 min], group 4 (‘◇’ symbol) exhibits nearly a constant variation. Higher than Δ~+70 min, similar fluctuations concern three groups, i.e., group 5 (‘□’ symbol), group 6 (‘x’ symbol), and group 7 (‘♦’ symbol). A prior precursor (PP) takes place between 23 January and 26 January, where similar variations are observed at Δ~−60 min and Δ~+120 min. This is designated by blue rectangles in the upper panel of Figure 7.

<figure><img src="image_9.png" /></figure>

**Figure 7.** Variation of Δ time difference (expressed in minutes) versus the day of the year (DOY) from 23 January to 9 February 2023 for the amplitude sunrise (upper panel) and for the amplitude, sunset (lower panel) anomalies. The green, magenta and red horizontal dashed lines design, respectively, the time difference for Graz facility (Δ = 0 min), for TBB transmitter station (Δ is close to +50 min), and for EQ epicenter (Δ is about +85 min). The red, blue and black colored symbols correspond to minima, maxima and inflexions spectral features. Δ has been estimated for seven intervals as defined in Section 3.7: ‘plus’ or ‘+’ symbol for group 1, ‘asterisk’ or ‘*’ symbol for group 2, ‘empty diamond’ or ‘◇’ symbol for group 3, ‘triangle’ or ‘Δ’ symbol for group 4, ‘square’ or ‘□’ symbol for group 5, ‘x’ symbol for group 6, and ‘fully diamond’ or ‘♦’ symbol for group 7. The inflexion anomaly is indicated by a ‘circle’ or ‘●’ symbol. The prior precursor (PP), the first and the second main precursors (MP1 and MP2) are displayed in blue rectangles.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** As in Figure 7, for the phase sunrise (**upper panel**) and for the phase sunset (**lower panel**) anomalies, the red, blue and black colored symbols correspond to minima, maxima and jumps spectral features. The green, magenta and red horizontal dashed lines design, respectively, the time difference for Graz facility (Δ = 0 min), for TBB transmitter station (Δ is close to +50 min), and for EQ epicenter (Δ is about +85 min). Δ has been estimated for seven intervals as defined in Section 3.7: ‘plus’ or ‘+’ symbol for group 1, ‘asterisk’ or ‘*’ symbol for group 2, ‘empty diamond’ or ‘◊’ symbol for group 3, ‘triangle’ or ‘Δ’ symbol for group 4, ‘square’ or ‘□’ symbol for group 5, ‘x’ symbol for group 6, and ‘fully diamond’ or ‘◆’ symbol for group 7. The jump anomaly is indicated by a ‘circle’ or ‘●’ symbol. The prior precursor (PP), the first and the second main precursors (MP1 and MP2) are indicated by blue rectangles.</figcaption></figure>

In the lower panel of Figure 7, the situation of the amplitude sunset anomalies is different from the previous example because the quasi-parallel lines are found to be a combination of all spectral features (i.e., minima, maxima and inflexions). However, two lines appear around −25 min and +50 min, which merge and combine from 31 January to 3 February. The blue rectangle in the lower panel of Figure 7 marks the beginning and the end of the second main precursor MP2. This sunset precursor signature is not evident when compared to the previous one, i.e., the sunrise precursor. One notes that, for both
cases, sunrise and sunset anomalies, the inflections (black circles in Figure 7) mainly emerge close to the Graz facility ($\Delta \sim -30$ min) and to the EQ epicenter ($\Delta \sim +80$ min). The main groups appear at two $\Delta$ intervals: (a) [-50 min, +20 min] with four groups, i.e., group 1 ('+' symbol), group 2 ('*' symbol), group 3 ('◊' symbol), and inflexions ('●' symbol), and (b) [+20 min, +180 min] with three groups, i.e., group 4 ('$\Delta'$ symbol), group 5 ('□' symbol), group 6 ('x' symbol) and group 7 ('♦' symbol). From 2 February to 7 February, we notice an abrupt increase of $\Delta$ from +80 min to about~+180 min in the case of group 7 ('♦' symbol). As for the amplitude sunrise fluctuations, one notes the presence of a prior precursor (PP) when $\Delta \sim -40$ min and $\Delta \sim +40$ min. In the lower panel of Figure 7, the blue rectangles show PP between 23 January and 25 January 2023.

### 3.7.2. Phase Sunrise and Sunset Spectral Anomaly Time Shifts

The phase sunrise anomalies show minima and maxima quasi-parallel lines, as displayed in the upper panel of Figure 8. Such lines are similar to those observed for the amplitude sunrise anomalies (upper panel of Figure 8), with alteration between maxima and minima lines. Irregularities in fluctuations take part in $\Delta$ interval [-180 min, +10 min], and are more organized when $\Delta$ is bigger than +30 min, and particularly at about $\Delta \sim +80$ min. Contrary to the inflection anomalies, the jumps are not at the $\Delta$ edges (i.e., minima and maxima of $\Delta$) but closer to the central part at around +50 min. The precursor signature occurs between 26 January and 1 February and principally from 27 to 31 January, where one can see that the first jump occurs on 29 January, and not before this date, and appears daily until 9 February. In the upper panel of Figure 8, the blue rectangle shows the first main precursor (MP1), as found for amplitude sunrise variations. In the period from 26 January to 1 February, one also remarks that, on the same day 30 January, three minima (red colored symbols) occurred, and three maxima (blue colored symbols), where the jump (black colored symbol) may play the role of a maxima. Lower than $\Delta \sim +10$ min, only group 1 ('+' symbol) and group 2 ('*' symbol) are observed with similar fluctuations after 1 February. Two other groups, i.e., group 3 ('◊' symbol) and group 7 ('●' symbol), are very close to each other, nearly overlapping, from 30 January to 9 February in the $\Delta$ interval [+40 min, +60 min]. Group 4 ('$\Delta'$ symbol) materializes for the first time (nearly quasi-absent in the amplitude anomalies of Figure 7) and corresponds to maxima (blue colored '$\Delta$'). This group has nearly quasi-parallel variations when compared to group 5 ('□' symbol), which is found to contain only minima (red colored '□') anomalies. Both groups are confined to the $\Delta$ interval [+70 min, +100 min]. As for the amplitude sunrise fluctuations, on notes the presence of a prior precursor (PP) when $\Delta \sim 0$ min and $\Delta \sim +80$ min. In the upper panel of Figure 8, the blue rectangles show where PP occurred.

The phase sunset anomalies exhibit a variation confined to Graz terminators, at about $\Delta \sim 0$ min, as shown in the lower panel of Figure 8. Jumps take place at the beginning, from 23 January to 31 January, when maximum anomalies appear from 31 January to 9 February. The second main precursor (MP2) occurs between 31 January and 3 February, as indicated by the blue rectangle in the lower panel of Figure 8. Above the TBB terminators (i.e., $\Delta$ bigger than +50 min), one can see a clear progressive increase of $\Delta$ from 45 min on 24 January, reaching a maximum of about 130 min on 3 February and ending at +90 min 6 days later. Group 1 (i.e., '+' symbol), group 2 (i.e., '*' symbol), group 4 (i.e., '$\Delta$' symbol) and jumps anomalies ('●' symbol) are mainly discovered around the Graz terminators between -10 min and +20 min. Group 5 (i.e., '□' symbol) is quasi-absent around TBB terminators at about +50 min, and only group 6 (i.e., 'x' symbol) and group 7 (i.e., '♦' symbol) principally emerge above +70 min and up to 130 min. In the time interval 23 January to 25 January, one notes a similar variation as found for the phase sunrise anomalies. The prior precursor (PP) happens when $\Delta \sim -10$ min and is shown in a blue rectangle in the lower panel of Figure 8.
# 4. Results and Analysis of VLF Sub-Ionospheric Spectral Anomalies

The main outcomes of our analysis are emphasized in this Section, where we insist on the efficiency of the multi-terminators method. Four aspects are highlighted: (a) the new spectral feature anomalies as the first stage; (b) the presence of the quasi-parallel lines as a second stage to justify further research into the approaching earthquake; (c) the role of jumps and inflexions as spectral anomalies as a third stage to enhance the prior precursor detections; and (d) the physical processes behind the reported precursors.

## 4.1. Minima, Maxima, Inflexions and Jumps as Spectral Feature Anomalies

Detailed analysis of the TBB amplitude and phase signal around the terminators is essential. The daily variations of the transmitter are well-known to be subject to the Earth's rotation, with an alternation between day and night. The crucial transition from night-to-day and day-to-night defines the terminators, respectively, for sunrise and sunset. Enlarging the time window around the terminators has allowed the monitoring of new spectral features, like maxima, inflections and jumps, beside minima features. All these spectral anomalies are the basic ingredients with which to extricate and distinguish between real/authentic and wrong/false precursor(s). Hence a total number of 349 anomalies have been derived from the TBB signal, 204 at sunrise and 145 at sunset, in the period from 23 January to 8 February 2023, i.e., 17 days. It is appropriate with such a high number of anomalies to use the statistical tool approach. This can be foreseen as a first stage in improving the efficiency of the multi-terminators method.

The anomaly occurrence probability quickly brings to light peaks not related to the Graz or TBB terminators but to the EQ epicenter region, as shown in Figure 6. Further, we have introduced groups, which leads to characterization of how the VLF signal propagates among two fundamental geographical points, the emitter (TBB transmitter) and the reception (Graz facility) stations. For this, we have distinguished the anomalies which appear: (a) before (i.e., group 1) and close to Graz terminator (i.e., group 2); (b) between the Graz and TBB terminators (i.e., group 3) and nearby TBB terminators (i.e., group 4); (c) between the TBB and EQ epicenter terminators (i.e., group 5) and near the TBB terminators (i.e., group 5); and (d) finally, before the EQ epicenter terminators (i.e., group 7). Group time intervals have been derived from occurrence probability in Figure 6, which makes evident when precursors start to be noticeable and visible. The Δ time difference linked to each group provided an opportunity to characterize the TBB signal propagation in the Earth's waveguide.

## 4.2. Seismic and Non-Seismic Disturbances of TBB Radio Signal

One needs to insist on the role of the TBB transmitter signal in this analysis. First, the VLF emitted signal is subject to pre-seismic ionospheric disturbances before its detection by the Graz facility. Such disturbances are used as a marker and a tracer to characterize the precursor features. Hence, the reported spectral anomalies are directly linked to these ionospheric disturbances. Second, the TBB transmitter cannot provide information on seismic precursors, despite its close distance to the epicenter. Such powerful emitters are exclusively used for navigation and communication services [3,8], and not for seismic studies. The emitted radio wave is only subject to pre-seismic disturbances when it propagates in the ionospheric waveguide. However, space observations may allow investigation the radio emission beam [53] of the transmitter station and its variability before EQ occurrence. Unfortunately, the altitudes and trajectories of satellites (e.g., DEMETER micro-satellite) cannot lead to continuous and exclusive observations above the transmitter station.

The following steps have been respected to distinguish and separate seismic and non-seismic disturbances of the TBB transmitter signal. Primarily, we have examined the solar
and geomagnetic activities from the beginning of January to the end of February 2023, which have been found to be quiet and moderate. Then, we have investigated a shorter period of 10 days (28 January to 6 February) to check if anomalies occurred before the EQ. We found that two precursors are evident, on 30–31 January and 2–3 February. Then, we agreed to examine the time evolutions of the anomalies before 28 January and after 6 February. We have extended the investigated time interval from 10 days (28 January to 6 February) to 17 days (23 January to 8 February 2023), enabling discovery of the following: (a) the prior precursor (23–25 January) which exhibits weaker signal disturbances when compared to those of the main precursor, and (b) the post-evolution (6 February to 8 February) of the anomalies, which were reduced but continued after the EQ occurrence. This denotes that the selected period of 17 days may be ideal. because it provides details on the main precursor phase and how it grows (23 January–25/26 January), develops (27 January–3 February), and progressively attenuates (3 February to 8 February).

We find that the variation in the Δ time difference versus the observation days exhibits the presence of quasi-regular parallel lines with inherent minima and maxima spectral features. This may be expected as a second stage in improving the efficiency of the multi-terminators method. Hence, one notes, particularly, lines for the amplitude and phase sunrise anomalies when, respectively, Δ is on average at about +35 min and +80 min. Such spectral features seem to be not related to seismic disturbances but to the usual and ordinary TBB wave propagation. However, on some occasions, in the parallel lines, irregularities occur, like drops in minima and absence of maxima. This is the case for the amplitude and phase sunrise anomalies between 27 January and 30 January. Simultaneously and on the same day interval, an increase and decrease in extrema are evident around the Graz terminators in the amplitude sunrise at Δ~±20 min, and for the phase sunrise at Δ~+30 min. The sunset anomalies also show clear irregularities four days later, i.e., 31 January to 3 February. The time difference is about Δ~−25 min and −50 min, respectively, for amplitude and phase anomalies.

Those two consecutive precursors on 27/28 January and 31 January/1 February take place, respectively, 8 days and 4 days before the earthquake occurrence. However, it is necessary to examine whether earlier precursors exist, and how they are linked to specific geographic location(s). This may be comparable to the analysis of Boudjada et al. [48], but using the Δ time difference. A more precise analysis of the spectral anomalies leads to the discovery of similar variations from 23 to 25/26 January: (a) in the amplitude sunrise at Graz (i.e., Δ~[-70 min, -50 min]) and the EQ epicenter (i.e., Δ~[+105 min, +120 min]), (b) in the phase sunrise at Graz (i.e., Δ~[0 min, +10 min]) and the EQ epicenter (i.e., Δ~[+90 min, +100 min]). This indicates the possibility of further investigation of observations which may be considered as of the prior percussors appearing on 23 January, earlier than the main precursors on 27/28 January and 31 January/1 February. Such prior percussors occurring 12-days before the EQ can be assumed as indicators which support extra investigation of an approaching seismic event.

### 4.3. *Inflexion and Jump Spectral Features as Prior Precursor Parameters*

Inflexion and jump anomalies have been observed in the amplitude and phase spectral features, respectively. We remark that only two inflexions occurred on one day, on 31 January, and the other on 1 February, in the amplitude sunset anomalies (i.e., lower panel of Figure 7). We also notice that jumps have only been found to be negative on two occasions, on 29 and 31 January. All these dates belong to the main precursors previously reported, i.e., 8 days and 4 days before the EQ occurrence. The inflexions and jump anomalies also display similar behaviors, i.e., daily occurrence from 23 to 25 January, as in the prior precursor period 23 to 25/26 January.
We have reported in sub-Section 3.4 the absences of inflexions and jump anomalies, respectively, on 27 and 31 January (Table 2), and 26 and 28 January (Table 3). We note that the inflexion amplitude anomalies occurred in the prior precursor period (23–26 Jan at sunset), and disappeared on 27 January at the beginning of the main precursor (e.g., MP1), and again vanished on 31 January (i.e., marking the end of the first main precursor phase). In both cases, we find that maxima occurred on 27 January and 31 January near the Graz facility (around −30 min), when inflexions were not observed. In the case of the jump phase anomalies, and as for inflexions, they occurred in the prior precursor period (23–26 January), vanished on 26 January, came back on 27 January, and again disappeared on 28 January. All these jumps were observed at sunset. With the beginning of the first main precursor, after 28 January, all jumps continuously occurred at sunrise until 8 February. The combination of the presence and the absence of anomalies (inflexions and jumps) leads us to conclude that a precursor process affecting the amplitude and phase reveals different behavior during the prior (e.g., PP) and main precursors (e.g., MP1 and MP2).

The analysis of the slopes for inflexions (sixth column in Table 2) and jumps (sixth column in Table 3) leads us to insist on common behavior. Hence at three occasions, 28 January, 1 February and 4 February for the inflexion slopes, and 27 January, 30 January and 2 February for the jump drifts, maxima have been recorded. We note that the maxima jump drifts occurred earlier with a difference of one/two day(s). These maxima can also be presumed as pointers of the end of the prior precursor (i.e., 25/26), the beginning of the first main precursor (i.e., on 27 January), and finally the end of the second main precursor (i.e., on 3 February).

It is possible to guess that both spectral features, the inflexions and the jumps, are originally minima or maxima anomalies. The inflection anomaly seems to be the result of a minimum, which in a very short time transforms to a maximum, and vice versa, with a difference in intensity level of less than 1 dB. Conversely, the jump anomaly appears as the consequence of the addition of ±180° to a minimum/maximum phase. These spectral features, like the irregular minima and maxima, are all signatures of the seismically disturbed ionosphere, which perturbates TBB transmitter propagation along its ray path before the EQ.

## 4.4. *Physical Processes Behind the Prior and Main Precursors*

The main precursors reported at the end of January and beginning of February are linked to distinct disturbances of the amplitude and the phase of the TBB signal due to ionospheric instabilities 10 days to one week before the EQ occurrence. These precursors have been reported by other studies (see Section 1.4). The prior precursors, as shown in Figures 7 and 8, occur 15 days before the seismic event. We believe that both the prior precursor and the main precursors are generated by the LAIC mechanism, with an increased growth rate starting two weeks before the EQ and reaching a maximum at the end of January 2023, i.e., about one week before the EQ.

However, we note the existence of non-seismic phases corresponding to Δ time intervals close to zero. This is the case for the following specific days for the following: (a) the amplitude sunrise in the upper panel of Figure 7: 27 DOY and 30 DOY; (b) the amplitude sunset in the lower panel of Figure 7: 23 DOY and 34 DOY; (c) the phase sunrise in the upper panel of Figure 8: 24 DOY and 25 DOY; and (d) the phase sunset in the lower panel of Figure 8: 23 DOY and 32 DOY. The presence of non-seismic and seismic days in the investigated period suggests a particular phase process in the LAIC mechanism.

Figure 9 displays the acoustic channel of the LAIC mechanism, where we refer to the model of Yoshida et al. [36] on the relationship between minima/maxima spectral anomalies linked to a drop in/growth of the ionosphere electron density. The left right
panel shows the main physical processes in the atmosphere (i.e., AGWs and GWs), which occur at the origin of the ionosphere pre-seismic disturbances. The right panel provides a connection between the variations in the ionospheric electron density *Ne* and the observed anomalies. In this context, the maxima and the jumps can be the results, respectively, of the growth and the extreme drops in *Ne*. Moderate and significant drop in *Ne* may be the consequence of the minima and inflection anomalies.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left panel shows a schematic representation of the acoustic channel as reported by Hayakawa [3], where the origin of the pre-seismic ionosphere disturbances is due to the seismic preparation zone in the lithosphere, which generates instabilities in the atmosphere, particularly acoustic gravity waves, then later gravity waves, ending in turbulence in the ionosphere electron density *Ne*. The right panel displays the behavior of the investigated spectral anomalies versus the ionospheric electron density fluctuations.</figcaption></figure>

We conclude that, as shown in Figure 9, the features of the Turkey Syria precursors are
as follows:

* the seismic preparation zone exhibits a coupling between the lithosphere, the atmosphere and the ionosphere.

* the atmospheric gravity waves amplify with altitude because of air density reduction.

* the disturbances of the ionosphere electron density by GWs generates anomalies in the transmitter signal detected by the ground-based facility.

* The low/high decrease/increase in ionospheric density produced the reported anomalies, principally extrema (~80%, i.e., 43% minima and 38% maxima), which are explained by Yoshida et al. [36].

* The inflexions and jumps (~20%) are clearly inherent in the same model [36].

**5. Conclusions**

In this study, we attempt to investigate the utility, the realism and the accuracy of
this method. We have computed the anomaly occurrences by separating sunrise and
sunset spectral anomalies. The characterization of the efficiency of the multi-terminators
method leads to forecast of the sunrise and sunset terminators linked to the EQ preparation
zone. This method enables estimation of the time difference Δ between the sunrise/sunset
terminators of the reception station (e.g., Graz facility) and the occurrence time of the
spectral anomalies generated by the pre-seismic areas. The time interval Δ, when converted into geographical longitude, indicates the probable location of the EQ preparation zone. In the case of the Turkey–Syria EQ, we have applied the MTM method and found the prior and main precursors, respectively, 11 days and 5 days before seismic activity in the Anatolian region, well-known as predisposed to seismic activity. We have used the TBB transmitter signal because of its ray path in the pre-seismic sensitive region. principally defined by the Dobrovolsky preparation zone and the fifth Fresnel elliptic area. Ray paths of other European transmitters, like the ITS and ICV Italian transmitters, are found outside this sensitive region. Besides minima and maxima spectral anomalies reported in Boudjada et al. [48], we have added two new features, inflections and jumps, appearing, respectively, in the amplitude and phase of the TBB transmitter signal.

Hence, a total of 349 anomalies have been derived from the MTM method, where 43.0% are minima, 38.4% maxima, 9.7% jumps, and 8.9% inflections. Such spectral anomalies provide a more precise idea of the ways the transmitter signal is perturbated before the EQ occurrence. These pre-seismic anomalies have been recorded in the time interval from 23 January to 8 February 2023. This has allowed us to use the statistical tool approach, where anomaly occurrence probability peaks are a first stage in restricting the location of the preparation EQ zone. The occurrence probability quickly brings to light peaks not related to the Graz or TBB terminators but to the EQ preparation zone. The computed occurrence has been estimated each 12 min, corresponding to a longitude interval of about 3° when sunrise and sunset terminators vary at 04.40–07.60 UT and 13.20–16.80 UT, respectively. The presence of several occurrence peaks provides the possibility of analyzing the performances of spectral anomalies in defined zones. Seven groups have been considered, taking into account the terminator time intervals. For this, we have distinguished anomalies which appear as follows: (a) before (i.e., group 1) and close to the Graz terminator (i.e., group 2); (b) between the Graz and TBB terminators (i.e., group 3) and near the TBB terminators (i.e., group 4); (c) between the TBB and EQ epicenter terminators (i.e., group 5) and near the TBB terminators (i.e., group 6); and (d) finally, before the EQ epicenter terminators (i.e., group 7). This second stage details the behavior of the anomalies in each group versus the investigated time interval, and the estimation of Δ time difference between the Graz terminators and the anomaly time occurrence. The enlarge of the observation time window around the Graz terminators have allowed us to recognize an increase the number of spectral anomalies and to provide an overview of the VLF wave propagation between the emitter station and the reception facility.

This stage has provided the possibility of studying the propagation effects which occur in the detected transmitter signal, where alternance of minima and maxima above the TBB station and the EQ epicenter region, have been observed. These extrema bring an overview of the combined effect of the sky and ground waves along ray paths. However, irregularities (like absence of extrema, or change in anomaly Δ time difference) reveal the presence of two main precursors, from 27 January to 30 January (e.g., MP1), and from 31 January–3 February (e.g., MP2). More important is the prior precursor (e.g., PP) detected from 23 January to 25/26 January, where anomaly fluctuations around the Graz terminators are found to be similar to those approximately situated at the EQ epicenter area. The investigation of the maximum slopes related to the inflexions and the jumps leads us to consider that such drifts are markers of the end of the prior precursor (i.e., 25/26), the beginning of the first main precursor (i.e., on 27 January), and finally the end of the second main precursor (i.e., on 3 February).

A forecasting model can be suggested within the frame of this work, where the following steps can be observed: (a) a choice of a powerful VLF transmitter with a preferable constant phase variation; (b) transmitter ray paths included in the pre-seismic sensitive
region defined by the Dobrovolsky preparation zone and Fresnel elliptic area; (c) a daily analysis of spectral features (minima, maxima, inflexions and jumps), particularly those with $\Delta$ time difference more than $\pm$12 min; (d) an anomaly occurrence probability based on sunrise and sunset terminators of the transmitter station and the VLF reception facility; (e) A study of $\Delta$ time difference versus day of the year where transmitter and reception terminators are indicated. These steps can provide, in the presence of spectral anomalies, first hints on the longitudinal locations of the seismic preparation zone, particularly when simultaneous anomalies like prior precursors (as shown in Figures 7 and 8) occur in the amplitude sunrise and sunset, and also in the phase sunrise and sunset. The daily change in the jump drift rates, from positive to negative is a good indicator for the near-future start of the main precursors.

Future perspectives will be devoted to the comprehension of the relationship between the prior and the main precursors. This can only be achieved by applying the multi-terminator method to other EQs. Of course, the new reception system (i.e., UltraMSK-2) used for this investigation and the powerful transmitter signal (i.e., TBB transmitter) play a crucial role in the spectral anomaly detection. All these ingredients may provide a clear hint on the origin of seismic precursors, and a possibility of earthquake forecast by radio observations in the very low frequency band. In the near future, new UltraMSK-2 stations are planned to be installed in Guyancourt (France), Reunion and Sri Lanka as reported by Galopeau et al. [51] permitting an enlargement of the VLF network. New collaborations are also intended with other networks, like AWESOME, SAVNET, teams in Japan and in India, and VLF communities like lightning location network (i.e., WWLLN). The synergy between the VLF ground-based networks (e.g., INFREP network) and space-based satellites (e.g., CSES1, CSES2, Swarm) can lead to a better comprehension of the ionospheric pre-seismic perturbation behaviors, particularly when simultaneous observations can be performed.

**Author Contributions:** Conceptualization, methodology, M.Y.B.; original draft preparation, M.Y.B., H.U.E. and G.N.; investigations, P.H.M.G., M.C. and P.F.B.; formal analysis, S.S.; data curation, H.U.E. and W.V., supervisions, H.L. and W.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data used in this work are available upon request to the corresponding author. The data are not publicly available due to privacy.

**Acknowledgments:** The authors are thankful to the editors and unknown reviewers for precious and valuable suggestions and constructive comments.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Gokhberg, M.B.; Morgounov, V.A.; Yoshino, T.; Tomizawa, I. Experimental measurement of electromagnetic emissions possibly related to earthquakes in Japan. *J. Geophys. Res.* **1982**, *87*, 7824–7828. [CrossRef]

2. Warwick, J.W.; Stoker, C.; Mayer, T.R. Radio emissions associated with rock fracture: Possible application to the great Chilean earthquake of May 22, 1980. *J. Geophys. Res.* **1982**, *87*, 2851–2859. [CrossRef]

3. Hayakawa, M. *Earthquake Prediction with Radio Techniques*; John Wiley and Sons: Singapore, 2015; 296p.

4. Molchanov, O.A.; Hayakawa, M. *Seismo Electromagnetics and Related Phenomena: History and Latest Results*; Terrapub: Tokyo, Japan, 2008; 189p.

5. Pulinets, S.A.; Boyarchuk, K. *Ionospheric Precursors of Earthquakes*; Springer: Berlin/Heidelberg, Germany, 2004; 315p.

6. Gokhberg, M.B.; Morgounov, V.A.; Pokhotelov, O.A. *Earthquake Prediction. Seismoelectromagnetic Phenomena*; Gordon and Breach Science Publishers: Amsterdam, The Netherlands, 1995; 287p.

7. Davies, K. *Ionospheric Radio*, 3rd ed.; Institution of Engineering and Technology: London, UK, 2008; 603p.
8. Hunsucker, R.D.; Hargreaves, J.K. *The High-Latitude Ionosphere and Its Effects on Radio Propagation*; Cambridge University Press: Cambridge, UK, 2003; 639p.

9. Gokhberg, M.B.; Gufeld, I.L.; Rozhnoy, A.A.; Marenko, V.F.; Yampolsky, V.S.; Ponomarev, E.A. Study of seismic influence on the ionosphere by super long-wave probing of the Earth-ionosphere waveguide. *Phys. Earth Planet. Inter.* **1989**, *57*, 64–67. [CrossRef]

10. Hayakawa, M.; Molchanov, O.A.; Ondoh, T.; Kawai, E. Precursory signature of the Kobe earthquake on VLF subionospheric signal. *J. Atmos. Electr.* **1996**, *16*, 247–257.

11. Molchanov, O.A.; Hayakawa, M. Subionospheric VLF signal perturbations possibly related to earthquakes. *J. Geophys. Res. Space Phys.* **1998**, *103*, 17489–17504. [CrossRef]

12. Molchanov, O.; Fedorov, E.; Schekotov, A.; Gordeev, E.; Chebrov, V.; Surkov, V.; Rozhnoi, A.; Andreevsky, S.; Iudin, D.; Yunga, S.; et al. Lithosphere-atmosphere-ionosphere coupling as governing mechanism for preseismic short-term events in atmosphere and ionosphere. *Nat. Hazards Earth Syst. Sci.* **2004**, *4*, 757–767. [CrossRef]

13. Kilimenko, M.V.; Kilimenko, V.V.; Karpov, I.V.; Zakharenkova, I.E. Simulation of seismo-ionospheric effects initiated by internal gravity wave. *Russ. J. Phys. Chem.* **2011**, *B5*, 393–401. [CrossRef]

14. Hayakawa, M.; Kasahara, Y.; Nakamura, T.; Hobara, Y.; Rozhnoi, A.; Solovieva, M.; Molchanov, O.A.; Korepanov, K. Atmospheric gravity waves as a possible candidate for seismo-ionospheric perturbations. *J. Atmos. Electr.* **2011**, *31*, 129–140. [CrossRef]

15. Korepanov, V.; Hayakawa, M.; Yampolski, Y.; Lizunov, G. AGW as a seismo-ionospheric coupling responsible agent. *Phys. Chem. Earth.* **2009**, *34*, 485–495. [CrossRef]

16. Molchanov, O.A.; Hayakawa, M.; Miyak, K. VLF/LF sounding of the lower ionosphere to study the role of atmospheric oscillations in the lithosphere-ionosphere coupling. *Adv. Polar Upper Atmos. Res.* **2001**, *15*, 146–158.

17. Hayakawa, M. On the fluctuation spectra of seismo-electromagnetic phenomena. *Nat. Hazards Earth Syst. Sci.* **2011**, *11*, 301–308. [CrossRef]

18. Yasuoka, Y.; Shinogi, M. Anomaly in atmospheric radon concentration: A possible precursor of the 1995 Kobe, Japan, earthquake. *Health Phys.* **1997**, *72*, 759–761. [CrossRef]

19. Yasuoka, Y.; Igarashi, G.; Ishikawa, T.; Tokonami, S.; Shinogi, M. Evidence of precursor phenomena in the Kobe earthquake obtained from atmospheric radon concentration. *Appl. Geochem.* **2006**, *21*, 1064–1072. [CrossRef]

20. Pulinets, S.; Ouzounov, D.; Ciraolo, L.; Singh, R.; Cervone, G.; Leyva, A.; Dunajecka, M.; Karelin, A.; Boyarchuk, K. Thermal, atmospheric and ionospheric anomalies around the time of Colima M7.8 earthquake of January 21, 2003. *Ann. Geophys.* **2006**, *24*, 835–849. [CrossRef]

21. Pulinets, S.; Ouzounov, D. Lithosphere-atmosphere-ionosphere coupling (LAIC) model: An unified concept for earthquake precursors validation. *J. Asian Earth Sci.* **2011**, *41*, 371–382. [CrossRef]

22. Parrot, M.; Tramutoli, V.; Liu, J.Y.; Pulinets, S.; Ouzounov, D.; Genzano, N.; Lisi, M.; Hattori, K.; Namgaladze, A. Atmospheric and ionospheric coupling phenomena associated with large earthquakes. *Eur. Phys. J. Special Topics* **2021**, *230*, 197–225. [CrossRef]

23. Pulinets, S. Physical bases of the short-term forecast of earthquakes. *Astron. Astrophys. Trans.* **2023**, *34*, 65–84. [CrossRef]

24. Hayakawa, M.; Hirooka, S.; Michimoto, K.; Potirakis, S.M.; Hobara, Y. Meteorological anomalies during earthquake preparation: A case study for the 1995 Kobe earthquake (M = 7.3) based on statistical and machine learning-based analyses. *Atmosphere* **2025**, *16*, 88. [CrossRef]

25. Ray, S.; Chakrabarti, S.K. A study of the behavior of the terminator time shifts using multiple VLF propagation paths during the Pakistan earthquake (M = 7.2) of 18 January 2011. *Nat. Hazards Earth Syst. Sci.* **2013**, *13*, 1501–1506. [CrossRef]

26. Chakrabarti, S.K.; Saha, M.; Khan, R.; Mandal, S.; Acharyya, K.; Saha, R. Possible detection of ionospheric disturbances during Sumatra-Andaman Islands. *Indian J. Radio Space Phys.* **2005**, *34*, 314–317.

27. Sasmal, S.; Chakrabarti, S.K.; Chakrabarti, S. Studies of the correlation between ionospheric anomalies and seismic activities in the Indian subcontinent. *AIP Conf. Proc.* **2010**, *1286*, 270–290.

28. Ghosh, S.; Chakraborty, S.; Sasmal, S.; Basak, T.; Chakrabarti, S.K.; Samanta, A. Comparative study of the possible lower ionospheric anomalies in very low frequency (VLF) signal during Honshu, 2011 and Nepal, 2015 earthquakes. *Geomat. Nat. Hazards Risk* **2019**, *10*, 1596–1612. [CrossRef]

29. Sasmal, S.; Chakrabarti, S.K.; Ray, S. Unusual behavior of VLF signals observed from Sitapur during the Earthquake at Honshu Japan on 11 March, 2011. *Indian J. Phys.* **2014**, *88*, 103–119. [CrossRef]

30. Wait, J.R. *Electromagnetic Waves in Stratified Media*; Pergamon Press: New York, NY, USA, 1962; 608p.

31. Politis, D.Z.; Sasmal, S.; Hayakawa, M.; Haralambous, H.; Datta, A.; Potirakis, S.M. A six-year (2014–2020) statistical correlation study of VLF terminator time shift with earthquakes in Japan. *Remote Sens.* **2024**, *16*, 4162. [CrossRef]

32. Sasmal, S.; Chakrabarti, S.K. Ionospheric anomaly due to seismic activities—Part 1: Calibration of the VLF signal of VTX 18.2 kHz station from Kolkata and deviation during seismic events. *Nat. Hazards Earth Syst. Sci.* **2009**, *9*, 1403–1408. [CrossRef]

33. Lowrie, W. *Fundamentals of Geophysics*; Cambridge University Press: Cambridge, UK, 2007; 393p.

34. Shalimov, S.L. Lithosphere-ionosphere relationship: A new way to predict earthquakes? *Episodes. Int. Geophys. Newsmag.* **1992**, *15*, 252–254. [CrossRef]
35. Garmash, S.V.; Lin'kov, E.M.; Petrova, L.N.; Shved, G.M. Generation of atmospheric oscillations by seismic gravitational oscillations of the Earth. *Izvessitya Atmos. Ocean. Phys.* **1989**, *25*, 1290–1299.

36. Yoshida, M.; Yamauchi, T.; Horie, T.; Hayakawa, M. On the generation mechanism of terminator times in subionospheric VLF/LF propagation and its possible application to seismogenic effects. *Nat. Hazards Earth Syst. Sci.* **2008**, *8*, 129–134. [CrossRef]

37. Hazelgrove, J. *Ray Theory and a New Method for Ray Tracing*; Physics of the Ionosphere; The Physical Society: London, UK, 1955; 355p.

38. Rapoport, Y.; Reshetnyk, V.; Grytsai, A.; Grimalsky, V.; Liashchuk, O.; Fedorenko, A.; Hayakawa, M.; Krankowski, A.; Błaszkiewicz, L.; Flisek, P. Spectral analysis and information entropy approaches to data of VLF disturbances in the waveguide Earth-Ionosphere. *Sensors* **2022**, *22*, 8191. [CrossRef]

39. Haken, H.; Fraser, A.M. Information and self-organization: A macroscopic approach to complex systems. *Am. J. Phys.* **1989**, *57*, 958–959. [CrossRef]

40. Potirakis, S.; Minadakis, G.; Eftaxias, K. Relation between seismicity and pre-earthquake electromagnetic emissions in terms of energy, information and entropy content. *Nat. Hazards Earth Syst. Sci.* **2012**, *12*, 1179–1183. [CrossRef]

41. De Santis, A.; Abbattista, C.; Alfonsi, L.; Amoruso, L.; Campuzano, S.A.; Carbone, M.; Cesaroni, C.; Cianchini, G.; De Franceschi, G.; De Santis, A.; et al. Geosystemics View of Earthquakes. *Entropy* **2019**, *21*, 412. [CrossRef]

42. Barbot, S.; Luo, H.; Wang, T.; Hamiel, Y.; Piatibratova, O.; Javed, M.T.; Braitenberg, C.; Gurbuz, G. Slip distribution of the February 6, 2023 Mw 7.8 and Mw 7.6, Kahramanmaraş, Turkey earthquake sequence in the East Anatolian Fault Zone. *Seismica 2023*, *2*, 1–17. [CrossRef]

43. Contadakis, M.E.; Pikridas, C.; Bitharis, S.; Scordilis, E. TEC variation over Europe during the intense tectonic activity in the area of Turkey on February of 2023. In *Proceedings of the European Geosciences Union General Assembly 2024*, Vienna, Austria, 14–19 April 2024. EGU24-6275.

44. Contadakis, M.E.; Arabelos, D.N.; Vergos, G.; Spatalas, S.D.; Scordilis, E.M. TEC variations over the Mediterranean before and during the strong earthquake (M = 6.5) of 12th October 2013 in Crete, Greece. *Phys. Chem. Earth* **2015**, *85*, 9–16. [CrossRef]

45. Vesnin, A.; Yasyukevich, Y.; Perevalova, N.; Sentürk, E. Ionospheric response to the 6 February 2023 Turkey–Syria earthquake. *Remote Sens.* **2023**, *15*, 2336. [CrossRef]

46. Akhoondzadeh, M.; Marchetti, D. Study of the preparation phase of Turkey's powerful earthquake (6 February 2023) by a geophysical multi-parametric fuzzy inference System. *Remote Sens.* **2023**, *15*, 2224. [CrossRef]

47. Liu, J.; Zhang, X.; Yang, M.; Yang, Y.; He, F.; Xue, L.; Yao, X.; Yang, X.; Wu, W.; Qiu, G. Pre-seismic anomaly analysis of the Turkey earthquakes on 6 February 2023 based on multi-source satellite observations. *Nat. Hazards* **2024**, *120*, 12491–12513. [CrossRef]

48. Boudjada, M.Y.; Biagi, P.F.; Eichelberger, H.U.; Nico, G.; Schwingenschuh, K.; Galopeau, P.H.M.; Solovieva, M.; Contadakis, M.; Denisenko, V.; Lammer, H.; et al. Unusual sunrise and sunset terminator variations in the behavior of sub-ionospheric VLF phase and amplitude signals prior to the Mw7.8 Turkey Syria earthquake of 6 February 2023. *Remote Sens.* **2024**, *16*, 4448. [CrossRef]

49. Boudjada, M.Y.; Biagi, P.F.; Eichelberger, H.U.; Nico, G.; Galopeau, P.H.M.; Ermini, A.; Solovieva, M.; Hayakawa, M.; Lammer, H.; Voller, W.; et al. Analysis of pre-seismic ionospheric disturbances prior to 2020 Croatian earthquakes. *Remote Sens.* **2024**, *16*, 529. [CrossRef]

50. Dobrovolsky, I.P.; Zubkov, S.I.; Myachkin, V.I. Estimation of the size of earthquake preparation zones. *Pure Appl. Geophys.* **1979**, *117*, 1025–1044. [CrossRef]

51. Galopeau, P.H.M.; Maxworth, A.; Boudjada, M.Y.; Eichelberger, H.U.; Meftah, M.; Biagi, P.F.; Schwingenschuh, K. A VLF/LF facility network for preseismic electromagnetic investigations. *Geosci. Instrum. Method. Data Syst.* **2023**, *12*, 231–237. [CrossRef]

52. Fuller, W.A. *Introduction to Statistical Time Series*; John Wiley and Sons: New York, NY, USA, 1976.

53. Boudjada, M.Y.; Eichelberger, H.U.; Al-Haddad, E.; Magnes, W.; Galopeau, P.H.M.; Zhang, X.; Pollinger, A.; Lammer, H. Case study of radio emission beam associated to very low frequency signal recorded onboard CSES satellite. *Adv. Radio Sci.* **2023**, *20*, 77–84. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.