Number of analyzed days in specified data groups according to F10.7 range and month. The last column shows mean number of days per F10.7 range.
Nowadays, the Global Navigation Satellite System (GNSS) is commonly used for positioning, navigation and timing. The GNSS based services are used in many areas such as maritime, aviation, agriculture, public transportation and geodesy. The GNSS receiver computes its position by trilateration using ranges between satellites and the receiver, where the ranges are calculated from measurements of time-of-arrival of satellite radio signals . However, the signals do not propagate ideally. Many factors can change signals’ propagation speed or trajectory and they can consequently cause incorrect determination of the receiver position.
One of the factors which affect GNSS signal propagation is the ionosphere. The ionosphere causes delay of radio signals, and if not mitigated, it can be the largest source of error (ionospheric error) in GNSS positioning and navigation . There are several possibilities to compensate for the ionospheric effect. First technique is to use multi-frequency satellite-receiver communication which takes advantages of dispersive nature of the ionosphere. This approach called ionosphere-free combination can remove about 99 % of the ionospheric error . In case the receiver uses only one frequency it can use a Satellite Based Augmentation Systems (SBAS) such as the Wide Area Augmentation System (WAAS) or European Geostationary Navigation Overlay Service (EGNOS). These systems determine condition of nearby ionosphere from a network of reference stations and send the information to the user via a geostationary satellite . If the SBAS service is not available or it is not supported by the receiver, an ionospheric model can be used to estimate the ionospheric error. Ionospheric models are also used for satellite and receiver inter-frequency bias estimation and Total Electron Content (TEC) calibration .
There are several empirical ionospheric models. Well known empirical models included in our study are the Klobuchar model, International Reference Ionosphere (IRI) and NeQuick. In addition, we evaluated accuracy of a relatively new model developed at DLR (Deutsches Zentrum für Luft-und Raumfahrt), the Neustrelitz TEC Model (NTCM) . As each of the models applies different modeling approach and was developed with different background data, we assume that their TEC modeling performance differ and the use of one model in particular condition would be better than use of another.
Some of the recent analysis evaluated TEC modeling performance of the Klobuchar model and the NeQuick 2 by comparison of modeled TEC with GNSS measurements , . In our study, we evaluated performance of four empirical ionospheric models comparing them with Global Ionospheric Maps (GIMs) produced at the Center for Orbit Determination in Europe (CODE), using a similar approach as in  for the Klobuchar model and the NeQuick 2. This approach has the advantage that all TEC data are in the zenith direction, therefore, we can avoid conversion between slant and vertical TEC which can produce additional error . On the other hand, GIMs do not contain direct measurements and we have to consider their data accuracy.
To estimate TEC mismodeling for each model, we compared modeled TEC data with CODE GIMs over three years (2010, 2011 and 2012). Such data were produced at the CODE using measurements from about 200 GPS/GLONASS stations. The evaluation will show variation in models' TEC modeling performance with respect to time of day, season, location and space weather condition. In addition, it will show in which cases the CODE GIMs' inaccuracy prevents performance evaluation.
2. The ionosphere
The ionosphere is the part of the atmosphere with large amount of charged particles (ions and electrons). A typical vertical profile of ionospheric electron density (Figure 1) is divided into the several layers according to the different ionization and recombination principles . The ionospheric structure significantly varies with geographical location, local time and with changes in solar-terrestrial environment.
We can divide ionospheric variations into two groups. The first group includes variations with periodic behavior that can be distinguished from empirical data. These variations are: daily variation, seasonal variation, dependence on the geomagnetic field and climatological dependence on space weather. These phenomena can be analytically described and modeled by empirical ionospheric models.
The second group of ionospheric variations includes sudden ionospheric disturbances (SID) or small rapid changes in the electron density causing scintillation. Even though these phenomena are often observed, they do not show any behavior pattern to the magnitude or period of occurrence . As these variations belong to ionospheric weather rather than to climatology they are not modeled by empirical ionospheric models.
Radio signal which propagates through the ionosphere experiences changes of its propagation speed and trajectory. These changes depend on the signal carrier frequency and electron density of the ionosphere . The ratio between the group propagation velocity and the speed of light in vacuum can be described as refractive index :
The refractive index of the ionosphere can be derived from the Appleton-Hartree formula. Usually, we considering only the first two terms of the equation as the rest contributes to the ionospheric error by less than 1% . The group refractive index then leads to
where is the electron density and is the signal carrier frequency. If we write the ionospheric propagation time as integration of the propagation speed over the propagation path and include (2) we get
The integration of the electron density along the path is usually referred to as Total Electron Content (TEC):
which is the total number of electrons in a tube of 1 m2 cross-section along the GNSS signal path through the ionosphere. If we know the frequency of the signal, the value of TEC allows us to compute the propagation delay introduced by the ionosphere, which is an error for our point of view. All ionospheric models considered in this study can provide TEC values.
There are two commonly used types of TEC values: vertical and slant TEC. Vertical total electron content (vTEC) at a certain geographical point stands for TEC in the direction of the zenith. TEC between a satellite and a receiver is usually referred to as the slant TEC (sTEC). This notation signifies that the TEC is at a different angle then the zenith. TEC is usually given in TEC Units (TECU) where 1 TECU=1016 electrons/m2.
3. Ionospheric models
3.1. Klobuchar model
The Klobuchar model was developed by John A. Klobuchar at the Air Force Geophysics Laboratory, U.S. The algorithm is used to correct ionospheric time-delay in GPS for single frequency communication. The model was developed in 1975 keeping in mind limited computation memory and capability of receivers, therefore, the model algorithm is very fast and has minimum complexity.
One of the main criteria of the algorithm design was to fit best the daily period with the largest TEC values, i.e. afternoon period. The Klobuchar model approximates daytime variation of ionospheric time delay as a half period of cosine function with maximum at 14 hours local time. The amplitude and period of the cosine are each calculated with 4 coefficients transmitted by GPS navigation message. The night time ionospheric delay is set as a constant value of 5 ns (Figure 2.) .
The first version of the International Reference Ionosphere was developed as a joint project of Committee on Space Research (COSPAR) and Union of Radio Science (URSI) in 1978. Since then, the model has been continuously improving.
The IRI is able to compute vertical electron density profile and vTEC as well as other ionospheric parameters such as ion densities and ion temperatures. The IRI divides the ionosphere into six sub-regions where each of them is described by several parameters. The model also uses lists of foF2 and M(3000)F2 parameters (foF2 is critical frequency of the ionospheric layer F2 which is usually the layer with Ne maximum, M(3000)F2 is ratio between maximum usable frequency for ionospheric radio link over the distance of 3000 km using layer F2 and foF2, more information can be found in ). The integration height for TEC computation is limited to 2000 km .
To calculate vertical el. density profile of the ionosphere, user can choose one from several models for each region. In this work, we used standard IRI setting. The IRI uses both space weather and geomagnetic indices, i.e. Solar Radio Flux (F10.7) index, International Sunspot Number (Ri) and magnetospheric Ap index which describes variations in the Earth's geomagnetic field.
Comparing to the Klobuchar model, the ionosphere structure modeled by the IRI is more complex including also the equatorial anomaly (Figure 3). Equatorial anomaly is the area with higher TEC distanced about 20° north and south from the equator. Source code of current IRI version is available on model's webpage (http://irimodel.org/).
3.3. NeQuick 2
The NeQuick is an empirical model based on the model introduced by Di Giovanni and Radicella in 1990. Its modified version is used in the GNSS Galileo to aid single-frequency positioning . The model has been also included into the ITU-R recommendation as a suitable method for TEC modeling. In addition, the IRI model uses NeQuick algorithm as a default option for the upper ionosphere computation.
The NeQuick is able to calculate electron density at any given location in the ionosphere. Therefore, it can provide TEC and electron density profile between any two given points . For the analyses, we used the NeQuick version 2 which we obtain at the International Centre for Theoretical Physics in Trieste, Italy. The model is driven by monthly-mean solar radio flux.
NeQuick is a complex electron density model. As well as in case of the IRI, we can distinguish equatorial anomaly on the NeQuick GIM (Figure 4).
Recently, a new global ionospheric model NTCM was developed at the Institute of Communications and Navigation, DLR in Neustrelitz, Germany. The model can provide values of vTEC at any given time and location. The core of the model consists of 12 coefficients which can be autonomously used for full solar cycle. The driver of the NTCM is the F10.7 index. The model does not use any integration of electron density profile, therefore, it is very simple and fast . The model analytically describes daily variation, seasonal variation, equatorial altitude anomaly and solar flux dependency as harmonic functions. All the formulas of the model algorithm can be found in .
The NTCM, as well as the Klobuchar, is a TEC model. The GIM structure is rather simple and more similar to the Klobuchar one than to the GIMs produces with the electron density models IRI and NeQuick (Figure 5).
4. Data and methodology
All the ionospheric models discussed here are climatological. For our comparison, we chose CODE as a form of climatological reference rather than real measurements which can be affected by local ionosphere weather.
The Center for Orbit Determination in Europe at the Astronomical Institute at University of Berne, Switzerland provides GIMs on daily bases from 1995. The maps are available in the IONosphere map EX change format (IONEX) from 1997. CODE GIMs cover area from 87.5° northern to 87.5° southern latitude and from 180° western to 180° eastern longitude. The grid point step is 2.5° in latitude and 5° in longitude. GIMs are produced with 2 hours interval and the maps are generated using data from about 200 GPS/GLONASS sites of the IGS (International GNSS Service) and other institutions. Each map grid point contains a vTEC value calculated from measurements from the GNSS sites .
The accuracy of the CODE vTEC values depends on the local density of the GNSS reference network. The reference stations are not equally distributed over the whole world and CODE has higher level of inaccuracy at places with lack of reference stations, typically over oceans. Along vTEC values, CODE IONEX files contain also RMS maps (Figure 6) which give value of RMS error for each GIM grid point.
As the goal is a climatological study a lot of data are required. In this study, we processed data from years 2010, 2011, and 2012. Considering 12 maps per day and 2012 as a leap year, we made comparisons with 13,152 CODE maps. As each map has 71x73 (±87.5 latitude with 2.5° step, ±180 longitude with 5° step) values, we processed 68,166,816 CODE values in total.
First, we created similar GIM databases as the CODE one with all 4 models so that, for each CODE map, we created 1 map for each ionospheric model (4 maps per 1 CODE map). The modeled maps have the same grid points as CODE maps and each grid point contains corresponding modeled vTEC value.
As the next step, we divided data into groups. First division was according to the universal time (UT). For example, maps for 12 UT were compared only with each other and not with maps for different UT. As the CODE time resolution is 2 hours, this division created 12 data groups. Secondly, we divided data according to months in order to analyze effect of seasonal variation. This division created 12 month sub-groups for each of the 12 UT data groups.
The last division criterion was solar activity. Solar radio flux F10.7 and international sunspot number Ri are considered to be two primary long-term solar indices. The F10.7 is often used as a proxy for Ri making these two indices interchangeable, however, a recent study showed that there is a disagreement during the last decade . We chose F10.7 over the Ri as we assume that the index derived from measurements on Earth surface corresponds more to the ionospheric behavior rather than index derived from measurements of a solar phenomena (sun spots). The Figure 7 shows the observed solar radio flux variation for years 2010, 2011, and 2012. The data for the 1st of January 2011 and 2012 were missing and were filled as linear interpolation of values of adjacent days.
Considering 3 years of data, each group has approximately 91 days (e.g. March group has 93 days as 3 years mean 3 Marches (3x31=93)). To keep the number of days within one F10.7 group high enough we decided to divide the groups into only 4 F10.7 sub-groups. Applying this division, each sub-group should have approximately 23 days where every sub-group includes only days within particular F10.7 range. We divided the dataset in order to keep the amount of days in groups as balances as possible, which led to non uniform distribution of F10.7 intervals (Figure 7). Even so, some groups have significantly less or more days then ideal 23 days (Table 1).
|F10.7 range [sfu]||Number of days per group||mean|
|0 – 82||32||13||11||30||32||30||19||20||14||13||13||12||19.9|
|83 – 100||30||37||28||9||17||27||40||30||19||20||23||25||25.4|
|101 – 122||6||34||38||43||30||15||12||30||20||17||10||26||23.4|
4.3. Comparison method
We estimated the TEC mismodeling of all models comparing the modeled ionospheric maps with the reference CODE maps for each data group. We calculated the mismodeling as mean absolute difference between model's and CODE's for each map grid point as
where is the reference CODE vTEC value and is the corresponding vTEC value modeled by the model m1. The number N stands for the number of analyzed maps and depends on particular data group. For example, for the month January and F10.7 range of 83 – 100 sfu the N is 30 (Table 1). We computed the mean absolute difference for each UT separately.
As it was mentioned, the accuracy of CODE data varies and should be considered. We calculated mean CODE RMS error for each grid point of each data sub-group as
where the N is again the number of analyzed maps for each data group.
The decision of the most accurate model was made for each map grid point of each data group. Referring to Figure 8, in case the for all models for particular grid point was higher than the corresponding , the model with the lowest is marked as decisively the best model of the grid point. Also, in case only one model has below corresponding this model is marked as decisively the best model. In case two or more models have lower than corresponding , the model with lowest is still identified as the best in average but not decisively with respect to the other models within the threshold, because the CODE accuracy for that particular grid point is not high enough.
Results for each field of Table 1 are represented in form of grid maps. As there are 48 groups and each of them has 12 UT-sub-groups we do not show all the results but only two examples. It should be noted that all the shown maps display interpolated data (1° x 1°) but the original grid resolution is 2.5° in latitude and 5° in longitude.
The first example (Figure 9) shows results for October, 14 UT and F10.7 range of 101 – 122 sfu. Figure 9a shows calculated by the best models for particular areas. The average value of is 3.7608 TECU and 63 % of the values are lower than the average. Higher values of are mostly at the area of equatorial anomaly. The spatial distribution of the best models over the globe is shown in Figure 9b. This maps shows which model has the lowest to CODE for particular location. If we apply the criterion of the (Figure 9c) we can identify the regions for which the insufficient accuracy of CODE data prevents identification of the decisively best model (marked with gray-white stripes, Figure 9d). In these areas two or more models have their lower then . The Figure 9e additionally shows regions where only two models have lower then . The areas are marked with stripes with colors of the particular models.
The second example of results is for January, 02 UT for 0 – 82 sfu data group (Figure 10). As it was expected, the values are lower during the periods with low solar radio flux. The average value of is 1.755 TECU and 59 % of the values are lower than the average (Figure 10a). The best models distribution can be seen in Figure 10b and the map considering the mean RMS error criterion in Figure 10c. The amount of the gray-white areas indicates that for this data group two or more models have their lower then in majority of cases. Additionally, Figure 10d shows that for about half of the globe three or all four models have their lower then (gray-white areas).
The results can be also expressed in form of a graph for particular location. Both Figure 11 and 12 show performance of the models for the location: 60° northern latitude and 15° eastern longitude. Results for October, 14 UT and 101 – 122 sfu is shown in Figures 11 and results for January, 02 UT and 0 – 82 sfu in Figure 12. The graph for October shows that in all cases there is no more than one model with lower than so that the best model was decisively identified for the whole day. On the other hand, for January, it was not possible to identify the best model decisively during the morning and evening and night as for these periods more models have lower then .
If we summarize results for all data groups we can show how much were models identified as the best with respect to different areas (Table 2). The EU region was chosen between 65° – 30° northern latitude and 10° western to 50° eastern longitude and the USA between 55° – 0° northern latitude and 130° – 50° western longitude. Similar results but with respect to F10.7 ranges are shown in Table 3. Table 4 summarizes the amount of cases for which the analysis was able to decisively determine the best model with respect to both F10.7 and different regions.
|Model||Portion cases [%]|
|All data [%]||Decisive EU [%]||Decisive USA [%]|
|Region||Portion of cases [%]|
|F10.7 range [sfu]|
|0 – 82||83 – 100||101 – 122||> 122|
|Region||Portion of cases [%]|
|F10.7 Range [sfu]|
|0 – 82||83 – 100||101 – 122||> 122|
One of the constrains of this research is the small amount of solar radio flux sub-groups. In our case, the same model is marked as the best one for particular location and time for both F10.7 of 83 sfu and 100 sfu. A finer solar flux scale would be more adequate. We used wide solar flux ranges to keep amount of days in data groups high enough. However, there are still 4 groups with number of days lower then 10 (one group with only 1 day) which we consider to be insufficient. Future research can overcome this issue by including more years into the analysis, preferably a whole solar cycle.
Considering Table 2 and 3, the Klobuchar and NTCM were marked as the best models in more cases than the IRI and NeQuick. Such results are not in agreement with results from  and  where Klobuchar always performs worse than the NeQuick. However, in our case, both the IRI and NeQuick are driven by averaged indices while NTCM and Klobuchar by daily indices. This provides the advantage for NTCM and Klobuchar of ability to respond on any rapid day-to-day variation of the ionosphere. In the studies  and  the NeQuick is driven by Ionization level, while we used monthly mean F10.7. It can be expected that applying Ionization level would significantly change our results, this will be done in future.
Good performance of the Klobuchar model in our test is surprise and it will be better investigated in the future work using larger database of reference data, different indices to drive the models and comparing results for different reference data sources.
The Table 4 shows that we were able to decisively identify the best model for most of the cases in the EU and USA region, especially during the middle and high solar activity. This is caused by the fact that during the periods with higher solar activity the difference between modeled TEC and the CODE's TEC rises while CODE RMS error stays roughly the same. In particular for these cases the information about the best model can be very important to minimize the potential impact of the ionosphere mismodeling in single frequency positioning. However, to verify this assumption the results should be tested on real TEC measurements.
We compared TEC data modeled by empirical ionospheric models: IRI2012, Klobuchar, NeQuick2 and NTCM to the CODE TEC data. We analyzed CODE GIMs for every two hours for years 2010, 2011 and 2012. The results show that the CODE RMS error values are low enough to identify the decisively best model in most cases above Europe and North America, especially during the days with higher solar radio flux. The ability to decisively recognize the best model decreases for lower solar radio flux values. For these periods more models have mean absolute TEC deference lower than corresponding mean CODE RMS error. The study shows that the Klobuchar and NTCM were marked as the best model in more cases than NeQuick 2 and IRI 2012. However, the performance of all models varies according to the time, location and solar flux which indicates that there is not one model which performs the best under all conditions. In addition, there is still a significant portion of cases for which the best model could not be identified decisively.
In the future work, we plan to analyze data for more years, improve the solar radio flux resolution and include the Galileo version of the NeQuick into our analysis. We also plan to perform the analysis on the different reference databases and test the results on real TEC measurements.
Pavel Najman's research work is undertaken in the scope of the TRANSMIT ITN, funded by the Research Executive Agency within the 7th Framework Program of the European Commission, People Program, Initial Training Network, Marie Curie Actions-GA no 264476.
We would like to thank to Astronomical Institute of the University of Bern for freely available CODE maps. We would also like to thank to Bruno Nava and Sandro M. Radicella for the NeQuick 2 source code.