Tomographic imaging of a large-scale travelling ionospheric disturbance during the Halloween storm of 2003

The most intense ionospheric storm observed in recent times occurred between 29 and 31 October 2003. The disturbances to the high-latitude regions set off several large-scale travelling ionospheric disturbances (LSTIDs), wave-like perturbations in the ionospheric electron density. This paper investigates one particular TID on 31 October 2003 using North American Global Positioning System (GPS) receiver network data and a tomographic imaging technique. The TID has an estimated period of 30 min and an estimated horizontal wavelength of 700 km and propagates south-westward over North America. The tomographic reconstruction of the wave is validated using a simulation of the observations and with independent observations from ionosondes and the CHAMP planar Langmuir probe. The results are discussed in the context of the magnetic and ionospheric conditions that may have contributed to the launch of the wave. Large-scale TIDs are challenging to study over large regions of the Earth, and the GPS network here is shown to offer a unique perspective on the spatial and temporal variation of the TID. The experimental results are backed up by simulations that show a denser network of receivers, as is available in more recent years, would produce improved accuracy in the TID imaging.

1 Introduction

Travelling ionospheric disturbances (TIDs) are ionospheric manifestations of atmospheric gravity waves (AGWs) occurring in the neutral atmosphere (Hines, 1960). AGWs are buoyancy waves in the atmosphere and can be observed as TIDs when they transfer momentum to ions in the ionosphere by collision. Large-scale TIDs (LSTIDs) are a common occurrence during geomagnetic storms. LSTIDs are wave-like perturbations in the ionospheric electron density, with typical wavelengths over 1000 km and periods between 0.5 and 3 h (Hocke and Schlegel, 1996) and that typically travel equatorwards from the auroral regions (Davis and da Rosa, 1969). LSTIDs perturb the electron density and hence the total electron content (TEC), the number of free electrons along a path through the ionosphere, on scales of up to several TEC units (1 TEC_u = 10¹⁶ free electrons per m²). TEC is proportional to the first-order ionospheric delay of transionospheric radio waves propagating in the ionosphere and is therefore a crucial parameter for the Global Navigation Satellite System (GNSS).

Between 29 and 31 October 2003 a series of coronal mass ejections (CMEs) – expulsions of plasma from the solar corona – reached the magnetosphere of the Earth, causing strong geomagnetic storms. These are often referred to as the Halloween Storm(s) of 2003. The CMEs caused two sudden storm onsets on 29 October 2003 and 30 October 2003 (e.g. Mannucci et al., 2005; Horvath and Lovell, 2010). The planetary K-index (Kp) peaked at 9 on 29 and 30 October 2003 and at 8 on 31 October. Kp remained above 4 throughout 31 October, which, although still disturbed, constituted the recovery phase corresponding to the second sudden onset. The auroral electrojet index (AE) reached a maximum of 1827 nT at 06:31 UTC on 31 October, which is plotted in Fig. 1. Change in AE is related to auroral ionospheric current activity, which has been correlated with the appearance of TIDs at mid-latitudes (Hajkowicz, 1991; Hunsucker, 1982; Hocke and Schlegel, 1996; Lewis et al., 1996). These TIDs are thought to be launched by Joule heating of the atmosphere caused by increased ionospheric currents. High variability in AE occurred several times throughout 31 October, as seen in Fig. 1. This variability in AE provides evidence for a potential TID generation mechanism being present.

Figure 1 The auroral electrojet index at 1 min intervals on 31 October 2003.

LSTIDs during the first two days of the October 2003 ionospheric storms have been studied extensively (e.g. Afraimovich et al., 2006; Ding et al., 2007; Perevalova et al., 2008; Valladares et al., 2009; Borries et al., 2009; Horvath and Lovell, 2010). This study focuses on the less intense third day of the storms, 31 October 2003, and specifically on a high-amplitude TID observed over North America in the local morning hours (16:00–20:00 UTC).

Section 2 covers the data instrumentation used for the study of the TID and shows examples of the GNSS slant TEC (sTEC) observed. In Sect. 3.1, observations from different instruments and techniques – GPS tomography, an ionosonde and a space-borne planar Langmuir probe (PLP) – are compared. To investigate the effects of using a sparse network of GPS receivers, an additional tomographic inversion using simulated data is performed in Sect. 4. Section 5 contains a short discussion on the results and generation of the TID and final conclusions.

2 Data and instrumentation

The primary data used in this study were sTEC measurements derived from phase delay observations by a network of ground-based dual-frequency GPS receivers. In addition to the GPS sTEC used to image the TID, independent ionosonde data and measurements from the Challenging minisatellite Payload (CHAMP) PLP were used to confirm the presence of a TID.

2.1 GPS TEC

The GPS receiver network is shown in Fig. 2 and includes 40 stations in North America (listed in Table ) which are part of the International GNSS Service (IGS) and UNAVCO networks.

Slant TEC values were calculated using the geometry-free combination. It should be noted that MIDAS (Sect. 2.1.1) uses time-differenced sTEC measurements, so satellite and receiver biases which change slowly over time have no effect on the accuracy of the inversion (Mitchell and Spencer, 2003).

Figure 2 Network of GPS receivers used (circles) and the location of the Dyess and Millstone Hill ionosondes (triangles).

Figure 3 shows an example of pseudorange-calibrated sTEC observations from one receiver station, tono, where wave-like perturbations can be seen in the sTEC of several satellites. The satellites with the clearest TID signatures, PRNs 3 and 31, had ionospheric pierce points (IPPs) moving north. It should be noted that the movement of the satellites relative to a TID may result in distortions to the apparent TID, as it introduces a Doppler-like shift in the apparent period of the TID perturbations (e.g. Wan et al., 1997; Hernández-Pajares et al., 2006; Penney and Jackson-Booth, 2015). Bolmgren et al. (2020) showed, using simulations, that MIDAS has the capacity to correctly image LSTIDs without explicitly taking this effect into account.

Figure 3 Biased sTEC from GPS receiver station tono on 31 October 2003.

2.1.1 MIDAS

Computerised ionospheric tomography is a method that can estimate the 2D or 3D ionospheric electron density over an area using integrated electron density measurements, such as TEC. In general, ionospheric tomography can be described as solving an inverse problem formulated by the relationship between the geometry, the observations and the discretised electron density distribution. For a historical review of different methods of ionospheric tomography, see Bust and Mitchell (2008).

In this study, the electron density was imaged using the Multi-Instrument Data Analysis Software (MIDAS) tomography algorithm (Mitchell and Spencer, 2003). MIDAS uses differential phase observations from a network of ground-based geodetic GNSS receivers and solves for an estimate of the ionospheric electron density. Empirical orthogonal functions (EOFs) are used as a change in basis in the height dimension; this constrains the problem by decreasing the degrees of freedom and by providing a basic structure to the variation of electron density with height. MIDAS has previously been tested as a TID imaging algorithm using a simulation approach in Bolmgren et al. (2020), which established that the algorithm can successfully reproduce LSTIDs using GNSS data. In this study we will show that this is possible with real data even in relatively challenging conditions.

2.2 Ionosondes

The first scientific observations of TIDs were made using ionosondes (Munro, 1948). Ionosondes are ground-based radio instruments that characterise the bottomside electron density of the ionosphere. Ionosondes work by generating signal pulses that sweep through a span of frequencies. The pulses reflected back to the Earth from close to the zenith are used to estimate the height distribution of the plasma frequency, which is proportional to the square root of the electron density, directly above the ionosonde. The highest plasma frequency is usually found in the F2 layer and is denoted foF2. Since electromagnetic waves with frequencies above foF2 pass through the ionosphere, ionosondes provide no information on the electron density above the height of the F2 layer (referred to as hmF2).

Ionosondes at Dyess (32.4^∘ N, 99.8^∘ W) and Millstone Hill (42.6^∘ N, 71.5^∘ W) were both active on 31 October 2003. Figure 2 indicates the locations of these two ionosondes. The Millstone Hill ionosonde is used as a reference when setting up the MIDAS EOFs, while measurements from the Dyess ionosonde are used in Sect. 3.2.

2.3 CHAMP planar Langmuir probe

The CHAMP satellite was active for 10 years between 2000 and 2010 and was equipped with atmospheric and ionospheric observation instruments. CHAMP has a near-circular polar orbit and had an altitude around 390 km at the time of the storm, which usually would be in the topside of the ionospheric F layer. This study makes use of electron density data from the CHAMP PLP, a planar Langmuir probe which was used to measure in situ electron temperature as well as electron density in the front of the spacecraft every 15 s. Details on the CHAMP PLP can be found in McNamara et al. (2007).

3 Results

Section 3.1, 3.2 and 3.3 present the results in terms of the tomographic GPS inversion, foF2 and hmF2 from the Dyess ionosonde, and CHAMP PLP in situ electron density respectively.

3.1 Tomographic inversion

Differential phase observations from the GPS receiver network were used with MIDAS to estimate the ionospheric electron density distribution on 31 October 2003. The reconstructions in MIDAS used voxels of $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }\times \mathrm{10}$ km in latitude, longitude and height respectively and time steps of 10 min. Two EOFs were generated using a set of Chapman profiles (Chapman, 1931), adjusted to fit the vertical profiles observed by the Millstone Hill ionosonde.

Figure 4 Series of vTEC from the MIDAS GPS inversion. Each frame is separated by 10 min. Arrows have been added to indicate identified wave crests.

Figure 4 shows six consecutive time frames between 17:10 and 18:00 of the inversion results, with electron density integrated vertically to give vTEC. Between two and four wave fronts aligned NW–SE can be observed in the figure, spanning latitudes between 45 and 30^∘. These features are also visible in the electron density viewed as a cross section spanning 100–1200 km in altitude along the direction of travel, shown in Fig. 5. The wave-like perturbations are presumed to be the result of a passing TID.

Figure 5 (b) NE–SW cross section of the full inverted electron density and (a) the path of the placement of the cross section in the vTEC map.

Using consecutive tomographic images from MIDAS, the TID parameters were estimated as follows: horizontal wavelength λ_h≈700 km, phase velocity v_ph≈390 m s⁻¹, and direction of travel $\approx \mathrm{195}{}^{\circ }$ S–W. The period T was estimated as $T={\mathit{\lambda }}_{\mathrm{h}}/v_{\mathrm{ph}}\approx \mathrm{30}$ min. These parameters would qualify the TID as medium scale, following the definitions in Hunsucker (1982). However, considering the high amplitude, geomagnetic conditions and equatorward direction of travel, we will consider it a LSTID.

3.2 Ionosonde observations

The Dyess ionosonde is located within the area that was visibly affected by the TID in the MIDAS images. There is an indication of a periodical signature in the F2 layer critical frequency (foF2) with a 30 min period between 18:00 and 19:30 UTC, which may be related to the TID visible in the GPS data. However, the 15 min sampling makes it impossible to detect potential shorter-period perturbations. In Fig. 6, foF2 and hmF2 from the Dyess ionosonde are plotted against the equivalent parameters calculated from the MIDAS result. In Table 4 of Bruno et al. (2020), MIDAS results were compared against ionosonde data, and for a setup close to what is used here Bruno et al. (2020) found errors of 0.55 MHz in foF2 and 40 km in hmF2. The discrepancies in Fig. 6 are of the same order. The other ionosonde with data readily available during this period, Millstone Hill (42.6^∘ N, 71.5^∘ W), does not show a similar indication of TID passage. This is expected, since it is located outside of the area visibly affected by the TID in the tomographic inversions.

Figure 6 foF2 (a) and hmF2 (b) observations for the Dyess ionosonde (32.4^∘ N, 99.8^∘ W) and MIDAS equivalent sampled at the same location on 31 October 2003.

3.3 CHAMP PLP observations

The CHAMP satellite had one north-to-south pass over North America between 17:00 and 19:00 UTC on 31 October 2003, when the TID was visible in the GPS TEC. The in situ electron density measured by the PLP at altitudes between 391 and 395 km for this pass over North America is plotted in Fig. 7. At 17:43 UTC, CHAMP passes North America at longitude 76^∘ W, i.e. east of the area where the TID is visible in the tomographic images. Two dips in electron density separated by an apparent latitudinal wavelength of around 700–825 km are visible in Fig. 7 between latitudes 15 and 30^∘ N.

The dotted line in Fig. 7 shows the electron density estimated by MIDAS at 17:40 UTC, sampled at the location of CHAMP. Apart from not displaying the same wave perturbations, the electron density at this altitude is overestimated by approximately 3×10¹¹ electrons per m³. This is the result of a mismatch between the in situ observation and integrated estimate of the vertical density distribution in this area.

The perturbations in Fig. 7 may indeed be caused by the passage of the TID seen further west in the tomographic images, but poor receiver coverage in the region may explain why the wave fronts do not appear to reach 76^∘ W in the tomography result. The effect of possible poor data coverage is further examined by testing the tomography procedure on simulated data in Sect. 4.

Figure 7 In situ electron density (a) measured by the CHAMP planar Langmuir probe (line) and the electron density sampled from the MIDAS inversion (dotted line) for 31 October 2003. The corresponding CHAMP satellite track is plotted on top of the MIDAS vTEC result for 17:40 UTC (b).

4 Method verification by simulation

The Dyess ionosonde and CHAMP PLP electron density both suggest the presence of a TID, but the wave-like features observed by these instruments are not clearly translated onto the same spatial and temporal coordinates in the MIDAS inversion results. The ionosonde suggests the presence of a TID with a period similar to that in the GPS inversion, but it appears later than it does in the inversion. The CHAMP satellite measurements suggest wave-like perturbations can affect the electron density as high up as 390 km in a region where the wave is not visible in the tomographic inversion. It is possible that these features are not visible in the tomographic images due to poor receiver coverage below 30^∘ latitude.

Figure 8 The modelled LSTID vTEC (a) and inversions from the simulated data, (b) and (c), with a 1 h running mean background subtracted. Panel (c) uses a denser network than the real inversion, where the receiver locations of the network are indicated with black markers.

To investigate the effect of data coverage and geometry used for the tomographic inversion, simulated TEC from a model ionosphere was inverted with MIDAS under the same geometric conditions (satellite geometry and receiver coverage) as the original inversion. Any discrepancies between the model and simulated inversion results can be used to identify where there may be issues in the results presented in Sect. 3.1. A second inversion of the simulated data using a denser, fictional network is used to identify the effect of receiver geometry.

The TID parameters estimated in Sect. 3.1 were used together with the Hooke (1968) TID model and the International Reference Ionosphere, IRI2016 (Bilitza et al., 2017), to generate a model ionosphere with TID, through which sTEC measurements were integrated (following Bolmgren et al., 2020). A single frame of the model ionosphere is shown in Fig. 8a.

The resulting inversion shows that while the reconstruction with the regular network (Fig. 8b) is able to conserve the main morphology of the TID, it does not correctly replicate the perturbations of the wave east of 105^∘ W and south of 30^∘ N. In addition, the wave fronts in Fig. 8b appear skewed when compared to the model in Fig. 8a. In all panels of Fig. 8, a 1 h running mean was subtracted from each voxel post inversion to minimise the background ionosphere and to better see the TEC perturbations caused by the modelled TID.

The wave is more accurately reproduced if a denser network of GPS receivers than was available in 2003 is used. Figure 8c shows the improved simulation result, which uses a larger number of receivers. The simulated receiver network is marked by points in the same sub-figure. This inversion more accurately reproduces the perturbations in Fig. 8a, including the direction of the wave fronts.

Figure 9 OVATION Prime 2013 particle flux for 31 October 2003 at 14:30 UTC.

5 Discussion and conclusions

In this paper, we have used GPS tomography to reconstruct the ionospheric electron density over North America for 31 October 2003, the third day of the Halloween storm of 2003, and to identify a LSTID. The presence of a LSTID was evidenced by other instrumentation. A potential discrepancy in the TID morphology was observed between the measurements of two other instruments and the large-scale MIDAS reconstructions. While indications of the TID were captured by the Dyess ionosonde and CHAMP PLP, this was in areas where the MIDAS reconstruction showed no clear wave pattern. However, this was identified from computer end-to-end simulation to be the result of poor receiver coverage available for the MIDAS inversion, as discussed in Sect. 5. The receiver network used has an approximate receiver density of 1 per $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$, compared to approximately 6 per $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$ for the denser synthetic network shown in Fig. 8c. For comparison, the modern North American network used by Bruno et al. (2020) has an average receiver density close to 15 per $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$.

The observed TID had an estimated phase velocity of 390 m s⁻¹, an estimated period of 30 min, a horizontal wavelength of 700 km and a south-westerly direction, suggesting a source in the auroral region. The high variability in AE occurring between 11:00 and 14:00 UTC (Fig. 1) may indicate a possible time of launch of the observed LSTID if it were launched by Joule heating resulting from variations in the auroral electrojets. Another possible source mechanism may have been heating by auroral particle precipitation. The auroral oval was centred at latitude 63^∘ N at this time, with the region experiencing strong energetic particle precipitation at 14:30 UTC, as estimated by OVATION Prime 2013 (Newell et al., 2014) as shown in Fig. 9. The highest levels of precipitation around the presumed launch time of the LSTID occurred between 08:00 and 10:00 magnetic local time (MLT), which coincides with northern North America at the presumed launch time of the LSTID (11:00–14:00 UTC) and with the increased levels of AE activity around the same time. However, further analysis of additional datasets would be needed to obtain a detailed understanding of the generation mechanisms responsible for this LSTID. Since TIDs are effectively relative changes in the background electron density, the enhanced storm density likely contributed to the high perturbation amplitudes.

The work discussed in this paper is built on that of Bolmgren et al. (2020), where MIDAS was demonstrated to be capable of imaging certain TIDs, and has shown that the tomographic algorithm is capable of imaging LSTIDs with relatively small spatial dimensions, provided that a sufficiently dense ground receiver network is available.

Tomographic maps like the ones produced here could be used in practical navigation systems to provide ionospheric delay corrections for GNSS positioning during LSTID activity, which may otherwise induce unmodelled TEC fluctuations impairing the quality of the solution.

Appendix A: List of receiver stations

A1 North American GPS receiver stations used for the tomographic inversion

ID	Lat (deg)	Long (deg)
alrt	82.4939	−62.34179
bake	64.3178	−96.00235
bcov	50.5443	−126.8426
bogt	4.65553	−74.10725
cags	45.5851	−75.80731
chur	58.7591	−94.08873
cvms	35.5414	−89.64351
dsl1	70.3334	−148.4728
flin	54.7256	−101.978
frdn	45.9335	−66.65992
garl	40.4165	−119.3555
gtrg	43.2441	−113.2412
guat	14.5904	−90.52018
holm	70.7363	−117.7613
hvlk	37.6515	−99.10675
invk	68.3062	−133.527
kely	66.9874	−50.94485
kuuj	55.2784	−77.74545
mig1	34.0383	−120.3514
mkea	19.8018	−155.456
modb	41.9023	−120.3028
nain	56.537	−61.68864
ormd	29.2982	−81.10889
prds	50.8714	−114.2935
ptal	49.2563	−124.8609
qaq1	60.7153	−46.04779
reso	74.6911	−94.8961
ross	48.8337	−87.5196
sask	52.1963	−106.3984
sch2	54.8321	−66.83255
sg00	47.9218	−97.08662
thu2	76.537	−68.82508
thu3	76.537	−68.82508
tono	38.0972	−117.184
tukt	69.4383	−132.9944
vald	48.097	−77.56419
wslr	50.1265	−122.9211
yell	62.4809	−114.4806
ztl4	33.3797	−84.29673
will	52.2369	−122.1679