© I. Tsagouri et al., Published by EDP Sciences 2018

1 Introduction

As the societal aspirations for a technologically sophisticated future are increasingly influenced by Earth's near-space environment, the modern society depends more and more on the accurate prediction of the space weather impact on it. In this context, the forecasting of the conditions in the Earth’s upper atmosphere receives special attention. This area plays an important role in ground-based and satellite radio communication and navigation and its density determines the orbit and the lifetime of low-Earth orbiting (LEO) satellites. Two parameters of special interest are the peak F region electron density (NmF2) or the critical frequency (foF2) and the total electron content (TEC). The NmF2/foF2 is related to the maximum usable frequency (MUF) for oblique propagation of radio waves, while TEC is a significant parameter for phase delay effects on the ground-to-satellite navigation signals (Fuller-Rowell et al., 2000).

A key challenge in forecasting the ionospheric behaviour concerns the ionospheric response to geomagnetic storm events. During periods of enhanced geomagnetic activity, the highly variable solar wind and the magnetospheric energy inputs to the Earth's upper atmosphere in the high latitudes activate a variety of mechanisms. These mechanisms include photochemical and chemical reactions, as well as dynamical and electrodynamical processes that drive exchange and transportation of mass, momentum and energy between the elements of the magnetosphere-ionosphere-thermosphere system. Moreover, the interconnection between southward interplanetary magnetic field (IMF) and the Earth's magnetic field leads to prompt penetration of zonal electric fields in the equatorial latitudes (Kelley et al., 2003) that cause redistribution of the ionospheric plasma upward and poleward (Tsurutani et al., 2004). As a result of the combined impact, the ionospheric structure is substantially altered at global scale to compose a complex sequence of ionospheric disturbances. Among others, large-scale increases and decreases in the peak electron density, the so called positive and negative storm effects respectively that appear with strong seasonal, latitudinal and local time dependence hold a critical role (see for example the review articles presented by Rishbeth, 1991; Prölss, 1995; Buonsanto, 1999 and Mendillo, 2006).

Despite the complexity, the operational requirements pushed forward the development of relevant ionospheric forecasting models and algorithms. Significant progress within operational environments has been achieved mainly through empirical approaches (e.g., Tsagouri et al., 2013). In Europe, TEC forecasting products and services are made available by the Ionosphere Monitoring and Prediction Center (IMPC) of German Aerospace Center (DLR, http://impc.dlr.de/), which is the successor of the Space Weather Application Center − Ionosphere (SWACI) service. These include regional and global TEC forecasting maps provided one hour in advance. The forecast is calculated as a sum of the actual European TEC map and a weighted sum of the temporal TEC gradient of the previous hour and the temporal gradient of the previous day at the same time. The IMPC services are also integrated with ESA/SSA Space Weather Service Network (SWE) (http://swe.ssa.esa.int/). On the other hand, a full range of regional foF2 forecasting products and services up to 24 h in advance are routinely provided by the European Digital Upper Atmosphere Server (DIAS, http://dias.space.noa.gr). These are supported by the implementation of the Solar Wind driven autoregression model for Ionospheric Forecast (SWIF − Tsagouri & Belehaki, 2006, 2008; Tsagouri et al., 2009; Tsagouri, 2011; Tsagouri & Belehaki, 2015) and include: warnings for upcoming storm-time ionospheric disturbances, single site forecasts for all DIAS locations and regional forecasting maps for the European area. DIAS forecasting services are also integrated with ESA/SSA/SWE to provide ionospheric forecasts over whole Europe for HF propagation users (http://swe.ssa.esa.int/web/guest/dias-federated), as part of the European Ionosonde Service (EIS).

The SWIF’s concept relies on the combination of two distinct algorithms available for ionospheric storm and non-storm conditions, respectively. The storm time component is based on the empirical Storm Time Ionospheric Model (STIM, Tsagouri & Belehaki, 2006, 2008), which assumes the IMF disturbances observed at the L1 point by ACE spacecraft as drivers for upcoming ionospheric foF2 storm-time disturbances. The forecasting efficiency of the model has been assessed in several occasions and particularly successful performance has been identified under the occurrence of ionospheric disturbances associated to interplanetary coronal mass ejections (ICMEs) detected at L1 point (see Tsagouri & Belehaki, 2015 for a comprehensive description of the SWIF model and the performance tests).

In its present formulation, SWIF has been designed and implemented to serve foF2 forecasting applications. In this paper, we explore SWIF’s capabilities, particularly the STIM component, in forecasting the TEC storm-time variations. This investigation is driven by: (i) the argument that most models of F region disturbances concentrate on the peak density are fully relevant to TEC (Mendillo, 2006); (ii) recent results that support potentially improved predictions of TEC deviations during space weather events based on solar wind measurements (Borries et al., 2015); (iii) promising evidence on the efficiency of the SWIF’s alert criteria in predicting TEC storm-time disturbances in long term investigations (Bergeot et al., 2013). To meet the goal, the ionospheric disturbances are comparatively analysed in conjunction to the solar wind conditions at L1 point as they are monitored by the ACE spacecraft. The ionospheric disturbances are determined through the comparison between actual observations of the foF2 critical frequency and GNSS-TEC estimates obtained over European locations with their corresponding median values during eleven storm events (see Data presentation in Section 2 for a detailed description of the datasets used here). The results reveal similarities and differences in the response of the two ionospheric characteristics with respect to the local time and the latitude of the observation point, driving the upgrade of the SWIF’s modeling technique in forecasting the storm-time TEC variation several hours in advance. This part of the analysis is discussed in Section 3. Preliminary considerations of the SWIF’s upgraded formulation and results on the evaluation of the model's efficiency in TEC prediction are presented in Sections 4 and 5, respectively. A comprehensive summary/discussion of the results and conclusive notes are provided in Sections 6 and 7, respectively.

2 Data presentation

This investigation exploits foF2 and hmF2 records of hourly resolution obtained by ground-based Digisondes and GNSS-derived TEC estimates, also of hourly resolution. The Digisonde observations were retrieved from the Global Ionosphere Radio Observatory (GIRO − http://giro.uml.edu/) data repository (Reinisch and Galkin, 2011) and they mainly come as the result of automatically scaled ionograms. TEC estimates used here are based on data from GNSS receivers located close to the Digisondes. In particular, at each station vertical TEC (vTEC) values are calculated from data extracted from Receiver Independent Exchange Format (RINEX) files with 30 s sampling, using the single station solution proposed by Ciraolo (2005) and Ciraolo et al. (2007). Dual-frequency carrier-phase and code-delay GNSS observations are combined to obtain ionospheric observables related to the slant total electron content (sTEC) along the satellite-receiver line-of-sight. Only data files received from GNSS satellites with elevation larger than 70° are considered in this analysis. The sTEC has been converted to vTEC at Ionosphere Pierce Point (IPP) by applying the ionospheric thin shell model. The conversion can be done using the relation given by where χ is the zenith angle of a GNSS satellite at IPP.

The ionospheric disturbances are determined through the comparison between observed and monthly median values that are obtained from actual observations recorded within each calendar month. The monthly medians were estimated for each hour/month/year.

The present analysis is carried out for the middle latitude European region with data obtained over the ionospheric locations that are listed in Table 1 (also mapped in Fig. 1). It applies to all intense storm events (min Dst ≤ −100 nT) occurred in the present solar cycle 24 that are associated with CME-related solar wind flows (e.g., sheath fields or the ejecta itself) in the near-Earth solar wind. The distinction is based on the entries of the ICMEs catalogue that is available at http://www.srl.caltech.edu/ACE/ASC/DATA/level/icmetable2.htm (Richardson and Cane, 2010). The storm events analysed here are listed in Table 2.

Fig. 1 A European map indicating the locations of the observational points listed in Table 1. The locations of the Digisondes and the GNSS receivers are noted with red and blue circles, respectively.

3 Analysis of the TEC storm-time response with respect to the SWIF’s concept

Of central importance within the SWIF’s concept is the determination of the storm onset time at L1 point. To this effect, SWIF exploits the STIM’s alert detection algorithm (ADA) that is able to identify storm (alert) conditions and therefore to activate the storm-time algorithm i.e., the STIM. Storm (alert) conditions are identified by applying quantitative criteria to IMF-total magnitude Bmag, its rate of change, dBmag/dt and the IMF-Bz component. The criteria were derived mainly through superposed epoch analysis that was carried out over several intense storm events occurred in solar cycle 23. The set of the criteria were initially introduced by Tsagouri & Belehaki (2008) and include: (i) an increase in the IMF–Bmag: time derivative values greater than 3.8 nT/h or absolute values greater than 13 nT, and (ii) a southward turning of the IMF–Bz component (Bz < −10 nT) for at least 3 h.

In this context, we investigate the interplanetary conditions for the specific set of storm events analyzed here as the first step in testing the validity of the SWIF’s concept in forecasting the vTEC storm-time response. In Figure 2, we present the results of the superposed epoch analysis that was applied again to the IMF-total magnitude, its rate of change and the Bz-component. The Dst index is also included to verify geomagnetic storm conditions. The results demonstrate that the SWIF’s criteria are still fully applicable.

The next key issue for SWIF is the dependence of the ionospheric response on the latitude and the local time (LT) of the observation point at the storm onset. Following the determination of the storm onset at L1 point, STIM algorithm determines the ionospheric storm onset and quantifies/formulates the storm-time disturbances over single locations. To this effect, the model exploits empirical expressions that take into account the latitude of the observation point and its LT at the storm onset. The model considers two latitudinal zones, one for latitudes greater than 45° N (middle-to-high, MtH) and one for latitudes between 30° and 45° N (middle-to-low, MtL), and four LT sectors: morning (0300–0600 in LT), prenoon (0700–1200 in LT), afternoon (1300–1800 in LT) and evening (1900–0200 in LT).

The validity of STIM's background assumptions are further explored below for vTEC forecasting purposes (see Figs. 3–6), through superposed epoch analysis that is applied to the vTEC storm-time variation for the set of the storm events presented in Table 2. Each of Figs. 3–6 includes the results for one LT sector (i.e., morning, prenoon, afternoon or evening) and both latitudinal zones (MtH and MtL in the top and bottom panel, respectively). For comparison purposes, we also include in the plots the corresponding variation of the foF2 for the same set of events and the empirical formulation of the storm-time response that the original SWIF’s expressions can provide. In this part of the analysis, the ionospheric response was investigated over all ionospheric locations reported in Table 1 except for Dourbes and San Vito. These two stations were kept out in order to serve model validation experiments described later as test stations. The ionospheric response is expressed by the ratios foF2/foF2median and vTEC/vTECmedian for foF2 and vTEC, respectively following the SWIF’s original formulation that is designed to provide a correction factor to the normal ionosphere.

The results verify significant disturbances in ionospheric plasma properties that are apparent in both foF2 and vTEC. According to the obtained trends one may argue that:

1. Concerning the LT dependence, the ionospheric storm-time response starts almost simultaneously in foF2 and vTEC. Moreover, the vTEC tends to follow qualitatively the trends present in the foF2 pattern. More precisely, positive storm effects tend to predominate in the daytime response of both foF2 and vTEC (see the trends in Figs. 3 and 4), while negative storm effects characterize the nighttime response in all cases. However, there are noticeable differences in the quantitative characteristics of the disturbances in each case. The most important one applies in the intensity of positive storm effects that appears significantly greater in the vTEC case. It is also interesting to note that positive storm effects at higher latitudes in case the observation point at storm onset time is located in the morning sector are present only for vTEC. The vTEC positive phase is more intense in the lower latitudes especially in the afternoon hours (see Figs. 3 and 4).

2. Concerning the SWIF’s prediction efficiency: the model in its original version is able to reproduce satisfactorily the intensity and the duration of significant negative storm effects in foF2 in all cases. With respect to vTEC, it fails to reproduce the significant negative response that follows the positive phase in the morning sector (see. Fig. 3), which is apparent only in the vTEC variation. With respect to positive effects, although it tends to underestimate daytime ionization increases observed in lower latitudes, the model seems to work satisfactorily well for foF2. Nevertheless, some adjustments should be anticipated in the vTEC formulation (see also point i above).

Fig. 2 The results of superposed epoch analysis applied to the IMF-total magnitude (Bmag), its rate of change (dBmag/dt) and the Bz-component. The Dst index is also included to verify geomagnetic storm conditions. The vertical red line indicates the storm onset time (0 h).

Fig. 3 The results of the superposed epoch analysis for the vTEC and foF2 storm-time response in case the observation point is located in the morning LT sector at the storm onset, for MtH (top panel) and MtL (bottom panel) locations. The predictions that the original SWIF’s expressions can provide is also included (green line). Error bars represent standard deviations. The vertical black line at 0 hrs indicates the storm onset time.

Fig. 4 Same as Figure 3 in case the observation point is located in the prenoon LT sector at the storm onset.

Fig. 5 Same as Figure 3 in case the observation point is located in the afternoon LT sector at the storm onset.

Fig. 6 Same as Figure 3 in case the observation point is located in the evening LT sector at the storm onset.

4 SWIF’s empirical formulation for the vTEC forecast

Following the findings of the previous sections, and with the focus on the forecast of the vTEC, we propose the adjustment of the SWIF’s empirical formulation as below.

1. Since the SWIF’s alert criteria are still valid (see Fig. 2), we keep SWIF’s original settings for the determination of the storm onset at L1 point (see Sect. 3 for the description of the alert criteria).

2. For the prediction of the vTEC storm-time response, we exploit the results of the superposed epoch analysis obtained in the previous section that are reproduced in Figure 7 for reasons of clearness. As in case of foF2, the model idea is to apply the determined variations to the normal ionosphere that can be expressed e.g., by monthly medians in order to predict vTEC values.

The vTEC response depends on the local time and the latitude of the observation point. It may be important to note that the vTEC variations in Figure 7 include also the time delay in the ionospheric storm onset in each local time sector (see x-axis that counts the time from the storm onset at L1 point). This approach is slightly different than the original version of the model that includes more detailed formulations for the local time dependence of the foF2 storm-time response in each latitudinal zone (Tsagouri & Belehaki, 2008). This simplification is made here for practical reasons, since the number of storm events with the onset in a specific local time sector is rather small to support a detailed analysis, which remains subject for future improvements in the proposed model. Nevertheless, the reader should note that the overall trends in the time difference between the storm onset at L1 point and the ionospheric storm onset at middle latitudes are still comparable in the two versions of the model (see Fig. 4 in Tsagouri & Belehaki, 2008). This was expected since the ionospheric storm-time response starts almost simultaneously in foF2 and vTEC: positive storm effects appear within 1–3 h after the storm onset, while the time difference for the appearance of negative storm effects ranges from about 2–3 h up to 12–13 h. The minimum is recorded when the observation point is located in the afternoon and evening sectors at the storm onset, while the maximum time delay is recorder when the observation point is located in the morning and prenoon sectors at the storm onset.

In principle, the model is designed to predict the storm-time response from the onset to full development and recovery. The results given in Figure 7 represent the model predictions for a storm event with one single onset. However, in practice the geomagnetic activity may be evolved through successive intensifications that can prolong both the geomagnetic and ionospheric activity for many days. In such cases, the model's predictions are superposed to anticipate all detected onsets. This is how SWIF operates in case of the foF2 in DIAS (see also the example presented in Fig. 9).

Fig. 7 SWIF’s empirical formulation for the vTEC storm time response taking into account the latitudinal zone (i.e., MtH and MtL) and the local time sector of the observation point at the storm onset (i.e., morning, prenoon, afternoon, evening).

5 Performance evaluation tests for the proposed model

In this section, we attempt an assessment of the proposed model's prediction capabilities. The assessment tests described below are carried out on the same set of storm events listed in Table 2. This compromise results from the limitations in the vTEC data availability (see Sect. 6 for more explanation) and in this respect, the validation plan may not be considered as fully independent. Nevertheless, one may argue that the results can still provide useful input regarding the achieved prediction efficiency.

As the SWIF’s (ADA) alert criteria remain the same in the new version of the model, we exploit the results of the DIAS automated algorithm for the determination of the storm onset time. To some extent, one may argue that this input serves the assessment of key components of the proposed model in operational mode. The storm onset time for each of the storm events analyzed here is presented in Table 2 and is also available through DIAS (at http://dias.space.noa.gr, under Alerts' Archive list after registration). Multiple onsets are detected for several of the storm events analyzed here. It should be noted that only one of the storm events, the one occurred in the time interval 16–20 March 2013, was “missed” by SWIF’s ADA (see Tsagouri & Belehaki, 2015 for the evaluation of the results of the SWIF’s ADA).

The accuracy of the upgraded SWIF’s performance in forecasting the vTEC response was then tested over two GNSS locations: DOUR, which is close to Dourbes in Belgium and USAL, which is close to San Vito in Italy. As mentioned above, these two stations were not included in the analysis for the determination of the model's empirical formulation to help, in the maximum possible extent, an unbiased evaluation of model's capabilities. For comparison purposes, in this part of the analysis we also include the vTEC output of the IRI (International Reference Ionosphere) model. In particular, we use the output of the IRI2012 version of the model (e.g. Bilitza et al., 2014), available for on line calculations at https://omniweb.gsfc.nasa.gov/vitmo/iri2012_vitmo.html. Key input parameters for the model include: Upper limit for Electron content = 2000 km; F peak model = CCIR; Ne Topside = IRI01-corr; foF2 Storm model = on; Bottomside Thickness = ABT-2009; F1 occurrence probability = Scotto-1997 no L; foE Storm model = on; D-Region Ne = IRI-95; Auroral boundary = off.

As indicative examples of the capabilities of the proposed SWIF’s vTEC version, we present in Figures 8 and 9 the model's predictions obtained for the time intervals 23–26 April 2012 and 21–25 June 2015, respectively. The corresponding GNSS-derived vTEC estimates, monthly medians and IRI2012 output are also included for comparison purposes. The results demonstrate the models efficiency in capturing the ionospheric storm onset and the ionospheric variation on the main storm day. The IRI2012 tends to reproduce only the median variation in the first case, while in the second event it tends to underestimate the observations with no clear sensitivity to the storm effects.

For a more comprehensive overview of the performance of the proposed model, in Figure 10, we present the scatter plots between vTEC modeled and estimated values obtained over the two test locations for all the events analyzed here. The results demonstrate satisfactory agreement between model’s predictions and GNSS-derived estimates. This agreement is also verified by the findings of Figure 11, which provides the distributions of modeled and estimated values for both test locations in a boxplot format. The distributions of the monthly median values are also included in each case for comparison purposes. The model's predictions reproduce successfully the distributions of the GNSS-derived estimates in both test locations and more effectively than the median values, indicating that the model can efficiently capture the observed storm-time disturbances. In Table 3, the performance of the proposed model is assessed through the estimation of several metric parameters: the Mean Error (ME), the Mean Absolute Error (MAE), the Mean Relative Error (MRE) and the Root Mean Square Error (RMSE) and the relative improvement over climatology. The latter is calculated by using the following formula (see for example Araujo-Pradere & Fuller-Rowell, 2002; Tsagouri, 2011; Tsagouri & Belehaki, 2015):

Based on the results presented there, one may argue that the model's predictions do not suffer by significant biases (see the ME over both locations), while the MRE is really small in all cases. In overall, the performance of the model tends to be more successful for DOUR. In this case, the model's error in terms of both MAE and RMSE is comparable to 3 TECu that in turn, is comparable to the TEC measurement error at middle latitudes (e.g. Zhang et al., 2010 and references therein). Significant improvement (>10%) over climatology is recorded in both cases. The same tests were also performed for IRI2012 for comparison purposes. The corresponding results are also included in Table 3. The calculated improvement over climatology demonstrates that the IRI2012 predictions are comparable to the climatological estimates (i.e., medians).

Fig. 8 SWIF’s predictions for the vTEC storm time response during the time interval 23–26 April 2012 for DOUR and USAL locations. The corresponding GNSS-derived vTEC estimates, monthly medians and IRI2012 output are also included for comparison purposes. The storm onset time was determined in the afternoon sector (23/04/2012, 14:00 UT), so that the vTEC response was forecasted 4 and 6 h before its onset at DOUR and USAL, respectively. In this example, the IRI2012 tends to reproduce only the median variation.

Fig. 9 Same as Fig. 8 for the time interval 21–25 June 2015. This was a storm event with multiple onsets determined in the afternoon, evening and morning sectors (22/06/2015 18:00 UT; 22/06/2015 23:00 UT; 23/06/2015 04:00 UT). The vTEC response was forecasted 6 and 8 h before the first storm onset at DOUR and USAL, respectively. The IRI2012 tends to underestimate the observations with no clear sensitivity to the storm effects.

Fig. 10 The scatter plots between vTEC modeled and estimated values obtained over the DOUR (top panel) and USAL (bottom panel) locations for all the events listed in Table 2. The linear regression is also plotted with the red line in both cases.

Fig. 11 Boxplots representing the distributions of vTEC estimates, modeled and median values for DOUR (top panel) and USAL (bottom panel) locations. The box has lines at the lower quartile, median (red line) and upper quartile values. Whiskers extend from each end of the box to the adjacent values in the data; in our case to the most extreme values within 1.5 times the interquartile range from the ends of the box. Outliers (e.g., data with values beyond the ends of the whiskers) are displayed with red crosses.

6 Summary and discussion

In this paper, we attempt the expansion of the capabilities of the SWIF model, particularly its storm component in forecasting the vTEC storm-time variation. The SWIF model runs in the DIAS system since 2009 and provides ionospheric forecasting products and services with respect to the foF2 critical frequency for the DIAS and ESA/SSA/SWE users. The model is driven by solar wind disturbances detected at L1 point through the analysis of the ACE real-time measurements. The present work was motivated by recent results that indicated potentially improved predictions of TEC disturbances during space weather events based on solar wind measurements (Borries et al., 2015), but also by the findings received through the exploitation of the SWIF’s output in the investigation of the space weather impact on TEC variation in long-term perspective. More precisely, Bergeot et al. (2013) demonstrated that the temporal and spatial evolution of TEC disturbances detected in solar cycle 23 and specified by SWIF’s storm onset conditions were consistent with well established phenomenological models. These findings provided a first indication of the potential efficiency of the SWIF model in predicting TEC storm-time deviations.

The expansion of the model's forecasting capabilities was based on the comparative analysis of foF2 and vTEC storm-time disturbances during eleven storm events occurred in the present solar cycle 24. Since there are some evidence available today that indicate dependence of the ionospheric response on the solar wind drivers of the storms (e.g. Tsagouri et al., 2017) and moreover, the performance of the SWIF model had been proved to be particularly successful in forecasting the impact of intense CME-related events (Tsagouri & Belehaki, 2015), the focus in the present analysis was given in this type of storm events in order to build progress on a solid and clear basis. The distinction of the interplanetary drivers of the storms is compatible with the list of the ICMEs that is available at http://www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm (Richardson & Cane, 2010).

The ionospheric disturbances were determined through the comparison between observed and monthly median values that were obtained from actual observations recorded within each calendar month. The vTEC estimates used here are based on data from GNSS receivers located close to the Digisondes' locations to ensure fair comparison between foF2 and vTEC disturbances. At each station, vTEC values were calculated using the single station solution proposed by Ciraolo (2005) and Ciraolo et al. (2007). This method assumes that the ionosphere is a thin layer around 300 km. The overall distributions of the hmF2 observations obtained in the Digisonde locations during the storm events analyzed here show that in the majority of the cases, the height of the peak density lies within the assumed range in all reference and test locations (see Fig. 12). This may be considered as a rough indication of the reliability of the vTEC estimates that were used here. Further, one may argue that the exploitation of results provided by a single site solution assigns a significant advantage to this work, since these results could be considered as more representative for local vTEC estimates, being free of extra interpolation errors that could be introduced by any mapping technique.

The vTEC storm-time disturbances were found to be compatible with SWIF’s specifications for the determination of the storm onset at L1 point, indicating that the SWIF’s ADA and the alert it can provide could be considered reliable tools in forecasting the vTEC storm-time variation. This variation was further investigated through superposed epoch analysis in combination with foF2 storm-time disturbances that evolve simultaneously, taking into account the latitude and the local time of the observation point at storm onset. The results revealed similarities and differences in the response of the two characteristics. As a first outcome, the evidence indicate that in principle the morphology of the ionospheric storm effects is maintained in both foF2 and vertical TEC. This is rather expected since the physical causes for NmF2/foF2 and vertical TEC storm-time variations are the same (Mendillo, 2006). Nevertheless, there are quantitative and qualitative differences in their response. Most probably, these differences reflect the ability of each parameter to characterize the ionosphere during storms (Mendillo, 2006), but also the relative importance of the background mechanisms with altitude. One has to keep in mind that vertical TEC represents the integral over the entire ionospheric electron density. Such an integral has its maximum contribution from the F2 layer, with approximately 2/3 of the total electron content coming from regions above the altitude of the peak density, hmF2 (e.g., Mendillo, 2006; Belehaki & Tsagouri, 2002).

Today, there are several mechanisms proposed in the literature for the explanation of the observed storm-time effects. According to the Prölss phenomenological model for the dependence of the F-layer disturbances on local time (e.g., Prölss, 1995, 2011), positive storm effects are attributed to travelling atmospheric disturbances (TADs) and meridional winds, while negative storm effects are attributed to changes in the neutral gas composition that follow the injection of solar wind energy to the polar upper atmosphere. Any new intensification in the solar wind energy input results in a new intensification of the ionospheric activity (Tsagouri et al., 2000). In this scenario, positive storm effects in the middle latitudes are expected to be observed in the daytime ionosphere shortly after the storm onset (i.e. within 1–3 h) following the propagation of the TADs from the auroral to lower latitudes and the imposed changes to the wind circulation. On the other hand, neutral composition changes are evolved in different time scale and are mainly detected in the post-midnight hours. As a result, the negative storm effects are observed when the observation point is located (or rotated) into the nighttime hemisphere that conditionally can provide further delays in the ionospheric response. This scenario may explain the differences in the time lag of the ionospheric response in different LT sectors that are evident in Figure 7.

More generally, the proposed mechanisms include upwelling/downwelling of the gas due to storm-induced thermospheric circulation, penetration of electric fields of magnetospheric or interplanetary origin into the ionosphere (e.g. the super-fountain mechanism proposed by Tsurutani et al., 2008) and plasma fluxes from the plasmasphere (Danilov, 2013). Presently, it is believed that the electron concentration at heights around the F2 layer maximum tends to be more sensitive to changes in neutral composition, thermospheric gas temperature and horizontal winds, whereas electrodynamics (i.e., penetrating electric fields and plasma fluxes from the plasmasphere) is more important at higher altitudes (400–800 km) (Danilov, 2013). Taking into account that the electrodynamical processes reported above result mainly in ionization increases, this argument may explain the detected differences in the occurrence and the intensity of positive storm effects: they appear more frequently and significantly greater in the vTEC case, maybe due to the greater contribution from the topside ionosphere in the estimation of vTEC; they are more intense in the lower latitudes, which is especially true in the afternoon hours (the so-called “dusk effect”) as a result of the combined effect of penetrated electric fields in the equatorial latitudes and downward plasma fluxes from the plasmasphere (Mendillo, 2006).

Although the ADA remains in place for the forecasting of vTEC, the identified differences between foF2 and vTEC response made it clear that one should look for a new empirical formulation in case of the vTEC response. This task was driven here by the results of the superposed epoch analysis. According to the new formulation, the vTEC storm-time response depends also on the latitude and the local time of the observation point at the storm onset time at L1 point and can be predicted from 1 to 12–13 h in advance (see Fig. 7), depending on the local time of the observation point at storm onset at L1.

The upgraded version of SWIF was also preliminary evaluated in performance during the storm events under study over two GNSS test locations: DOUR in Belgium and USAL in Italy. Although the validation tests that were carried out here may not be considered as totally independent and such an approach remains a strong requirement for any future developments, we believe that our findings provide indicative results regarding the forecasting capabilities of the proposed model. In particular, the performance of the model has been proven to be more successful in the higher latitudes, giving a prediction error that is comparable to 3 TECu, which in turn is comparable to the TEC measurement error at middle latitudes. Moreover, the results demonstrate significantly improved performance with respect to standard prediction approaches such as the median estimates. Concerning the RMSE, the improvement was 17% and 37% for USAL and DOUR, respectively. More successful performance was also determined in comparison to the standard IRI2012 model: while the proposed model was able to capture the ionospheric disturbances during the main storm day, IRI2012 provided predictions that were rather comparable to the median estimates. At this point, one may note that the upper limit for the estimation of the total electron content by IRI is only 2000 km. This limitation may partially explain some of the discrepancies. Nevertheless, taking into account that the maximum contribution in the vTEC estimate comes from the F2 layer, the results obtained here may still be considered as indicative of the general trends.

As a shortcoming of the present work, one may indicate the relatively small number of the storm events analyzed for the development of the new model, but in this respect, one should also consider the limitations that are imposed by the availability of vTEC data. The selected set of the events include all storm events that occurred in the present solar cycle 24 and meet the requirements of the present analysis (i.e., they are intense and CME-driven events). As it was pointed out previously, here we used GNSS-vTEC estimates that were provided by a single station solution as most appropriate for the analysis of storm effects in local perspective. However, the extraction of these results is not a straightforward task, since the code execution is heavily time consuming. On the other hand, GNSS-vTEC estimates provided by regional or global maps come with significant discrepancies as they depend strongly on their calculation method. In addition, they usually provide a smoothed representation of the ionospheric conditions through the background interpolation that could not effectively support the analysis of storm-time effects. In any case, the further we go back in time, the less dense is the network of the available receivers and the lower is the quality of the data and the accuracy of the GNSS-vTEC estimates.

Further on the previous point, the availability of the vTEC data is also a big challenge for the operational implementation for the proposed storm-time model. As it requires reliable GNSS-vTEC estimates available in real time to provide an estimation of the background ionosphere (e.g. by running medians as in case of foF2), key future plans apply to the adjustment of the algorithm used here towards its real time exploitation.

Another point of further consideration regarding the potential implementation of the proposed model in operational mode may be the availability of real-time solar wind measurements at L1 point. Following the settings of the SWIF’s operational version for foF2, the present analysis was based also on the analysis of IMF observations from ACE spacecraft that is still able to provide valuable input. Future upgrades may include the adjustment of the model's settings to new monitoring capabilities at L1, such as the NOAA's DSCOVR (Deep Space Climate Observatory), which comes as the operational replacement for the real-time solar wind instruments on the ACE satellite (Knipp & Biesecker 2015; Constable, 2016).

Future plans include also the upgrade of the SWIF’s non-storm component for forecasting vTEC during quiet ionospheric conditions. This may be possible through the implementation of autoregressive technique, following the philosophy of the TSAR (Time Series AutoRegressive) algorithm (Koutroumbas et al., 2008) that presently supports SWIF’s forecasts provided by DIAS for the foF2 in non-alert conditions at middle latitudes. For the forecast of vTEC, the time series that will be used are sampled every 30 s and the methodology relies on the use of autoregressive (AR) modeling. In the sequel, we describe briefly the proposed general forecasting strategy and then we discuss the way it can be utilized in the framework of vTEC forecasting problem. Assuming that the process that generates an observed time series {x(n), n = 0,1,…} is described by an AR model, the estimation of the value of the time series s steps ahead, x(n+s), based on the M measurements before the current time step is expressed as follows where M is the order of the AR model, w = [w₀,w₁,...,w_M]^T is its parameter vector (see e.g., Kalouptsidis, 1997) and x_n = [x(n), x(n-1),..., x(n−M)]^T. For fixed M, the best w results from the solution of a linear system of equations (for more details, see Kalouptsidis, 1997; Koutroumbas et al., 2008). Note that both M and w are crucial for the complete determination of an AR model. Based on a time series

of length l, with l >> M, the estimation of the order of the model, M, as well as the parameter vector w, can be carried out in the following two candidate ways: (a) estimate first the best (in the mean square error sense) M and then, for the best M, determine the corresponding parameter vector w and (b) estimate M and w simultaneously, employing the concept of sparsity (Themelis et al., 2012). Let us describe the above methodologies in a little more detail. First, we split Y into the following two sets, i.e., X₁ = {x(0), x(1), …, x(l/2)} and X₂ = {x(l/2 + 1),…, x(l)}¹. According to the first methodology, we apply for a range of values of M, [M_min, M_max], the following procedure: we determine the parameter vector w, based on the set X₁ and then we measure the associated mean square error, MSE_i, i = M_min,…,M_max, based on the set X₂. Then, among all the MSE_i's, we select the minimum one and adopt the corresponding parameter vector w.

According to the second methodology, we consider a parameter vector w, of relatively large size and we estimate it based on the set X₁, adopting specific algorithms (e.g., Themelis et al., 2012) that impose sparsity on the entries of w; that is, they push several values of w to zero. Clearly, this is an implicit way of determining the order of the model. In addition, in contrast to the previous case, the entries of w that are pushed to zero need not be consecutive.

Taking into account that ionosphere changes dynamically, the above methodology is modified as follows, in order to cope with the vTEC forecasting. The order M and the parameter vector w of the AR model are re-estimated (with either of the two methodologies) at the beginning of each calendar month, taking into account the vTEC values of the previous couple of months.

Fig. 12 The distributions of the hmF2 values that were recorded in the Digisonde locations included in Table 1 during the storm events analyzed here.

7 Conclusions

In this paper, a new model for the short-term forecasting of the vTEC parameter was proposed, with the focus in the forecast of the storm-time effects. The development of the new model was based on the concept of the SWIF model, which is operationally implemented in the DIAS system and provides a full set of ionospheric forecasting products and services over Europe in terms of the foF2 critical frequency. The model is triggered by detected CME-related disturbances in the solar wind at L1 point and estimates the ionospheric storm-time response 1–13 h in advance, expanding significantly current forecasting capabilities for vertical TEC. The predictions are based on empirical formulations that take into account the latitude and the local time of the observation point at the storm onset (at L1). The performance of the model was preliminary evaluated over European locations in comparison to GNSS-vTEC estimates, to climatological estimates (i.e., monthly medians) and to IRI2012 predictions. The results are encouraging, demonstrating significant improvement over climatology and relatively small prediction error, especially in the higher latitudes, where the improvement was found up to 37% and the error was found comparable to 3 TECu. Future improvements are mainly considered towards the operational implementation of the proposed model. Key challenges include the availability of vertical GNSS-TEC estimates over single locations in real time and the development of an algorithm for the vTEC forecast during quiet conditions in support of the storm-time model. Improvements in the availability of the GNSS-vTEC data would also help the analysis and the evaluation of the model's performance in more storm events in support of future developments.

Acknowledgments

The Dst index are provided by the World Data Center for Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp/dstdir/). Special acknowledgements are due to Luigi Ciraolo for providing the source code for the calculation of TEC values from RINEX files. RINEX files from the European IGS stations have been downloaded from the IGS ftp archive. The authors thank the ACE MAG instrument teams and the ACE Science Center for providing the ACE data. This investigation uses ionospheric data from Global Ionospheric Radio Observatory (GIRO, http://spase.info/SMWG/Observatory/GIRO). The authors are grateful to Anna Belehaki for many useful discussions. This work is supported by the NATO Science for Peace and Security project #984894 and partly by the European Commission Horizon 2020 project #776011 TechTIDE “Warning and Mitigation Technologies of Travelling Ionospheric Disturbances Effects”. I.T. was in part supported by ESA SSA Contract 4000113184/15/D/MRP. The editor thanks two anonymous referees for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 A European map indicating the locations of the observational points listed in Table 1. The locations of the Digisondes and the GNSS receivers are noted with red and blue circles, respectively.
In the text

	Fig. 2 The results of superposed epoch analysis applied to the IMF-total magnitude (Bmag), its rate of change (dBmag/dt) and the Bz-component. The Dst index is also included to verify geomagnetic storm conditions. The vertical red line indicates the storm onset time (0 h).
In the text

Fig. 3 The results of the superposed epoch analysis for the vTEC and foF2 storm-time response in case the observation point is located in the morning LT sector at the storm onset, for MtH (top panel) and MtL (bottom panel) locations. The predictions that the original SWIF’s expressions can provide is also included (green line). Error bars represent standard deviations. The vertical black line at 0 hrs indicates the storm onset time.

In the text

	Fig. 4 Same as Figure 3 in case the observation point is located in the prenoon LT sector at the storm onset.
In the text

	Fig. 5 Same as Figure 3 in case the observation point is located in the afternoon LT sector at the storm onset.
In the text

	Fig. 6 Same as Figure 3 in case the observation point is located in the evening LT sector at the storm onset.
In the text

	Fig. 7 SWIF’s empirical formulation for the vTEC storm time response taking into account the latitudinal zone (i.e., MtH and MtL) and the local time sector of the observation point at the storm onset (i.e., morning, prenoon, afternoon, evening).
In the text

Fig. 8 SWIF’s predictions for the vTEC storm time response during the time interval 23–26 April 2012 for DOUR and USAL locations. The corresponding GNSS-derived vTEC estimates, monthly medians and IRI2012 output are also included for comparison purposes. The storm onset time was determined in the afternoon sector (23/04/2012, 14:00 UT), so that the vTEC response was forecasted 4 and 6 h before its onset at DOUR and USAL, respectively. In this example, the IRI2012 tends to reproduce only the median variation.

In the text

Fig. 9 Same as Fig. 8 for the time interval 21–25 June 2015. This was a storm event with multiple onsets determined in the afternoon, evening and morning sectors (22/06/2015 18:00 UT; 22/06/2015 23:00 UT; 23/06/2015 04:00 UT). The vTEC response was forecasted 6 and 8 h before the first storm onset at DOUR and USAL, respectively. The IRI2012 tends to underestimate the observations with no clear sensitivity to the storm effects.

In the text

	Fig. 10 The scatter plots between vTEC modeled and estimated values obtained over the DOUR (top panel) and USAL (bottom panel) locations for all the events listed in Table 2. The linear regression is also plotted with the red line in both cases.
In the text

Fig. 11 Boxplots representing the distributions of vTEC estimates, modeled and median values for DOUR (top panel) and USAL (bottom panel) locations. The box has lines at the lower quartile, median (red line) and upper quartile values. Whiskers extend from each end of the box to the adjacent values in the data; in our case to the most extreme values within 1.5 times the interquartile range from the ends of the box. Outliers (e.g., data with values beyond the ends of the whiskers) are displayed with red crosses.

In the text

	Fig. 12 The distributions of the hmF2 values that were recorded in the Digisonde locations included in Table 1 during the storm events analyzed here.
In the text