Regular paper 16 Mar 2020
Regular paper  16 Mar 2020
Global total electron content prediction performance assessment of the IRI2016 model based on empirical orthogonal function decomposition
 ^{1}School of Land Science and Technology, China University of Geosciences (Beijing), Beijing 100083, China
 ^{2}Shanxi Key Laboratory of Resources, Environment and Disaster Monitoring, Jinzhong 221116, China
 ^{1}School of Land Science and Technology, China University of Geosciences (Beijing), Beijing 100083, China
 ^{2}Shanxi Key Laboratory of Resources, Environment and Disaster Monitoring, Jinzhong 221116, China
Correspondence: Shuhui Li (li.shuhui@163.com)
Hide author detailsCorrespondence: Shuhui Li (li.shuhui@163.com)
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 (IRI2016) model in 2013. Results showed that the main spatial patterns and timevarying 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 IRI2016 and GIM TECs were analyzed by extracting the same timevarying 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 timevarying characteristics between the IRI2016 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 IRI2016 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 IRI2016 and GIM TECs.
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 highfrequency 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 ionosphericbased 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 IRI2016 version (Bilitza et al., 2017). The model has been improved by ingesting all available data from worldwide groundbased and satellite observations to enhance the model capacity. IRI2016 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 GNSSderived 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 IRI2016 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 IRI2016. 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 IRI2016 and the GIM TECs at a global scale.
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 GIMTEC hereafter).
2.2 IRI2016
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 IRI2016 (Bilitza et al., 2016, 2017). After IRI2012, IRI2016 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 IRI2016 can be downloaded from http://irimodel.org/ (last access: 12 March 2019). The IRI software package contains FORTRAN subroutines, model coefficients, index files for IRI2016 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 IRI2016 will be called IRITEC hereafter. IRITEC 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
The EOF base functions E_{i}, with i=1, 2, 3, …, N, are the eigenvectors of the covariance matrix and obtained by solving
where λ_{i} is the associated eigenvalues. Once the EOF base functions are known, the EOF coefficients A_{k} are obtained using
The original data set X can be decomposed in terms of the EOF base functions and associated coefficients in accordance with
The percentage of the total variance in the data set accounted for by the ith EOF component is given as follows:
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 1year 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 timevarying features, can be obtained by arranging IRITEC and GIMTEC according to the same number of columns. Accordingly, comparing the two data sets' spatial variation features represented by the base functions is feasible.
If IRITEC and GIMTEC 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.
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:
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 IRI2016 and GIMs or the values of base functions or coefficients of base functions. The mean relative bias (Bias_rel) can be calculated as follows:
The RMS error of the bias can be calculated using the following expression:
The 2D linear correlation coefficient was used to investigate the similarity of the spatial pattern of IRITEC and GIMTEC. The 2D linear correlation coefficient ρ for two matrices A and B with M×N dimension is calculated as
where $\stackrel{\mathrm{\u203e}}{A}$ and $\stackrel{\mathrm{\u203e}}{B}$ are the mean values of matrices A and B, respectively, and they are written as
3.1 GIMTEC and IRITEC in 2013
Figure 1 shows the season averages of global GIMTEC and IRITEC 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.
The figure illustrates that the spatial distribution characteristics, which change with the latitude and longitude exhibited by IRITEC and GIMTEC, have good consistency. However, the equatorial ionospheric anomaly of IRITEC is more pronounced than that of GIMTEC. 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 IRITEC and GIMTEC at 12:00 UT are all negative. This result indicates that the TEC level predicted by the IRI2016 model is lower than that of the GIM. This characteristic can also be seen in Fig. 1. The IRI2016 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 IRITEC and GIMTEC 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 GIMTEC and IRITEC 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. GIMTEC values over the Equator and low latitudes are much larger than IRITEC 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 lowlatitude 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 IRImodelderived 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 GIMTEC and IRITEC 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 IRI2012 model. However, their conclusions are based on singlestation time series data. In this article, we will further analyze the IRI model for spatiotemporal data.
The gridded values of the global IRITEC and GIMTEC 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 demonstrates that the daily RMS of the differences between global IRITEC and GIMTEC 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 IRI2016 model presents a strong correlation with solar activity.
3.2 Differences of spatial patterns between IRITEC and GIMTEC based on the same timevarying characteristics
We combined the IRITEC and GIMTEC data to obtain the same TEC timevarying characteristics using Eq. (6) and analyzed their differences in terms of spatial patterns.
The timevarying 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.
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 GIMTEC, whereas those in the right column of Fig. 5 depict the base functions of IRITEC.
The first base function E_{1} of GIMTEC and IRITEC 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_{1} of the two models are very consistent. However, the peak GIMTEC value over the geomagnetic equator is greater than that of the IRITEC. The ionospheric trough near the geomagnetic equator is evident in Fig. 5b. The daily mean A_{1} 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_{1} and F10.7 index is 0.61. Solar activity is the primary determinant of the first base function E_{1}.
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_{2} and E_{3} 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 GIMTEC and IRITEC is generally consistent; their main difference is reflected in the peak region of the Equator, and GIMTEC shows large peak values. The EOF coefficients A_{2} and A_{3} 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_{2} and A_{3} during equinox seasons are larger than those during solstice seasons.
The fourth base function E_{4} reflects interhemispheric asymmetry, which is mainly caused by the seasonal variation of the inclination angle of the Earth's orbit. A_{4} 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_{4} in the Southern Hemisphere is positive. In the Northern Hemisphere, the maximum value of the E_{4} component is on DOY150, whereas that in the Southern Hemisphere is on DOY347.
Similar to E_{2} and E_{3}, the fifth and sixth base functions E_{5} and E_{6} 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_{5} and E_{6} represent a semidiurnal variation that changes with LT, and their coefficients A_{5} and A_{6} 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 GIMTEC and IRITEC. 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.
We showed the difference between the six base functions of GIMTEC and IRITEC in Fig. 7 to have an intuitive understanding of the difference between the IRI and the GIM base functions.
Figure 7 shows that the differences of other modes exhibit a large deviation in the equatorial and lowlatitude regions, except for the interhemispheric asymmetry feature E_{4}. The magnitudes of the spatial distribution changes of the IRITEC for all six base functions are significantly smaller than those of GIMTEC.
The mean relative bias statistics of the base functions of GIMTEC and IRITEC in Table 2 are negative. This finding indicates that the spatial variations of the base functions of IRITEC are generally underestimated compared with those of GIMTEC. Here, the mean relative bias of E_{4} reached −56 %, and the underestimation is serious. This outcome is consistent with the statistical results in Table 1.
3.3 Differences of timevarying characteristics between IRITEC and GIMTEC based on the same spatial patterns
Equation (7) shows that the same EOF base functions are extracted for GIMTEC and IRITEC. The differences of the corresponding coefficients of the EOF base functions between GIMTEC and IRITEC 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.
The coefficients of the different base functions of GIMTEC and IRITEC obtained in accordance with Eq. (7) are shown in Fig. 9.
The timevarying 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_{1} are mainly related to solar activity, and solar activity is the primary determinant of the first base function E_{1} in Fig. 8a, which describes the overall average of global TEC. From Fig. 9c–f, the EOF coefficients A_{2} and A_{3} of GIMTEC and IRITEC 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_{4} values in Fig. 9g and h indicate a strong annual cycle variation of the interhemispheric asymmetry of the TEC. A_{5} and A_{6} show a semiannual period of the base functions E_{5} and E_{6}, which represent a longitudinal variation that changes with LT. The EOF coefficients of GIMTEC and IRITEC have consistent annual, semiannual, diurnal, and semidiurnal variations. Therefore, Fig. 9 shows that GIMTEC and IRITEC 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_{1}–A_{6} of the two types of data and the bias statistics of such coefficients are shown in Table 3.
The magnitudes of coefficients A_{1}–A_{6} of IRITEC are generally smaller than those of the GIMTEC, especially for A_{4}. The maximum and minimum values of GIMTEC A_{4} in Fig. 9g are 302.27 and −431.47, respectively. The variation range of the IRITEC A_{4} in Fig. 9h is −138.99 to 165.13. Results in Table 3 indicate that A_{4} exhibits the largest mean relative bias.
Figure 9 shows that A_{1}–A_{6} reflect the timevarying characteristics of different scales. We conducted EOF decomposition on A_{1}–A_{6} according to the following equation to divide their diurnal and seasonal variation characteristics:
where A_{i} represents the coefficient of the ith order the EOF base function. This part is the secondlayer EOF decomposition in this study.
Equation (13) shows that the timevarying feature E_{ik} depending on UT and seasonal variation A_{ik} can be obtained. According to the percentage variance of the secondlayer 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 firstorder result of the secondlayer EOF decomposition in this study. The decomposed first base function E_{i1} and associated coefficient A_{i1} are shown in Fig. 10.
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 secondlayer 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_{1} in Fig. 8a describes the overall average global TEC, and Fig. 10a shows E_{11}, the diurnal variation characteristic of E_{1}. GIMTEC and IRITEC have similar magnitudes, whereas the diurnal variation of IRITEC is insignificant. A_{11} of GIMTEC and IRITEC in Fig. 10b shows a pronounced semiannual period. However, A_{11} values of GIMTEC on most days are larger than those of IRITEC, and the correlation between the F10.7 index and A_{11} of GIMTEC is evidently observed.
As shown in Fig. 10c, e, and g, the diurnal variations of the second, third, and fourth base functions E_{2}–E_{4} of GIMTEC and IRITEC show minimal discrepancy. Hence, the IRI2016 model accurately captures the diurnal variations of the solar radiation according to LT and interhemispheric asymmetry.
A_{21} and A_{31} of GIMTEC and IRITEC are shown in Fig. 10d and f. These functions evidently demonstrate a semidiurnal variation period. A_{21} and A_{31} of IRITEC during the equinox season are lower than those of GIMTEC. The correlation between the F10.7 index and A_{21} and A_{31} of GIMTEC is also observed. A_{41} of GIMTEC and A_{41} of IRITEC in Fig. 10h exhibit an evident annual period variation of interhemispheric asymmetry. However, the summertowinter annual variation of GIMTEC is much larger than that of IRITEC.
The fifth and sixth base functions E_{5} and E_{6} in Fig. 8e and f reflect the spatial distribution characteristics along the longitude due to LT. E_{51} and E_{61} in Fig. 10i and k represent a semidiurnal variation. However, shifts in the peak value time between GIMTEC and IRITEC are detected in E_{51} and E_{61}. A_{51} and A_{61} in Fig. 10j and l exhibit a semiannual variation, and A_{51} and A_{61} of GIMTEC are relatively consistent with those of IRITEC.
We calculated the correlation coefficients between A_{i1} of GIMTEC and the solar F10.7 index. Results are shown in Table 4. Coefficients A_{11}, A_{21}, and A_{31} are highly related to solar activity.
A_{11}–A_{61} in Fig. 10 show that IRITEC mainly reflects the annual and semiannual variations of the ionospheric TEC. The monthly and shortterm variations with solar activity are unrepresented by IRITEC.
Although the IRITEC will be smaller than the GIMTEC because of the missing plasmaspheric content, A_{11} of IRITEC in Fig. 10b shows quite a large underestimation compared with that of GIMTEC. The strong correlation between A_{11} of GIMTEC and solar activity is unrepresented by A_{11} of IRITEC. The diurnal variation of the first base function of GIMTEC represented by E_{11} is partially represented by E_{11} of IRITEC. The variance contribution rate of the first EOF component reaches 79.03 %; thus, the influence of its coefficient is large for the deviation of IRITEC and GIMTEC.
In this study, the global TEC prediction performance of the IRI2016 model was evaluated. The EOF decomposition method was introduced to compare the global TEC data from the IRI2016 model and GIMs in 2013. The prediction performance of the IRI2016 model could be evaluated from two perspectives, namely, spatial pattern and temporal variation. The main conclusions are as follows:

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

The six base functions extracted by performing EOF decomposition on the global TEC data from IRI2016 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 IRITEC and GIMTEC data have good consistency. The IRI2016 model follows the variation patterns of the observed GIMTEC.

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

The temporal variation characteristics of IRITEC and GIMTEC are extracted and compared on the basis of the same EOF base functions. The degree of IRITEC changes with time is weaker than that of GIMTEC. 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 IRITEC are consistent with those of GIMTEC. However, the change with solar activity is unrepresented by IRITEC. The diurnal variation of the first base function and the annual variation of the fourth base function have a relatively large deviation between IRITEC and GIMTEC.

Results of the spatial and temporal variation characteristic analyses show that the deviation of the first and fourth EOF components between IRITEC and GIMTEC are the two main influencing factors.
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).
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.
The authors declare that they have no conflict of interest.
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).
This paper was edited by Ana G. Elias and reviewed by two anonymous referees.
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/angeo2915012011, 2011.
Abdu, M. A.: Electrodynamics of ionospheric weather over low latitudes, Abdu Geosci. Lett., 3, 11, https://doi.org/10.1186/s4056201600436, 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 lowlatitude region using empirical orthogonal functions, Ann. Geophys., 37, 65–76, https://doi.org/10.5194/angeo37652019, 2019.
Atici, R.: Comparison of GPS TEC with modelled values from IRI 2016 and IRIPLAS over Istanbul, Turkey, Astrophys. Space Sci., 363, 231, https://doi.org/10.1007/s1050901834570, 2018.
Bilitza, D.: International Reference Ionosphere 2000, Radio Sci., 36, 261–275, https://doi.org/10.1029/2000RS002432, 2001.
Bilitza, D.: The International Reference IonosphereStatus 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.: IRI2016: 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 realtime 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: GIMTEC data, available at: ftp://cddis.gsfc.nasa.gov/, last access: 10 April 2019.
COSPAR and URSI: Fortran source code of IRI2016, 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 GPSTEC 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 lowlatitude 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 GPSTEC 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., JacksonBooth, N., Jakowski, N., Hoque, M., HernándezPajares, M., Aragón, Á., María, Á., and OrúsPé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/s0019001811756, 2019.
Kenpankho, P., Supnithi, P., and Nagatsuma, T.: Comparison of observed TEC values with IRI2007 TEC and IRI2007 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 IRI2012 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 GPSTEC 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 IRI2012 and IRI2016 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 F2layer peak height based on ionospheric radiooccultation and groundbased 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/angeo3411092016, 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 MultiGNSS, Satellite Altimetry and Formosat3/COSMIC Data, GPS Solut., 22, 104, https://doi.org/10.1007/s1029101807706, 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/angeo2336152005, 2005.