Ionosonde total electron content evaluation using International Global Navigation Satellite System Service data

In this work, a period of 2 years (2016–2017) of ionospheric total electron content (ITEC) from ionosondes operating in Brazil is compared to the International GNSS (Global Navigation Satellite System) Service (IGS) vertical total electron content (vTEC) data. Sounding instruments from the National Institute for Space Research (INPE) provided the ionograms used, which were filtered based on confidence score (CS) and C-Level flag evaluation. Differences between vTEC from IGS maps and ionosonde TEC were accumulated in terms of root mean squared error (RMSE). As expected, we noticed that the ITEC values provided by ionosondes are systematically underestimated, which is attributed to a limitation in the electron density modeling for the ionogram topside that considers a fixed scale height, which makes density values decay too rapidly above ∼ 800 km, while IGS takes in account electron density from GNSS stations up to the satellite network orbits. The topside density profiles covering the plasmasphere were re-modeled using two different approaches: an optimization of the adapted α-Chapman exponential decay that includes a transition function between the F2 layer and plasmasphere and a corrected version of the NeQuick topside formulation. The electron density integration height was extended to 20 000 km to compute TEC. Chapman parameters for the F2 layer were extracted from each ionogram, and the plasmaspheric scale height was set to 10 000 km. A criterion to optimize the proportionality coefficient used to calculate the plasmaspheric basis density was introduced in this work. The NeQuick variable scale height was calculated using empirical parameters determined with data from Swarm satellites. The mean RMSE for the whole period using adapted α-Chapman optimization reached a minimum of 5.32 TECU, that is, 23 % lower than initial ITEC errors, while for the NeQuick topside formulation the error was reduced by 27 %.


Introduction
The understanding of ionospheric behavior provides important information about the space weather. In addition, the electron content affects group and phase delays of radio waves passing through the ionosphere and impacts, among others, the Global Navigation Satellite System (GNSS). Different instruments are capable of evaluating electron density in the ionosphere, and validations among different sources of data can lead to interesting conclusions. While ionosonde instruments provide "ground truth" measures for the bottom side of the ionospheric profile and estimate the topside using an exponential decay function, ground GNSS stations receiving radio signals from orbiting satellites can provide largescale details of the entire ionosphere structure and even plasmasphere Jakowski, 2005;Reinisch and Galkin, 2011;Jin and Jin, 2011). The analysis proposed in this work is based on comparisons between total electron content (TEC) estimated using density profiles derived from ionograms and vertical total electron content (vTEC) from the International GNSS Service (IGS). While the IGS has its own intrinsic quality control through a ranking system (Hernández-Pajares et al., 2009), ionosonde data are evaluated by its auto-scaling system, and quality scores are assigned to each ionogram. The study was performed in the Brazilian region for a 2-year period (2016)(2017), where ionosonde data from National Institute for Space Research (INPE) were available. Indeed, the plasmaspheric electron density has been considered using two different models: an adapted α-Chapman function (Jakowski, 2005) with a simple optimization and a corrected version of the NeQuick topside formulation (Pezzopane and Pignalberi, 2019).

IGS vTEC maps
IGS vTEC maps are considered to be a reliable ionospheric information product, which was achieved from integrating scientific community efforts (Hernández-Pajares et al., 2009). Such maps are generated by a combination of data from different research institutions within a method that is based on ranking different vTEC maps to compose the final product (Hernández-Pajares et al., 2009). The process begins with raw data from the GNSS ground network being acquired and sent to Ionospheric Associate Analysis Centers (IAACs) so the vTEC maps can be generated using the IONosphere map EXchange (IONEX) format ). To achieve a high level of quality, these vTEC maps are evaluated, and its ability to reproduce corresponding slant TEC (STEC) maps is checked. Next, a combination process takes place using a weighted mean of the available vTEC maps. The final step before making the maps available for access on the IGS server is a validation process. It compares the vTEC maps to an independent source: dual-frequency altimeters onboard the TOPEX, JASON, and ENVISAT satellites.

Ionosonde data
An ionosonde measures the returning echoes of pulse signals at a fixed location to estimate ionospheric characteristics, and the ionogram trace can be processed to result in a vertical electron density profile. The bottom side profile starts with measures at ∼ 90 km up to the peak of the F2 layer (foF2), around 350 km. The ionosphere topside profile, instead, is modeled using an exponential decay function. The integration of electron density in height produces an estimate for the TEC value.
Ionograms can be interpreted either manually by experts or automatically using software. The auto-scaling ionosonde data availability, rather than manual scaling, contributes to meeting practical applications (Jiang et al., 2015). Different systems were created and concentrated efforts have been applied to improve auto-scaling (Reinisch and Xueqin, 1983;Pezzopane, 2002, 2007;Reinisch et al., 2005). Also, a standard archiving output (SAO) format was created by the initiative of the Ionospheric Informatics Working Group (IIWG) to store and disseminate auto-scaled data. Initially, SAO format considered only ionograms scaled by an automatic real-time ionogram scaler with true height (ARTIST); however, it evolved to hold scaled data from other sounder systems (Galkin, 2006).

Ionogram quality
Several attempts had been made to verify ionogram autoscaling system quality (Reinisch et al., 2005;Enell et al., 2016;Pezzopane et al., 2017). Early comparisons between manual and automatic scaled ionospheric parameters revealed limitations in ARTIST system performance due to the absence of quality metrics (Gilbert and Smith, 1988). Recent versions of the ARTIST system improved qualityproof methods, enhancing their results (Bamford et al., 2008;Galkin and Reinisch, 2008) and facing problems related to auto-scaling Stankov et al., 2012).   Ionosonde data scaled by ARTIST 5 have a quality metric called the confidence score (CS). Such a metric is based on quality criteria supported by concepts of ionogram interpretation and algorithms that specify the uncertainty and confidence of scaled results (Galkin et al., 2013). The CS metric includes quality checking solutions introduced by confidence calculation schemes developed since the late 1980s: the auto-scaling confidence level (ACL) quality flag, the twodigit confidence level (C-Level), and the QualScan quality control (McNamara, 2006;Galkin et al., 2013). The estimation of the confidence score occurs during ionogram processing. Ionogram interpretation criteria consider not only analysis of extracted trace shapes, but also ionospheric conditions to compute per-point-error reduction. The CS starts with a value of 100. If an interpretation criterion is found, its perpoint-error value is subtracted from the CS. To be considered acceptable for further use, auto-scaling records need to reach a CS above a predefined threshold value, which is generally 40 (Galkin et al., 2013).
The SAO format version 4 does not have the CS on its specifications, but has the C-Level representation. The two digits' range goes from 11 (highest confidence) to 55 (lowest confidence). The CS produced by ARTIST 5 can be converted to C-Level representation using Table 1 (Galkin et al., 2013).

Methodology
Ionosonde data were obtained from the INPE database using files in SAO format (version 4) and scaled by the ARTIST auto-scaling system (version 5). Data from up to five instruments (see Fig. 1) were available at the same time for the period considered (2016)(2017). Although ionosondes can generate ionograms in less than 10 min intervals, a 1 h interval between soundings was considered in this work, except for the comparisons with IGS data. In that case, a 2 h interval is used to match IGS data availability.
C-Level values were extracted from all ionograms, and despite auto-scaled ionosonde data with a CS above 40 able to be considered acceptable (Galkin et al., 2013), we chose to use only those achieving a CS above 60, corresponding to C-Levels 11 and 22. Figure 2 shows the total number of C-Level flag occurrences for the available data. Figure 3 shows the same distribution, however, considering each ionosonde station separately. It can be seen from Figs. 2 and 3 that the majority of C-Level flag occurrences are in classification levels 11 or 22. Also, more than half of C-Level flags achieve a CS above 80 (C-Level 11). Considering the daily variation in ionosphere electron density, it would be interesting to analyze also the ionosondes' data quality variation within daytime hours taking into account all available data. In Fig. 4 the increase in the occurrence of C-Level 11 after sunset can be seen, while C-Level 22 decreases. The aggregation of C-Level flags 11 and 22 ionograms (used in this work) provides ionosonde data with over 1000 samples even for the period with low occurrences -see the red curve in Fig. 4.
TEC values from IGS maps and ionosondes were compared at the same date or time using the closest geographic correspondence as shown in Table 2, considering the IGS data grid (5 • in longitude per 2.5 • in latitude, every 2 h). The analysis is mainly based on the accumulation of TEC differences by applying the root mean squared error (RMSE) as defined in Eq. (1) (Chai and Draxler, 2014): where errors e i are the differences (with i = 1, 2, . . ., n) between TEC values from the IGS and ionosonde, and n is the total number of values considered.

Experiments and results
Ionospheric total electron content (ITEC) daily variability for each ionosonde C-Level flag is shown in variation, they are coherent. On the other hand, it can be seen that the vTEC values are consistently higher than ITEC for the whole period and for every ionosonde. In Fig. 5, a noisy and incoherent TEC variation for flags greater than 22 in both vTEC and ITEC is noticeable. Obviously, since the results for flags greater than 22 have low confidence, they may have errors, but this reason can not be used to explain the noise in vTEC values. The data representative of higher flags is low; i.e., there is a reduced number of points and heterogeneous distributions during daytime and nighttime. Such unbalanced distribution can produce daily mean TEC representing only daytime or nighttime that during consecutive days leads to noisy curves. In Fig. 6 we can observe that some ionosondes presented a lack of data for a few days or even entire months. The seasonal variation in ITEC was similar for all stations. During the autumn and winter seasons in the Southern Hemisphere we can notice a decrease in ITEC values for both years evaluated.
All the panels in Fig. 7 present daily mean values, considering the density profiles of all ionosondes. In panel (a) the RMSE when comparing ITEC and IGS vTEC is shown. The seasonal variability in TEC differences seems highly correlated with the ionization distribution along the analyzed period. Figure 7b shows the peak of plasma frequency (foF2), and we can observe the periods of high foF2 values corresponding to high RMSE. The maximum altitude used for electron density integration in ionosondes (Fig. 7c) does not change significantly, rarely reaching 900 km. Figure 7d shows the plasma frequency at the maximum altitude of the density profile, which indicates the level from where it is necessary to extend the ionosphere structure evaluation to higher altitudes. Considering the fixed scale height used in digisonde topside profile modeling, such a contribution has not been included in the ITEC calculation, since values of electron density decay too rapidly above ∼ 800 km, and the simple extension of maximum integration altitude is insufficient for proper comparisons to IGS vTEC. This is the main reason why the ITEC values from ionosondes are underestimated when compared to IGS values. It is well known that vTEC values from the IGS data represent the integrated electron density along the signal path between the receiver and the satellite altitude (∼ 20 000 km). Thus, this analysis is in agreement with what is shown in Figs. 5 and 6. Different analytical functions have been used to model the topside ionospheric density profile (e.g., exponential, Epstein, Chapman) (Nsumei et al., 2012;Pignalberi et al., 2018a;Reinisch et al., 2007). These functions and their variations may adopt fixed or variable scale height. In this  work, the ionogram topside profiles were re-modeled using two different approaches that consider fixed and variable scale heights: an adapted α-Chapman exponential decay (Jakowski, 2005) that includes an optimized transition function between the F2 layer and plasmasphere and the NeQuick topside formulation with modeled scale height as a function of a corrected version of the empirical parameter H 0 (Pezzopane and Pignalberi, 2019).

Adapted α-Chapman
The adapted α-Chapman introduced by Jakowski (2005) defines the topside profile N T as The ionograms provided F2 scale height H T , the electron peak density NmF2 that can be derived from measured critical frequency foF2 using NmF2 = (1/80.6)·(foF2) 2 , and peak height hmF2. According to Jakowski (2005), the plas-maspheric scale height H p can be defined as 10 000 km and the plasmaspheric basis density N P 0 is assumed to be proportional to NmF2, i.e., N P 0 = K · NmF2. Using the topside reconstruction of the density profile shown above, the maximum integration height used to estimate ionosonde TEC values was defined as 20 000 km, corresponding to an approximation for the satellite orbit.
In this work, different values for the proportionality coefficient K are examined, and the optimal factor that minimizes the global RMSE is used. Figure 8 shows the ionosonde TEC differences to IGS vTEC in terms of the mean RMSE in the whole period, considering all ionosondes. When K is set to zero, Eq. (2) is reduced to regular α-Chapman decay, and the plasmaspheric slowly decaying exponential term is ignored. As K increases, the underestimated ionosonde TEC values move closer to IGS, hence reducing RMSE. However, after an optimal K, in this experiment equal to 1/175, the plasmaspheric contribution is exceeded, increasing again RMSE.

NeQuick topside formulation
The NeQuick topside analytical formulation (Nava et al., 2008;Pezzopane and Pignalberi, 2019) is based on a semi-Epstein layer describing the topside electron density profile N T as In this approach, the modeled scale height H T is dependent on height h, hmF2, and also on the empirical parameter H 0 : H 0 can be calculated using NmF2, foF2, the propagation factor (M(3000)F2), hmF2, and the smoothed sunspot number (R 12 ) as presented by Nava et al. (2008). Pezzopane and Pignalberi (2019) proposed a new formulation for H 0 based on electron density measurements made by the Swarm satellite constellation. The formulation is where two 2-D grids provide the values of H 0,AC and H 0,B as a function of foF2 and hmF2. Specifically, the grids have been calculated as the median values obtained by using the NeQuick topside formulation, IRI UP modeled values (Pignalberi et al., 2018b, c), Swarm A and C (for H 0,AC ), or Swarm B (for H 0,B ) electron density measurements (Pezzopane and Pignalberi, 2019). In our experiments, hmF2 and foF2 were obtained from the post-processed ionograms and H 0,AC and H 0,B grids were generously provided by Michael Pezzopane and Alessio Pignalberi from the Istituto Nazionale di Geofisica e Vulcanologia, Italy. Figure 9 presents a comparison, considering all ionosondes, among the daily mean ITEC (in blue), IGS vTEC (in red), and density profile integration up to 20 000 km with topside reconstruction using an adapted α-Chapman function (in orange) and a NeQuick formulation (in green). The results are for the years 2016 (top panel) and 2017 (bottom panel). All calculated values follow the reference, i.e., the IGS vTEC variations showing a semi-annual dependency with a minimum in June solstice as well as the day-to-day variability. Indeed, the TEC values obtained with the optimization of an adapted α-Chapman are very similar to the ones from the NeQuick topside procedure, and both are much closer to IGS vTEC than ITEC values.

Comparative evaluation
To identify the best methodology to estimate the plasmaspheric TEC, the daily mean RMSE variation shown in Fig. 10 can be assessed. The differences to IGS vTEC using adapted α-Chapman (orange triangle) and NeQuick (green cross) approaches can be compared to the ITEC errors (blue circle) as shown in Fig. 7a. In general, the RMSE values calculated with adapted α-Chapman and NeQuick are similar. However, the lowest RMSE values along the 2 years 2016 and 2017 belong to the NeQuick criterion, as shown in Fig. 10. In addition, the mean RMSE for the whole period, using adapted α-Chapman reconstruction with the proportionality coefficient optimization, has a minimum value of 5.32 TECU, while using NeQuick topside reconstruction the error is 5.05 TECU.

Conclusions
This paper presented a 2-year period validation of ionosonde data, using IGS vTEC as a reference. Ionogram electron density profiles were first selected based on the confidence score and then integrated in height. As expected, ITEC values were systematically underestimated, which is consistent with an ionospheric topside modeling limitation that uses a fixed scale height, which almost neglects plasma above ∼ 800 km, while IGS data consider electron densities from the GNSS stations up to the satellite. This claim was supported by the examination of Figs. 5, 6, and 7. The ionogram topside profiles were re-modeled using two different approaches: optimization of an adapted α-Chapman exponential decay and the NeQuick topside formulation, based on a semi-Epstein layer with modeled scale height as a function of a corrected version of the H 0 empirical parameter. The electron density integration height was extended to an approximation of satellite orbits. Hence, as expected, for both topside reconstructions the plasmaspheric ionization contribution brought ionosonde TEC values closer to IGS observations. In our experiments, the improvement was significant for determining TEC using ionosonde data, as shown in Figs. 9 and 10. Although both procedures for calculating plasmaspheric TEC yield similar results, the NeQuick cri-terion shows lower RMSE values, as we can clearly see in Fig. 10.
Author contributions. TdSK assisted in conceiving the study, led the implementation of data processing and analysis, and actively contributed to the discussion of results and paper writing. AP proposed the study, assisted with data processing and analysis, and actively contributed to the discussion of results and paper writing. JRdS assisted in conceiving the study and actively contributed to the discussion of results and paper writing. GSF assisted with data acquisition and manuscript writing. ERdP and HFdCV assisted in discussing the results and reviewing the manuscript.
Competing interests. The authors declare that they have no conflict of interest.
Special issue statement. This article is part of the special issue "7th Brazilian Meeting on Space Geophysics and Aeronomy". It is a result of the Brazilian Meeting on Space Geophysics and Aeronomy, Santa Maria/RS, Brazil, 5-9 November 2018.
Acknowledgements. The authors would like to acknowledge the Crustal Dynamics Data Information System (CDDIS), NASA Goddard Space Flight Center, Greenbelt, MD, USA, for providing online access to IGS vTEC data. The authors thank the Embrace/INPE program from MCTIC for providing ionosonde data for this work. The authors thank Michael Pezzopane and Alessio Pignalberi from the Istituto Nazionale di Geofisica e Vulcanologia, Italy, for providing a deeper understanding regarding the implementation of the NeQuick topside formulation and access to median values of H 0,AC and H 0,B as a function of foF2 and hmF2. Review statement. This paper was edited by Inez Batista and reviewed by two anonymous referees.