Journal topic
Ann. Geophys., 38, 331–345, 2020
https://doi.org/10.5194/angeo-38-331-2020
Ann. Geophys., 38, 331–345, 2020
https://doi.org/10.5194/angeo-38-331-2020

Regular paper 16 Mar 2020

Regular paper | 16 Mar 2020

# Global total electron content prediction performance assessment of the IRI-2016 model based on empirical orthogonal function decomposition

Global total electron content prediction performance assessment of the IRI-2016 model based on empirical orthogonal function decomposition
Shuhui Li1,2, Jiajia Xu1, Houxiang Zhou1, Jinglei Zhang1, Zeyuan Xu1, and Mingqiang Xie1 Shuhui Li et al.
• 1School of Land Science and Technology, China University of Geosciences (Beijing), Beijing 100083, China
• 2Shanxi Key Laboratory of Resources, Environment and Disaster Monitoring, Jinzhong 221116, China

Correspondence: Shuhui Li (li.shuhui@163.com)

Abstract

In this study, the empirical orthogonal function (EOF) decomposition technique was utilized to analyze the similarities and differences of the spatiotemporal characteristics between the total electron content (TEC) of the International GNSS Service global ionospheric map (GIM) and that derived from the International Reference Ionosphere 2016 (IRI-2016) model in 2013. Results showed that the main spatial patterns and time-varying features of the data set have good consistency. The following four main spatiotemporal variation features can be extracted from both data sets through EOF decomposition: the variation with the geomagnetic latitude reflecting the daily averaged solar forcing, the diurnal and semidiurnal periodic changes with longitude due to local time, and the interhemispheric asymmetry caused by the annual variation of the inclination angle of the Earth's orbit. The differences between the spatial patterns represented by the EOF base functions of IRI-2016 and GIM TECs were analyzed by extracting the same time-varying coefficients. The deviations of the interhemispheric asymmetry component between the two data sets showed roughly equal values throughout the Southern or Northern Hemisphere, whereas those of the other spatial modes were mainly concentrated on the equatorial region. The differences of the time-varying characteristics between the IRI-2016 and GIM TECs were also compared by extracting the same EOF base functions. Although the EOF coefficients of the two data sets presented consistent seasonal variations, the magnitude of IRI-2016 TEC changes over time was less than that of GIM TEC. The diurnal variation of the daily averaged solar forcing component and the annual variation of the interhemispheric asymmetry component exhibited relatively large deviations between the two data sets. Considering the variance contribution of the different EOF components and their average relative deviations, both analyses showed that the daily averaged solar forcing and interhemispheric asymmetry components were the main factors for the deviation between the IRI-2016 and GIM TECs.

1 Introduction

The ionosphere is a shell of electrons and electrically charged atoms and molecules that surrounds the Earth and stretches from a height of approximately 60 km to more than 1000 km. The variations in the ionosphere should be accurately measured, modeled, or estimated because the ionosphere critically affects high-frequency satellite communication and navigation system signals. Total electron content (TEC), which is the number of free electrons along the path where the signal is traveling, is a critical quantity that describes the ionosphere and its variability. Modeling and predicting temporal and spatial variations in ionospheric TEC are crucial to ionospheric physics research and ionospheric-based applications (Yao et al., 2018).

Many attempts have been made to specify ionospheric parameters using empirical approaches, because an empirical model can describe the general condition of the ionosphere without actual measured data (Feltens et al., 2011). Several ionosphere empirical models, such as Klobuchar, NeQuick, the Standard Plasmasphere Ionosphere Model (SPIM), and International Reference Ionosphere (IRI; Bilitza, 2001), are currently available. The IRI is one of the most accepted standard global empirical ionosphere models. This model can be used to estimate the values of electron density and temperature, ion temperature and composition, and TEC at altitudes ranging from approximately 50 to 2000 km at a particular location, at a particular time, and on a particular day. The IRI model is continuously improved when new data and techniques become available. This model was recently upgraded to the IRI-2016 version (Bilitza et al., 2017). The model has been improved by ingesting all available data from worldwide ground-based and satellite observations to enhance the model capacity. IRI-2016 includes two new model options for the F2 peak height hmF2 and an enhanced representation of topside ion densities at low and high solar activities. Several small changes were made concerning the use of solar indices and the speedup of the computer program (Bilitza et al., 2017).

The performance of the previous versions of the IRI model in terms of predicting TEC have been investigated to improve the model effectively and provide reference for the application (Maltseva et al., 2012; Scidá et al., 2012; Kenpankho et al., 2013; Okoh et al., 2013; Zakharenkova et al., 2015; Li et al., 2016). Comparative studies with GNSS-derived TEC have validated the performance of different IRI versions over years of varied solar activity in diverse regions. Given the predictability of the diurnal variation of TEC, deficiencies have varied with local time (LT), season, and latitude. After the release of IRI-2016 as the recent version, its performance in predicting TEC has attracted the attention of many researchers (Atici, 2018; Sharma et al., 2018; Tariku, 2018; Jiang et al., 2019). Most existing studies for ionospheric models aimed at the low and middle latitudes. Studies on the TEC prediction performance of different IRI versions worldwide are relatively sparse. Most comparative studies are based on the contrast between TEC derived from the IRI model and that derived from the global ionospheric map (GIM) or GNSS. The variations of diurnal and seasonal changes and those in different solar activity years on certain sites have been investigated from several aspects, such as bias, root mean square (RMS) error, and correlation coefficients. Although several assessments of the IRI models have been conducted, few studies on the comprehensive evaluation of the temporal and spatial distribution prediction performance of the IRI model are available. The predictive performance of the IRI model for ionospheric temporal and spatial changes should be evaluated using efficient analytical methods.

Many scholars have recently used the empirical orthogonal function (EOF) decomposition method to analyze the spatial patterns and temporal variations of the TEC and their relationships with influencing factors (Zhao et al., 2005; Mao et al., 2008; A et al., 2011; Zhang et al., 2011, 2013; Bouya et al., 2012; Uwamahoro and Habarulema, 2015; Talaat and Zhu, 2016; Dabbakuti and Ratnam, 2016, 2017; Chang et al., 2017; Andima et al., 2019; Li et al., 2019). The spatial patterns and temporal variations of the TEC are separated by EOF decomposition and can be properly represented by the base functions and associated coefficients, respectively. The data analysis results of a single station and the regional or global TEC indicated that the EOF method is a potentially useful tool for data compression and separation of different physical processes. The EOF method contributes to the comprehensive analysis of the overall spatiotemporal variations in ionospheric TEC.

In this work, GIM TEC data in 2013 were selected as reference values, and the EOF method was introduced to analyze the global TEC prediction performance of IRI-2016. A comparison between the modeled TEC and the reference values was conducted from the perspective of spatial patterns and time variation characteristics. Results provide a reference for the further understanding of the differences between the IRI-2016 and the GIM TECs at a global scale.

2 Data and method

## 2.1 GIM TEC

The GIM TEC used in this study is the official IGS combined final product provided by the Crustal Dynamic Data Information System (ftp://cddis.gsfc.nasa.gov, last access: 10 April 2019). Final GIMs are regular products of the International GNSS Service (IGS) since 1998. These GIMs are provided in the ionosphere exchange format with a spatial resolution of 2.5× 5 in geographic latitude and longitude and a temporal resolution of 2 h.

In this study, we downloaded and extracted the 2013 global TEC data from GIMs (referred to as GIM-TEC hereafter).

## 2.2 IRI-2016

The IRI is the international standard empirical model for the terrestrial ionosphere and recommended for international use by the Committee On Space Research and International Union of Radio Science (Bilitza, 2001; Bilitza and Reinisch, 2008; Chauhan and Singh, 2010). The first version was released in 1978, followed by several steadily improved ones in 1986, 1990, 1995, and 2012 (Rawer et al., 1978; Bilitza, 2015). The most recent version of this model is IRI-2016 (Bilitza et al., 2016, 2017). After IRI-2012, IRI-2016 exhibits the latest improvement in the model by introducing two new F2 peak height hmF2 modeling options with their data sources from ionosonde measurements (Altadill et al., 2013) and COSMIC radio occultations (Shubin, 2015). Hence, this version is independent of the propagation factor M(3000)F2 (Bilitza et al., 2017).

The software package of IRI-2016 can be downloaded from http://irimodel.org/ (last access: 12 March 2019). The IRI software package contains FORTRAN subroutines, model coefficients, index files for IRI-2016 models, README files, and license files. The user can calculate relevant parameters by inputting location, time, height range, model selection, and certain parameters. The global TEC data calculated by using IRI-2016 will be called IRI-TEC hereafter. IRI-TEC can also be calculated online in accordance with https://ccmc.gsfc.nasa.gov/modelweb/models/iri2016_vitmo.php (last access: 1 March 2020).

## 2.3 EOF decomposition

The EOF decomposition analysis method was originally invented by Pearson (1901). This method is performed by using an orthogonal transformation to decompose the original data set into a set of uncorrelated and ordered base functions and associated coefficients.

If an original data matrix X with the dimension M×N is present, then the covariance matrix is determined from the data matrix X in accordance with

$\begin{array}{}\text{(1)}& \mathrm{\Sigma }={\mathbf{X}}^{\mathrm{T}}\mathbf{X}.\end{array}$

The EOF base functions Ei, with i=1, 2, 3, …, N, are the eigenvectors of the covariance matrix and obtained by solving

$\begin{array}{}\text{(2)}& \mathrm{\Sigma }{E}_{i}={\mathit{\lambda }}_{i}{E}_{i},\end{array}$

where λi is the associated eigenvalues. Once the EOF base functions are known, the EOF coefficients Ak are obtained using

$\begin{array}{}\text{(3)}& {A}_{k}=\mathbf{X}{E}_{k}.\end{array}$

The original data set X can be decomposed in terms of the EOF base functions and associated coefficients in accordance with

$\begin{array}{}\text{(4)}& \mathbf{X}=\sum _{k=\mathrm{1}}^{N}{E}_{k}{A}_{k}.\end{array}$

The percentage of the total variance in the data set accounted for by the ith EOF component is given as follows:

$\begin{array}{}\text{(5)}& {r}_{i}=\mathrm{100}×\frac{{\mathit{\lambda }}_{i}}{\sum _{j=\mathrm{1}}^{N}{\mathit{\lambda }}_{j}}\mathit{%},\end{array}$

where N denotes the total number of the EOF components accounting for the total variance in the original data set.

Talaat and Zhu (2016) reported that the effectiveness of the EOF technique for TEC is nearly insensitive to the horizontal resolution and length of the data records. We analyzed the global TEC over a 1-year time period (2013) with a 2 h temporal resolution and 37×36 spatial grids.

We first organized the data set $\mathrm{TEC}\left(\mathrm{Lat},\mathrm{Lon},\mathrm{UT},\mathrm{DOY}\right)$ used in this study into a 2D matrix according to location and time epoch, that is, TEC(epoch,grid), where grid is a grid point arranged according to the latitude and longitude, and its total number is $\mathrm{37}×\mathrm{36}=\mathrm{1332}$; epoch is arranged according to Universal Time (UT), with an interval of 2 h. The total epoch number of the study period was $\mathrm{12}×\mathrm{365}=\mathrm{4380}$. After performing EOF decomposition, the base function Ek(grid) expressing a spatial pattern and the associated coefficient Ak(epoch) varying with time are obtained.

The EOF method can separate the temporal and spatial variation characteristics. If the IRI TEC and GIM TEC are decomposed separately, it is difficult to directly compare their EOF base functions and coefficients in magnitude. Therefore, we combined the data to form a whole data set for EOF decomposition and compared the two data sets.

The same coefficients of the EOF base function, that is, the same time-varying features, can be obtained by arranging IRI-TEC and GIM-TEC according to the same number of columns. Accordingly, comparing the two data sets' spatial variation features represented by the base functions is feasible.

$\begin{array}{}\text{(6)}& \begin{array}{rl}\left[\begin{array}{c}{X}_{\mathrm{GIM}}\\ {X}_{\mathrm{IRI}}\end{array}\right]& =\sum _{k=\mathrm{1}}^{N}\left[\begin{array}{c}{E}_{k,\mathrm{GIM}}\\ {E}_{k,\mathrm{IRI}}\end{array}\right]\cdot {A}_{k}\\ & =\left[\begin{array}{c}{\sum }_{k=\mathrm{1}}^{N}{E}_{k,\mathrm{GIM}}\cdot {A}_{k}\\ {\sum }_{k=\mathrm{1}}^{N}{E}_{k,\mathrm{IRI}}\cdot {A}_{k}\end{array}\right]\end{array}\end{array}$

If IRI-TEC and GIM-TEC are arranged in the same number of rows, then the same spatial variation features represented by EOF base functions will be obtained. Accordingly, the time variation characteristics of the two data sets can be compared.

$\begin{array}{}\text{(7)}& \begin{array}{rl}\left[{X}_{\mathrm{GIM}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{X}_{\mathrm{IRI}}\right]& =\sum _{k=\mathrm{1}}^{N}{E}_{k}\cdot \left[{A}_{k,\mathrm{GIM}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{A}_{k,\mathrm{IRI}}\right]\\ & =\left[\sum _{k=\mathrm{1}}^{N}{E}_{k}\cdot {A}_{k,\mathrm{GIM}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\sum _{k=\mathrm{1}}^{N}{E}_{k}\cdot {A}_{k,\mathrm{IRI}}\right]\end{array}\end{array}$

## 2.4 Evaluation indicators

In this study, the mean bias was calculated to represent the difference between two data sets. The equation is shown as follows:

$\begin{array}{}\text{(8)}& \mathrm{Bias}=\frac{\mathrm{1}}{n}\sum _{i=\mathrm{1}}^{n}\left({Y}_{i}-{Y}_{i}^{\prime }\right),\end{array}$

where n is the total number of sample data, and Yi and ${Y}_{i}^{\prime }$ are sample data for two different data sets. These variables can be TEC from IRI-2016 and GIMs or the values of base functions or coefficients of base functions. The mean relative bias (Bias_rel) can be calculated as follows:

$\begin{array}{}\text{(9)}& \mathrm{Bias}\mathrm{_}\mathrm{rel}\mathit{%}=\frac{\mathrm{1}}{n}\sum _{i=\mathrm{1}}^{n}\frac{\left({Y}_{i}-{Y}_{i}^{\prime }\right)}{{Y}_{i}^{\prime }}×\mathrm{100}.\end{array}$

The RMS error of the bias can be calculated using the following expression:

$\begin{array}{}\text{(10)}& \mathrm{RMS}=\sqrt{\frac{\mathrm{1}}{n}\sum _{i}^{n}\left({Y}_{i}-{Y}_{i}^{\prime }{\right)}^{\mathrm{2}}}.\end{array}$

The 2D linear correlation coefficient was used to investigate the similarity of the spatial pattern of IRI-TEC and GIM-TEC. The 2D linear correlation coefficient ρ for two matrices A and B with M×N dimension is calculated as

$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathit{\rho }& =\sum _{m=\mathrm{1}}^{M}\sum _{n=\mathrm{1}}^{N}\left({A}_{mn}-\stackrel{\mathrm{‾}}{A}\right)\left({B}_{mn}-\stackrel{\mathrm{‾}}{B}\right)\cdot {\left[\sum _{m=\mathrm{1}}^{M}\sum _{n=\mathrm{1}}^{N}\left({A}_{mn}-\stackrel{\mathrm{‾}}{A}{\right)}^{\mathrm{2}}\right]}^{-\frac{\mathrm{1}}{\mathrm{2}}}\\ & \cdot {\left[\sum _{m=\mathrm{1}}^{M}\sum _{n=\mathrm{1}}^{N}\left({B}_{mn}-\stackrel{\mathrm{‾}}{B}{\right)}^{\mathrm{2}}\right]}^{-\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$

where $\stackrel{\mathrm{‾}}{A}$ and $\stackrel{\mathrm{‾}}{B}$ are the mean values of matrices A and B, respectively, and they are written as

$\begin{array}{}\text{(12)}& \stackrel{\mathrm{‾}}{A}=\frac{\mathrm{1}}{MN}\sum _{m=\mathrm{1}}^{M}\sum _{n=\mathrm{1}}^{N}{A}_{mn};\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\stackrel{\mathrm{‾}}{B}=\frac{\mathrm{1}}{MN}\sum _{m=\mathrm{1}}^{M}\sum _{n=\mathrm{1}}^{N}{B}_{mn}.\end{array}$
3 Results and analysis

## 3.1 GIM-TEC and IRI-TEC in 2013

Figure 1 shows the season averages of global GIM-TEC and IRI-TEC at 12:00 UT in 2013. The months are divided into the following four seasons: March equinox (February, March, and April), June solstice (May, June, and July), September equinox (August, September, and October), and December solstice (November, December, and January). The global level of ionospheric TEC at 12:00 UT is lowest during the June solstice compared with that during other seasons. By contrast, the ionospheric TEC reaches the highest level during the December solstice.

Figure 1Season averages of global TEC obtained from GIM and IRI at 12:00 UT in 2013. (a) GIM-TEC in the March equinox; (b) IRI-TEC in the March equinox; (c) GIM-TEC in the June solstice; (d) IRI-TEC in the June solstice; (e) GIM-TEC in the September equinox; (f) IRI-TEC in the September equinox; (g) GIM-TEC in the December solstice; and (h) IRI-TEC in the December solstice.

The figure illustrates that the spatial distribution characteristics, which change with the latitude and longitude exhibited by IRI-TEC and GIM-TEC, have good consistency. However, the equatorial ionospheric anomaly of IRI-TEC is more pronounced than that of GIM-TEC. The 2D correlation coefficients of the two types of TEC data are shown in Table 1. The correlation coefficients of the four seasons are at least 0.924.

Table 1 reveals that the mean biases between the season averages of global IRI-TEC and GIM-TEC at 12:00 UT are all negative. This result indicates that the TEC level predicted by the IRI-2016 model is lower than that of the GIM. This characteristic can also be seen in Fig. 1. The IRI-2016 model provides ionospheric parameters of up to 2000 km and is expected to be lower than the TEC up to GNSS satellites located at an altitude of approximately 20 000 km because of the missing plasmaspheric content. The mean bias and mean relative bias between IRI-TEC and GIM-TEC during the December solstice are larger than those in other seasons.

Considering the different levels of ionospheric activities at different latitudes, mean and RMS values of the discrepancies between seasonal averages of GIM-TEC and IRI-TEC over different latitudinal regions in 2013 were calculated. Results are shown in Fig. 2. From Fig. 2, the mean and RMS values over the area near the Equator generally exhibit peak values. GIM-TEC values over the Equator and low latitudes are much larger than IRI-TEC values, especially over the ionospheric trough near the magnetic Equator shown in Fig. 1. Due to high solar radiation in the equatorial region and Earth electric and magnetic field, the ionosphere over the equatorial region is at a high ionization level and its changes are complex. There are also anomalies such as an equatorial ionization anomaly (EIA) characterized by two low-latitude ionization crests of global maximum of plasma densities (Abdu, 2016). The IRI model has been reported to underestimate the ionospheric TEC at the equatorial station by Shreedevi et al. (2018), and a comparison of IRI-model-derived TEC and GPS TEC showed a wide departure, with ∼60 % deviation in their study. The mean and RMS values over the Southern Hemisphere during the December solstice are significantly large, and they are also very large over the Northern Hemisphere during the June solstice. Therefore, there are large discrepancies between GIM-TEC and IRI-TEC over the summer hemisphere. The large deviation of the ionospheric TEC estimated by the IRI model in the summer hemisphere indicates that the model cannot fully reflect the periodic seasonal variation in the ionosphere. As discussed by Li et al. (2016), solar activity component and periodic components are supposed to be the main reasons which account for the difference between the GIM TEC and the TEC from the IRI-2012 model. However, their conclusions are based on single-station time series data. In this article, we will further analyze the IRI model for spatiotemporal data.

Figure 2Mean and RMS values of the discrepancies between GIM-TEC and IRI-TEC at different latitudes during four seasons.

The gridded values of the global IRI-TEC and GIM-TEC at different UTs for each day of the year 2013 were used to calculate the daily RMS value. Results are shown in Fig. 3, which also displays the daily solar F10.7 index and daily average of geomagnetic AE index in 2013. The solar F10.7 and geomagnetic AE indexes are available at https://omniweb.gsfc.nasa.gov/form/dx1.html (last access: 21 April 2019).

Figure 3Daily (a) RMS of the differences between global IRI-TEC and GIM-TEC, (b) solar F10.7 index, and daily (c) average geomagnetic AE index in 2013.

Figure 3 demonstrates that the daily RMS of the differences between global IRI-TEC and GIM-TEC is in good agreement with the daily solar F10.7 index. The correlation coefficients between the RMS and the solar F10.7 or geomagnetic AE index are 0.78 and −0.19, respectively. Results indicate that the ionospheric TEC prediction error of the IRI-2016 model presents a strong correlation with solar activity.

Table 1Correlation coefficient and bias statistics among the season averages of global IRI-TEC and GIM-TEC at 12:00 UT in 2013.

## 3.2 Differences of spatial patterns between IRI-TEC and GIM-TEC based on the same time-varying characteristics

We combined the IRI-TEC and GIM-TEC data to obtain the same TEC time-varying characteristics using Eq. (6) and analyzed their differences in terms of spatial patterns.

The time-varying characteristics are reflected in the coefficient Ak of the EOF decomposition. Given that the TEC data are in accordance with the 2 h time interval, coefficient Ak is also the data that vary with the 2 h time interval. We described the coefficients of the base function according to the changes in UT and day of year (DOY) in Fig. 4 to reflect the seasonal changes effectively.

Figure 4Associated coefficients A1A6 of the first six orders of EOF base functions based on Eq. (6), and A1A6 plotted against UT and DOY.

The main EOF base functions extracted from Eq. (6) are shown in Fig. 5. The graphics in the left column of Fig. 5 exhibit the first six base functions Ei of GIM-TEC, whereas those in the right column of Fig. 5 depict the base functions of IRI-TEC.

Figure 5First six orders of EOF base functions E1E6 extracted on the basis of Eq. (6). The figures in the left column are the base functions of GIM-TEC, and those in the right column are the base functions of IRI-TEC.

The first base function E1 of GIM-TEC and IRI-TEC in Fig. 5a and b describe the overall average of global TEC. This function reflects the daily average effect of solar forcing and offset magnetic field (Talaat and Zhu, 2016). The TEC over the area near the geomagnetic equator exhibits a peak value. The TEC value decreases with the increase in geomagnetic latitude. The spatial distribution characteristics of E1 of the two models are very consistent. However, the peak GIM-TEC value over the geomagnetic equator is greater than that of the IRI-TEC. The ionospheric trough near the geomagnetic equator is evident in Fig. 5b. The daily mean A1 and solar F10.7 index are illustrated in Fig. 6, which shows that these two data sets demonstrate a consistent trend. The correlation coefficient between daily mean A1 and F10.7 index is 0.61. Solar activity is the primary determinant of the first base function E1.

Figure 6Daily mean first EOF coefficient A1 and daily solar F10.7 index.

Figure 5c–f show that the second and third base functions reflect the spatial distribution that varies in the longitude direction. The two base functions E2 and E3 approximately have the same magnitude and show a phase shift of π∕2, which is consistent with the results of Talaat and Zhu (2016). These functions reflect the change of diurnal solar radiation as it changes with the LT. This change of GIM-TEC and IRI-TEC is generally consistent; their main difference is reflected in the peak region of the Equator, and GIM-TEC shows large peak values. The EOF coefficients A2 and A3 corresponding to Fig. 4b and c show the change of the diurnal variation, and a change characteristic of the semiannual period is observed. The levels of A2 and A3 during equinox seasons are larger than those during solstice seasons.

The fourth base function E4 reflects interhemispheric asymmetry, which is mainly caused by the seasonal variation of the inclination angle of the Earth's orbit. A4 in Fig. 4d indicates the seasonal variation of the interhemispheric asymmetry of the TEC and a strong annual cycle. The TEC component corresponding to base function E4 in the Southern Hemisphere is positive. In the Northern Hemisphere, the maximum value of the E4 component is on DOY150, whereas that in the Southern Hemisphere is on DOY347.

Similar to E2 and E3, the fifth and sixth base functions E5 and E6 also reflect the spatial distribution characteristics along the longitude (Fig. 5i to l). In conjunction with Fig. 4e and f, these two base functions have semidiurnal period changes, and the phases of the two base functions differ by π∕4 and are of approximately equal magnitude. Base functions E5 and E6 represent a semidiurnal variation that changes with LT, and their coefficients A5 and A6 show a semiannual period. The intensity of the semidiurnal variation is strong during the equinox season and weak during the June solstice.

We calculated the variances, correlation coefficients, biases, and their relative biases to analyze the spatial distribution characteristics of GIM-TEC and IRI-TEC. The statistical results are shown in Table 2, which indicates that the base functions of the two data sets are correlated and present good consistency with Fig. 5.

Table 2Variances of the base function, correlation coefficient, and bias statistics among the base functions of GIM-TEC and IRI-TEC.

We showed the difference between the six base functions of GIM-TEC and IRI-TEC in Fig. 7 to have an intuitive understanding of the difference between the IRI and the GIM base functions.

Figure 7Differences of the first six orders of the base functions of GIM-TEC and IRI-TEC.

Figure 7 shows that the differences of other modes exhibit a large deviation in the equatorial and low-latitude regions, except for the interhemispheric asymmetry feature E4. The magnitudes of the spatial distribution changes of the IRI-TEC for all six base functions are significantly smaller than those of GIM-TEC.

The mean relative bias statistics of the base functions of GIM-TEC and IRI-TEC in Table 2 are negative. This finding indicates that the spatial variations of the base functions of IRI-TEC are generally underestimated compared with those of GIM-TEC. Here, the mean relative bias of E4 reached −56 %, and the underestimation is serious. This outcome is consistent with the statistical results in Table 1.

## 3.3 Differences of time-varying characteristics between IRI-TEC and GIM-TEC based on the same spatial patterns

Equation (7) shows that the same EOF base functions are extracted for GIM-TEC and IRI-TEC. The differences of the corresponding coefficients of the EOF base functions between GIM-TEC and IRI-TEC are then compared, and those of their time variation characteristics can be analyzed.

Figure 8 shows the six EOF base functions extracted in accordance with Eq. (7). Similar to the EOF base function extracted in Fig. 5, the first base function is consistent with the average variation of the TEC, varying with geomagnetic latitude. The second and third base functions are related to the diurnal variation of solar radiation change with longitude due to the LT. The fourth base function reflects the interhemispheric asymmetry caused by the seasonal variation of the inclination angle of the Earth's orbit. The fifth and sixth base functions reflect the characteristics of the semidiurnal variation with longitude due to the LT.

Figure 8Six EOF base functions E1E6 extracted in accordance with Eq. (7).

The coefficients of the different base functions of GIM-TEC and IRI-TEC obtained in accordance with Eq. (7) are shown in Fig. 9.

The time-varying characteristics of the coefficients in Fig. 9 are very consistent with the results shown in Fig. 4. From Fig. 9a and b, the variations of A1 are mainly related to solar activity, and solar activity is the primary determinant of the first base function E1 in Fig. 8a, which describes the overall average of global TEC. From Fig. 9c–f, the EOF coefficients A2 and A3 of GIM-TEC and IRI-TEC all obviously exhibit a diurnal period and a semiannual period. They reflect the diurnal variation of solar radiation change with longitude due to the LT. A4 values in Fig. 9g and h indicate a strong annual cycle variation of the interhemispheric asymmetry of the TEC. A5 and A6 show a semiannual period of the base functions E5 and E6, which represent a longitudinal variation that changes with LT. The EOF coefficients of GIM-TEC and IRI-TEC have consistent annual, semiannual, diurnal, and semidiurnal variations. Therefore, Fig. 9 shows that GIM-TEC and IRI-TEC have highly consistent temporal variation characteristics based on the same spatial distribution modes Ek according to Eq. (7). The variance and correlation coefficients of A1A6 of the two types of data and the bias statistics of such coefficients are shown in Table 3.

Figure 9Associated coefficients A1A6 of the EOF base functions extracted in accordance with Eq. (7).

Table 3Variances of base function, correlation coefficient, and bias statistics among coefficients A1A6 of GIM-TEC and IRI-TEC.

The magnitudes of coefficients A1A6 of IRI-TEC are generally smaller than those of the GIM-TEC, especially for A4. The maximum and minimum values of GIM-TEC A4 in Fig. 9g are 302.27 and −431.47, respectively. The variation range of the IRI-TEC A4 in Fig. 9h is −138.99 to 165.13. Results in Table 3 indicate that A4 exhibits the largest mean relative bias.

Figure 9 shows that A1A6 reflect the time-varying characteristics of different scales. We conducted EOF decomposition on A1A6 according to the following equation to divide their diurnal and seasonal variation characteristics:

$\begin{array}{}\text{(13)}& {A}_{i}\left(\mathrm{UT},\mathrm{DOY}\right)=\sum _{k=\mathrm{1}}^{N}{E}_{ik}\left(\mathrm{UT}\right)×{A}_{ik}\left(\mathrm{DOY}\right),\end{array}$

where Ai represents the coefficient of the ith order the EOF base function. This part is the second-layer EOF decomposition in this study.

Equation (13) shows that the time-varying feature Eik depending on UT and seasonal variation Aik can be obtained. According to the percentage variance of the second-layer EOF decomposition, the first EOF component has already explained more than 99 % of the total variance of Ai. Therefore, the first EOF component is the most significant, and we will only present the first-order result of the second-layer EOF decomposition in this study. The decomposed first base function Ei1 and associated coefficient Ai1 are shown in Fig. 10.

Figure 10First base function Ei1 and associated coefficient Ai1 of the six coefficients A1A6 according to Eq. (12). The monthly smoothed Ai1 of GIM-TEC and daily solar F10.7 index are shown together with Ai1.

The left column of Fig. 10 shows base function Eik, which represents the diurnal variation characteristic of the base function Ei. The coefficients of the second-layer EOF decomposition Ai1 represent the variations in long timescales. Ai1 is shown in the right column of Fig. 10. Previous studies have shown that the long timescale variations of TEC are mainly influenced by solar and geomagnetic activities and periodical variation. The solar F10.7 index is also shown in the right column of Fig. 10 together with Ai1.

The first base function E1 in Fig. 8a describes the overall average global TEC, and Fig. 10a shows E11, the diurnal variation characteristic of E1. GIM-TEC and IRI-TEC have similar magnitudes, whereas the diurnal variation of IRI-TEC is insignificant. A11 of GIM-TEC and IRI-TEC in Fig. 10b shows a pronounced semiannual period. However, A11 values of GIM-TEC on most days are larger than those of IRI-TEC, and the correlation between the F10.7 index and A11 of GIM-TEC is evidently observed.

As shown in Fig. 10c, e, and g, the diurnal variations of the second, third, and fourth base functions E2E4 of GIM-TEC and IRI-TEC show minimal discrepancy. Hence, the IRI-2016 model accurately captures the diurnal variations of the solar radiation according to LT and interhemispheric asymmetry.

A21 and A31 of GIM-TEC and IRI-TEC are shown in Fig. 10d and f. These functions evidently demonstrate a semidiurnal variation period. A21 and A31 of IRI-TEC during the equinox season are lower than those of GIM-TEC. The correlation between the F10.7 index and A21 and A31 of GIM-TEC is also observed. A41 of GIM-TEC and A41 of IRI-TEC in Fig. 10h exhibit an evident annual period variation of interhemispheric asymmetry. However, the summer-to-winter annual variation of GIM-TEC is much larger than that of IRI-TEC.

The fifth and sixth base functions E5 and E6 in Fig. 8e and f reflect the spatial distribution characteristics along the longitude due to LT. E51 and E61 in Fig. 10i and k represent a semidiurnal variation. However, shifts in the peak value time between GIM-TEC and IRI-TEC are detected in E51 and E61. A51 and A61 in Fig. 10j and l exhibit a semiannual variation, and A51 and A61 of GIM-TEC are relatively consistent with those of IRI-TEC.

We calculated the correlation coefficients between Ai1 of GIM-TEC and the solar F10.7 index. Results are shown in Table 4. Coefficients A11, A21, and A31 are highly related to solar activity.

Table 4Correlation coefficients between Ai1 of GIM-TEC and solar F10.7 index.

A11A61 in Fig. 10 show that IRI-TEC mainly reflects the annual and semiannual variations of the ionospheric TEC. The monthly and short-term variations with solar activity are unrepresented by IRI-TEC.

Although the IRI-TEC will be smaller than the GIM-TEC because of the missing plasmaspheric content, A11 of IRI-TEC in Fig. 10b shows quite a large underestimation compared with that of GIM-TEC. The strong correlation between A11 of GIM-TEC and solar activity is unrepresented by A11 of IRI-TEC. The diurnal variation of the first base function of GIM-TEC represented by E11 is partially represented by E11 of IRI-TEC. The variance contribution rate of the first EOF component reaches 79.03 %; thus, the influence of its coefficient is large for the deviation of IRI-TEC and GIM-TEC.

4 Conclusion

In this study, the global TEC prediction performance of the IRI-2016 model was evaluated. The EOF decomposition method was introduced to compare the global TEC data from the IRI-2016 model and GIMs in 2013. The prediction performance of the IRI-2016 model could be evaluated from two perspectives, namely, spatial pattern and temporal variation. The main conclusions are as follows:

1. A general underestimation of the IRI-2016 model can be observed compared with the season averages of global GIM-TEC in 2013, and the RMS of the global TEC deviation is strongly correlated with the solar activity F10.7 index.

2. The six base functions extracted by performing EOF decomposition on the global TEC data from IRI-2016 and GIMs include the following: the variation with the geomagnetic latitude reflecting the daily averaged solar forcing, the diurnal and semidiurnal periodic changes with longitude due to local time, and the interhemispheric asymmetry caused by the annual variation of the inclination angle of the Earth's orbit. The spatiotemporal features extracted from IRI-TEC and GIM-TEC data have good consistency. The IRI-2016 model follows the variation patterns of the observed GIM-TEC.

3. The spatial variation characteristics of IRI-TEC and GIM-TEC can be extracted for comparison on the basis of the same EOF coefficients. Results show that the spatial distribution fluctuation of the IRI-TEC is smaller than that of GIM-TEC. The average relative deviation of the base function representing the interhemispheric asymmetry reaches −56.7 %. The interhemispheric asymmetry presents a relatively stable deviation between IRI-TEC and GIM-TEC. The other spatial distribution variations have large deviations at the Equator and low latitudes.

4. The temporal variation characteristics of IRI-TEC and GIM-TEC are extracted and compared on the basis of the same EOF base functions. The degree of IRI-TEC changes with time is weaker than that of GIM-TEC. The average relative deviation of the fourth base function coefficient reaches −52.83 %. Most diurnal, annual, and semiannual variations of the six base functions of IRI-TEC are consistent with those of GIM-TEC. However, the change with solar activity is unrepresented by IRI-TEC. The diurnal variation of the first base function and the annual variation of the fourth base function have a relatively large deviation between IRI-TEC and GIM-TEC.

5. Results of the spatial and temporal variation characteristic analyses show that the deviation of the first and fourth EOF components between IRI-TEC and GIM-TEC are the two main influencing factors.

Data availability
Data availability.

The data used in this study were downloaded from ftp://cddis.gsfc.nasa.gov/ (CDDIS, 2019), http://irimodel.org/ (COSPAR and URSI, 2019), and https://omniweb.gsfc.nasa.gov/ (CCMC, 2019).

Author contributions
Author contributions.

SL contributed to the conception of the study. SL and JX contributed significantly to the data analysis and manuscript preparation. HZ and JZ performed the model validation and wrote part of the manuscript. ZX and MX contributed to some data analysis work.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Financial support
Financial support.

This research has been supported by the Fundamental Research Funds for the Central Universities (grant no. 2652017105) and the National Natural Science Foundation of China (grant no. 41574011).

Review statement
Review statement.

This paper was edited by Ana G. Elias and reviewed by two anonymous referees.

References

A, E., Zhang, D.-H., Xiao, Z., Hao, Y.-Q., Ridley, A. J., and Moldwin, M.: Modeling ionospheric foF2 by using empirical orthogonal function analysis, Ann. Geophys., 29, 1501–1515, https://doi.org/10.5194/angeo-29-1501-2011, 2011.

Abdu, M. A.: Electrodynamics of ionospheric weather over low latitudes, Abdu Geosci. Lett., 3, 11, https://doi.org/10.1186/s40562-016-0043-6, 2016.

Altadill, D., Magdaleno, S., Torta, J. M., and Blanch, E.: Global empirical models of the density peak height and of the equivalent scale height for quiet conditions, Adv. Space Res., 52, 1756–1769, https://doi.org/10.1016/j.asr.2012.11.018, 2013.

Andima, G., Amabayo, E. B., Jurua, E., and Cilliers, P. J.: Modeling of GPS total electron content over the African low-latitude region using empirical orthogonal functions, Ann. Geophys., 37, 65–76, https://doi.org/10.5194/angeo-37-65-2019, 2019.

Atici, R.: Comparison of GPS TEC with modelled values from IRI 2016 and IRI-PLAS over Istanbul, Turkey, Astrophys. Space Sci., 363, 231, https://doi.org/10.1007/s10509-018-3457-0, 2018.

Bilitza, D.: International Reference Ionosphere 2000, Radio Sci., 36, 261–275, https://doi.org/10.1029/2000RS002432, 2001.

Bilitza, D.: The International Reference Ionosphere-Status 2013, Adv. Space Res., 5, 1914–1927, https://doi.org/10.1016/j.asr.2014.07.032, 2015.

Bilitza, D. and Reinisch, B. W.: International Reference Ionosphere 2007: Improvements and new parameters, Adv. Space Res., 42, 599–609, https://doi.org/10.1016/j.asr.2007.07.048, 2008.

Bilitza, D., Watanabe, S., Truhlik, V., and Altadill, D.: IRI-2016: Description and Introduction, 41st COSPAR Scientific Assembly, July, Istanbul, Turkey, 2016.

Bilitza, D., Altadill, D., Truhlik, V., Shubin, V., Galkin, I., Reinisch, B., and Huang, X.: International Reference Ionosphere 2016: from ionospheric climate to real-time weather predictions, Space Weather, 15, 418–429, 2017.

Bouya, Z., Terkildsen, M., Francis, M., and Neudegg, D.: EOF Analysis applied to Australian Regional Ionospheric TEC modelling, AOSWA, 22–24 February 2012, Chiang Mai, Thailand, 2012.

CCMC: Solar activity data, available at: https://omniweb.gsfc.nasa.gov/, last access: 21 April 2019.

CDDIS: GIM-TEC data, available at: ftp://cddis.gsfc.nasa.gov/, last access: 10 April 2019.

COSPAR and URSI: Fortran source code of IRI-2016, available at: http://irimodel.org/, last access: 12 March 2019.

Chang, X., Zou, B., Guo, J., Zhu, G., Li, W., and Li, W.: One sliding PCA method to detect ionospheric anomalies before strong Earthquakes: Cases study of Qinghai, Honshu, Hotan and Nepal earthquakes, Adv. Space Res., 59, 2058–2070, https://doi.org/10.1016/j.asr.2017.02.007, 2017.

Chauhan, V. and Singh, O. P.: A morphological study of GPS-TEC data at Agra and their comparison with the IRI model, Adv. Space Res., 46, 280–290, https://doi.org/10.1016/j.asr.2010.03.018, 2010.

Dabbakuti, J. R. K. K. and Ratnam, D. V.: Characterization of ionospheric variability in TEC using EOF and wavelets over low-latitude GNSS stations, Adv. Space Res., 57, 2427–2443, https://doi.org/10.1016/j.asr.2016.03.029, 2016.

Dabbakuti, J. R. K. K. and Ratnam, D. V.: Modeling and analysis of GPS-TEC low latitude climatology during the 24th solar cycle using empirical orthogonal functions, Adv. Space Res., 60, 1751–1764, https://doi.org/10.1016/j.asr.2017.06.048, 2017.

Feltens, J., Angling, M., Jackson-Booth, N., Jakowski, N., Hoque, M., HernándezPajares, M., Aragón, Á., María, Á., and Orús-Pérez, R.: Comparative testing of four ionospheric models driven with GPS measurements, Radio Sci., 46, RS0D12, https://doi.org/10.1029/2010RS004584, 2011.

Jiang, H., Liu, J., Wang, Z., An, J., Ou, J., Liu, S., and Wang, N.: Assessment of spatial and temporal TEC variations derived from ionospheric models over the polar regions, J. Geodesy, 93, 455–471, https://doi.org/10.1007/s00190-018-1175-6, 2019.

Kenpankho, P., Supnithi, P., and Nagatsuma, T.: Comparison of observed TEC values with IRI-2007 TEC and IRI-2007 TEC with optional foF2 measurements predictions at an equatorial region, Chumphon, Thailand, Adv. Space Res., 52, 1820–1826, https://doi.org/10.1016/j.asr.2013.08.012, 2013.

Li, S., Li, L., and Peng, J.: Variability of Ionospheric TEC and the Performance of the IRI-2012 Model at the BJFS Station, China, Acta Geophys., 64, 1970–1987, 2016.

Li, S., Zhou, H., Xu, J., Wang, Z., Li, L., and Zheng, Y.: Modeling and analysis of ionosphere TEC over China and adjacent areas based on EOF method, Adv. Space Res., 64, 400–414, https://doi.org/10.1016/j.asr.2019.04.018, 2019.

Maltseva, O. A., Mozhaeva, N. S., Poltavsky, O. S., and Zhbankov, G. A.: Use of TEC global maps and the IRI model to study ionospheric response to geomagnetic disturbances, Adv. Space Res., 49, 1076–1087, https://doi.org/10.1016/j.asr.2012.01.005, 2012.

Mao, T., Wan, W., Yue, X., Sun, L., Zhao, B., and Guo, J.: An empirical orthogonal function model of total electron content over China, Radio Sci., 43, RS2009, https://doi.org/10.1029/2007RS003629, 2008.

Okoh, D., McKinnell, L., Cilliers, P., and Okeke, P.: Using GPS-TEC data to calibrate VTEC computed with the IRI model over Nigeria, Adv. Space Res., 52, 1791–1797, https://doi.org/10.1016/j.asr.2012.11.013, 2013.

Pearson, K.: On lines and planes of closest fit to systems of points in space, Philos. Mag., 2, 559–572, 1901.

Rawer, K., Ramakrishnan, S., and Bilitza, D.: International Reference Ionosphere 1978. International Union of Radio Science, URSI Special Report, 75 pp., Bruxelles, Belgium, 1978.

Scidá, L. A., Ezquer, R. G., Cabrera, M. A., Mosert, M., Brunini, C., and Buresova, D.: On the IRI 2007 performance as a TEC predictor for the South American sector, J. Atmos. Sol.-Terr. Phy., 81–82, 50–58, https://doi.org/10.1016/j.jastp.2012.04.001, 2012.

Sharma, S. K., Ansari, K., and Panda, S. K.: Analysis of Ionospheric TEC Variation over Manama, Bahrain, and Comparison with IRI-2012 and IRI-2016 Models, Arab. J. Sci. Eng., 43, 3823–3830, 2018.

Shreedevi, P. R., Choudhary, R. K., Yadav, S., Thampi, S., and Ajesh, A.: Variation of the TEC at a dip equatorial station, Trivandrum and a mid latitude station, Hanle during the descending phase of the solar cycle 24 (2014–2016), J. Atmos. Sol.-Terr. Phy., 179, 425–434, 2018.

Shubin, V. N.: Global median model of the F2-layer peak height based on ionospheric radio-occultation and ground-based Digisonde observations, Adv. Space Res., 56, 916–928, https://doi.org/10.1016/j.asr.2015.05.029, 2015.

Talaat, E. R. and Zhu, X.: Spatial and temporal variation of total electron content as revealed by principal component analysis, Ann. Geophys., 34, 1109–1117, https://doi.org/10.5194/angeo-34-1109-2016, 2016.

Tariku, Y. A.: Assessment of improvement of the IRI model over Ethiopia for the modeling of the variability of TEC during the period 2013–2016, Adv. Space Res., 63, 1634–1645, https://doi.org/10.1016/j.asr.2018.11.014, 2018.

Uwamahoro, J. and Habarulema, J. B.: Modelling total electron content during geomagnetic storm conditions using empirical orthogonal functions and neural networks, J. Geophys. Res.-Space, 120, 11000–11012, https://doi.org/10.1002/2015JA021961, 2015.

Yao, Y., Liu, L., Kong, J., and Zhai, C.: Global Ionospheric Modeling Based on Multi-GNSS, Satellite Altimetry and Formosat-3/COSMIC Data, GPS Solut., 22, 104, https://doi.org/10.1007/s10291-018-0770-6, 2018.

Zakharenkova, I. E., Cherniak, Iu. V., Krankowski, A., and Shagimuratov, I. I.: Vertical TEC representation by IRI 2012 and IRI Plas models for European mid latitudes, Adv. Space Res., 55, 2070–2076, https://doi.org/10.1016/j.asr.2014.07.027, 2015.

Zhang, S., Foster, J. C., Coster, A. J., and Erickson, P. J.: East–West Coast differences in total electron content over the continental US, Geophys. Res. Lett., 38, L19101, https://doi.org/10.1029/2011GL049116, 2011.

Zhang, S., Chen, Z., Coster, A. J., Erickson, P. J., and Foster, J. C.: Ionospheric symmetry caused by geomagnetic declination over North America, Geophys. Res. Lett., 40, 5350–5354, https://doi.org/10.1002/2013GL057933, 2013.

Zhao, B., Wan, W., Liu, L., Yue, X., and Venkatraman, S.: Statistical characteristics of the total ion density in the topside ionosphere during the period 1996–2004 using empirical orthogonal function (EOF) analysis, Ann. Geophys., 23, 3615–3631, https://doi.org/10.5194/angeo-23-3615-2005, 2005.