© D.G. Markowski et al., Published by EDP Sciences 2025

1 Introduction

Quantifying the day-to-day bottom-side ionospheric variability at both spatial and temporal scales has been a long-standing problem in the scientific community (Munro & Appleton, 1950; Hines, 1960; Forbes et al., 2000). From early studies by Hines (1960) and Hunsucker (1982) it has been well established that traveling ionospheric disturbances (TID) are a major contributor to that variability. This is due in part to their ability to create wave-like disturbances in the ionospheric electron density. As a consequence, there are fluctuations in the reflection height of high-frequency (HF) radio waves (Davies, 1990).

TIDs have been broadly categorized into two groups (Hunsucker, 1982), large-scale TIDs (LSTIDs) and medium-scale TIDs (MSTIDs). LSTIDs usually have periods between 30 and 120 min, have horizontal velocities of 300–1000 m/s, and have horizontal wavelengths of approximately 300–3000 km (Hocke & Schlegel, 1996). One of the main drivers of LSTIDs is believed to be geomagnetic activity, which can create significant atmospheric gravity waves, eventually causing LSTIDs (Richmond, 1978; Crowley et al., 1987; Hunsucker, 1982). While the basic physics of atmospheric gravity wave propagation seems to be understood, possible sources for LSTIDs remains unresolved.

It is well-documented that MSTIDs usually have periods between 10 and 60 min, horizontal velocities of 50–300 m/s (Jacobson et al., 1995) and horizontal wavelengths of 100–300 km (Hocke & Schlegel, 1996). MSTIDs are thought to be to be caused by either atmospheric gravity waves propagating up to ionospheric heights (Hines, 1960; Hooke, 1968) or electrodynamics processes typically associated with the Perkins instability (Perkins, 1973; Kelley, 2011; Chou et al., 2017). However, similar to LSTIDs, establishing the driving mechanism is still a priority (Hocke & Schlegel, 1996; Zhang et al., 2019).

As stated above, understanding the spatial and temporal propagation characteristics of TIDs remains an open question in the community, and as a result, several studies have explored these characteristics using a range of different instruments. Two radar systems used for these investigations are incoherent scatter radars (Kirchengast et al., 1995; Galushko et al., 2003; Nicolls & Heinselman, 2007; Vlasov et al., 2011; Negale et al., 2018), and coherent scatter radars (Fukao et al., 1991; Frissell et al., 2014). Both techniques are powerful tools for making observations of electron density as a function of altitude and time. However, ISRs have very limited 2-D spatial coverage, although these instruments have excellent altitude resolution. Coherent scatter radars, i.e., SuperDARN, have much better spatial coverage but require more assumptions to determine TIDs. GPS TEC measurements (Valladares & Hei, 2012; Crowley et al., 2016; Kelley et al., 2023) and airglow imagers (Mendillo et al., 1997; Shiokawa et al., 2009) are also commonly used tools since they have excellent spatial coverage. However, TEC and airglow are both height integrated quantities. In addition, airglow imagers are great for studying the bottom-side variability; however, they rely on clear skies to obtain measurements, which can greatly limit data availability in certain locations. GPS TEC measurements are less sensitive to the bottom-side ionosphere but have much more data availability.

HF radio wave sounding techniques have also been exploited extensively to understand the propagation characteristics of TIDs. Two examples are HF Doppler Systems and ionosondes. To get TID characteristics with HF Doppler systems (Georges, 1968; Waldock & Jones, 1986; Chum et al., 2014; Chilcote et al., 2015; Oinats et al., 2016; Trop et al., 2025), the frequency-and-angular sounding (FAS) technique, which was first developed by Beley et al. (1995) and later expanded upon by Galushko et al. (2003) is commonly used. Given angular deflections of angle of arrival, delay, and Doppler measurements it is possible to determine the propagation characteristics when the ionosphere is treated as a reflecting mirror. This method has typically been applied using a single HF transmitter. Ionosondes have also been used extensively to understand the temporal evolution of isodensity contours due to the passing of a TID (Morgan et al., 1978; Tedd & Morgan, 1985; Galushko et al., 2003; Amorim et al., 2011; Cervera & Harris, 2014; Reinisch et al., 2018; Emmons et al., 2020). Traditionally, single-station ionosondes were used, however, more recent studies have started using networks of ionosondes for better spatial resolution (Reinisch et al., 2018). HF techniques are able to probe the electron density of the bottom-side ionosphere in contrast to the integrated measurements from GPS TEC and airglow imagers. Therefore, in order to better understand the 2-D spatial and temporal variations of the bottom-side ionosphere additional coverage from HF instrumentation is necessary.

Recently, Kaeppler et al. (2022), from here on referred to as K22, demonstrated that the Coastal Ocean Dynamics and Applications Radars (CODARs) could be used as bistatic single-frequency oblique ionospheric sounders. The frequency modulated continuous wave (FMCW) waveforms of these radars can be used to extract group delay measurements from which the virtual height of an ionospheric layer can be estimated. The cadence of observations is approximately 2 min for this investigation. An advantage of using this technique is the potential increase in coverage in areas with limited ground-based infrastructure. This is made possible by the widespread locations of the CODAR transmitters, the mobility of the individual receivers, and their relatively low cost (K22).

The purpose of this investigation is to expand upon the previous results from K22 by including more CODAR transmitters to estimate TID propagation speed and direction. We present a case study from 06 October 2020, in which we observe fluctuations in the ionospheric virtual height over a large spatial region. From these fluctuations, we derive horizontal velocities and propagation directions of the disturbance using bistatic CODAR observations. We demonstrate a generalized approach to include n stations rather than the standard three station triangle, which has been used extensively in several other studies (Mitra, 1949; Ratovsky et al., 2008; Medvedev et al., 2013; Aksonova et al., 2024). The advantage of using the CODAR systems is the increased coverage, which can be added to already existing networks of ionosondes, such as those described by Reinisch et al. (2018). To validate our findings we then compare our results to those derived using two Digisonde 256 digital ionosondes and one Digisonde-Portable-Sounder-4D (Reinisch & Galkin, 2011) located near the estimated ionospheric reflection point.

2 Methodology

2.1 Measurement techniques

In previous work done by K22, the authors showed several events with oscillations in virtual height; however, they did not attempt to estimate any TID propagation parameters. Here we use their same method for estimating the ionospheric virtual height from the CODAR observations. As noted in K22, the CODAR transmitters are co-channel but separated by different time delay offsets relative to GPS PPS. We refer to the time delay that includes the GPS offset as the pseudo-time delay or if including the speed of light, the pseudo-group range. This offset can be estimated by using the virtual height associated with the relatively stable E-region reflection. Using this method, we are able to calibrate out the GPS offset and produce the group delay. We used a standard matched filter of the known transmitted waveform template, which is cross-correlated with the received signal, that was measured using an Ettus USRP N200 with GPSDO. Details regarding the data acquisition, calibration methodology, and basic signal processing can be found in K22.

In order to first extract the time series of the pseudo-group range, we used the SciPy peak finding algorithm (Virtanen et al., 2020) for both O (left-hand circular polarization in the northern hemisphere) and X (right-hand circular polarization in the northern hemisphere) modes of propagation. The incoming elliptically polarized waves were decomposed into O and X modes based on the calculated Stokes parameter from the IQ data; see K22 for more details.

Given that the majority of our ground distances exceed 600 km, the virtual height, h, can then be calculated by assuming mirror-like reflection of our signal at the midpoint between transmitter and receiver, using the following curved Earth approximation derived from Davies (1990) and compared with Fabrizio (2013).

$$h=\sqrt{{\left(\frac{P}{2}\right)}^2-\left(a^2{\sin }^2\frac{D}{2a}\right)}-a\left(1-\cos \frac{D}{2a}\right)$$(1)

where P is the calculated group range (P = cΔT, where c is the speed of light and ΔT is the time delay of ionospheric reflection) and D is the great circle distance between the transmitter and the receiver. a is the radius of the Earth.

2.2 CODAR transmitters

The CODAR transmitters are located along the Eastern and Western U.S. Coasts, and the Gulf of Mexico. These transmitters are typically grouped by region and in most cases transmission frequency. In this study we used four groups of CODAR stations that were selected based on data quality, geographic location, and transmission frequency. Due to the location of our receiver, from here on referred to as “CARL” (black star in Fig. 1), we selected CODARs located along the Eastern Coast of the U.S. with a transmission frequency of approximately 4 MHz. Using nearly identical transmission frequencies, each with very similar equivalent vertical frequencies, allows us to assume that we are probing similar isodensity contours when using the transmitters located in New Jersey, North Carolina, and Florida. Within each of the four groups, there were multiple transmitters with very similar midpoint locations; therefore, the station with the strongest return signal for each group was used in this study. The four stations selected have the following callsigns: BRIG, CORE, NAPL, and PINS. Table 1 shows the frequency, geographic location, distance to CARL, and distance from the CODAR-CARL midpoint to the BRIG-CARL midpoint. It should be noted that all ground distances presented in this paper were calculated using the Haversine formula to account for the curvature of the Earth. We additionally use data from several ionosondes, including Wallops Island, VA (WP 937), Eglin Air Force Base, FL (EG 931), and Austin, TX (AU 930). These stations were chosen based on the great circle CODAR-CARL midpoint, as seen in Figure 1 by the blue squares.

Figure 1 A map of each CODAR and ionosonde station used in the study. The black dots represent the CODAR transmitter locations, and the X’s represent their respective estimated ionospheric reflection point. The dotted line is the assumed path of propagation and the star labeled “CARL” is the location of the receiver for all CODAR transmitters. The blue squares represent the three ionosonde stations used.

2.3 Estimation of the period of oscillation

To determine the periods of oscillation in the ionospheric virtual heights, a Lomb-Scargle periodogram (Lomb, 1976; Scargle, 1982) was used as shown in Figure 2. We evaluated the Lomb-Scargle for periods between 5 min and 8 h on the CODAR-derived virtual heights after applying a Gaussian filter. The inputted data ranged from 0100 UTC – 0900 UTC (2100 – 0500 Local Time (LT) for BRIG, CORE, and NAPL and 2000 – 0400 LT for PINS). These hours were selected in order to avoid any possible effects caused by the passing of the terminators. We note that, generally, we have observations only during local night due to D-region absorption during the day, where a link between the CODAR and CARL cannot be established. An example of our methodology for one CODAR-CARL link on 06 October 2020 is shown in Figure 2. Figure 2A is the initial periodogram, showing a dominant 8-hour trend. To remove this trend, an 8-hour period was used to reconstruct the inputted waveform (black dots in Fig. 2B) by fitting the parameters of a sinusoidal expansion (red curve in Fig. 2B). Figure 2C shows the residual virtual heights after the background 8-hour trend has been subtracted from the original virtual heights, leaving only the smaller periods of oscillation (black dots). The red line shows the reconstructed waveform, as described above, using all dominant periods as the input frequencies. Figure 2D is the Lomb-Scargle periodogram calculated again using the same frequency input on the residual virtual heights. Figure 2D also shows that only the peaks in the normalized power spectral density (PSD) with false alarm probabilities (VanderPlas, 2018) much less than 0.5% were considered. While we showed the example from the BRIG-CARL link, similar Lomb-Scargle results, with an approximately 2.5-hour period peak, were observed across all four stations.

Figure 2 Here we show the process as described in Section 2.3 for the BRIG transmitter only. Panel A is the Lomb-Scargle periodogram performed on the CODAR-derived virtual heights. Panel B is the CODAR-derived virtual heights (black dots) with the reconstructed waveform based on the dominant 8-hour period of oscillation plotted over top (red curve). Panel C is the CODAR-derived virtual heights with the dominant 8-hour trend removed and the reconstructed waveform plotted over top (red curve). Panel D is the Lomb-Scargle periodogram based on the data shown in panel C. The dashed line refers to a significance level of 0.01.

2.4 Estimation of time delay

After initial inspection, we determined that the disturbance appears to impact the BRIG station first; therefore, our t₀ value is based upon the BRIG station. t₀ here represents the timestamp of our structure crossing the BRIG-CARL midpoint, which is then used to find the time lag between each remaining station. This time lag was determined by calculating the maximum cross-correlation value between each station and BRIG (Rees & Mobbs, 1988), using data sets with the 8-hour trend removed. This process is similar to that used by Morgan et al. (1978) and Tedd & Morgan (1985). The same process was done for the three ionosondes using WP 937 as the t₀. For the calculations, including all seven stations, BRIG was used as the t₀ for all CODAR and ionosonde stations.

2.5 Determination of TID velocity and uncertainty quantification

Once the time lag was determined, the horizontal velocity was calculated from first principles using the assumed ionospheric reflection points and the calculated time delay from the cross-correlation analysis. It should be noted that the ionospheric reflection point is assumed to be the midpoint along the great circle path between each CODAR transmitter and the receiver, as marked in Figure 1. This methodology of determining the velocity is a similar process to several past investigations (Mitra, 1949; Nappo, 2002; Ratovsky et al., 2008; Aksonova et al., 2024). However, the method we develop here is generalized to include n stations, instead of the more typical 3-stations in a triangle method. Previously, Chum et al. (2014) has shown a different method for deriving velocities using multiple stations by calculating the slowness vector.

For the technique used in this study, it was assumed that we observed a single coherent planar structure moving across the four midpoints. This is justified because the correlation coefficient of each signal relative to BRIG exceed 0.5. The correlation coefficients are as follows: BRIG & CORE = 0.928, BRIG & PINS = 0.619, and BRIG & NAPL = 0.720.

We derived the speed and direction of propagation through the following model. We first consider the previous 3-midpoint example as a starting point, in which we define that the total delay we expect to observe between two stations can be written as,

$$\begin{array}{l}{t}_{12}=\frac{x_{12}}{V_x}+\frac{y_{12}}{V_y}=x_{12}{S}_x+y_{12}{S}_y,\\ t_{13}=\frac{x_{13}}{V_x}+\frac{y_{13}}{V_y}=x_{13}{S}_x+y_{13}{S}_y.\end{array}$$(2)

The t₁₂ here corresponds to the time lag between the midpoints of stations 1 and 2, where midpoint 1 is considered the reference location, t₁₃ follows a similar definition. The x₁₂ and y₁₂ correspond to x₁ – x₂ and y₁ – y₂, respectively, which are the distances in the x and y directions between station midpoints. The velocity components are represented by V_x and V_y, with the slowness being defined as S_(x/y) = 1/V_(x/y). Finally, positive x and y in this configuration correspond to eastward and northward directed, respectively.

We can rewrite equation (2) in the following matrix form as,

$$\left[\begin{array}{l}{t}_{12}\\ t_{13}\end{array}\right]=\left[\begin{array}{ll}{x}_{12}& y_{12}\\ x_{13}& y_{13}\end{array}\right]\left[\begin{array}{l}{S}_x\\ S_y\end{array}\right]$$(3)

which can be generalized as,

$$\mathbf{T}=\mathbf{AS}$$(4)

where A is a known matrix given the locations of the midpoints, T is the calculated time-delay vector, and S is the slowness vector. Solving equation (3) is straightforward, provided that the matrix A is not singular. We can solve for the velocities and we obtain

$$\begin{array}{l}{S}_x=\frac{1}{V_x}=\frac{y_{13}{t}_{12}-y_{12}{t}_{13}}{x_{12}{y}_{13}-y_{12}{x}_{13}}\\ S_y=\frac{1}{V_y}=\frac{x_{12}{t}_{13}-x_{13}{t}_{12}}{x_{12}{y}_{13}-y_{12}{x}_{13}}\end{array}$$(5)

which is consistent with equation (1) in Aksonova et al. (2024).

The problem can now be generalized for n stations which can be solved in matrix form as,

$$\left[\begin{array}{ll}{x}_{12}& y_{12}\\ x_{13}& y_{13}\\ ⋮& ⋮\\ x_{\mathrm{jn}}& y_{\mathrm{jn}}\end{array}\right]\left[\begin{array}{l}{S}_x\\ S_y\end{array}\right]=\left[\begin{array}{l}{t}_{12}\\ t_{13}\\ ⋮\\ t_{\mathrm{jn}}\end{array}\right]$$(6)

where j ≠ n. Equation (6) can be recast into the form of a matrix, i.e., equation (4). This corresponds to an overdetermined problem and can be solved using linear least squares, where the optimal solution can be written as,

$$\mathbf{S}={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathbf{T}$$(7)

which can be solved explicitly as,

$$V_x=\frac{\left(\sum _{i=2}^n\mathrm{}{x}_{1i}^2\right)\left(\sum _{i=2}^n\mathrm{}{y}_{1i}^2\right)-(\sum _{i=2}^n\mathrm{}{x}_{1i}{y}_{1i}{)}^2}{\left(\sum _{i=2}^n\mathrm{}{y}_{1i}^2\right)\left(\sum _{i=2}^n\mathrm{}{x}_{1i}{t}_{1i}\right)+(-\sum _{i=2}^n\mathrm{}{x}_{1i}{y}_{1i})\left(\sum _{i=2}^n\mathrm{}{y}_{1i}{t}_{1i}\right)}$$(8)

and

$$V_y=\frac{\left(\sum _{i=2}^n\mathrm{}{x}_{1i}^2\right)\left(\sum _{i=2}^n\mathrm{}{y}_{1i}^2\right)-(\sum _{i=2}^n\mathrm{}{x}_{1i}{y}_{1i}{)}^2}{\left(\sum _{i=2}^n\mathrm{}{x}_{1i}^2\right)\left(\sum _{i=2}^n\mathrm{}{y}_{1i}{t}_{1i}\right)+(-\sum _{i=2}^n\mathrm{}{y}_{1i}{x}_{1i})\left(\sum _{i=2}^n\mathrm{}{x}_{1i}{t}_{1i}\right)}.$$(9)

The horizontal phase velocity (V_h), hereafter referred to as horizontal velocity, and the propagation direction (ϕ) are therefore,

$$V_h=\frac{V_x{V}_y}{(V_x^2+V_y^2{)}^{1/2}}\to \varphi ={\tan }^{-1}\left(-\frac{V_y}{V_x}\right)+180°$$(10)

which, as stated above, is consistent with the previous results from Aksonova et al. (2024). We note that several texts discuss the inversion of linear least squares problems, both from a frequentist perspective (Aster et al., 2013) and a Bayesian approach (Tarantola, 2005; Heinselman & Nicolls, 2008). These methods can be additionally applied to better condition the matrix A, but exploring this topic is outside of the scope of this investigation.

The uncertainty in both the velocity and propagation direction arises from the time delay of our signals. Therefore, to calculate the total uncertainty, a Monte Carlo simulation was performed. First, a Gaussian function was fit to the largest peak of each cross-correlation output, and its standard deviation was then calculated. This value was used to define the range of random noise added to the time delay, and the process was repeated 10,000 times for each required transmitter pair. From this, each method of velocity and propagation direction estimation was then performed on each output of the time delay. The errors shown in Tables 2 and 3 represent the 68% confidence interval of those velocity and propagation direction estimates.

3 Results

Here we present the results from a case study of 06 October 2020, selected based on the October 2020 data set first examined by K22. As an extension to that previous study, we show several time series of virtual heights from each of the transmitters discussed in Section 2.2 for the hours of 0100 – 0900 UTC.

Figure 3A shows the CODAR-derived virtual heights color-coded by station and therefore transmission frequency. For each stations the averaged virtual height error was calculated using,

$$\Delta h^2=\frac{\Delta P^2}{(4h/P{)}^2}$$(11)

where ΔP is the uncertainty of the group range, found to be 11.7 km (K22), and h is the calculated virtual height. The uncertainty in the ground distance is assumed to be zero and therefore is neglected in this calculation. This uncertainty is represented in Figure 3 as the grouping of points located in the upper left-hand corner of Figure 3A. Each point is color-coded by the corresponding signal for which the error was calculated for. It is important to note that this panel is showing the CODAR-derived virtual height traces before the 8-hour period of oscillation was removed. In Figure 3, it is clear to see the slight time delay between each CODAR trace that was discussed in Section 2.4. This example is best shown during the hours of 0300 – 0500 UTC.

Figure 3 A time series of virtual height oscillations from October 06, 2020. Panel A shows the four CODAR stations: BRIG, CORE, NAPL, and PINS color coded by transmission frequency. Panel B is the Digisonde station WP 937, located in Virginia. Panel C is the Digisonde station AU 930, located in Texas, and Panel D is the Digisonde station EG 931, located in Florida. All virtual height traces shown in each panel are for O-mode propagation only.

For validation of our results, we also present a time series of vertical incidence Digisonde-derived virtual heights for the same time span as the CODAR measurements. This is shown in the bottom three panels of Figure 3. Figure 3B is the virtual heights derived from WP 937. Figure 3C shows the virtual heights derived from AU 930. Figure 3D shows the virtual heights derived from EG 931. For all three stations, each line is color-coded by transmission frequency, ranging from 2.5 MHz to 4.00 MHz.

It is important to note that the frequencies selected for the comparison between CODAR-derived horizontal velocities and Digisonde-derived horizontal velocities were chosen based on the calculated equivalent vertical frequency of each CODAR station, using the modified secant law adapted from Davies (1990).

$$k{f}_v\sec \varphi =f_o$$(12)

f_o is the frequency observed at the receiver, and ϕ = sin⁻¹d/P with d being the length of the chord between transmitter and receiver and P is the derived group range as discussed in Section 2.1. k is the correction factor used to account for curvature of both the Earth and ionosphere; in this case, a correction factor of 1.1 is used (McNamara, 1994). Thus, the average equivalent vertical frequency calculated across the hours of 0100 – 0900 UTC is as follows: CORE = 2.85 MHz, BRIG = 2.17 MHz, NAPL = 2.29 MHz, and PINS = 1.42 MHz. Therefore, only the following ionosonde frequencies were used in the calculations shown in Table 2, AU 930 = 1.80 MHz, WP 937 = 2.65 MHz, and EG 931 = 2.65 MHz. The frequencies were selected based on their proximity to the calculated equivalent vertical frequencies of the four CODAR transmitters as well as the data quality of the ionosonde; for example, a frequency of 2.65 MHz was chosen for WP 937 because it is approximately the average frequency between the BRIG (2.17 MHz) and CORE (2.85 MHz) stations. Table 2 shows three different outcomes of equation (10), each with a different input. The top row labeled CODAR, is the horizontal velocity and propagation direction of the presumed TID, only using the four CODAR measurements. The middle row labeled Digisonde is the horizontal velocity and propagation direction of the presumed TID only using the three Digisonde stations. Finally, the third row labeled CODAR + Digisonde, shows the calculated speed and direction using the four CODAR stations in combination with the three Digisonde stations for a total of seven stations. Using the four CODAR stations the horizontal velocity was found to be 580 m/s, the three Digisonde stations gave a horizontal velocity of 674 m/s, and the CODAR plus Digisonde stations gave a horizontal velocity of 555 m/s. Along with propagation directions for the three examples of 178°, 180°, and 184° respectively. The propagation direction is assuming North is 0°. This further supports our original theory that this structure is an LSTID traveling equatorward.

For completeness, we also show the Ionogram-derived parameters of the peak height of the F2-layer (hmF2) and the critical frequency of the F2 layer (foF2) as a function of time (Fig. 4), for each of the three Digisonde stations used. Note that these values are from automatically scaled ionosonde data obtained from the GIRO online database (Reinisch & Galkin, 2011). It can be seen that the foF2 values remain relatively consistent throughout the night from 0200 UTC forward. However, the hmF2 values follow similar trends in oscillation as the CODAR and Digisonde measurements. This could suggest that the oscillations we are seeing in virtual height could be caused by a layer height change rather than a density change in the F2 layer.

Figure 4 Ionogram-derived parameters of foF2 (left panels) and hmF2 (right panels) for the night of October 06, 2020 at stations WP 937, AU 930, and EG 931.

It is important to note that, similar to K22, we see the same oscillations in virtual height across all three Digisondes that we do from the four CODAR observations, which provides further evidence that the CODAR stations can be used as oblique ionospheric sounders. This correlation can be clearly seen during the first large dip in virtual height across all stations between the hours of 0400 UTC and 0500 UTC. This can also be seen with the two smaller oscillations during the hours of 0500 UTC and 0700 UTC. However, we do still see some discrepancies in virtual height oscillation in both height and time. This is specifically seen during the hours of 0200 – 0300 UTC for the 3.4425 MHz trace of AU 930 and the 4.000 MHz trace of EG 931. The first is likely due to the operating frequency of the ionosonde at that time being very close to the critical frequency of the F2 layer, as seen in Figure 4. The EG 931 trace may have encountered a similar issue; however, this could also be due to the data quality at that frequency.

4 Discussion

4.1 Calculation comparison

To further support this technique and the derived numerical values of horizontal velocity and propagation direction, we also compare our results with two other techniques that use n = 3 stations. The first method is discussed in several studies, including Mitra (1949), Ratovsky et al. (2008), and Aksonova et al. (2024). These investigations all have a similar setup to the 3-midpoint technique described in Section 2.5, and use a slightly modified version of equation (5) for their initial setup. The second method uses more of a trigonometric approach and is discussed in Section 7.9 of Davies (1990). Both sets of equations can be found in Appendix.

For simplicity, we will only be showing the comparison of methods using three CODAR stations, BRIG, NAPL, and PINS. Based on the values shown in Table 3, it can be seen that using our method for n = 3 stations yields results very similar to those from the other two methods. These values are also within approximately 2 m/s of the horizontal velocity calculated using our new generalized method on all seven stations and within 5 degrees for propagation direction, further validating our technique.

4.2 Possible drivers

From the accepted categorization of TIDs by Georges (1968) and Hocke & Schlegel (1996), the derived values of horizontal velocity, propagation direction, and period of oscillation would suggest that this disturbance is an LSTID traveling equatorward from higher latitudes. Therefore, a possible driver of this event could be geomagnetic activity as suggested by several other studies (Lyons et al., 2019; Frissell et al., 2022). We test this hypothesis by examining the geomagnetic activity level preceding the day of the observations. Figure 5 shows the OMNI 2 three-hour Kp values (King & Papitashvili, 2005) and the SuperMAG electrojet SME index (Gjerloev, 2012) with a one-minute sampling rate for 5 days before and 7 days after 06 October 2020. From Figure 5, it is clear across both indices that there is an increase in both Kp and SME indices the day prior to this event, highlighted in yellow, with values greater than 4 and 810 nT for Kp and SME respectively. Beginning on October 3, 2020 there is a decrease in both Kp and SME indices, until the sharp increase late in the day on October 5, 2020. After this event, both indices are quiet again.

Figure 5 Top: Kp values for Oct 01, 2020 – Oct 13, 2020. Bottom: SME index for Oct 01, 2020 – Oct 13, 2020.

We also consider ground-based magnetometer outputs from the SuperMAG database (Gjerloev, 2012), beginning on October 5, 2020, at 0000 UTC and ending on October 6, 2020, at 1000 UTC. Based on the calculated propagation direction of our disturbance, we choose to investigate the magnetometer located at Rankin Inlet (RAN) in Nunavut, Canada (62.82°N, 92.11°W). It is important to note that SuperMAG reports its data in a coordinate system decomposed into north, east, and vertically down components (N, E, and Z), with E being similar to the standard H component. For simplicity, we will continue using this naming convention. More information about the coordinate system SuperMAG uses can be found here: https://supermag.jhuapl.edu/info/. Based on Figure 6, we observe very little fluctuation in each of the components until 0500 UTC on October 5, 2020. From that point on, we see small fluctuations in each component until 2000 UTC on October 5, 2020. Beginning at 2000 UTC, we observe very large fluctuations in all components until approximately 2300 UTC, indicating an increase in geomagnetic activity. This first event could have been the cause of our recorded virtual height fluctuations; however, based on our calculated horizontal velocities, this is rather unlikely. At approximately 0100 UTC on October 6, 2020, we observe another onset of activity, suggesting another increase in geomagnetic activity, though slightly weaker. Therefore, based on our calculated horizontal velocity of approximately 555 m/s, the timing of the increases in Kp and SME, and the fluctuations in magnetic field components, it is quantitatively possible that this geomagnetic activity was the driver of the oscillations in virtual height beginning at approximately 0300 UTC. As another visual representation of our event, we show Figure 7, which is an estimation of the wavefront locations every 20 min, assuming our t₀ value is the time the wavefront crosses the BRIG-CARL link. This estimation assumes that our wavefront is traveling equatorward at speeds of approximately 550 m/s, and we assume that the phase and group velocities have similar propagation directions and speeds. However, a further modeling investigation of the TID phase and group velocities is outside the scope of this paper.

Figure 6 Magnetometer output from RAN station. Starting at 0000 UTC on October 05, 2020 and going until 1000 UTC on October 06, 2020. The portion highlighted in yellow is possible cause of our presumed LSTID.

Figure 7 Estimated locations of a wavefront originating from high latitudes and propagating equatorward. The black dots represent the CODAR transmitters, the black X’s are the assumed ionospheric reflection point, the blue squares are the ionosonde stations, and the black star is the receiver. The red triangle is the magnetometer station “RAN”.

4.3 Limitations

One of the main limitations of this technique is that we are considering virtual height rather than real height, as done in Emmons et al. (2020). This could cause some difficulties in the interpretation of our data due to our inability to distinguish between the oscillations in Figure 3 being height variations or density variations. Additional observations from redline all-sky imagers and/or GPS TEC would provide potential assistance to understand this further; however that investigation remains outside the scope of this paper. Also, we were limited in the transmission frequencies we could use for this study since we have no control over the transmission frequencies or locations of the CODARs. This caused us to only consider nighttime disturbances because the lower frequencies are absorbed more strongly by the D-region during the day.

5 Summary and conclusion

The purpose of this investigation was to show that CODARs can be used as an additional network of transmitters that can augment existing ionosonde networks for the estimation of the propagation velocity and direction of traveling ionospheric disturbances. We presented a case study from 06 October 2020 comparing CODAR-derived virtual heights from four stations with large spatial distribution to three Digisonde-derived virtual heights at similar geographic locations. We generalized the previous 3-station method for the estimation of TID direction to a more general approach that recasts this problem as a linear least-squares estimation problem. Additionally, we present uncertainty estimates for the velocity parameters. Given the more generalized nature of our derivation, it is possible to solve the linear least squares problem both in a frequentist and Bayesian approach to obtain the velocity and direction of propagation. We have also compared our generalized technique to two other previously accepted methods to show this technique gives nearly the same values for two TID parameters, providing validation of our technique. In conclusion, our results have shown the capabilities of the HF oceanographic radars to capture a presumed TID and estimate its propagation speed and direction.

Acknowledgments

We gratefully acknowledge the SuperMAG collaborators (https://supermag.jhuapl.edu/info/?page=acknowledgement). We also wish to thank the creators of https://github.com/breid-phys/dgsraw for writing dgsraw, which was used to process the digisonde data. As well as, Lowell digisonde for collecting ionosonde data as part of the Global Ionospheric Radio Observatory (GIRO), available at https://doi.org/10.5047/eps.2011.03.001. We also acknowledge use of NASA/GSFC’s Space Physics Data Facility’s OMNIWeb service, and OMNI data. The OMNI data were obtained from the GSFC/SPDF OMNIWeb interface at https://omniweb.gsfc.nasa.gov. The editor thanks Sergey Fridman and an anonymous reviewer for their assistance in evaluating this paper.

Funding

This work is supported by Office of Naval Research (grant no. N0001424C2201) with a contract to Clemson University.

Conflicts of interest

The authors declare no conflict of interest.

Data availability statement

Processed and analyzed CODAR observations from CARL that were used in this investigation can be found at https://doi.org/10.5281/zenodo.16904730. The data for AU 930, EG 931, and WP 937 digisondes can be found at https://data.ngdc.noaa.gov/instruments/remote-sensing/active/profilers-sounders/ionosonde/mids12/.

Appendix

Equations used in the method comparison:

Aksonova et al. (2024) equation 1

$$\begin{array}{cc}{u}_x=\frac{x_{12}{y}_{13}-x_{13}{y}_{12}}{y_{13}{\tau }_{12}-y_{12}{\tau }_{13}},& u_y=\frac{x_{12}{y}_{13}-x_{13}{y}_{12}}{x_{12}{\tau }_{13}-x_{13}{\tau }_{12}}\to V_h=\frac{u_x{u}_y}{\sqrt{u_x^2+u_y^2}}.\end{array}$$(A1)

Davies (1990) equations (7.71) and (7.72)

$$\tan \varphi =-\frac{V_{\mathrm{EU}}\cos {\mathrm{\Psi }}_{\mathrm{EU}}-V_{\mathrm{ES}}\cos {\mathrm{\Psi }}_{\mathrm{ES}}}{V_{\mathrm{EU}}\sin {\mathrm{\Psi }}_{\mathrm{EU}}-V_{\mathrm{ES}}\sin {\mathrm{\Psi }}_{\mathrm{ES}}}$$(A2)

$$\frac{{\sin }^2({\mathrm{\Psi }}_{\mathrm{ES}}-{\mathrm{\Psi }}_{\mathrm{EU}})}{V_H^2}=\frac{1}{V_{\mathrm{EU}}^2}+\frac{1}{V_{\mathrm{ES}}^2}-\frac{2\cos ({\mathrm{\Psi }}_{\mathrm{ES}}-{\mathrm{\Psi }}_{\mathrm{EU}})}{V_{\mathrm{EU}}{V}_{\mathrm{ES}}}$$(A3)

where ϕ is the propagation angle from North of the wavefront, Ψ_ES is the angle from North between stations E and S, and V_ES is the horizontal speed along ES (V_ES = ES/t_ES, t_ES = time delay between stations E and S). V_H is the horizontal velocity.

References

All Tables

All Figures

Figure 1 A map of each CODAR and ionosonde station used in the study. The black dots represent the CODAR transmitter locations, and the X’s represent their respective estimated ionospheric reflection point. The dotted line is the assumed path of propagation and the star labeled “CARL” is the location of the receiver for all CODAR transmitters. The blue squares represent the three ionosonde stations used.

In the text

Figure 2 Here we show the process as described in Section 2.3 for the BRIG transmitter only. Panel A is the Lomb-Scargle periodogram performed on the CODAR-derived virtual heights. Panel B is the CODAR-derived virtual heights (black dots) with the reconstructed waveform based on the dominant 8-hour period of oscillation plotted over top (red curve). Panel C is the CODAR-derived virtual heights with the dominant 8-hour trend removed and the reconstructed waveform plotted over top (red curve). Panel D is the Lomb-Scargle periodogram based on the data shown in panel C. The dashed line refers to a significance level of 0.01.

In the text

Figure 3 A time series of virtual height oscillations from October 06, 2020. Panel A shows the four CODAR stations: BRIG, CORE, NAPL, and PINS color coded by transmission frequency. Panel B is the Digisonde station WP 937, located in Virginia. Panel C is the Digisonde station AU 930, located in Texas, and Panel D is the Digisonde station EG 931, located in Florida. All virtual height traces shown in each panel are for O-mode propagation only.

In the text

	Figure 4 Ionogram-derived parameters of foF2 (left panels) and hmF2 (right panels) for the night of October 06, 2020 at stations WP 937, AU 930, and EG 931.
In the text

	Figure 5 Top: Kp values for Oct 01, 2020 – Oct 13, 2020. Bottom: SME index for Oct 01, 2020 – Oct 13, 2020.
In the text

	Figure 6 Magnetometer output from RAN station. Starting at 0000 UTC on October 05, 2020 and going until 1000 UTC on October 06, 2020. The portion highlighted in yellow is possible cause of our presumed LSTID.
In the text

	Figure 7 Estimated locations of a wavefront originating from high latitudes and propagating equatorward. The black dots represent the CODAR transmitters, the black X’s are the assumed ionospheric reflection point, the blue squares are the ionosonde stations, and the black star is the receiver. The red triangle is the magnetometer station “RAN”.
In the text