© A.A. Ferreira et al., Published by EDP Sciences 2025

1 Introduction

Geomagnetic storms are closely related to the occurrence of disturbances in the thermosphere and ionosphere. The list of disturbances associated with geomagnetic storms includes for example, changes in the neutral and electron densities (Prölss, 1980; Mendillo, 2006; Codrescu et al., 2022), enhancement of ionospheric conductivities and currents (Buonsanto, 1999), and variations in the neutral winds (Fuller-Rowell & Codrescu, 1996; Jonah et al., 2020). During such events, a considerable amount of energy and momentum can be deposited in the auroral zone via enhanced electric fields, precipitation of particles, and Joule heating (Aa et al., 2019). The deposition of energy and momentum contributes to the triggering of wave-like perturbations in the ionospheric electron density that exhibit an apparent motion with horizontal velocities ranging from 400 up to 1000 m/s, periods varying from 30 min to 3 h and wavelengths greater than 1000 km. These perturbations are called Large Scale Travelling Ionospheric Disturbances (LSTIDs) and are the manifestation in the ionosphere of Atmospheric Gravity Waves (AGWs) generated due to auroral and geomagnetic activity (Hines, 1960; Hunsucker, 1982; Cherniak & Zakharenkova, 2018). For decades, TIDs and their theoretical understanding have been the subject of several investigations (Munro, 1948; Hines, 1960; Francis, 1975; Hocke & Schlegel, 1996) and different studies have been conducted in order to model, detect, and statistically characterize their occurrence using different instruments. Such investigations have included, for example, the use of ionosondes (Morgan et al., 1978; Maruyama et al., 2004; Reinisch et al., 2018), All-Sky imagers (Kubota et al., 2000; Otsuka et al., 2004) and Global Navigation Satellite Systems (GNSS) data (Saito et al., 1998; Tsugawa & Saito, 2004; Hernández-Pajares et al., 2006; Otsuka et al., 2013; Zakharenkova et al., 2016). In addition, significant efforts towards the modelling/prediction of TIDs using semi-empirical modelling approaches and general circulation models have been made (e.g., Mayr et al., 1990; Millward et al., 1993; Ridley et al., 2006; Fedorenko et al., 2013; Yokoyama, 2014; Sheng et al., 2020; Vadas et al., 2023; Bukowski et al., 2024). Despite the considerable progress in modelling/understanding TIDs, predicting them is still a challenging task due to the complexity of the coupling mechanisms in the thermosphere-ionosphere system. Recently, a new methodology for predicting TID activity over mid-latitude Europe during geomagnetic storms has been presented by Borries et al. (2023). This methodology is based on the use of a linear regression approach in order to predict the TID activity index (A_TID). The index is derived from GNSS Total Electron Content (TEC) measurements and is also well-suited for statistical analysis and investigation of the source mechanisms of LSTIDs. The model uses solar wind data derived from Lagrangian point L₁ as input, more specifically the Kan–Lee merging electric field. Results presented by Borries et al. (2023) have shown that the proposed approach outperforms the persistence model (PM) for predictions beyond 30 min for strong levels of TID activity (i.e. A_TID ≥ 0.5 TECU).

The methodology proposed by Borries et al. (2023) is based on a linear regression approach and it has been proven to have a good performance in predicting LSTIDs activities at mid-latitudes. Different prediction methodologies have been developed throughout the years and applied in Space Weather research. Such methodologies may include, for example, neural networks (NNS, Haykin, 2009; Camporeale, 2019) and multi-model ensembles (Murray, 2018). In this context, we investigate to what extent state-of-the-art methodologies can improve the prediction of LSTID activity at mid-latitudes. For this purpose, we develop additional models for predicting the level of LSTID activity using NNs and multi-model ensemble approaches. The performance of the methodologies is evaluated in a set of different geomagnetic storms (see Appendix) and the geomagnetic storm observed on November 20, 2003 is used to illustrate how the methodology works.

This work is structured as follows: Section 2 presents the data used in this investigation and the methodologies proposed for predicting the LSTID activity. The results are presented in Section 3. The discussion of the results and the conclusions are presented in Sections 4 and 5, respectively.

2 Data and methodology

2.1 TID activity index (A _TID )

As previously mentioned, TIDs can be detected from GNSS data. The slant TEC (sTEC) is obtained by combining carrier phase/code measurements at two frequencies f₁ and f₂, according to equations (1) and (2):$$\mathrm{sTEC}=\frac{f_1^2{f}_2^2}{40.3(f_1^2-f_2^2)}\left[({\mathrm{\Phi }}_1-{\mathrm{\Phi }}_2)+B_{\mathrm{amb}}+{\epsilon }_{{\mathrm{\Phi }}_1-{\mathrm{\Phi }}_2}\right],$$(1) $$\mathrm{sTEC}=\frac{f_1^2{f}_2^2}{40.3(f_1^2-f_2^2)}\left[\left({\mathrm{\Psi }}_2-{\mathrm{\Psi }}_1\right)+{\epsilon }_{{\mathrm{\Psi }}_2-{\mathrm{\Psi }}_1}\right],$$(2)where Φ and Ψ correspond to carrier and code pseudoranges, respectively. Noises (e.g. thermal noise) are represented by ${\epsilon }_{{\mathrm{\Phi }}_1-{\mathrm{\Phi }}_2}$ and ${\epsilon }_{{\mathrm{\Psi }}_1-{\mathrm{\Psi }}_2}$. B_amb is the carrier-phase ambiguity. Additional terms, such as inter-frequency biases and multipath are not included for the sake of simplicity (Hoque & Jakowski, 2012).

For the TID detection, the slant TEC at an elevation angle ε is then converted to an equivalent vertical TEC (vTEC) in a thin-shell model of the ionosphere according to equation (3):$$\mathrm{vTEC}=\sqrt{1-{\left(\frac{R_e\cos \epsilon }{R_e+h_i}\right)}^2}\times \mathrm{sTEC},$$(3)where R_e is the Earth’s radius, h _i is the height of the thin-shell model, assumed here to be 350 km. After removing the regular trend of the background vTEC of the ionosphere, one can obtain a perturbation TEC for each satellite-receiver link, which allows the identification of the TIDs. The removal of the regular trend of the vTEC can be done in different ways (e.g. polynomial fitting, moving average or other trend estimation approaches), with no significant difference between the methodologies (Borries et al., 2023).

Recently, a new TID index was proposed by Borries et al. (2023), which is obtained as an extension of the already established ways of detecting TIDs from GNSS measurements. The index can be used to monitor the TIDs activity and also to perform statistical analysis and prediction of such disturbances. In this case, the TEC variability obtained by a bandpass filtering (TEC_bp) in the range of 30–60 min is extracted by calculating the difference of moving averages with 60 and 30 min window size, as presented in equation (4) $${\mathrm{TEC}}_{\mathrm{bp}}\left(t\right)=\frac{1}{T_{30}}\sum _{n=t-0.5T_{30}}^{t+0.5T_{30}}\mathrm{}\mathrm{v}\mathrm{TEC}\left(n\right)-\frac{1}{T_{60}}\sum _{m=t-0.5T_{60}}^{t+0.5T_{60}}\mathrm{}\mathrm{v}\mathrm{TEC}\left(m\right),$$(4)where vTEC corresponds to the vertical TEC; T₃₀ = 30 min/Δt and T₆₀ = 60 min/Δt, with Δt corresponding to the sampling time. For each time step t, the TID activity index is obtained by taking the difference between the maximum TEC_bp and the minimum TEC_bp in a 60-minute window size centred on t and multiplying by 0.5, as presented in equation (5) $${\mathrm{A}}_{\mathrm{TID}}\left(t\right)=0.5\left(\max \right(\mathrm{TE}{\mathrm{C}}_{\mathrm{bp}}\left(x\right)\times F(x-t))-\min (\mathrm{TE}{\mathrm{C}}_{\mathrm{bp}}\left(x\right)\times F(x-t)\left)\right),$$(5)where t − (30 min/Δt) ≤ x ≤ t + (30 min/Δt) and F corresponds to a Gaussian function applied to reduce the impact of the values that are away from the centre of the 60-minute window.

The index can be used for individual GNSS stations, but can also be employed on the mapping of TIDs activity. It also allows the investigation of LSTIDs and MSTIDs. It is important to highlight that although a 30–60 minute bandpass filter is used to obtain the A_TID index, other disturbances may also be reflected in the A_TID, especially at high and low latitudes. The source of such disturbances cannot be distinguished from the index itself. However, at mid-latitudes, the impact of other disturbances on the A_TID index is expected to be smaller than at high and low latitudes.

Correlation studies with several solar and geomagnetic indices, such as the Kan–Lee merging electric field (E_KL) and Auroral Electrojet (AE) index, have shown that the magnitude of the LSTIDs at mid-latitudes is well correlated with them. Based on this information, in this work we propose different methodologies for predicting the TIDs activity using solar wind data from Lagrangian Point L1 as input. The description of each methodology is given in the following sections.

2.2 November 20, 2003 geomagnetic storm

In this investigation, we analyze the prediction of the TID activity index over a set of more than 50 geomagnetic storms (see Appendix). To illustrate the methodology, the geomagnetic storm that occurred on November 20, 2023 is used. Figure 1 presents the evolution of SYM-H, Bz, solar wind speed and AE index during the storm. This was a strong geomagnetic storm event, which was related to the occurrence of an Interplanetary Coronal Mass Ejection (ICME) (Zhang et al., 2007), which reached the Earth on the morning of November 20, 2003 and led to a SYM-H index of approximately −490 nT at around 18 UT. Figure 1b shows the Bz component of the interplanetary magnetic field, which increased to about 33 nT and turned southward for more than 12 h, reaching a minimum value of about −52 nT. A significant increase in the solar wind speed is observed at about 11:21 UT, reaching 766.3 km/s. A significant increase in the auroral electrojet index is also observed, reaching over 2000 nT during this period, which indicates a high level of auroral activity. This geomagnetic storm and the ionospheric and thermospheric responses associated with this event have been reported in several investigations (Meier et al., 2005; Crowley et al., 2006; Becker-Guedes et al., 2007; Borries et al., 2017).

Figure 1 (a) SYM-H index, (b) Bz component of the interplanetary magnetic field, (c) solar wind speed and (d) AE index for 19–21 November, 2003.

2.3 Persistence Model (PM)

The PM, although simple, is considered a useful prediction approach which assumes that the future value of a predicted variable is equal to the most recent observation (Reikard, 2018; Paulescu et al., 2021). In this case, the prediction of the TID activity A_TID in a horizon prediction n, in minutes, is assumed to be equal to the most recent observation, i.e. ${\widehat{\mathrm{A}}}_{\mathrm{TID}}(t+n)={\mathrm{A}}_{\mathrm{TID}}\left(t\right)$. This type of approach has been used for the prediction of different space weather parameters, such as solar radio flux, solar flare activity and the geomagnetic K-index (Devos et al., 2014). In this work, the PM is used as a baseline to evaluate the prediction performance of the proposed methodologies.

2.4 Linear Regression Models (LR)

The main goal of this analysis is to use the A_TID to forecast TID activity in mid-latitudes, based on solar wind-derived parameters. It is therefore necessary to know which geophysical parameters or indices reflect the driving mechanisms for the generation of such disturbances. This analysis has been performed by Borries et al. (2023) by computing the cross-correlation of the A_TID with different potential candidate parameters. The study separates two different types of geophysical parameters: geomagnetic indices and solar-wind magnetosphere coupling functions. In both cases, different parameters show a good correlation with the A_TID index.

Although the investigation presented by Borries et al. (2023) analyses the correlation with 15 geophysical parameters and the TID activity index, in this work only the AE index and the Kan–Lee merging electric field are going to be used in the predictions, due to their high correlation shown in the previous study. It is important to highlight that other indices presented comparable correlations, such as the intermediate function (E_WAV) and the modified version of the Akasofu index (ε₃), but for simplicity, only the AE index and E_KL are used here. In the LR models, the prediction of the TID activity is performed using a linear regression approach, with one of the aforementioned indices as input (AE index or E_KL). A description of each one of these cases is given below.

2.4.1 LR_EKL

Borries et al. (2023) proposed LR_EKL and presented good results in predicting the TIDs activity using data from Lagrangian point L1. The methodology is based on the Kan–Lee merging electric field (E_KL), which has been chosen due to its good performance in the study conducted in Newell et al. (2007), which investigated the geoeffectiveness of several parameters. In addition, the correlation studies conducted by Borries et al. (2023) have shown a correlation of 0.79 between the maximum A_TID and the maximum E_KL registered in the 18 h prior to the occurrence of the maximum A_TID for different geomagnetic storms. Equation (6) describes the E_KL coupling function$${\mathrm{E}}_{\mathrm{KL}}=v{B}_T\mathrm{si}{\mathrm{n}}^2({\theta }_c/2),$$(6)where v is the solar wind speed, $B_T=(B_y^2+B_z^2{)}^{1/2}$ is the transverse component of the Interplanetary Magnetic Field (IMF), with B _y and B _z corresponding to y and z components of the IMF and θ _c  = arctan(B _y /B _z ) is the IMF clock angle.

The investigation presented by Borries et al. (2023) was performed for a set of stations in the European–African sector to verify how the latitude impacts the correlation coefficients. It was observed that good correlations are observed mainly at mid-latitude stations and the prediction results for the GNSS station GLSV (50° N, 30° E) were presented. In our study, the GLSV station is also selected to investigate the TID activity at mid-latitudes.

In order to derive the linear regression model, the leave-one-out cross-validation is used. In this approach, one of the N geomagnetic storm events is used for validation and the other N − 1 remaining events are used to derive a linear regression model to be used for A_TID prediction on the validation event. This allows to use all the available geomagnetic storms for the validation of the proposed approach. In this study, for the GLSV station, 56 linear regression models are generated, each with another storm event kept for validation. To illustrate the approach, the geomagnetic storm registered on November 20, 2003 is selected. Figure 2 shows the linear regression fit line obtained for the prediction of the A_TID during the event.

Figure 2 Scatter plot of the maximum TID activity index of the GNSS ground station GLSV during each one of the investigated storm events versus the maximum Kan-Lee electric field EKL in the 18 h ahead of the maximum TID activity index (Borries et al., 2023).

For the period of 19–21 November 2003, which includes the main phase of the geomagnetic storm of November 20, 2003, for example, the following simple linear equation is obtained from the data plotted in Figure 2, y = 0.152 + 0.021x. This equation serves as the basis for estimating the A_TID for this event.

Since the correlation study presented in Borries et al. (2023) was carried out by analysing the maxima values, the maximum E_KL in a 2-hour interval (max_2h E_KL, Fig. 3b) is used to predict the A_TID index. Different intervals could be used, but we have empirically chosen 2 h, aiming at maximizing the prediction performance.

Figure 3 (a) EKL and (b) max2 h EKL for one day before, the day of the onset and one day after November 20, 2003 geomagnetic storm.

The following linear regression equation derived from Figure 2 is then used for the prediction:$${\widehat{\mathrm{A}}}_{\mathrm{TID}}(t+n)=0.152+0.021\times \mathrm{ma}{\mathrm{x}}_{2\mathrm{h}}{\mathrm{E}}_{\mathrm{KL}}\left(t\right),$$(7)where n represents the 0, 30, 60, … 180-minute prediction time. The ${\widehat{\mathrm{A}}}_{\mathrm{TID}}$ estimates are then compared with the A_TID reference values at t, t + 30, …, t + 180.

This procedure is performed for each storm event. Here we have used a leave-one-out cross-validation scheme to make the predictions. In this case, assuming that the maxima from N geomagnetic storm events are available, we use the maxima from N − 1 events to generate the scatter plot and then derive the linear regression equation, which is used to predict the A_TID for the Nth storm event following the procedure described before.

The results of this methodology for predicting the A_TID for the November 20, 2003 geomagnetic storm can be found in Section 3.1.

2.4.2 LR_NNAE

Previous investigations have shown that the Auroral Electrojet (AE) index presents a high correlation with the amplitude of LSTIDs (Borries et al., 2009, 2023). The recent study presented in Borries et al. (2023) has found a correlation coefficient of 0.72 between the maximum A_TID registered at each event, and the maximum AE index registered in the 18 h prior to the occurrence of the maximum A_TID. For our predictions purpose, however, instead of using the AE index, a predicted AE index is used. Although the World Data Center (WDC) for Geomagnetism Kyoto provides plots of the real-time AE index on its website, the real-time data is currently not available to the public, and the provisional data may take several days to be officially released (WDC, 2022). Therefore, to evaluate the prediction of the A_TID for an operational scenario, we use the AE index derived from solar wind measurements at Lagrangian Point L1 and predicted by an NN model. The NN predicted AE (NNAE) is based on the model described in Ferreira & Borges (2021), where the AE index is predicted using the solar wind velocity, the B _y and B _z components of the IMF together with the day of the year and temporal information. The idea of including the NNAE index is to verify if it can provide better performance than using the E_KL coupling function as input. It is important to note that the predictions of the A_TID obtained by using the NNAE may differ from those that would be obtained by using the AE index.

The prediction of the A_TID index is then made in a similar way to the procedure presented in Section 2.4.1. In this case, the max_2h NNAE is used as the input to the model, as shown in equation (8):$${\widehat{\mathrm{A}}}_{\mathrm{TID}}(t+n)=0.098+0.0002\times \mathrm{ma}{\mathrm{x}}_{2\mathrm{h}}\mathrm{NNAE}\left(t\right).$$(8)

Figure 4 presents the maximum values of A_TID and NNAE together with the linear regression equation, which is used to perform the A_TID index prediction for the November 20, 2003 geomagnetic storm.

Figure 4 Scatter plot of the maximum TID activity index of the GNSS ground station GLSV during each of the 56 storm events versus the maximum NNAE in the 18 h ahead of the maximum TID activity index (Borries et al., 2023).

The results of this methodology and their comparison with the other investigated methodologies are presented in Section 3.1.

2.5 Neural Network model (NN)

As an alternative to the PM and the LR models presented in the previous sections, a model based on artificial neural networks (NNs), more specifically, the multi-layer perceptron (Haykin, 2009) is proposed. Different investigations in space weather and ionospheric research have shown that multi-layer perceptron NNs can be a useful tool for prediction tasks (Orus-Perez, 2018; Camporeale, 2019). These networks are based on single units called neurons, which can be represented as$$a=g(W_{x}X+B),$$(9)where a is the output of the neuron; g is the activation function, which is usually a non-linear function such as hyperbolic tangent or rectified linear unit (ReLu); W _x is the synaptic weight of the neuron; X is the input value; and b is a bias term (Haykin, 2009; Orus-Perez, 2018). These neurons are massively interconnected and the network can be structured in various layers. The network can then be trained iteratively using input–output pairs.

The NN model proposed in this work was implemented based on the TensorFlow open source library¹ and takes as inputs the Solar Zenith Angle (SZA), to account for the daily variability observed in the LSTIDs activity described in Borries et al. (2009), and the historical values of A_TID and E _KL . An illustrative schematic of the NN used in this work to predict the A_TID is presented in Figure 5.

Figure 5 Schematic representation of the NN Model used on the ATID index prediction.

The NN model proposed in this work consists of one input layer, one output layer and two hidden layers of neurons. For the hidden layers, the ReLu activation function was used and a linear activation function was used in the output layer. The network was trained with 600 epochs, a learning rate of 1.5 × 10⁻⁴ and a batch size equal to 2048, and the training/validation procedure was performed using the leave-one-out cross-validation scheme in the same way as presented in Section 2.4. In this case, the performance of the model is tested for all available events, which differs from the procedure presented in Kim et al. (2021), which presents the prediction of the ionospheric state during geomagnetic storms using LSTM neural networks. In that study, from a set of 70 geomagnetic storms, 60 events are selected for training, 7 events for validation and 3 events for testing. Since in the proposed approach the model is tested in all events, it may provide a more comprehensive assessment of the model’s performance during different parts of the solar cycle and for storms with different characteristics.

The results of the NN model and the comparison with the other proposed methodologies are presented and discussed in the following sections.

3 Results

This section presents the evaluation results for each one of the proposed models, taking into consideration 56 geomagnetic storm events observed from 2001 to 2017. The storms from 2001 to 2007 were selected based on the list of geomagnetic storms presented in Borries et al. (2009). In addition, geomagnetic storm events from 2008 to 2017 have been included based on SYM-H index values below −50 nT. For each geomagnetic storm event, 3 days are usually included in the analysis: the day on which the SYM-H index reached the minimum value, and the days before and after. It is important to highlight that the list of selected events presented in Appendix does not include all geomagnetic storm events with SYM-H index values below −50 nT in the period from 2001 to 2017. Geomagnetic storm events with missing or significant gaps in GNSS data, AE index and solar wind measurements were not included in the analysis. Similarly to the investigation presented in Borries et al. (2023), the geomagnetic storm event observed on November 20, 2003 is used to illustrate the results of each model on the prediction activity.

Figures 6a and 6c show the results of the 1-hour prediction using the LR models. The A_TID computed for the GLSV station during this event reached magnitudes of about 0.8 TECU, indicating a moderate/strong TID activity over the station. During this period one can observe that both proposed LR based methodologies were able to clearly reproduce the increase of the A_TID. However, the rapid fluctuations in the A_TID are not well reproduced by the models.

Figure 6 Left panels: One-hour prediction of the TID activity index during the geomagnetic storm event registered on the 20th of November 2003 using: (a) LREKL, (c) LRAE, (e) PM and (g) NN model. Right panels: The same predictions presented in the left panel, but evaluated as a classification problem, showing True Positives (TP), True Negatives (TN), False Positives (FP) and False Negatives (FN) for a threshold of 0.5 TECU obtained using: (b) LREKL, (d) LRAE, (f) PM and (h) NN model.

Figure 6e shows the prediction results obtained from the PM and Figure 6g the results of the NN model described in Section 2.5. It can be observed that the NN model also clearly reproduces the A_TID increase. Similar to the results obtained with the LR_EKL and LR_AE models, one can observe that the NN model also overestimated the A_TID index in the early morning of November 21, 2023. The right column of Figure 6 presents the results of the model when it is evaluated from the perspective of a classification problem, which is described in the following section.

3.1 LSTID activity detection performance

As presented in Figure 6, all the proposed methodologies are able to predict fairly well the A_TID index increases observed during the storm. The rapid fluctuations observed in the A_TID index are, on the other hand, not very well depicted by the methodologies. Similar to what is presented in Borries et al. (2023), the models are evaluated from the perspective of a classification task. According to this approach, any predicted value of A_TID greater than or equal to the threshold level is considered as an LSTID activity event and a predicted value below this threshold is considered as a non-LSTID activity event. This type of classification leads to the possibilities of True Positive (TP, when both the reference and predicted values are above the threshold), True Negative (TN, when both the reference and predicted values are below the threshold), False Positive (FP, when the reference value is below the threshold, but the predicted value is above the threshold) and False Negative (FN, when the reference value is above the threshold, but the predicted value is below the threshold) by comparing the reference and estimated values (Borries et al., 2023). Figure 6 (right panels) shows the predictions of the proposed models, based on this classification methodology. Here the chosen threshold to identify the occurrence or not of LSTIDs is equal to 0.5 TECU, which corresponds to the 98% percentile of the A_TID considering all values and events. One can observe that all methodologies present good estimates of TNs and different performances in the identification of TPs. In this specific case, one can note that the LR_EKL present the best performance in the detection of TP values.

The metric used in the evaluation of the predictions is the True Skill Statistic (TSS), which is computed based on the difference between the True Positive Rate (TPR) and the False Positive Rate (FPR). The latter corresponds to the probability of false alarm and is computed as the ratio of False Positives to the number of negative events (N). The former indicates the ability to find all positive events and corresponds to the ratio of true positives to all positive events (P), as shown in equation (10).$$\mathrm{TSS}=\mathrm{TPR}-\mathrm{FPR}=\frac{\mathrm{TP}}{\mathrm{P}}-\frac{\mathrm{FP}}{\mathrm{N}}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}-\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}.$$(10)

The TSS has the advantage of being unbiased with respect to class imbalance (Detman & Joselyn, 1999; Camporeale, 2019; Borries et al., 2023). It ranges from −1 to 1, where −1 is interpreted as always wrong predictions (Bobra & Couvidat, 2015), 0 or less mean no better than random performance (Detman & Joselyn, 1999) and 1 means perfect forecasts. The TSS is obtained by considering the predictions of all the investigated events for the different models and different prediction times, which vary from 0 to 180 min. Figure 7a shows the TSS obtained by each model considering the different lead times for prediction.

Figure 7 (a) TSS for the different models under investigation, i.e. regression (obtained from EKL and AE), NN and the PM and (b) the TSS obtained from the combination of some of such models using a weighted average scheme.

As expected, for nowcast activities, the PM reaches the maximum achievable TSS of 1, since it is identical to the ground truth. For longer prediction times, one can observe that the TSS for the PM significantly drops when increasing the leading time of prediction, reaching about 0.15 for 180-minute predictions. The other models presented lower TSS for nowcasting activities, but one can notice that the reduction in performance with the increase in the leading times is not as dramatic as the one observed for the PM.

For lead times greater than 30 min, the TSS of LR_EKL is greater than the TSS of all the models under analysis. One can also observe that in this case the TSS peaks at 30 min (TSS = 0.8) and decreases continuously for longer lead times, reaching a TSS of 0.5 for 180 min predictions. The LR_NNAE model outperforms the PM for lead times equal to or greater than 60 min, but has a lower performance when compared to the LR_EKL and a similar performance when compared to the NN model. It is important to highlight that the performance of the LR_NNAE model is affected by the performance of the AE index prediction. The NN-based model shows an overall TSS higher than the PM and the LR_NNAE for prediction lead times equal to or greater than 60 min. One can observe from Figure 7a that all proposed models presented higher TSS when compared to the baseline PM for predictions beyond 60 min.

To further investigate the performance of each model, Figure 8 shows the TPR, FPR and the TSS for a one-hour prediction. From the results, one can clearly note that all the models present low FPR, with the NN and PM models presenting the lowest levels (below 1%). The main difference in performance is led by the TPR, which is the highest for the predictions obtained from the LR_EKL. In this case, the observed TSS was 0.76. The other proposed models showed similar performance among them.

Figure 8 Statistics obtained from the different models for 60 min prediction: (a) True Positive Rate, (b) False Positive Rate and (c) True Skill Statistic.

3.2 Multi-Model Ensemble (MME)

An additional methodology investigated in this work is the multi-model ensemble (MME), which corresponds to a combination of different models in order to increase the performance of a single-model prediction. The main focus here is to investigate whether such a methodology can improve the performance of the model with the highest performance obtained so far, i.e. LR_EKL. A basic approach is the linear combination of different methods in order to create a better ensemble forecast (Murray, 2018). Taking probability forecasts as an example, one can combine multiple forecast probabilities outputs from n different models according to P_ens = ∑ _n ω _n p _n , where p denotes a forecast probability value and ω is the weight value (often chosen in such a way that the sum of all weights equals to 1 (Murray, 2018; Zhou, 2012).

Different approaches can be used to set the weights for each model, which may go from simple ensemble averages (which tend to be used in operational forecasting methods, for example) up to nonlinear weighting schemes. Other weighting schemes based on performance metrics can be used in order to improve the forecasts based on the user-end requirements (Murray, 2018). In this work, a weighting scheme that takes into account the performance of each model in predicting the A_TID at different levels is adopted. The weights were heuristically chosen in this case. The weighting procedure is performed according to equation (11) $${\mathrm{A}}_{\mathrm{TID}}\left(t\right)=\{\begin{array}{ll}0.8{\mathrm{L}\mathrm{R}}_{\mathrm{EKL}}\left(t\right)+0.2{\mathrm{M}}_2\left(t\right),& \mathrm{if} \mathrm{L}{\mathrm{R}}_{\mathrm{EKL}}\left(t\right)\ge 0.5 \mathrm{and} M_2\left(t\right)\le 0.5\\ 0.2{\mathrm{L}\mathrm{R}}_{\mathrm{EKL}}\left(t\right)+0.8{\mathrm{M}}_2\left(t\right),& \mathrm{if} \mathrm{L}{\mathrm{R}}_{\mathrm{EKL}}\left(t\right)\le 0.5 \mathrm{and}{ M}_2\left(t\right)\ge 0.5\\ 0.5{\mathrm{L}\mathrm{R}}_{\mathrm{EKL}}\left(t\right)+0.5{\mathrm{M}}_2\left(t\right),& \mathrm{otherwise},\end{array}$$(11)where LR_EKL(t) is the prediction obtained from the LR_EKL model step t, and M₂ (t) is the prediction obtained from one of the other models under investigation (i.e. PM or NN model). The LR_EKL model was included in all multi-model combinations due to its good performance when compared to the other models, in an attempt to investigate if such combinations could improve the performance of the single model. Many other approaches for assigning the weights could be adopted, such as defining the weight according to the prediction time (i.e., higher weights for PM for short-term predictions and LR_EKL for long-term predictions). In this work, for the sake of simplicity, we only demonstrate the performance of the approach shown in equation (11).

The TSS for the ensembles investigated in this work are presented in Figure 7b. The results indicate that the combination of the LR_EKL and PM models can improve by about 22% and 8% for the nowcast and 30-minute predictions, respectively, when compared to the predictions obtained by the LR_EKL alone. For a prediction lead time of 60 min, no significant changes were observed. Beyond this prediction horizon, the results indicate that the combination of models did not lead to an improvement in performance, but rather to a slight degradation. The ensemble of LR_EKL with the NN model shows an improvement in the nowcasting performance, with the performance of the combination exceeding the individual performance of either methodology. For predictions equal to or longer than 30 min, no improvement was observed, but a slight decrease in performance.

4 Discussion

The analysis presented in the previous section shows different approaches for the prediction of LSTID activities, which were based on linear regression, neural networks and multi-model ensembles. The methodologies are based on the use of solar wind data from the Lagrangian point L1 to predict the level of TID activity. The model proposed by Borries et al. (2023), which is based on a linear regression with the E_KL as input, has already shown good performance in predicting the A_TID index. By considering the solar wind propagation time from L1 to the Earth’s bow shock, an immediate impact of the solar wind on the TID generation and the time taken by the LSTID to travel from the border of the auroral oval to the mid-latitudes, one can see a good agreement with the lead times that presented good performance, confirming the applicability of the models for predicting LSTID activity moments ahead of the storm commencement. In order to further advance and in an attempt to improve such predictions, a new linear regression-based model taking the NNAE index as input has been investigated. The results obtained have shown that the performance of the LR_NNAE has not surpassed the performance of the LR_EKL in any of the prediction horizons investigated, which suggests that, although it outperforms the baseline model (PM) for predictions beyond 60 min, the NNAE is not an optimal parameter for estimating the LSTIDs activity index via LR.

By analysing the performance of the NN model presented in Figure 7a, one can note that, except for nowcasting purposes, no improvement was observed when using a more complex model to predict the A_TID index. For predictions beyond 60 min, the NN model outperforms the PM model and has a similar performance to the LR_NNAE model, but the model presents a performance of 35% lower than LR_EKL on average. These results, although unexpected at first glance given the recent successes reported for the use of machine learning algorithms in various applications, represent a good opportunity to discuss the idiosyncrasies and challenges of the machine learning methods that have to be carefully considered when using such a methodology. Such features may include data quality, over/underfitting, network topology, training/validation, testing and metric methodologies (Camporeale, 2019). In this case study, although a larger number of geomagnetic storms were registered from 2001 to 2017, a limited amount of geomagnetic storms was used to develop and test the model. Several geomagnetic storm events were removed from the analysis, due to the occurrence of data gaps, either in solar wind data, geomagnetic data (AE index) or GNSS data, which limits the amount of representative data to be used to train the NN model. In addition, each one of the observed geomagnetic storm events has its own features in terms of time of occurrence, background ionization levels, and coupling from below, for example. As reported by Camporeale (2019), the too often too quiet problem makes the space weather datasets quite unbalanced, which poses a serious challenge to any machine learning algorithm trying to find patterns in the data. Such complexity summed up with a reduced dataset and high diversity of the geomagnetic events, poses a challenge for the prediction using the proposed NN model. Similar reduced performance of artificial neural network models when compared to simple methodologies has been already observed in other investigations, such as Reikard (2018), where logistic regression consistently outperformed the artificial neural network in predicting solar irradiance, sunspot number, and the Aa and Am geomagnetic indices. Also, in some cases presented in Wrench et al. (2022), NN models performed worse than simple models (e.g. linear interpolation) in interpolation tasks for predicting statistics of sparse IMF data series obtained from the Parker Solar Probe. These studies evidence that a gain in performance from using such a methodology may not be always the case. However, it is important to highlight that in this study the NN model outperformed the LR_EKL for nowcasting and outperformed both the LR_NNAE and the PM for the predictions beyond 60 min.

When comparing all the investigated models with the LR_EKL, one can note a clear difference in the performance. For the 1-hour predictions, the LR_EKL model is followed in performance by the NN, LR_NNAE (both with similar performance) and PM models, respectively. Although all the models presented a low FPR, the LR_EKL presented the highest TPR, which in this case translated into a higher TSS. Different factors may contribute to this.A low signal-to-noise ratio summed up with a limited number of geomagnetic events poses significant challenges for the development of prediction models. Improving the data pre-processing for the generation of the A_TID and increasing the database of events could be a good way to improve the preconditions for the development of the prediction models. It is important to highlight that all models presented higher performance than the baseline PM model for predictions beyond 60 min. The NN model (described in Sect. 2.5) presented better results than almost all the other models for predictions beyond 60 min. Its lower performance when compared to the LR_EKL, however, may be attributed to the higher level of underestimation in this model when higher levels of A_TID are observed, which leads to a decrease in TPR. One of the possible causes for the underestimation could be the scarce dataset showing a low number of storms with significantly high A_TID values, which could then affect the NN estimates. Also, due to the different characteristics and features of each geomagnetic storm event, the amount of geomagnetic storm events used to train the model may not be representative enough. The MMEs did not significantly improve the prediction results beyond 60 min. However, improvements were observed for prediction lead times shorter than 60 min, which suggests that the combination of the LR_EKL and the PM could be a useful solution for short-term predictions/nowcasting of the LSTID activity level.

4.1 Longitudinal variability

As presented in the previous sections, the analysis has so far been performed considering the GLSV station only. In order to investigate whether there is a difference in performance for different longitudinal sectors, we include the stations BOR1, HERS, and NVSK in the analysis, whose locations are presented in Figure 9.

Figure 9 Location of the GLSV, BOR1, HERS and NVSK GNSS stations. Geomagnetic latitudes of 30, 40, 50 and 70° N of are indicated by the dashed blue lines.

Given the better performance of the LR_EKL model presented for the GLSV station, we have selected this model to investigate how it performs for the other stations located at a similar geomagnetic latitude. The prediction performance is illustrated for the November 20, 2003 geomagnetic storm, and more than 50 storms are used to obtain the regression lines, which are shown in Figure 10a for this geomagnetic storm. The complete list of storms used in this analysis is presented in Appendix. In this case, we again use a leave-one-out approach to test the model. With this approach, the model is tested for predicting the A_TID index for more than 50 stations, and the results are summarized in the TSS plots for each station in Figure 10c. Similar to what was observed for GLSV, the A_TID predicted by the LR_EKL reflects the increase in the reference A_TID observed during the storm. The TSS considering all the storms for each station is quite similar to the results obtained for the GLSV station.

Figure 10 (a) Regression lines for predicting ATID index for the November 20, 2003 geomagnetic storm, (b) prediction of the ATID index for the November 20, 2003 geomagnetic storm and TSS (c) for the BOR1, HERS and NVSK stations obtained after evaluating the model on the storms described in the Appendix.

5 Conclusions

In this study, we investigated the performance of different methodologies in predicting the LSTIDs activity via the A_TID index, which is a useful index developed for statistical analysis of TIDs activity. The set of models investigated in this work includes: two models based on linear regression (LR_EKL and LR_NNAE), a neural network-based model ($\mathrm{NN}$), and the persistence model (PM), which is taken as the baseline model. The results show that all proposed models outperformed the baseline model PM for predictions beyond 1 h. For 1-hour predictions, all models presented a rather low FPR, which indicates that the difference in TSS for each model is associated with the TPR in this case. The investigation of the performance for other GNSS stations (BOR1, HERS and NVSK) has shown similar results to those obtained for the GLSV station, demonstrating the capabilities of the LR_EKL model to predict TID activity during geomagnetic storms for other longitudes at mid-latitudes in the European region. As already presented, the main goal of this study is to understand if and how much other state-of-the-art methodologies can improve the linear regression-based methodology for predicting the A_TID index at mid-latitude Europe during geomagnetic storms presented in Borries et al. (2023). The results have shown that more complex models do not necessarily lead to an increase in performance, as may be sometimes expected. For the GLSV station, the LR_EKL model outperformed the NN model in almost all lead times of prediction. Although the NN model provides the lowest FPR, the TSS obtained with this methodology could not surpass the TSS of the LR_EKL, except for nowcast activities.

It is important to note that although the results of the NN model for this particular study did not outperform the LR_EKL, it performed similarly to the LR_NNAE and better than the PM for predictions beyond one hour. This suggests that the NN model should not be completely discarded, but could also be improved in terms of hyperparameters optimization, number of neurons and layers, and selection of input parameters. In addition, the quality of the dataset in terms of coverage, noise, and cause-effect relationship has to be carefully taken into consideration when using such a methodology. The study also showed that multi-model ensembles may be a good alternative to improve the estimates of different models. Even the combination of simple models (LR_EKL and PM) led to an increase in performance for short-term predictions. In this specific case, this methodology provided a good improvement in performance for 30-minute prediction. The recent developments in machine learning algorithms and multi-model ensembles provide new possibilities for space weather prediction and modelling applications, and different investigations have shown their advantages and limitations. The already established simple methodologies, should not be disregarded and can, depending on the situation, still be considered powerful tools for approaching complex problems.

Acknowledgments

For generating the TID activity index, we used RINEX files from IGS for the dates before 2015 and processed line of sight TEC data from Madrigal database. For predicting the AE index, we used ACE solar wind data from ACE-MAG, ACE-SWICS and ACE-SWEPAM instruments www.srl.caltech.edu/ACE/ASC/level2/, and the auroral electrojet indices AE from https://wdc.kugi.kyoto-u.ac.jp/aeasy/index.html. The data is publicly available and we like to thank the data providers for sharing their data. The authors would like to thank the Brazilian agencies Federal District Research Support Foundation (FAPDF) and the National Council for Scientific and Technological Development (CNPq) that partially supported this work. The editor thanks three anonymous reviewers for their assistance in evaluating this paper.

References

Appendix

All Tables

All Figures

	Figure 1 (a) SYM-H index, (b) Bz component of the interplanetary magnetic field, (c) solar wind speed and (d) AE index for 19–21 November, 2003.
In the text

	Figure 2 Scatter plot of the maximum TID activity index of the GNSS ground station GLSV during each one of the investigated storm events versus the maximum Kan-Lee electric field EKL in the 18 h ahead of the maximum TID activity index (Borries et al., 2023).
In the text

	Figure 3 (a) EKL and (b) max2 h EKL for one day before, the day of the onset and one day after November 20, 2003 geomagnetic storm.
In the text

	Figure 4 Scatter plot of the maximum TID activity index of the GNSS ground station GLSV during each of the 56 storm events versus the maximum NNAE in the 18 h ahead of the maximum TID activity index (Borries et al., 2023).
In the text

	Figure 5 Schematic representation of the NN Model used on the ATID index prediction.
In the text

Figure 6 Left panels: One-hour prediction of the TID activity index during the geomagnetic storm event registered on the 20th of November 2003 using: (a) LREKL, (c) LRAE, (e) PM and (g) NN model. Right panels: The same predictions presented in the left panel, but evaluated as a classification problem, showing True Positives (TP), True Negatives (TN), False Positives (FP) and False Negatives (FN) for a threshold of 0.5 TECU obtained using: (b) LREKL, (d) LRAE, (f) PM and (h) NN model.

In the text

	Figure 7 (a) TSS for the different models under investigation, i.e. regression (obtained from EKL and AE), NN and the PM and (b) the TSS obtained from the combination of some of such models using a weighted average scheme.
In the text

	Figure 8 Statistics obtained from the different models for 60 min prediction: (a) True Positive Rate, (b) False Positive Rate and (c) True Skill Statistic.
In the text

	Figure 9 Location of the GLSV, BOR1, HERS and NVSK GNSS stations. Geomagnetic latitudes of 30, 40, 50 and 70° N of are indicated by the dashed blue lines.
In the text

	Figure 10 (a) Regression lines for predicting ATID index for the November 20, 2003 geomagnetic storm, (b) prediction of the ATID index for the November 20, 2003 geomagnetic storm and TSS (c) for the BOR1, HERS and NVSK stations obtained after evaluating the model on the storms described in the Appendix.
In the text