© K. Themelis et al., Published by EDP Sciences 2025

1 Introduction

Large-Scale Travelling Ionospheric Disturbances (LSTIDs) are generally considered to be the manifestation of the Atmospheric Gravity Waves (AGWs) in the ionosphere and are associated with auroral and geomagnetic activity (Hunsucker, 1982). Joule heating and Lorentz forces caused by the enhancement of the auroral electrojet or intense precipitation of charged particles are the primary excitation sources of the LSTIDs. Ferreira et al. (2020) reported more frequent and higher-amplitude LSTIDs during periods of intensified auroral electrojet currents. This general correspondence – that LSTIDs are most prevalent immediately following peaks in auroral current – suggest that the ionospheric convection pattern variations are pivotal in exciting TIDs. Essentially, when ion convection intensifies (raising the electric field and plasma flow speed in the auroral zone), the resulting currents and energy dissipation increase, launching gravity waves. LSTIDs have typical horizontal velocities between 300 and 1000 m/s, horizontal wavelengths of more than 1000 km, and periods in the range between 30 min and 3 h. In general, in the Northern Hemisphere, LSTIDs have a south to south-east propagation direction, since they are generated in the northern auroral zone. There are rare cases of poleward propagating LSTIDs observed either in connection with energy dissipation from Medium Scale Travelling Ionospheric Disturbances (MSTIDs) of lower atmosphere origin (Ding et al., 2013) or in connection with interhemispheric constructive interference (Guo et al., 2014).

LSTIDs can impose disturbances with amplitudes of up to ~20% of the ambient electron density, a Doppler frequency shift of the order of 0.5 Hz on HF signals (Reinisch et al., 2018), and perturbations in the Total Electron Content (TEC) from less than 1 TEC unit (TECU, which is equal to 10¹⁶ el/m²) up to 10 TECU (Hernández-Pajares et al., 2006). These disturbances in ionospheric characteristics may cause effects on the operation of systems using transionospheric propagation or those that use the ionosphere as a propagation medium. For example, LSTIDs affect the performance of Global Navigation Satellite Systems (GNSS), and in particular, the Satellite-Based Augmentation Systems (SBAS), such as the European SBAS European Geostationary Navigation Overlay Service (EGNOS) (Pintor & Roldán, 2015), as they can produce variations in TEC of several TECUs, resulting in a decrease of the observed accuracy and a limitation of the availability of these navigation systems for the different applications that they support (mainly aviation). The presence of LSTIDs introduces difficulties to the correct interpretation of radio astronomy observations, such as the Low Frequency Array (LOFAR), which operates at frequencies from 10 to 250 MHz (Van Haarlem et al., 2013). LSTIDs are also the principal remaining mechanism limiting the performance of target detection algorithms for geolocation systems, because the associated range/azimuth deflections spread the target return in the algorithm results (Nickisch et al., 2006). The effect of LSTIDs on the Maximum Usable Frequency for a 3000 km radio path via reflection from the F2 layer (MUF(3000)F2) and consequently on High Frequency (HF) propagation paths is also demonstrated by Altadill et al. (2020a, 2020b).

Within the framework of the European Project “Warning and Mitigation Technologies for Travelling Ionospheric Disturbances Effects” (TechTIDE) funded by the Horizon 2020 Programme, new methodologies have been designed, developed, and deployed to detect the occurrence of travelling ionospheric disturbances in real-time and their potential impact on ground-based and aerospace applications, using activity metrics derived from degradation data of affected applications (Belehaki et al., 2020). The analysis of the TID effects in the applications concerned has confirmed the need not only for nowcasting information but also for early warnings, short-term and long-term forecasts, that can support mitigation technologies capable of preventing degradation of the concerned applications. However, the forecast of LSTIDs remains a requirement for researchers and users, with several issues still needing to be resolved, mainly arising from the complex triggering mechanisms and the multiparametric propagation pattern. The follow-up project “Travelling Ionospheric Disturbances Forecasting System” (T-FORS), funded by Horizon Europe, proposes for the first time LSTID forecasting models based on Machine Learning (ML) classification algorithms. In the framework of the T-FORS project, a first ML model was proposed to forecast the occurrence of LSTIDs over the European continent up to 3 h in advance (Ventriglia et al., 2025). The model provides a single forecast representative for the whole European region. The model, based on Categorical Boosting (CatBoost) – a gradient boosting library – and trained on a human-validated LSTID catalogue (Segarra et al., 2024), uses a diverse set of physical drivers, ranging from geomagnetic indices, solar wind observations, solar activity data, to ionosonde measurements. The model results are compared with LSTID characteristics derived with independent methods, demonstrating a high predictive robustness. The quantitative assessment is performed for distinct operating modes, namely balanced, high precision, and high accuracy. Accuracy is defined as the percentage of correct predictions for the test data. Precision is defined as the fraction of relevant examples. The validation for the balanced mode identifies 51% of the LSTIDs reported in the source catalogue, while the high-precision mode correctly identifies 36% of the reported LSTIDs. Finally, with high-accuracy mode, the model correctly identifies 60% of the reported LSTIDs. This level of performance might be due to the high variability of LSTID characteristics at different European sites, which cannot be correctly represented with a single number for the whole European region, indicating the need for single-site LSTID forecasting.

This contribution proposes a new forecasting model to estimate the level of LSTIDs activity at specific locations at middle latitudes in Europe. We restrict the forecasting scheme to high-latitude sources, where the most likely triggering candidates are the Joule heating and Lorentz forces associated with the auroral electrojet and intense particle precipitation events. Among these mechanisms, the surges in the auroral electrojet are the most important sources of Atmospheric Gravity Wave (AGW) driven TIDs, as clearly explained in the Hunsucker (1982) review paper. The localized heating caused by intense precipitation of charged particles is the second most effective source mechanism of AGWs at high latitudes. Lyons et al. (2019) found evidence that LSTIDs detection is more likely near the longitude range of the initiating disturbance than further away. This is consistent with the expectation that it should be more likely to detect LSTIDs near the longitude range of the initiating auroral disturbance than several hours away in Magnetic Local Time (MLT). The forecasting modeling framework presented here uses as drivers of LSTIDs the Auroral Electrojet (AE) Indices derived from the International Monitor for Auroral Geomagnetic Effects (IMAGE) network, instead of the global AE indices, to introduce more accurately the concept of the localized longitudinal range of the excitation sources. A supplementary source is required to detect AGW propagating and interacting with the middle latitude ionosphere, and this is provided by the TEC gradients. Finally, a third set of data provides computed characteristics of LSTIDs over Digisonde locations at European middle latitudes, acquired from the TechTIDE database, as described in detail in the sections below. It is therefore possible to build an ML forecasting framework using features derived from both the auroral oval and mid-latitude regions, whether observed directly or calculated. Following recent trends in ML for time series forecasting, we aspire to employ novel tools to forecast the occurrence of LSTIDs and compare their performance against more classical ML approaches.

Section 2 describes the data used for training and validation of the ML classifiers. Section 3 provides the concept underpinning the development of ML classifiers tested in this contribution. Section 4 presents the results, which are further analysed and discussed in Section 5. Finally, the conclusions and outlook are presented in Section 6.

2 Data

The data series used in this study are drivers of LSTIDs and their detected characteristics. Considering that the proposed model is restricted to the forecast of LSTIDs occurrence having their sources in the high-latitude ionosphere, the drivers are magnetic field disturbances recorded in the European sector of the auroral region, and TEC gradients inferred from TEC maps covering Europe.

Magnetometer data are obtained from the IMAGE network, a dense array of magnetometers located at subauroral and auroral latitudes (approximately 55° to 80° geomagnetic latitude), providing good regional coverage across European longitudes. The derived electrojet indicators (Kallio et al., 2000) IL (IMAGE Lower), IU (IMAGE Upper), and IE (IMAGE Electrojet), are simple estimates of the total eastward and westward currents crossing the magnetometer network, based on the deviation of the H-component (the horizontal northward magnetic field component) of the magnetic field from its quiet time baseline, recorded at each station. The IL index is the minimum value of the H-component (negative perturbations) from all stations and represents the intensity of the westward electrojet. Similarly, the IU index takes the maximum value of the H-component (positive perturbations) from all the auroral latitude stations and represents the eastward electrojet.

Kauristie et al. (1997) demonstrated that in some cases, the auroral oval is either so expanded or contracted that the latitudinal coverage of the AE magnetometer chain is not adequate. The coverage provided by the IMAGE magnetometers network follows the transient variations of the oval much better and this is an additional argument in favor of the use of these indicators as representatives of the electrojets’ intensity in the European sector.

The second driver considered in this model is TEC gradient inferred from TEC maps covering Europe (Borries et al., 2017). The estimation of TEC gradients is based on TEC maps provided by the German Aerospace Centre (DLR) for the European region. TEC is first converted into equivalent delay units (in millimeters, mm) because it relates directly to the impact on GNSS signal paths. Then, the spatial gradient is calculated, i.e., the rate of change of TEC delay with distance across the ionosphere. Thus, the units become millimeters per kilometer (mm/km). Borries et al. (2017) showed that sharp gradients in TEC mark regions where LSTIDs are likely generated. Therefore, TEC gradients are considered precursors of LSTID activity. Steep TEC gradients often form near the equatorward edge of the auroral oval during storms, indicating substantial ionospheric convection changes. The European maps of TEC gradients are provided in real-time through the TechTIDE service of the ESA Space Weather Network.

Since the generation of TEC maps averages out steep TEC gradients, relatively low thresholds must be assumed for the indication of the probability of LSTID generation. The statistical analysis of TEC gradients shows that during quiet conditions, TEC gradients are small (about 0.2 mm/km), but during storm events, gradients exceed 2 mm/km. In this contribution, the TEC gradient index, which is the average TEC gradient over Europe, is used as an equivalent of a LSTID precursor index.

Finally, the third time series that is used as an input feature to the forecasting model is associated with the LSTID detections. This is the Spectral Energy Contribution (SEC) of the detected Travelling Ionospheric Disturbances, extracted from the results of the HF-Interferometry (HF-INT) method. The method is designed to detect LSTIDs in near real-time using monostatic measurements from a network of HF sensors (i.e., ionosondes), as described in detail by Altadill et al. (2020a). The method is transitioned to operations as part of the TechTIDE project (Belehaki et al., 2020), and the results are provided in real-time from the TechTIDE service of the European Space Agency Space Weather Network. The HF-INT method utilizes ionospheric characteristics, specifically the Maximum Usable Frequency, MUF(3000)F2, obtained from autoscaled ionograms of ten European ionosondes, available through the Digital Ionogram Data Base (DIDBase) of the Global Ionospheric Radio Observatory (GIRO)¹. To effectively identify LSTIDs, the spatial distribution of the ionosonde network must be sufficiently dense, with the distance between measuring sites not exceeding 1000 km.

The HF-INT method applies spectral analysis (Scargle, 1982; Hocke, 1998) to the detrended MUF(3000)F2 time series, ΔMUF, after removing the main daily harmonics. HF-INT computes the amplitudes (A(T)) of the periodic components (T) of the time series for a 6-h long time interval and estimates the dominant period (TAM) of the LSTID-like variations, ensuring a high confidence level (>0.95) and coherent TAM values within the different Digisonde sensors. The method requires data from at least three sensors that exhibit coherent periodicities and high cross-correlation values. If the above conditions are met, HF-INT estimates the propagation velocity vector of LSTID by measuring time delays of the disturbance across different sensor sites, assuming plane wave propagation. However, if the above conditions are not met, HF-INT returns no LSTID detected at that time, resulting in output values marked as NaN. The characteristics of the detected LSTIDs are referenced to the last time instant of the 6-h analysis window. For example, LSTIDs features obtained in a window from 0300UT to 0900UT are referenced to 0900UT. HF-INT estimates the contribution of the TID-like variation to the total variability of the time series of the different Digisonde sensors, the so-called Spectral Energy Contribution (SEC) of the disturbance, which is estimated by the percentage ratio of the spectral energy of the periodic range attributed to the LSTID to the total spectral energy (Eq. (1))

$$\mathrm{S}\mathrm{EC}=\frac{\sum _{T=T_{{\mathrm{TID}}_{\mathrm{S}}}}^{T=T_{{\mathrm{TID}}_{\mathrm{E}}}}{A\left(T\right)}^2}{{\sum }_{T=T_{\mathrm{S}}}^{T=T_{\mathrm{E}}}{A\left(T\right)}^2}\times 100$$(1)

where $T_{{\mathrm{TID}}_{\mathrm{S}}}$ and $T_{{\mathrm{TID}}_{\mathrm{E}}}$ are the starting and ending bounds of the periodic range of the TID-like variation, respectively, and T_S and T_E are the starting and ending periods of the total periodic range under analysis.

Note that the spectral energy for a given periodic range can be calculated as a sum of squared amplitudes A(T). Moreover, the quantitative evaluation of the contribution of the variations of an arbitrary periodic range to the total variation of the data series can be expressed as the percentage ratio of the signal energy in this periodic range to the total signal energy (e.g., Altadill et al., 1998; Apostolov et al., 1998), and according to the Parseval’s relation, the percentage ratio of the signal energy is equivalent to the percentage ratio of the spectral energy of this periodic range to the total spectral energy. In the HF-INT method, T_S = 30 min and T_E = 360 min, and $T_{{\mathrm{TID}}_{\mathrm{S}}}=0.8$TAM and $T_{{\mathrm{TID}}_{\mathrm{E}}}=1.4$TAM, where TAM is the period of the Maximum Amplitude. The classification of TID activity in the HF-INT method is based on the SEC of the detected TID. A higher SEC indicates a greater impact on variability. Altadill et al. (2019, 2020a) performed a statistical analysis of events detected in 2018 in the European region, categorizing levels of activity based on the SEC distribution. The activity levels are defined as Insignificant (SEC < 18%), Weak (18% ≤ SEC < 65%), Moderate (65% ≤ SEC < 80%), Strong (80% ≤ SEC < 86%), and Very Strong (SEC ≥ 86%). An example of SEC calculation, from the MUF data time series to SEC results in percentage is presented in Figure 1, for data obtained from Dourbes Digisonde (DB049) on April 23, 2023, at 1800UT.

Figure 1 On the top panel, MUF (MUF(3000)F2) data series from the 24 previous hours, from April 22 at 1800UT to April 23 at 1800UT. The red curve corresponds to the fit applied to remove the main daily trends for the last 6 hours, from April 23 at 1200UT to April 23 at 1800UT. The second panel shows the residuals, ΔMUF, for the last 6  h, the difference between the original data series and the fit. The third panel shows the periodogram of the last 6 h, with a dominant period around 130 min. Painted in magenta, the area defined by 0.8 TAM and 1.4 TAM, where TAM is the maximum amplitude period, in this case around 130 min. The final value of SEC is obtained from the ratio between the area painted and the total area under the curve.

An example of SEC calculation, from the MUF data time series to the final value of SEC in percentage, is presented in Figure 1, for data obtained from Dourbes Digisonde DB049 on April 23, 2023, at 1800UT.

HF-INT SEC calculations for four Digisonde locations (Table 1) are used in this contribution.

Figure 2 shows an example of the variations of the characteristics of the LSTID for April 23, 2023, over Dourbes Digisonde DB049 provided by the HF-INT method. The sharp peaks in the SEC are in fact spikes in the time series derived. The latter is because HF-INT data is a real-time estimate subject to the real-time MUF(3000)F2 data availability and quality. In addition, while SEC calculation originates from individual sensors, the identification of an LSTID at a reference sensor at a specific time relies on the coherency of HF-INT estimates between that sensor and its neighbors. Note that in Figure 2, all LSTID characteristics (period, SEC, azimuth, and velocity) are relatively stable in time from 1500UT to 2000UT, whereas large scattering in magnitude is observed outside that time interval. Some drift in the magnitude of these parameters is also observed at the beginning and ending of the detection interval. This is because the spectral analysis requires that the window catches at least “half-period” of the wave-like disturbance.

Figure 2 Time variations of the characteristics of the LSTID detected by the HF-INT at Dourbes, Belgium, for April 23, 2023. The top plot depicts the activity levels of LSTIDs (green, yellow, orange, red, and magenta for Insignificant, Weak, Moderate, Strong, and Very Strong activity, respectively). The middle plot depicts the dominant period (black diamonds) and SEC (green diamonds), and bottom plots depict the magnitude of the velocity (red diamonds) and azimuth direction (blue squares) of propagation of the detected LSTID.

Figure 3 presents the response of the TEC gradients index, inferred from TEC maps that are based on ground based GNSS receiver data (and refer to hereafter as TEC Grad index) and that of the HF-INT SEC parameter over Juliusruh, Dourbes and Sopron Digisonde stations to a moderate geomagnetic storm occurred from 26 to 28 February, 2023, and the associated disturbances recorded by the IMAGE magnetometer chain. This storm triggered an intensification in segments of the westward and the eastward electrojet circulating in the auroral oval at the Northern European longitudes, and recorded by the IL and IU indices, respectively. The TEC Grad index increased around 2000UT on 26 February 2023, almost simultaneously with the IU and IL indices, but more gradually, indicating intense precipitation of charged particles and generation of AGWs near the equatorial edge of the auroral oval. LSTID ws detected by the HF-INT SEC at Juliusruh, Dourbes, and Sopron Digisonde stations several hours later and mainly from 0600UT until 2400UT on 27 February 2023. Some weaker LSTID activity is detected on 28 February 2023 after 1000UT as the result of a weak auroral currents intensification depicted in the plots of IU, IL, and TEC Grad indices from 0400UT to 0900UT the same day. The time sequence of the evolution of the phenomenon, from its triggering at higher latitudes to its detection at lower latitudes, cannot be specified precisely due to the 30° longitudinal spread of the observing points of Digisonde locations; however, a large-scale consistency is observed.

Figure 3 The response of the LSTID drivers (IL and IU indices) and precursors (TEC Grad index) during the geomagnetic storm that occurred from 26 to 28 February 2023. The SEC calculations, using ionospheric characteristics recorded at Juliusruh, Dourbes, and Sopron Digisonde stations, are presented in the bottom panels.

The dataset used for model training and testing spans 18 months. Specifically, the period from January 1, 2022, to April 30, 2023, was used for model training, while data from May 1, 2023, to June 30, 2023, served as the testing set to evaluate the model’s performance. While the IL and IU data series are complete and free of outliers, the TEC gradient index and SEC calculations from Digisonde stations contain a significant proportion of missing values, as detailed in Table 2. Within the training dataset, these missing values were addressed through interpolation (as analyzed in Sect. 3). However, for the testing dataset, the missing values were explicitly disregarded during the performance evaluation.

Figure 4 displays the distribution of the whole SEC time series calculated over the Ebro Digisonde station. The histogram counts the SEC occurrences exceeding 50% (blue), and it is superimposed by the corresponding histogram of SEC exceeding 70% (orange), using 24 bins each representing an hour of the day. “No data” values are disregarded. The whole time series data spans a period from January 1, 2022, to June 30, 2023. In this histogram, we see a peak of LSTID occurrence during late evening and early morning. This agrees with the climatological characteristics reported by Tsugawa et al. (2004), Ferreira et al. (2020), and Nykiel et al. (2024). During daytime hours, the dayside cusp (especially near magnetic noon) can be a strong daytime LSTID source, driven by particle precipitation (Nishimura et al., 2020). In the distribution presented in Figure 4, an increase is observed in the pre-noon sector, followed by a minimum around noon. This is consistent with the neutral wind-based storm scenario, which drives LSTIDs from auroral to lower-latitude regions. As long as auroral energy input is present, positive and negative effects in the ionosphere can be generated and drive TIDs (e.g., Prölss, 1993; Fuller-Rowell et al., 1994). Paznukhov et al. (2020) studied the local time dependence of the ionospheric storm effects, showing distinctly different local time patterns for negative and positive ionospheric storms. These results suggest that the composition disturbance zone responsible for the negative ionospheric storms tends to form around the sector of 0300–0500 LT and that the positive ionospheric storms prefer to begin around local noontime. Considering the latency time needed to detect an LSTID by the HF-INT (at least the time equivalent to half the period of the detected LSTID), the results of Figure 4 are consistent with this scenario.

Figure 4 Distribution of the occurrence of the SEC over the hours of the day (in UT) for Ebro Digisonde station data (January 1, 2022–June 30, 2023), illustrating the counts where SEC_Ebro > 50% (marked with blue) and the counts where SEC_Ebro>70% (orange). In this plot, “no data” values are disregarded.

The same type of distribution for the TEC gradient index is shown in Figure 5 for data collected between January 1, 2022 and June 30, 2023. Here, we consider the TEC gradient index that exceeds the 1.2 mm/km threshold (blue color) and the superimposed histogram of TEC gradient index that exceeds the 2 mm/km threshold (orange color). At the same time, “no data” is disregarded. The 1.2 mm/km threshold corresponds to moderate to intense activity, while the 2 mm/km threshold corresponds to intense activity. Considering that TEC is not affected by ionosphere irregularities, such as the sporadic Es layer, the strength of TEC gradient index peaks at noontime and decreases at night. TEC gradient intensification is more detectable during daytime thanks to the following mechanisms: i) the higher ambient TEC from solar ionization which allows sharper spatial gradients to develop and be measured, ii) the cusp particle precipitation during local noon which often produces abrupt TEC increases and strong localized gradients, iii) the ionospheric trough boundary near the equatorward edge of the auroral oval which is typically more pronounced in the afternoon sector, creating steep TEC drop-offs.

Figure 5 The number of TEC grad index computed at each hour of the day (in UT) between January 1, 2022, and June 30, 2023, for TEC grad exceeding the threshold of 1.2 mm/km (low activity) and the superimposed histogram for TEC grad exceeding the threshold of 2.0 mm/km (moderate activity).

The distribution of the IL and IU indices per hour of the day is shown in Figure 6. The IU index corresponds to the intensity of the eastward auroral electrojet and reaches its peak value at 1600UT, which is local afternoon hours in the European sector. The IL index peaks around midnight, where the westward electrojet intensity becomes stronger. Superposing the strength of the two electrojet indices, it is concluded that the intensity of the electrojets and consequently the generation mechanism of LSTID is in place at any time of the day. This agrees with the conclusion extracted from the behavior of the TEC gradient index distribution presented in Figure 5.

Figure 6 The mean values of the IU index (Left) and the IL index (Right) by hour of the day (in UT).

In summary, considering the complex landscape, we should also account for the time lag in the propagation of LSTIDs from higher to middle latitudes, which depends on local time. According to a statistical study performed by Tsugawa et al. (2004), longer delays are found during local daytime and faster, steeper wavefronts at night. According to Nykiel et al. (2024), the longer delays in the daytime sector are due to slower atmospheric wave propagation from the cusp. This complex behavior justifies the choice made in this contribution to approach the LSTID forecasting classification problem with a ML technology.

3 Methods

In data-driven ML, the idea that historical values are primary predictors when building time series forecasting models is widely supported (Wang et al., 2017; Montero-Manso & Hyndman 2021). From a signal processing perspective, this is also a common practice (e.g., Hua et al., 2019; Lim & Zohren, 2021). This approach is essential because the patterns and trends observed in the past data are used to make accurate predictions about future values.

In this contribution, the forecast of the occurrence of LSTIDs is considered a binary classification problem, with a pre-defined threshold theta of SEC, which is computed with the HF-INT method. More specifically, a value of SEC greater than theta at a specific timestamp, indicates the occurrence of a LSTID (corresponding to class 1), while a SEC value less than theta indicates no LSTID (corresponding to class 0). This section outlines the methodological framework employed to formulate and solve the binary classification problem for the LSTID forecast.

First, significant effort has focused on preparing the time series of data for analysis through pre-processing. As mentioned in the previous paragraphs, our approach leverages time series from three distinct data sources: IU-IL indices, TEC gradients inferred from TEC maps covering Europe, and HF-INT SEC computations over Digisonde stations. The pre-processing steps in our methodology include data cleaning, data integration, and data transformation. To achieve this, we implemented a series of steps to ensure the data’s consistency and reliability.

To ensure consistency, a standard format is adopted for all the time series. All data points were indexed chronologically based on their timestamps, and missing values were then addressed using interpolation techniques. Specifically, the time series data were resampled at 5-minute intervals. Any missing values resulting from upsampling were filled using the nearest valid measurement. TEC grad indices have undergone a temporal adjustment of 30 s to achieve alignment of timestamps at a minute level. To enrich the data for the classifier, the time of day and month of year is extracted as additional features. Finally, each data timeseries is normalized so that all its elements lie in the range from zero to one.

Following the data preparation process, the prepared data are fed into three different classifiers: two traditional ones, namely a k-Nearest Neighbors (k-NN) classifier, and a Feedforward Neural Network (FNN) (see e.g., Bishop, 2006; Hastie et al., 2009; Theodoridis & Koutroumbas, 2009), and a more advanced one, namely, the Temporal Fusion Transformer (TFT) (Lim et al., 2021).

k-NN is a fundamental and straightforward algorithm used in ML for classification tasks. Unlike complex models that require extensive training, k-NN leverages a simpler approach, which requires no learning at all and is based on the idea that similar data points tend to be categorized into the same class. More specifically, the classifier uses a reference set of vectors whose class is known. When presented with a new, unseen data vector, k-NN identifies its k nearest neighbors in the reference set based on a chosen distance metric (like Euclidean distance). And then, it assigns the data point to the class where most of its k nearest neighbors in the reference set belong. For the binary classification task at hand, the k-NN classifier’s input is solely the SEC time series, and its output is a binary value indicating the presence or absence of LSTID occurrences.

FNNs are one of the fundamental types of artificial neural networks. They are based on a simple structure processing element (usually called a neuron), which takes several input values, weighs them through a set of weights, computes the weighted sum (plus a bias term), and passes the summation through a nonlinear function. This result is the output of the neuron. Each neuron is parameterized by its weights and bias. The architecture of an FNN is composed of a set of neurons arranged in layers. Specifically, it consists of an input layer of nodes (which perform no processing), followed by one or more layers of neurons. The last of these layers is called the output layer (the output of the neurons in this layer constitutes the output of the FNN), while the previous ones are called hidden layers. The whole network is defined by the parameters and biases of the neurons that constitute it. In the classification framework, the training of an FNN (i.e., the learning of the input-output mapping of vectors to their correct classes) is carried out through the estimation of its parameters using a given set of data vectors (called the training set) for which the class is known. After training, the performance of the FNN is assessed using a separate set of vectors that were not used for training. If the performance is satisfactory, the FNN can be used operationally to classify other unknown vectors. FNNs are versatile tools, finding applications in various domains such as image recognition and speech classification. Their strength lies in their ability to learn complex relationships between inputs and outputs, making them a popular choice for numerous ML tasks. To train the FNN in the framework of the present application, time series data from IU-IL indices, TEC grad indices over Europe, and HF-INT SEC calculations from Digisonde stations were used.

The output of the FNN model is a binary value indicating LSTID occurrence (1) or not (0).

The Temporal Fusion Transformer (TFT), a novel approach for time series forecasting, was adapted – in our case – for binary classification of LSTID occurrence. Unlike traditional forecasting models, TFTs prioritize interpretability, allowing users to understand the reasoning behind their predictions (Lim et al., 2021). In general, TFT is most appropriate for time series forecasting, especially when dealing with complex, multivariable, and temporally dependent data. TFT handles sequential/temporal data, uses attention for long-term context, is designed for forecasting multiple future time steps simultaneously, includes built-in interpretable mechanisms, handles different types of inputs (e.g., calendar events) and dynamically weights inputs at each time step. In this application, the TFT was trained using a combination of inputs: time series data from IU-IL indices, the TEC grad indices across Europe, and HF-INT SEC calculations from Digisonde stations. By combining recurrent layers (capturing short-term trends) and self-attention mechanisms (learning long-term patterns), the TFT effectively processed these diverse inputs. This architecture enables multi-horizon forecasting, which, in this context, translates to understanding the temporal evolution of LSTID occurrence. The model also offers interpretability by highlighting the influence of specific input features, aiding in the identification of key factors contributing to LSTID events. Ultimately, the TFT, leveraging LSTM Sequence-to-Sequence and Transformer’s Self-Attention mechanisms, enriched its temporal representation with static contextual information, providing accurate predictions alongside valuable insights into the underlying dynamics of the phenomena. The model output is a binary value indicating LSTID occurrence (1) or not (0).

4 Results

For the needs of the experiments performed to find the best forecasting model, time series of HF-INT computed SEC over four Digisonde stations are collected: Juliusruh, Dourbes, Sopron, and Ebro. These datasets are combined with IU-IL indices and TEC grad indices over Europe, and have been preprocessed, as described in the previous section. The goal is to forecast the occurrence of LSTIDs. To achieve this, it is assumed that an LSTID occurs when the SEC is over a predefined threshold, theta, which is set to 50% or to 70%. The forecasting horizon varies from 5 min to 2 h ahead. Three different classifiers are implemented: a k-NN, an FNN, and a TFT classifier. Classification accuracy is measured in terms of the F1-score, which takes values between 0 and 1, with values close to 0 indicating degraded performance and values close to 1 indicating superior performance. F1-score is defined as the harmonic mean of two performance indices (usually competing with each other), namely, the Recall and the Precision. Recall is the fraction of real LSTID disturbance cases that are correctly classified by the classifier, while Precision is the fraction of cases classified as LSTID disturbances and that are indeed LSTID disturbances.

The F1-score curves versus the forecast horizon for the three classifiers are plotted in Figure 7. These curves show the performance of the classifiers on SEC time series from the Digisonde stations Juliusruh, Dourbes, Sopron, and Ebro, respectively, as the forecast horizon increases from 5 min to 120 min, for theta = 50% and theta = 70%. The F1-score is determined by averaging the scores across all non-missing samples within the two-month test period, which extends from May 1 to June 30, 2023.

Figure 7 The performance of the proposed TFT classifier compared to state-of-the-art tools, such as a FNN (a three-layer network with 30, 10, and 1 nodes in the 1st, 2nd, and 3rd layers, respectively) and a k-NN classifier (for k = 25). In our experiment, LSTID forecasts are performed for a time period up to two hours ahead for the Digisonde locations in Juliusruh (Germany, JR055), Dourbes (Belgium, DB049), Ebro (Spain, EB040), and Sopron (Hungary, SO148). The performance is measured in terms of the F1-score for the cases where theta takes the values 50% and 70%.

The forecasting accuracy for two hours ahead decreases from 0.95 to 0.65 approximately with increasing forecasting horizon, for theta = 70%, while the success depends on the chosen neural network, the combination of input features, and the geographic latitude of the Digisonde.

5 Discussion

Figure 8 presents the TFT classifier results obtained for an ionospheric storm that occurred during the test period, using training data from January 1, 2022, to June 30, 2023. The forecasted SEC and the calculated SEC are overplotted for three Digisonde locations, in Juliusruh, Dourbes, and Ebro.

Figure 8 TFT classifier results obtained for an ionospheric storm that occurred within the test period. The forecasted SEC and the calculated SEC are overplotted at three Digisonde locations, in Juliusruh, Dourbes, and Ebro. Grey horizontal lines in the last three panels correspond to theta = 70%. Predictions are generated with a 40-min forecast horizon.

The TFT classifier performance for Juliusruh and Dourbes locations is quite successful, considering it predicts all the disturbances with SEC > 70%. For the Ebro location, the TFT shows several false positive predictions, while two peak values at the beginning and at the end of this time series are correctly forecasted. The large number of “no data” points, especially for the Ebro station, might have an influence on the reliability of the model’s performance.

Several experiments are performed for the following two distinct scenarios: (a) values of SEC greater than 50% indicate moderate and strong LSTID activity, and (b) values of SEC greater than 70% indicate strong LSTID activity. The performance is assessed through the F1-score metric.

The confusion matrices² for the whole test period (May and June 2023), the value are given in Figure 9, for two values of theta (50% and 70%) and for two forecast windows (5 min and 40 min), and for model execution with the data obtained at the Juliusruh Digisonde location.

Figure 9 Confusion matrices with the TFT technology for the Juliusruh Digisonde station.

Regarding the results obtained for the Juliusruh location, with a theta value of 70% and a 5-min forecast horizon, the model achieved a high degree of accuracy, correctly identifying 964 out of 1037 positive cases and 920 out of 946 negative cases. When the forecast horizon increased from 5 min to 40 min, the mode’s performance deteriorated only for the positive class. It misclassified 382 out of 933 positive cases and only 80 out of 946 negative cases. Overall, for the easiest case (a theta value of 50% and a 5-min forecast horizon), the performance is close to excellent, indicating the possibility of using this result to issue a warning. For the most difficult forecast scenario (a theta value of 70% and a 40-min forecast horizon), the model performs with a relatively high number of false alarms, which calls for some improvements. Such improvements could consider the quality of input data and the use of additional features (see also the end of this section). However, it should be noted that the recall measure about the presence of LSTIDs remains high in all scenarios. This is very important in operational forecasting, where the cost of a missed event is usually far greater than the cost of a false alarm.

Regarding the results obtained for the Ebro location, the same type of performance is noted from the analysis of the confusion matrices presented in Figure 10. Considering the latitudinal difference between Ebro and Juliusruh, the similarity of results shown in the corresponding confusion matrices underlines the robustness of the model regarding the choice of the specific classifier and of the specific features.

Figure 10 Confusion matrices with TFT technology for the Ebro Digisonde station.

The unique characteristics inherent in the calculation of the SEC time series preclude a comparison of the proposed model’s performance to some of the standard benchmark models such as the persistent estimator. The SEC, derived from the HT-INT method, does not solely depend on data from a single sensor but also requires specific conditions to be met across multiple sensors to confirm the detection of LSTIDs. These conditions include specific period ranges (37–150 min), high confidence levels (above 0.975) for the dominant period, coherence across at least three neighboring sensors (differing by less than ±30%), a cross-correlation greater than 0.5 among these sensors’ time series, and an estimated velocity below 1800 m/s. Additionally, LSTIDs with SEC less than 10% are discarded. If any of these requirements are not met, the HF-INT method returns no estimates for all parameters, including SEC. Consequently, the detection of LSTIDs and the corresponding SEC can change significantly over short, overlapping time intervals, as a sudden fulfillment or non-fulfillment of these criteria can transition from a “nothing detected” state to a state with significant SEC percentage. Due to this inherent discontinuity and the multi-sensor dependency in SEC calculation, direct comparison with a persistent estimator – which assumes a smooth and predictable evolution of values – lacks validity. Although the persistent estimator serves as a baseline benchmark in many modeling contexts, its application here would misrepresent the dynamics and complexity of the SEC dataset.

To further assess the TFT model’s performance, the importance of the features’ contribution is provided in Table 3. For the highest latitude stations, the SEC has a dominant role. On the contrary, at the Ebro location, which is at lower latitudes, the importance of the SEC time series to the model results is inferior. The difference in the importance of features to the model performance might be connected to the relative location of the three locations with respect to the auroral oval and the trough. In other words, this difference might be the result of the decreasing LSTIDs amplitude as they travel equatorward and could highlight the dissipative nature of the ionospheric medium, as reported by Cherniak and Zakharenkova (2018) and Ding et al. (2012). This result provides evidence for future improvements with the consideration of additional features in the deep learning models, such as the propagation pattern, the damping effect, and the interhemispheric propagation, in addition to the drivers considered here.

The above methodology can also be considered in an adaptive framework, using as a performance criterion the F1-score. When the value of the F1-score decreases over time, retraining is required using the most recent data. The retraining can be executed in a single CPU computer and takes approximately 2 h; however, the use of a GPU could accelerate the training up to 50 times. The results of the retrained model using data collected up to June 30, 2024, are presented in Figure 11, for Dourbes Digisonde for July 30, 2024. In this case, a substorm is recorded in the first morning hours of the day, and the TEC grad index shows an increase with a delay of almost 3 h after the substorm onset. At the lower latitudes, the ionospheric observed characteristics from the Digisonde in Dourbes indicate high SEC in the evening of the same day. However, the “no SEC calculation” result is possibly due to the absence of data from at least three sensors that exhibit coherent periodicities and does not necessarily imply the absence of LSTIDs. The indication provided from the three drivers, i.e., IU, IL, and TEC grad indices, forms the condition for LSTID generation during daylight hours as a result of the auroral currents intensification in the early morning hours. This is successfully calculated with the TFT classifier with a 40-min forecasting horizon, as indicated in Figure 11.

Figure 11 The results of the forecast of the occurrence of LSTID with 40 min horizon, with the TFT classifier (bottom panel red line) obtained for an SEC threshold 70% on July 30, 2024, provided together with the SEC calculated with the HF-INT method. The TEC grad index and the IU and IL indices from the IMAGE magnetometer network, used by the TFT classifier as drivers of the LSTID activity, are also provided in the top and middle panel.

Several strategies can help to detect and overcome data or concept drift in the proposed SEC-based model:

• Performance monitoring: Regularly evaluate the TFT classifier’s performance on unseen data (holdout set) or live data (production). This helps to identify drops in accuracy, precision, recall, or other relevant metrics that might signal model drift.

• Data and concept drift detection tools: There are tests specifically designed for detecting data drift, such as the Page-Hinkley test or the Kolmogorov-Smirnov test, as well as algorithms that detect concept drift, such as the adaptive windowing (ADWIN) algorithm (Bifet & Gavalda, 2017). These tools can automatically identify changes in the input data or task that may indicate model drift.

• Data and concept drift prevention techniques: these techniques can help build ML models that are more robust against changes in the underlying data. For instance, using data augmentation or creating synthetic data can expose the model to a broader and more diverse set of examples. This wider exposure is expected to make the model more adaptable to shifts in real-world data patterns. Similarly, techniques like transfer learning and multitask learning can help the model adjust to evolving tasks or goals.

• Retraining and fine-tuning: If data drift is detected, retraining or fine-tuning the model on new data can help to overcome it. This can be done periodically or in response to significant changes in the data or task.

6 Conclusions

In this contribution, the forecast of the occurrence of LSTIDs is considered as a binary classification problem, for a pre-defined threshold theta of SEC, which is computed with the HF-INT method. Three different features are used. Specifically, the IU-IL indices time series and the TEC gradient Activity index time series are used as precursors of the LSTIDs, while the HF-INT SEC computations are used as a binary indicator for LSTID occurrence (1) or no LSTID occurrence (0). The SEC computations provide information for detected LSTIDs over specific Digisonde locations.

Three different classifiers are tested: the more advanced TFT classifier, which is the proposed forecasting model, and two traditional ones, (k-NN and FNN). Several experiments are performed for the following two distinct scenarios: (a) values of SEC greater than 50% indicate moderate and strong LSTID activity, and (b) values of SEC greater than 70% indicate strong LSTID activity. The performance is assessed through the F1-score metric. The forecasting accuracy decreases from 0.9 to 0.6 approximately with increasing forecasting horizon up to 2 h ahead, for the proposed TFT model, while FNNs have the next best performance and k-NN has the most inferior performance. The performance also depends on the geographic latitude of the Digisonde and especially its location with respect to the auroral surge intensification. The missing data is also an important issue, especially for the SEC computations that are performed by the HF-INT method in real time and the results depend on the availability and quality of the ionograms automatic scaling and on physical mechanisms that disturb considerably the identifications of the various ionospheric layers, such as sporadic E layers and negative ionospheric storm effects.

Even though, for a very short forecasting horizon (e.g., 5 min ahead), TFT may have worse performance compared to simpler models if the available dataset is small or lacks diversity (i.e., too many gaps), TFT exhibits a very satisfactory performance for longer forecasting horizon, which is the requirement for the operational application of the model. To improve the performance of the TFT model, a much longer series of data should be used for its training. In addition, noisy data should also be eliminated to avoid overfitting or misinterpretation, since TFT’s adaptability might lead it to adjust to minor fluctuations that are not indicative of significant changes, potentially reducing its forecast accuracy in stable periods.

Regarding future improvements, additional features may be tested, such as the time derivative of the auroral electrojet indicators and the SEC detected at stations of higher latitude to forecast the LSTID occurrence at lower latitudes.

The dependence of the TFT importance of the feature variables from the geographic location implies a future need for major restructuring of the LSTID forecasting framework, considering the propagation pattern, the ionospheric damping and the interhemispheric circulation.

The model currently operates offline. The source code is available in Zenodo (Themelis, 2023) with open access, and it is available for use and reuse. Development and exploitation plans include the operational application of the model in the PITHIA-NRF eScience Centre (Belehaki et al., 2025), the design of high-level data products with the forecasted SEC over the European Digisonde stations with a minimum forecast horizon of two hours, and the real-time provision of these products following the standards of the Space Weather Network of the European Space Agency.

Acknowledgments

The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Funding

This research has received funding from the Horizon Europe T-FORS project (Grant Agreement 101081835) of the European Commission. Partial support is received by the project PITHIA-NRF, funded by the Horizon 2020 Programme of the European Commission (Grant Agreement 101007599) and by the SWESNET/EIS project funded by the European Space Agency.

Data availability statement

Auroral electrojet indices used in this study were obtained by the FMI-managed IMAGE Network, publicly available at https://space.fmi.fi/image/www/index.php?page=real_time. The SEC time series and TEC grad indices were retrieved through the TechTIDE service of the ESA Space Weather Network (https://swe.ssa.esa.int/techtide). The specific preprocessed data required to reproduce the above findings are available to download from the GitHub repository: https://github.com/themelis/LSTIDs_forecasting. The Python code and its documentation supporting this research are available at Themelis K. 2023. LSTIDs Forecasting with the Temporal Fusion Transformer. Zenodo. https://doi.org/10.5281/zenodo.10442655.

References

All Tables

All Figures

Figure 1 On the top panel, MUF (MUF(3000)F2) data series from the 24 previous hours, from April 22 at 1800UT to April 23 at 1800UT. The red curve corresponds to the fit applied to remove the main daily trends for the last 6 hours, from April 23 at 1200UT to April 23 at 1800UT. The second panel shows the residuals, ΔMUF, for the last 6  h, the difference between the original data series and the fit. The third panel shows the periodogram of the last 6 h, with a dominant period around 130 min. Painted in magenta, the area defined by 0.8 TAM and 1.4 TAM, where TAM is the maximum amplitude period, in this case around 130 min. The final value of SEC is obtained from the ratio between the area painted and the total area under the curve.

In the text

Figure 2 Time variations of the characteristics of the LSTID detected by the HF-INT at Dourbes, Belgium, for April 23, 2023. The top plot depicts the activity levels of LSTIDs (green, yellow, orange, red, and magenta for Insignificant, Weak, Moderate, Strong, and Very Strong activity, respectively). The middle plot depicts the dominant period (black diamonds) and SEC (green diamonds), and bottom plots depict the magnitude of the velocity (red diamonds) and azimuth direction (blue squares) of propagation of the detected LSTID.

In the text

	Figure 3 The response of the LSTID drivers (IL and IU indices) and precursors (TEC Grad index) during the geomagnetic storm that occurred from 26 to 28 February 2023. The SEC calculations, using ionospheric characteristics recorded at Juliusruh, Dourbes, and Sopron Digisonde stations, are presented in the bottom panels.
In the text

	Figure 4 Distribution of the occurrence of the SEC over the hours of the day (in UT) for Ebro Digisonde station data (January 1, 2022–June 30, 2023), illustrating the counts where SEC_Ebro > 50% (marked with blue) and the counts where SEC_Ebro>70% (orange). In this plot, “no data” values are disregarded.
In the text

	Figure 5 The number of TEC grad index computed at each hour of the day (in UT) between January 1, 2022, and June 30, 2023, for TEC grad exceeding the threshold of 1.2 mm/km (low activity) and the superimposed histogram for TEC grad exceeding the threshold of 2.0 mm/km (moderate activity).
In the text

	Figure 6 The mean values of the IU index (Left) and the IL index (Right) by hour of the day (in UT).
In the text

Figure 7 The performance of the proposed TFT classifier compared to state-of-the-art tools, such as a FNN (a three-layer network with 30, 10, and 1 nodes in the 1st, 2nd, and 3rd layers, respectively) and a k-NN classifier (for k = 25). In our experiment, LSTID forecasts are performed for a time period up to two hours ahead for the Digisonde locations in Juliusruh (Germany, JR055), Dourbes (Belgium, DB049), Ebro (Spain, EB040), and Sopron (Hungary, SO148). The performance is measured in terms of the F1-score for the cases where theta takes the values 50% and 70%.

In the text

	Figure 8 TFT classifier results obtained for an ionospheric storm that occurred within the test period. The forecasted SEC and the calculated SEC are overplotted at three Digisonde locations, in Juliusruh, Dourbes, and Ebro. Grey horizontal lines in the last three panels correspond to theta = 70%. Predictions are generated with a 40-min forecast horizon.
In the text

	Figure 9 Confusion matrices with the TFT technology for the Juliusruh Digisonde station.
In the text

	Figure 10 Confusion matrices with TFT technology for the Ebro Digisonde station.
In the text

Figure 11 The results of the forecast of the occurrence of LSTID with 40 min horizon, with the TFT classifier (bottom panel red line) obtained for an SEC threshold 70% on July 30, 2024, provided together with the SEC calculated with the HF-INT method. The TEC grad index and the IU and IL indices from the IMAGE magnetometer network, used by the TFT classifier as drivers of the LSTID activity, are also provided in the top and middle panel.

In the text