Capturing the Storm-Time F-Region Ionospheric Response

in an Empirical Model

by

T.J. Fuller-Rowell, M.V. Codrescu, and E.A. Araujo-Pradere

NOAA, Space Environment Center and CIRES, University of Colorado, Boulder, Colorado 80305, USA



Abstract

Understanding the response of the thermosphere and ionosphere to geomagnetic disturbances has reached a point where some of the F-region ionospheric characteristics can be captured in an empirical model. The integrated effect of magnetospheric energy injection during a storm disturbs the thermospheric temperature, global circulation, and neutral composition.  The composition modifications accumulate over tens of hours, have a particular seasonal/latitude structure, and are very slow to recover. The long-lived reaction of the neutral upper atmosphere produces an ionospheric response that is consistent from storm to storm, and so enables the characteristics to be captured by an empirical model. The goal of this first empirical storm-time model is to establish a correction to the F-region peak density, or critical frequency, as a function of season and latitude for any given ap time history of a storm. Guided by the knowledge gained from previous data analysis, and from simulations with a physically-based model, observations of the F-region peak density from available sites and from many storms have been sorted by latitude and season, and by the magnitude of the storm. A coherent picture begins to emerge particularly in the summer and equinox mid-latitudes. Several features are still unable to be included in the empirical model, although they are clearly important, and can be simulated in physical models. These include the local-time dependence, and the dynamic response to transient large-scale gravity waves. The latter in particular requires accurate knowledge of the spatial and temporal variation of the geomagnetic sources.



Introduction


The maturity of an area of science can often be reflected in the ability to capture the physics within an empirical model. The knowledge of the physics enables the appropriate choice of sorting parameters. As this maturity grows so does the accuracy with which the empirical models are able to reproduce observations. For the thermosphere, the neutral component of the upper atmosphere, the Mass Spectrometer and Incoherent Scatter (MSIS) empirical model (Hedin et al., 1987) is the undisputed standard for specifying neutral temperature and composition around the globe as a function of season, latitude, longitude, Universal Time, solar and geomagnetic activity. The model has evolved from the early versions in the 1970s to the current version of the 1990s, by including more data and by incorporating new physical understanding of the controlling parameters. The International Reference Ionosphere (IRI) can be regarded as the equivalent standard for the ionosphere (Bilitza, 1990). The global distribution of ionospheric parameters can be specified as a function of a similar set of sorting variables. The one exception is that, currently, the IRI has no geomagnetic activity dependence. The reason is simple; the characteristics of the ionospheric response to storms had not been understood to a level where the development of an empirical model was feasible.

Physical understanding of the ionospheric response to storms has now sufficiently matured to enable the first rudimentary characterization of the storm-time response of the peak F-region ionospheric density to be made as a function of season and latitude. This step has been made possible by a combination of analyses of extensive networks of ionospheric data and detailed simulations using coupled thermosphere ionosphere models.

Geomagnetic storms result when high-speed plasma injected into the solar wind from coronal mass ejections or coronal holes impinges upon Earth’s geomagnetic field. If the arriving solar wind plasma has a southward magnetic field energy is coupled efficiently into Earth’s magnetosphere and upper atmosphere. The magnitude of the ensuing geomagnetic storm has come to be defined by the strength of the low latitude magnetic index, Dst, which is a measure of the magnetospheric ring current. Although the ring current is not the main driver of the upper atmosphere, the sources for both are related; Dst can reflect the level of magnetospheric energy input to the upper atmosphere.

A geomagnetic storm, from the perspective of the upper atmosphere, is a period of intense energy input from the magnetosphere for a period of several hours to days. The manifestations of a storm are well defined and include the following effects. Auroral electron precipitation increases in magnitude and expands to lower latitudes than normal, and these particles heat and ionize the gas and increase the conductivity of the atmosphere. The magnetospheric electric field mapped to the atmosphere intensifies and expands in concert with the aurora, and combines with the increased conductivity to produce large enhancements of Joule heating, which is the dominant atmospheric energy source during a storm. Joule heating can increase from tens of gigawatts during quiet times, to hundreds of gigawatts during severe geomagnetic disturbances. The combined input can dump thousands of Terajoules of energy into the upper atmosphere during the course of a storm, and can raise the temperature of the gas by hundreds of degrees Kelvin. Thermal expansion of the atmosphere raises neutral density and can have significant effects on satellite drag. Ionospheric ions drift in response to the electric field and, by colliding with the atmosphere, can drive winds in excess of 1 km/s at high latitudes. Many interesting science problems center around these high latitudes processes, and the combined energy source is the driver of the global storm-time changes in the upper atmosphere.

The energy injection drives global thermospheric storm winds and composition changes, and many of the ionospheric changes at midlatitude can be understood as a response to these thermospheric perturbations. Wind surges propagate from high latitudes around the globe and transport plasma along magnetic field lines to regions of altered chemical composition, changing ion recombination rates. Many of the increases in NmF2 and TEC are thought to result from these “traveling atmospheric disturbances.” The divergent nature of the wind causes upwelling of molecular rich thermospheric gas from lower altitudes. These regions of enhanced molecular species at F region altitudes, or composition “bulges,” can be transported by the pre-existing background wind fields and by the storm winds [Fuller-Rowell et al., 1996]; the changed composition again feeds back to the ionosphere. The regions of upwelling (increases of N2 and O2) cause the ionosphere to decay faster and create negative phases of ionospheric storms.

The dynamical changes are complicated because they are driven by the highly variable magnetospheric energy sources. The wind surges propagate and interact around the globe and often appear as a random mixture of waves. Exactly where a composition bulge will be created is also difficult to determine; composition changes are created by persistent divergence of the wind field, in areas of significant energy injection (mainly Joule heating). Accurate knowledge of the spatial and temporal distribution of the magnetospheric sources is required to predict where and when these composition changes will manifest themselves.

Much of the interest in understanding the response of the upper atmosphere to geomagnetic storms has stemmed from the need to predict the ionospheric response. The need arises for practical reasons; the requirement for ground-to-ground communication via the ionosphere using HF radio propagation and from ground-to-satellite through the ionosphere at higher frequencies. The parameters that have received then most attention are the peak F region electron density (NmF2), or critical frequency (foF2), which is related to the maximum usable frequency (MUF) for oblique propagation of radio waves, and the total electron content (TEC), which is significant for the phase delay of high frequency ground-to-satellite navigation signals.



Physical Understanding from Data and Models


In spite of the complexities in the observed response of the thermosphere and ionosphere to a geomagnetic storm, systematic features are apparent. The breakthrough in understanding of the storm-time ionosphere has come from analysis of extensive ionospheric observations and from interpretation of the data by physical models. Rodger et al. (1989) showed that the response of the ionosphere reveals both seasonal and local-time (LT) dependencies during a geomagnetic storm. They demonstrated that at a southern magnetic mid-latitude station, a consistent local time signature in the ratio of disturbed to quiet NmF2 existed throughout the year, with a minimum in the morning hours around 06 LT and a maximum in the evening hours around 18 LT (see Figure 1).The local-time "AC" variation was superimposed on a "DC" shift of the mean level that varied with season, being most positive in winter (May-July) and most negative in summer (October-February). The data supported the widely held belief that “positive phases” of storms (increases in electron density) are more likely in winter mid-latitudes, and “negative phases” of storms (decreases in electron density) are more likely in summer. Field and Rishbeth (1997) showed that these same characteristics are true for other longitude sectors. Rodger et al (1989) stressed the point that individual storms show large deviations from the average behavior.

A cause of the LT variation was suggested by Fuller-Rowell et al. (1994), who extended the theory of Prölss (1993). Prölss suggested that negative storm effects are due to regions in which the neutral gas composition is changed — the ratio of molecular gas concentration (N2 + O2) to the atomic oxygen concentration is increased. Such a region, which we call a "composition bulge" because it represents a region of increased mean mass, is originally produced through heating and upwelling of air by the magnetospheric energy inputs at auroral latitudes. The bulge can be moved to middle latitudes by the nightside equatorward winds and brought onto the dayside as Earth rotates.

Fuller-Rowell et al. (1994) made computational simulations of this process for storms at equinox. They showed that once the composition bulge is created it can be transported by a wind field, either that of the background quiet-day thermospheric circulation, or of the storm circulation driven by the high latitude heat input. They attributed the local time AC effect derived by Rodger et al. (1989) to an oscillation in latitude of the composition bulge in response to the diurnally varying winds. Skoblin and Förster (1993) also showed a case where steep gradients in thermospheric composition could be advected by the meridional wind.

The work of Fuller-Rowell et al. (1996) extended this idea to suggest an explanation of the seasonal variations. Numerical computations suggest that the prevailing summer-to-winter circulation at solstice transports the molecular rich gas to mid- and low-latitudes in the summer hemisphere over the day or two following the storm (see Figure 2a). In the winter hemisphere, poleward winds restrict the equatorward movement of composition (see Figure 2b). The altered neutral-chemical environment in summer subsequently depletes the F-region mid-latitude ionosphere to produce a negative phase. In winter mid-latitudes a decrease in molecular species, associated with downwelling, persists and produces a positive phase. The seasonal migration of the bulge is superimposed on the diurnal oscillation.

This seasonal dependence of the response of the neutral atmosphere composition as a function of season is depicted in Figure 2. The top panel shows summer, the middle is winter conditions, and the lower panel is the equinox case. Each one shows a snap-shot of the northern hemisphere from 10º latitude to the pole, at about 300 km altitude, 24 hours after a 12 hour storm. The figure is from numerical simulations using the CTIM model (Fuller-Rowell et al., 1996). In the upper panel the “bulge” of increased mean molecular mass (equivalent to an increase in the N2/O ratio) has been transported by the wind field to low latitudes. In the middle panel, in the winter solstice, the composition bulge has been constrained to high latitudes. The lower panel for equinox is the intermediate case.

Fuller-Rowell et al. (2000) demonstrated that numerical simulations are able to capture the seasonal variation in the ionospheric response to storms and indicated clear “scientific success”. Figure 3 shows a numerical simulation of a 25-day period during November and December of 1997, for sites in a) the winter hemisphere (Rome) and b) the summer hemisphere (Grahamstown, SA). The top panels of the Figures 3a and 3b show the hourly

F-region peak critical frequency (NmF2), compared with the monthly mean. The ratio of the two, in the middle panel, shows the increase or decrease during the geomagnetic disturbance on days 326 and 327. An increase occurs in the winter hemisphere at Rome (3a), and a decrease occurs in the summer hemisphere at Grahamstown (3b). The results from a numerical simulation using the CTIM model (Fuller-Rowell et al., 1996), where variability of the high-latitude electric field had been included (Codrescu et al., 2000), show good agreement in both cases. The lower panels in Figure 3 show the variation in the geomagnetic index ap during the interval.

An important point to note is that the physical understanding that has emerged does not translate automatically to operational value. In contrast to the “visual” success of the model, Fuller-Rowell et al. (2000) showed that detailed quantitative comparisons of the physical model with data, which are necessary for space weather applications, were less impressive. The accuracy, or value, of the model was quantified by evaluating the daily standard deviation, the root-mean-square error, and the correlation coefficient between the data and model predictions. During the storm periods, the RMS error of the model improved on IRI only slightly, and there were occasionally false-alarms. Using unsmoothed data over the full interval, the correlation coefficients between the model and data were low, between 0.3 and 0.4. Isolating the storm intervals increased the correlation to between 0.43 and 0.56, and by smoothing the data the values increased up to 0.65. The study illustrated the substantial difference between scientific success and a demonstration of value for space weather applications.



Empirical Storm-Time Correction Model


The physical understanding that has emerged from the numerical simulations and the data analysis has set the stage for developing an empirical storm-time ionospheric correction model. Specifically, we would like to be able to correct the monthly median, or climatological value (IRI, for instance), of the F-region density during a geomagnetic storm, as a function of season and latitude. One advantage of the empirical approach is that data is used to establish the correct relationship between the parameters and the ionosphere, providing the correct sorting parameters are selected. The physical model can be used to determine the qualitative relationship, but we do not have to rely on the model to provide the quantitative dependence. The question is, do we have the knowledge and understanding to include a geomagnetic dependence in the International Reference Ionosphere, that demonstrates a real improvement? The goal of this first empirical storm-time model is to establish a correction to the F-region peak density, or critical frequency, as a function of season and latitude for any given ap time history of a storm.

One of the keys to developing a successful algorithm is to determine the relationship between the parameter representing the magnitude of the geomagnetic forcing and the ionospheric response (the function F(t) in equation 1). For example, the upper panel of Figure 4 illustrates the time history of the ap index over a five-day interval in June, 1981. In the lower panel, the evolution of the ionosphere at mid-latitude in the southern hemisphere is shown. The parameter shown is the ratio of the disturbed-to-quiet F region peak critical frequency averaged over a number of sites. The first task in developing an empirical model is to determine the relationship between the driving parameter, ap, and the ionospheric response.

The knowledge gained from the past data analysis and numerical simulations suggests an ionospheric storm-time correction algorithm of the form:

F = {a0 + a1 X (t0) + a2 X2 (t0) + a3 X3 (t0) } { 1 + a4 sin ( LT + a ) } (1)

where X(t0) = ò F(t) ap(t0 - t) dt

The target parameter F, is the value used to scale the quiet-time ionospheric F-region critical frequency, as a function of latitude, day-of-year, and local time. X(t0) is a new index representing the integrated effect of geomagnetic activity over the previous day or two. The index is calculated using a filter function, F(t), that weights the time history of ap. which is described further below.“a”, “a1”, “a2”, and “a3” are coefficients for the polynomial fit to the non-linear relationship between the filtered ap and the ionospheric response, which are dependent on magnetic latitude and season, and “a4” and “a are the amplitude and phase of the local-time dependence.

The optimum shape and length of the filter F(t) defining the relationship between the time history of ap and the ionospheric response at mid-latitude (the impulse response function) was determined by a linear regression technique (Detman and Vassiliadis, 1997). Ideally, the function should be derived for each latitude and season, but with the limited data volume this was not possible. Data from summer mid-latitudes were used to develop the filter, by minimizing the mean square error between the filter input (ap) and filter output (the ionospheric ratios). Figure 5 shows the shape of the unsmoothed filter (dashed line), together with a smoothed profile (solid line) that was adopted for the algorithm. During the early time history of ap, 0 to 7 hours into the past, the filter shape is quite variable probably due to the ionospheric effects of gravity wave propagation and penetration electric fields. The predominant negative response is due to the development of the composition changes, as described earlier, and the slow recovery over the next day or so, reflects the time-scale for molecular diffusion that gradually returns the neutral atmosphere to its pre-storm configuration.

The theory described previously suggests that the ionospheric response will depend strongly on both season and latitude. Using the new index X(t0), ionosonde data from 75 ground-based stations covering latitudes from 79°S to 83°N for a number of storms over the last twenty years have been sorted as a function of latitude and season (see Araujo-Pradere et al. 2001 and 2000, for further details). The results of this analysis are shown in Figure 6. Each panel shows a scatter plot of the storm-to-quiet NmF2 ratio as a function of the filtered ap index, X(t0). This has been done for four magnetic latitude bands 0-20°, 20-40°, 40-60°, and 60-80°, and five seasonal bands centered on the summer and winter solstices, equinox, and the two intermediate times. The data show a consistent negative response in summer mid-latitudes, while in the winter hemisphere the response is not so well defined, showing a boundary around 400. The consistent response in summer is due to the prevailing summer-to-winter circulation. In the winter hemisphere, theory suggested a boundary exists in the prevailing circulation and in the composition response. Such a boundary also exists in the sorted data producing a negative phase in latitudes greater than 400, while in lowest latitudes a decrease in molecular species, associated with downwelling, persists and produce the characteristic positive storm.

Another important difference between summer and winter hemispheres is the variability in both sets of data. Summer hemisphere and equinox mid latitudes shown a very coherent behavior, with a narrow band of variability around the fit following the negative phase, while winter hemisphere shows a highest dispersion around the fit. In each panel a polynomial cubic fit to the data has been determined to provide the set of coefficients “a0”, “a1”, “a2”, and “a3” required in Equation 1.



Testing the Empirical Model


The empirical storm-time correction model has been tested on a twenty-five day interval towards the end of 1997, between November 12 and December 6. During this period, a significant disturbance occurred on November 22/23. Figure 7 illustrates the response of the ionosphere to the event and the prediction of the empirical model at two sites, Rome in the northern winter midlatitudes (top panel), and Grahamstown, SA, in the southern summer midlatitudes (lower panel). The filtered ap index, using the function in Figure 5, shows a large increase on the two storm days. Both panels show the ionospheric response depicted as the storm-to-quiet ratio of foF2. At Rome, the ionospheric response is positive, consistent with expectations in winter midlatitudes. At Grahamstown, the ionospheric F region decreases, again consistent with expectations in summer midlatitudes. In both cases, the empirical model captures the direction of the change, and the magnitude is particularly good in summer.



Outstanding Issues


At this stage in the development of the empirical algorithm, we have not included the local time dependence represented by coefficient a4 in Equation 1. The analysis by Rodger et al. (1989) showed a strong local time signature with a variation of about 40% in NmF2, but we have been unable to show such a strong dependence in the present analysis. The results of our analysis are shown in Figure 8. The figure shows the residuals after performing the cubic fit to the summer midlatitude data in Figure 6, as a function of local time. The phase of the local time variation is in agreement with the Rodger et al. (1989) study and the theory described by Fuller-Rowell et al. (1996), with the maximum in the negative phase at dawn, and the peak in the positive phase at dusk. The surprising result is the very small amplitude of the local-time response, which does not agree with the previous study.

At the present time it is also not possible to capture the transient effect during the early phase of a geomagnetic storm. The numerical simulation and observations clearly show the propagation of large-scale gravity waves that are known to have a significant influence on the ionosphere at mid latitudes. At present their ionospheric consequences are represented by the structure in the new index based on the time filtered geomagnetic index, represented in Figure 5.

The response at low latitudes is also poorly represented by the empirical model. This is due to our current rudimentary understanding of the temporal evolution of the ionosphere at low latitude during geomagnetic storms. Due to the nature of the magnetic field configuration, the low latitude ionosphere is strongly influenced by electric fields.  Changes in the low-latitude electric field can arise from the penetration of magnetospheric fields and through the development of polarization fields from the changing thermospheric wind fields. The effect on the ionosphere at low latitudes has yet to be understood. When the response of the equatorial region has reached the same level of maturity as at mid-latitudes, then maybe a relatively simple algorithm that captures the response will become apparent.



Conclusion


Our understanding and knowledge of the way the upper atmosphere and ionosphere responds to an increase in geomagnetic activity has matured considerably over the last ten years. An empirical storm-time ionospheric correction model has been developed that shows promise for improving IRI. The algorithm is designed to scale climatological F-region peak critical frequencies, foF2, or electron densities, NmF2, and possibly total electron content, TEC, every hour through a geomagnetic event. This model is to be included in a new version of IRI in an effort to parallel the development of MSIS by including a dependence on geomagnetic activity.

The prospect for improvement of the model in the future is uncertain. The most promising area is to redo the analysis of the local-time variation. The analysis may require the amplitude of the local-time response to be normalized by the mean response in a given latitude and seasonal bin (as given in Figure 6). Work in this direction appears to be showing promise, which will hopefully further improve the fit to the data, and reduce the residual error.

Another avenue to pursue is the use of the state-space reconstruction method, to separate the response characteristics as a function of the phase of the disturbance. This approach may be a way of capturing the affect of the propagating wave features or penetrating electric field, but will have to await further analysis. The index ap was chosen to capture the time-history of the geomagnetic events in order to maximize the use of the model by the community, for operational use, and to enable the correction model to be implemented in the next generation IRI. The limitations of a single index are obvious, but the availability of a simple alternative, suitable for implementation in an empirical model, is not clear. The time-history of solar wind parameters is one possible option to replace ap in the future.



References


Araujo-Pradere, E.A., T.J. Fuller-Rowell, and M.V. Codrescu, STORM: An empirical storm-time ionospheric correction model. Radio Science, submitted, 2001.

Araujo-Pradere, E.A. and T.J. Fuller-Rowell, A model of the perturbed ionosphere using auroral power as the input. Geof. Int., 39, 1, 2000.

Bilitza, D., International reference ionosphere 1990, National Space Science Data Center, WDC A for Rockets and Satellites, Greenbelt, MD, 1990.

Codrescu, M.V., T.J. Fuller-Rowell, J.C. Foster, J.M. Holt, S.J. Cariglia, Electric field variability associated with Millstone Hill electric field model, J. Geophys. Res., 105, 5265-5273, 2000.

Detman, T.R. and D. Vassiliadis, Review of techniques for magnetic storm forecasting, AGU Monograph, 98, 253-266, 1997.

Field, P. and H. Rishbeth, The response of the ionospheric F2-layer to geomagnetic activity: an analysis of worldwide data, J. Atmos. Solar-Terr. Phys., 59, 163-180, 1997.

Fuller-Rowell, T.J., M.V. Codrescu, R.J. Moffett, and S. Quegan, Response of the thermosphere and ionosphere to geomagnetic storms, J. Geophys. Res., 99, 3893-3914, 1994.

Fuller-Rowell, T.J., M.V. Codrescu, R.J. Moffett, and S. Quegan, On the seasonal response of the thermosphere and ionosphere to geomagnetic storms, J. Geophys. Res., 101, 2343-2353 1996.

Fuller-Rowell, T.J., E.A. Araujo-Pradere, and M.V. Codrescu, An empirical ionospheric storm-time correction model. Adv. Space Res., 25, 1, 139-146, 2000.

Hedin, A.E., MSIS-86 thermospheric model. J. Geophys. Res., 92, 4649, 1987.

Prölss, G.W., Common origins of positive ionospheric storms at middle latitudes and the geomagnetic effect at low latitudes, J. Geophys. Res., 98, 5981-5991, 1993.

Rodger, A.S., G.L. Wrenn, and H. Rishbeth, Geomagnetic storms in the Antarctic F region, II, physical interpretation, J. Atmos. Terr. Phys., 51, 851-866 1989.

Skoblin, M.G., and M. Forster, An Alternative Explanation of Ionization Depletions in the Winter Night-time Storm Perturbed F2-Layer, Ann. Geophys., 11, 1026-1032, 1993.