Abstract

[1] Two assimilative models of the ionosphere are described: the Electron Density Assimilative Model (EDAM), developed by QinetiQ, and the IonoNumerics model developed by Fusion Numerics, Inc. Output from each technique has been compared to independent validation data measured by oblique and vertical ionosondes. Results indicate that when tested against vertical ionosondes, both models can reduce RMS errors compared to a median model. For example, the daytime RMS errors in the F region critical frequency above the Eglin Air Force Base ionosonde are 1.2 MHz for the Parameterised Ionospheric Model, 0.7 MHz for EDAM, and 1.0 MHz for IonoNumerics. Testing conducted against an oblique ionosonde has helped to expose a potential problem with ionospheric hmF2 due to poorly specified empirical driver models within IonoNumerics. The testing shows the usefulness of ray-tracing experiments for improving numerical models and the limitations of only testing models against independent total electron content measurements.

1. Introduction

[2] Comprehensive, global, and timely specifications of the Earth's atmosphere (particularly refractivity profiles of the troposphere and ionosphere) are required to ensure the effective operation, planning, and management of many radio frequency systems. There are currently a number of groups worldwide who are actively working on systems to provide three-dimensional (3-D) models of ionospheric electron density based on a combination of an ionospheric model with measured data. Data for such assimilative ionospheric modeling may be provided by a range of measurement techniques. Of particular note are data from the International GNSS Service (IGS) receiver network [Beutler et al., 1999]. Some IGS stations provide data in hourly files with low latency (∼15 min), which can be used to calculate near-real-time slant total electron content (TEC). Radio occultation (RO) methods are also being investigated. These measurements are made by monitoring transmissions from GPS satellites using receivers on low Earth orbiting (LEO) satellites and provide the potential for measuring refractivity profiles in regions where ground-based sensors cannot easily be located, such as deep seawaters.

[3] As different ionospheric assimilative modeling techniques are developed, it is increasingly important to undertake performance comparisons between different models and against independent ionospheric data. This paper describes two independently developed ionospheric modeling systems: the Electron Density Assimilative Model (EDAM), developed at QinetiQ under funding from the UK Ministry of Defence, and an early version of the IonoNumerics model developed by Fusion Numerics, Inc., under a contract with the U.S. Air Force Research Laboratory (http://fusionnumerics.com/ionosphere). Both these models remain under active development.

[4] EDAM uses measured slant TEC to adjust an empirical 3-D climatological distribution of ionospheric electron density. Such climatologies have been developed over many years and have been shown to capture average characteristic features of the ionosphere very well. Since EDAM does not rely on prognostic equations describing the evolution of ionospheric plasma, it can only provide simple forecasting by means of a persistence model.

[5] IonoNumerics, on the other hand, utilizes a set of time-dependent differential equations describing the conservation of mass, momentum, and energy for seven ion species and electrons. It also contains a comprehensive first-principles description of chemical transformations and energy exchange between these particles and energy inputs due to solar heating. IonoNumerics therefore does not utilize any climatological information on ionospheric electron densities and does allow forecasting via time integration of the prognostic equations. Both models can utilize measured slant TEC from a network of ground-based dual-frequency GPS receivers to constrain their results and to achieve better ionospheric nowcasting.

[6] Results from the two systems have been compared with data measured by vertical and oblique ionosondes. In the latter case, ray tracing has been used, in conjunction with the 3-D ionospheric electron density grids, to synthesize oblique ionograms. These were compared with measured oblique ionograms to assess the fidelity of the ionospheric images.

2. Assimilative Modeling Techniques

2.1. Introduction

[7] The last 15 years have seen a rapid growth in the development and use of ionospheric modeling systems. Traditional tomographic techniques represent the ionosphere by means of a two-dimensional (2-D) grid of pixels [Bernhardt et al., 1998]. Then, in the case of an integrated measurement such as TEC, the ith observation can be modeled by the sum of the electron density in the jth pixel (xj) multiplied by the ray length within the jth pixel (Hij) (Figure 1):

or in matrix notation,

where H is known as the observation operator, x is the rasterized representation of the ionosphere, and y is the TEC measurement.

Figure 1.

Diagram illustrating the construction of the observation operator that relates a slant TEC measurement to a pixel representation of the ionosphere.

[8] It can be seen that this formulation of the problem is entirely general and can therefore be extended to a 3-D grid of voxels. Furthermore, there is no requirement that the slant TEC must be measured by a ground station, and consequently, it is also simple to directly include measurements of TEC made in low Earth orbit. This generalization does require that the raypaths from the transmitter to the receiver remain straight. In practice, for GPS measurements, the departure of the rays from straight line propagation in the ionosphere can reach a few kilometers at worst, which is much smaller than the typical vertical scale height in the F region ionosphere [Schreiner et al., 1999]. Therefore bending can be neglected for this measurement geometry.

[9] The disadvantage of moving from 2-D to 3-D is that a very large amount of data is required in order to employ the traditional tomographic methods. In practice, and especially for GPS data (where the satellites are slow moving), the required data are not available, and other, more sophisticated, methods must be used. Such methods use a priori data to constrain the inversion from slant TEC to the electron density grid. This a priori data may be a model of the ionosphere [Wang et al., 2004] or may take the form of assumptions about how the ionosphere can be decomposed into a set of basis functions [Fremouw and Secan, 1992; Fremouw et al., 1994]. However, to make most effective use of a priori information, it is necessary to estimate the error covariance matrices of both the a priori information and the data. This allows the a priori information and the data to be combined in an optimal way [Rodgers, 2000].

2.2. Electron Density Assimilative Model

[10] The Electron Density Assimilative Model (EDAM) has been developed to assimilate measurements into a background ionospheric model [Angling and Cannon, 2004]. This background model is provided by the Parameterised Ionospheric Model (PIM) [Daniell et al., 1995], and the majority of the input data is GPS TEC derived from IGS stations. The assimilation is based on a weighted, damped least mean squares estimation algorithm. This is a form of minimum variance optimal estimation (also referred to as best linear unbiased estimation) that provides an expression for an updated estimation of the state (known as the analysis) that is dependent upon an initial estimate of the state (the background model) and the differences between the background model and the observations [Menke, 1989; Twomey, 1977]. The error covariance matrices of the background model and the observations are also included to control the relative contributions of the background and the observations to the analysis. The algorithm can be expressed algebraically thus:

where xa is the analysis, xb is the background model, K is the weight matrix, y is the observation vector, B is the background error covariance matrix, and R is the error covariance matrix of the observations [Rodgers, 2000]. As before, H is the observation operator that relates the measurements to the state:

where ɛ is the observation error. The analysis error covariance matrix (S) may also be calculated thus:

[11] The EDAM assimilations are conducted using a geomagnetic coordinate system that remains fixed in space with respect to the Sun. An assimilation time step of 15 min has been used, and the electron density differences between the voxels of the analysis and the background model are propagated from one time step to the next by assuming persistence combined with an exponential decay. The time constant for this decay is set at 4 hours. Thus, if the data feed is interrupted, the analysis will decay back to the background model. The variances of the background model are also propagated using a persistence model with the same time constant.

2.3. IonoNumerics

[12] IonoNumerics [Khattatov et al., 2004a, 2004b, 2004c, 2005] is a numerical global model of the ionosphere loosely based on descriptions given by Bailey and Balan [1996], Fuller-Rowell et al. [1996], Millward et al. [1996], and Huba et al. [2000]. The model computes the spatial distribution and temporal evolution of H+, O+, He+, O2+, NO+, N2+, and N+ ions as well as electrons. The model solves momentum and mass conservation equations for all seven ion species and electrons and the energy conservation equation for the three major ions and electrons. It also includes chemical interactions with neutrals and ions, recombination, ion-ion and ion-neutral collision rates, photoionization, and different types of heating.

[13] The model domain covers all latitudes and longitudes; however, implementation of polar transport and high-latitude effects are still in an experimental stage and are not considered reliable. As is customary in ionospheric applications, the dynamic equations are solved in magnetic coordinates since, in the absence of electric fields, plasma moves predominantly parallel to the direction of the magnetic field. A detailed discussion of the coordinate transformation and related equations can be found, for instance, in work by Millward et al. [1996] and Bailey and Balan [1996]. Here a brief overview is given for the benefit of readers not familiar with this subject.

[14] The Earth's magnetic field is approximated as that of a tilted eccentric dipole. Model equations are given in dipole coordinates, along magnetic field lines. The first magnetic coordinate is magnetic longitude. For each magnetic longitude, a “vertical stack” of magnetic field lines is characterized by the distance of the apex of each line from the Earth's center at the magnetic equator. This distance, normalized by the Earth's radius, is the second magnetic coordinate. Finally, for each field line, the distance from the apex to a particular point along the magnetic field line gives the third coordinate.

[15] Plasma motion along the magnetic flux tubes is strongly influenced by the ionospheric neutral wind speeds. Depending on their direction, and on the magnetic dip angle, neutral winds force plasma to move up or down field lines, strongly affecting the altitude of the maximum of the F2 layer (hmF2). IonoNumerics uses the empirical HMW-93 model for specifying these winds.

[16] Once the field-aligned transport is solved, the plasma evolution due to cross-field transport is computed. Cross-field transport is forced by electric fields either imposed externally from the magnetosphere or generated internally from the action of the neutral wind. In the lower thermosphere, the mobility of the ions is inhibited by collisions with the neutral atmosphere. The dynamo action of the neutral winds drives currents that create polarization electric fields through continuity. The ions respond to these electric fields by drifting perpendicular to both the electric and magnetic fields. This is often referred to as E × B transport. In non–fully coupled models, the E × B plasma velocity is specified from external empirical models. Once this velocity is known, solving the plasma advection equation is relatively straightforward. To model E × B drift at low latitudes, IonoNumerics uses coefficients of the Scherliess and Fejer [1999] model. In the polar regions, E × B drift is modeled using coefficients of the Weimer [1996] model.

[17] As well as providing a physical background model, IonoNumerics uses a Kalman filter approach to achieve the assimilation of data. The approach is derived from that described by Khattatov et al. [2000]. Since the model contains about 106 points, the direct implementation of a Kalman filter is impossible even with modern computing capabilities. As is customary in large-scale Kalman filter implementations in numerical weather prediction and other areas of atmospheric sciences, only the evolution of the diagonal of the background error covariance matrix is tracked, while parameterizations are used for the off-diagonal elements. If it is assumed that correlations between variations of electron densities at two different points become negligible at certain distances, the matrices become sparse, and the above-mentioned calculations become possible. Since the plasma equations are solved in magnetic coordinates, IonoNumerics specifies the decorrelation lengths and distances between points in geomagnetic rather than geographic coordinates.

3. Assimilation Scenarios

[18] For the purposes of the model comparison described in this paper, both EDAM and IonoNumerics have been limited to using data from ground-based GPS receivers. Both models have used an identical set of 66 IGS stations chosen to provide a global distribution of data (Figure 2). For both models, assimilations were conducted for 19–27 December 2003.

Figure 2.

4. Validation Data

[19] The EDAM and IonoNumerics results have been validated by comparison with data from vertical and oblique sounders (Figure 2). In the case of comparison with vertical sounders, the F2 critical frequency (foF2) above each sounder has been estimated from the maximum value of the electron density profile (i.e., NmF2) interpolated from the assimilative models' output grids to the sounders' locations:

This can then be directly compared to the value of foF2 reported by the vertical ionosondes. The height of the electron density peak (hmF2) can also be extracted from the assimilative models' output grids. This can then be compared with hmF2 values from the ionosondes. For this study, IIWG format data from 43 ionosondes have been obtained from the National Geophysical Data Center (NGDC) Space Physics Interactive Data Resource (SPIDR) database.

[20] In order to provide oblique validation data, an Improved Radio Ionospheric Sounder (IRIS) oblique chirp sounder was operated between Inskip, United Kingdom (54°N, 3°W), and Rome, Italy (42°N, 12°E). Data from this instrument were manually scaled to extract the channel's O-mode maximum usable frequency (MUF). Then a ray trace program, Segmented Method of Analytical Ray Tracing (SMART) [Norman and Cannon, 1997], was used to synthesize ionograms from the electron density grids provided by EDAM and IonoNumerics. The MUFs of these synthetic ionograms were then compared to the measured values.

5. Results

5.1. Comparison to Vertical Sounding

[21] The foF2 data from EDAM (dashed line), PIM (dotted line), and the Eglin AFB ionosonde (solid line) are shown in Figure 3 for the period 19–27 December 2003. Eglin is 824 km from the nearest IGS station used in this test. PIM has been run with the smoothed monthly sunspot number from NGDC. This results in PIM producing a repeated pattern of foF2 for each day of the test. The enhanced performance of EDAM is most apparent on 20, 21, and 27 December when the measured foF2 is well above the PIM estimate. While the EDAM results were available every 15 min, the IonoNumerics results were made available only every hour. The IonoNumerics results for Eglin are shown in Figure 4. The performance is comparable to that of EDAM, but the dawn increase in foF2 consistently lags the measurements. For this ionosonde, the daytime RMS errors in foF2 are 1.2, 0.7, and 1.0 MHz for PIM, EDAM, and IonoNumerics, respectively. The foF2 RMS errors for all the test ionosondes are shown in Figures 5 and 6 as a function of magnetic latitude. The data have been split into two groups by local time (0600–1759 and 1800–0559) at each ionosonde site. As might be expected, the errors are worst at low and high latitudes where ionospheric variability is likely to be greatest. Conversely, the data do not show a clear dependence on the distance of each ionosonde to the nearest IGS station (Figures 7 and 8).

Figure 7.

Figure 8.

Nighttime foF2 RMS error as a function of the distance to the nearest IGS station.

[22] Twenty-two of the ionosondes also reported values for hmF2. It should be noted that hmF2 is not a quantity measured by the ionosondes, but rather it is derived from a true height inversion of the ionogram [e.g., Titheridge, 1988]. Consequently, hmF2 is sensitive to errors in the (often automatic) scaling of the ionograms and tends to be a rather noisy parameter. Figure 9 shows a typical example of a noisy hmF2 measurement as well as hmF2 estimates from PIM, EDAM, and IonoNumerics. In this example, it can be seen that IonoNumerics is biased high, while the EDAM hmF2 remains close to PIM during daytime but is forced up at night. However, the uplift of the EDAM hmF2 at night does not closely track the ionosonde data.

Figure 9.

[23] RMS results from all 22 ionosondes reporting hmF2 are shown in Figures 10 and 11. Again the data have been split into local day and night and are plotted as a function of distance to the nearest IGS station. Daytime results show that the EDAM RMS errors remain close to those produced by PIM and that IonoNumerics performs somewhat more badly than the other two models. Performance of EDAM and IonoNumerics is better at night, especially at ranges of 300–600 km from an IGS station. The high RMS errors shown by IonoNumerics during the day can be attributed to a positive bias in the IonoNumerics hmF2 that is present during the day but not at night (Figures 12 and 13).

Figure 12.

Figure 13.

Nighttime hmF2 mean error as a function of the distance to the nearest IGS station.

5.2. Comparison to Oblique Sounding

[24]Figure 14 shows the measured MUFs from the IRIS oblique ionosonde (solid line), synthesized MUFs from PIM (dotted line), synthesized MUFs from EDAM (dashed line), and synthesized MUFs from IonoNumerics (dash-dotted line) for 23–25 December 2003. This date range is representative of the whole test period. The improvement of the EDAM grid over that provided by PIM (the EDAM background model) is evident during the day. However, EDAM does not always model the short-term MUF variability exhibited by the data. Furthermore, the nighttime performance of EDAM is poor, when a positive MUF bias can be observed. IonoNumerics exhibits the opposite behavior: Nighttime results are good, while MUF are underestimated during the day.

Figure 14.

[25] The mean, standard deviation, and RMS of the MUF errors are given in Table 1. As with the vertical ionosondes, the data have been divided into day and night, though in this case the local time is calculated at the path's midpoint.

6. Discussion

[26] Examination of the test results raises a number of issues that are of practical importance to developers and users of ionospheric assimilation systems.

6.1. Performance as a Function of Distance From Input Data Points

[27] Intuitively, it is expected that the performance of an assimilative model should increase as the test point moves closer to an input station. However, this does not appear to be the case for either EDAM or IonoNumerics. Performance does, however, show a clear dependence on the location of the test point. This suggests that the ionospheric variability experienced at low and high latitudes may be the limiting factor of performance for sparse networks of input stations as used in this study. In order to improve performance, it is likely that it will be necessary to increase the density of the input data and to use a mix of sensor types.

6.2. Vertical Structure of the Ionosphere

[28] It has been demonstrated that during the daytime, IonoNumerics exhibits a positive bias in hmF2 as compared to the values derived from vertical ionosondes. For practical HF applications of an assimilative model, a positive hmF2 bias can result in a negative bias in MUF even if foF2 is correctly modeled. The effects of this are seen in the poor daytime performance of IonoNumerics as compared to the oblique ionosonde. The hmF2 bias is caused by inaccuracies in the empirical neutral wind and E × B models that are used as inputs to the physical ionospheric model. This leads to the conclusion that either fully coupled physical models should be used in assimilative models or a mechanism should be included to adjust the empirical inputs during the assimilation process.

[29] It is also instructive to compare hmF2 from PIM with that of EDAM. It can be seen that during the day, PIM provides a reasonable estimate of hmF2 for this test period and that the EDAM hmF2 remains close to it. However, at night, the PIM hmF2 is biased low, and EDAM makes some adjustment to the hmF2 during the assimilative process. However, the small magnitude of the change is indicative of the difficulty in modifying the vertical structure of the model using only ground-based TEC measurements. In order to improve hmF2 performance, it is likely that additional sensors will be required that can provide higher vertical resolution, i.e., radio occultation.

[30] Another area where EDAM performs poorly is at nighttime in Europe. This is shown by the large overestimation of MUF on the United Kingdom to the Italy HF link at night and can be confirmed by an examination of the vertical ionosonde foF2 results. It is probable that this is caused by a mismodeling of the contributions to the TEC from the ionosphere and plasmasphere. The problem therefore becomes apparent at night when the relative contribution of the plasmasphere becomes much larger. Ideally, this problem would be resolved by having independent measurements of the plasmasphere TEC, i.e., from zenith-pointing GPS receivers on low Earth orbit/medium Earth orbit satellites.

6.3. Assimilating Vertical Ionosonde Data

[31] It is clear that assimilative models will require data sources that can provide information about the vertical structure of the ionosphere and that vertical ionosondes can potentially provide such information. However, in order to meet the requirements for near-real-time operation, it is unlikely that such data will be manually scaled before being assimilated. The resultant noisy hmF2 values will make effective assimilation difficult. It is probable therefore that effective assimilation of vertical ionosonde data will only be achieved by removing the intermediate step of true height inversion and assimilating the ionogram directly.

7. Conclusions

[32] Ionospheric electron density maps have been produced using two independently developed data assimilation systems (EDAM and IonoNumerics). It has been demonstrated that in comparison to a median model (PIM, which is also used as the EDAM background model), both systems are capable of reducing the RMS foF2 errors in relation to vertical HF sounding. It has also been shown that problems in effectively retrieving height information from the assimilated data can have a significant impact on the models' performance.

[33] This study has highlighted that testing against sounder information can provide important validation of the many ionospheric imaging and data assimilation techniques that are currently emerging. Indeed, it was only by testing against oblique ionosonde measurements that the hmF2 problems in IonoNumerics were exposed and can now be addressed. This has shown that only using other TEC measurements to provide validation, be they from GPS, TOPEX, or some other source, may not provide a sufficiently rigorous test of a model that is itself largely driven by TEC inputs.

Acknowledgments

[34] EDAM has been developed under funding from the Operating Environment domain of the United Kingdom Ministry of Defence Corporate Science and Technology program. IGS data were obtained from the IGS, JPL, and SOPAC data canters. Differential code biases were obtained from the Centre for Orbit Determination in Europe. Vertical ionosonde data were obtained via SPIDR. BK is grateful to the Air Force Research Laboratory personnel for useful discussions and guidance and to T. Fuller-Rowell, M. Codrescu, M. Murphy, M. Gnedin, J. Wright, N. Zabotin, and F. Berkey. AFRL has funded development of IonoNumerics by Fusion Numerics, Inc.

Index terms:

Supporting Information

Please note: Wiley-Blackwell is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.

3M. Pietrella, Achievement of a short term three dimensional electron density mapping of the ionosphere in the European sector: Comparisons with the IRI model for quiet-moderate geomagnetic-ionospheric conditions, Advances in Space Research, 2016, 58, 7, 1242CrossRef