Polar tongue of ionisation during geomagnetic superstorm

During the main phase of geomagnetic storms, large positive ionospheric plasma density anomalies arise at middle and polar latitudes. A prominent example is the tongue of ionisation (TOI), which extends poleward from the dayside storm-enhanced density (SED) anomaly, often crossing the polar cap and streaming with the plasma convection flow into the nightside ionosphere. A fragmentation of the TOI anomaly contributes to the formation of polar plasma patches partially responsible for the scintillations of satellite positioning signals at high latitudes. To investigate this intense plasma anomaly, numerical simulations of plasma and neutral dynamics during the geomagnetic superstorm of 20 November 2003 are performed using the Thermosphere Ionosphere Electrodynamics Global Circulation Model (TIE-GCM) coupled with the statistical parameterisation of high-latitude plasma convection. The simulation results reproduce the TOI features consistently with observations of total electron content and with the results of ionospheric tomography, published previously by the authors. It is demonstrated that the fast plasma uplift, due to the electric plasma convection expanded to subauroral mid-latitudes, serves as a primary feeding mechanism for the TOI anomaly, while a complex interplay between electrodynamic and neutral wind transports is shown to contribute to the formation of a mid-latitude SED anomaly. This contrasts with published simulations of relatively smaller geomagnetic storms, where the impact of neutral dynamics on the TOI formation appears more pronounced. It is suggested that better representation of the high-latitude plasma convection during superstorms is needed. The results are discussed in the context of space weather modelling.

1 Introduction

In the course of a geomagnetic storm, large amounts of solar wind energy and momentum are deposited into the high-latitude ionosphere through the Joule dissipation of magnetosphere/ionosphere currents and auroral particle precipitation (Rodger et al., 2001). During the storm main phase a large positive dayside ionospheric plasma anomaly, known as the storm-enhanced density (SED), arises at sub-auroral mid-latitudes (Mendillo et al., 1972; Buonsanto, 1999; Immel and Mannucci, 2013). The morphologies of dayside SEDs have strong seasonal, local time, longitudinal, and other dependencies (Borries et al., 2015). A formation of the SED anomaly is largely (though not exclusively) attributed to the storm-time changes in plasma transport (Prölss, 1995, 2008; Immel and Mannucci, 2013), especially to the uplift of plasma to higher altitudes with longer recombination times. Storm-time changes in plasma/neutral composition and chemistry play greater roles in the formation of negative plasma anomalies, which are more common during the storm recovery phase (Rishbeth et al., 1987; Prölss and Werner, 2002).

The key physical mechanisms contributing to the storm-time plasma uplifts include (a) equatorward thermospheric neutral winds driven by the storm-time Joule dissipation (Anderson, 1976; Rishbeth, 1998) and (b) a vertical component of the electric E×B plasma convection expanded equatorward to mid-latitudes (Deng and Ridley, 2006; Heelis et al., 2009). Also, a horizontal plasma transport due to the poleward expansion of the equatorial plasma anomaly (Tsurutani et al., 2004) or due to the westward plasma drift caused by subauroral polarisation streams (Foster et al., 2007) has been invoked to explain the SED anomaly. However, the importance of the last two mechanisms, which involve substantial horizontal plasma transport over mid-latitudes, has been downplayed by Rishbeth et al. (2010) and Fuller-Rowell (2011), respectively, based on considerations of plasma transport times and global plasma density distributions. In this study the vertical uplifts due to (a) neutral winds and (b) expanded E×B convection are considered to be the key competing mechanisms for the generation of SED anomalies. We also note that this study does not aim to explain the formation of an SED anomaly.

The focal point of this study is the tongue of ionisation (TOI), which is a storm-time plasma density anomaly originating at the poleward edge of the SED anomaly, spreading anti-sunward across the polar cap and reaching the nightside auroral zone (Knudsen, 1974; Foster et al., 2005). The TOI anomaly has been observed during large geomagnetic storms using multiple radar systems (Foster et al., 2005) emphasising the role of cross-polar plasma transport by the enhanced E×B plasma convection flow. Using tomographic inversions of total electron content (TEC) observations, the 3D structure of the TOI anomaly has been revealed (Mitchell et al., 2008) and the role of dayside plasma uplift has been demonstrated (Yin et al., 2006). In situ satellite observations using ion drift instruments during the 20 November 2003 storm (Pokhotelov et al., 2008) suggested that the uplift can be attributed to the equatorward expansion of E×B convection flow. Storm-time observations of cross-polar plasma convection and plasma density using polar cap digital ionosondes (Pokhotelov et al., 2009) and SuperDARN radars (Thomas et al., 2013) demonstrated that sudden changes in the convection regime (e.g. due to rapid changes in the interplanetary magnetic field) can effectively disrupt the formation of a TOI anomaly. The fragmentation of a TOI anomaly is considered to be one of the mechanisms producing polar patches responsible for radio scintillations (e.g. Moen et al., 2013).

Earlier numerical simulations of the SED anomaly demonstrated competing roles of the plasma uplift mechanisms due to neutral winds and electric fields (e.g. Lin et al., 2005; Crowley et al., 2006; Swisdak et al., 2006). Since the mid-latitude SED anomaly provides a source of the uplifted dense plasma for the TOI anomaly, it is reasonable to assume that the same two mechanisms may control the formation of a TOI anomaly. However, the SED anomaly covers the entire local day–evening sector and often persists throughout the storm main phase and through an early part of the recovery phase, while the TOI anomaly is relatively narrow in longitude and persists for shorter times during the main phase. With recent developments of higher-resolution ionospheric circulation models (e.g. Maute, 2017), it became possible to simulate the dynamics of the TOI across the polar cap. Recently Liu et al. (2016) modelled the development of TOI anomalies during two comparable geomagnetic storms of March 2013 and March 2015 using the new release of the Thermosphere-Ionosphere Electrodynamics General Circulation Model (TIE-GCM) with 2.5^∘ horizontal resolution. Based on the simulations, they concluded that the uplift and horizontal transport due to the E×B drifts generally dominate over other possible drivers, such as neutral winds and compositional/chemical changes. Klimenko et al. (2019) modelled the TOI dynamics during the March 2015 geomagnetic storm, concluding that neutral dynamics and compositional changes may contribute to the suppression of a TOI anomaly beyond the geomagnetic North Pole. Using an ultra-high-resolution (0.6^∘) version of the TIE-GCM model, Dang et al. (2019) modelled the separation of the TOI anomaly into “double tongues” associated with morning and evening convection cells during the March 2013 storm event. However, these recent modelling studies simulated relatively moderate geomagnetic storms. The storm of March 2015, the largest in solar cycle 24, has a disturbance storm time (Dst) index minimum of −226 nT. Such storms are not in the category of “great storms”, commonly defined as having a Dst below $\sim -\mathrm{300}$ nT (Kamide et al., 1997). The current study is an attempt to model the TOI formation with a physics-based ionospheric model during a great storm (superstorm) event. The magnetosphere–ionosphere interactions in general (Kamide et al., 1997) and the formation of SED/TOI anomalies in particular (Pokhotelov et al., 2008) can both be quantitatively and qualitatively different during great storms.

In this study we use an example of the 20 November 2003 geomagnetic superstorm to analyse key mechanisms responsible for the formation and evolution of the TOI anomaly. The 20 November 2003 storm provides an advantage of being an isolated event driven by a single coronal mass ejection (Zhang et al., 2007). It is among the largest geomagnetic storms observed by modern space/ground instrumentation, including Global Navigation Satellite Systems (GNSS). Early studies of this storm using radars (Foster et al., 2005) and GNSS tomography (Pokhotelov et al., 2008) revealed the dynamics and 3D morphology of the TOI anomaly. However, self-consistent numerical simulations of the TOI anomaly were not possible at that time due to resolution limits of the existing ionospheric models and other factors. In this study the high-resolution version of the TIE-GCM model is used to model the TOI anomaly, with the analysis focusing on possible roles of the E×B drifts and neutral winds. A comparison with earlier GNSS tomography reconstructions is presented as well as with TEC distributions using conventional geometric TEC mapping. Limitations of other ionospheric circulation models in reproducing the TOI anomaly are also discussed, including the models currently used by space weather services.

2 Geomagnetic storm of 20 November 2003

In terms of the equatorial ring current disturbance magnitude, the geomagnetic storm of 20 November 2003 was the largest storm of solar cycle 23 and one of the largest storms recorded by modern instruments, with the Dst index reaching a value of −422 nT (Zhang et al., 2007). The storm was an isolated event preceded by a ∼ 20 d period of relatively quiet geomagnetic activity. Following the interplanetary shock arrival at 08:35 UT on 20 November 2003, the main phase of the storm lasts until ∼ 19:00 UT. During the main phase, the north–south IMF component (B_z) turns strongly negative, reaching to below −50 nT, while the dawn–dusk IMF component (B_y) increases to +50 nT in the beginning of the main phase and then goes down and turns negative after 18:00 UT (Fig. 1). With the solar wind speed (V_SW) exceeding 700 km/s, this IMF configuration should lead to a very strong two-cell plasma convection pattern. The observed IMF B_y change from positive to negative is expected to alter the east–west orientation of the “throat” (the entry region) of the cross-polar convection channel throughout the main phase (e.g. Sojka et al., 1994). During the main phase, the two-cell convection pattern expands dramatically to lower latitudes (to $\sim \mathrm{35}{}^{\circ }$ magnetic latitude), as also confirmed by in situ plasma drift measurements using the Defence Meteorological Satellite Program (DMSP) spacecraft (Pokhotelov et al., 2008). This expanded convection is expected to cause an anomalous vertical plasma transport at subauroral latitudes due to the resulting vertical component of E×B drift.

Figure 1 Solar wind speed (a), interplanetary magnetic field components (b), and symmetric horizontal component disturbance index (c) during the 20 November 2003 storm.

3 Simulations of the storm

To analyse the ionospheric dynamics during the 20 November 2003 storm, simulations have been performed using the TIE-GCM (Richmond et al., 1992; Qian et al., 2014). The TIE-GCM is a first-principle model simulating energy and momentum equations in the coupled thermosphere–ionosphere system. The current high-resolution version TIE-GCM v2.0 (Maute, 2017) uses the hydrostatic grid with 57 logarithmically spaced pressure levels ($\mathrm{1}/\mathrm{4}$ scale height resolution), covering geopotential heights from ∼ 97 to ∼ 600 km, with uniform horizontal 2.5^∘ grid resolution in longitude and latitude.

To facilitate the thermosphere/ionosphere forcing from above and below, the TIE-GCM should be coupled with external models. Mean horizontal neutral winds at the lower simulation boundary can be specified according to the Horizontal Wind Model (HWM07; Drob et al., 2008) and atmospheric tides are specified according to the Global Scale Wave Model (GSWM; Hagan and Forbes, 2002, 2003). Most relevant to the high-latitude storm dynamics, the plasma convection pattern is specified according to the statistical parameterisations of Heelis et al. (1982) or Weimer (2005). In this study, the Weimer parameterisation (Weimer, 2005) is used with the electrostatic potential expressed as a function of solar wind and IMF parameters measured upstream of the Earth's magnetosphere and time-shifted to the bow shock according to King and Papitashvili (2005).

The TIE-GCM simulation is performed throughout the 20 November 2003 storm after the 20 d initialisation run to reach the model equilibrium. Following the methodology of Liu et al. (2016), the simulated outputs of the 19 November 2003 quiet day are subtracted from the simulated 20 November 2003 storm day outputs. The resulting relative ΔTEC¹ anomalies for the 20 November 2003 day are shown as a snapshot at 15:00 UT in Fig. 2 and as an animated sequence for the interval 10:00–23:00 UT in the Supplements (movie01.avi). For reference, absolute values of TEC are also shown in Fig. 2.

Figure 2 Modelled TIE-GCM distributions of relative ΔTEC (a), plasma density at the F2 peak (b), absolute TEC (c) and the height of the F2 peak (d) at 15:00 UT 20 November 2003.

The following parameters relevant to storm-time plasma dynamics are also extracted from the TIE-GCM simulation. The height of the ionospheric F2 peak (hmF2) and plasma density at the F2 peak (NmF2) are shown in Fig. 2. Using the electrostatic potential (ϕ) given by the Weimer convection model, horizontal and vertical components of the plasma electric drift (${\mathbf{U}}_{E\times B}=-(\mathrm{\nabla }\mathit{\varphi }\times \mathbf{B})/B^{\mathrm{2}}$) are computed as vector products with the Earth's internal magnetic field (B). Expressed in geographic coordinates, northward meridional (V_E×B) and vertical (W_E×B) components of the electric drift are shown in Fig. 3 (a snapshot at 15:00 UT) and as an animated sequence for the interval 10:00–23:00 UT in the Supplements (movie02.avi). Meridional neutral winds (V) at pressure levels corresponding to ∼120 and ∼400 km geopotential heights and the Joule heating per unit mass (Q_Joule) are shown in Fig. 4 (a snapshot at 15:00 UT) and in the Supplements (movie03.avi).

Figure 3 Modelled TIE-GCM distributions of relative ΔTEC (a), electrostatic potential (b), and relative horizontal (c) and vertical (d) components of the E×B convection flow at 15:00 UT 20 November 2003.

Figure 4 Modelled TIE-GCM distributions of relative ΔTEC (b), Joule heating, and relative meridional neutral wind different levels (c)–(d) at 15:00 UT 20 November 2003.

Simulations of the 20 November 2003 storm with the Coupled Thermosphere Ionosphere Plasmasphere electrodynamics (CTIPe) model (Fuller-Rowell et al., 1996; Millward et al., 2001) have also been used in this study to compare to the TIE-GCM simulations described above. CTIPe is the first-principle model solving plasma and neutral dynamics on the hydrostatic grid with a resolution of $\mathrm{2}{}^{\circ }\times \mathrm{18}{}^{\circ }$ in latitude and longitude, respectively, and 15 pressure levels in the vertical direction going from the lower boundary at ∼ 80 to ∼ 400 km in altitude. The atmospheric forcing is specified according to the Whole Atmosphere Model (WAM) (Akmaev et al., 2008). WAM fields (neutral temperature, zonal and meridional neutral winds) are averaged in every local hour sector of a given month and thus contain the monthly-averaged mean winds and tides. The high-latitude electrodynamic forcing is specified according to the statistical parameterisations of Weimer (2005). The capability of CTIPe to reproduce SED anomalies during the main phase of the 20 November 2003 superstorm has been demonstrated in Fernandez-Gomez et al. (2019) for the European sector. In this study, the extended CTIPe run has been used to analyse anomalies in the North American sector, with the results discussed in Sect. 5.4. The purpose is to see whether the operational CTIPe model reproduces similar features to the higher-resolution research model (TIE-GCM).

4 Total electron content from GNSS mapping and tomography

The radio signals transmitted by GNSS can be used to retrieve information about ionospheric plasma anomalies. Due to dispersive properties of the ionospheric plasma, GNSS signals carry information about the TEC along the signal trajectory. Using thin ionospheric shell approximation (e.g. Horvath and Crozier, 2007) and taking proper care of the receiver and transmitter biases, slant TEC observations by ground GNSS receivers can be converted into the 2D distributions of vertical TEC. GNSS-based maps of TEC, using the thin shell transformation, are available from the International GNSS Services (IGS) with typical grid resolutions of $\mathrm{2.5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ in latitude × longitude (Hernández-Pajares et al., 2009). IGS TEC maps can be directly compared to the results of numerical simulations, though one needs to be careful with artefacts caused by sparse/uneven distributions of ground GNSS receivers. Examples of IGS TEC maps during the 20 November 2003 storm are presented in Fig. 5, also showing a comparison with the distributions of absolute TEC simulated by the TIE-GCM. An animated sequence of absolute TEC maps from TIE-GCM simulations and from IGS services for the interval 11:00–23:00 UT 20 November 2003 is included in the Supplements (movie04.avi).

Figure 5 TEC distributions obtained from TIE-GCM simulations (a, b) and from IGS services (c, d) at 16:00 UT (a, c) and 18:00 UT (b, d) during the 20 November 2003 storm.

A tomographic inversion of multiple slant TEC observations is also possible, yielding the 3D distribution of plasma density. The 3D time-dependent algorithm of ionospheric plasma tomography is described by Mitchell and Spencer (2003), and it has been previously applied to reconstruct the high-latitude plasma anomalies during the 20 November 2003 storm (Pokhotelov et al., 2008) using the network of 60 ground IGS receivers. Additional information about the E×B plasma drifts has been included in the tomographic algorithm using the Kalman filters with the Weimer convection model as a priori information (Spencer and Mitchell, 2007). The distributions of TEC obtained from the tomographic reconstructions were previously published and are presented in Fig. 6 of Pokhotelov et al. (2008) in the same format and at the same time moments (16:00 and 18:00 UT) as TEC distributions from TIE-GCM simulations and from IGS services shown here in Fig. 5.

5 Discussion

5.1 Total electron content

Total electron content maps provide global coverage showing the morphology of plasma anomalies on ionospheric mesoscales comparable to the horizontal resolution of TIE-GCM simulations presented here ($\mathrm{2.5}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$). However, the TEC mapping experiences potential problems at high latitudes due to (a) sparse/uneven distribution of ground GNSS receivers, (b) singularities of the latitude–longitude grid at the geographic poles, and (c) configuration of the GNSS satellite orbits. The network of ground IGS receiver stations available at polar latitudes during the 20 November 2003 storm is presented in Fig. 2 of Pokhotelov et al. (2008), showing separations between some of the polar cap receivers far greater than the desired horizontal resolutions. The inclination of GPS satellite orbits of $\sim \mathrm{55}{}^{\circ }$ (Samama, 2008) also contributes to the deficiencies of TEC reconstructions in the polar cap region. The tomographic reconstruction algorithm (Mitchell and Spencer, 2003; Spencer and Mitchell, 2007) partially mitigates these deficiencies by using rotated tomographic grids without the polar singularity and by including a priori information about plasma convection in the polar cap. Thus the tomographic reconstruction has advantages over the thin shell IGS TEC mapping, providing a more homogeneous solution across the polar cap region.

Taking into account the above limitations, we can compare the simulated TEC distributions with the results of TEC mapping and tomography. As shown in Fig. 5 and in the animation (movie04.avi), the TOI anomaly is visible in TEC maps (both from IGS and from the TIE-GCM model) starting from ∼ 12:00 UT, though the poleward extension of the SED anomaly appears at first some 20^∘ further westward in the TIE-GCM simulations relative to the IGS TEC maps (southern tip of Greenland in the simulations vs. east of Iceland in IGS maps). The reasons for this mismatch in the local time/location of the TOI formation are not clear and will be discussed further in Sect. 5.4. One has to note that TEC reconstructions are not reliable over the Atlantic Ocean sector due to poor GNSS receiver coverage.

The main development of the TOI anomaly (from ∼ 13:00 to 18:00 UT) is seen over the eastern sector of the United States–Canada, spreading further anti-sunward over the geomagnetic North Pole and northern tip of Greenland. In simulations and in TEC observations, the TOI anomaly develops in the same longitudinal sector (60–90^∘ W), though the simulated TOI appears narrower in longitude and more homogeneous in latitude. After 19:00 UT the TOI anomaly starts to disintegrate and disappear and the remains of the plasma are transported across the polar cap, merging into the nightside auroral TEC enhancement seen over the European sector. Overall, the location and general morphology of the simulated TOI anomaly is remarkably close to the IGS TEC observations and the tomography, except for the difference in the TOI onset time/location mentioned earlier. The amplitudes of modelled TEC anomalies (both SED and TOI) appear somewhat higher relative to the observations, confirming the assessment of Liu et al. (2016) that the TIE-GCM generally overestimates the magnitude of positive storm anomalies at high latitudes, though the specific reasons for this overestimation cannot be addressed here. The overestimation of TEC in TIE-GCM simulations, relative to IGS TEC maps, may reach up to 30–40 TEC units in the vicinity of the geomagnetic North Pole (Fig. 5). For the smaller geomagnetic storm of March 2015, TIE-GCM simulations overestimate the TOI anomaly by 10–15 TEC units relative to IGS TEC maps (see Fig. 4 in Liu et al., 2016, for comparison). IGS TEC maps are also expected to suffer from the lack of ground GNSS receivers in the Arctic Ocean sector. Thus, the comparison between the modelled TEC and the observed IGS TEC is questionable in the nightside region beyond the North Pole. The TIE-GCM grid singularity at the geographic North Pole may also lead to numerical problems in cross-polar plasma transport and continuity.

5.2 Plasma uplift dynamics

At first we analyse the dynamics of plasma uplift without looking into specific uplift mechanisms. As most of the ionospheric plasma is expected to be confined in the vicinity of the F2 peak, it is instructive to compare TEC distributions to the height and density of the F2 peak. The comparison (see Fig. 2 and the Supplements) confirms that the F2 peak plasma density (NmF2) largely mimics the behaviour of TEC. In contrast, the change in F2 peak height (ΔhmF2) shows a more complex behaviour. Substantial enhancements of hmF2 (up to 300 km) appear in the following longitude sectors (as referred to 15:00 UT 20 November 2003): (a) central part of the mainland USA west of 80^∘ W, westward of the main SED anomaly; (b) eastern coast of Canada and towards the geomagnetic North Pole at 45–65^∘ W, corresponding to the TOI location; and (c) eastern part of Europe and Central Asia east of 20^∘ E, in the post-sunset sector. Out of these three major hmF2 enhancements, only the TOI-related enhancement (b) is accompanied by a clear increase in plasma density and TEC, while the other two enhancements are accompanied by negative density anomalies. The post-sunset enhancement in hmF2 (c) is considered to be related to a sudden significant increase in hmF2 reported in Borries et al. (2017), which is accompanied by an extreme increase in the equivalent slab thickness. The authors consider intensive plasma transport with strong vertical components at this period of time over the respective region. The most westward enhancement in hmF2 (a) is due to the early formation of an SED anomaly in that sector and does not have a clear connection to the TOI anomaly. Some secondary positive/negative anomalies in hmF2 are seen in conjunction/alignment with the auroral TEC anomalies and will not be discussed here. The main focus here is the clear enhancement of hmF2 coinciding with the positive anomaly in NmF2 and TEC at the poleward edge of the SED anomaly and the throat of the TOI anomaly, lasting from ∼ 14:00 to 19:00 UT.

5.3 Electrodynamic vs. neutral wind transport

We first focus on the comparison between the modelled relative TEC distributions and the electrodynamic transport parameters shown in Fig. 3 and in the Supplements. As indicated by the electric potential distributions, the high-latitude plasma convection pattern greatly expands equatorwards and develops the characteristic two-cell pattern following the southward IMF turn at 11:00–12:00 UT. The expanded two-cell convection pattern persists through the storm's main phase, reaching the maximum expansion at 17:00–18:00 UT around the minimum of the SYMH index. This is consistent with the DMSP satellite observations of E×B convection during this storm (see Figs. 5 and 6 in Pokhotelov et al., 2008). The comparison shows that the Weimer model used in TIE-GCM underestimates the degree of equatorward expansion. Due to IMF By being strongly positive in the early main phase (11:00–15:00 UT), the convection “throat” is initially oriented NW–SE, later changing its orientation to NE–SW, when the IMF B_y becomes negative around 18:00 UT. This change in orientation of the convective channel is clearly reflected in the shape of the TOI TEC anomaly. The influence of east–west convection asymmetry on the TOI anomaly due to the IMF B_y dynamics has been reported before (e.g. Sojka et al., 1994), and it requires further analysis, which is outside the scope of this study. The important feature of electrodynamic plasma transport is the enhancement in the vertical electric drift component (W_E×B) seen at latitudes from 60^∘ N down to 40–45^∘ N, which accompanies the equatorward expansion of plasma convection. The vertical drift component arises from the E×B convection expanded to latitudes where dipolar magnetic field lines are far from vertical (e.g. Swisdak et al., 2006). The vertical electric drift maximises in the same longitudinal sector as the TOI anomaly. It maximises at the poleward edge of the SED anomaly and in the throat of the cross-polar convection channel (∼ 70^∘ W, 50–60^∘ N) but has a larger E–W extension (∼ 30–100^∘ W) than the TOI anomaly itself. Additionally, enhanced vertical drifts are seen at ∼ 40 ^∘ N in a broader range of longitudes extending into the central–western USA sector (west of 90^∘ W). The amplitudes of vertical drifts of ∼ 200 m/s appear to be very large, but they are generally consistent with occasional storm-time measurements of large vertical plasma drifts by the mid-latitude Millstone Hill incoherent scatter radar (Yeh and Foster, 1990; Erickson et al., 2010; Zhang et al., 2017) and the uplifts of the F2 peak by ∼ 400 km within 1 h estimated from tomographic reconstructions during the main phase of the 30 October 2003 superstorm (Yin et al., 2006).

During storms, Joule heating in the auroral region changes thermospheric winds and generates so-called storm wind cells (Volland, 1983). The model results show (Fig. 4) the enhanced Joule heating near the throat region at about 60^∘ N, but the amplitude is small compared to the Joule heating in the nightside auroral region ($\sim \mathrm{140}{}^{\circ }$ W$--\mathrm{120}{}^{\circ }$ E). The heating-induced equatorward neutral winds are expected to cause plasma uplift at subauroral latitudes, contributing to the formation of SED (Rishbeth, 1998; Swisdak et al., 2006) and possibly TOI anomalies. The modelled distributions of meridional neutral winds (see Fig. 4 and the Supplements) clearly show enhancements of winds (200–300 m/s at 120 km height and up to 500 m/s at 400 km height) in the longitudinal sector of the TOI anomaly blowing in the anti-sunward (cross-polar) direction, even partially equatorwards of the heating region (60–80^∘ W). Enhanced equatorward neutral winds are primarily seen in the central–western USA sector (west of 90^∘ W). We also notice that at the early stage of the TOI formation (∼ 13:00 UT) the meridional neutral winds are nearly zero at the poleward edge of SED and at the throat of the convective channel but become polewards later on and appear at higher latitudes. This is an indication that at the poleward edge of SED and in the throat region, forcing from the enhanced E×B convection flow is stronger than the forcing from heating-induced winds. The cross-polar neutral wind is mainly driven by the plasma convection, thus forming the polar cap neutral tongue anomaly (Burns et al., 2004).

5.4 Relations to other modelling efforts and space weather applications

After comparing the electrodynamics and neutral wind dynamics, we conclude that the uplift due to the vertical component of enhanced E×B convection is the dominant mechanism forming the TOI anomaly. This is generally consistent with the conclusions of Liu et al. (2016), based on TIE-GCM modelling of two relatively smaller storms (with Dst minima of −132 and −226 nT) driven with the Weimer convection model. Liu et al. (2016) concluded that around the F2 region peak and above, the electric field transport is the dominant driver of positive SED/TOI anomalies, while at lower altitudes (∼ 280 km) neutral winds could play a major role in producing the positive anomalies. The dominant role of electrodynamic uplift/transport is also confirmed by Huba et al. (2017), who used the SAMI3 model driven with the Rice Convection Model (RCM), showing that the realistic TOI anomaly can be reproduced even without including the neutral wind dynamo. The dominant role of electrodynamic plasma uplift in the formation of the TOI anomaly does not rule out a complex interplay between electric convection, neutral winds, and other possible mechanisms responsible for the formation of a mid-latitude SED anomaly (e.g. Swisdak et al., 2006; Crowley et al., 2006), which is outside the scope of this study. It is also possible that during relatively smaller geomagnetic events, such as the March 2013 and March 2015 storms (Liu et al., 2016; Klimenko et al., 2019), the effects of neutral winds and compositional changes are more pronounced compared to the superstorm case presented here. For instance, we do not observe such clear suppression of the TOI anomaly beyond the geomagnetic North Pole, as noticed by Klimenko et al. (2019) during the March 2015 storm. In contrast to Liu et al. (2016), showing that at lower altitudes (∼ 280 km) the neutral wind effects may cause the enhancement of TOI density (see Fig. 7 in Liu et al., 2016), we obtain strongly poleward neutral winds at all altitudes down to ∼ 120 km (see Fig. 4), implying negative or no contribution of the neutral wind effects to the TOI formation.

The conclusions above are subject to the right choice of high-latitude E×B plasma convection model. The Weimer parameterisation (Weimer, 2005) used here to drive the TIE-GCM simulations and also for the earlier tomographic reconstructions (Spencer and Mitchell, 2007; Pokhotelov et al., 2008) should provide a realistic response to the rapid changes in solar wind/IMF conditions, which could be missing in the case of Heelis parameterisation (Heelis et al., 1982) based on the 3 h resolution planetary index (Kp). Our TIE-GCM simulations repeated for the 20 November 2003 storm using the Heelis convection parameterisation (not shown here but available on request) demonstrated relatively poor agreement with IGS TEC maps and tomography. Pokhotelov et al. (2008) demonstrated that the statistical Weimer parameterisation may not be able to capture the true extent of equatorward expansion of the E×B convection pattern during the superstorm. The mismatch between the times/longitudes of the early TOI formation (the TOI anomaly appears earlier in time and more eastward in IGS TEC maps relative to the TIE-GCM simulations, as noted in Sect. 5.1) is likely due to this underestimation of the E×B expansion. Simulations driven with more realistic convection patterns obtained from e.g. radar network observations during a specific storm (Wu et al., 2015) or from assimilative models (Lu et al., 2016) may be needed to overcome these deficiencies.

While it is clear that the numerical setup of the CTIPe model (namely, the coarse resolution of 18^∘ in longitude) is not ideal for analysing the TOI anomaly, it is beneficial to discuss the results of this model in the context of space weather applications as the CTIPe is currently used for operational analysis and forecast by the US National Oceanic and Atmospheric Administration Space Weather Prediction Center – https://www.swpc.noaa.gov/models (last access: 21 September 2021) (Codrescu et al., 2012). The CTIPe simulation of the 20 November 2003 storm by Fernandez-Gomez et al. (2019) extended to the North American sector does not show clear TOI developments, though the CTIPe reproduces enhanced neutral wind patterns in the polar cap (Fig. A1) similar to those modelled by the TIE-GCM. On the other hand, Pryse et al. (2009) demonstrated that the CTIP model (Millward et al., 1996) was able to reproduce some features of the TOI anomaly consistent with ionospheric tomography when the simulation was driven by the SuperDARN radar observations of plasma convection. The use of SuperDARN data for driving the simulations was not addressed here but should be exploited in the future.

A fragmentation of the TOI anomaly due to IMF dynamics and other mechanisms has long been attributed to the formation of polar cap plasma patches (Sojka et al., 1994; Carlson et al., 2004). Climatological studies of ionospheric GNSS scintillations at high latitudes (e.g. Prikryl et al., 2015) demonstrate strong correlations with the plasma patches, especially near noon in the cusp region and near midnight, i.e. near the exit from a cross-polar convection channel². One has to note that polar patches are formally defined as drifting F-region plasma irregularities with horizontal scales ∼ 100 km and densities 2–10 times above the background and could also be formed during geomagnetically quiet times (Moen et al., 2013). Nevertheless, the TOI anomaly is expected to be a dominant source of the high-latitude GNSS disruptions during geomagnetic storms, and it needs to be addressed in space weather applications.

6 Summary and conclusions

The feeding mechanisms of the TOI anomaly have been analysed using the simulations of the geomagnetic superstorm of 20 November 2003, which have been conducted using the high-resolution version of the TIE-GCM ionospheric circulation model with the Weimer parameterisation of high-latitude E×B plasma convection. The simulation results are compared to the IGS TEC maps and to the results of ionospheric GNSS tomography for this storm event, published earlier by the authors (Pokhotelov et al., 2008). The main conclusions are summarised as the following.

• a.

  The TIE-GCM simulations reproduce the development of the polar TOI anomaly consistently with the IGS TEC maps and the tomographic TEC reconstructions. Differences between the model and observations are seen in the early formation of the TOI anomaly and in the magnitude/longitudinal extent of the TEC anomaly across the polar cap. The results of TIE-GCM simulations are qualitatively consistent with earlier modelling of less severe geomagnetic storms with the TIE-GCM and other ionospheric models (Liu et al., 2016; Huba et al., 2017). The TIE-GCM substantially overestimates the amplitude of the TOI anomaly in the polar cap (up to 40 %–50 % overestimation in TEC units) relative to IGS TEC maps, which is generally consistent with the TIE-GCM simulations of smaller storms (Liu et al., 2016). The large uplift velocities shown by the TIE-GCM near the poleward edge of the SED anomaly and in the convection throat agree with earlier ionospheric tomography results and with radar observations of vertical drifts during large storms. The noted differences between the modelled TEC and IGS TEC maps can be attributed to the model deficiencies (especially the E×B convection parameterisations during storms) and to poor GNSS data coverage in the polar cap. More rigorous data–model comparisons using more recent storms with better GNSS coverage are needed.

• b.

  Simulated distributions of the plasma and neutral dynamics demonstrate that the plasma uplifts of ∼ 200 m/s due to the high-latitude E×B plasma convection expanded to mid-latitudes appear to be the dominant mechanism responsible for the formation of the TOI anomaly. The neutral winds, enhanced during the storm, show the pattern which is not able to actively contribute to the TOI formation. This contrasts with the published simulations of relatively smaller geomagnetic storms (Liu et al., 2016; Klimenko et al., 2019), when neutral winds play a substantial role in forming and/or suppressing the TOI anomaly, especially at lower altitudes (Liu et al., 2016). We also show that the SED anomaly at mid-latitudes is likely to be influenced by both neutral wind and electrodynamic transport mechanisms, which is consistent with earlier simulations.

• c.

  Comparisons between the TIE-GCM and CTIPe model show that the lower-resolution CTIPe model, currently used for space weather operations, is not able to reproduce the TOI anomaly correctly. On the other hand, TIE-GCM simulation of the TOI anomaly also has clear deficiencies. Better model representation of the E×B plasma convection during extreme geomagnetic storms is needed.

Appendix A

Figure A1 Modelled distributions of relative ΔTEC (a), Joule heating (b), and relative meridional neutral winds at different levels (c)–(d) at 15:00 UT 20 November 2003 obtained from the CTIPe simulations.