© X. Meng et al., Published by EDP Sciences 2016

1. Introduction

The terrestrial ionosphere is significantly influenced by solar and interplanetary conditions. Disturbances on the Sun and structures in the solar wind may cause geomagnetic activities and consequently ionospheric storms, which perturb the neutral composition, temperature, wind, electric fields, and electron densities in the ionosphere (Buonsanto 1999; Mendillo 2006). Such perturbations can affect many technological systems that human activities rely on, including telecommunication systems operating on the ground and in space, power grids and oil pipelines, satellite operations, global navigation satellite systems, and communication with aircrafts. Therefore, predicting ionospheric conditions, especially during geomagnetically disturbed periods, is crucial to reduce or even avoid damages in the technological systems. The ultimate goal would be long-term operational forecasts of ionospheric weather.

Numerical models of the ionosphere offer opportunities of ionospheric forecasts. Currently, both data-assimilative models (Schunk et al. 2004) and fully physics-based models (Richmond et al. 1992; Roble & Ridley 1994; Millward et al. 1996; Huba et al. 2000; Ridley et al. 2006) are capable of predicting ionospheric conditions. While data-assimilative models, including physics-based data-assimilative models, are mainly applicable for nowcasts and near-term forecasts (Keil 2007; Schunk et al. 2005, 2014), fully physics-based models can provide forecasts with a lead time of several days, or more when coupled with predictive solar and solar wind models. Despite their uncertainties and missing physics (Schunk et al. 2012), the current generation of physics-based global ionospheric models are adequately advanced to permit an analysis of forecasting (Mannucci et al. 2015), given that these models have successfully reproduced the overall ionospheric dynamics and have been validated against measurements under various solar conditions. This is the primary motivation for our study. We aim to understand modeling capabilities through a forecasting analysis so that future forecasts can be improved.

A complete ionospheric forecast consists of forecasting both interplanetary drivers and the corresponding ionospheric conditions. Given that the current predictive heliospheric models are not mature enough to accurately reproduce the solar wind at 1 AU, our forecast study is carried out in two steps. In the first step, we evaluate the forecast performance of ionospheric models and explore their uncertainties. In the second step, we bring in solar wind forecasts to drive ionospheric models and investigate how uncertainties in the solar wind affect ionospheric forecasts. This paper focuses on the first step. In particular, we present a forecasting-oriented analysis that employs the physics-based Global Ionosphere-Thermosphere Model (GITM; Ridley et al. 2006) and produces ionospheric total electron content (TEC) distributions for corotating interaction region/high-speed-stream (CIR/HSS) events in a forecast mode. We call the modeled TEC output a forecast, even though it is generated retrospectively.

TEC is an excellent quantity to forecast, as it is a representative variable of the ionosphere and can be directly validated against measurements. CIRs and associated high-speed solar wind streams (Tsurutani et al. 2006) are known drivers of ionospheric storms especially during the declining phase of a solar cycle. Unlike previous modeling studies of ionospheric and thermospheric responses to CIR/HSS events (Pedatella & Forbes 2011; Burns et al. 2012; Solomon et al. 2012), which utilize fully physics-based models as tools to examine the storm dynamics, we emphasize the forecasting perspective. We perform forecast-mode GITM simulations for four selected CIR/HSS events. In the forecast mode, the simulations are only driven by inputs that are forecastable given a suitable solar wind forecast, and none of the model parameters are tuned from default values. We define a metric to characterize TEC variations relative to pre-storm quiet-time conditions. The TEC metric is displayed on hourly global maps as a possible type of forecast product. The forecasted TECs are compared to the TECs from Global Ionospheric Maps (GIM), which are based on Global Positioning System (GPS) satellite observations. We investigate examples of model-data agreements and disagreements to better interpret GITM and evaluate its forecasting capabilities.

The following content of the paper is divided into four sections. Section 2 introduces the CIR/HSS events to forecast, the GIM TEC data, GITM simulations, and the TEC metric. Section 3 presents the forecast-mode results and compares them to the GIM data. Section 4 discusses possible reasons for model-data agreements and disagreements. Section 5 concludes the paper.

2. Techniques

2.1. CIR/HSS events

For our forecast study, we selected four typical CIR/HSS events in January 2007, April 2011, May 2012, and June 2012, respectively. These events span from the declining phase of Solar Cycle 23 (2007 event) to the ascending phase of Solar Cycle 24 (2011 and 2012 events), and from northern hemisphere winter to summer. They are typical in the sense of conforming to the schematic picture for CIR/HSS events outlined in Tsurutani et al. (2006).

Figure 1 displays the solar wind conditions as well as the SYM-H and AE indices for the four events. The data are obtained from the OMNI 1-minute resolution dataset (http://omniweb.gsfc.nasa.gov). For each event, four days of data are shown. The January 2007 event is represented by the black line in each panel. The CIR, marked by the shaded area, is characterized by the gradual increase of the solar wind speed from around 300 km/s to about 700 km (panel (a)) on January 29th, the increase of the solar wind density just before the speed increase (panel (b)), the increases of the solar wind temperature (panel (c)) and the interplanetary magnetic field (IMF) magnitude (panel (d)), and the highly varying IMF B_z (panel (e)) on January 29th. The SYM-H index drops to about −50 nT, and the AE index reaches about 1500 nT on the same day.

Fig. 1. Solar wind conditions and geomagnetic activity indices for the four CIR/HSS events.

For the other three events shown in Figure 1, the OMNI data are time-shifted such that the first peaks of the IMF intensity within the CIRs line up with the one of the January 2007 event, marked by the vertical dotted line in panel (d). This allows comparisons between different events. In terms of the overall trend of the solar wind conditions, the three CIR/HSS events in 2011 and 2012 are very similar to the 2007 event. However, comparing to the 2007 event, notable differences include a broader CIR in the May 2012 event, a less intensified IMF across the CIR and weaker geomagnetic responses in the June 2012 event, and less intensified solar wind densities ahead of the CIRs in the April 2011 and June 2012 events.

For each of the HSS events, we also choose several quiet days just before event to establish quiet states of the ionosphere for GITM simulations (Sect. 2.2) and the GIM data (Sect. 2.3). Our criterion for a quiet day is based on the daily Ap index. We selected six or seven consecutive days with the Ap index less than 8 as the quiet days for each event. The specific quiet days and storm days that we study are given in Table 1.

2.2. GITM simulations

The TEC forecasts are made with GITM, a widely used first-principle model of the global ionosphere and thermosphere. GITM is governed by three-dimensional (3D) transport equations of magnetized plasmas (Ridley et al. 2006), and its computational grid structure is based on geographic longitude, latitude, and altitude. GITM has a rich variety of options in specifying the solar irradiance, high-latitude electric fields, and auroral precipitation. This flexibility encourages us to configure simulations in a forecasting mode. In a forecasting-mode simulation, the solar irradiance is determined by the F10.7 flux, the high-latitude electric field is obtained from the Weimer 2005 empirical model (Weimer 2005a, 2005b) that is driven by solar wind conditions, and the auroral precipitation is obtained from the OVATION (Oval Variation, Assessment, Tracking, Intensity, and Online Nowcasting) Prime model (Newell et al. 2009) that is also driven by the solar wind. Both the Weimer 2005 and OVATION Prime models are embedded within GITM. To summarize, the input parameters of a forecasting-mode GITM simulation only include the F10.7 flux and solar wind conditions, both considered to be predictable. The F10.7 radio flux has been routinely forecasted at the Space Weather Prediction Center (SWPC) (http://www.swpc.noaa.gov/forecasts). The solar wind conditions can be predicted with physics-based heliospheric models (Odstrcil 2003; Riley et al. 2013; van der Holst et al. 2014). Moreover, in the forecasting mode, none of the model parameters are tuned specifically for a simulation.

For each CIR/HSS event, we perform a forecasting-mode GITM simulation, starting from the first quiet day until the end of the storm days in Table 1. A complete forecast requires predicted F10.7 flux and solar wind to drive GITM. However, these predictions would bring in additional forecasting uncertainties that are difficult to separate from the uncertainties of GITM itself. This is beyond the scope of the paper. Since we would like to evaluate the performance of GITM solely, we drive GITM with measured F10.7 flux and solar wind. The simulation input includes the daily F10.7 flux from the Space Physics Interactive Data Resource (http://spidr.ngdc.noaa.gov/spidr/home.do) and the time-varying solar wind from the OMNI 1-minute resolution dataset. The simulation is initialized with neutral atmospheric states from MSIS (Mass Spectrometer and Incoherent Scatter) (Hedin 1991) and HWM (Horizontal Wind Model) (Drob et al. 2008), as well as ionospheric states from IRI (International Reference Ionosphere) (Rawer et al. 1978).

We extract global distributions of TECs every hour from the simulations. The spatial resolution is 3.33° in longitude and 1° in latitude. Since GITM covers from 100 km to about 600 km in altitude, the TECs are produced from vertically integrating the electron density between 100 km and 600 km altitudes. Since the TEC contributions below 100 km altitude and above 600 km altitude are not modeled by GITM, the simulated TECs should be constantly lower than the vertical TECs in reality. This is partially accounted for in our TEC metric 2.4, such that the simulated TECs may be used to forecast TEC variations.

2.3. Global ionospheric maps (GIM)

We utilize the ionospheric TEC data provided by the JPL GIM (Mannucci et al. 1998; Iijima et al. 1999) as reference TECs to validate the forecasted TECs. GIM are essentially global TEC maps retrieved from GPS-based measurements. The measurements are collected from a global GPS receiver network and interpolated in space and time. The spatial resolution of the TEC maps is uniformly 1° in both longitude and latitude. We obtain the TEC maps for all quiet and storm days specified in Table 1. These TEC maps are taken from the average TECs during the first 15 min of each hour within the entire time interval. This is somewhat different from the GITM TEC maps taken at every hour, but we can still make comparisons between the TEC maps from GITM and GIM, if we assume that the TEC changes within 15 min are small compared to the TEC changes on an hourly basis during geomagnetic storms Afraimovich et al. 2013.

Figure 2 shows examples of TEC maps from both GITM and GIM. The TEC maps are for two successive hours during the January 2007 event. The overall TEC patterns from GITM and GIM are similar. However, the GIM TEC is significantly larger than the GITM TEC. This could be partially due to the fact that GITM does not model the atmosphere below 100 km and above 600 km altitudes. Besides the difference in the TEC magnitude, the GITM and GIM TECs are peaked at different local times. The GITM TEC is peaked around local noon: near 0° longitude at 12:00 UT and shifted to the west at 13:00 UT, while the GIM TEC is peaked at post-noon.

Fig. 2. TEC maps from the GITM simulation (left) and GIM data (right) at two successive hours for the January 2007 event.

Although GIM TECs are good representations of the actual TECs, the accuracy of GIM depends on the distribution of GPS receivers. For certain areas like oceans, where GPS receivers are few or even absent, the GIM-provided TECs are less accurate due to a local time interpolation and smoothing method. Therefore, it is worthwhile to consider GPS receiver distribution when evaluating the GITM forecasts with GIM.

2.4. TEC metric

To quantitatively characterize ionospheric TEC responses to the CIR/HSS events, we develop a TEC metric that describes TEC changes normalized to quiet-time TEC variations. In this way the metric partially compensates for the different altitude regions covered by GITM and GIM, and permits GITM forecasts to be informative regarding TEC effects. We apply this TEC metric to TECs from both the GITM forecasts and GIM data, in order to make comparisons between the two.

We divide each of the hourly TEC maps into geographically registered grid boxes of size 30° (longitude) × 15° (latitude). This resolution is much coarser than the resolutions of the original GITM and GIM TEC maps, but it still reveals the longitudinal and latitudinal dependence of TEC variations, and is useful for model evaluation and analysis. The total number of grid boxes for a coarsened TEC map is 144. For a given universal time (UT) on a storm or quiet day, we define the TEC perturbation of the grid box centered at [lon, lat] as(1)where TEC(UT, lon, lat) is the mean TEC in this grid box at UT. is the N-day-averaged quiet-time TEC for the same grid box at UT, written as(2)σ_TEC,quiet(UT, lon, lat) is the standard deviation of the quiet-time TEC:(3) N is the number of quiet days. For January 2007, April 2011, and May 2012 events, N = 7. For June 2012 event, N = 6. σ_TEC,quiet(UT, lon, lat) location. Typical values are between 0.5 TECU and 4 TECU for GIM, and between 0.3 TECU and 2 TECU for GITM.

The TEC metric, expressed by Eq. (1), is analogous to a signal-to-noise ratio. We compute the metric for all quiet and storm days during the four events. The storm-time dTEC is expected to have a larger amplitude than the quiet-time value. Again, for GITM we have hourly dTEC snapshots, while for GIM we have dTEC based on average TECs during the first 15 min of each hour.

The metric is fundamental to our forecasting study, because it is meant to reveal TEC variations that are due to storm-time physics in particular. We note that this metric is specifically constructed to elucidate physical conditions during storms versus quiet time. We make the distinction here between a user-oriented metric and a metric meant for scientific analysis of the forecast. For the latter, variations that exceed quiet-time variability are key to identify storm-time physics.

3. Results

Taking the 2011 event as an example, we visualize the simulated dTEC over the global map in Figure 3. dTEC maps computed from the GIM data are also shown for comparison. Four different times before and during the storm are shown in four rows. The grid boxes are displayed on individual maps, bordered by dotted lines. The value of dTEC for each grid box is written in the box. To identify TEC perturbations during the storm, we define five levels for dTEC:

1. dTEC ≤ −8.0

2. −8.0 ≤ dTEC ≤ −4.0

3. −4.0 ≤ dTEC ≤ 4.0

4. 4.0 ≤ dTEC ≤ 8.0

5. dTEC ≥ 8.0

Fig. 3. dTEC maps from the GITM simulation (left) and GIM data (right) for four different times during the April 2011 event.

These levels are represented by purple, blue, green, yellow, and red, respectively. A grid box is filled with a certain color corresponding to the dTEC value in this grid box. The entire map is thus color-coded by the TEC disturbance level. In this way the regions with severe or moderate TEC perturbations, whether positive or negative, are distinctive on the map.

We examine the GITM-simulated dTEC maps in the left column of Figure 3 first. At 23 UT on 28 April 2011, the CIR has not arrived yet. The ionosphere is considered to be quiet, since all regions are covered by green in both the GITM simulation and GIM data. The dTEC values provided in the grid boxes are mostly around 1.0, indicating that the TEC perturbations are close to the quiet-time variability. When the storm hits the ionosphere at 19 UT on 29 April 2011, the TEC in the southern hemisphere high-latitude region is positively perturbed, as the yellow color shows up in this region. Four hours later, at 23 UT, the ionosphere is highly perturbed in all latitudes. This marks the main phase of the storm. Interestingly, the southern hemisphere seems much more perturbed than the northern hemisphere. Moreover, the TEC perturbations in the southern hemisphere are mainly positive, while the TEC perturbations in the northern hemisphere are mostly negative. A snapshot at 7 UT on 1 May shows the ionospheric response during the storm recovery phase. Most of the regions have returned to a quiet state, except some moderate activity in the middle-low latitudes indicated by the yellow and blue boxes. The simulated dTEC maps show similar trends as the ones from the GIM data, shown in the right column of Figure 3. However, the geographic locations and magnitudes of the TEC perturbations from the GIM data are not well reproduced by the GITM simulation.

We propose that the simulated dTEC maps shown in the left column of Figure 3 can be employed as ionospheric weather forecasting products. Such a dTEC map could be produced every hour, or more frequently if needed, generating a time sequence of dTEC maps. It would be very straightforward to tell when and where a TEC perturbation occurs from these dTEC maps. In operational ionospheric weather forecasting, certain alarms could be raised according to the disturbance level on the dTEC maps.

As seen earlier, the GIM data helps validate and evaluate the GITM forecasts. We make comparisons between dTEC from GITM and GIM for all grid boxes and for all events. Examples from the May 2012 event are displayed in Figure 4. These comparisons are for the grid boxes in the longitudinal bin 120° E–150° E. This particular longitudinal bin is selected because of the excellent coverage of GPS receivers in this bin, which ensures the quality of the GIM data. From top to bottom, the panels in Figure 4 show dTEC in the grid boxes distributed from northern to southern hemispheres. The latitudinal ranges are labeled on the panels. All panels share the same horizontal axis representing the time, which covers the quiet days from May 1 to May 7 and the storm days from May 8 to May 10. During the quiet days, dTEC values from both the GITM forecast and the GIM data are close to zero with very small oscillations, though the agreement between GITM and GIM results is not perfect. During the storm days, large magnitudes of dTEC appear in both the forecast and the data. At some latitudes, like 30° S–45° S and 45° S–60° S, the forecasted peaks and dips of dTEC match with the ones in the data, especially on May 9th. However, the forecast overestimates or underestimates the storm-time dTEC at some other latitudes and other times. For example, in 0°–15° S and 15° S–30° S, the forecast significantly overestimates the dTEC values on May 9; in 15° N–0°, the forecast clearly fails to predict the dTEC dip on May 8, while in 75° S–90° S, the forecast misses the dTEC peaks on May 8.

Fig. 4. The GITM dTEC (black) and GIM dTEC (blue) as functions of time for each grid box within the longitudinal bin 120° E–150° E during the May 2012 event.

Below we statistically evaluate the performance of the dTEC forecasts for all grid boxes and for all events. We set a threshold for the magnitude of dTEC, i.e., |dTEC|, as |dTEC|_threshold. For a given grid box and at a given time, we define the following categories: if both the GIM and GITM dTEC magnitudes reach |dTEC|_threshold, we have a “hit”; if the GIM |dTEC| reaches |dTEC|_threshold but GITM |dTEC| does not, we have a “miss”; if the GIM |dTEC| is below |dTEC|_threshold but GITM |dTEC| reaches, we have a “false alarm”. We do not consider the sign of dTEC here for simplicity, as we are more concerned about the magnitude of a TEC perturbation than its sign at the present time. We think this would be also true for operational ionospheric forecasting. Table 2 summarizes the amount of hits, misses, and false alarms for the storms days. Two different thresholds, 4.0 and 8.0, are applied for each event. According to the table, the forecasting performance is merely acceptable. Although several hits are achieved, more misses and false alarms are found for both thresholds in all events.

To find out the forecasting performance for various disturbance levels, we plot Hits/(Hits + Misses) as a function of |dTEC|_threshold in Figure 5. The total number of hits and misses is essentially the amount of “hits” in the GIM data, i.e., how many times that the GIM |dTEC| reaches the threshold. Thus, Hits/(Hits + Misses) is the ratio between the GITM hits and GIM “hits”. This ratio is a measure of how well the forecasts reproduce the TEC perturbations in the GIM data (therefore false alarms are not taken into account). Hits/(Hits + Misses) is always between 0 and 1. The larger the value, the better the prediction is. For a perfect prediction, Hits/(Hits + Misses) equals 1. In Figure 5, Hits/(Hits + Misses) for eight different thresholds: 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, and 10.0 are shown. The four CIR/HSS events are represented by different colors as explained in the legend. The figure reveals several interesting features. First of all, Hits/(Hits + Misses) is below 0.5 for all cases, which implies that less than half of the GIM TEC disturbances are reproduced by the GITM forecasts. Second, a common trend is that the prediction gets worse with an increasing |dTEC|_threshold. In other words, larger TEC disturbances during the storms are not so well captured by the GITM forecasts as are moderate and small TEC disturbances. An exception is the January 2007 event, for which the prediction slightly improves when |dTEC|_threshold increases from 7.0 to 10.0. Third, while the prediction performances are similar for the January 2007, May 2012, and June 2012 events, the performance for the April 2011 event drops well below the others. When |dTEC|_threshold ≥ 8.0, there are even no hits for the April 2011 event. The reason for the worse prediction of hits is not clear. However, we do notice that for a given threshold, the amount of false alarms for the April 2011 event is remarkably less than the one for the other three events. This distinction is apparent in the last column of Table 2. Hence, a more complete evaluation of forecasting performance should also include false alarms.

Fig. 5. The performance of the GITM forecasts for various |dTEC| thresholds.

The performance of the GITM forecasts depends on a variety of factors including: the accuracy of the solar wind conditions measured at L1 point, the accuracy of the solar irradiance derived from the F10.7 proxy, the accuracies of empirical models Weimer 2005 and OVATION Prime, and the sufficiency of the coupled ionosphere-thermosphere physics in GITM. All these factors can affect the forecasting performance and thus are important to consider when interpreting the forecasting performance.

4. Discussion

The validation of the GITM forecasts with the GIM data presented in the previous section indicates both model-data agreements and disagreements. To understand model behavior, we look into model calculations of TECs within individual grid boxes. Here we analyze two neighboring middle-latitudinal grid boxes during the initial and main phases of the January 2007 event. This is a US region of good coverage by GPS receivers measuring TEC. The specific time interval we choose is 13 UT–20 UT on January 29. During the time and within one of the grid boxes, [90° W–60° W, 30° N–45° N], the forecasted dTEC agrees with the GIM data well, while within the other grid box to the west, [120° W–90° W, 30° N–45° N], the forecast greatly underestimates the dTEC. More details about the two grid boxes can be found in Figure 6. Panels (a) and (e) show the GIM and GITM dTEC over the entire event period for the two grid boxes respectively. We mark the time interval of interest with colored vertical lines on both panels. From purple to red, the colors represent every hour between 13 UT and 20 UT on January 29. Since the two grid boxes cover a wide range of longitudes, the local time varies. We take the local time at the center of each grid box as a reference time. For grid box I, [90° W–60° W, 30° N–45° N], the local time at 75° W is from 8 AM to 3 PM. For grid box II, [120° W–90° W, 30° N–45° N], the local time at 105° W is from 5 AM to 12 PM.

Fig. 6. Analysis of two neighboring grid boxes for 13 UT–20 UT on 29 January 2007: panel (a): dTEC from the GIM data (gray) and GITM simulation (black) for the entire event period for grid box I, and the vertical colored lines mark the time interval of interest; panel (b): storm-time (triangles) and 7-quiet-day average TECs (squares) from the GIM data for grid box I; panel (c): storm-time (triangles) and 7-quiet-day average TECs (squares) from the GITM simulation for grid box I; panel (d): simulated hourly electron density profiles for grid box I; panels (e)–(h): same as panels (a)–(d) but for grid box II.

As dTEC incorporates the quiet-time variability of the TEC, it might not sufficiently reflect the absolute TEC changes during the storm time. We plot the absolute GIM TECs in panel (b) and GITM TECs in panel (c) for grid box I. Similarly, the GIM TECs and GITM TECs for grid box II are in panels (f) and (g), respectively. These plots not only show the TECs during 13 UT to 20 UT on January 29, marked by the triangles, but also show the TECs averaged over the seven quiet days (January 22 to January 28) for the same UT time interval, marked by the squares. The colors correspond to various UT times, from 13 UT to 20 UT, as well as various local times as noted just below panels (a) and (e). In panels (b) and (c), although the values of TECs from GIM and GITM are different, both GIM and GITM TECs are elevated significantly during the storm time in grid box I. For the GIM data, the storm-time TEC reaches its peak at 17 UT, which is about 10 TECU larger than the quiet-time average at 17 UT. For the GITM simulation, the storm-time TEC peaks at 16 UT, which is about 5 TECU larger than the quiet-time value. For grid box II, panels (f) and (g) show that the storm-time GIM TEC is elevated by about 7 TECU at its peak (18 UT) from the quiet-time average, while the storm-time GITM TEC has a very little increase from the quiet-time average. Therefore, in grid box II, the GITM simulation does not reproduce an adequate TEC enhancement during 13 UT–20 UT on January 29 seen in the GIM data. As the TECs are height-integrated electron densities, we also plot electron density profiles from the GITM simulation in panel (d) for grid box I and panel (h) for grid box II. There are eight colored lines in each panel, corresponding to hourly electron density profiles between 13 UT and 20 UT on January 29. These are electron densities averaged over longitudes and latitudes within grid box I or II. January 2007 (average F10.7 flux 80 sfu) is at a solar minimum, as a result the ionospheric electron density profiles only exhibit the F-layer instead of F1 and F2 layers during a solar maximum. For grid box I, the variation of the F-layer peak density and height follows closely with the TEC change. At 16 UT when the TEC peaks, the F-layer peak reaches the maximum height. The local time is 11 AM at 75° W. For grid box II, the variation of the electron density profiles from 13 UT to 16 UT is due to the diurnal change. From 16 UT to 19 UT, the electron density increases, but the F-layer peak stays at about the same height or even drops. Comparing the electron density profiles from the two grid boxes at the same local times, for example at 11 AM, the maximum electron density in grid box I is almost 20% larger than the maximum electron density in grid box II. In summary, in terms of the simulated electron density, GITM behaves very differently for the two neighboring grid boxes.

From the continuity equations that GITM solves, we identify two primary contributors to the electron density. One is the transport term, which is the spatial derivative of the momentum. The other is the source term, which includes the production and loss of electrons due to chemical reactions. Our following analysis focuses on the F layer, since it supplies the most electron density to the TEC. We consider the electron density change rate at the F-layer peak heights found from the electron density profiles of grid boxes I and II. To address the electron density change rate, the transport term and the source term due to chemical reactions are evaluated separately. For the transport term, we take the velocity of O⁺ as the electron velocity, because GITM does not directly compute the electron velocity, and O⁺ is the dominate ion species in the F-layer. The electron density change rate due to the east-west transport is calculated as the difference between the electron momentum fluxes at the east and west boundaries of the grid box, divided by the distance between the two boundaries. The calculation of the electron density change rate due to the north-south transport is the same, except replacing the east and west boundaries by the north and south boundaries. Similarly, the vertical transport contribution is calculated. For the chemical reactions, the major electron production mechanism in the F-layer is the photoionization of O, and the electrons are lost through dissociative recombinations with and NO⁺, preceded by the atom-ion interchange of O⁺ + O₂ and O⁺ + N₂ (Prölss 1995; Schunk & Nagy 2000). The photoionization coefficient of O depends on the solar irradiance. In GITM, this coefficient is computed through an EUV flux model. We obtain the time-varying photoionzation coefficient from the GITM simulation. The electron loss rate is controlled by the atom-ion interchange reactions. These reaction rates are temperature-dependent (Schunk & Nagy 2000) and are extracted from the GITM simulation.

Figure 7 illustrates the electron density and its change rate at F-layer peak altitudes for the two grid boxes. Each UT (or local) time is represented by a color in the same way as in Figure 6 for easy identification. The electron density at the F-layer peak height contributes a significant portion to the TEC. Therefore, the electron density variations, shown in the top row of Figure 7, are very similar to the TEC variations represented by the triangles in panels (c) and (g) of Figure 6. The electron density change rates, from the second to the last rows in Figure 7, provide insights of the physics behind the simulation. For grid box I, the electron density increase from 13 UT to 16 UT is due to the transport, or more specifically, the transport in the east-west and vertical directions. The contributions of these two directions are comparable. After 16 UT, the chemical reactions dominate over the transport and reduce the electron density. This reduction is mostly caused by the enhanced electron loss through the O⁺ + O₂ process. For grid box II, the electron density increase from 13 UT to 17 UT is due to the transport, similar to the grid box I case. However, the north-south transport contributes positively, while the east-west transport contributes negatively to the electron density in grid box II, which is the opposite to the situation in grid box I. The behavior of the east-west transport is expected, since the two grid boxes are neighbors in the east-west direction. The vertical transport also plays a role in increasing the electron density during 13 UT and 17 UT, yet its contribution is less than the contribution of the north-south transport after 15 UT. After 17 UT, the positive contribution from the transport roughly balances with the negative contribution from the chemical reactions, so that the electron density variation is small. The chemical reactions are dominated by the O⁺ + O₂ process after 16 UT. According to the GIM TEC data (Fig. 6f), the storm-time electron density in grid box II should start to increase right after 14 UT. The simulated storm-time electron density increase, however, is insufficient especially during 14 UT and 16 UT according to Figure 6g. In comparison, during 14 UT and 16 UT, the storm-time electron density increase in grid box I is reproduced by the simulation (Figs. 6b and 6c). Examining contributors to the electron density changes between 14 UT and 16 UT from Figure 7, we notice that the rate of the electron loss due to the O⁺ + O₂ process increases during the time interval and substantially exceeds the electron production rate at 16 UT in grid box II. On the contrary, for grid box I the electron loss due to the O⁺ + O₂ process is at about the same rates as the electron production during the time. Hence, we suspect that the insufficient storm-time TEC enhancement in grid box II could be attributed to the electron density loss through the O⁺ + N₂ process.

Fig. 7. GITM-simulated electron density and its change rate at F-layer peak height for grid box I (left column) and grid box II (right column) during 13 UT–20 UT on 29 January 2007.

Our analysis implies that the transport process dominates over chemical reactions in increasing the storm-time TEC in grid box I. The lack of the storm-time TEC enhancement in grid box II might be due to the electron loss via the intensified reaction between O⁺ and O₂ and the resultant dissociative recombination process. This may be caused by composition changes that increase the abundance of O₂ and/or temperature changes that accelerate the reaction. Further investigation is required to understand why the O⁺ + O₂ reaction is preferred in the simulation.

5. Conclusions

We have conducted an ionospheric forecasting study with a first-principle model of the global ionosphere and thermosphere, GITM. We apply GITM to simulate the ionospheric TEC responses to four typical CIR/HSS events in a forecasting mode. We quantify the forecasted TECs through a well-defined metric, which gives the TEC perturbation normalized to the quiet-time TEC variability within mesoscale-sized grid boxes, represented by dTEC. We show that dTEC can be visualized on the global map and can be potentially utilized as a product for operational ionospheric weather forecasting.

Comparing the forecasting results with the GIM data based on GPS observations, we find both agreement and disagreement. According to the magnitudes of GITM-forecasted dTEC and GIM dTEC, we categorize the forecasts into three groups: hits, misses, and false alarms. The statistical results from all grid boxes and four CIR/HSS events suggest that: the GITM forecasts successfully reproduce less than half of the storm-time TEC perturbations seen in the GIM data; the forecasts capture smaller TEC disturbances better than larger TEC disturbances during the storms; the forecasts also generate a large amount of TEC disturbances that do not exist in the GIM data; for unknown reasons, the forecasting performance for the April 2011 event is different from the other events.

Understanding the forecasting performance requires a detailed analysis of the TEC computation in GITM. We present such an analysis for two neighboring grid boxes during the January 2007 event. For the specific 8-hour time interval we chose, the simulation captures the storm-time TEC enhancement in one of the grid boxes well, yet underestimates the TEC enhancement in the other grid box. By tracing down the electron density calculation in GITM, we find that the simulated TEC enhancement in the “agreement” grid box can be attributed to the ion transport process in the east-west direction and the chemical production/loss of electrons. The underestimation of storm-time TEC in the other grid box is likely due to the enhanced electron loss through the O⁺ + N₂ reaction and the subsequent dissociative recombination of electrons with NO⁺.

In this paper we present our initial effort on formulating ionospheric forecasting with coupled ionosphere-thermosphere physics-based models. To achieve the final goal of long-term operational ionospheric weather forecasting, we need to further explore model capabilities and interpret simulation results. Achieving this goal also requires improved physical understanding of the storm-time ionospheric-thermospheric response. Future work may include but is not limited to: analyzing more individual grid boxes in detail from the existing CIR/HSS event simulations, extending the forecasting study to Coronal Mass Ejection (CME)-driven ionospheric storms, applying data-assimilative models for validation and insight, and evaluating state-of-the-art models of the ionosphere other than GITM.

Acknowledgments

This work was done at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. The authors would like to acknowledge Prof. Aaron Ridley at the University of Michigan for providing help in setting up GITM simulations. The authors also thank J.H. King, N. Papatashvilli at AdnetSystems, NASA GSFC and CDAWeb for providing the OMNI data (http://cdaweb.gsfc.nasa.gov/istp_public/). The computing resources were provided by the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center. The editor thanks two anonymous referees for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1. Solar wind conditions and geomagnetic activity indices for the four CIR/HSS events.
In the text

	Fig. 2. TEC maps from the GITM simulation (left) and GIM data (right) at two successive hours for the January 2007 event.
In the text

	Fig. 3. dTEC maps from the GITM simulation (left) and GIM data (right) for four different times during the April 2011 event.
In the text

	Fig. 4. The GITM dTEC (black) and GIM dTEC (blue) as functions of time for each grid box within the longitudinal bin 120° E–150° E during the May 2012 event.
In the text

	Fig. 5. The performance of the GITM forecasts for various |dTEC| thresholds.
In the text

Fig. 6. Analysis of two neighboring grid boxes for 13 UT–20 UT on 29 January 2007: panel (a): dTEC from the GIM data (gray) and GITM simulation (black) for the entire event period for grid box I, and the vertical colored lines mark the time interval of interest; panel (b): storm-time (triangles) and 7-quiet-day average TECs (squares) from the GIM data for grid box I; panel (c): storm-time (triangles) and 7-quiet-day average TECs (squares) from the GITM simulation for grid box I; panel (d): simulated hourly electron density profiles for grid box I; panels (e)–(h): same as panels (a)–(d) but for grid box II.

In the text

	Fig. 7. GITM-simulated electron density and its change rate at F-layer peak height for grid box I (left column) and grid box II (right column) during 13 UT–20 UT on 29 January 2007.
In the text