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.

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 asym-metry 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.

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., , 2013Bouya et al., 2012;Uwamahoro and Habarulema, 2015;Talaat and Zhu, 2016;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.

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).

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(Rawer et al., 1978Bilitza, 2015). The most recent version of this model is IRI-2016 (Bilitza et al., 2016(Bilitza et al., , 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 subrou-tines, 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).

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 = X T X. (1) 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 1-year time period (2013) with a 2 h temporal resolution and 37 × 36 spatial grids.
We first organized the data set TEC(Lat, Lon, UT, 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 37×36 = 1332; epoch is arranged according to Universal Time (UT), with an interval of 2 h. The total epoch number of the study period was 12×365 = 4380. After performing EOF decomposition, the base function E k (grid) expressing a spatial pattern and the associated coefficient 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.
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.

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 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: 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 IRI-TEC and GIM-TEC. The 2D linear correlation coefficient ρ for two matrices A and B with M × N dimension is calculated as where A and B are the mean values of matrices A and B, respectively, and they are written as 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 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-modelderived 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.
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 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.

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.
The main EOF base functions extracted from Eq.   The first base function E 1 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 1 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 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 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 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 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.
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 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 4 . 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 4 reached −56 %, and the underestimation is serious. This outcome is consistent with the statistical results in Table 1.

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.
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 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 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 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 longitudi-nal 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 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 IRI-TEC are generally smaller than those of the GIM-TEC, especially for  Table 3 indicate that A 4 exhibits the largest mean relative bias. Figure 9 shows that A 1 -A 6 reflect the time-varying 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 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.
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 1 in Fig. 8a describes the overall average global TEC, and Fig. 10a shows E 11 , the diurnal variation characteristic of E 1 . GIM-TEC and IRI-TEC have similar magnitudes, whereas the diurnal variation of IRI-TEC is insignificant. A 11 of GIM-TEC and IRI-TEC in Fig. 10b shows a pronounced semiannual period. However, A 11 values of GIM-TEC on most days are larger than those of IRI-TEC, and the correlation between the F10.7 index and A 11 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 2 -E 4 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 21 and A 31 of GIM-TEC and IRI-TEC are shown in Fig. 10d and f. These functions evidently demonstrate a semidiurnal variation period. A 21 and A 31 of IRI-TEC during the equinox season are lower than those of GIM-TEC. The correlation between the F10.7 index and A 21 and A 31 of GIM-TEC is also observed. A 41 of GIM-TEC and A 41 of IRI-TEC in Fig. 10h exhibit an evident annual period variation of interhemispheric asymmetry. However, the summerto-winter annual variation of GIM-TEC is much larger than that of IRI-TEC.
The fifth and sixth base functions E 5 and E 6 in Fig. 8e  gitude 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 GIM-TEC and IRI-TEC 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 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 11 , A 21 , and A 31 are highly related to solar activity.
A 11 -A 61 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 11 of IRI-S. Li et al.: Global TEC prediction performance assessment of IRI-2016 model  Figure 10. First base function E i1 and associated coefficient A i1 of the six coefficients A 1 -A 6 according to Eq. (12). The monthly smoothed A i1 of GIM-TEC and daily solar F10.7 index are shown together with A i1 . TEC in Fig. 10b shows quite a large underestimation compared with that of GIM-TEC. The strong correlation between A 11 of GIM-TEC and solar activity is unrepresented by A 11 of IRI-TEC. The diurnal variation of the first base function of GIM-TEC represented by E 11 is partially represented by E 11 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.

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. 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. The authors declare that they have no conflict of interest.
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. This paper was edited by Ana G. Elias and reviewed by two anonymous referees.