# Automatic Detection of Large-Scale Traveling Ionospheric Disturbances Using GNSS Data and Image Processing Techniques

Marco Guerra, Claudio Cesaroni, Nicola Fiorentino, Federico Tosone, Vincenzo Ventriglia, Emanuele Pica, Elvira Astafyeva, Boris Maletckii, Luca Spogli

► To cite this version:

Marco Guerra, Claudio Cesaroni, Nicola Fiorentino, Federico Tosone, Vincenzo Ventriglia, et al. Automatic Detection of Large-Scale Traveling Ionospheric Disturbances Using GNSS Data and Image Processing Techniques. Space Weather: The International Journal of Research and Applications, 2025, 23 (9), pp.e2025SW004423. 10.1029/2025SW004423. hal-05382202

HAL Id: hal-05382202

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

Submitted on 26 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Space Weather®

RESEARCH ARTICLE

10.1029/2025SW004423

**Key Points:**

* An automatic method is developed to detect large-scale traveling ionospheric disturbances (LSTIDs) from GNSS dTEC keograms using image processing and RANdom SAmple Consensus fitting

* A 4-year LSTID database over Europe is created, enabling detailed statistical analysis and comparison with past studies

* Clear relationships are found between LSTID properties and geomagnetic indices, confirming known generation mechanisms

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. Guerra,
marco.guerra@ingv.it

**Citation:**

Guerra, M., Cesaroni, C., Fiorentino, N., Tosone, F., Ventriglia, V., Pica, E., et al. (2025). Automatic detection of large-scale traveling ionospheric disturbances using GNSS data and image processing techniques. *Space Weather*, 23, e2025SW004423. https://doi.org/10.1029/2025SW004423

Received 7 MAR 2025
Accepted 10 SEP 2025

**Author Contributions:**

**Conceptualization:** Marco Guerra,
Claudio Cesaroni, Vincenzo Ventriglia,
Luca Spogli

**Data curation:** Marco Guerra,
Boris Maletckii

**Formal analysis:** Marco Guerra,
Emanuele Pica

**Funding acquisition:** Luca Spogli

**Investigation:** Marco Guerra,
Claudio Cesaroni

**Methodology:** Marco Guerra,
Claudio Cesaroni, Nicola Fiorentino,
Federico Tosone, Vincenzo Ventriglia,
Luca Spogli

© 2025. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

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

# Automatic Detection of Large-Scale Traveling Ionospheric Disturbances Using GNSS Data and Image Processing Techniques

Marco Guerra¹,² ID, Claudio Cesaroni¹ ID, Nicola Fiorentino³ ID, Federico Tosone³ ID, Vincenzo Ventriglia¹ ID, Emanuele Pica¹ ID, Elvira Astafyeva⁴ ID, Boris Maletckii⁴,⁵ ID, and Luca Spogli¹ ID

¹Istituto Nazionale di Geofisica e Vulcanologia, Rome, Italy, ²Sapienza University of Rome, Rome, Italy, ³Spindox, Rome, Italy, ⁴Institut de Physique du Globe de Paris, Paris, France, ⁵Earthquake Research Institute, The University of Tokyo, Tokyo, Japan

**Abstract** Large-scale traveling ionospheric disturbances (LSTIDs) are wave-like structures propagating through the ionosphere from auroral to lower latitudes, affecting radar-based technologies and satellite navigation. Due to their disruptive effects, timely detection is crucial. This work presents an autonomous detection technique based on travel-time diagrams of detrended Total Electron Content (TEC). Our method employs a network of GNSS receivers across Europe, selected using a K-Means clustering algorithm to ensure comprehensive spatial coverage. Applying image processing techniques, we identify propagating features of positive and negative amplitude in detrended TEC keograms. These features are fitted using the RANdom SAmple Consensus method to estimate propagation speed along the North-South direction. Following feature identification, we refine our results by merging positive and negative perturbations and filtering out false positives based on physical constraints. This methodology allowed us to compile a catalog of LSTID with periods ranging from 30 to 90 min, spanning 2021–2024. To validate our approach, we conducted case studies and statistical analyses, confirming accuracy and reliability. Additionally, testing in North America and Japan demonstrated the algorithm broader applicability. Our findings, obtained with user-defined parameters tailored to examine LSTIDs with periodicities in the 30–90 min range, show an average speed of 629 ± 215 m/s and period of 62 ± 14 min. A direct relation with auroral electrojet indices confirms that stronger auroral currents produce faster, higher amplitude and shorter period LSTIDs. Furthermore, our analysis reveals a link between LSTID occurrences and geomagnetic activity, reinforcing previous findings.

**Plain Language Summary** Large Scale Traveling Ionospheric Disturbances (LSTIDs) are wave-like disturbances in the Earth's upper atmosphere that can disrupt GPS and communication systems. They are mainly caused by solar storms and geomagnetic activity, making their detection crucial for space weather monitoring. This study presents an automated technique to detect LSTIDs using data from GNSS receivers. By analyzing variations in Total Electron Content, we identify LSTIDs through image processing and machine learning techniques. Our method was applied to 2021–2024 data, showing that using user-defined parameters tailored to reproduce previous findings our LSTIDs results in Europe show an average speed of 629 m/s and a period of around one hour. We highlight a strong direct relation between LSTID activity and geomagnetic storms, confirming that intense auroral currents generate faster and more energetic disturbances. Testing on one case event in the United States and Japan demonstrated that the method is reliable across different regions. This fully automated approach improves ionospheric monitoring and could be integrated into real-time space weather forecasting, helping mitigate the effects of ionospheric disturbances on satellite and radar-based technologies.

## 1. Introduction

The Earth's ionosphere is a highly dynamic region, constantly shaped by solar activity, geomagnetic forces, and atmospheric interactions. While long-term variations in solar flux and thermospheric composition drive gradual changes (Forbes et al., 2000; Rishbeth & Mendillo, 2001; Spogli et al., 2024), short-term disturbances frequently occur due to various geophysical phenomena (see e.g., Alfonsi et al., 2013; Astafyeva, 2019; Mendillo, 2006; Spogli et al., 2019; Takahashi et al., 2020). Among the ionospheric disturbances, Traveling Ionospheric Disturbances (TIDs) stand out as a key manifestation of ionospheric variability. As of now, TIDs have been studied for many decades through many different sensing technologies (Afraimovich et al., 2000; Chan & Villard, 1962;
**Project administration:** Luca Spogli

**Resources:** Marco Guerra,

Claudio Cesaroni

**Software:** Marco Guerra,

Nicola Fiorentino, Federico Tosone,

Emanuele Pica

**Supervision:** Claudio Cesaroni,

Elvira Astafyeva, Luca Spogli

**Validation:** Marco Guerra

**Visualization:** Marco Guerra

**Writing – original draft:** Marco Guerra,

Claudio Cesaroni, Nicola Fiorentino,

Federico Tosone, Vincenzo Ventriglia,

Emanuele Pica, Elvira Astafyeva,

Boris Maletckii, Luca Spogli

**Writing – review & editing:**

Marco Guerra, Claudio Cesaroni,

Nicola Fiorentino, Federico Tosone,

Vincenzo Ventriglia, Emanuele Pica,

Elvira Astafyeva, Boris Maletckii,

Luca Spogli

Davis, 1971; Francis, 1974; Georges, 1968; Habarulema et al., 2013; Haralambous et al., 2023; Hines, 1960; Hocke & Schlegel, 1996; Hunsucker, 1982; Munro, 1948; Shiokawa et al., 2007; Verhulst et al., 2022). TIDs are normally caused by acoustic-gravity waves propagating in the thermosphere that displace the ionospheric plasma through particle collisions. Three families of TIDs have been identified: large-scale, medium-scale (MSTIDs), and small-scale TIDs (SSTIDs), with their parameters (period, wavelength, and speed) being the discerning factor. While SSTIDs are the least studied and only recently started to catch the attention of the scientific community (see, e.g., Boyde et al., 2022), much work has been done to investigate the generation and propagation of MSTIDs and LSTIDs. LSTIDs are usually defined by periods ranging from 45 to 180 min, speeds of 400 to 1,000 m/s, and corresponding wavelengths that can be larger than 1,000 km (Borries et al., 2009; Ding et al., 2008; Habarulema et al., 2013; Tsugawa et al., 2004). MSTIDs are the main phenomenon during quiet periods. They are predominantly linked to lower atmospheric sources, natural hazards (e.g., Astafyeva, 2019), variations in the ionizing radiation from the Sun and electrodynamical processes, such as Perkins instability associated with E- and F-region coupling processes (Makela & Otsuka, 2012). In contrast, LSTIDs are mainly linked to geomagnetic storms and substorms. Originating primarily within the auroral region during increased geomagnetic activity, LSTIDs tend to propagate equatorward along geomagnetic field lines (Borries et al., 2009; Cesaroni et al., 2017; Frissell et al., 2022; Hunsucker, 1982; Jing & Hunsucker, 1993; Zakharenkova et al., 2016). Although the full chain of physical mechanisms behind LSTID generation remains complex and not fully understood, the primary cause is likely ascribed to AGWs produced by thermospheric heating associated with Joule heating, particle precipitation and Lorentz forces (Banks, 1977). Oyama and Watkins (2012) demonstrated that these forces impact wind velocity and induce temperature fluctuations with different magnitudes. Lorentz forces are responsible for increased wind velocities, while Joule heating causes a rise in thermospheric temperature. As the excitation phenomena are linked to the auroral current system, auroral electrojet indices (AE) should proxy the LSTID generation potential. Borries et al. (2009) studied around 40 geomagnetic storms and showed a correlation between AE and LSTID amplitude. Some examples of poleward and quiet LSTIDs were recently identified (Habarulema et al., 2022). While the first case can also be connected to the typical generation process, as storm LSTIDs can propagate down/up to the opposite hemisphere, the quiet time LSTIDs might be related to tropospheric phenomena such as deep convective plumes (Vadas & Crowley, 2010; Vadas & Liu, 2009).

In the early years, most studies of TIDs relied on ionosondes and Doppler sounders. However, with the advent of GNSS technology, Total Electron Content (TEC) measurements have become a widely used tool for studying TIDs. The extensive global network of GNSS receivers enables precise reconstruction of the horizontal propagation of LSTIDs. Despite its advantages, GNSS TEC measurements have a key limitation compared to HF-based techniques: being an integrated quantity, it cannot resolve the vertical structure of traveling waves. TEC measurements are typically detrended to emphasize wave-like features and remove diurnal trends. The resulting detrended TEC (dTEC) data can then be used to produce dTEC maps and keograms. These keograms, which visualize disturbances over space and time, are especially useful for studying LSTIDs. Since LSTIDs propagate preferably along magnetic field lines, approximated as parallel to the North-South direction at mid-latitudes, latitude is commonly used as the spatial variable in these analyses. As the magnetic field declination is around 4°–6° for the European area, the meridional velocity estimation error is negligible when mapping the magnetic N-S movement to the geographic N-S direction (Borries et al., 2009).

Several studies have addressed the automatic detection of LSTIDs using ionosonde data. For example, the TechTIDE project developed various methodologies for the real-time detection, characterization, and tracking of TIDs based on ionosonde and other data sources (Belehaki et al., 2020). One of these methods specifically leverages ionosonde measurements to identify and track LSTIDs in real-time, as Altadill et al. (2020) demonstrated. More recently, Segarra et al. (2025) presented a climatological analysis of LSTID activity over Europe from 2014 to 2023 using the HF-INT method, which enables near-real-time detection from a European ionosonde network. Their Catalog of LSTID events includes onset times, durations, dominant periods, and propagation vectors, highlighting the dependence of LSTID activity on local time, season, and geomagnetic conditions. These studies confirm the viability of ionosonde-based automatic detection systems and their relevance in operational space weather services. Other independent TEC-based methods not originally designed for LSTID detection have also shown the capability to capture such disturbances under specific conditions. For example, a recent study successfully detected TIDs associated with the Lamb wave generated by the 2022 Tonga eruption in the far field, which, despite not being geomagnetically induced, can be classified as large-scale TIDs due to their propagation characteristics (Maletckii & Astafyeva, 2022). Furthermore, another approach based on high dTEC/dt values was
<table><caption>**Table 1**
Summary of the Characteristics and Results of Statistical Analysis Performed in Previous Studies</caption><thead><tr><th></th><th>Geographical area</th><th>Period considered</th><th>LSTID speed [m/s]</th><th>LSTID period [min]</th></tr></thead><tbody><tr><td>Tsugawa et al. (2004)</td><td>Japan</td><td>April 1999–December 2002</td><td>475 ± 171</td><td>80 ± 29</td></tr><tr><td>Ding et al. (2008)</td><td>North America</td><td>Major magnetic storms from 2003 to 2005</td><td>300</td><td>108</td></tr><tr><td>Borries et al. (2009)</td><td>Europe</td><td>Major magnetic storms from 2001 to 2007</td><td>684 ± 195</td><td>59 ± 11</td></tr></tbody></table>

able to detect storm-related LSTID activity, as discussed in (Maletckii & Astafyeva, 2024). These results demonstrate that flexible detection systems, even when not explicitly tailored for LSTIDs, can still identify them effectively under favorable conditions.

This work introduces an autonomous detection algorithm designed to identify LSTID signatures in TEC data derived from GNSS observations. The method operates in both geomagnetically quiet and disturbed conditions, and it is optimized to extract key LSTID parameters such as the horizontal apparent speed (along the N-S direction), period, and wavelength from dTEC keograms. To assess its performance, the algorithm was applied to the period from January 2021 to December 2024 over the European area. This enabled a comprehensive statistical analysis of LSTID occurrence, allowing comparison with known climatological patterns reported in the literature.

In particular, Tsugawa et al. (2004) visually inspected dTEC maps over Japan to identify LSTID signatures from April 1999 to December 2002. Ding et al. (2008) performed a similar study by visually inspecting dTEC maps of North America for major storms from 2003 to 2005. Lastly, Borries et al. (2009) performed an automated analysis of LSTIDs derived from dTEC maps and keograms of the European area for 37 geomagnetic storms. The results of the aforementioned studies are summarized in Table 1.

Previous studies that leveraged TEC either required visual data inspection or focused only on major geomagnetic storms. In contrast, our technique no longer requires visual inspection and was applied to identify LSTIDs under quiet and disturbed conditions. The automatization of our technique raises some difficulties, such as the robustness against false positives; on the other hand, it does not include the visual inspection bias and also allows identifying LSTIDs that are not directly related to geomagnetic storms, as we applied it to all the days in the considered period.

The paper is structured as follows: the Data and Methods section describes the steps of the automatic detection algorithm and the creation of detrended TEC maps. This section also details the image processing techniques used for identifying regions of interest (ROIs) and extracting LSTID parameters such as speed and period. The Results and Discussion section presents the outcomes of applying the detection algorithm to data from 2021 to 2024. It also includes two case studies (one during a major geomagnetic storm and another during a moderately disturbed period) to illustrate the algorithm's performance under different conditions. For the major geomagnetic storm, we also analyzed the performance of the detection algorithm in other geographical regions, specifically Japan and North America. Then, we present a statistical analysis of LSTID occurrences, propagation speeds and periods, highlighting their dependence on geomagnetic indices like Hp30, SymH, and the IMAGE auroral electrojet index (IE). The strengths and weaknesses of the proposed methodology are also discussed. The Conclusions and Future Work section summarizes the key findings, emphasizing the effectiveness of the automated detection technique and its agreement with previous studies. It discusses limitations, such as difficulties in detecting overlapping disturbances or events spanning two consecutive keograms. We also suggest potential improvements, including integrating machine learning for forecasting LSTIDs and developing real-time detection capabilities.

## 2. Data and Methods

Figure 1 outlines the automated LSTID detection process. It begins with Data Gathering and Station Selection, where dual-frequency GNSS data are retrieved from an optimized subselection of stations belonging to the EUREF network. Next, in Detrended TEC and keogram Images, phase measurements are used to derive slant TEC, which is then verticalized and detrended with a Savitzky-Golay filter. These detrended TEC (dTEC) values are organized into latitude-time keograms for further analysis. The ROI Identification Process identifies potential LSTID features using image processing. Binary masking and segmentation techniques, including Bayesian Gaussian Mixture Models (BGMM), isolate coherent positive and negative deviations. The ROI Feature
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Flowchart of the large-scale traveling ionospheric disturbances detection method.</figcaption></figure>

Extraction step estimates propagation speed and other parameters using the RANSAC algorithm to fit lines to the ROIs, providing metrics like start time, latitude range, and velocity. Following this, the algorithm joins opposite ROIs to retrieve LSTID occurrences, filtering out false detections based on criteria like speed, period, and latitudinal overlap.

## 2.1. Data Gathering and Station Selection

Our detection technique leverages dTEC measurements to identify LSTID features in latitude-time images. Image retrieval starts from uncalibrated verticalized TEC derived from GNSS phase measurements. Thus, the first necessary step is to retrieve dual-frequency phase measurements from the area of interest. We exploit multi-constellation RINEX files in 30 s format gathered by the EUREF Permanent GNSS Network receivers (Bruyninx et al., 2019). As keograms are made by selecting observations pertaining to a certain longitudinal range, we collect data from receivers located between 10° and 25°E and in the latitudinal range 30° to 60°N.

To reduce computational time and manage data size when deriving keograms, we downsized the original EUREF GNSS network by selecting a subset of stations. Instead of randomly choosing stations, we employed a clustering algorithm to ensure adequate spatial coverage using the K-Means algorithm (Lloyd, 1982; MacQueen, 1967). K-Means is an unsupervised learning algorithm that groups data into a predefined number of clusters, 100 in our case. Using 100 stations, we obtained a good compromise between the Ionospheric
<figure><img src="image_4.png" /><figcaption>**Figure 2.** EUREF GNSS permanent network (black dots) and our selection of GNSS stations performed by applying K-Means clustering (red dots). The blue shaded area represents the geographical region considered for generating the keograms.</figcaption></figure>

Pierce Point (IPP) coverage of the area of interest and the density of IPP per pixel of the keogram image. More details are provided in Text S5 in Supporting Information S1, where we show IPP density keograms for different numbers of stations. The station selection process begins by randomly selecting 100 stations as initial cluster centroids. Each GNSS station from the full network is assigned to the nearest centroid, forming clusters. The centroids are iteratively recalculated as the mean position of the stations in their respective clusters. This process continues until the assignments stabilize, meaning no stations change clusters during an iteration. The primary goal of this clustering was to eliminate redundant data from densely covered regions while preserving the spatial distribution necessary to capture ionospheric variations. To mitigate the impact of randomly initialized centroids, we repeated the K-Means algorithm 5,000 times with different starting points. The solution minimizing the average distance between stations and their cluster centroids was selected among these iterations. This ensured a robust and well-optimized distribution of the selected stations. Additionally, each year, the station selection was repeated to account for temporal changes in the network configuration. For example, the selected stations for 2024 are shown in Figure 2. Our keogram generation process operates by analyzing data one day at a time. Thus, to remove the impacts of measurements gathered around 00:00 and 24:00 UT of the given day, we include the last/first 2 hr of the previous/following day in our process.

To analyze the feasibility of applying the detection algorithm outside Europe, we also applied it to two different geographical regions for the major geomagnetic storm of 11–12 August 2024. The data used belong to GEONET (Sagiya, 2004) for Japan and CORS (National Geodetic Survey, 2025) for North America. As GEONET and CORS comprise thousands of GNSS stations, we applied the same station selection process that was applied to EUREF. Figure 3 shows the spatial distribution of the 100 stations selected for the two regions.

## 2.2. Detrended TEC and Keogram Images

Once the RINEX files for the optimized network have been collected, we calculate the detrended verticalized uncalibrated TEC. To begin, dual-frequency phase measurements belonging to the four GNSS constellations (GPS, GLONASS, Galileo and BeiDou) are combined to obtain the Geometry-Free Linear Combination (GFLC) of phase measurements *s*TECph according to Equation 1:
<figure><img src="image_5.png" /><figcaption>Figure 3. GNSS stations considered for the performance analysis of our algorithm for two different geographical regions. Red dots show the 100 stations selected from the CORS network, while blue dots are the GEONET ones.</figcaption></figure>

$$sTEC_{ph} = \frac{1}{A} \frac{f_1^2 f_2^2}{f_1^2 - f_2^2} (L_1 \lambda_1 - L_2 \lambda_2) \quad (1)$$

$A = 40.308 \text{ m}^3\text{s}^{-2}$, $L_1$ and $L_2$ are the phase measurements, $\lambda_1$ and $\lambda_2$ are the wavelengths of the two considered frequency bands $f_1$ and $f_2$ (e.g., E1 and E5a from Galileo). A thin-shell ionosphere assumption is made to derive the coordinates of the IPPs related to the GFLC measurements (Mannucci et al., 1998). The IPPs correspond to the intersection between the receiver-satellite line of sight and the ionospheric 2D layer. As the plasma density distribution is mainly located in a thin layer, it is fair to assume the whole ionosphere is a 2D layer in latitude and longitude. We set the ionospheric shell height to 300 km and derived the IPP location accordingly. Then, we verticalize the GFLC of each continuous observational arc after setting the initial point of observation to 0 TECu. According to Guerra et al. (2024), this should allow the best compromise between amplitude error after the observable is detrended and computational speed. The mapping function represents the relationship between the slant and verticalized observables:

$$vTEC_{ph} = sTEC_{ph} \cos \left( \arcsin \left( \frac{R_e \cos \theta}{R_e + H_{IPP}} \right) \right) \quad (2)$$

where $R_e$ is the Earth's radius, $\theta$ is the elevation angle, and $H_{IPP}$ is the ionospheric shell height. To reduce the impact of cycle slips related to multipath, only observations with an elevation higher than 20° are considered at this stage. Once the TEC observable is obtained, a detrending technique is applied to highlight periodic features ascribable to LSTID. Again, following the findings of Guerra et al. (2024), we implemented a Savitzky-Golay filter (Schafer, 2011), as it was found to ensure the best performance with respect to detrending errors and computational weight. Specifically, we employ a second-order filter with a time window of 95 min. The rationale behind choosing the 95-min window is as follows: Guerra et al. (2024) demonstrated that optimal performance was achieved with a 120-min window when identifying a modeled LSTID with a period of 75 min. According to Borries et al. (2009), we expect European LSTIDs to have, on average, a period of approximately 60 min. Therefore, by maintaining the same window-to-period ratio, we deduce that the optimal window for our case would be 95 min. This choice of window size imposes a limitation on the LSTIDs that our technique can identify: a larger time window would enable the identification of LSTIDs with longer periods. However, since we expect most LSTIDs to have periods shorter than approximately 90 min, we accept this trade-off. Moreover, our
<figure><img src="image_6.png" /><figcaption>**Figure 4.** An example of the result of the Region of Interest identification step for the 15th of October 2024. Panel (a) shows the keogram image colored according to dTEC amplitude, with red/blue area representing positive/negative deviations from the background trend. Panel (a) also shows the output of the RANdom SAmple Consensus fitting algorithm in gray, a straight line that models the North-South propagation of Total Electron Content disturbances. Panel (b) shows the different ROIs with their pixels in different colors.</figcaption></figure>

technique could be applied with different parameters (such as adjusting the detrending window) to allow targeted studies of longer-period LSTIDs. In addition to Savitzky-Golay, we apply a moving average filter on a time window of 15 min to reduce the high-frequency noise. As detrending techniques are known to suffer from border effects, we now apply a further elevation cutoff of 30°.

Starting from the dTEC measurements for a given day, it is possible to create the keogram by using the data pertaining to the latitude and longitude range of 37.5° to 57.5°N and 10° to 25°E. After this, the median of all measurements falling in each 0.05° bin of latitude is calculated and assigned to the relative latitude/time pixel. The data from the last/first 2 hours of the previous/following day are removed, thus obtaining a keogram free of boundary effects. As an example, Figure 4a shows the keogram image obtained for the 15th of October 2024. The final step before the ROI identification process (Section 2.3) is an interpolation of dTEC on a grid equal to the keogram image size. By doing this, the dTEC values for pixels not covered by observations are filled with linear interpolation and we obtain a processable 2D image. Once the database of keogram images is generated for January 2021 to December 2024, it is possible to proceed with the ROI identification process.

## 2.3. Region of Interest (ROI) Identification Process

Being wave-like structures, LSTIDs can be identified as oblique bands of aligned positive or negative dTEC. Therefore, to identify an LSTID in a keogram, the first step would be identifying connected regions characterized
by alternating negative/positive amplitude. The regions that are potential LSTID features are hereafter referred to as ROIs. Before applying the ROI identification algorithm, some preprocessing is performed.

First, we perform a local average downsampling of size 4 in time. As the original keogram image is 401 × 2880 pixels, the downsampled one will be 401 × 720. This preprocessing has two major advantages: the first is that ROI identification is faster due to the reduced image size, and the second is that the averaging step reduces image noise, which is detrimental for the instance segmentation algorithm that will be later described and applied. A step-by-step process is used to identify the ROIs in the preprocessed keogram images. First, all the values with an absolute amplitude higher than 0.015 TECu are used to create one negative and one positive binary mask. The two masks are defined as follows:

$$B_{+}(x,y) = \begin{cases} 1, & M(x,y) > T \\ 0, & M(x,y) \le T \end{cases} \qquad (3)$$
$$B_{-}(x,y) = \begin{cases} 1, & M(x,y) < T \\ 0, & M(x,y) \ge T \end{cases}$$

Where $B_+$ and $B_-$ are the two binary masks, $M$ is the preprocessed matrix, and $T$ is the amplitude threshold set to 0.015 TECu. ROIs are clusters of spatially connected pixels within the binary masks that satisfy certain conditions. Two pixels are considered connected if they satisfy a 2-connected adjacency criterion, meaning they can be linked through a path consisting of one or two orthogonal (horizontal or vertical) steps. The connected regions that have a latitudinal range smaller than 10° and a solidity lower than 0.35 are discarded. Here, solidity refers to the ratio between the pixels belonging to the connected region and the total number of pixels inside the convex hull. The connected regions that are characterized by a minor axis of the fitted ellipse longer than 25 pixels (50 min) are checked for the size of their largest concavity. To fit the ellipse, we build a covariance matrix from the pixel coordinates. This matrix captures how the pixel positions vary along the horizontal, vertical, and diagonal directions. The ellipse parameters are then derived by performing an eigenvalue decomposition of this covariance matrix. The eigenvectors of the matrix indicate the orientation of the ellipse, while the eigenvalues correspond to the squared lengths of the major and minor axes. The axis with the larger eigenvalue represents the direction in which the region is most elongated (the major axis), and the smaller eigenvalue corresponds to the minor axis. To evaluate the size of the largest concavity, we compare the convex hull with the alphashape derived concave boundary (Edelsbrunner et al., 1983) and retrieve the region concavities. If the largest concave area is bigger than 700 pixels, the connected region undergoes a second step of instance segmentation using the BGMM clustering algorithm applied to the normalized coordinates of points. Other clustering algorithms have been evaluated, with BGMM showing the best performance. However, the comparison between different techniques is not included in the text. The BGMM (Roberts et al., 1998) is a probabilistic model that assumes data points come from a mixture of K Gaussian distributions. The parameters of the Gaussian distributions (such as means and covariances) are treated as random variables in BGMM, and Bayesian inference methods are used to estimate their values. The number of clusters $K$ is gradually increased, and the algorithm stops when the performance improvement from increasing $K$ by one is lower than in the previous step. This unsupervised approach ensures that the segmentation is not overly fragmented and helps avoid creating too many clusters when the additional subdivision does not improve the model significantly. The irregular region is then divided into $K$ subregions by disconnecting the pixels belonging to the intersection of each couple of the $K$-connected subregions. The new $K$ subregions are checked again based on the solidity and latitudinal elongation criteria. This secondary segmentation step tackles the scenario typically related to two consecutive LSTID bands of the same magnitude connected by noise at higher latitudes that appear as a single region. More details are provided in Text S4 in Supporting Information S1. All the connected regions obtained at this stage are now filtered based on the inclination of the major axis of the fitted ellipse. The connected regions are discarded if the inclination is smaller than 75° (corresponding to a meridional speed of ~150 m/s). All the remaining connected regions are now considered ROIs. The output of the ROI identification process is shown in Figure 4b. Figure 4 shows how the technique can detect most of the band-like structures for the 15th of October 2024. Panel (b) of Figure 4 shows a large number of ROIs in different colors; nevertheless, some are not ascribable to LSTID occurrences. For instance, feature 5 is likely related to the
passage of the solar terminator, while features 6, 7, 8 and 9 are likely connected to daytime MSTIDs. These are discarded in the following steps of the detection algorithm, where ROI pairing is performed.

## 2.4. ROI Feature Extraction

Once the ROIs of a given image are identified, it is possible to derive some of their parameters. Since a keogram is a space-time representation of dTEC, the slope of a given ROI corresponds to the meridional component of the horizontal velocity. To derive the slope of ROIs, a straight line is fitted to the ROI pixels with the RANSAC algorithm (Fischler & Bolles, 1981), which excludes outliers in the fitting process to ensure a higher level of robustness with respect to other fitting techniques.

RANSAC is an iterative algorithm employed to derive the parameters of a given model. In our case, the model considered is a weighted linear regression fit (Draper & Smith, 2014). The weights employed in the fitting process are set according to the pixel dTEC intensity. First, a random subselection of the points S is chosen and the model parameters $\theta_s$ are estimated:

$$ \theta_s = \operatorname{argmin}_{\theta} \sum_{x \in S} \operatorname{loss}(x, \theta) \quad (4) $$

Then, the consensus is evaluated by calculating the number of points $x_i$ of the full set X characterized by an error lower than $\epsilon$:

$$ C_s = \{x_i \in X : \operatorname{loss}(x_i, \theta_s) < \epsilon\} \quad (5) $$

After the consensus is evaluated for a given number of subselections, the set of $\theta_s$ that maximizes consensus is chosen. Thanks to RANSAC, it is possible to derive the ROI start time, end time, latitude range, slope/velocity and the fitting quality metrics (mean absolute error, root mean squared error, and R2 score). As shown in panel (a) of Figure 4, RANSAC can correctly fit the middle line of all the band-like features related to LSTIDs.

## 2.5. From ROIs to LSTIDs

In our assumptions, LSTIDs are wave-like structures propagating at a velocity above 300 m/s with periods longer than 30 min. Thus, our ROIs cannot be considered LSTIDs, as a single ROI represents only positive or negative variations. Therefore, the ROIs obtained from daily keograms are merged into one single database and undergo further filtering processes to ensure that they are related to LSTID occurrences and not to noise or other ionospheric processes that are not of interest. A given ROI is considered an LSTID ROI if an ROI of the opposite magnitude exists and meets the following requirements:

* The difference in time propagated at 47.5°N is smaller than 45 min

* The speed is in the range of 300–1,500 m/s

* The slope angle difference is smaller than 2°

* The latitudinal overlap is at least 7.5°

All the ROIs that meet the requirements are considered LSTID signatures. The high threshold for the velocity range (1,500 m/s) was chosen due to the high sensitivity of the slope-based velocity retrieval, as working with angles close to 90°, the tangent grows dramatically for small angle differences. Moreover, we acknowledge the fact that we discard all potential LSTIDs defined by a period longer than 90 min. However, as specified in the previous sections, the technique parameters are not considered to be fixed. In contrast, they allow users to tailor the methodology to different periodicities or velocities. We show how different parameters impact the methodology results on our 4-year database in Supporting Information S1. After this final filtering process, an LSTID catalog is obtained, enabling statistical characterization analyses.

## 2.6. From LSTID ROI Identification to Statistical Analysis

LSTID ROI occurrences are investigated in relation to the hour of the day and the month to highlight potential influences of the background ionization and relative solar wind/geomagnetic field geometry. Moreover, the relationship between geomagnetic indices and the LSTID detection occurrences and rates is investigated.
The magnetic indices considered are Hp30, SymH and IMAGE auroral electrojet indices (IE). The geomagnetic Hp30 index, described by Yamazaki et al. (2022), is similar to the widely adopted Kp index, but offers higher time resolutions, specifically 30 min. Moreover, unlike the Kp index, which is capped at a maximum value of 9, the Hp30 index is open-ended, allowing it to better represent the intensity of severe geomagnetic storms in greater detail than the three-hourly Kp index. The two indexes, Hp30 and Kp, are strongly correlated. For example, for the 11–12 August 2024 storm, the maximum Kp reached 8 between 12 and 15 UT, while Hp30 had an average value for the same period of 7.88.

The Symmetric H-component (SYMH) index is a key indicator of ring current disturbances during geomagnetic storms, utilizing data from a longitudinal array of low-latitude ground-based magnetometers (Iyemori, 1990). SYMH is similar to the Disturbance Storm Time (Dst) index as it characterizes the ring current disturbances; however, it offers a refined 1-min temporal resolution, making it particularly valuable for investigating short-term fluctuations in geomagnetic disturbances.

Lastly, IE is an alternative to the classic Kyoto AE index (Tanskanen, 2009) provided by the Finnish Meteorological Institute (FMI). The main difference is that the IE encompasses a larger network of magnetometers distributed across the northern European area, thus allowing a more significant representation of the area where European LSTIDs are normally generated. In addition to the statistical study and geomagnetic activity influences, thanks to the RANSAC fitting procedure, it is also possible to estimate the velocity and period distribution along with their dependence on geomagnetic conditions. Before moving to the statistical results, the definition of the LSTID ROI speed and period must be clarified. Any given LSTID ROI has a related period and speed. While the speed is derived by converting the slope to a meridional velocity, the LSTID period is twice the distance between the given ROI and the adjacent and opposite amplitude ROIs. For instance, ROI 10 and 11 of Figure 4 will have different periods since ROI 10 will be characterized by a period twice the time distance between ROI 10 and 11. In contrast, the ROI 11 period will be twice the average time distance between ROI 11 and ROIs 10 and 12.

## 3. Results and Discussion

This section will present two case events, one severe storm and one moderately disturbed period, that will help show the reliability of our technique and its ability to recognize LSTID of different characteristics and in different geographical areas. Further, to get a broader perspective, the results over the whole period are analyzed by looking at the detection occurrences and probability with respect to month, hour of the day and geomagnetic conditions. Moreover, the statistical behavior of LSTID speed, period and amplitude is presented. Lastly, we present a statistical analysis to evaluate the relationship between the geomagnetic indices and the LSTID detection rate.

### 3.1. Examples of LSTID ROIs Identification

The first example of LSTID ROI identification in a keogram image corresponds to the extreme geomagnetic storm that occurred on 11–12 August 2024, as shown in Figure 5. During this event, the SymH index reached a minimum of approximately –220 nT, IE grew to 1,940 nT, and the Hp30 index peaked at 9. Figure 5 highlights that the generation of numerous LSTIDs characterized this storm.

The first notable LSTID activity is a low-amplitude wave train propagating at 500–600 m/s, observed around 21:00 on the 11th of August. This train consists of two negative and one positive deviation from the background. Following this, after 22:00, many continuous fluctuations of varying amplitudes and speeds propagate consistently until 03:00 on the 12th of August, and the technique correctly detects these. However, at approximately 02:00, our method misidentifies a positive ROI due to high noise levels between latitudes 45° and 50°. This results in the ROI being incorrectly split into two parts: the lower segment is too short in latitude, and the upper segment fails to overlap with adjacent positive ROIs. Consequently, both parts are discarded. A similar issue arises with the broad negative ROI around 04:00, which is discarded due to a failure in the BGMM step of the segmentation process. Here, the model fails to separate two consecutive overlapping negative disturbances. Around 05:00, an ambiguous detection occurs with a negative dTEC feature propagating at a speed of 568 m/s. Additionally, before the correctly detected LSTID ROIs between 12:00 and 17:00 on the 12th of August, several northward-propagating features were identified between 11:00 and 12:00. These are discarded due to their negative speed values. Although these features are likely LSTIDs crossing the equator, they were excluded based on our velocity criteria. Adjusting this criterion would allow the detection of poleward-propagating LSTIDs, which could be useful for studying such phenomena. Later, around 19:00, two consecutive features (a negative deviation
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The latitude-time evolution of dTEC disturbances for the 11–12 August 2024 storm, with lines representing the RANdom SAmple Consensus fitted line. White dashed lines indicate detected regions of interest (ROIs) that are not considered Large-scale traveling ionospheric disturbance (LSTID) ROIs, while white solid lines are ROIs that satisfy the LSTID conditions and are plotted along with their corresponding speed in m/s shown in the white boxes.</figcaption></figure>

followed by a positive one) appear to form an LSTID but are not fully detected. This is because only the negative feature is identified as a connected region, likely due to lower noise levels (as indicated by the white dashed line around 19:00), especially at latitudes above 45°. The identified negative feature is discarded because no corresponding positive feature is detected within the preceding or following 45 min. Lastly, the white dashed line visible around 00:00 highlights a limitation of our technique. Because our method processes daily keogram images, ROIs that extend across two keograms are split and identified in both images with incorrect velocities. This limitation explains the high-speed white solid line just before 00:00 showing a speed of 1,376 m/s and the subsequent line, which is discarded for unrealistically high velocity. Here, we specify that Figure 5 was obtained by the concatenation of two daily keograms available in the database.

Figure 6 presents the LSTID ROI identification process for the moderately disturbed period on 1–2 April 2024. During this event, the SymH index reached a minimum value of −44 nT, IE grew to 390 nT, and Hp30 peaked at 4.667. As anticipated, in contrast to the strong storm depicted in Figure 5, the weaker geomagnetic activity and

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Keogram belonging to the weak geomagnetic storm of 1–2 April 2023 with lines representing the RANdom SAmple Consensus fitted line. White solid lines are Region of Interest that satisfy the Large-scale traveling ionospheric disturbance conditions and are plotted along with their corresponding speed in m/s, which is shown in the white boxes.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** dTEC keogram for the 11–12 August 2024 storm over North America. The latitudinal extend goes from 28° to 48° North. The longitudinal range considered to generate the keogram is from 85° to 105° West. White lines represent the RANdom SAmple Consensus fitted line. White dashed lines indicate detected regions of interest (ROIs) that are not considered large-scale traveling ionospheric disturbance (LSTID) ROIs, while white solid lines satisfy the LSTID conditions and are plotted along with their corresponding speed in m/s shown in the white boxes.</figcaption></figure>

reduced energy injection at high latitudes resulted in fewer visible LSTIDs. This observation is particularly
significant, as previous statistical studies on LSTID for the European region (e.g., Borries et al., 2009) have
primarily focused on major storms, whereas moderately disturbed periods occur more frequently and are thus
critical for our study. In Figure 6, no white dashed lines are visible, indicating that all identified ROIs satisfied the
LSTID conditions. A solitary LSTID propagating at approximately 650 m/s appears around 15:00 and is followed
by intermittent LSTID activity between 18:00 and 02:00 on the morning of the 2nd of April. Notably, no major
misdetections are evident in this figure. However, two negative features with speeds of 903 and 365 m/s represent
false detections of LSTIDs. These false detections highlight the trade-off involved in setting the filtering pa-
rameters for the ROI identification technique. Increasing the strictness of the criteria reduces the number of false
detections but at the cost of more missed detections. We show the rationale behind the selection of the parameters
in Supporting Information S1. Since no ground truth for LSTIDs is available, it was not feasible to conduct a more
rigorous quantitative evaluation of the impact of filtering parameters on ROI identification performance in terms
of false positives and negatives.

After evaluating the detection algorithm's performance in two case studies, we now examine its results for the
same geomagnetic storm depicted in Figure 5 but across different geographical regions. Figure 7 presents the
detection output for the GNSS stations shown in the left panel of Figure 3. The keogram spans a 20° latitudinal
range, from 28° to 48° North, maintaining the same extent as EUREF but with different upper and lower limits.
Additionally, the longitudinal slice used to generate the keograms is 5° wider. A key challenge in these regions is
the limited availability of satellite observations. The CORS and GEONET networks do not provide BeiDou data,
and only a small number of stations collect Galileo signals. As a result, the total number of available line-of-sight
measurements is significantly reduced, even though the number of stations remains the same. To mitigate this
issue, we expanded the longitudinal extent for both North America and Japan, including more data in the keogram
generation process. This broader longitudinal coverage helps lower noise levels and ensures sufficient data
coverage across all latitudes and periods.

In the case of North America, the keogram reveals a high noise level, making LSTID detection more irregular and challenging. The first wave trains appear between 22:00 and 04:30, but increased ionospheric activity introduces overlapping disturbances beyond this period, complicating the detection. The lower latitudinal boundary also increases the likelihood of interhemispheric LSTIDs, which can interfere with image segmentation. Around 11:30, two ROI are filtered out (marked with dashed white lines) because the detected red feature exceeds the maximum velocity threshold (1,550 m/s). A similar issue arises between 19:00 and 20:00, where two more ROIs
<figure><img src="image_10.png" /><figcaption>**Figure 8.** dTEC keogram for the 11–12 August 2024 storm over Japan. The latitudinal extension goes from 30° to 45° North. The longitudinal range considered to generate the keogram is from 125° to 145° East. White lines represent the RANdom SAmple Consensus fitted line. White dashed lines indicate detected regions of interest (ROIs) that are not considered Large-scale traveling ionospheric disturbance (LSTID) ROIs, while white solid lines satisfy the LSTID conditions and are plotted along with their corresponding speed in m/s shown in the white boxes.</figcaption></figure>

are discarded due to their speeds exceeding 2 km/s. Lastly, at 15:00, an ROI is missed by the identification process, causing the adjacent positive and negative ROIs to be discarded as well since they lack a properly paired opposite-amplitude feature with matching physical characteristics. Despite these challenges stemming from high noise levels and overlapping ionospheric disturbances, our algorithm maintains an acceptable level of performance.

Figure 8 presents the results for the Japanese region, where the keogram spans a 15° latitudinal range, narrower than the European and North America cases. Overall, when comparing the detection results shown in Figure 7 with those of Figure 8, it is possible to note how Figure 8 shows more coherent parallel structures and thus LSTID ROIs over Japan are more easily identified. However, two notable misdetections occur around 15:00 and 17:00. The first misdetection arises from underestimating the negative band's inclination due to noise at low latitudes, which prevents the required speed match with the preceding positive ROI (both fitted lines appear dashed). The second misdetection occurs due to an incorrect segmentation of the irregular region during the BGMM processing step. This results in an underestimated speed for the positive ROI, preventing it from being merged with the following correctly classified LSTID ROIs. Despite these minor misclassifications, the algorithm demonstrates strong detection capabilities, particularly in the Japanese sector, where noise levels are comparatively lower than in the North America case, even if the latitudinal extension is 5° smaller.

## 3.2. Statistical Analysis

Having demonstrated the technique's capability to reconstruct LSTIDs in two distinct scenarios, we now focus on its statistical characterization when applied to the entire database of keogram images.

As LSTIDs need time to propagate from the auroral zone to the area where we detect them, a time shift was introduced to the indices to better understand the relationship between geomagnetic parameters and LSTID ROI occurrences: 2 hr for Hp30 and SymH and 1 hr for IE. The time shift values were selected to maximize the detection rates.

Additionally, the hourly and monthly distributions of LSTID ROI occurrences are shown in Figure 9. Figure 9a shows that LSTIDs are notably more frequent around the equinoxes, with a preference for the autumn equinox. The distribution exhibits a clear bimodal pattern. This behavior can be attributed to the increased strength and frequency of geomagnetic storms during equinoxes (Gonzalez et al., 1994; Tsugawa et al., 2004), driven by the geomagnetic field's orientation relative to the solar wind. Moreover, the monthly distribution shows that LSTIDs are more likely to happen in autumn than in spring in Europe; this difference is partially explained by the higher
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Large-scale traveling ionospheric disturbance (LSTID) regions of interest (ROIs) parameters statistics. Panels (a) and (b) show the monthly and hourly distribution of LSTID Region of Interest occurrences. Moreover, (a) also shows the monthly mean of the Hp30 index in orange. Similarly (b) shows the hourly mean of the 1-hr shifted IE index in orange. Panel (c) shows the period distribution, with the mean and standard deviation on top of the panel. Panel (d) shows the mean and standard deviation for the LSTID ROIs speed.</figcaption></figure>

value of average Hp30 during autumn showed by the orange curve. However, some other phenomena must play a crucial role, as the difference in the number of occurrences for autumn and spring is not completely justified by the difference in Hp30.

Furthermore, as LSTIDs represent displaced plasma, their higher amplitude during elevated background plasma density periods facilitates the detection. A similar explanation applies to the hourly distribution in Figure 9b, which shows that LSTIDs are primarily detected between 12:00 UT and 01:00 UT. The orange line in Figure 9b represents the hourly mean of the one-hour shifted IE index, revealing similarities between the two distributions. LSTIDs are thought to be generated by plasma displacement caused by neutral waves. However, post-midnight discrepancies arise due to the reduced background plasma density, which limits the detectability of LSTIDs, even if electrojet current enhancements are still generating neutral waves. Turning to LSTID ROI periods, Figure 9c shows a mean period of 62 min with a standard deviation of 14 min, closely matching the results of Borries et al. (2009), who reported an average period of 59 min with a similar standard deviation. The speed distribution, presented in Figure 9d, has a mean of 629 m/s and a standard deviation of 215 m/s, aligning well with Borries et al. (2009), who reported values of 684 m/s and 195 m/s, respectively. It is worth noting, however, that the distributions for period and speed differ in one key aspect. While the mean period corresponds closely to the mode, the mode of the speed distribution is lower than its mean (525 m/s compared to 629 m/s). This discrepancy arises from the velocity estimation process, which relies on the slope of ROIs in the keogram. The technique is susceptible to overestimating the average velocity because small errors of different signs in the slope of quasi-vertical lines result in different errors in the calculated velocities. In particular, positive deviations produce larger errors due to the tangential relationship between slope and speed, which can lead to significant over-estimations of the calculated speed.

As LSTID generation is closely linked to Joule heating, a relationship between the LSTID wave parameters and the IMAGE IE index is expected, as it measures the disturbances in the horizontal component of the geomagnetic field at the ground around the auroral oval caused by the intensification of ionospheric currents responsible for
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Violin plot of the period (a), speed (b) and mean Large-scale traveling ionospheric disturbance (LSTID) Region of Interest amplitude (c) with respect to the one-hour shifted IE index. The black diamonds represent the median, the gray horizontal lines is the mean, the darker shaded area represents the interquartile range, and the shape of the half-violin represents the distribution shape. Every violin shows the statistics for the LSTID ROIs related to a 1-hr shifted IE corresponding to the value shown on the x-axis.</figcaption></figure>

Joule heating. Figure 10 illustrates how the LSTID ROI period, speed, and mean LSTID ROI amplitude distribution vary with the shifted IE index. Here, mean LSTID ROI amplitude refers to the average value of dTEC for the pixels belonging to the considered LSTID ROI. Figure 10a reveals an almost linear decrease in the LSTID period as IE increases. For an IE of 0–100 nT, the median period is approximately 67 min, but it decreases to around 54 min for IE values exceeding 1,000 nT. While the period exhibits a clear inverse relationship with IE, the relation to speed, shown in Figure 10b, is more complex. Speed retrieval, being the most error-prone aspect of the technique, has a higher variability, particularly for low IE values. Nevertheless, a slight upward trend in speed is apparent, particularly for IE values above 200 nT. The median speed rises from approximately 566 m/s at an IE of 300 nT to around 628 m/s for IE values greater than 1,000 nT. Finally, the median of the mean LSTID ROI amplitude, shown in Figure 10c, demonstrates a direct relationship with the IE index. The amplitude represents the average absolute dTEC value for pixels within each ROI. The amplitude increases from about 0.05 TECu for an IE of 0–100 nT to approximately 0.12 TECu for IE values exceeding 1,000 nT. This relationship is consistent with findings by Borries et al. (2009). In summary, a stronger energy injection into the auroral oval leads to LSTIDs with shorter wavelengths and higher amplitudes. These results suggest that higher IE values generate stronger atmospheric disturbances mapped to the ionosphere as faster, shorter-period, and higher-amplitude LSTIDs. This relation highlights the link between increased IE and intensified wind, temperature, and pressure disturbances in the thermosphere.

After examining the parameters of LSTID ROIs and their hourly and monthly dependencies, we focus on their occurrence and detection rate relative to shifted geomagnetic indices. Figure 11 presents the results, where the blue bars show the distribution of events with respect to Hp30 (Figure 11a), SymH (Figure 11b), and IE (Figure 11c). Since geomagnetic disturbances are less frequent with higher magnitudes, the raw occurrence distribution of LSTID ROIs alone does not provide a complete picture of their relationship to geomagnetic activity. To address this, the green squares in Figure 11 depict the probability distribution of geomagnetic indices, while the orange dots show the detection rate expressed as a percentage. The detection rate is the number of LSTID ROIs normalized by the corresponding probability distribution of the geomagnetic index values along the x-axis. Figure 11a indicates that the highest occurrence of LSTID ROIs (blue histogram) corresponds to Hp30 values of 3. However, the probability distribution of Hp30 (green line) shows that the most frequent values are 1 and 2, with Hp30 values above five occurring each in less than 1% of cases. This implies that most LSTIDs are generated during mildly disturbed conditions (Hp30 between 3 and 5). Notably, the detection rate increases almost linearly with Hp30, reaching a maximum of 46% at Hp30 = 7 before decreasing. This behavior aligns with findings by Tsugawa et al. (2004), who observed a maximum occurrence rate of 75% for Kp = 9. The lower detection rate in our study is due to the finer temporal resolution of Hp30 (30 min vs. 3 hr for Kp), which increases
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Large-scale traveling ionospheric disturbance (LSTID) regions of interest (ROIs) detection occurrences and rates with respect to shifted geomagnetic parameters along with their probability distribution. Panels (a), (b) and (c) show in blue the number of LSTID ROIs detected when the corresponding shifted Hp30, SymH and IE index corresponded to the values shown on the x-axis. The orange dots and line represent the LSTID ROI detection rate in percentage normalized with respect to the distribution of occurrences of the corresponding shifted geomagnetic indices shown on the x-axis. The green squares and lines instead represent the percentage probability distribution of the geomagnetic indices shown on the x-axis.</figcaption></figure>

the likelihood of time bins without detectable LSTIDs. Additionally, the geographic coverage in our study introduces higher noise levels compared to studies where denser receiver networks were used, such as in Japan (Tsugawa et al., 2004).

For extreme storms, interference between equatorward and poleward LSTIDs at midlatitudes and noise from auroral oval displacement and other phenomena likely reduce detection rates. Another possible explanation is that during extreme storms, the assumption of a thin layer ionosphere might fail due to the increase in slab thickness (Pignalberi et al., 2022), causing an increased error in the projection of IPP and in the integration of disturbances along the line of sight. Figure 11b shows the distribution of LSTID ROIs with respect to the SymH index. Similar to Hp30, most LSTIDs are associated with quiet to mildly disturbed periods, with the highest occurrences for SymH values between -25 and -15 nT. The detection rate increases steadily as SymH decreases, peaking at 58% for SymH ≈ -120 nT. However, a saturation effect is observed as detection rates decline for extreme disturbances (SymH < -120 nT). An intriguing feature in Figure 11b is the spike in detection rate to 87.5% for SymH ≈ 60 nT. Despite the low number of occurrences, this anomaly suggests a unique generation mechanism for LSTIDs under such conditions, as enhanced auroral currents are unlikely for highly positive SymH values. Investigating the specific events responsible for this anomaly and their physical causes is beyond the scope of this study and will be addressed in future work. Figure 11c examines the relationship between LSTID ROIs and the IMAGE IE index, which directly measures the phenomena responsible for gravity wave excitation. Most LSTID occurrences correspond to IE values between 175 and 225 nT, representing mildly disturbed conditions (7% probability distribution of IE). As with Hp30 and SymH, the detection rate increases with IE, reaching 40% for values between 875 and 925 nT before saturating or declining to around 35% at higher IE values. In summary, the highest number of LSTID ROIs are associated with medium-low geomagnetic activity, as extreme storms are relatively rare. However, the detection probability increases dramatically during highly disturbed periods, reaching 46%. It is important to note that the detection rate in this study is based on 30-min time bins. Larger time bins, such as the 3-hr bins characteristic of the Kp index, would yield higher detection rates. For example, under 3-hr binning, the detection rate for a Kp of 7 would increase to approximately 75%.

## 4. Conclusions and Future Work

This study presents a novel methodology for automatically detecting LSTIDs derived from detrended TEC data. By leveraging GNSS networks and advanced image processing techniques, the method provides robust identification of LSTIDs, extracting key parameters such as propagation speed and period from keogram images. The approach successfully identifies LSTIDs during quiet and disturbed periods, as demonstrated through case studies and comprehensive statistical analysis for the European sector spanning 2021–2024. Moreover, we showed that the same technique can be applied to other geographical areas with comparable performances, requiring only
some parameter tuning to tackle differences in the GNSS networks. The key findings of our work are the
following:

* LSTIDs with periods falling in the 30–90 min range are more frequently observed during equinoxes and are most commonly detected during mildly disturbed geomagnetic conditions (e.g., Hp30 values between 3 and 5 or SymH values between –25 and –15 nT). The probability of detection increases significantly during periods of high geomagnetic activity, with detection rates reaching up to 50%.

* The results confirm a direct relationship between LSTID occurrence and geomagnetic indices such as Hp30, SymH, and IMAGE IE. Increased geomagnetic activity leads to LSTIDs characterized by shorter wavelengths, higher speeds, and greater amplitudes, consistent with known generation mechanisms tied to Joule heating and associated atmospheric disturbances.

* The results obtained using our methodology are consistent with previous studies conducted over different regions and periods using dTEC measurements (Borries et al., 2009; Ding et al., 2008; Tsugawa et al., 2004). In particular, our findings closely align with those of Borries et al. (2009), who reported an average period of 59 min and a phase velocity of 684 m/s, values that are in good agreement with our observed period of 62 min and phase speed of 629 m/s. Here, we acknowledge that by setting the Savitzky-Golay filter detrending window to 95 min, we tailored our user-definable algorithm to reproduce the detection capabilities of Borries et al. (2009). When compared with results from other geographic areas, some differences emerge. Tsugawa et al. (2004) identified a bimodal seasonal distribution, with a peak in LSTID occurrence during the spring equinox, whereas our analysis shows a maximum during the autumn equinox. Similarly, Ding et al. (2008) reported different average values for period and velocity. Their study also revealed a different dependence on local time. The discrepancies in speed and period across these studies can be attributed to several factors, including differences in the inclination and declination of the geomagnetic field, variations in geomagnetic latitude ranges, and, most importantly, differences in detection methodology. While our approach and that of Borries et al. (2009) rely on automated algorithms, the analyses by Tsugawa et al. (2004) and Ding et al. (2008) are based on manual visual inspection of dTEC maps. These methodological differences likely account for much of the variation in the reported LSTID parameters.

Despite the promising results, the study also identifies key limitations, such as the inability to fully resolve LSTIDs during periods of high noise or overlapping disturbances. One possible improvement is the implementation of a detection confidence metric to flag uncertain ROI. A human operator could review these flagged cases, offering a hybrid approach that balances reliability, speed, and objectivity. This scientist-in-the-loop approach would be particularly valuable in building a high-quality LSTID database for machine learning models or operational forecasting systems.

The obtained catalog of LSTID provides a foundation for machine learning-based forecasting models to predict LSTID occurrences based on geomagnetic indices and other environmental factors (Ventriglia et al., 2025). Such a possibility is encouraged by the demonstrated relation between the shifted geomagnetic indices and the LSTID detection rates and parameters. Moreover, by applying the required adjustments, the algorithm could be applied for near-real-time LSTID detection and could benefit operational systems that rely on ionospheric monitoring. Adjusting the detrending window size, the ROI definitions and detection thresholds we could apply the method to study poleward-propagating LSTIDs (negative speeds), dTEC responses to solar flares (extremely high speeds), or longer-period disturbances. Additionally, parameter customization allows users to tailor detection sensitivity to specific operational needs. Extending the temporal range will allow for broader statistical studies and more comprehensive method validation. Further analysis of unique phenomena, such as the anomalously high detection rates for specific SymH values, may provide new insights into LSTID generation mechanisms. This work lays the groundwork for improved understanding and monitoring of LSTIDs, contributing to developing predictive capabilities and mitigation strategies for ionospheric disturbances affecting critical technologies.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.
## Data Availability Statement

GNSS data were provided by the IAG (International Association of Geodesy) Regional Reference Frame sub-commission for Europe, EUREF (https://www.epncb.oma.be/). The Japan GNSS data are available from the GeoSpatial Authority of Japan (GSI; ftp://terras.gsi.go.jp/data; requires registration at https://www.gsi.go.jp/ENGLISH/geonet_english.html. More information: https://www.gsi.go.jp/ENGLISH/page_e30286.html). Hp30 values were gathered from the GFZ Postdam website (https://kp.gfz-potsdam.de/en/hp30-hp60/data) while SymH was retrieved from OMNIWeb (https://omniweb.gsfc.nasa.gov/form/omni_min.html). The IMAGE IE index was sourced via web services (https://space.fmi.fi/image/www/wget_instructions.php) provided by the Finnish Meteorological Institute at https://space.fmi.fi/image/www/index.php?page=home. The MATLAB library developed to retrieve GNSS vTEC from RINEX files is available at Guerra, 2025b. The python and MATLAB codes related to the generation of keograms, their analysis and the ROI pairing process are available at Guerra, 2025a.

## Acknowledgments

We thank the institutes who maintain the IMAGE Magnetometer Array: Tromsø Geophysical Observatory ofUiT the Arctic University of Norway (Norway), Finnish Meteorological Institute (Finland), Institute of Geophysics Polish Academy of Sciences (Poland), GFZ German Research Centre for Geosciences (Germany), Geological Survey of Sweden (Sweden), Swedish Institute of Space Physics (Sweden), Sodankylä Geophysical Observatory of the University of Oulu (Finland), DTU Technical University of Denmark (Denmark), and Science Institute of the University of Iceland (Iceland). The provisioning of data from AAL, GOT, HAS, NRA, VXJ, FKP, SIN, BOR, SCO and KUL is supported by the ESA contracts number 4000128139/19/D/CT as well as 4000138064/22/D/KS. This work has been funded by the Artificial Intelligence based Forecasting of Large-Scale Travelling Ionospheric Disturbances over Europe (ARISTIDE, https://progetti.ingv.it/it/aristide) project, funded by the Lazio Region, by means of Lazio Innova S.P.A (Management body). Marco Guerra's PhD is funded by an INGV project named "Rete Multiparametrica" under the topic "Space Weather: PECASUS." Vincenzo Ventriglia's current research fellowship is funded by the Travelling Ionospheric Disturbances Forecasting System (T-FORS, https://t-fors.eu/) project, which is funded by the European Union (GA-101081835). Open access publishing facilitated by Istituto Nazionale di Geofisica e Vulcanologia, as part of the Wiley - CRUI-CARE agreement.

## References

Afraimovich, E. L., Kosogorov, E. A., Leonovich, L. A., Palamartchouk, K. S., Perevalova, N. P., & Pirog, O. M. (2000). Determining parameters of large-scale traveling ionospheric disturbances of auroral origin using GPS-arrays. *Journal of Atmospheric and Solar-Terrestrial Physics*, 62(7), 553–565. https://doi.org/10.1016/S1364-6826(00)00011-0

Alfonsi, L., Spogli, L., Pezzopane, M., Romano, V., Zuccheretti, E., De Franceschi, G., et al. (2013). Comparative analysis of spread-F signature and GPS scintillation occurrences at Tucumán, Argentina. *Journal of Geophysical Research: Space Physics*, 118(7), 4483–4502. https://doi.org/10.1002/JGRA.50378

Altadill, D., Segarra, A., Blanch, E., Juan, J. M., Paznukhov, V. V., Buresova, D., et al. (2020). A method for real-time identification and tracking of traveling ionospheric disturbances using ionosonde data: First results. *Journal of Space Weather and Space Climate*, 10, 2. https://doi.org/10.1051/swsc/2019042

Astafyeva, E. (2019). Ionospheric detection of natural hazards. *Reviews of Geophysics*, 57(4), 1265–1288. https://doi.org/10.1029/2019RG000668

Banks, P. M. (1977). Observations of joule and particle heating in the auroral zone. *Journal of Atmospheric and Terrestrial Physics*, 39(2), 179–193. https://doi.org/10.1016/0021-9169(77)90112-X

Belehaki, A., Tsagouri, I., Altadill, D., Blanch, E., Borries, C., Buresova, D., et al. (2020). An overview of methodologies for real-time detection, characterisation and tracking of traveling ionospheric disturbances developed in the TechTIDE project. *Journal of Space Weather and Space Climate*, 10, 42. https://doi.org/10.1051/swsc/2020043

Borries, C., Jakowski, N., & Wilken, V. (2009). Storm induced large scale TIDs observed in GPS derived TEC. *Annales Geophysicae*, 27(4), 1605–1612. https://doi.org/10.5194/angeo-27-1605-2009

Boyre, B., Wood, A., Dorrian, G., Fallows, R. A., Themens, D., Mielich, J., et al. (2022). Lensing from small-scale travelling ionospheric disturbances observed using LOFAR. *Journal of Space Weather and Space Climate*, 12, 34. https://doi.org/10.1051/swsc/2022030

Bruyninx, C., Legrand, J., Fabian, A., & Pottiaux, E. (2019). GNSS metadata and data validation in the EUREF Permanent Network. *GPS Solutions*, 23(4), 106. https://doi.org/10.1007/s10291-019-0880-9

Cesaroni, C., Alfonsi, L., Pezzopane, M., Martinis, C., Baumgardner, J., Wroten, J., et al. (2017). The first use of coordinated ionospheric radio and optical observations over Italy: Convergence of high-and low-latitude storm-induced effects. *Journal of Geophysical Research: Space Physics*, 122(11), 11794–11806. https://doi.org/10.1002/2017JA024325

Chan, K. L., & Villard, O. G., Jr. (1962). Observation of large-scale traveling ionospheric disturbances by spaced-path high-frequency instantaneous-frequency measurements. *Journal of Geophysical Research*, 67(3), 973–988. https://doi.org/10.1029/JZ067i003p00973

Davis, M. J. (1971). On polar substorms as the source of large-scale traveling ionospheric disturbances. *Journal of Geophysical Research*, 76(19), 4525–4533. https://doi.org/10.1029/JA076i019p04525

Ding, F., Wan, W., Liu, L., Afraimovich, E. I., Voeykov, S. v., & Perevalova, N. p. (2008). A statistical study of large-scale traveling ionospheric disturbances observed by GPS TEC during major magnetic storms over the years 2003–2005. *Journal of Geophysical Research*, 113(A3). https://doi.org/10.1029/2008JA013037

Draper, N. R., & Smith, H. (2014). Applied regression analysis (3rd ed.). https://doi.org/10.1002/978118625590.ch15

Edelsbunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of a set of points in the plane. *IEEE Transactions on Information Theory*, 29(4), 551–559. https://doi.org/10.1109/TIT.1983.1056714

Fischler, M. A., & Bolles, R. C. (1981). Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. *Communications of the ACM*, 24(6), 381–395. https://doi.org/10.1145/358669.358692

Forbes, J. M., Palo, S. E., & Zhang, X. (2000). Variability of the ionosphere. *Journal of Atmospheric and Solar-Terrestrial Physics*, 62(8), 685–693. https://doi.org/10.1016/S1364-6826(00)00029-8

Francis, S. H. (1974). A theory of medium-scale traveling ionospheric disturbances. *Journal of Geophysical Research*, 79(34), 5245–5260. https://doi.org/10.1029/JA079i034p05245

Frissell, N. A., Kaeppeler, S. R., Sanchez, D. F., Perry, G. W., Engelke, W. D., Erickson, P. J., et al. (2022). First observations of large scale traveling ionospheric disturbances using automated amateur radio receiving networks. *Geophysical Research Letters*, 49(5), e2022GL097879. https://doi.org/10.1029/2022GL097879

Georges, T. M. (1968). HF Doppler studies of traveling ionospheric disturbances. *Journal of Atmospheric and Terrestrial Physics*, 30(5), 735–746. https://doi.org/10.1016/S0021-9169(68)80029-7

Gonzalez, W. D., Joselyn, J. A., Kamide, Y., Kroehl, H. W., Rostoker, G., Tsurutani, B. T., & Vasyliunas, V. M. (1994). What is a geomagnetic storm? *Journal of Geophysical Research*, 99(A4), 5771–5792. https://doi.org/10.1029/93JA02867

Guerra, M. (2025a). Automatic detection of large-scale travelling ionospheric disturbances using GNSS data and image processing techniques [Software]. Zenodo. https://doi.org/10.5281/zenodo.14988317

Guerra, M. (2025b). My_Ionosphere library [Software]. Zenodo. https://doi.org/10.5281/zenodo.15753170
Guerra, M., Cesaroni, C., Ravanelli, M., & Spogli, L. (2024). Travelling ionospheric disturbances detection: A statistical study of detrending techniques, induced period error and near real-time observables. *Journal of Space Weather and Space Climate, 14*, 17. https://doi.org/10.1051/swsc/2024017

Habarulema, J. B., Katamzi, Z. T., & McKinnell, L.-A. (2013). Estimating the propagation characteristics of large-scale traveling ionospheric disturbances using ground-based and satellite data. *Journal of Geophysical Research: Space Physics, 118*(12), 7768-7782. https://doi.org/10.1002/2013JA018997

Habarulema, J. B., Thaganyana, G. P., Katamzi-Joseph, Z. T., Yizengaw, E., Moldwin, M. B., & Ngwira, C. M. (2022). A statistical study of poleward traveling ionospheric disturbances over the African and American sectors during geomagnetic storms. *Journal of Geophysical Research: Space Physics, 127*(4), e2021JA030162. https://doi.org/10.1029/2021JA030162

Haralambous, H., Guerra, M., Chum, J., Verhulst, T. G. W., Barta, V., Altadill, D., et al. (2023). Multi-Instrument observations of various ionospheric disturbances caused by the 6th of February 2023 Turkey earthquake. *Journal of Geophysical Research: Space Physics, 128*(12), e2023JA031691. https://doi.org/10.1029/2023JA031691

Hines, C. O. (1960). Internal atmospheric gravity waves at ionospheric heights. *Canadian Journal of Physics, 38*(11), 1441-1481. https://doi.org/10.1139/p60-150

Hocke, K., & Schlegel, K. (1996). A review of atmospheric gravity waves and travelling ionospheric disturbances: 1982-1995. *Annales Geophysicae, 14*(9), 917-940. https://doi.org/10.1007/s00585-996-0917-6

Hunsucker, R. D. (1982). Atmospheric gravity waves generated in the high-latitude ionosphere: A review. *Reviews of Geophysics, 20*(2), 293-315. https://doi.org/10.1029/RG020i002p00293

Iyemori, T. (1990). Storm-time magnetospheric currents inferred from mid-latitude geomagnetic field variations. *Journal of Geomagnetism and Geoelectricity, 42*(11), 1249-1265. https://doi.org/10.5636/jgg.42.1249

Jing, N., & Hunsucker, R. D. (1993). A theoretical investigation of sources of large and medium scale atmospheric gravity waves in the auroral oval. *Journal of Atmospheric and Terrestrial Physics, 55*(13), 1667-1679. https://doi.org/10.1016/0021-9169(93)90171-T

Lloyd, S. (1982). Least squares quantization in PCM. *IEEE Transactions on Information Theory, 28*(2), 129-137. https://doi.org/10.1109/TIT.1982.1056489

MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. In *Proceedings of the fifth Berkeley symposium on mathematical statistics and probability* (Vol. 1, pp. 281-297). Retrieved from http://projecteuclid.org/euclid.bsmsp/1200512992

Makela, J. J., & Otsuka, Y. (2012). Overview of nighttime ionospheric instabilities at low- and mid-latitudes: Coupling aspects resulting in structuring at the mesoscale. *Space Science Reviews, 168*(1), 419-440. https://doi.org/10.1007/s11214-011-9816-6

Maletckii, B., & Astafyeva, E. (2022). Near-real-time analysis of the ionospheric response to the 15th of January 2022 Hunga Tonga-Hunga Ha'apai volcanic eruption. *Journal of Geophysical Research: Space Physics, 127*(10), e2022JA030735. https://doi.org/10.1029/2022JA030735

Maletckii, B., & Astafyeva, E. (2024). Near-real-time identification of the source of ionospheric disturbances. *Journal of Geophysical Research: Space Physics, 129*(11), e2024JA032664. https://doi.org/10.1029/2024JA032664

Mannucci, A. J., Wilson, B. D., Yuan, D. N., Ho, C. H., Lindqwister, U. J., & Runge, T. F. (1998). A global mapping technique for GPS-derived ionospheric total electron content measurements. *Radio Science, 33*(3), 565-582. https://doi.org/10.1029/97RS02707

Mendillo, M. (2006). Storms in the ionosphere: Patterns and processes for total electron content. *Reviews of Geophysics, 44*(4). https://doi.org/10.1029/2005RG000193

Munro, G. H. (1948). Short period changes in the F-region of the ionosphere. *Nature, 162*(4127), 886-887. https://doi.org/10.1038/162886a0

National Geodetic Survey. (2025). NOAA continuously operating reference stations (CORS) network. Retrieved from https://www.fisheries.noaa.gov/inport/item/39662

Oyama, S., & Watkins, B. J. (2012). Generation of atmospheric gravity waves in the polar thermosphere in response to auroral activity. *Space Science Reviews, 168*(1), 463-473. https://doi.org/10.1007/s11214-011-9847-z

Pignalberi, A., Pietrella, M., Pezzopane, M., Nava, B., & Cesaroni, C. (2022). The ionospheric equivalent slab thickness: A review supported by a global climatological study over two solar cycles. *Space Science Reviews, 218*(4), 37. https://doi.org/10.1007/s11214-022-00909-z

Rishbeth, H., & Mendillo, M. (2001). Patterns of F2-layer variability. *Journal of Atmospheric and Solar-Terrestrial Physics, 63*(15), 1661-1680. https://doi.org/10.1016/S1364-6826(01)00036-0

Roberts, S. J., Husmeier, D., Rezek, I., & Penny, W. (1998). Bayesian approaches to Gaussian mixture modeling. *IEEE Transactions on Pattern Analysis and Machine Intelligence, 20*(11), 1133-1142. https://doi.org/10.1109/34.730550

Sagiya, T. (2004). A decade of GEONET: 1994-2003 — the continuous GPS observation in Japan and its impact on earthquake studies. *Earth Planets and Space, 56*(8). https://doi.org/10.1186/BF03353077

Schafer, R. (2011). What is a Savitzky-Golay filter? [Lecture notes]. *IEEE Signal Processing Magazine, 28*(4), 111-117. https://doi.org/10.1109/MSP.2011.941097

Segarra, A., Altadill, D., Paula, V. de, & Navas-Portella, V. (2025). Climatology of large scale traveling ionospheric disturbances above Europe during the 2014-2023 period. *Journal of Space Weather and Space Climate*. https://doi.org/10.1051/swsc/2025024

Shiokawa, K., Lu, G., Otsuka, Y., Ogawa, T., Yamamoto, M., Nishitani, N., & Sato, N. (2007). Ground observation and AMIE-TIEGCM modeling of a storm-time traveling ionospheric disturbance. *Journal of Geophysical Research, 112*(A5), A05308. https://doi.org/10.1029/2006JA011772

Spogli, L., Piersanti, M., Cesaroni, C., Materassi, M., Cicone, A., Alfonsi, L., et al. (2019). Role of the external drivers in the occurrence of low-latitude ionospheric scintillation revealed by multi-scale analysis. *Journal of Space Weather and Space Climate, 9*, A35. https://doi.org/10.1051/swsc/2019032

Spogli, L., Sabbagh, D., Perrone, L., Scotto, C., & Cesaroni, C. (2024). Investigating the drivers of long-term trends in the upper atmosphere over Rome across four decades. *Journal of Space Weather and Space Climate, 15*, 8. https://doi.org/10.1051/swsc/2024040

Takahashi, H., Wrasse, C. M., Figueiredo, C. A. O. B., Barros, D., Paulino, I., Essien, P., et al. (2020). Equatorial plasma bubble occurrence under propagation of MSTID and MLT gravity waves. *Journal of Geophysical Research: Space Physics, 125*(9), e2019JA027566. https://doi.org/10.1029/2019JA027566

Tanskanen, E. I. (2009). A comprehensive high-throughput analysis of substorms observed by IMAGE magnetometer network: Years 1993-2003 examined. *Journal of Geophysical Research, 114*(A5), A05204. https://doi.org/10.1029/2008JA013682

Tsugawa, T., Saito, A., & Otsuka, Y. (2004). A statistical study of large-scale traveling ionospheric disturbances using the GPS network in Japan. *Journal of Geophysical Research, 109*(A6). https://doi.org/10.1029/2003JA010302

Vadas, S. L., & Crowley, G. (2010). Sources of the traveling ionospheric disturbances observed by the ionospheric TIDDBIT sounder near Wallops Island on the 30th of October 2007. *Journal of Geophysical Research, 115*(A7). https://doi.org/10.1029/2009JA015053

Vadas, S. L., & Liu, H.-l. (2009). Generation of large-scale gravity waves and neutral winds in the thermosphere from the dissipation of convectively generated gravity waves. *Journal of Geophysical Research, 114*(A10). https://doi.org/10.1029/2009JA014108
Ventriglia, V., Guerra, M., Cesaroni, C., Spogli, L., Altadill, D., Segarra, A., et al. (2025). An explainable machine learning model for large-scale travelling ionospheric disturbances forecasting. *Journal of Space Weather Space Climate, 15*, 25. https://doi.org/10.1051/swsc/2025020

Verhulst, T. G. W., Altadill, D., Barta, V., Belehaki, A., Burešová, D., Cesaroni, C., et al. (2022). Multi-instrument detection in Europe of ionospheric disturbances caused by the 15th of January 2022 eruption of the Hunga volcano. *Journal of Space Weather and Space Climate, 12*, 35. https://doi.org/10.1051/swsc/2022032

Yamazaki, Y., Matzka, J., Stolle, C., Kervalishvili, G., Rauberg, J., Bronkalla, O., et al. (2022). Geomagnetic activity index Hpo. *Geophysical Research Letters, 49*(10), e2022GL098860. https://doi.org/10.1029/2022GL098860

Zakharenkova, I., Astafyeva, E., & Cherniak, I. (2016). GPS and GLONASS observations of large-scale traveling ionospheric disturbances during the 2015 St. Patrick's Day storm. *Journal of Geophysical Research: Space Physics, 121*(12), 12138-12156. https://doi.org/10.1002/2016JA023332