Assessment of Global Ionospheric Maps Performance by Means of Ionosonde Data

This work presents a new method for assessing global ionospheric maps (GIM) using ionosonde data. The method is based on the critical frequency at the F2 layer directly measured by ionosondes to validate VTEC (vertical total electron content) values from GIMs. The analysis considered four different approaches to using foF2. The study was performed over one of the most challenging scenarios, the Brazilian region, considering four ionosondes (combined in six pairs) and thirteen GIM products available at CDDIS (Crustal Dynamics Data Information System). Analysis was conducted using daily, weekly, one year (2015), and four years (2014–2017) of data. Additional information from the ionosphere was estimated to complement the daily analysis, such as slab thickness and shape function peak. Results indicated that slab thickness and shape function peak could be used as alternative indicators of periods and regions where this method could be applied. The weekly analysis indicated the squared frequency ratio with local time correction as the best approach of using foF2, between the ones evaluated. The analysis of one-year data (2015) was performed considering thirteen GIMs, where CODG and UQRG were the two GIMs that presented the best performance. The four-year time series (2014–2017) were analyzed considering these two products. Regional and temporal ionospheric influences could be noticed in the results, with expected larger errors during the solar cycle peak in 2014 and at locations with pairs of ionosondes with the larger distance apart. Therefore, we have confirmed the viability of the developed approach as an assessment method to analyze GIMs quality based on ionosonde data.


Introduction
Vertical total electron content (VTEC) based on GNSS (Global Navigation Satellite System) measurements has been one of the main parameters used to describe the ionosphere. Given such valuable information, many efforts have been made in the last few decades to provide reliable VTEC values in terms of global ionospheric maps (GIMs). Since 1998, for instance, several Ionosphere Associated Analysis Centers (IAACs) from IGS (International GNSS Service) have developed methods to provide GIMs [1]. Some of the main analysis centers are established by CODE (Center for Orbit Determination in Europe), ESA (European Space Agency), JPL (Jet Propulsion Laboratory), UPC (Universitat Politècnica de Catalunya), and WHU (Wuhan University). Despite all methods are based on VTEC inputs provided by GNSS measurements, the quality of each product varies, depending mainly on the modeling process, DCB (Differential Code Bias) computation, and shell-model specification [2].
The quality of ionospheric models is commonly defined based on comparisons with reference values from consolidated ionospheric models such as IRI (International Reference Ionosphere [3]) [4][5][6], IGS GIMs [7][8][9], ionosonde data [10][11][12], and VTEC values from altimeter and GNSS data [13][14][15]. In the case of GIMs, mostly of the validation and assessment methods are based on comparison with two complementary sources of reference data: VTEC measurements from altimeters and STEC (slant TEC) variability from independent GNSS receivers. Previous studies related to GIMs assessment concerning this approach can be found in Hernández-Pajares et al. [14] and Roma-Dollase et al. [2], indicating the quality and consistency of IGS GIMs.
A complimentary assessment to benefit the quality indicators of GIM's accuracy would be the incorporation of ionosonde measurements as a reference to specify the "true" ionosphere. The electron density peak in the F 2 layer (N m F 2 ) is commonly measured by ionosondes in terms of the associated critical frequency (foF 2 ), which is specially trustable when they are directly provided from manually calibrated ionosonde datasets [16]. Some studies related to the estimation of foF 2 can be found in works such as Altinay et al. [17] and Pietrella et al. [18].
Several investigations regarding the comparison of VTEC with foF 2 have been performed. Spalla and Ciraolo [19] investigated the calibration of Faraday Rotation measurements by the slab thickness. In this work, the correlation between VTEC and foF 2 is analyzed with a limitation of the assumption of constant slab thickness, used during nighttime over long periods (long enough to allow sufficient points to statistics, but short enough to avoid long periods effects caused by variations of slab thickness). Leitinger et al. [20] present an analysis of the relations between VTEC and foF 2 at regular and extreme behaviors, highlighting the influence of slab thickness for estimating VTEC based on foF 2 or vice versa, considering mid-latitudes. Kouris et al. [21] present a study of the variation of VTEC and foF 2 values, showing a clear correlation between them, specifically when low/medium solar activity periods are considered. A study correlating VTEC and foF 2 , and using foF 2 to predict VTEC, by means of linear regression, was conducted by Cander [22]. Some papers applied the relation between VTEC and foF 2 to perform ionospheric analysis at global [23] and regional approaches [24,25]. In addition, VTEC and foF 2 relations have also been used for improving climatological models, see for instance Habarulema et al. [11] and Pignalberi [26]. However, there is still room for the definition of a systematic GIMs quality assessment based on foF 2 measurements.
Considering the relevance of assessing GIMs and the consistency of foF 2 from ionosonde data, we present a new GIM assessment method. The assessment is based on how much VTEC from GIM can help in the interpolation of directly measured and manually scaled foF 2 considering periods and regions with variable slab thickness (Section 2). We have applied this strategy in one of the most challenging ionospheric scenarios, which comprehends the most recent solar maximum (2014) at low latitudes. Results of the proposed method are presented in Section 3, the discussion in Section 4, and Section 5 shows the conclusions.

Materials and Methods
The assessment of the IGS GIMs has been carried out using ionospheric products available at the Crustal Dynamics Data Information System (CDDIS). The VTEC global ionospheric maps selected for this study are indicated in Table 1. They are provided in IONEX (ionosphere map exchange) format by IGS [27], in two different latencies and correspondingly different availability of receivers: rapid (1-2 days) and final (1-2 weeks). More information about GIM products can be found for instance in Hernandez-Pajares et al. [1] and Roma-Dollase et al. [2].  [28]. In this study, only manually scaled foF2 data were used with a 10-min time resolution. The peculiarity of our validation method is that we form a pair of ionosonde stations. The foF2 variability of the ionosphere between both stations is considered as relatively close to the VTEC variability. To this end, Table 2 presents the selected pairs of ionosonde stations with their respective distances.  The peculiarity of our validation method is that we form a pair of ionosonde stations. The foF 2 variability of the ionosphere between both stations is considered as relatively close to the VTEC variability. To this end, Table 2 presents the selected pairs of ionosonde stations with their respective distances. We have applied the assessment method at different time series: one week; one year and four years. The first analysis aimed to select the assessment strategy with better performance, using one-week data from 2015 with the six possible ionosonde pair combinations ( Table 2) and all the rapid and final GIMs (Table 1). In the second assessment, after the selection of the best strategy, we extended the application of this method for a one-year time series (2015), assessing all GIMs previously used. The third assessment was conducted with two of the products with the best performance in the 2015 time series, CODG and UQRG, considering a four-year time series (2014-2017).
The method of assessment proposed in this paper is based on the correlation between the critical frequency at the F 2 layer and the VTEC [19][20][21]. Figure 2 presents plots comparing the behavior of foF 2 measured by two pairs of ionosondes and the respective VTEC values estimated using two GIM products (CODG and UQRG) for the ionosonde positions. Values presented correspond to the day of year (DOY) 21 of 2015. It can be noticed that foF 2 and VTEC in general present a clear similarity. Some specific behaviors of foF 2 are not always well represented by VTEC but we assume that the quality depends on the GIM product.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 18 and final GIMs (Table 1). In the second assessment, after the selection of the best strategy, we extended the application of this method for a one-year time series (2015), assessing all GIMs previously used. The third assessment was conducted with two of the products with the best performance in the 2015 time series, CODG and UQRG, considering a four-year time series (2014-2017).
The method of assessment proposed in this paper is based on the correlation between the critical frequency at the F2 layer and the VTEC [19][20][21]. Figure 2 presents plots comparing the behavior of foF2 measured by two pairs of ionosondes and the respective VTEC values estimated using two GIM products (CODG and UQRG) for the ionosonde positions. Values presented correspond to the day of year (DOY) 21 of 2015. It can be noticed that foF2 and VTEC in general present a clear similarity. Some specific behaviors of foF2 are not always well represented by VTEC but we assume that the quality depends on the GIM product. In the proposed GIM assessment method, data from one ionosonde ( 1) is estimated based on data from the other ionosonde of the considered pair ( 2). The VTEC variation between both places is used to transport the foF2 from ionosonde 2 to the location of the ionosonde 1. In this way, estimated foF2 can be compared with the real foF2 measured at the first ionosonde ( 1), and, implicitly, the GIM quality can be assessed in front of a manually scaled foF2 measurement. The GIM assessment was based on the root mean square (RMS) of the differences between the values of foF2 measured at the first station ( F ) and the values estimated with GIMs ( F ): In the proposed GIM assessment method, data from one ionosonde (Iono1) is estimated based on data from the other ionosonde of the considered pair (Iono2). The VTEC variation between both places is used to transport the foF 2 from ionosonde 2 to the location of the ionosonde 1. In this way, estimated foF 2 can be compared with the real foF 2 measured at the first ionosonde (Iono1), and, implicitly, the GIM quality can be assessed in front of a manually scaled foF 2 measurement. The GIM assessment was based on the root mean square (RMS) of the differences between the values of foF 2 measured at the first station ( f oF 2Iono1 ) and the values estimated with GIMs ( f oF 2 iono1 ): where n is the number of values estimated. The assessment also considered the normalized root mean square (NRMS): to confirm the results obtained. The rates of improvement (RoI) were also calculated: where f oF 2Iono2 represents the critical frequencies at the F 2 layer measured at the second ionosonde. Two equations were analyzed to carry out foF 2 transportation. The first approach is based on the simple ratio between critical frequency and VTEC values, i.e., on the assumption of a direct linear correlation: where VTEC Iono1 and VTEC Iono2 represent the vertical total electron content values interpolated for the ionosondes positions using a GIM product. The second approach is based on the ratio of the squared critical frequency and VTEC values. Considering the relationship between frequency and VTEC in the ionospheric delay estimation: we are implicitly assuming a linear correlation between VTEC and the peak density, N m F 2 . When taking into account the first-principles relationship between N m F 2 and foF 2 (see for instance Kelley [29]), we formulate the second approach as: In addition to the approaches presented, we also evaluated the use of local time (LT) for the time-collocation vs. the usage of plain Universal Time (UT). In this way, we had analyzed four approaches, considering equation (4) with LT and UT and equation (6) using LT and UT for time-collocation.
In the method proposed, we consider that the behavior of foF 2 is correlated to VTEC, assuming constant slab thickness. One limitation of this assumption is related to regions with significantly distinct ionosphere (and consequently different slab thickness in the two regions considered) [19][20][21]. To analyze this influence, in the first assessment section, we also estimated other ionospheric indicators. First, we analyzed the slab thickness (τ) values [30]: Remote Sens. 2020, 12, 3452 6 of 18 with N m F 2 representing the maximum electron density at the F 2 layer, and, second, we have analyzed the inverse value of the slab thickness, or the shape function peak (S m ): which presents a higher spatial correlation than the electron density profiles [31].

Results
We present the main results divided into three sections. The first one concerns the selection of the best strategy to transport foF 2 based on VTEC values. The second section presents further analysis with daily and weekly data, considering error versus time and distance and additional discussion concerning the estimation of two ionospheric indicators (slab thickness and shape function peak). In the third section, we present two time series analyses, the first with one-year data (2015) and thirteen GIMs and the last one with four-year data (2014-2017) and two ionospheric products.

Strategy Selection
To characterize the state of the ionosphere in the week used for the selection of the best strategy, Figure 3 presents two indices: Dst (Disturbance storm time), on the left, and F 10.7 (solar radio flux at 10.7 cm), on the right. The period corresponds to days 20 to 26 of 2015, the indices were obtained at OMNIWEB [32]. It can be verified that no significant magnetic storms occurred in this week, the minimum Dst obtained (−34 nT) occurred on day 26, it was also the only value to reach the scale of a weak storm (−30 nT to −50 nT) [33]. Regarding the F 10.7 , the values are above 110, indicating the active influence of the Sun.

Strategy Selection
To characterize the state of the ionosphere in the week used for the selection of the best strategy, Figure 3 presents two indices: Dst (Disturbance storm time), on the left, and F10.7 (solar radio flux at 10.7 cm), on the right. The period corresponds to days 20 to 26 of 2015, the indices were obtained at OMNIWEB [32]. It can be verified that no significant magnetic storms occurred in this week, the minimum Dst obtained (−34 nT) occurred on day 26, it was also the only value to reach the scale of a weak storm (−30 nT to −50 nT) [33]. Regarding the F10.7, the values are above 110, indicating the active influence of the Sun.  Each plot corresponds to the results of a distinct strategy, being the first strategy related to the Equation (4), named here as freq 1 , with UT time-colocation. The second strategy is related to freq 1 using LT time-colocation. The third strategy is related to the squared frequency (freq 2 ) of Equation (6) in UT; and the last is related to the squared frequency (freq 2 ) in LT. The left panels show the ionosonde pair with the smallest distance and the right ones are referred to the pairs with the longest distance.
Visually, the influence of the distance between the ionosondes and the error obtained is clear. The results on the right (longest distance) present larger differences and higher variability. It is also possible to notice that the approach considering the ratio with squared frequency of Equation (6) presents better performance for both pairs. Regarding the use of UT or LT, the selection of the best strategy is not so clear from Figure 4 only, except for the larger distance with the frequency ratio, where the use of LT provided better results.   Each plot corresponds to the results of a distinct strategy, being the first strategy related to the Equation (4), named here as freq 1 , with UT time-colocation. The second strategy is related to freq 1 using LT time-colocation. The third strategy is related to the squared frequency (freq 2 ) of Equation (6) in UT; and the last is related to the squared frequency (freq 2 ) in LT. The left panels show the ionosonde pair with the smallest distance and the right ones are referred to the pairs with the longest distance.
Visually, the influence of the distance between the ionosondes and the error obtained is clear. The results on the right (longest distance) present larger differences and higher variability. It is also possible to notice that the approach considering the ratio with squared frequency of Equation (6) presents better performance for both pairs. Regarding the use of UT or LT, the selection of the best strategy is not so clear from Figure 4 only, except for the larger distance with the frequency ratio, where the use of LT provided better results.
Despite of presenting detailed results for two pairs, the other four ionosonde pairs also presented similar behavior. To summarize the values obtained and to allow the selection of the best strategy, we present in Table 3 the mean RMS errors obtained with the same dataset for the six ionosondes pairs. The values from the four approaches were classified from the smaller (lightest gray) to the largest (darkest gray) differences for each pair independently for each GIM. For instance, analyzing the first Remote Sens. 2020, 12, 3452 7 of 18 column of Table 3, i.e., the first pair (FZA0M and SAA0K) of CORG GIM, 0.87 (lightest gray) is the smallest value and 1.17 (darkest gray) is the largest value. In each row, the smaller values between the thirteen GIMs are highlighted in italic, indicating the product with better performance.
Each plot corresponds to the results of a distinct strategy, being the first strategy related to the Equation (4), named here as freq 1 , with UT time-colocation. The second strategy is related to freq 1 using LT time-colocation. The third strategy is related to the squared frequency (freq 2 ) of Equation (6) in UT; and the last is related to the squared frequency (freq 2 ) in LT. The left panels show the ionosonde pair with the smallest distance and the right ones are referred to the pairs with the longest distance.
Visually, the influence of the distance between the ionosondes and the error obtained is clear. The results on the right (longest distance) present larger differences and higher variability. It is also possible to notice that the approach considering the ratio with squared frequency of Equation (6) presents better performance for both pairs. Regarding the use of UT or LT, the selection of the best strategy is not so clear from Figure 4 only, except for the larger distance with the frequency ratio, where the use of LT provided better results.    From this analysis, it can be noticed that the fourth strategy provided the lowest RMS values, which means that the strategy selected to assess GIMs was the ratio considering the squared frequency (Equation (6)) in LT collocation. Indeed, the influence of time transformation on LT is more significant with the simple ratio strategy than when the squared frequency is used in agreement with the typical stationarity of the electron content distribution on a sun-fixed reference frame motivated by the main solar ionization source. It can also be noticed that the larger differences in the time correction to LT are related to larger distances in longitude, once the correction is based on the longitude of the stations. For instance, for the ionosondes CAJ2M and SAA0K, which present similar longitudes (45.0W and 44.3W), the time correction is the same, as are the results in UT and LT, consequently. The results also showed that CODG and UQRG presented the best performance between the products analyzed for the week selected. The influence of the distance between the ionosondes is observed in both foF 2 and VTEC values. The influence of VTEC can be noticed, for instance, around 10 h LT for the second pair (CAJ2M and BVJ03), where the foF 2 measured at BVJ03 is larger than that at CAJ2M. In this case, UQRG better represented the ionospheric variability, leading to a closer estimation of foF 2 to the values measured by the ionosonde. Figure 6 presents the differences between the estimated and measured foF 2 for one day of the same two pairs considered before (with the smallest (a) and largest (b) distances). For these same data, we also present the difference in slab thickness (Figure 7) and the shape function peak (Figure 8). The results obtained in the estimation of foF 2 are correlated to the behavior of the ionosphere reflected by slab thickness and shape function peak. The shape function peak demonstrated to be a good indicator to characterize the ionospheric similarity between two regions. Values presented in Figure 8a are much more similar than those from the most distant pair (Figure 8b). Considering the values for the first pair, the highest variability of the RMS (Figure 6a) occurs close to 0 h, 05 h, and 20 h LT, when similar periods with larger differences were observed in the shape function peak estimation. The higher differences between FAZ0M and SAA0K are associated with the geomagnetic locations and the effects of the equatorial ionization anomaly (EIA). SAA0K is located closer to the magnetic equator, which leads to lower ionization in comparison to FAZ0M (closer to the EIA crest). The peak of difference after 18 LT Remote Sens. 2020, 12, 3452 9 of 18 in Figure 6a is associated with the post-sunset enhancement of ionization, where EIA becomes more evident due to the transition between nighttime and daytime [12]. In the case of CAJ2M and BVJ03 (Figure 6b), they are located at a distinct magnetic hemisphere, with a large distance apart from them, and several reasons justify large foF 2 differences, such as the EIA asymmetries, recombination rates, solar zenith angle, geomagnetic field configurations, and seasonal variabilities. Figure 9a shows the mean RMS differences of estimated foF 2 and measured for one-week data versus the distance of the ionosondes. This indicates that most of the GIM products provide better results predicting the foF 2 spatial variability than directly comparing data from both ionosondes (red line), even when great distances are considered. To verify the influence of ionospheric electron content on error behavior, we present the results from two more weeks in Figure 9: The results indicate that in periods of lower ionospheric electron content, the use of any GIM provides better results than the direct comparison of data from both ionosondes, and the worst performance of some of the products occurred in periods of high ionospheric electron content. However, in general, the performance of the approach with GIMs was better, for most of the products available, than not using any GIM.  Figure 6 presents the differences between the estimated and measured foF2 for one day of the same two pairs considered before (with the smallest (a) and largest (b) distances). For these same data, we also present the difference in slab thickness (Figure 7) and the shape function peak ( Figure  8). The results obtained in the estimation of foF2 are correlated to the behavior of the ionosphere reflected by slab thickness and shape function peak. The shape function peak demonstrated to be a good indicator to characterize the ionospheric similarity between two regions. Values presented in Figure 8a are much more similar than those from the most distant pair (Figure 8b). Considering the values for the first pair, the highest variability of the RMS (Figure 6a) occurs close to 0 h, 05 h, and 20 h LT, when similar periods with larger differences were observed in the shape function peak estimation. The higher differences between FAZ0M and SAA0K are associated with the geomagnetic locations and the effects of the equatorial ionization anomaly (EIA). SAA0K is located closer to the magnetic equator, which leads to lower ionization in comparison to FAZ0M (closer to the EIA crest). The peak of difference after 18 LT in Figure 6a is associated with the post-sunset enhancement of ionization, where EIA becomes more evident due to the transition between nighttime and daytime [12]. In the case of CAJ2M and BVJ03 (Figure 6b), they are located at a distinct magnetic hemisphere, with a large distance apart from them, and several reasons justify large foF2 differences, such as the EIA asymmetries, recombination rates, solar zenith angle, geomagnetic field configurations, and seasonal variabilities.      Figure 9 clearly shows that the increase of the error is not only affected by the distance, but also by other factors, such as the region of the ionosondes considered. Figure 10 presents the same results now divided into two sections: one with the pairs formed with CAJ2M ionosonde (right), located at different latitudes, and one with the other three pairs (left). It can be noticed that the influence of distance has a different impact on the two sections. In the section of pairs at similar latitude (left), the increase of the error related to the distance happens in the period with low electron content (b), for the other section the relation is valid for the period with high electron content (c). Considering the pairs without CAJ2M, the decrease occurs with the pair of ionosondes FZA0M and BVJ03, this behavior could be related to the location of the stations. SAA0K is located closer to the magnetic equator than the other two ionosondes, so in the period with high electron content, the proximity with the equator could lead to a more significant difference in the ionospheric influence, compared to the other station. Considering the pairs with CAJ2M, as previously mentioned, several sources of influence are involved, such as EIA asymmetries, seasonal variabilities, and recombination rates. In addition, this figure also shows the significant improvement in the results with the use of GIMs for the period with lower electron content considering pairs with CAJ2M (b). content. The results indicate that in periods of lower ionospheric electron content, the use of any GIM provides better results than the direct comparison of data from both ionosondes, and the worst performance of some of the products occurred in periods of high ionospheric electron content. However, in general, the performance of the approach with GIMs was better, for most of the products available, than not using any GIM.  Figure 9 clearly shows that the increase of the error is not only affected by the distance, but also by other factors, such as the region of the ionosondes considered. Figure 10 presents the same results now divided into two sections: one with the pairs formed with CAJ2M ionosonde (right), located at different latitudes, and one with the other three pairs (left). It can be noticed that the influence of distance has a different impact on the two sections. In the section of pairs at similar latitude (left), the increase of the error related to the distance happens in the period with low electron content (b), for the other section the relation is valid for the period with high electron content (c). Considering the pairs without CAJ2M, the decrease occurs with the pair of ionosondes FZA0M and BVJ03, this behavior could be related to the location of the stations. SAA0K is located closer to the magnetic equator than the other two ionosondes, so in the period with high electron content, the proximity with the equator could lead to a more significant difference in the ionospheric influence, compared to the other station. Considering the pairs with CAJ2M, as previously mentioned, several sources of influence are involved, such as EIA asymmetries, seasonal variabilities, and recombination rates. In addition, this figure also shows the significant improvement in the results with the use of GIMs for the period with lower electron content considering pairs with CAJ2M (b).       Figure 12 shows the RMS of foF2 differences for 2015 considering the thirteen GIMs and the six ionosonde pairs. Considering the general results, the seasonal influence can be noticed for all the pairs, with higher values from October to March, a period with higher ionospheric electron content in the Southern Hemisphere. In addition, the three pairs with ionosonde CAJ2M presented a higher improvement with the use of GIMs from March to October. It is worth highlighting that CAJ2M is the ionosonde with larger distances in latitude compared to the others, so the behavior of the ionosphere is less similar in this region. The results indicated that for the period with less intense ionospheric electron content, the use of GIMs led to significant improvement in the representation of the different regions considered, confirming the behavior observed in Figure 10b. That indicates this scenario like the one where the use of the GIMs presented the strongest improvements. On the other hand, when a period with higher ionospheric electron content is considered, the error estimated with the help of the GIMs is closer to the simple difference between the values of the ionosondes. These results indicate that besides being a challenging region, when periods with lower ionospheric electron content are considered, GIMs can be able to represent the Brazilian ionospheric variability. However, when periods with higher electron content are considered, the performance of GIMs representing the ionosphere in this region may not be so realistic. The worse performance of GIMs in this period can also be related to the increase of the ionospheric gradients that affect foF2 (local measurement) in a different way than VTEC (integrated measurements). The different influence in the foF2 and VTEC can be provided by the slab thickness, that in the case of the presented method is assumed to be constant. In this way, a worse performance of the method was expected, when regions and periods with different ionospheric gradients are considered. More details on the influence of slab  Figure 12 shows the RMS of foF 2 differences for 2015 considering the thirteen GIMs and the six ionosonde pairs. Considering the general results, the seasonal influence can be noticed for all the pairs, with higher values from October to March, a period with higher ionospheric electron content in the Southern Hemisphere. In addition, the three pairs with ionosonde CAJ2M presented a higher improvement with the use of GIMs from March to October. It is worth highlighting that CAJ2M is the ionosonde with larger distances in latitude compared to the others, so the behavior of the ionosphere is less similar in this region. The results indicated that for the period with less intense ionospheric electron content, the use of GIMs led to significant improvement in the representation of the different regions considered, confirming the behavior observed in Figure 10b. That indicates this scenario like the one where the use of the GIMs presented the strongest improvements. On the other hand, when a period with higher ionospheric electron content is considered, the error estimated with the help of the GIMs is closer to the simple difference between the values of the ionosondes. These results indicate that besides being a challenging region, when periods with lower ionospheric electron content are considered, GIMs can be able to represent the Brazilian ionospheric variability. However, when periods with higher electron content are considered, the performance of GIMs representing the ionosphere in this region may not be so realistic. The worse performance of GIMs in this period can also be related to the increase of the ionospheric gradients that affect foF 2 (local measurement) in a different way than VTEC (integrated measurements). The different influence in the foF 2 and VTEC can be provided by the slab thickness, that in the case of the presented method is assumed to be constant. In this way, a worse performance of the method was expected, when regions and periods with different ionospheric gradients are considered. More details on the influence of slab thickness in the relations between foF 2 and VTEC can be found in Spalla and Ciraolo [19], Leitinger et al. [20], and Kouris et al. [21]. Regarding the general performance of GIMs, it is worth mentioning the consistency of most of the products, presenting similar behavior, mainly if only the final products are considered, similar to other GIM assessment techniques [2].  Table 4 summarizes the mean RMS of the values from 2015 for each product. Table 5 presents the mean improvement percentages compared with the simple difference of the ionosondes foF2, without the representation of the ionospheric variation from GIM. Products with the best improvement rates are highlighted in gray. The lowest rates of improvements were observed in the pair with the smaller distance between the ionosondes, and the highest in the pair with a higher distance. This indicates the relevance of representing the ionospheric variability, especially when larger distances are considered. Analysis using NRMS values presented similar behavior to that obtained with RMS values, confirming the RMS analyses.   Figure 12. RMS values of the differences between foF 2 values estimated with GIMs and measured at the ionosondes for the six pairs of ionosondes for one-year data (2015). Table 4 summarizes the mean RMS of the values from 2015 for each product. Table 5 presents the mean improvement percentages compared with the simple difference of the ionosondes foF 2 , without the representation of the ionospheric variation from GIM. Products with the best improvement rates are highlighted in gray. The lowest rates of improvements were observed in the pair with the smaller distance between the ionosondes, and the highest in the pair with a higher distance. This indicates the relevance of representing the ionospheric variability, especially when larger distances are considered. Analysis using NRMS values presented similar behavior to that obtained with RMS values, confirming the RMS analyses.  In Figure 13, we present the mean RMS from 2015 for each product versus distance. The results show an expected increase in the mean error related to the distance of the ionosonde and that even the products with a worse performance presented improvements in the estimation of critical frequency, compared with the approach without GIM information. These results confirm what was previously indicated with one-week data, and show that, on average, the results with GIM data presented a better performance than the simple differences between the measured values.

Time Series Analysis
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 18 frequency, compared with the approach without GIM information. These results confirm what was previously indicated with one-week data, and show that, on average, the results with GIM data presented a better performance than the simple differences between the measured values. The products with the highest rates of improvement were CODG, IGSG, JPLG, UQRG, WHRG, and WHUG (Table 5). From those products, we selected two (CODG and UQRG) to analyze the RMS of the differences in foF2 estimation considering four years ( Figure 14): 2014 (a), 2015 (b), 2016 (c), and 2017 (d). The same behavior noticed with all the GIMs in 2015 can be observed for CODG and UQRG for the other considered years, except the one considering the ionosonde CAJ2M in 2014, due to lack of data. In addition, with the four-year time series, it is also possible to notice the influence of the solar cycle. The results show that higher RMS occurred closer to the last solar cycle peak (between 2013 and 2014). The products with the highest rates of improvement were CODG, IGSG, JPLG, UQRG, WHRG, and WHUG (Table 5). From those products, we selected two (CODG and UQRG) to analyze the RMS of the differences in foF 2 estimation considering four years ( Figure 14): 2014 (a), 2015 (b), 2016 (c), and 2017 (d). The same behavior noticed with all the GIMs in 2015 can be observed for CODG and UQRG for the other considered years, except the one considering the ionosonde CAJ2M in 2014, due to lack of data. In addition, with the four-year time series, it is also possible to notice the influence of the solar cycle. The results show that higher RMS occurred closer to the last solar cycle peak (between 2013 and 2014). Table 6 summarizes the results presenting the mean RMS values, and Table 7 presents the percentages of improvement compared to the simple difference of the values measured at the ionosondes. In general, as noticed for 2015, the lowest rates of improvement were observed for the pair with the closest ionosondes, at similar latitudes and the highest rates for the most distant pair, involving ionosondes at quite different latitudes. As observed in the RMS values from Figure 13, the rates of improvement also do not follow a regular increase related to the distance. Considering the proximity of the solar cycle peak, it was not possible to notice a clear pattern of improvement for all the pairs. For instance, in the closest pair, the rate improvement decreases, and in the most distant pair, the rate tended to increase. In general, the pairs with ionosondes with a larger difference in latitude (CAJ2M_FZA0M, CAJ2M_SAA0K, and CAJ2M_BVJ03) tended to present an increase in the rates of improvement for data closer to solar minimum, while the pairs with ionosondes with closer latitudes (FZA0M_SAA0K, SAA0K_BVJ03, and FZA0M_BVJ03) tended to present a decrease in the rates of improvement. In fact, by Table 6, considering columns 'off' (indicating the simple difference between the foF 2 of both ionosondes), gradients from each pair can be compared. This is consistent with the gradient for ionosondes in similar latitudes (spatio-longitudinal), which decrease from solar maximum (2014) to the solar minimum (2017); with the ionosondes in different latitudes (spatio-latitudinal) for which the gradient slightly varies (increase or decrease). 2017 (d). The same behavior noticed with all the GIMs in 2015 can be observed for CODG and UQRG for the other considered years, except the one considering the ionosonde CAJ2M in 2014, due to lack of data. In addition, with the four-year time series, it is also possible to notice the influence of the solar cycle. The results show that higher RMS occurred closer to the last solar cycle peak (between 2013 and 2014).

Discussion
We have presented a new way of assessing GIMs, by looking at how much they can help in the interpolation of directly measured and manually scaled foF 2 , applied up to four years of data in one of the most challenging ionospheric scenarios that we can have: low latitude and around the most recent solar maximum. The best strategy, linear interpolation of squared foF 2 based on VTEC ratio, is in agreement with the good behavior of the shape function (electron density divided by the corresponding VTEC) for modeling radio occultations [31].
The performance of the different GIMs presented consistency, especially considering the final products. The highest errors and smaller rates of improvement were obtained with rapid products. The agreement of IGS GIMs was also observed in the assessment performed by Roma-Dollase et al. [2]. Concerning the performance over the Brazilian region, the main improvements were observed in the period of lower ionospheric electron content, considering pairs with ionosondes in different latitudes. This highlights the relevance of the use of GIMs information to approximate the variability of the ionosphere. It could also indicate some limitations of the ionospheric modeling in the Brazilian region in periods with larger electron content. That performance can also be related to the limitation of the method when periods and regions with different slab thickness values are considered. This in agreement with the behavior presented in Spalla and Ciraolo [19], Leitinger et al. [20], and Kouris et al. [21].
Many efforts have been made to model and monitor the ionosphere, especially due to its impact on activities as the ones related to GNSS. Many ionospheric assessment methods rely on ionosonde data, which can be tricky when a limited number of ionosondes are available, mainly considering regions with high ionospheric variability. In this way, the method presented may also be an alternative approach to assessing ionospheric models, by using ionosonde data aided by GIM. A limitation observed in the method is related to periods with significant differences in the ionospheric behavior of the considered regions. A possible improvement of the method could be the use of filtering based on an ionospheric indicator, such as the shape function peak estimation, to identify periods and regions when the method can be applied.

Conclusions
This paper presented a new method for assessing global ionospheric maps (GIM) by means of ionosonde data. Four strategies were presented, the best one was the linear interpolation of squared foF 2 based on the VTEC ratio. The assessment was performed over one of the most challenging scenarios, the Brazilian region and near the last solar cycle peak. The analyses are considered daily, weekly, one-year, and four-year time series and thirteen of the main GIMs products available. The analysis was based on the RMS values and additional ionospheric indicators, slab thickness, and shape function peak were used to complement the daily analysis.
Regarding the application of the method in the Brazilian region, error increase related to distance and seasonal influence could be noticed, with larger errors near 2014 (solar cycle peak) and especially considering pairs with ionosondes in different latitudes. Considering the GIMs selected, CODG, IGSG, JPLG, UQRG, WHRG, and WHUG provided the best improvement in 2015, with mean rates of improvement between 5% and 42% in comparison to not using any GIM. With data from 2014-2017, CODG and UQRG provided improvement rates from 5% up to 49%. Funding: This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES)-Finance Code 001 with the support as well of the UPC.