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

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 E_i, 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 A_k 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}\times {\displaystyle \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}(\mathrm{Lat},\mathrm{Lon},\mathrm{UT},\mathrm{DOY})$ 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}\times \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}\times \mathrm{365}=\mathrm{4380}$. After performing EOF decomposition, the base function E_k(grid) expressing a spatial pattern and the associated coefficient A_k(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}}{X}_{\mathrm{IRI}}\right]& =\sum _{k=\mathrm{1}}^N{E}_k\cdot \left[A_{k,\mathrm{GIM}}{A}_{k,\mathrm{IRI}}\right]\\ & =\left[\sum _{k=\mathrm{1}}^N{E}_k\cdot A_{k,\mathrm{GIM}}\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}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(Y_i-Y_i^{\prime }),\end{array}$$

where n is the total number of sample data, and Y_i 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{\%}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\displaystyle \frac{(Y_i-Y_i^{\prime })}{Y_i^{\prime }}}\times \mathrm{100}.\end{array}$$

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

$$\begin{array}{}\text{(10)}& \mathrm{RMS}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _i^n(Y_i-Y_i^{\prime }{)}^{\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(A_{mn}-\stackrel{\mathrm{‾}}{A})(B_{mn}-\stackrel{\mathrm{‾}}{B})\cdot {\left[\sum _{m=\mathrm{1}}^M\sum _{n=\mathrm{1}}^N(A_{mn}-\stackrel{\mathrm{‾}}{A}{)}^{\mathrm{2}}\right]}^{-\frac{\mathrm{1}}{\mathrm{2}}}\\ & \cdot {\left[\sum _{m=\mathrm{1}}^M\sum _{n=\mathrm{1}}^N(B_{mn}-\stackrel{\mathrm{‾}}{B}{)}^{\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}={\displaystyle \frac{\mathrm{1}}{MN}}\sum _{m=\mathrm{1}}^M\sum _{n=\mathrm{1}}^N{A}_{mn};\mathrm{and}\stackrel{\mathrm{‾}}{B}={\displaystyle \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 1 Season 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 2 Mean 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 3 Daily (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 1 Correlation 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 A_k of the EOF decomposition. Given that the TEC data are in accordance with the 2 h time interval, coefficient A_k 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 4 Associated coefficients A₁–A₆ of the first six orders of EOF base functions based on Eq. (6), and A₁–A₆ 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 E_i of GIM-TEC, whereas those in the right column of Fig. 5 depict the base functions of IRI-TEC.

Figure 5 First six orders of EOF base functions E₁–E₆ 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 E₁ 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 E₁ 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 A₁ 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 A₁ and F10.7 index is 0.61. Solar activity is the primary determinant of the first base function E₁.

Figure 6 Daily mean first EOF coefficient A₁ 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 E₂ and E₃ 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 A₂ and A₃ 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 A₂ and A₃ during equinox seasons are larger than those during solstice seasons.

The fourth base function E₄ reflects interhemispheric asymmetry, which is mainly caused by the seasonal variation of the inclination angle of the Earth's orbit. A₄ 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 E₄ in the Southern Hemisphere is positive. In the Northern Hemisphere, the maximum value of the E₄ component is on DOY150, whereas that in the Southern Hemisphere is on DOY347.

Similar to E₂ and E₃, the fifth and sixth base functions E₅ and E₆ 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 E₅ and E₆ represent a semidiurnal variation that changes with LT, and their coefficients A₅ and A₆ 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 2 Variances 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 7 Differences 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 E₄. 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 E₄ 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 8 Six EOF base functions E₁–E₆ 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 A₁ are mainly related to solar activity, and solar activity is the primary determinant of the first base function E₁ in Fig. 8a, which describes the overall average of global TEC. From Fig. 9c–f, the EOF coefficients A₂ and A₃ 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. A₄ values in Fig. 9g and h indicate a strong annual cycle variation of the interhemispheric asymmetry of the TEC. A₅ and A₆ show a semiannual period of the base functions E₅ and E₆, 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 E_k according to Eq. (7). The variance and correlation coefficients of A₁–A₆ of the two types of data and the bias statistics of such coefficients are shown in Table 3.

Figure 9 Associated coefficients A₁–A₆ of the EOF base functions extracted in accordance with Eq. (7).

Table 3 Variances of base function, correlation coefficient, and bias statistics among coefficients A₁–A₆ of GIM-TEC and IRI-TEC.

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

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

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

where A_i 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 E_ik depending on UT and seasonal variation A_ik 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 A_i. 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 E_i1 and associated coefficient A_i1 are shown in Fig. 10.

Figure 10 First base function E_i1 and associated coefficient A_i1 of the six coefficients A₁–A₆ according to Eq. (12). The monthly smoothed A_i1 of GIM-TEC and daily solar F10.7 index are shown together with A_i1.

The left column of Fig. 10 shows base function E_ik, which represents the diurnal variation characteristic of the base function E_i. The coefficients of the second-layer EOF decomposition A_i1 represent the variations in long timescales. A_i1 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 A_i1.

The first base function E₁ in Fig. 8a describes the overall average global TEC, and Fig. 10a shows E₁₁, the diurnal variation characteristic of E₁. GIM-TEC and IRI-TEC have similar magnitudes, whereas the diurnal variation of IRI-TEC is insignificant. A₁₁ of GIM-TEC and IRI-TEC in Fig. 10b shows a pronounced semiannual period. However, A₁₁ values of GIM-TEC on most days are larger than those of IRI-TEC, and the correlation between the F10.7 index and A₁₁ 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 E₂–E₄ 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.

A₂₁ and A₃₁ of GIM-TEC and IRI-TEC are shown in Fig. 10d and f. These functions evidently demonstrate a semidiurnal variation period. A₂₁ and A₃₁ of IRI-TEC during the equinox season are lower than those of GIM-TEC. The correlation between the F10.7 index and A₂₁ and A₃₁ of GIM-TEC is also observed. A₄₁ of GIM-TEC and A₄₁ 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 E₅ and E₆ in Fig. 8e and f reflect the spatial distribution characteristics along the longitude due to LT. E₅₁ and E₆₁ 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 E₅₁ and E₆₁. A₅₁ and A₆₁ in Fig. 10j and l exhibit a semiannual variation, and A₅₁ and A₆₁ of GIM-TEC are relatively consistent with those of IRI-TEC.

We calculated the correlation coefficients between A_i1 of GIM-TEC and the solar F10.7 index. Results are shown in Table 4. Coefficients A₁₁, A₂₁, and A₃₁ are highly related to solar activity.

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

A₁₁–A₆₁ 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, A₁₁ of IRI-TEC in Fig. 10b shows quite a large underestimation compared with that of GIM-TEC. The strong correlation between A₁₁ of GIM-TEC and solar activity is unrepresented by A₁₁ of IRI-TEC. The diurnal variation of the first base function of GIM-TEC represented by E₁₁ is partially represented by E₁₁ 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.