2.1 WACCM-X

We use Whole Atmosphere Community Climate Model with thermosphere and ionosphere extension (WACCM-X) version 2.0 (Liu et al., 2018), a configuration of the NCAR Community Earth System Model (CESM 2.0; Hurrell et al., 2013) to perform model simulations in this study. WACCM-X extends from the surface to the upper thermosphere with its top boundary, depending on the solar and geomagnetic activity, lying between 500 and 700 km. The vertical resolution of WACCM-X above the stratosphere is one fourth of a scale height, and the horizontal resolution is 1. x 2. in latitude and longitude, respectively. WACCM-X is built upon the Whole Atmosphere Community Climate Model (WACCM) (Marsh et al., 2013) and Community Atmosphere Model version 4 (Neale et al., 2013). The new version of WACCM-X has a coupled ionosphere and incorporates self-consistent low-mid-latitude ionospheric electrodynamics adapted from the Thermosphere-Ionosphere Electrodynamics General Circulation Model (TIE-GCM). At high-latitudes, WACCM-X uses an empirical electric potential pattern (Heelis et al., 1982), which is parameterized by the 3-h geomagnetic Kp index and an auroral precipitation oval based on the formulation described by Roble and Ridley (1987). More details about the physical processes included in WACCM-X 2.0 can be found in Liu et al. (2018).

WACCM-X provides a comprehensive tool to study the entire atmosphere-ionosphere system. The impact of lower atmospheric forcing on the upper atmospheric variability can be studied using WACCM-X during specific time periods by constraining the tropospheric and stratospheric dynamics to meteorological reanalysis fields. In the present study, we use the specified dynamics (SD) set up in WACCM-X to simulate the SSWs by constraining the winds and temperatures from 0 to 50 km toward the National Aeronautics and Space Administration Modern Era Retrospective Analysis for Research and Applications (MERRA) Version 2 (Gelaro et al., 2017). Using the approach described in Kunz et al. (2011), the WACCM-X model fields are constrained to the MERRA2 meteorological fields at every model time step (i.e., 5 min). The default SD-WACCM-X set up does not include forcing from the M2 tide but we implement it in our simulations based on the method described by Pedatella, Liu, and Richmond (2012) because the M2 forcing becomes an important source of MLT and ionospheric variability during SSWs. Hourly outputs of winds, temperature and geopotential height are obtained for the 2009 and 2019 SSWs through the SD-WACCM-X runs in this study.

2.2 TIE-GCM

TIE-GCM is a three dimensional, self-consistent numerical model of the coupled thermosphere-ionosphere system that has been developed at the High Altitude Observatory at the National Center for Atmospheric Research (NCAR). The model spans from 97 km to about 500–700 km depending on the solar cycle activity. In this study, we use TIE-GCM version 2.0 with a horizontal resolution of 2. by 2. in geographic longitude and latitude and a vertical resolution of 0.25 times the scale height. The input parameters for TIE-GCM include the solar XUV, EUV and FUV spectral fluxes that are defined by the EUVAC model (Richards et al., 1994) using the F10.7 index. The global electric potential due to the wind dynamo is solved by the self-consistent TIE-GCM ionospheric electrodynamo at low- and mid-latitudes. At high-latitudes, however, the electric potential is prescribed through empirical convection electric field patterns using the Heelis (Heelis et al., 1982) or Weimer (Weimer, 2005) models. TIE-GCM also uses an analytical auroral model to account for high-latitude auroral particle precipitation in the default set up. The upper part of WACCM-X is based on TIE-GCM and it uses this same auroral model.

The effect of lower atmospheric tidal forcing can be specified in TIE-GCM using the tidal perturbations at its lower boundary (LB). The amplitudes and phases of upward propagating atmospheric tides specified at the TIE-GCM LB in the default setup are based on the Global Scale Wave Model (GSWM) (Hagan et al., 1999). The default TIE-GCM LB assumes constant neutral temperature (T = 181 K), geopotential height (Z = 96.37 km) and zero horizontal winds (e.g., Maute, 2017). For a more realistic LB conditions, TIE-GCM includes the option to specify hourly inputs at its LB from any other source. In this study, we specify the hourly WACCM-X outputs for 2009 and 2019 SSWs at the TIE-GCM LB and carry out three simulations for each of the SSWs in order to isolate the effects of geomagnetic and lower atmospheric forcings on the TEC variability. In the first simulation setup (hereafter referred to as S1), the TIE-GCM forced by WACCM-X is run in its default mode and the obtained day-to-day ionospheric variability from this run includes the effects of both geomagnetic and lower atmospheric forcings. In the second simulation setup (hereafter referred to as S2), we turn-off the geomagnetic forcing and carry out a similar run for both SSWs. The day-to-day ionospheric variability resulting from the second run arises solely from lower atmospheric forcing during both SSWs. The geomagnetic forcing is turned off from the TIE-GCM runs by reducing the hemispheric power from 18 GW in the first setup to 0.1 GW and the cross-polar cap potential from 30 kV in the first setup to 0.1 kV. Additionally, the Heelis convection model and the analytical auroral model have also been turned off in the second setup to remove the magnetospheric energy input. Using S1 and S2, the contribution of geomagnetic forcing to the TEC variability during both SSWs is estimated. We also perform a third TIE-GCM simulation (hereafter referred to as S3) to assess the importance of SW2 and M2 tides to the TEC variability during both SSWs. In S3, we utilize the settings of S2 and additionally remove the SW2 and M2 tides from the WACCM-X input to TIE-GCM LB. Through comparison of S2 and S3, the contribution of these semidiurnal tides to the TEC changes during both SSWs is estimated. The experiment setups used in this study have been summarized in Table 1. We used TIE-GCM instead of WACCM-X to carry out these three simulations because of the former being computationally less expensive than the latter.

The use of simulations to study the ionospheric effects solely due to lower atmospheric forcing have been implemented in some earlier studies. Pedatella and Maute (2015) simulated the ionospheric effects during the 2013 SSW event solely due to lower atmospheric forcing by running the TIME-GCM simulations under constant solar and geomagnetic activity levels. Yamazaki et al. (2014) used a similar method by running the TIE-GCM under constant solar and magnetospheric energy inputs to study the day-to-day variability of the equatorial electrojet due to lower atmospheric forcing. The occurrence of an SSW event during quiet-time ionospheric conditions is a rarity and there have been many SSWs that occur under active geomagnetic conditions, which complicates the separation of SSW driven ionospheric variability. Our method builds upon the techniques that are presented in Pedatella and Maute (2015) and Yamazaki et al. (2014) and provides a tool to separate the geomagnetic and lower atmospheric forcing effects on the ionospheric variability.

4.1 Zonal Mean and PW Variability During 2009 and 2019 SSWs

In Figure 1, the zonal mean temperatures averaged between and N during January–February 2009 (Figures 1a and 1c) and December 2018–January 2019 (Figures 1b and 1d) from Aura MLS observations (top panels) and SD-WACCM-X simulations (bottom panels) are presented. The vertical dashed black lines in the figure mark the day of polar vortex weakening (PVW). As an alternative to the classical definition of SSW provided by the WMO, PVW has been used to correlate the tidal enhancements in the MLT and ionosphere with the magnitude of the reversal of stratospheric ZMZW (e.g., Chau et al., 2015; Siddiqui et al., 2017; Zhang & Forbes, 2014). A PVW day is identified by locating the earliest and most extreme reversal of ZMZW at N and 48 km altitude (1 hPa) that occurs simultaneously with the increase in zonal mean temperature at North Pole and 40 km altitude (3 hPa) between December and February.

It can be seen that there is a reasonable agreement between the observed and simulated zonal mean temperatures during the considered time intervals. With the onset of each of the SSWs, the warm stratopause descends from its climatological position near 60 km (0.2 hPa) toward lower altitudes, which results in warming at these heights (e.g., Labitzke, 1981) and later breaks down completely in the peak stages of the SSW. The stratopause then reappears at a higher altitude before slowly returning to its original position. This altitudinal shift and reemergence of the stratopause is called an elevated stratopause event and is often associated with major SSW events (e.g., Manney et al., 2008; Siskind et al., 2007). In WACCM-X simulations, the elevated stratopause events are being reproduced for both SSWs but their appearances tend to occur slightly earlier in comparison with the Aura MLS observations. The PVW date for the 2009 SSW event occurs on day 23 (Jan 23) and for the 2019 SSW event on day −3 (Dec 28). The warming of the stratosphere is accompanied by cooling in the mesosphere for both these events. In case of the 2009 SSW, mesospheric cooling can be observed to start around day 20 above 0.1 hPa in Figure 1a and in case of the 2019 SSW around day −8 in Figure 1b. The model temperature simulations are able to reproduce these features to a large extent and are overall consistent with the Aura MLS observations but some discrepancies can also be seen especially at altitudes where SD-WACCM-X stops being constrained to MERRA2, that is, above 50 km (1 hPa). The differences between SD-WACCM-X and Aura MLS observations start to become apparent above this altitude and are more pronounced in the MLT region. These temperature differences could be related to gravity wave forcing, which is parameterized in SD-WACCM-X and may be contributing to the discrepancy between the observations and model simulations (e.g., Smith, 2012).

The left panels of Figure 2 present the ZMZW at N from SD-WACCM-X for the time intervals that include the 2009 and 2019 SSWs while the right panels show the Kp index and daily solar flux (s.f.u) conditions during these SSWs. The dashed black lines show the PVW days in all the panels. For the 2009 SSW, the ZMZW between 1 and 10 hPa in the stratosphere changes from eastward to westward direction starting around day 20 in Figure 2a. At 10 hPa, the ZMZW remains in the reversed direction until day 54 while at 1 hPa, where the PVW is defined, the ZMZW reaches a peak reversal on day 23 with a value of about −52 m/s. For the 2019 SSW, the reversal of the ZMZW between 1 and 10 hPa is seen toward the end of December around day −8 in Figure 2c. The westward direction of ZMZW at 10 hPa remains till day 21 while at 1 hPa, the ZMZW is found to reach a peak reversal on day −3 with a value of −15 m/s. From the comparison of the ZMZW in Figures 2a and 2c, it can be clearly seen that the 2009 SSW event was stronger in terms of ZMZW reversal and more prolonged than the 2019 SSW event.

From the Kp indices and solar flux values in Figures 2b and 2d, it can be inferred that both the 2009 and 2019 SSWs were recorded under periods of relatively low solar and geomagnetic activities. The solar flux levels during 2009 and 2019 SSWs remained below 75 s. f.u for both the events. The geomagnetic activity during the onset and peak phase of the 2009 SSW hovered mostly around Kp with an exception on day 19 when Kp values reached during the 3-hourly intervals. The 2019 SSW event showed higher geomagnetic activity levels as compared to the 2009 SSW event with brief periods of spike in Kp values that reached up to on day −3 (Dec 28) and on day 5 but overall the period during the 2019 SSW remained geomagnetically quiet.

Figure 3 presents the variability of PWs with wave number 1 (PW1) and 2 (PW2) at 10 hPa in temperature from SD-WACCM-X for the 2009 and 2019 SSWs. In this figure, the colorbar scales of PW1 are chosen to be twice as large in magnitude as compared to those for PW2. The dashed black lines show the PVW days in all the panels. In Figures 3a and 3c, the amplitudes of PW1 and PW2, respectively, are presented for the 2009 SSW event. The enhancement of PW1 is seen in NH high-latitudes particularly around N with peaks on days 6 and 23 while in the case of PW2, the enhancement begins in the second week of January and peaks on day 19. Based on the enhanced amplitudes of PW2 in Figure 3c, we find that our results are consistent with earlier studies (e.g., Manney et al., 2009) that have classified the 2009 SSW as PW2 forced split SSW event. The amplitudes of PW1 and PW2 for the 2009 SSW are also similar to the results shown using four different whole atmosphere models by Pedatella et al. (2014).

The amplitude of PW1 and PW2 in temperature are presented in Figures 3b and 3d, respectively, for the 2019 SSW event. The PW1 amplitudes show enhancement poleward of N during December with a maximum on day −7. The peak PW1 amplitudes for this event are almost twice as large as compared to that of peak PW1 for the 2009 SSW. The PW2 amplitudes are smaller in magnitude as compared to that of PW1 for this event but enhancement in PW2 is also seen poleward of N toward the end of December with a maxima centered on day −9. The 2019 SSW has been classified in literature as neither a typical displaced nor a typical split SSW event (Rao et al., 2019) but rather a mixed type event, which was initially a displaced SSW and later became a split SSW. From Figure 3b, it can be inferred by the dominance of PW1 amplitudes that the 2019 SSW must have started as a displaced SSW event.

4.2 Migrating Semidiurnal Tides During 2009 and 2019 SSWs

It is well established now that the primary reason for the variability in the ionosphere during SSWs is due to the modulation of atmospheric tides. In particular, the variability of SW2 and M2 tides based on modeling and observational studies have been found to be the most significant driver of ionospheric variabilities (e.g., Chau et al., 2009; Fejer et al., 2010; Forbes & Zhang, 2012; Goncharenko et al., 2012; Lin et al., 2013; Pedatella & Forbes, 2010; Pedatella et al., 2014; Siddiqui et al., 2018; Vineeth et al., 2009; Yamazaki et al., 2012). Atmospheric tides refer to planetary-scale oscillations of the atmosphere that are mainly excited by gravitational forces of the moon and by thermal forcing from the sun (e.g., Forbes & Garrett, 1978; Forbes, 1982; Lindzen & Chapman, 1969). These oscillations have periods and sub-periods of a solar or a lunar day. The solar tides form the dominating component of the tidal oscillations and are predominantly thermally forced. The excitation mechanisms include daily periodic absorption of solar energy by tropospheric water vapor and stratospheric ozone (e.g., Forbes & Wu, 2006; Zhang et al., 2010) while the relatively smaller lunar tides are mainly forced due to the lunar gravitational effects on the Earth's atmosphere. In this section, we investigate the variability of the SW2 and M2 tides during both SSW events.

The hourly outputs of neutral temperature from SD-WACCM-X simulations are used to extract the components of the solar and lunar tides by performing a least squares fit with a moving window of 15 days at each latitude. The following equation based on Pedatella, Liu, and Richmond (2012) has been used to make the fit:

(1)

where , is universal time in hours, denotes the harmonics of a solar day, is the longitude, is the zonal wave number, and are the amplitude and phase of the solar tidal components. For a wave propagating westward , while for a wave propagating eastward . and represent the amplitude and the phase of the semidiurnal lunar tide, respectively.

Figure 4 presents the amplitudes of the SW2 (top panels) and M2 (bottom panels) tides in neutral temperature at  hPa (110 km altitude) from SD-WACCM-X simulations for the 2009 and 2019 SSWs. The vertical dashed black lines mark the days of PVW for each of these SSWs. In Figure 4a, the SW2 tidal variability for the 2009 SSW shows enhanced amplitudes of 21 K at SH mid-latitudes on day 15 and amplitudes 25 K starting around day 35. In the NH, the SW2 enhancements are only seen after day 25 and enhanced amplitudes of 24 K are found from simulations. Following the PVW day, there is a sudden decrease in the amplitude of SW2 at SH mid-latitudes. This feature of SW2 was also reported in a study by Pedatella et al. (2014) where they compared the temporal variability of SW2 during the 2009 SSW event using four different whole atmosphere models. Figure 4c presents the amplitude of M2 tide where we notice its enhancement to 10 K in both hemispheres a few days after the PVW day. Another enhancement of M2 reaching 9 K is seen in the NH after day 40. The M2 variability from SD-WACCM-X during this SSW event is consistent with the results of Zhang and Forbes (2014). They also reported similar M2 enhancements using neutral temperature measurements at 110 km altitude from the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument onboard the Thermosphere Ionosphere Mesosphere Energetics Dynamics (TIMED) satellite. Figures 4b and 4d shows the amplitudes of SW2 and M2 tides, respectively, in neutral temperature between December 2018 and January 2019. Similar to the SW2 variability during the 2009 SSW event, the SW2 amplitude in Figure 4b shows enhancement at SH mid-latitudes on either side of the PVW day with a reduction of SW2 amplitude in between. The SW2 amplitude reaches a peak of 17 K on day −12 and 20 K on day 8 in the SH. In the NH, the SW2 amplification starts close to the PVW day and it reaches a value of 18 K on day 9. In Figure 4d, the first M2 enhancement in SH with peak amplitudes of 6.5 K happens a few days after the PVW day. The second M2 enhancement in the SH with peak amplitude of 7.5 K is seen around day 11. Following Chau et al. (2015) and Conte et al. (2017), we also performed tidal analyses using a 21-day running window (see Supporting Information S1) to further reduce any artifacts or ambiguity between the determination of SW2 and M2 tides but we found the amplification of M2 to be similar to that shown in Figures 4c and 4d. We find that the M2 amplitude during the 2019 SSW is smaller (2.5 K) as compared to that during the 2009 SSW event. One reason for the difference in amplification could be related to the timing of the SSW event relative to the phase of the moon. Pedatella and Liu (2013) have shown from simulation results that the lunar tidal response in the ionosphere is dependent upon the phase of the moon relative to the timing of the SSW event. Further, Fejer et al. (2010) have found that the lunar effects in the ionosphere during SSWs amplify close to the new or the full moon days. The new moon day for the 2009 SSW occurred 3 days after the PVW day while for the 2019 it was recorded 6 days before the PVW. This difference in the timing of the new moon days relative to the PVW days may also be contributing to the reduced amplification of M2 during the 2019 SSW.

4.3 Comparison of Observed and Simulated M2 Variability During 2009 and 2019 SSWs

Figure 5 presents the M2 tidal amplitudes for the 2009 and 2019 SSW events obtained using the V2.0 temperature measurements from the SABER instrument onboard the TIMED satellite at an altitude of 110 km. We employ the least squares fitting method mentioned in Zhang and Forbes (2014) to determine the M2 amplitudes from zonally averaged SABER temperature residuals. As the period of the M2 tide from the frame of the TIMED satellite is 11.85 days (Forbes et al., 2013), we use a 12-day moving window and fit only the M2 tide to the daily zonally averaged temperature residuals. Figures 5a and 5b show the amplitude of the M2 tide for the 2009 and 2019 SSWs, respectively. The vertical dashed black lines show the days of PVW. The enhancement of M2 following the PVW can be seen at low- and mid-latitudes for the 2009 SSW event in Figure 5a. This plot of M2 amplitude is similar to the one shown in Zhang and Forbes (2014) (see Figure 1). The M2 tides in neutral temperature from SD-WACCM-X simulations for the 2009 SSW, shown in Figure 4c, match very well with the M2 from SABER temperature observations in terms of the timing of the M2 enhancement. In the NH, the day-to-day variability of the M2 amplitudes from SD-WACCM-X simulations is slightly more consistent with those obtained from SABER temperature observations as the M2 peaks are located at similar latitudes. In the SH, the M2 peaks from SABER observations on day 27 are slightly more equatorward as compared to those from SD-WACCM-X simulations. The M2 amplitudes for the 2019 SSW presented in Figure 5b shows that the level of M2 enhancement for this SSW was clearly lower than that for the 2009 SSW. The M2 enhancements can be seen to occur toward the end of December 2018 between days −7 and 0 poleward of S and around day 22 above N. The enhancement in the SH toward the end of December is also captured in M2 from SD-WACCM-X simulations in Figure 4d. But the M2 amplitudes from SD-WACCM-X do not exactly reproduce the observations for the 2019 SSW event as the M2 peak around day 22 in the NH is not seen in the simulations. There is similarity between the observed and simulated M2 amplitudes in the sense that the enhancements were particularly weaker for the 2019 SSW as compared to the 2009 SSW event, which is confirmed from both the simulated and observed M2 amplitudes in the simulations. There is similarity between the observed and simulated M2 amplitudes in the sense that the enhancements were particularly weaker for the 2019 SSW as compared to the 2009 SSW event, which is confirmed from both the simulated and observed M2 amplitudes.

4.4 GNSS TEC Variability Over Europe During 2009 and 2019 SSWs

In this section of the paper, we analyze the variability of GNSS based TEC over Europe during the 2009 and 2019 SSWs. We compute the hourly TEC data between 35°N and N geographic latitudes and W-E, in the European sector by taking the mean vertical TEC values within 30 min for each local time (LT). The standard deviation for each hourly bin lies typically between 0 and 1.5 TECu. In order to examine the SSW associated TEC perturbations, we first determine the hourly quiet background TEC using median values from a 27-day moving window that is centered on a given day. The 27-day median TEC () has been used as an effective method to estimate the quiet background TEC for comparison against disturbances (e.g., Borries et al., 2015; Mendillo, 2006). The perturbations in TEC, TEC, are then calculated during the 2009 and 2019 SSW events by using,

(2)

The typical uncertainties in the 27-day median values, estimated using the median of absolute deviations (MAD) (Müller, 2000) comes out to be between 0.5 and 1 TECu.

We looked at the TEC variability over Europe for all the LTs during the 2009 and 2019 SSWs and found that the maximum variability is observed at around 12 LT. In Figures 6 and 7, we therefore present TEC (in TECu) observed at 12 LT over Europe for 2009 and 2019 SSWs. For both these events, we cover the pre-, peak-, and post-SSW periods. TEC is presented on alternate days between January 15, 2009 and February 12, 2009 for the 2009 SSW. For the 2019 SSW, TEC is similarly presented between December 18, 2018 and January 15, 2019. In Figure 6, TEC for the 2009 SSW event shows minor changes at European latitudes with magnitudes remaining less than 4 TECu over the entire period. Compared to TEC at low-latitudes, where magnitudes of more than 10 TECu were observed in the American sector at similar LTs during the 2009 SSW (e.g., Goncharenko et al., 2010), the TEC values at mid-latitudes are relatively minor in absolute terms. In Figure 7, TEC for the 2019 SSW is presented. In the pre-SSW period, TEC values are 3 TECu but on December 28, 2018, during the peak SSW phase, a relatively large enhancement is observed in TEC, which covers the majority of European latitudes. In the South-West parts of Europe, TEC values reach 10 TECu on this day. A reduction in TEC values by 4 TECu is seen on January 3, 2019 between 35°N and N that is followed by an enhancement of a similar amount between 40°N and N on January 5, 2019. The TEC magnitudes during the 2019 SSW event in Figure 7 are clearly stronger than for the 2009 SSW event in Figure 6.

It is to be noted that the timing of TEC enhancement on December 28, 2018 for the 2019 SSW event coincides with its PVW day. Moreover, from Figure 2d we have seen that the daytime Kp values increased to 40 on this day. The enhancement in TEC could therefore be a result of processes involving both geomagnetic and lower atmospheric forcings. The individual contribution of both these mechanisms to the TEC enhancements needs to be ascertained in order to quantify their relative importance for the respective SSW events. We investigate the TEC enhancement in more detail in the next section with TIE-GCM simulations to find the dominant forcing mechanism that caused the observed TEC variability over Europe.

5.1 Comparison Between Observed and Simulated TEC During 2009 and 2019 SSWs

Based on the GNSS TEC observations in Figure 7, it seems that the TEC variability over Europe during the 2019 SSW may also be affected by moderately enhanced geomagnetic activity levels as noted by the increase in Kp values. To identify the main driver for this TEC variability, we isolate the individual contributions of geomagnetic and lower atmospheric forcings to TEC variability using TIE-GCM simulations. Additionally, in order to check the veracity of the simulated TEC from TIE-GCM simulations during SSWs, we performed a comparison with the observed low-latitude TEC variations that was reported by Goncharenko et al. (2010) over the American sector during the 2009 SSW (see Supporting Information S2). We see the characteristic semi-diurnal signature in simulated TEC around the central day of 2009 SSW over the American sector that was reported by Goncharenko et al. (2010). As the TIE-GCM simulations are able to qualitatively reproduce the TEC observations at low-latitudes during the 2009 SSW, it gives us the confidence that they can also be used to study the mid-latitude TEC variability during SSWs.

In this section, we first present the observed GNSS TEC results in Figures 8a and 9a before comparing it with the simulated TEC from three different TIE-GCM setups. The readers may note that there has been a change in the sequential order with the 2019 SSW event being presented prior to the 2009 SSW event in order to facilitate the topical discussion. In Figure 8a, we take the averaged value of GNSS TEC over the European region shown in Figure 7 for each LT and present the temporal evolution of averaged TEC (Figure 8a) and TEC (Figure 8e) as a function of local time for the 2019 SSW. The dashed black lines mark the PVW day. The diurnal variation of TEC is evident in Figure 8a with the TEC values increasing gradually in the morning hours before reaching a maximum in the afternoon hours and then declining gradually post sunset. The day-to-day variability of TEC, which is subject to solar, geomagnetic and lower atmospheric forcings, can also be seen in Figure 8a. The TEC enhancements reach values of 12.5 TECu at 12 LT on day −3 and 8.5 TECu at 12 LT on day 5. TEC in Figure 8e is calculated using Equation 2 as before and shows the perturbations of TEC from the quiet background TEC values. The spike in both averaged TEC and TEC values on days −3 (Dec 28) and 5 (Jan 5) between 10 and 12 LT can also be seen in this figure, with TEC reaching 6 TECu on day −3 and to 2.5 TECu on day 5 at 12 LT.

Figure 8b shows the average of the simulated TEC over Europe for the S1 run. We note that the TEC derived from TIE-GCM simulations are able to reproduce the TEC spikes on days −3 (Dec 28) and 5 (Jan 5) and are qualitatively similar in comparison to the averaged GNSS TEC in Figure 8a. The modeled and observed TEC may only be compared in a qualitative sense owing to the upper boundary limit of TIE-GCM, which depending on the solar activity extends only up to 500–700 km in altitude, thus hindering a quantitative comparison with GNSS TEC observations. The modeled TEC results demonstrate that the TIE-GCM simulation includes the various forcing mechanisms that are responsible for the observed TEC variability in Figure 8a and can be used to filter out contributions of the mid-latitude sources of TEC variability. The averaged TEC enhancements over Europe reach values of 8.5 TECu at 12 LT on day −3 and 6.5 TECu at 13 LT on day 5. From the simulated TEC in Figure 8f, we can see that the TEC spikes on day −3 at 12 LT reach 4 TECu and on day 5 to 2 TECu at 12 LT.

In Figure 8c, the averaged TEC over Europe is presented for the S2 run. It can be noticed that there is an apparent reduction of the averaged TEC values after turning off geomagnetic forcing and the peak TEC values have been reduced to less than 5 TECu. The major TEC enhancements in this plot can be seen on days −3 and 4 with peak TEC values reaching 5 and 4.5 TECu, respectively, between 10 and 14 LT. In Figure 8d, the averaged TEC over Europe is presented for the S3 run. For this setup, the averaged TEC over Europe shows further reduction in the TEC peaks and on days −3 and 4 it reduces to the level of background TEC values between 10 and 14 LT, which suggests that the contribution of SW2 and M2 tides to the TEC variability for this SSW is not negligible and therefore needs to be assessed in more detail.

In Figure 8g, the difference of TEC values from S1 and S2 setups ( = ), is presented. The filled contour lines are plotted when absolute value of exceeds 0.25 TECu. It can be clearly noted that the major difference is seen on days −3 and 5 when peak values between 12 and 14 LT reach 3 and 2 TECu, respectively. Compared to the simulated background TEC values of 4 TECu at 12 LT, we find that the geomagnetic forcing is responsible for increasing the background TEC values by 75% on day −3 (Dec 28) and by 50% on day 5 (Jan 5) at 12 LT. It can therefore be easily inferred that the major processes causing the TEC spikes on days −3 and 5 between 10 and 14 LT are related to geomagnetic forcing. In Figure 8h, the difference of TEC values from S2 and S3 setups ( = ), is presented. Through this plot, the contribution of the SW2 and M2 tides to TEC variability can be assessed. The SW2 and M2 tides have been used in the analysis as they have been reported to be a major source of ionospheric variability during an SSW event (e.g., Fejer et al., 2010; Goncharenko et al., 2010; Jin et al., 2012; Pedatella et al., 2012). The filled contour lines are again plotted when absolute value of exceeds 0.25 TECu. We find that the TEC variability due to SW2 and M2 tides start close to the SSW onset on day −10 and continues up to day 12. The SW2 and M2 tidal changes associated with the 2019 SSW seem to modulate the daytime electric field through the E-region dynamo process, which results in the TEC variability as seen in this figure. The SW2 and M2 tides increase the background TEC values by approximately 1–1.5 TECu, which corresponds to 20%–25% increase in the background TEC.

The individual contributions of the geomagnetic and lower atmospheric forcings on the TEC variability over Europe is also similarly assessed for the 2009 SSW event. In Figure 9a, the averaged GNSS TEC over Europe as a function of local time is presented for the 2009 SSW event. In addition to the diurnal and day-to-day variability of TEC, major TEC enhancements are observed in Figure 9a between 10 and 12 LT on days 40 and 45. In Figure 9e, TEC shows the perturbations of TEC from the quiet median TEC values. Here, large TEC perturbations particularly on days 40 and 45 are clearly noticeable with TEC being 2 TECu between 10 and 13 LT. The averaged TEC derived from the TIE-GCM run with S1 setup is presented in Figure 9b. A qualitative comparison with GNSS TEC observations in Figure 9a suggests that the primary features of the averaged TEC variability are similar in the simulations. The comparatively lower levels of averaged TEC before day 30 between 10 and 15 LT and the moderately enhanced averaged TEC levels after this day have been correctly reproduced in the simulation. The simulated TEC enhancements on days 40 and 45 at 13 LT are smaller in comparison to observed TEC. We present the TEC from TIE-GCM simulations in Figure 9f. Compared to the GNSS TEC, we find that the TEC enhancements on days 40 and 45 between 10 and 14 LT are smaller in the simulated TEC with values reaching 1 TECu.

In Figures 9c and 9d, the averaged TEC over Europe from TIE-GCM runs with S2 and S3 setups are presented. We notice that there is a decrease in TEC values when the geomagnetic forcing has been turned off but we also find that most of the features of the averaged TEC from S1 setup have been reproduced in the averaged TEC from S2 setup, which suggests that the lower atmospheric forcing processes are the dominant drivers of TEC variability during the 2009 SSW event. We calculate the difference of TEC values from S1 and S2 setups () in Figure 9g and plot the filled contour lines when absolute values of exceed 0.25 TECu. We find that the TEC variability solely due to geomagnetic forcing lies mostly below 0.25 TECu in this period and exceeds this limit noticeably only on day 45 between 12 and 14 LT. In Figure 9h, the difference of TEC values from S2 and S3 setups, , is presented. As seen earlier for the case of 2019 SSW, the contribution of SW2 and M2 tides to TEC variability can be assessed through this plot. The filled contour lines are plotted when exceeds 0.25 TECu. The TEC variability due to SW2 and M2 tides starts around day 20 when the onset of the 2009 SSW event occurs and continues onward till day 50. The TEC variability reaches 1 TECu at 10 LT between days 23 and 28 and at 11 LT around day 40 before dropping to 0.75 TECu at 13 LT around day 45. For the 2009 SSW event, the SW2 and M2 tides are responsible for 20%–25% increase in the background TEC, which is similar to the levels seen during the 2019 SSW event.

5.2 Possible Reasons for the Observed TEC Variability

Most of the studies that documented the ionospheric effects of SSWs, especially during the 2009 SSW, focused on the variability at equatorial and low-latitudes (e.g., Chau et al., 2012; Yiğit & Medvedev, 2015, and references therein). While it is now accepted that the mechanisms causing the variability at these latitudes are driven by the changes in the vertically propagating semidiurnal solar and lunar tides, the mechanisms responsible for the mid-latitude ionospheric variability during SSWs are not as well understood. Simulation results by Pedatella and Maute (2015) have shown that the variability of the mid-latitude ionospheric F-region peak height (hmF2) during SSWs is predominantly driven by the field-aligned neutral winds, which is modulated by the M2 tidal enhancements. Yue et al. (2010) observed the global ionospheric response using Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) satellites during the 2009 SSW and suggested the changes in the neutral wind and composition due to direct propagation of tides as another mechanism for the ionospheric mid-latitude variability during SSWs. It is also known from observations and modeling studies that the influence of SSWs at mid- and high-latitudes ionosphere is generally smaller as compared to that at low-latitudes. Oyama et al. (2014) used the FORMOSAT-3/COSMIC peak ionospheric electron density (NmF2) data and found that changes in mid-latitude NmF2 to be only between 20% and 30% during the 2009 SSW, which was comparably lower than the changes at low-latitudes. Their observation results were found to be consistent with the simulations shown by Pedatella and Maute (2015). For the 2009 and 2019 SSW events, our results also show similar numbers as the background TEC values over Europe increase by 20%–25% due to SSW associated SW2 and M2 tidal changes. The amplification of SW2 and M2 tides in the peak SSW phase were relatively larger for the 2009 SSW in comparison with the 2019 SSW but their contributions to the TEC variability, in terms of TEC enhancement, remain similar. However, the SW2+M2 contribution to TEC variability lasts for a week longer during the 2009 SSW event as compared to the 2019 SSW event. For the 2019 SSW, the SW2+M2 influence on the TEC variability are seen between days −10 and 10 (3 weeks) and for the 2009 SSW between days 20 and 50 (4 weeks). This could be a result of a second M2 enhancement seen around day 45 for the 2009 SSW event, which contributes to a longer influence of SW2+M2 on TEC for this SSW.

During the 2019 SSW event, along with the increase in Kp values to 4 on December 28, the meridional component of the interplanetary magnetic field (IMF) in Geocentric Solar Magnetospheric (GSM) coordinate system, turned southward and reached down to −7.5 nT, and the Auroral Electrojet (AE) index reached to 1,100 nT. The symmetric disturbance field in H (SYM-H) index declined from 26 nT on December 27 to −30 nT on December 28. Based on the statistical measures, the geomagnetic activity parameters resemble the conditions of a weak geomagnetic storm (e.g., Yokoyama & Kamide, 1997). From the results in Figure 8, it is clear that the sudden surge in daytime TEC observations over Europe on December 28 is due to geomagnetic forcing. To explore the mechanisms behind the daytime TEC enhancement on this day, we present hourly GNSS TEC maps over Europe between 10 and 14 LT constructed with a resolution of 30 min in the top panels of Figure 10. The TEC values are seen to increase across Europe after 10 LT and reach a maxima at 12 LT before declining at later hours. An equatorward movement in the TEC can be observed between 12 and 14 LT, which suggests that these TEC fluctuations could be due to large-scale traveling ionospheric disturbances (LSTIDs). LSTIDs are ionospheric fluctuations associated with traveling atmospheric disturbances (TADs), which are generated in the auroral region when the polar ionosphere is subjected to energy input from the magnetosphere (e.g., Prólss, 1995). As the TADs propagate to middle- and low-latitudes, they induce perturbations in the neutral wind, which can lead to the upward/downward movement of plasma along the magnetic field lines leading to ionospheric fluctuations commonly known as LSTIDs (e.g., Lei et al., 2008). The propagation direction of LSTIDs is usually in the equatorward direction and they have a horizontal scale of 1,000 km and a period between 30 min and 3 h (e.g., Hunsucker, 1982). GNSS TEC maps are an excellent tool to identify the structure and temporal evolution of LSTIDs (e.g., Tsugawa et al., 2004). LSTIDs can be detected by detrending the observed TEC data with a 1 h running mean to retrieve the perturbation component of TEC (TECP) (e.g., Tsugawa et al., 2004). Below the top panels in Figure 10, we present sequence of TECP maps at an interval of every 20 min starting from 10 LT. It is possible to recognize the LSTID signatures in the TECP maps after 11 LT with the crest of LSTID starting around N and E and propagating in the equatorward direction. This crest continues to progress equatorward at later times and when it reaches N at 12:20 LT, a second crest appears at N. The LSTID signatures moved below N after 13 LT. The LSTID features seen in Figure 10 is consistent with the earlier observations of LSTIDs made over Europe (e.g., Zakharenkova et al., 2016). Therefore, we believe that major contribution to the TEC enhancement on December 28, 2018 could be caused by geomagnetic storm associated LSTIDs.

Along with LSTIDs, the increase in TEC over Europe could also be contributed by positive storm effect mechanism. During geomagnetic storms, the relative increase in the ionospheric plasma density with respect to quiet-time conditions is referred to as positive storm effect (e.g., Astafyeva et al., 2016; Goncharenko et al., 2007; Matsushita, 1959). The storm-time increase in TEC in Figure 10 can result from several competing physical processes that are mainly driven by variations in electric fields of different origins and changes in thermospheric neutral winds composition and circulation. During geomagnetic disturbances, electric fields of magnetospheric origin map to high-latitude ionosphere along the geomagnetic field lines and generate various disturbance electric fields that include prompt-penetration electric field (PPEF) (e.g., Kelley et al., 1979), disturbance dynamo electric field (DDEF) (Blanc & Richmond, 1980) and the subauroral polarization stream (SAPS) electric field (e.g., Foster & Burke, 2002). The effects of these electric fields along with the mechanical effects of the thermospheric wind (Balan et al., 2010) are responsible for the redistribution of ionospheric plasma during a geomagnetic storm and could be contributing to the TEC enhancements seen over Europe.

It has also been observed that minor geomagnetic storms can cause anomalous positive storm effects on the mid-latitude ionosphere under extremely low solar activity conditions that are comparable with the effects of strong geomagnetic storms occurring under higher solar activity conditions (e.g., Buresova et al., 2014; Cander & Haralambous, 2011; Sunda et al., 2013). In the study by Buresova et al. (2014) it was observed that under low solar activity conditions between 2007 and 2009, the changes in F2 layer critical frequency (fof2) could reach up to 100% as compared to its 27-day running mean at mid-latitude stations during minor geomagnetic disturbances. In a recent study by Rajesh et al. (2020), extreme and unexpected TEC enhancements were observed globally during a weak geomagnetic storm on August 5, 2019. Their results showed that despite of Kp values remaining at about 5+, the TEC values over the European region on August 05 almost tripled with respect to the reference value on the pre-storm day (see Figures 4 and 5 of Rajesh et al. (2020)). They reported that the mid-latitude TEC enhancements over Europe on August 05, 2019 seems to be unrelated to the PPEF driven enhanced equatorial fountain but appears to result from the plasma accumulation associated with the storm-induced equatorward surge of neutral winds. This hypothesis is supported in our results from the observations of equatorial propagating LSTIDs. The mechanisms behind this unexpected and intense positive storm effect feature at mid-latitudes should be investigated further in more detail. However, understanding the role of each mechanism behind the positive storm effect feature at mid-latitudes remain the most unpredictable feature of geomagnetic storms and is beyond the scope of the current study at present.