© R. Ghidoni et al., Published by EDP Sciences 2025

1 Introduction

The LOw Frequency ARray (LOFAR) is an extensive network of radio telescopes primarily designed for astronomical observations (van Haarlem et al., 2013). However, in recent years, its capabilities have expanded beyond astronomy, generating considerable interest in its potential to probe the Earth’s ionosphere (Mevius et al., 2016; Fallows et al., 2020). Significant contributions have been made by LOFAR to ionospheric physics, particularly in identifying and analysing multi-scale ionospheric irregularities (Boyde et al., 2024) and, with a great level of detail, those embedded in Travelling Ionospheric Disturbances (TIDs) (Boyde et al., 2022; Dorrian et al., 2023).

Ionospheric irregularities are plasma density variations with respect to an ambient, uniform ionosphere, which vary on a large range of scale sizes, from the order of 1,000 kilometres down to the centimetre level (Wernik et al., 2003). The presence of these irregularities affects the propagation on trans-ionospheric radio waves through both diffraction and refraction effects and poses challenges to systems that depend on stable and precise signal transmission. Refraction is triggered by the full range of scales of the irregularities and causes deterministic fluctuations, while diffraction occurs when the scale size of the irregularities is below the Fresnel’s scale for every specific wavelength and is featured by a stochastic nature (McCaffrey & Jayachandran, 2019; Ghobadi et al., 2020). Both effects result in fluctuations in the signal amplitude and phase whenever received at ground. Those fluctuations are referred to as “scintillation” (Yeh & Liu, 1982).

One system affected by those ionospheric irregularities is the Global Navigation Satellite System (GNSS), which provides essential services in positioning, navigation, and timing across various sectors, including aviation, maritime navigation, surveying, and location-based services. LOFAR and GNSS operate across different probing frequencies and employ distinct observational geometries, allowing for a comprehensive examination of ionospheric irregularities on multiple scales. By studying these two systems together, we can uncover the intricate multiscale properties of ionospheric disturbances and gain deeper insights into their spatial and temporal evolution, as recently highlighted in the literature (Błaszkiewicz et al., 2021; Flisek et al., 2023) . Within this scope, we demonstrate how a plethora of phenomena can be highlighted simultaneously when LOFAR and GNSS are used in synergy to understand ionospheric response in the European sector. Thus, we focus on the geomagnetic storm that occurred on January 14th–15th, 2022, which was also characterised by the availability of LOFAR data for ionospheric studies using Cassiopeia A as a radio source. LOFAR is still a relatively new tool for ionospheric purposes, and the use of GNSS here adopted is also intended to step up our knowledge about the LOFAR realm for ionospheric studies.

The paper is structured as follows: Section 2 describes the geomagnetic conditions featuring the January 2022 storm, which triggered the ionospheric irregularities analysed in the paper; while Section 3 details the instruments, the related variables used and the adopted methods in the analysis. Section 4 provides the results and discusses them, with a particular focus on the determination of the velocity of the identified irregularities. Conclusions are given in Section 5, in which we summarise the results obtained in this paper, discuss their implications, address limitations, and suggest directions for future research.

2 The January 2022 storm

The arrival of an Interplanetary Coronal Mass Ejection (ICME) on January 14th, 2022, caused an increase in the main interplanetary magnetic field (IMF) and solar wind parameters measured at 1 AU (Fig. 1, panels a and b), resulting in a geomagnetic storm (Vanlommel, 2022; Pedersen et al., 2024) . The storm development was monitored using data collected by the Deep Space Climate Observatory (DSCOVR) spacecraft (Burt & Smith, 2012), located at the Lagrangian point L1. DSCOVR provides real-time data on solar wind (SW) parameters and the IMF (Loto’aniu et al., 2022). It collects plasma parameters (density, velocity, and temperature) at 1-min resolution, whereas the IMF measurements are at a 1-s resolution. The disturbance is also visible from the geomagnetic disturbance indices measured at ground (Fig. 1, Panels c, d and e). In fact, the Disturbance Storm Time (Dst, Panel e) reached a minimum of -91 nT on January 14th at around 23:00 UT. Figure 1d shows the three-hourly Kp index, in which the colour code indicates the NOAA Space Weather Scale for geomagnetic storms.¹

Figure 1 Total intensity |B| and Bz component of the IMF (Panel a), Solar Wind speed (Vsw) and pressure (Psw) (Panel b). Data in panel a and in panel b are measured at 1 AU. Other panels report the electrojet indicators IE, IL, IU in the European sector (Panel c), Kp (Panel d) and Dst (Panel e) indices during the storm. The colour code in panel d indicates the NOAA Space Weather Scale for geomagnetic storms. The blue shaded area in all panels represents the period between 16:00 UT of January, 14th and 08:30 UT of January 15th, the focus period of our analysis.

Related to the increment of the Kp index in the night between 14th and 15th, there was an intensification of the Auroral Electrojet in the European Sector. In Figure 1c, we observe the electrojet intensity indicators IU, IL, and IE crossing the International Monitor for Auroral Geomagnetic Effects (IMAGE) magnetometer network of the Finnish Meteorological Institute (Viljanen & Häkkinen, 1997; Kauristie et al., 2017). The increment in the electrojet indexes after 20:00 on January 15th is a trace of another geomagnetic storm that occurred concurrently with the main Hunga Tonga eruption (Verhulst et al., 2022), which occurred at around 4:00 UT on the same day. For this study, we focus on ionospheric variations between 16:00 UT on January 14th and 08:30 UT on January 15th, 2022, a period during which LOFAR data were available. The next section details the datasets and methodologies used to analyse these disturbances.

3 Data and methods

In this section, we provide details on the LOFAR telescope and GNSS receiver data. Specifically, we will mostly leverage the capability of the two instruments in detecting the presence of ionospheric irregularities. We exploit the diffraction effects on LOFAR received signals from Cassiopeia A and mostly refractive effects on GNSS signals. Additionally, we estimate the velocity of mid-latitude ionospheric irregularities as detected by LOFAR. A companion Defense Meteorological Satellite Program (DMSP) dataset of auroral radiance is introduced to characterise the position of the auroral oval boundaries.

3.1 LOFAR: Scintillation index (S4) and fresnel scale

LOFAR consists of an extensive network of 24 core stations in Exloo, the Netherlands, 14 remote stations in the Netherlands, and 16 international stations all over Europe. These are reported as red dots in Figure 2a, with a zoom on the core stations reported in Figure 2b. Although LOFAR operates within the theoretical frequency ranges of 10–90 MHz, for the Low Band Antennas (LBA), and 110–240 MHz, for the High Band Antennas (HBA), in this study, we only analyse data from the LBA. Specifically, we focus on the 25–63 MHz range, as data from the lower and upper boundaries of the LBA band are typically more contaminated by noise and instrumental artefacts. Thus, this selection ensures a more reliable analysis of ionospheric irregularities.

Figure 2 GNSS receivers and LOFAR stations map. Panel a) shows the locations of GNSS receivers across Europe, LOFAR international stations, and the CS002 core station, as representative of all CSs. Panel b) provides a zoomed-in view of the LOFAR CSs for greater detail.

For the purposes of our study, we analyse the normalised intensity of the signal coming from Cassiopeia A measured by the core stations and the international stations. A normalisation of the power data shown in Figure 3 is performed using a running median filter applied separately to each frequency channel. The window size for the median filter is 30 min, which suppresses slow variations in received signal power. The data have been processed to achieve a time resolution of ∼1 s and a frequency resolution of approximately ∼195 kHz. This reduction in resolution was implemented to decrease the volume of data relative to the raw observations initially recorded by the antennas, facilitating more efficient storage and processing while preserving essential scientific information.

Figure 3 Normalised power measured by the LBA antenna of CS002. In red, we highlighted a 2 MHz wide band around 60 MHz, the frequency we will use to investigate the variability induced by the ionosphere on the signal measured by LOFAR.

To quantify ionospheric diffraction effects, the standard metric is the amplitude scintillation index (S4) (Fremouw et al., 1978; Yeh & Liu, 1982). In this study, we compute this index from the dynamic spectra of the signal intensity measured by LOFAR:

$$S4=\sqrt{\frac{\langle I^2\rangle -{\langle I\rangle }^2}{{\langle I\rangle }^2}},$$(1)

in which I represents the normalised intensity of the signal. In this study, we measure the S4 on a sliding window of 60 s, considering the median normalised intensity, at each time step, in a 2 MHz-wide band around 60 MHz (red shade and dashed line in Fig. 3). The choice of 60 MHz is because it represents a frequency at which the signal is not strongly affected by ionospheric reflection of sub-horizon radio frequency interference (RFI) back toward the ground (van Haarlem et al., 2013).

The 30-min window for the median intensity normalisation does not fully eliminate long-term trends, introducing two primary sources of error. First, it may underestimate gradual scintillation increases, as the running median assumes a quasi-static reference level, reducing sensitivity to long-duration scintillation events. Consequently, when irregularities evolve on timescales comparable to or longer than the normalisation window, S4 values may be underestimated.

Second, it may overestimate short-lived fluctuations. Since the 30-min median filtering removes large-scale background variations, transient fluctuations can appear more pronounced in the S4 calculations, potentially inflating values, particularly under highly dynamic ionospheric conditions.

The scale dependence of S4 in LOFAR observations is primarily governed by the Fresnel’s scale, which determines the dominant irregularity size responsible for scintillation at a given frequency. Beyond frequency, factors such as propagation geometry and ionospheric altitude also play crucial roles in S4 calculation. This effect is particularly significant for low-elevation observations, where oblique propagation alters the effective scale of the irregularities. If the actual irregularity layer altitude differs from the modelled one, the expected Fresnel’s scale shifts, potentially leading to a discrepancy between expected and observed S4 variations.

Despite these limitations, S4 remains a valuable metric for characterising ionospheric scintillation. Its ability to quantify fluctuations in signal intensity makes it a crucial tool for studying small-scale plasma irregularities and their evolution over time.

For the LOFAR LBAs, the Fresnel’s frequency is of the order of a few tens of MHz while the distance to the irregularity layer varies with the elevation angle (θ).

For oblique incidence, the distance between the antenna and the irregularity layer is d = h/sin(θ), in which h is the altitude of the irregularity layer. Additionally, the cross-section of the irregularity, as seen by the antenna, becomes an ellipse, not a circle. Therefore, the length of the major axis of the irregularity can be approximated by $\frac{\sqrt{2\lambda h}}{\sin \left(\theta \right)}$, in which λ is the wavelength at which the observation is performed.

To conclude, as reported by Teunissen & Montenbruck (2017), the expression for the Fresnel’s scale for an oblique observation is:

$$d_F=\frac{\sqrt{2\lambda d}}{\sin \left(\theta \right)}=\frac{\sqrt{2\lambda h}}{\left(\sin \right(\theta ){)}^{\frac{3}{2}}}.$$(2)

To provide an estimate of the scale sizes involved in the LOFAR scintillation, we assume a single irregularity layer located at an altitude of 350 km. Such an assumption agrees with a reasonable value of the F-layer altitude, at which ionospheric irregularities affecting the mid-latitude ionosphere mainly occur.

Thus, to better corroborate the scale sizes causing LOFAR scintillation, Figure 4 reports the Fresnel’s scale at a frequency of 60 MHz, considering the Ionospheric Pierce Point (IPP) altitude at 350 km as seen by the LOFAR international stations, and the core station CS002LBA, as a function of time.

Figure 4 The figure represents the LOFAR Fresnel’s scale at a frequency of 60 MHz, considering the altitude of a single ionospheric layer at 350 km.

3.2 Irregularity velocity estimate with LOFAR

Using LOFAR core stations, it is possible to estimate the velocity of ionospheric irregularities passing across the field of view of the antennas (Fallows et al., 2020). Specifically, from the normalised power measured by each antenna, the velocity can be expressed in the following form:

$$\overrightarrow{V_{\mathrm{irr}}}=\frac{{\mathbf{B}}^{†}\vec{T}}{\left|\right|{\mathbf{B}}^{†}\vec{T}|{|}^2},$$(3)

where B^† represents the inverse of the matrix B containing the distances between the IPPs of the LOFAR core stations, while $\vec{T}$ is a vector containing the time lag between the signals of every pair of antennas. The element (i, j, t) of matrix B stores, at each time-step t, the East–West and North–South components of the distance between the IPP where the disturbance was seen by station i and the IPP where it was seen by station j. After constructing the matrix B, it is then transformed into a vectorized version, by stacking its columns vertically. This is done to obtain, as a final result, a vectorised version of $\overrightarrow{V_{\mathrm{irr}}}$.

The algorithm is a two-step process. In the first part, signals from each baseline are analysed to find the best cross-correlation times of signal power time series. This is done by analysing, for each pair of core stations, the time lag between the signals received by the two stations, identifying the lag maximising the cross-correlation. In our case, we consider Fast Fourier Transform (FFT) convolutions performed on 6-min windows for the first station and a one- minute window for the second station. This was done every 30 s. These parameters follow from a good balance between the accuracy of the results and the computation time. Using those parameters, the FFT convolution method results are robust against the presence of fast, unphysical data points observed in the signal, which appear as sudden drops or sharp increments over short time intervals and across all frequencies. In this first part, the algorithm thus provides the time lags and the cross-correlation for each baseline every 30 s for the entire duration of the analysis. To mitigate the impact of the data resolution, we apply a method to obtain more precise lag estimates. For each time step, we extract the entire cross-correlation matrix and identify the peak cross-correlation and a window of 11 s around it. To refine the lag estimation, we apply a cubic spline interpolation on the selected cross-correlation values and determine the maximum using the derivative of the spline. This approach allows for more precise lag estimates despite the one-second time integration of the data. This method, however, does not completely remove the errors related to the lag estimations. This is true in particular for small lags, where even a small absolute lag error can result in a large relative error in velocity estimation.

In the second part, the algorithm minimises the function ${\left|\right|\mathit{B} \frac{\overrightarrow{V_{\mathrm{irr}}}}{\overrightarrow{V_{\mathrm{irr}}^2}}-\vec{T}\left|\right|}^2$, at 30-s intervals for the variable $\overrightarrow{V_{\mathrm{irr}}}$. A linear least squares fit seeks the value of the velocity for which the time lags are as close as possible to the product of the normalised velocity and the length of each baseline.

An alternative method for estimating the velocity of ionospheric irregularities involves the Fresnel’s frequency and the Fresnel’s zone. The absolute value of the relative velocity between the irregularity velocity and the ionospheric pierce point (IPP) velocity can be estimated via (Spogli et al., 2021; Forte & Radicella, 2002):

$$v_F=d_F\cdot f_F,$$(4)

where:

• v_F is the Fresnel’s velocity,

• d_F is the Fresnel’s zone of the selected frequency (Eq. (2)),

• f_F is the Fresnel’s frequency, calculated as the roll-off frequency of the Fourier transform of the dynamic power spectra (P) every 5 min of data (Eq. (5)).

It is important to remember that this method estimates the relative speed between plasma and the IPP. While the IPP velocity is significantly lower than the plasma velocity, this calculation yields only the speed, not the velocity vector.

The Fresnel’s frequency is determined as the roll-off (f_R) frequency of the Power Spectral Density (PSD) based on the Fourier Transform, identified through a fitting procedure at five-minute intervals for each station.

To determine the roll-off frequency, a Fourier transform is applied to a five-minute time series sampled at 60 MHz. Then a piecewise model is fitted to the spectrum to identify the transition point corresponding to the roll-off frequency. The power spectrum (P) is ideally described as a function of frequency f:

$$\log \left(P\right)=\{\begin{array}{ll}C,& \log \left(f\right)\le \log \left(f_R\right)\\ s\cdot \log \left(f\right)+b,& \log \left(f\right)>\log \left(f_R\right)\end{array},$$(5)

where, C represents the constant power level at low frequencies, while the high-frequency region follows a power-law decay with slope s and intercept b. The roll-off frequency is indicated by f_R. A least-squares fitting method is applied to optimise the parameters. The roll-off frequency is identified as the point where the two segments of the piecewise fit intersect. Confidence intervals at 95% for the roll-off frequency and the slope are estimated to account for uncertainties in the fit.

Figure 5 reports an example of this process for data recorded between 23:22:30 and 23:27:30 UT on January 14th, where the red line represents the fitted function. The blue line represents the log-log spectrum. The roll-off frequency is indicated by the dashed blue line while its confidence interval is represented by the light blue box around it. In the gray box, we report the best fit values and their standard deviation for the constant, the slope and the intercept of equation (5).

Figure 5 Example of roll-off frequency estimation from the PSD based on FFT of 5 min of data between 23:22:30 UT and 23:27:30 UT of January 14th. The blue line represents the log-log spectrum, while the red line shows the fitted function following equation (5). The roll-off frequency is indicated by the dashed blue line, while its standard deviation is represented by the light blue box around it. In the gray box, we report the best fit values and their standard deviation for the constant, the slope and the intercept of equation (5).

This process is repeated for each core station, and the resulting velocities are averaged to obtain the mean Fresnel’s velocity. Identifying the Fresnel’s frequency with the roll-off frequency of the amplitude spectrum introduces some caveats. Within the weak scattering theory, it has been proposed that the amplitude power spectrum bends over when Fresnel-scale structures are present (Yeh & Liu, 1982). In this case, the roll-off frequency is assumed to approximate the Fresnel’s frequency. However, recent findings (Song et al., 2025) highlight significant differences between the roll-off and Fresnel’s frequencies, suggesting that the phase spectrum method may provide a more accurate determination, particularly in high-latitude scintillation studies. Bearing in mind these limitations, we retrieved the velocity from Fresnel’s frequency estimates, which can be compared with velocity estimates obtained via the cross-correlation method.

While plasma density irregularities in both the ionosphere (with moderate-to-high relative drift) and the inner heliosphere (with moderate-to-high drift) can produce these Fresnel’s frequency values, the variations observed between LOFAR stations suggest that the analysed signals are most likely of ionospheric origin (Forte et al., 2022).

3.3 GNSS: ROTI and travelling ionospheric disturbances

Besides being widely used for positioning applications, GNSS receivers serve also as a precious tool for providing ionospheric information. We leverage on dense networks of geodetic GNSS receivers located all around the European sector. Specifically, blue dots in Figure 2 report the position of all the considered GNSS receivers from various networks. All the considered GNSS receivers are dual-frequency so they can be exploited to retrieve the geometry-free linear combination (GFLC) and, conversely, estimate the Rate of Total Electron Content (TEC) change Index (ROTI) (Pi et al., 1997), as measured from:

$$\mathrm{ROTI}=\sqrt{\langle \mathrm{RO}{\mathrm{T}}^2\rangle -{\langle \mathrm{ROT}\rangle }^2}.$$(6)

It is defined as the standard deviation of the change of Rate of TEC (ROT) over a given time window and it is measured in TECU/min. ROT is defined as the variation of TEC in time (dTEC/dt) between two consecutive epochs. In this study, we consider a time window of 5 min, to measure the standard deviation of 10 points, considering the TEC data are extracted at 30-s intervals (see Spogli et al. (2023) and references therein).

As dTEC/dt is used in the ROTI calculation, and receiver and satellite biases can be considered constant between two epochs, therefore ROTI measurements do not face the usual TEC calibration issues (Ciraolo et al., 2007; Cesaroni et al., 2021; Tornatore et al., 2021). As a consequence, estimations derived from GFLC can be properly exploited to provide reliable ROTI. It’s important to clarify that ROTI is a general index for ionospheric activity and is not exclusive to GNSS; it can be calculated from any instrument measuring TEC. However, for this analysis, the dense GNSS network provides the ideal data source for its derivation. ROTI mainly accounts for the refraction effects induced by the ionosphere on the signal, and its capability of probing given scale sizes is linked to the relative velocity of satellite IPPs and the ionosphere and by the Nyquist frequency of the GNSS observations. Usually, ROTI values larger than the noise floor suggest the presence of ionospheric irregularities at scale lengths of a few kilometres, if the drift velocity of the irregularities in the direction perpendicular to the GNSS ray path is of the order of 100 m/s (Ma & Maruyama, 2006). Several factors can introduce errors in ROTI estimates. Specifically, these are mainly caused by the occurrence of cycle slip detection, carrier phase thermal noise, receiver oscillator noise, and multipath (see Zhang et al. (2025) and references therein). In this study, we firstly mitigate these single-receiver errors excluding observations below a 20° elevation angle, to minimise the multipath effect. We also account for cycle-slips in the GNSS phase signal in the following way: We compute the difference between consecutive phase observations and then we determine the standard deviation of these differences to estimate the normal variation. A relative threshold is set at 10 times the standard deviation, while an absolute threshold is fixed at 100 TECU. Any phase difference exceeding either of these thresholds is marked as a sudden jump. When a jump is detected, the GFLC data after the jump is flagged and removed from the analysis. However, this correction process leads to a reduction in the number of available samples for the ROTI calculation. As ROTI is a statistical measure, it requires a sufficiently large dataset for reliable estimation. Thus, we discard ROTI values derived from epochs containing fewer than 10 valid samples, considering that this represents the maximum number of available samples within the 5-min chosen window (considering the 30-s interval of the GNSS observables). In the end, we also remove the results of the ROTI coming from a specific satellite-reciever couple that presents unphysical behaviour, like constant ROTI equals zero. Due to the huge number of GNSS receivers across Europe the number of usable ROTI values after this screaning is still great enough to properly cover the ionospheric region observed by LOFAR.

Additionally, to make use of the time mean in the ROTI definition, we must assume an ergodic process (Carrano et al., 2019), i.e., that the rate of change of TEC fluctuations is a zero-mean random process and no additional detrending is needed on the TEC measurements from which the given ROTI originates. Additionally, when considering ROTI between different receivers, errors can arise from variations in satellite tracking algorithms, tracking loop designs, firmware, and multipath mitigation techniques used in different receiver models (see Sui et al. (2024) and references therein). The diverse tracking techniques employed by various GNSS receiver types have been highlighted as likely to be the major contributor to the inconsistency between ROTI derived from multiconstellation measurements (Liu et al., 2019)]. In the work by Liu et al. (2019), who compared four GNSS receiver types (Javad, Leica, Septentrio, and Trimble) at high-, low- and mid-latitudes, they reported differences between 1-s ROTI up to 6 TECU/min, which can be reduced to 0.2–0.5 TECU/min when adopting 15- and 2-Hz noise bandwidths in the tracking loop. Bearing this in mind, the use of ROTI in our work has a twofold aim:

• using the larger network of geodetic receivers providing 30 s data to increase the granularity of the ROTI information,

• identifying the presence of medium-scale (km to tens of km scale) irregularities from the GNSS perspective to complement and understand the information retrieved from LOFAR.

Bearing this in mind and for the purposes of our work, we need to identify the presence of medium scale irregularities with ROTI enhancements, i.e., values of ROTI that are above a threshold of 0.2 TECU/min. This value was chosen as it effectively accounts for the overall noise level of the observation both prior to 18:00 UT on January 14th and subsequent to 06:00 UT on January 15th, as visible in Figure 6. This selection allowed for a better distinction between ROTI increments attributable to the storm and those from other effects. Additionally, according to Cherniak et al. (2015), the value of 0.2 TECU/min is a reasonable threshold for differentiating between the absence of significant TEC gradients and the presence of the trough walls of the high-latitude ionosphere.

Figure 6 ROTI and S4, measured at 60 MHz, time variability depending on latitude. All the points are projected at an IPP altitude of 350 km. Only data between 45° N and 65° N are shown. This is the same observational range of LOFAR data during this storm.

The same GNSS raw measurements are also used to reconstruct possible wave-like structures that propagate as a result of the geospace forcing. Specifically, we make use of the vertical TEC (vTEC) estimates retrieved by applying the vTEC_Ne technique, described in Guerra et al. (2024), to the GFLC data used for the ROTI calculation. This conversion, to a vertical path, allows to create consistent global maps and compare ionospheric conditions over time and location, regardless of the satellite’s position.

Once the vTEC is derived, a third-order Savitsky–Golay filter (Schafer, 2011) on a 90-min window is used to remove all the low-frequency components together with a moving average smoothing on a window of 15 min to remove all the high-frequency ones. This filtering will highlight the expected signatures of large-scale Travelling Ionospheric Disturbances (TIDs), which are normally defined by periods between 45 and 90 min. To investigate the occurrence of TIDs, we exploited the widely used keograms, a latitude-time plot of δTEC variations to the mean value. The longitudinal sector of interest was between 0° E and 20° E. To reduce the noise, we smoothed the δTEC measurements by averaging over bins of 2.5 min in time and 5 km in space. In a keogram, TIDs typically appear as parallel diagonal bands of different colours (in our case, blue/red for negative/positive δTEC). The slope of these parallel bands provides insight into the meridional velocity and north-south direction of the TIDs, with steeper slopes indicating faster-moving disturbances and shallower slopes corresponding to slower-moving ones. The amplitude of the δTEC variations, as seen in the keogram, reveals the strength of the TIDs, with larger amplitude fluctuations indicating more significant perturbations in the ionosphere’s electron content.

3.4 Complementary satellite data

In our study, additional information is provided by the auroral oval reconstruction, modelled by leveraging on the data from the Special Sensor Ultraviolet Spectrographic Imager (SSUSI) on board the DMSP F17 and F18 satellites, operated by the U.S. Air Force (Gussenhoven et al., 1981; Lane et al., 2015) . The SSUSI instruments are highly sophisticated ultraviolet spectrographic imagers designed to observe and measure the Earth’s upper atmosphere and ionosphere by capturing ultraviolet emissions from the auroras. The use of auroral oval boundary location, as retrieved from the DMSP auroral radiance data, has been proven to be effective in decoupling the upper atmosphere phenomena occurring in the auroral and sub-auroral regions (Griffin et al., 2012; D’Angelo et al., 2021; Papini et al., 2023; Quattrociocchi et al., 2024). . Additionally, we leverage the in-situ informations provided by the ESA’s Swarm satellite mission (Friis-Christensen et al., 2008) . It consists of a constellation of three satellites named Alpha (A), Bravo (B), and Charlie (C), which operate in quasi-Sun-synchronous near-polar low-Earth orbits at different altitudes. Swarm A and C form the lower pair of satellites, flying in proximity throughout the mission, while Swarm B is in a higher orbit. During the Mother’s Day storm, Swarm A and C were orbiting at approximately 480 km above the Earth’s surface, while Swarm B was at a higher altitude of about 510 km. These satellites are equipped with various instruments. One of these is constituted by a pair of Langmuir Probes (LPs), which measure the plasma density and temperature at a 2 Hz rate, revealing interesting features of the topside ionosphere the topside ionosphere (Wood et al., 2022).

4 Results and discussion

Figure 6 shows the comparison between the time profiles of ROTI, from the GNSS stations and of S4 AT 60 MHz from the LOFAR stations. The S4 data are measured considering the median power of all the frequencies in a 2 MHz wide band around the centre for every temporal interval. The colour represents the value of the latitude of the IPP projected at an altitude of 350 km, which is a reasonable altitude for the ionospheric F layer and, thus, in this figure, we identify that layer as the irregularity layer. Discussion about this assumption is provided later in this section. The GNSS data were dismissed if at an elevation of less than 20°, to minimise the effect of the multipath, which can mimic ionospheric effects (Romano et al., 2013; D’Angelo et al., 2015; Li et al., 2022). ROTI data only refer to the IPP latitude range between 45° and 65° N, to match the same range of LOFAR data. It is worth noticing, however, how the GNSS receivers cover a larger latitudinal range, providing useful insights into the irregularities moving from the Arctic regions toward the equator. In the S4 plot, it is important to note that the brief increment in the S4 value around 18:00 UT of January 14th corresponds to a few-second increase in a single German station (DE603) that affected all frequencies. Therefore, we do not consider this increment to be due to ionospheric effects and have opted to exclude a detailed analysis of it from the manuscript.

To investigate the ionospheric sectors involved in the formation of the irregularities, we investigate geographical maps of ROTI and S4 at 5-min intervals. An animation covering the full storm development is provided in the Supplementary material A. Specifically, the animation shows the period between 16:00 UT of January 14th, 2022, to 6:00 UT of January 15th, 2022, with a time step of 5 min. Blue dots show the LOFAR data. We show the data from all the international stations and only the core station CS002. S4 data from other core stations are not reported as they do not differ from CS002, due to their proximity. The blue shade intensity represents the S4 magnitude for the LOFAR selected frequency. If the S4 measured by a station is below the threshold of 0.05, the data is represented with a red”x”. On the other hand, the yellow/black points represent the ROTI value from GNSS above 0.2 TECU/min. Values below this threshold are indicated with a small grey dot. The red lines represent the edges of the modelled auroral oval, obtained from the closest (in time) imaging data from the DMSP satellite. Both LOFAR and GNSS data are projected at an IPP of 350 km, if outside the auroral oval, or 120 km, if inside of it. This is because, the maximum height of the particle precipitation causing ionospheric irregularities within the oval is expected to be located at around 120 km, while outside, F-layer irregularities are more likely to occur. The validity of the assumption of putting the IPP of LOFAR at 350 km below the auroral oval is confirmed by the Juliusruh ionosonde observations, while the assumption of 120 km as the height within the oval are a reasonable value for the height of the irregularity has been verified by using the ionosonde located in Sodankylä, Northern Finland.

The title of each map reports the time of the ROTI data. Concerning the LOFAR data, we show the mean S4 value of the data collected by each station in the five minutes before the time shown in the title. We report the time of the nearest auroral oval position estimation, labelled as “Auroral time”. For the sake of clarity, we also report in Figure 7 the map at 20:05 UT of January 14th, when we begin observing irregularities measured by both instruments.

Figure 7 Map of the ionospheric irregularities measured around 20:05 UT of January14th over the European sector at latitudes between 45° and 80° N and longitudes between −15° and 25° E. The red lines represent the boundaries of the auroral oval modelled by the DMSP. The ROTI data are represented in the yellow–black scale if over the threshold of 0.2 TECU/min or as simple gray dots if under it. The LOFAR S4 follow logic with a threshold of 0.05 and a scale green/blue and red “x”.

In Figure 8, we report the results at 22:00 UT of January 14th, where we have stronger irregularities. It is interesting to notice how inside the auroral oval we can clearly observe irregularities measured by both GNSS and LOFAR. This was expected, considering the auroral oval is the region of main ionospheric activity during a geomagnetic storm. The same is true near the lower boundary of the auroral oval, where we are observing the effects of the ionospheric trough, in agreement with what found by Cherniak and co-authors found during the 2015 St.Patrick’s Day Storm (Cherniak et al., 2015) . The presence of a steep ionospheric trough is also visible in the electron density data collected by ESA Swarm Bravo satellite from 23:07 UT to 23:16 UT on January 14th. Plot a) of Figure 9 shows the map of the ionospheric irregularities measured around 23:10 UT of January14th over the European sector at latitudes between 45° and 80° N and longitudes between −15° and 25° E. The dark green line represents the trajectory of the Swarm Bravo satellite passing over Europe, from South to North, between 23:07 UT and 23:16 UT of January 14th. The light green line indicates the in-situ plasma density of Swarm Bravo, with a higher longitudinal positive deviation from the ground track, the higher the plasma density is. The other variables follow the same colour description as Figures 7 and 8.

Figure 8 Map of the ionospheric irregularities measured around 22:00 UT of January14th over the European sector at latitudes between 45° and 80° N and longitudes between −15° and 25° E. The red lines represent the boundaries of the auroral oval modelled by the DMSP. The ROTI data are represented in the yellow–black scale if over the threshold of 0.2 TECU/min or as simple gray dots if under it. The LOFAR S4 follow the same logic with a threshold of 0.05 and a scale green/blue and red “x”.

Figure 9 Plot a) shows the map of the ionospheric irregularities measured around 23:10 UT of January14th over the European sector at latitudes between 45° and 80° N and longitudes between −15° and 25° E. The dark green line represents the trajectory of the Swarm Bravo satellite passing over Europe, from South to North, between 23:07 UT and 23:16 UT on January 14th. The light green line indicates the in-situ plasma density of Swarm Bravo, with a higher longitudinal positive deviation from the ground track, the higher the plasma density is. The other variables follow the same colour description as Figures 7 and 8. Plot b) shows the in-situ electron density measured by Swarm Bravo (red line) between 23:07 UT and 23:16 UT of January 14th over the European sector. The blue dots represent the ROTI values measured by the GNSS. The yellow band represents the enhancement of ROTI between 0.2 and 0.5 TECU/min due to the presence of the poleward boundary of the ionospheric trough.

Plot b) shows the in-situ electron density measured by Swarm Bravo (red line) between 23:07 UT and 23:16 UT of January 14th over the European sector. The blue dots represent the ROTI values measured by the GNSS. The yellow band represents the enhancement of ROTI between 0.2 and 0.5 TECU/min, from which it is possible to see how values in such a range nicely correspond to the poleward boundary of the ionospheric trough.

However, the most intriguing region in Figure 8 is the one below the auroral oval, in which we have the presence of irregularities measured by LOFAR, but not by GNSS. This is probably due to the scales of the irregularities that are active for LOFAR S4 and not effective in increasing the GNSS ROTI values.

To investigate the ROTI/S4 mismatching at sub-auroral latitudes, we analyse the normalised power data from the LOFAR core stations. An example of this analysis is reported, in the period between 21:20:00 UT and 21:24:59 UT on January 14th, in the Supplementary material B. During this period, ionospheric structuring increases, due to the response to the storm. In the animation in Supplementary material B, IPPs of all the core stations are reported, considering an altitude of 350 km. The colour of the dots represents the normalised power measured by each station at 60 MHz. During this period, we can observe a wave-like perturbation crossing all IPPs. An example of the reconstructed velocity of the ionospheric irregularities (Eq. (3)) passing across the core stations is reported in Figure 10. Specifically, Figure 10 shows the normalised power measured by the LOFAR core stations projected at an IPP of 350 km between four consecutive instants, i.e., from 21:23:02 UT to 21:23:05 UT on January 14th. The grey arrow indicates the reconstructed direction of the velocity of the perturbation using the cross-correlation method. The velocity is measured every 30 s so, for every frame, we report the nearest available velocity, for which the time is in the plot title. The trajectory of the LOFAR observations (blue line) is shown in the upper left corner, with the position at the considered time indicated by a red dot. We also report the amplitude and direction of this measured velocity.

Figure 10 Normalised power measured by the LOFAR core stations projected at an IPP of 350 km for four consecutive instants, i.e., 21:23:02 UT to 21:23:05 UT of January 14th. The grey arrow reports the direction of the velocity of the perturbation. We report the nearest velocity available, for which the time is in the plot title. The trajectory of the LOFAR observations (blue line) is shown in the upper left corner, with the position at the specified time indicated by a red dot.

To measure the velocity of the irregularities seen by the LOFAR core stations, we first used Fresnel’s method (Eq. (4)). We then compare the results to the velocities obtained using the cross-correlation method (Eq. (3)). Both algorithms refer to LOFAR intensity data at 60 MHz.

Figure 11 shows the velocities obtained by the two methods. The blue line represents the mean velocity of all the core stations obtained using the Fresnel’s algorithm (Eq. (4)) while the black line represents the velocity obtained by the cross-correlation algorithm (Eq. (3)). The colour of the dots represent the mean cross-correlation coming from every baseline of the LOFAR core stations at an interval of 30 s. The similarity between those two algorithms, for the most part of the entire observation, supports the robustness of our results. The main discrepancies between the two methods occur at the beginning of the observation (between 16:00 UT and 18:00 UT of January 14th), in regions where the cross-correlation is low, which in turn reduces the capability of the method to provide reliable values of the velocity. In such cases, the cross-correlation method is less accurate, and the velocity estimation becomes less certain. Additionally, within the same time window, the S4 values are also low, which further limits velocity estimates from the PSD due to the absence of Fresnel’s filtering mechanism acting on the recorded signals. Another time window, characterised by significant differences in the retrieved velocity with the two methods, is between 20:00 UT and 21:00 UT of January 14th. One possible reason is that this period corresponds to the highest velocity values (≥800 m/s), leading to a greater likelihood of errors associated with the previously described spline-based method for lag estimation in the cross-correlation algorithm. However, the method appears to properly reconstruct the lag, as evidenced by the very high cross-correlation values, shown in red (indicating cross-correlation close to 1) in Figure 11. However, it is important to consider that the Fresnel’s method, using 1-s resolution data, is limited to a finite range of possible roll-off frequency. The greater one is e ^−0.5 Hz. This, combined with the Fresnel’s scale of around 2 km, between 20:00 and 21:00 UT (Fig. 4 in the manuscript), results in a maximum velocity reconstruction from LOFAR data of around 1,200 m/s. These results are below the velocities measured using the cross-correlation algorithm. For this reason, we have decided, for the data between 20:00 UT and 21:00 UT of January 14th, to consider the velocities obtained from the cross-correlation algorithm as more reliable than the velocities obtained from the Fresnel’s algorithm.

Figure 11 Velocity estimates during the storm by using cross-correlation (CC) and velocity from Fresnel’s Frequency identification. The blue line represents the mean velocity of all the core stations obtained using the Fresnel’s algorithm (Eq. (4)), while the black line represents the velocity obtained by the cross-correlation algorithm (Eq. (3)). Both algorithms refer to LOFAR intensity data at 60 MHz.

Additionally, velocity estimates derived from the cross-correlation algorithm allow us to compute both the North–South and East–West components of the velocity, providing essential directional information for understanding the ionospheric flow. In contrast, velocity retrieval based on Fresnel’s frequency only measures the magnitude of the relative velocity between the ionospheric plasma and IPPs.

Figure 12 shows the magnitude and the East–West (E–W) and North–South (N–S) components of the velocity measured by the LOFAR core stations. The colour code follows from the mean cross-correlation evaluated from all the baselines every 30 s. The cross-correlation provides a measure of the reliability of the velocity estimates: The larger the cross-correlation, the more precise is the measurement. To this scope, we consider the mean cross-correlation of all baselines, evaluated every 30 s. A negative sign in the velocity components indicates a wave propagating from East to West and from North to South. The speed is of the order of some hundreds of m/s, with slightly higher values for the E–W direction. However, before about 20:00 UT and after about 01:30 UT, the reconstruction of the velocity is not reliable, presenting speed jumps and low cross-correlation values. Between 00:00 UT and 01:30 UT, there are no meaningful speed jumps; however, the cross-correlation is lower (ranging from about 0.5 to 0.7) compared to the period between 20:00 UT and 00:00 UT.

Figure 12 Components (blue for North–South, light blue for East–West) and magnitude (black) of the velocity measured by the LOFAR core stations. The negative sign represents a wave propagating from East to West and from North to South. The colour code represented the mean cross-correlation coming from all the baselines every 30 s.We applied the algorithm to LOFAR intensity data at 60 MHz.

Another notable feature shown by Figure 12 is that the two velocity components, when the cross-correlation is high, tend to exhibit very similar behaviour, regardless of their numerical values. To investigate this behaviour, we analyse the azimuthal angle of the velocity, reported in Figure 13. In the figure, the black line represents the azimuth of the velocity as a function of time between 16:00 UT on January 14th and 08:00 UT on January 15th. Similarly to Figure 12, the colour represents the mean cross-correlation. A regular trend appears particularly during the period with high cross-correlation. Actually, we found that the magenta lines represent the azimuthal angle of the vector $\vec{B}\times \overrightarrow{\mathrm{LoS}}$ projected along the plane tangent to the Earth’s spheroid at 350 km. For the sake of brevity, in the remaining part of the text, we will refer to the vector $\vec{B}\times \overrightarrow{\mathrm{LoS}}$ project along the plane tangent to the Earth’s spheroid at 350 km simply as $\vec{B}\times \overrightarrow{\mathrm{LoS}}$.

Figure 13 Direction of the ionospheric velocity measured using the LOFAR power data in terms of the azimuthal angle of the velocity. The colour represents the mean cross-correlation of all the baselines every 30 s while the magenta lines represent the possible azimuthal angle of $\vec{B}\times \overrightarrow{\mathrm{LoS}}$. vector. We applied the algorithm to LOFAR intensity data at 60 MHz.

This phenomenon can be explained by considering the nature of LOFAR observations and the adopted methods. First, it is important to note that the LOFAR antennas simultaneously observe Cassiopeia A. Given the proximity of the core stations relative to the distance from Cassiopeia A, the lines of sight ($\overrightarrow{\mathrm{LoS}}$) vectors connecting the antennas to the star can be assumed to be parallel, i.e., we can assume the geometrical optics approximation. Along these lines, the signal originating from Cassiopeia A and directed towards the LOFAR core stations undergoes variations caused by the ionosphere. However, because of the observational geometry, determining the exact location of these variations is not possible. This is because LOFAR is sensitive only to the integral of the effects along the $\overrightarrow{\mathrm{LoS}}$ and we can estimate velocities by assuming an irregularity layer like a thin layer, located at 350 km.

Moreover, we must consider the role of the Earth’s magnetic field, B, along whose lines the ionospheric irregularities tend to elongate. To the scope, we use the International Geomagnetic Reference Field (IGRF)-13th Generation (Alken et al., 2020). Hence, we evaluate the cross-correlation between every couple of core stations, considering no time delay (dt = 0) between signals, as a function of IPP distances.

Figure 14 shows an example of cross-correlation of all baselines by considering no delay (dt = 0) among received signals at 21:02 UT of January 14th. The cyan line represents the component tangent to the Earth’s spheroid at 350 km of the magnetic field line at the IPP projected on the plane perpendicular to $\overrightarrow{\mathrm{LoS}}$. It is noteworthy that the points with maximum cross-correlation align almost perfectly with the cyan line. The magenta lines reported in Figure 13 corresponds to the direction perpendicular to such a cyan line. The good match with the azimuth of the velocity, when the mean cross-correlation is high, means that the measured velocity refers to the projection of those structures, elongated the magnetic field lines, on the plane perpendicular to $\overrightarrow{\mathrm{LoS}}$. In other words, LOFAR core stations measure a component of the velocity projected on a direction that varies with the observational geometry and, then, with time.

Figure 14 Cross-correlation of all baselines at dt = 0. The cyan line (B × LoS) represents the direction of the component tangent to the Earth’s spheroid at 350 km of the magnetic field line at the IPP projected onto the plane perpendicular to $\overrightarrow{\mathrm{LoS}}$. The cross-correlations refer to LOFAR intensity data at 60 MHz.

This provides evidence that the ionospheric structures seen by LOFAR, the velocity of which we are measuring, are elongated along the aforementioned magnetic field line projection, whose orientation changes with the observational geometry. This comes with a price: No information can be obtained by LOFAR about the velocity component along these projections. Therefore, using the cross-correlation method to retrieve the irregularity velocity, the only perturbations that LOFAR core stations can jointly distinguish are those perpendicular to both $\vec{B}$ and $\overrightarrow{\mathrm{LoS}}$.

Another feature to investigate in Figure 13 is the presence of a stable flip in the azimuth of the velocity, observed around 02:00 UT of January 14th and 06:00 UT on January 15th. These correspond to small absolute values of the E–W and N–S components of the measured velocity. Therefore, those flips probably happen when the ionospheric velocity is fairly aligned with the magnetic field line and, consequently, the algorithm struggles to reconstruct the projection of the velocity on the $\vec{B}\times \overrightarrow{\mathrm{LoS}}$ direction.

Bearing these considerations in mind, and although the technique to retrieve the velocity suffers some limitations, it still provides useful information about the dynamics of the identified ionospheric structures. For instance, we can use this velocity estimate to verify the hypothesis that the S4/ROTI mismatch in the sub-auroral region is due to the different scales of the measurable irregularities.

During the period of main activity of the storm, around 22:00 UT of January 14th, the velocity of the perturbation can be estimated from Figure 12. It equals about $\sqrt{{\left(400\frac{\mathrm{m}}{\mathrm{s}}\right)}^2+{\left(700\frac{\mathrm{m}}{\mathrm{s}}\right)}^2}=800\mathrm{ }\mathrm{m}/\mathrm{s}$. From this and from the Nyquist frequency of the ROTI observation (2 × 30 s), we can obtain a rough estimate of the minimum dimension of the measurable irregularity from GNSS observation, under the approximation of the perpendicular velocity between the ray path and irregularity velocities. This is equal to 800 m/s × 60 s = 48 km, meaning that for GNSS data at 30 s temporal resolution it is impossible to measure irregularities caused by irregularities below about 10 km, as seen by LOFAR (Fig. 4). For LOFAR, the minimum dimension of the perturbation is given by its Fresnel’s scale. The difference in measurable scale between LOFAR and GNSS is likely the main cause of the different perturbations seen by LOFAR and GNSS in the sub-auroral region in Figure 8. However, the nature of these wave-like irregularities needs further investigation, bearing in mind what we already found: (i) they are irregularities with a typical scale size below 10 km; (ii) they have a speed ranging between from about 1,500 m/s and 100 m/s; (iii) they have a E–W component whose speed is larger than the N–S component.

To further investigate the nature of these wave-like irregularities, we compare the velocity reconstructed by LOFAR with the one retrieved by means of GNSS data. Specifically, we consider detrended TEC (dTEC) measurements from GNSS and the related keogram. The keogram in Figure 15 shows the dTEC variations in a time window between 15:30 UT and 09:00 UT of the night of between January 14th and 15th, and a latitudinal range between 35°N and 65° N. The data were processed using a 15-min bandpass filter, as described in Section 3. The alternating positive and negative variations in the dTEC values, represented by diagonal red and blue stripes, are a clear signature of Large-Scale TIDs (LSTIDs) propagating in the N–S direction. The yellow line represents the average position of the equatorward boundary of the auroral oval in the European sector, as retrieved by the DMSP-based model. The black line represents the latitudinal IPP position of the LOFAR LBA core station CS002 at an altitude of 350 km, and the gray area below that line is the corresponding N–S component of the velocity measured every 30 s. The green line highlights a LSTID propagating N–S at around 630 m/s and intersecting the IPPs of the core stations right after 20:00 UT. Numerous other LSTID signatures are evident in the keogram; however, this particular one clearly crosses the LOFAR field of view, making a direct comparison with LOFAR data and associated velocity possible. When the TID crosses the LOFAR position, the algorithm measures a velocity in the North–South direction of around 930 m/s. Then, the speeds of the irregularities reconstructed with LOFAR and with GNSS, although numerically different, present the same order of magnitude.

Figure 15 Keogram reporting the latitudinal distribution of the dTEC from 35° to 65° N for the time interval 15:30 to 09:00 between January 14th and 15th, 2022. The yellow line represents the position of the lower boundary of the auroral oval. The black line represents the latitude of the LOFAR observation projected at an IPP of 350 km. The gray segments represent the Latitudinal component of the velocity measure at the corresponding time; one degree corresponds to 100 m/s. The green segment highlights the alternations of positive and negative dTEC, showing the presence of a TID. The green dot represents the velocity measured by LOFAR at the same time and latitude as the TID seen by GNSS.

Between 22:00 UT and 00:00 UT, we obtain the most reliable reconstruction of the velocity using LOFAR. During this period, the IPPs of the core stations were close to the boundary of the oval. As a result, we reconstruct the velocity only within a narrow band of a few degrees equatorward of the auroral oval. When the IPPs are at lower latitudes or inside the auroral oval, the cross-correlation weakens, indicating a loss of coherence in the structures.

A possible explanation about the nature of the small-scale (<10 km) wave-like disturbances is that they are related to the triggering process of LSTIDs. In fact, we reconstruct the coherent wavelike structure only in “far field” relatively to the auroral precipitation, likely capturing TID-related disturbances in which LSTID are forming. In the “near-field” region, i.e., when IPPs are crossing the selector in which particle precipitation occurs, the irregularities are neither large nor persistent enough to coherently traverse all baselines. This is also true when the IPPs fall below the previously mentioned narrow latitudinal band below the oval, i.e., before 20:00 UT. Despite this, that period is marked by IPPs crossing LSTID signatures in the keogram, particularly around 18:30 UT–19:00 UT. Then, in this case, no apparent relationship seems to exist between the LOFAR-detected wave-like irregularities and the LSTIDs.

Due to the fact that the LOFAR irregularities have an E–W component whose speed is larger than the N–S component, we cannot be completely sure that disturbances observed in LOFAR may have triggered the occurrence of TIDs seen by GNSS, which are conversely featured by larger values on the N–S component.

Another option deals with the association of these wave-like structures with a series of phenomena occurring in the ionospheric sub-auroral region. This is region is highly dynamic, especially during storm and substorm activity. In our case, the larger values of cross-correlation occur in correspondence with the larger activity of the westward currents, as demonstrated by the IL behaviour (Orange line in Fig.1c), especially between 20:00 and 22:00 UT. This marks the possible intensification of enhanced westward flows (<1 km/s), such as polarization jets, sub-auroral ion drifts (SAIDs), and sub-auroral polarization streams (SAPS) (Horvath & Lovell, 2021). The very large speed featuring these classes of sub-auroral phenomena matche with the velocity of the irregularities reconstructed by LOFAR and with the fact that these structures, whenever coherent, are confined in a narrow latitudinal band of a few degrees, in agreement with the literature. In other words, this would justify why, right after 20:00 UT, we obtain a velocity of ≃1,500 m/s with such a good mean cross-correlation right below the auroral oval.

5 Conclusions

In this study, we successfully leverage the combined capabilities of LOFAR and GNSS to investigate ionospheric irregularities during the geomagnetic storm of January 14th–15th, 2022. While the LOFAR and GNSS observations reveal distinct aspects of the ionospheric response to the geomagnetic storm, the specific nature of a direct causal link between small-scale irregularities detected by LOFAR and large-scale TIDs remain an open question that warrants further investigation. By integrating data from these two systems, however, we provided a comprehensive analysis of ionospheric disturbances, focusing particularly on the sub-auroral sector. Our analysis identified three distinct phenomena associated with the geomagnetic storm: The enhancement of direct particle precipitation within the auroral oval, the steepening of the equatorward edge of the high-latitude ionospheric trough and the propagation of wave-like structures in the sub-auroral latitudes, characterised by scales of a few kilometres and velocities of several hundred metres per second.

Thanks to the higher resolution of LOFAR data, we set up two methods to reconstruct the velocity of these wave-like structures and we discussed the advantages and disadvantages of them. We provided evidence that with LOFAR core stations, it is possible to measure the component of this velocity perpendicular to both $\vec{B}$ and $\overrightarrow{\mathrm{LoS}}$ and then projected on the plane tangent to the Earth’s spheroid at an altitude of 350 km. From the same algorithm, we also obtain a measure of the goodness of the algorithm itself. These observations contribute to a deeper understanding of the dynamic processes occurring in the ionosphere during geomagnetic storms.

However, we are able to reconstruct a reliable velocity when the IPPs are positioned sufficiently far from the auroral oval to have stable irregularities, yet close enough for the perturbations in the LOFAR measurements to remain strong, allowing for efficient cross-correlation of signal variations induced by the ionosphere between each pair of stations. We speculate on the possible nature of these wave-like structures, and identify two potential explanations (which may concur):

• they are related to the triggering process of LSTIDs and are wave-like modulations jointly travelling with them;

• they are related to intensification of enhanced westward flows in the sub-auroral ionosphere.

Further research is needed to explore the nature of such irregularities. This includes the relationship between velocities derived from LOFAR and those from GNSS, during more active periods where both measurements are available and a statistical assessment, which are currently under realisation.

Acknowledgments

The editor thanks P.T. Jayachandran and an anonymous reviewer for their assistance in evaluating this paper.

Funding

This work is framed under the activities of the Plasmasphere Ionosphere Thermosphere Integrated Research Environment and Access services: A Network of Research Facilities (PITHIA-NRF) project, which has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101007599. The authors gratefully acknowledge financial support from the Monitoring Earth’s Evolution and Tectonics (MEET) project, in the frame of NextGenarationEU.

Data availability statement

This paper is based on data obtained with the International LOFAR Telescope (ILT) under project code LT16_002 available through LOFAR Long Term Archive (LTA): https://lta.lofar.eu/Lofar.

The DMSP SSUSI data were obtained from The Johns Hopkins University Applied Research Laboratory and are available at: https://doi.org/10.5281/zenodo.17277629.

GNSS data are provided by the National Land Survey Finland (NLS) (https://www.maanmittauslaitos.fi/en/finpos/rinex), SWEPOS Sweden (https://www.lantmateriet.se/en/geodata/gps-geodesy-and-swepos/lantmateriets-doi-objects/swepos-rinex-data/), Norwegian Mapping Authority (Kartverket, https://www.kartverket.no/en/on-land/posisjon/guide-to-etpos), and by the “Archivio Dati GNSS Centralizzato” of INGV (ftp://gnssgiving.int.ingv.it), which collects GNSS data from eighty networks in the Euro-Mediterranean region and Africa.

The DSCOVR data are available at National Centers for Environmental Information of NOAA: https://www.ngdc.noaa.gov/dscovr/portal/index.html#/vis/summary/1d/1642114800000.

The Kp data are available at: https://doi.org/10.5880/Kp.0001.

Dst index data are available at: https://wdc.kugi.kyoto-u.ac.jp/dstdir/.

The FMI-IMAGE electrojet indicators are available at https://space.fmi.fi/image/www/index.php?page=il_index.

The data from the Juliusruh ionosonde are available at: https://giro.uml.edu/ionoweb/#.

The data from the Sodankyla ionosonde are available at https://www.sgo.fi/pub_ion/IonoDaysView/IonoView_2022/IonoDay_2022_01_14.pdf and https://www.sgo.fi/pub_ion/IonoDaysView/IonoView_2022/IonoDay_2022_01_15.pdf.

Supplementary material

Supplementary material A: Evolution of geographical maps of ROTI and S4 in the period between 16:00 UT of January 14th, 2022 to 6:00 UT of January 15th, 2022, with a time step of 5 minutes. Blue dots show the LOFAR data from all the international stations and the core station CS002. The blue shade intensity represents the S4 magnitude for the LOFAR selected frequency. If the S4 measured by a station is below the threshold of 0.05, the data is represented with a red 'x'. The yellow/black points, represents the ROTI value from GNSS above 0.2 TECU/min. Values below this threshold are indicated with a small grey dot. The red lines represent the edges of the modelled auroral oval, obtained from the closest (in time) imaging data from the DMSP satellite. Both LOFAR and GNSS data are projected at an IPP of 350 km, if outside the auroral oval, or 120 km, if inside of it. The title of each map reports the time of the ROTI data and the time of the nearest auroral oval position estimation, labelled as “Auroral time”.

Supplementary material B: Evolution of normalised power measured by the LOFAR core stations projected at an IPP of 350 km, in the period between 21:20:00 UT to 21:24:59 UT of January 14th. The grey arrow reports the direction of the velocity of the perturbation reconstructed at the nearest available time, as reported in the title. The trajectory of the entire LOFAR observations (blue line) is shown in the upper left corner, with the position at the specified time indicated by a red dot.

References

All Figures

Figure 1 Total intensity |B| and Bz component of the IMF (Panel a), Solar Wind speed (Vsw) and pressure (Psw) (Panel b). Data in panel a and in panel b are measured at 1 AU. Other panels report the electrojet indicators IE, IL, IU in the European sector (Panel c), Kp (Panel d) and Dst (Panel e) indices during the storm. The colour code in panel d indicates the NOAA Space Weather Scale for geomagnetic storms. The blue shaded area in all panels represents the period between 16:00 UT of January, 14th and 08:30 UT of January 15th, the focus period of our analysis.

In the text

	Figure 2 GNSS receivers and LOFAR stations map. Panel a) shows the locations of GNSS receivers across Europe, LOFAR international stations, and the CS002 core station, as representative of all CSs. Panel b) provides a zoomed-in view of the LOFAR CSs for greater detail.
In the text

	Figure 3 Normalised power measured by the LBA antenna of CS002. In red, we highlighted a 2 MHz wide band around 60 MHz, the frequency we will use to investigate the variability induced by the ionosphere on the signal measured by LOFAR.
In the text

	Figure 4 The figure represents the LOFAR Fresnel’s scale at a frequency of 60 MHz, considering the altitude of a single ionospheric layer at 350 km.
In the text

Figure 5 Example of roll-off frequency estimation from the PSD based on FFT of 5 min of data between 23:22:30 UT and 23:27:30 UT of January 14th. The blue line represents the log-log spectrum, while the red line shows the fitted function following equation (5). The roll-off frequency is indicated by the dashed blue line, while its standard deviation is represented by the light blue box around it. In the gray box, we report the best fit values and their standard deviation for the constant, the slope and the intercept of equation (5).

In the text

	Figure 6 ROTI and S4, measured at 60 MHz, time variability depending on latitude. All the points are projected at an IPP altitude of 350 km. Only data between 45° N and 65° N are shown. This is the same observational range of LOFAR data during this storm.
In the text

Figure 7 Map of the ionospheric irregularities measured around 20:05 UT of January14th over the European sector at latitudes between 45° and 80° N and longitudes between −15° and 25° E. The red lines represent the boundaries of the auroral oval modelled by the DMSP. The ROTI data are represented in the yellow–black scale if over the threshold of 0.2 TECU/min or as simple gray dots if under it. The LOFAR S4 follow logic with a threshold of 0.05 and a scale green/blue and red “x”.

In the text

Figure 8 Map of the ionospheric irregularities measured around 22:00 UT of January14th over the European sector at latitudes between 45° and 80° N and longitudes between −15° and 25° E. The red lines represent the boundaries of the auroral oval modelled by the DMSP. The ROTI data are represented in the yellow–black scale if over the threshold of 0.2 TECU/min or as simple gray dots if under it. The LOFAR S4 follow the same logic with a threshold of 0.05 and a scale green/blue and red “x”.

In the text

Figure 9 Plot a) shows the map of the ionospheric irregularities measured around 23:10 UT of January14th over the European sector at latitudes between 45° and 80° N and longitudes between −15° and 25° E. The dark green line represents the trajectory of the Swarm Bravo satellite passing over Europe, from South to North, between 23:07 UT and 23:16 UT on January 14th. The light green line indicates the in-situ plasma density of Swarm Bravo, with a higher longitudinal positive deviation from the ground track, the higher the plasma density is. The other variables follow the same colour description as Figures 7 and 8. Plot b) shows the in-situ electron density measured by Swarm Bravo (red line) between 23:07 UT and 23:16 UT of January 14th over the European sector. The blue dots represent the ROTI values measured by the GNSS. The yellow band represents the enhancement of ROTI between 0.2 and 0.5 TECU/min due to the presence of the poleward boundary of the ionospheric trough.

In the text

Figure 10 Normalised power measured by the LOFAR core stations projected at an IPP of 350 km for four consecutive instants, i.e., 21:23:02 UT to 21:23:05 UT of January 14th. The grey arrow reports the direction of the velocity of the perturbation. We report the nearest velocity available, for which the time is in the plot title. The trajectory of the LOFAR observations (blue line) is shown in the upper left corner, with the position at the specified time indicated by a red dot.

In the text

Figure 11 Velocity estimates during the storm by using cross-correlation (CC) and velocity from Fresnel’s Frequency identification. The blue line represents the mean velocity of all the core stations obtained using the Fresnel’s algorithm (Eq. (4)), while the black line represents the velocity obtained by the cross-correlation algorithm (Eq. (3)). Both algorithms refer to LOFAR intensity data at 60 MHz.

In the text

	Figure 12 Components (blue for North–South, light blue for East–West) and magnitude (black) of the velocity measured by the LOFAR core stations. The negative sign represents a wave propagating from East to West and from North to South. The colour code represented the mean cross-correlation coming from all the baselines every 30 s.We applied the algorithm to LOFAR intensity data at 60 MHz.
In the text

	Figure 13 Direction of the ionospheric velocity measured using the LOFAR power data in terms of the azimuthal angle of the velocity. The colour represents the mean cross-correlation of all the baselines every 30 s while the magenta lines represent the possible azimuthal angle of $\vec{B}\times \overrightarrow{\mathrm{LoS}}$. vector. We applied the algorithm to LOFAR intensity data at 60 MHz.
In the text

	Figure 14 Cross-correlation of all baselines at dt = 0. The cyan line (B × LoS) represents the direction of the component tangent to the Earth’s spheroid at 350 km of the magnetic field line at the IPP projected onto the plane perpendicular to $\overrightarrow{\mathrm{LoS}}$. The cross-correlations refer to LOFAR intensity data at 60 MHz.
In the text

Figure 15 Keogram reporting the latitudinal distribution of the dTEC from 35° to 65° N for the time interval 15:30 to 09:00 between January 14th and 15th, 2022. The yellow line represents the position of the lower boundary of the auroral oval. The black line represents the latitude of the LOFAR observation projected at an IPP of 350 km. The gray segments represent the Latitudinal component of the velocity measure at the corresponding time; one degree corresponds to 100 m/s. The green segment highlights the alternations of positive and negative dTEC, showing the presence of a TID. The green dot represents the velocity measured by LOFAR at the same time and latitude as the TID seen by GNSS.

In the text