© D. Altadill et al., Published by EDP Sciences 2020

1 Introduction

Travelling Ionospheric Disturbances (TIDs) are the ionospheric manifestation of internal Atmospheric Gravity Waves (AGWs) in the neutral atmosphere (e.g., Hines, 1974; Hocke & Schlegel, 1996). TIDs are detected as plasma density fluctuations that propagate through the ionosphere with velocities and amplitudes that depend on the source of the AGW (e.g., Fedorenko et al., 2013). The source of excitation can be anthropogenic (e.g., nuclear explosions) or natural (geomagnetic storms, auroral activity intensification, tropospheric convection, etc.). TIDs result in changes in the Total Electron Content (TEC) and ionospheric electron density distribution affecting the operation of HF geolocation, HF communication, and the Global Navigational Satellite System (GNSS); the TID can impose Doppler frequency shifts of the order of 0.5 Hz on HF signals. TIDs are usually classified into two groups depending to their wave velocity and period (Hocke & Schlegel, 1996). Large-scale TIDs (LSTIDs) have fluctuating periods (T) of 50 min to 3 h and horizontal velocities between 400 and 1000 m s⁻¹. Medium-scale TIDs (MSTIDs) have horizontal velocities between 100 and 250 m s⁻¹ and T in the range of 15–60 min. LSTIDs are related to auroral geomagnetic activity (e.g., Tsugawa et al., 2004), and have wavelengths (Λ) greater than 1000 km. During geomagnetically disturbed periods the high latitude thermosphere is heated via the Joule effect, and consequently energy is transferred towards lower latitudes in the form of thermospheric waves that interact with the ions in the ionospheric plasma (Prölss & Očko, 2000; Paznukhov et al., 2009; Karpachev et al., 2010). The MSTIDs are generally related to meteorological phenomena, and have wavelengths of 100–300 km (e.g., Vadas & Crowley, 2010). The energy dissipation of the MSTID can excite secondary LSTID and explain poleward propagation of LSTIDs as reported by Ding et al. (2013). The TID fluctuations can be detected with various ionospheric sensors such as ionosondes (e.g., Hajkowicz, 1991; Galushko et al., 2003), coherent/incoherent scatter radars (e.g., Evans et al., 1983), and with the analysis of TEC data from GNSS measurements (e.g., Hernandez-Pajares et al., 2006).

Studying the TID phenomena expands the understanding of the dominant energy distribution and momentum transfer mechanisms in the ionosphere and thermosphere, in turn enabling the development of advanced warning and mitigation applications for systems relying on predictable ionospheric radiowave propagation. Recently, in the framework of the TechTIDE project (Warning and Mitigation Technologies for Travelling Ionospheric Disturbances Effects; http://www.tech-tide.eu/), different methodologies to detect TIDs have been developed using data provided by Digisonde ionospheric sounders, by Continuous Doppler Sounding System and by GNSS receivers (Galkin et al., 2018). This manuscript focusses on the so called HF-Interferometry method, which was initially developed in the framework of the Net-TIDE project (Pilot Network for the Identification of Travelling Ionospheric Disturbance; https://sites.google.com/site/spsionosphere/) and finally developed and improved in the framework of the TechTIDE project.

HF-Interferometry, hereafter HF-Int, is a technique for the identification of LSTIDs in near-real time based on observations made by the European and South African networks of ionosondes. The method detects quasi-periodic oscillations of ionospheric characteristics, identifies coherent oscillation activity at different measuring sites of the network, and calculates the propagation parameters, velocity and azimuth of the perturbation. Due to the characteristics of the sensor network, especially the large distances between Digisondes and the relatively coarse cadence of measurements, the HF-Int method detects mainly large scale TIDs, or LSTIDs. However, the method could detect MSTIDs if the networks of Digisondes were denser (separated by less than typical wavelengths of MSTIDs) and if cadences of measurements were shorter than 5 min.

This paper is organized as follows. Section 2 presents a description of the method, Sections 3 and 4 are respectively devoted to the verification and validation of the method, including results of the application of the method for well-established LSTIDs events. Section 5 presents the summary and discussion of the results, the limitations of the technique for detection of LSTIDs and potential applications.

2 HF interferometry method: identification and tracking of LSTIDs using Digisonde data

HF-Int is a method to identify and track LSTIDs over regions with a suitable network of HF sensors; e.g., Digisondes (Reinisch et al., 2009). HF-Int has been developed to support the warning services of the TechTIDE project. HF-Int determines the dominant period of oscillation and the amplitude of the LSTIDs using spectral analysis, and estimates the propagation parameters of the LSTIDs from the measured time delays of the disturbance detected at different sensor sites. HF-Int uses the ionospheric characteristics MUF(3000)F2 or foF2 provided by 10 Digisondes in Europe and 4 Digisondes in South Africa (Fig. 1), which make routine vertical incidence ionograms. MUF(3000)F2 is the maximum usable frequency for a 3000 km radio-link with a single hop in the F2 region of the ionosphere; foF2 is the ordinary wave critical frequency of the F2 layer, which is proportional to the square root of the maximum electron density. Both characteristics, foF2 and MUF(3000)F2, are obtained in the ionograms. Table 1 lists details of the ionospheric stations providing data to the HF-Int. Data is obtained in near real time from the Global Ionospheric Radio Observatory (GIRO) DIDBase Fast Chars database (http://giro.uml.edu/didbase/scaled.php).

Fig. 1 Geographic locations of the Digisonde stations used for the HF-Int method in (a) Europe and (b) South Africa.

To identify and distinguish TIDs from other irregular disturbances in the time series we apply criteria based on the fact that TIDs are quasi-periodic fluctuations according to the classification of TIDs (Hocke & Schlegel, 1996). We identify the significant periodicities in the time series within the periodic range of the TIDs using spectral analysis of the data collected at different sites of the network. For this purpose, we compute the periodograms of the above time series using the Lomb-Scargle method (Scargle, 1982; Hocke, 1998). Then we make sure that the detected periodic variations are coherent for different sites of the network. Finally, we obtain the time delays of the disturbance at different sites of observation by cross-correlation analysis, and we estimate the propagation velocity of the TID assuming that the disturbance propagates as a plane wave. This is a first order approximation that can be done with the available data, and a more precise estimation of the velocity would require a higher cadence of measurements. We note that both Digisonde networks, in Europe and South Africa, have measuring sites separated by about 500 km on average, and typically make observations with 15 min cadence, which allows detecting only LSTIDs (Λ ≥ 1000 km and T ≥ 45 min). For the HF-Int, we take into account that the ionospheric characteristics contain strong daily trends that may differ significantly at stations separated by large distances. It is expected that any perturbation that propagates over the network of Digisondes can be seen in the residuals of the time series after removing the aforementioned trends. Moreover, as mentioned, most of the Digisondes collect data with 15 min (900 s) cadence (see Table 1), while a higher sampling rate would of course be desirable for a more precise estimation of the disturbance velocity. For example, a perturbation propagating with a velocity of 500 m s⁻¹ should be seen at the Digisondes with a baseline of 500 km (typical separation in our network) with a delay of 1000 s. Therefore, the current application of HF-Int for the European and South African Digisonde networks is limited to the detection of LSTIDs.

The data processing of the HF-Int to detect LSTIDs and to obtain the velocity of the disturbance consists of three steps: removal of the daily trend from the measured Digisonde data, spectral analysis to detect possible TID-like variations, and cross-correlation to obtain time delays of the disturbance between different sites in order to estimate the vector velocity.

2.1 Data preprocessing

Before the application of the HF-Int method the raw data is preprocessed in order to remove the daily trend and to increase the time resolution. The original data sampling rate (5–15 min) is rather coarse and results into a rough estimate of the time delays which is used for the calculation of the disturbance parameters. Preprocessing consists of Discrete Fourier Transform (DFT) interpolation and high-pass filtering. First, the raw data is DFT interpolated to increase the sampling rate (upsampling). This procedure consists of calculating the DFT spectrum of the raw data, and then doing the inverse DFT calculation using the obtained spectral amplitudes and specifying the sine and cosine spectral components with finer time step (higher sampling rate) than the raw data. Note that in the preprocessing, the real DFT (Smith, 2002) calculation is used, although the algorithm is equivalent to spectral zero-padding interpolation made with complex DFT calculation (Oppenheim et al., 2001). In the preprocessing, the data is typically upsampled to obtain a 1 min sampling rate which is kept the same for all stations (the raw data may be collected with different sampling rates at different stations). Then, the DFT spectrum is high-pass filtered in order to eliminate the slowly varying daily trend from the record. The remaining spectrum produces high frequency residuals which are associated with the wavelike ionospheric disturbances. Typically, the data cutoff frequency is set corresponding to the period of 3 h, which allows preserving different types of ionospheric waves for further analysis.

HF-Int calculations are made in near real time, using the most recent 24 h of the collected data in the preprocessing to optimally remove the daily trends. The actual HF-Int calculation of the disturbance characteristics uses the most recent 6-hour long data record, with each new calculation made every 5 min.

2.2 Detection of TID-like variation

The second step of the HF-Int method consists of applying spectral analysis in order to detect possible TID-like variations (periodicities between 38 and 160 min) in the 6-h long time interval of the high-pass filtered residual data. To do so, HF-Int computes periodograms (Scargle, 1982; Hocke, 1998) for the residuals from the individual measurement sites of the network. Then, the two periods with the largest amplitudes and with confidence level above 0.975 (if any) for each of the measuring sites are selected as TID-like variation candidates for which the periods observed at different sites are compared. If they are coherent for at least three different sites of the network (the candidate periods differ by less than 30%) we estimate the amplitude, the dominant period (T _AM) and the Spectral Energy Contribution (SEC) of the possible TID that satisfy the above criteria for each of the measuring sites. Of course, if the candidate periods have a confidence level lower than 0.975, or are not coherent for at least three different sites, no TID activity is reported. HF-Int also provides information about the contribution of the TID-like variation to the total variability of the analyzed time series. Applying Parseval’s relation, similarly to Altadill et al. (1998), HF-Int estimates the SEC of the periodic range of the TIDs to the total energy which is equivalent to the contribution of the TIDs to the total variability of the given time series:

$$\sum _{n=-\infty }^{\infty }{ \left|x\left[n\right]\right|}^2=\frac{1}{2\pi }\underset{-\pi }{\overset{\pi }{\int }}{\left|X\left(\omega \right)\right|}^2\sim \sum _{T=T_{\mathrm{S}}}^{T=T_{\mathrm{E}}}{A\left(\omega \right)}^2$$(1)

$$\mathrm{SEC}\left(\%\right)=\frac{{\sum }_{T=T_{{\mathrm{TID}}_{\mathrm{S}}}}^{T=T_{{\mathrm{TID}}_{\mathrm{E}}}}{A\left(T\right)}^2}{{\sum }_{T=T_{\mathrm{S}}}^{T=T_{\mathrm{E}}}{A\left(T\right)}^2}$$(2)where ω is the angular frequency of the period T, T _TIDS and T _TIDE are the starting and ending bounds of the periodic range of the TID-like variation, respectively, and T _S and T _E are the starting and ending periods of the total periodic range under analysis. In the HF-Int method, T _S = 30 min and T _E = 360 min, and T _TIDS = 0.8T _AM and T _TIDE = 1.4T _AM.

2.3 Estimation of the velocity of propagation of the disturbance

Finally, to verify whether a detected TID-like variation can be considered a TID, we compute the propagation velocity (speed and azimuth) of the disturbance. To do so, we have adapted a technique used for calculation of the cross-correlation between different time series of GNSS data (Hernandez-Pajares et al., 2006) to the processing of the multi-station Digisonde data. The method, which is very similar to Wan et al. (1997), consists of reconstructing disturbance propagation parameters from the measured time delays of the signal at different sensor sites, i.e., HF-Int assumes a pure plane wave propagation of the disturbance as a test hypothesis.

Therefore, for each ith Digisonde station, the time delay ∆t _i of the disturbance with respect to the reference station can be written as:

$$∆t_{\mathrm{i}}=\vec{s}\bullet ∆{\vec{r}}_{\mathrm{i}}$$(3)where $\vec{s}=\vec{v}/v^2$ is the slowness vector of the disturbance propagating at velocity $\vec{v}$, which is an unknown variable; $∆{\vec{r}}_{\mathrm{i}}$ is the relative position vector of the ith Digisonde with respect to the reference one. For a network of stations with i > 3, this system of equations is overdetermined and can be solved by a least squares fitting method.

At this point, it is possible to compare the perturbations from different sites with respect to a reference one. The time delays of disturbance Δt _i are obtained by calculating the time-lag for maximum cross-correlation (CCM_i). In this process it is important to select the central part of the filtered signals in order to mitigate the edge effects. For this reason, only a 360 min long record is selected for the reference signal, and the correlation maximum with the signals from other receivers is searched within ±60 min. Then, if Δt _i = ±60 this station is discarded.

Once the time delays Δt _i are estimated for each site “i”, and given the relative vector position (Δx _i, Δy _i), the following relationship must be fulfilled according to equation (3),

$$∆t_{\mathrm{i}}-\vec{s}\bullet ∆{\vec{r}}_{\mathrm{i}}=0.$$(4)

Equation (4) is solved by least square error fitting to estimate the slowness vector, from which one can obtain the propagation velocity according to the relationship,

$$\vec{v}=\frac{\vec{s}}{s^2}.$$(5)

We note that in the estimation process, data from the reference site are assigned higher weight than data from the other sites, i.e., data are weighted according to the maximum correlation coefficient (CCM_i). Therefore, the propagation parameters should be considered as “linked” to the reference site. Changing the reference station in the fitting procedure could result in different results, if the hypothesis of a pure plane wave propagation is not satisfied over such large distances.

As stated above, for estimating each Δt _i, the cross-correlation maximum is searched within ±60 min into a 360 min time interval, which is sufficient taking into account the expected TIDs periods. However, the large baselines between ionosondes and the oscillating nature of the TIDs can result into observing more than a single CCM_i providing, for instance, at a positive and a negative Δt _i value. This would represent an ambiguity in the solution of equation (4). However, this is not the case: the later apply for pure periodic variations and TIDs quasi-periodic (Hunsucker, 1982; Hocke & Schlegel, 1996); and one has to take into account that equation (4) is a linear relationship that, in fact, can be solved separately for s _x and s _y components of the slowness vector. Then, the relevant issue is not the distance with respect to the reference ionosonde but the coordinate differences (Δx _i and Δy _i) between ionosondes. In addition, Δx _i and Δy _i range from tens up to thousands of km and should have consistent values of Δt _i for all baselines. In this sense, selecting a wrong Δt _i for a specific baseline (for instance, a negative value instead of a correct positive one) would introduce an inconsistent value for the other baselines. In summary, the redundancy of ionosonde baselines helps to break any ambiguity in the estimation of Δt _i.

Therefore, to detect and estimate the propagation velocity of a TID, the HF-Int method applies the above cross-correlation technique to the reference sites that comply with the restrictions described in Section 2.2, and makes sure that the CCM_i > 0.5 for at least three sites of the network. Thus, if the candidate disturbance detected in the reference site does not sufficiently correlate for at least three sites of the network, then no TID activity is reported. This way, HF-Int ensures that the disturbance is coherent for a given region of the network. However, the latter might be difficult to accomplish for the Digisonde network in South Africa which has only four measuring sites in total. That is why HF-Int makes sure that the CCM_i > 0.5 for at least two sites in the South African network.

3 Verification of the HF Interferometry method

To verify the capability of the HF-Int method to detect TID events, we have created a set of synthetic data using the following analytic function that represents the propagation of a plane wave:

$$F=F_0\left(1+A{e}^{{-\left(t-t_0\right)}^2/b}\cos \left(\vec{k}\vec{x}-\omega t\right)\right)$$(6)where F ₀ is the background MUF(3000)F2 value obtained from IRI model (Bilitza et al., 2017) for different epochs of the year, the$\cos \left(\vec{k}\vec{x}-\omega t\right)$gives the oscillation of the wave and its spatial propagation, and $\vec{k}$ is the wave vector ($\vec{k}\vec{v}=\omega$). Exponential coefficients (A and b) allow us to control the amplitude and the duration of the perturbation respectively. The value t ₀ is the time when the amplitude of the perturbation reaches the maximum. Figure 2 shows the two simulated TIDs and the results of the application of the method to these synthetic data. Figure 2a shows the background calculated by IRI for EB040 in winter and two perturbations added using equation (6). The first perturbation is centered at 6 UT, has a period of 120 min, and propagates with a speed of 500 m s⁻¹ and azimuth of 270° (simulating a TID linked to sunrise effect). The second TID is centered at 21 UT, has a period of 80 min, and propagates with speed of 1500 m s⁻¹ and azimuth of 180° (simulating a TID of auroral origin). In addition to simulating data for the EB040 Digisonde station, we have computed a synthetic 24-hour long data set with a sampling interval of 15 min for each station of the European network (Fig. 1) to be used as input for the HF-Int method verification. Figures 2b and 2c shows the results of applying HF-Int method for DB049 and EB040 stations respectively. Red dots show the calculated TID velocity, blue dots show TID azimuth, TID period is shown with black dots, and spectral contribution in shown in green. TID period, velocity and azimuth reconstructed with the HF-Int method are in a good agreement with the simulated TID parameters. The larger discrepancies are observed at the beginning and the end of the detection. These edge discrepancies are due to the fact that the analysis window does not exactly match the duration of the simulated disturbance to precisely detect a “periodic” variation which is optimally one full cycle. Note also that the SEC maximizes when the 6-h analysis window overlaps the full wave disturbance exactly.

Fig. 2 Verification of the HF-Int method with synthetic data. (a) Simulated time variations of the MUF(3000)F2 for EB040 on 15 January 2017, with two TIDs events generated analytically according to equation (6) (red line) with the indicated characteristics. Black dots correspond to the background calculated by IRI. (b) and (c) Results of the HF-Int method for DB049 and EB040 stations respectively.

4 HF interferometry results and comparison with other techniques

This section analyzes the 12–13 October 2015 event as a case study, presenting results obtained with the HF-Int method and comparing them with the results from other methods to validate the HF-Int method performance in the detection of LSTIDs. Quiet to minor storm levels were observed on 12–14 October 2015 as Earth’s geomagnetic field was under the influence of a positive polarity Coronal Hole High Speed Streams (CH HSS). During this time frame, solar wind speeds reached a peak of about 615 km s⁻¹ early on 13 October before declining to about 450 km s⁻¹ by the end of the 14th. The geomagnetic activity on 12–13 October 2015 reached a value 5 of the Kp index during these 2 days, and the AATR ionospheric activity indicator (Juan et al., 2018) estimated moderate activity at auroral region for the night of 12–13 October 2015. For this night, some TID signatures were observed above Europe (Verhulst et al., 2017).

4.1 TID identification with FAS technique

The Frequency-and-Angular-Sounding (FAS) technique is based on bistatic Digisonde-to-Digisonde (D2D) “skymap” measurements that probe the ionospheric midpoint between transmitter and receiver using the 1-hop F2 (1F2) HF signal path (Reinisch et al., 2018). The DPS4D at the receiving end of the fixed-frequency HF link measures the signal observables: time-of-flight, Doppler frequency, and the elevation and azimuth angles of the arriving ray, selecting either ordinary or extraordinary wave polarization. The FAS analysis derives the amplitude, A _N, period, T, velocity, V, and propagation azimuth, Θ, of the fundamental TID wave from the measured observables (e.g., Beley et al., 1995; Galushko et al., 2003; Paznukhov et al., 2012).

On 12–13 October 2015, the D2D FAS technique, operating under the NetTIDE project website, detected a LSTID event with the 1082-km D2D link from Dourbes (Belgium) to Roquetes (Spain) (Verhulst et al., 2017). For this event, a TID wave with a period of 80 min propagating with a velocity of 170–330 m s⁻¹ in the direction 165–200° from North was determined using the FAS method discussed in Reinisch et al. (2018).

4.2 Time variations of MUF(3000)F2

Figure 3 shows the time variations of the MUF(3000)F2 high frequency residuals (see Sect. 2.1) during the night of 12–13 October 2015 event for the European Digisonde stations located near the 15 °E meridian (left panel), and 40 °N latitude (right panel) making it possible to estimate the longitude and latitude time delays. The MUF(3000)F2 characteristic was selected because of its sensitivity to variations in both the F2 layer height and peak density (e.g., Davies, 1990). Note that the stations operate at different measurement intervals; EB040, DB049, and AT138 operate at 5 min sampling, RL052 at 10 min and the others at 15 min sampling (see Table 1). A vertical offset separates the data from the different stations to better show the time delay of irregularity signatures arriving at the different sites.

Fig. 3 (a) Time variation of the MUF(3000)F2 recorded at several Digisonde stations arranged by latitude; and (b) arranged by longitude.

The MUF(3000)F2 curves depict quasi-periodic oscillations with a period of about 75–80 min (Fig. 3). The tilted dashed lines are drawn to help following the main valleys of the oscillations, and these lines indicate delayed arrival times at the different sites. Thus the left panel shows the latitudinal (North–South) delays, and the right panel the longitudinal (West–East) delays. Note that JR055 does not observe clear oscillation activity. The irregularity signatures propagate southward from PQ052 to VT139 (10° latitude difference) in ~50 min, and eastward from EB040 to AT138 (25° longitude difference) in ~40 min. Using equation (4) it is possible to estimate the velocity vector of the propagating disturbance. The s _x and s _y components of the slowness vector are 1/370 s m⁻¹ and 1/875 s m⁻¹, respectively, indicating NW to SE propagation with a velocity of about 340 m s⁻¹. Using the numbers for the observed period (T = 80 min) and velocity (V = 340 m s⁻¹), we can estimate the disturbance wavelength of Λ = 1632 km, which classifies as a LSTID.

4.3 TID identification with GNSS data

To deduce the TID signatures for the 12–13 October event from the GNSS data, we ran the Fast Precise Point Positioning (Fast-PPP) ionospheric model (e.g., Rovira-Garcia et al., 2016), over a global set of GNSS data. The Fast-PPP model allows (in well sampled areas like Europe) separating ionospheric delays occurring in the bottom layers (VTEC1) from those occurring in the top layers (VTEC2), thus making it especially suitable for observing wave-like behavior in the lower ionosphere, and comparing it to the Digisonde bottomside measurements.

Figure 4 shows the vertical electron content estimation of the bottom layer (VTEC1) for the European region for particular ionospheric grid points, enabling to better observe latitudinal and longitudinal differences of the time variations. Similar to the arrangements of the ionospheric MUF(3000)F2 curves in Figure 3, the VTEC1 data plots in Figure 4a are arranged north–south at longitude meridian 0 °E, and those in Figure 4b are arranged east–west in the latitude parallel 43 °N. Figure 4 clearly depicts three cycles of VTEC1 oscillations during a time interval of 4 h, with a period of about 80 min, which is in a good agreement with the Digisonde observations. Figure 4a shows a phase shift (time delay) in the time variations of the VTEC1 of the southernmost grid point compared to those of the northernmost grid point, indicating a North–South propagation of approximately 6° (from latitude 46–40) in 25 min (tilted grey lines). Moreover, Figure 4b shows a West–East propagation of 10° (from longitude 355–0) in approximately 20 min (tilted grey lines). Using the equations given above, we can roughly estimate the velocity vector of the disturbance in the VTEC1, and the s _x and s _y components of the slowness vector as 1/440 s m⁻¹ and 1/675 s m⁻¹ respectively, resulting in NW to SE propagation with a velocity of about 370 m s⁻¹. These results obtained from VTEC1, agree fairly well with the results deduced from the Digisonde data, and the period and velocity of this wave-like disturbance support the suggestion that a large scale TID is observed.

Fig. 4 (a) Time variation of the VTEC1 estimated for given latitudes but for fixed longitude of 0°; and (b) time variation of the VTEC1 estimated for given longitudes but for fixed latitude of 43°.

The higher density of GNSS receivers enables us to observe the fine structure of this wavelike disturbance. However, variations in TEC values with latitude and local time might hide the relatively small TEC fluctuations produced by the TID. That is why the perturbation components of TEC are obtained by removing the average (Tsugawa et al., 2003). To highlight the TID effects on the TEC for the same region analyzed with the Digisondes we have applied the method described in Borries et al. (2009). Figure 5 depicts snapshots of the perturbation TEC maps which clearly show a band-like structure and how the main disturbance progresses from North–West to South–East, indicating an LSTID. It is also possible to visually estimate the wavelength (Λ) of this structure of about 12–15°, which is very close to the 1632 km wavelength estimated from the results presented in Figure 3.

Fig. 5 Snapshots of the perturbation TEC maps estimated over Europe (a) 22:30 UT; (b) 22:45 UT; (c) 23:00 UT; and (d) 23:15 UT of 12 October 2015. Tilted solid and dashed lines in the plots indicate the valleys and crests of the TEC disturbance respectively.

The propagation of the LSTID can also be observed in the time-latitude-plot of TEC perturbation shown in Figure 6 which is centered at 0 °E longitude.

Fig. 6 Time–latitude-plot of TEC perturbations for the 12–13 October 2015 18:00–03:00 UT at 0 °E. The arrows indicate the valleys of the disturbance TEC.

The wave fronts of the LSTID in Figure 6 show several tilted signatures of positive and negative disturbances which are indicated with arrows in Figure 6. These signatures are separated by 70–85 min from each other, indicating the period of the disturbance, in agreement with the previous period estimates from Digisonde measurements. The slope of the disturbance pattern indicates the southward propagation of TEC perturbation which agrees with the results shown in Figure 4, where crests and valleys of the disturbances are clearly seen to be shifted in time for different locations. Using the numbers estimated from the GNSS data (Fig. 4) for the period (T = 80 min) and velocity of propagation (V = 370 m s⁻¹) for this event and assuming a plane wave propagation (Λ = V · T), we can infer a wavelength of Λ = 1776 km. These results compare quite favorably with those derived from the Digisonde measurements discussed above.

4.4 TID identification with HF Interferometry method

We have run the HF-Int method for 12 and 13 October 2015 for both Europe and South Africa networks and using MUF(3000)F2 time series. As mentioned earlier, the MUF(3000)F2 characteristic was selected because of its sensitivity to variations in both the F2 layer height and peak density (Davies, 1990). Figure 7 shows the results of detection of the LSTID activity for AT139 in Europe (left plots) and GR13L in South Africa (right plots). The top plots of Figure 7 depict the time variations of the period and SEC, and the bottom plots depict the time variations of the velocity magnitude and the propagation azimuth of the disturbance. The results indicate that the TID activity starts at around 1400UT over Europe, then breaks at 1800UT, and an additional structure is observed for the night-time, with a period of about 85 min, a weak SEC (<50%), and velocity of about 400 m s⁻¹ and azimuth of propagation of about 205°. Results for South Africa observe similar features, but the LSTID activity event, which is more stable and lasts for longer time than that in Europe, shows a period of about 120 min, a strong SEC (>80%), and velocity of about 500 m s⁻¹ and azimuth of propagation of about 345°.

Fig. 7 Results of detection of the LSTID activity by the HF-Int over two stations located in (a) Europe and (b) South Africa for the time interval from 12UT on October 12 to 12UT on October 13 of 2015. Red dots show velocity, blue dots show azimuth, black dots show period and green dots show spectral energy contribution (SEC). See text for details.

Table 2 provides more details of the values of the correlation coefficient maximum (CCM_i) of the data series, time delays (Δt _i) for the CCMs, for the European network with Roquetes (EB040) serving as a reference station. Note that Table 2 refers to the 6-h time interval from 17 to 23 UT. For each station (first column), the components of the relative position vector, Δx _i (East) and Δy _i (North), are shown in the second and third column, respectively. Note that the components Δx _i and Δy _i have been obtained at a reference ionospheric altitude of 350 km above the Earth’s surface. The CCM_i and Δt _i between the reference site EB040 and the other sites in the network are shown in the fourth and fifth column.

Note that the method first checks whether the individual period of the perturbation detected for each station is coherent with the dominant period of the reference site. Those stations with non-coherent period are excluded from the calculation, this is why in Table 2 for some stations CCM_i and Δt _i values are not given. Also note that the cases for which Δt _i = ±60 or CCM_i < 0.5 are discarded too. According to these considerations, we have four stations with coherent period, which is sufficient for the velocity and azimuth calculations.

According to equations (4) and (5) and values in Table 2, the slowness vector s _x and s _y for the MUF(3000)F2 data are 4.4 s km⁻¹ and −2.1 s km⁻¹ respectively with errors of 0.9 s km⁻¹ and 0.7 s km⁻¹ respectively, meaning a relative error of 20% and 34% respectively. By applying equations (4) and (5) to the results shown in Table 2, HF-Int estimated the propagation velocity and azimuth of the TID event at 23UT above EB040 station to be 553 m s⁻¹ and 205°, respectively.

For a detailed view of the disturbance detection in both Europe and South Africa, Figure 8 shows maps of the estimated velocity vectors of the LSTID event from 2130UT to 2300UT on 12 October 2015. Results shown in Figure 8 indicate that the two networks in Europe and South Africa observed a clear equatorward propagation of the LSTID, except at 2130UT when no significant TID activity is detected in Europe. It is noticeable that both networks observe similar velocity of propagation of about 500 m s⁻¹ on average. However, the European network observes a dominant Southward propagation with a West component, and the South African network observes a dominant Northward propagation with a West component, suggesting an auroral origin of the perturbations in both hemispheres.

Fig. 8 Maps of the estimated velocity vectors of the LSTID event detected over (a) Europe and (b) South Africa on 12 October 2015 for the indicated times. Length of arrows is proportional to the velocity values which is show in the labels.

The results obtained by the HF-Int method for this case study compare very well with the results obtained by other methods, as shown in previous Sections 4.1–4.3. Figure 7 shows that the European network observes a disturbance with T = 80 min and V = 400 m s⁻¹, and the South African network observes a disturbance with T = 120 min and V = 500 m s⁻¹. Assuming a planar wave propagation, one can infer that HF-Int detects LSTIDs with wavelengths of about 1920 km and 3000 km in Europe and South Africa respectively, indicating LSTIDs of a likely auroral origin. Therefore, it is demonstrated that the HF-Int method is capable of detecting and tracking TIDs, primarily LSTIDs (due to the configuration of the Digisonde network and to the time sampling of the measurements).

5 Summary and conclusion

The HF-Int method has been developed in the framework of the TechTIDE project with the aim to detect TIDs over Europe and South Africa in near real time. The method works with the data from a network of Digisondes (10 systems in Europe and 4 in South Africa). HF-Int processing runs every 5 min and searches for quasi-periodic oscillations of ionospheric characteristics, determines the dominant period of the oscillation through spectral analysis, checks if this oscillation is coherent at different observing stations, and in a positive case, calculates velocity and azimuth of propagation of the perturbation by cross-correlation analysis. This analysis is performed automatically in near real time, as soon as new data from the Digisonde network becomes available. A non-real time version of the method has also been developed to be used for a user-specified interval of time, allowing to make retrospective analyses. We have used this version for verification and validation purposes.

In this paper, the HF-Int method has been verified by comparing simulated TID parameters (period, velocity, and propagation azimuth) with those reconstructed using the HF-Int technique. Overall, the agreement is very good, however, the HF-Int results contain some discrepancies at the beginning and end of the detection intervals. Although, in case the analysis window sufficiently overlaps with the “disturbance” to clearly detect a “periodic” variation (i.e., the window contains one full cycle of the disturbance period) the reconstructed parameters are sufficiently accurate.

In addition, the HF-Int method has been validated by comparing its performance in detecting TIDs with the results obtained by others methods for well-established TID events. For this purpose, we have focused on the LSTID event observed on the night of 12–13 October 2015 when a minor storm occurred (Kp = 5) and moderate ionospheric activity was reported at high latitudes by the AATR indicator (Juan et al., 2018). The results obtained by the HF-Int method for this case study compare very well with the results obtained by other methods, as shown in Section 4. A summary of the comparison is presented in Table 3 where the TID identification with FAS technique (Reinisch et al., 2018) is listed as FAS, the TID identification with direct observation of the time variations of MUF(3000)F2 is given as DO-MUF, the TID identification with GNSS data is shown as GNSS, and TID identification with HF-Int is referred to as HF-Int. The different results obtained by the different methods are from the different accuracy of the methods and different type of measurements. Thus, it has been demonstrated that the HF-Int methods is capable of detecting and tracking TIDs, primarily LSTIDs (because of the configuration of the networks of Digisonde and time sampling intervals of the measurements).

Systematic observations with the HF-Int method started in May 2018, and since April 2019 the results are available in near real time at the TechTIDE web site (http://www.tech-tide.eu/), where codes are also available in open access. The comparison of HF-Int results with other methodologies of TID identification is conducted routinely in the framework of TechTIDE project and the results for different methodologies make the detection and identification of TIDs available in near real time. Although not shown here, starting from April 2019, HF-Int method detected (through real time operation) several events of LSTIDs of auroral origin primarily over the European region. These events occurred when certain geomagnetic and auroral activity was also observed. In particular, TID events occurred for the nights of 23–24 April, 4–5 June, 8–9 June, the afternoon and evening of 5 August 2019, and several events occurred in 1, 3 and 4 of September 2019 as the most recent. These events were also observed with other independent TID detection techniques which are running in the warning services of the TechTIDE project such as the AATR indicator, Continuous Doppler Sounding System and by GNSS receivers techniques (Galkin et al., 2018), proving again the capability of HF-Int methods to detect and track LSTIDs.

The HF-Int method assumes a uniform plane TID wave over the entire region of analysis, and the results are therefore an estimate. Nevertheless, the general good agreement between the results provided by the HF-Int and those provided by independent techniques suggests the reasonableness of the presented methodology. These results also prove the utility of the Digisonde network as a sensor for the TID activity in well covered regions. Therefore, HF-Int can be used also for other ionosonde networks provided they are sufficiently dense.

However, current Digisonde networks in Europe and South Africa can be used to detect and track LSTID only, i.e., disturbances with Λ ≥ 500–1000 km and T ≥ 45 min, being insufficient for detecting MSTIDs. Thus, HF-Int method detects mainly large scale TIDs in Europe and South Africa. However, HF-Int could detect MSTIDs if the networks of Digisondes were denser (separated by less than typical wavelengths of MSTIDs) and if cadences of measurements were shorter than 5 min.

One possibility to overcome the above limitation would be to synchronize all the compatible sounders (like DPS4D) to obtain oblique ionospheric soundings and provide information in the middle point of two stations and increasing the time sampling of measurements.

LSTIDs affect the phase and the amplitude of trans-ionospheric or ionospherically reflected radio signals creating problems in the operation of systems operating in HF and very high frequency to ultra-high frequency (VHF-UHF) frequency bands (Davies, 1990). The presence of LSTIDs also introduces difficulties to correct interpretation of radio astronomy observations (10–250 MHz, e.g., LOFAR; Van Haarlem et al., 2013) and HF geolocation radars (Nickisch et al., 2006). The output of HF-Int method, which provides LSTID amplitudes, periods and velocities can help calibrating/postprocessing to improve radio astronomy observations or to discard observations when a large TID is present. This will make it possible to produce validated LOFAR observation and HF geolocation data.

Acknowledgments

This research has been mainly funded by H2020 TechTIDE project (GA 776011). However, it has received support also by NATO Net-TIDE project (SPS #984894), by grant PGC2018-096774-B-I00 (MCIU, AEI, FEDER), and by Universitat Ramon Llull projects 2019-URL-Proj-002 and 2019-URL-IR1rQ-002, with funds provided by “Secretaria d’Universitats i Recerca de la Generalitat de Catalunya” and “Obra Social la Caixa”. The paper uses ionospheric data from the Global Ionosphere Radio Observatory, http://giro.uml.edu/, which includes the USAF NEXION Digisonde network, and we thank the national institutions and data providers that support them. The editor thanks Eduardo Araujo-Pradere and an anonymous reviewer for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 Geographic locations of the Digisonde stations used for the HF-Int method in (a) Europe and (b) South Africa.
In the text

Fig. 2 Verification of the HF-Int method with synthetic data. (a) Simulated time variations of the MUF(3000)F2 for EB040 on 15 January 2017, with two TIDs events generated analytically according to equation (6) (red line) with the indicated characteristics. Black dots correspond to the background calculated by IRI. (b) and (c) Results of the HF-Int method for DB049 and EB040 stations respectively.

In the text

	Fig. 3 (a) Time variation of the MUF(3000)F2 recorded at several Digisonde stations arranged by latitude; and (b) arranged by longitude.
In the text

	Fig. 4 (a) Time variation of the VTEC1 estimated for given latitudes but for fixed longitude of 0°; and (b) time variation of the VTEC1 estimated for given longitudes but for fixed latitude of 43°.
In the text

	Fig. 5 Snapshots of the perturbation TEC maps estimated over Europe (a) 22:30 UT; (b) 22:45 UT; (c) 23:00 UT; and (d) 23:15 UT of 12 October 2015. Tilted solid and dashed lines in the plots indicate the valleys and crests of the TEC disturbance respectively.
In the text

	Fig. 6 Time–latitude-plot of TEC perturbations for the 12–13 October 2015 18:00–03:00 UT at 0 °E. The arrows indicate the valleys of the disturbance TEC.
In the text

	Fig. 7 Results of detection of the LSTID activity by the HF-Int over two stations located in (a) Europe and (b) South Africa for the time interval from 12UT on October 12 to 12UT on October 13 of 2015. Red dots show velocity, blue dots show azimuth, black dots show period and green dots show spectral energy contribution (SEC). See text for details.
In the text

	Fig. 8 Maps of the estimated velocity vectors of the LSTID event detected over (a) Europe and (b) South Africa on 12 October 2015 for the indicated times. Length of arrows is proportional to the velocity values which is show in the labels.
In the text