Towards a Real-Time Description of the Ionosphere: A Comparison between International Reference Ionosphere (IRI) and IRI Real-Time Assimilative Mapping (IRTAM) Models

This paper focuses on a detailed comparison, based on the F2-layer peak characteristics foF2 and hmF2, between the International Reference Ionosphere (IRI), which is a climatological empirical model of the terrestrial ionosphere, and the IRI Real-Time Assimilative Mapping (IRTAM) procedure, which is a real-time version of IRI based on data assimilation from a global network of ionosondes. To perform such a comparison, two different kinds of datasets have been considered: (1) foF2 and hmF2 as recorded by 40 ground-based ionosondes spread all over the world from 2000 to 2019; (2) foF2 and hmF2 from space-based COSMIC/FORMOSAT-3 radio occultation measurements recorded from 2006 to 2018. The aim of the paper is to understand whether and how much IRTAM improves IRI foF2 and hmF2 outputs for different locations and under different diurnal, seasonal, solar and magnetic activity conditions. The main outcomes of the study are: (1) when ionosonde observations are considered for validation, IRTAM significantly improves the IRI foF2 modeling both in accuracy and precision, while a slight improvement in the IRI hmF2 modeling is observed for specific locations and conditions; (2) when COSMIC observations are considered for validation, no noticeable improvement is observed from the IRTAM side for both foF2 and hmF2. Indeed, IRTAM can improve the IRI foF2 description only nearby the assimilated ionosonde locations, while the IRI hmF2 description is always more accurate and precise than IRTAM one.

real-time measurements in a background empirical model is one of the most applied and fruitful methodologies for the real-time specification of the ionospheric electron density. In this context, the knowledge of the large-scale climatological behavior of the ionosphere provided by the underlying background empirical model is complemented with the small-scale weather information provided by real-time assimilated data. The effectiveness and quality of such data-assimilation procedures is critically dependent on the applied algorithm, on the quality, spatial distribution, and availability of assimilated data, and of course on the underlying background empirical model.
Empirical models, such as IRI, are based on analytical formulations whose numerical coefficients are obtained on the basis of the underlying datasets; as a consequence, when new datasets are released, it is of utmost importance to validate the model against new data and eventually recalculate the model's coefficients with the inclusion of the newest data. This validation and recalculation scheme is an ongoing process for empirical models and leads to the continuous improvement of the model itself. Over the years, IRI underwent many validation studies and comparisons with other ionospheric models [21][22][23][24][25][26][27][28]; on the contrary, validations of the IRTAM model are restricted to the works by Vesnin [29] and Galkin et al. [18] for specific locations and conditions. Due to the ever-growing importance that IRTAM is gaining as the most used and affirmed real-time specification of the ionosphere, it is important to validate its performances against large and different datasets to quantify the improvement made by IRTAM in the description of ionospheric weather when compared to the climatological representation made by IRI.
In the present paper, a global validation of the ionosphere F2-layer peak characteristics as modeled by IRI and IRTAM is presented. Specifically, the IRI and IRTAM models, the latter assimilating both the F2-layer ordinary critical frequency (foF2) and the F2-layer peak height (hmF2) from ground-based ionosondes, have been validated according to two different datasets: (1)  Corresponding results are represented in the form of grids of values as a function of the local time (LT) and month of the year, for three different levels of solar activity, for the different ionosonde locations. Moreover, the spatial variation of the calculated statistical values has been investigated through the COSMIC dataset. Comprehensive statistical results are provided for the entire ionosonde and COSMIC datasets as distribution of residuals, density plots, and residuals deviation ratio values, allowing us to draw a complete picture of IRI and IRTAM performances in the description of the F2-layer peak characteristics. As far as we know, it is the first time that IRI and IRTAM are crossvalidated on the basis of such a large dataset covering very different conditions and locations. Moreover, the use of foF2 and hmF2 datasets from different measurement techniques, such as ionosonde and radio occultation, represents an added value in the validation process because it allows us to validate IRTAM against independent data (i.e., COSMIC RO data) and evaluate how much IRTAM is tied to the assimilated data from ionosondes.
A brief description of both IRI and IRTAM models will be provided in Section 2. An overview of the two different datasets used for validation and some information about the runs of IRI and IRTAM models are given in Section 3. The statistics metrics, the binning procedures, and the graphical representation of the results are the subject of Section 4. The validation results for foF2 are described in Sections 5 and 6, while those for hmF2 are described in Sections 7 and 8; the validation shown in Sections 5 and 7 is based on ionosonde data, while that shown in Sections 6 and 8 is based on COSMIC RO data. Final analyses and considerations are the subject of Section 9, while the conclusive remarks are outlined in Section 10.

IRI
IRI is a project started in 1968 by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI) with the aim to develop an international standard for the terrestrial ionosphere. In April 2014, IRI became the official International Standardization Organization (ISO) standard for the ionosphere [20]. The IRI model is empirical, based on ground and space data, and provides monthly medians of the electron density, electron temperature, ion temperature, and ion composition in the altitude range of 60 km to 2000 km [19]. Additionally, it provides the vertical total electron content (vTEC) from the lower boundary to a user-specified upper boundary. Other IRI outputs include the vertical ion drift near the magnetic equator, the F1-layer and spread-F occurrence probability, and the representation of auroral boundaries.
The F2-layer peak plasma frequency foF2, and the corresponding maximum electron density NmF2, are two of the most important parameters when modeling the ionosphere and are related by the formula NmF2 = 1.24·10 10 (foF2) 2 , where NmF2 and foF2 are, respectively, expressed in m −3 and MHz. Both parameters are very important for a reliable characterization of the ionosphere for both scientific and applicative purposes [2,30].
IRI proposes two options for modeling them: the model recommended by the Consultative Committee on International Radio (CCIR) of the International Telecommunication Union (ITU) and the model developed by a special URSI working group [31]. Both models use the same mathematical functions but are based on different datasets of ionosonde data and different methods to fill data gaps, especially over the ocean areas. The CCIR model is recommended mainly for land regions, while URSI is recommended when the region under investigation includes large ocean areas. Specifically, both models are based on a procedure of numerical mapping proposed by Jones and Gallet [32][33][34], which is based on a Fourier time series describing the diurnal variation of monthly medians of foF2 observed at each of the ionosonde stations considered to develop the model (about 150 in total). Legendre special functions (see [35] for further details) are then used to represent the variation of the Fourier coefficients with geographic coordinates and modip, the modified dip coordinate introduced by Rawer [36] to better describe the magnetic field dependence of ionospheric parameters. As driver index both models use the 12-month running mean (IG12) of the ionosonde-based Ionospheric Global (IG) index introduced by Liu et al. [37]. To describe the global daily behavior of foF2, both CCIR and URSI models require 988 coefficients because the diurnal trend is described through a 6th order Fourier time series, then 13 time coefficients, and each time coefficient undergoes a harmonic spatial expansion to the 9th order, then 76 spatial coefficients. This reasoning is valid for a specific month. The seasonal variability is added by sorting data as a function of the month of the year; then, 12 datasets of 988 coefficients are obtained. Finally, the solar activity variation is described by calculating these coefficients for two levels of solar activity, IG12 = 0 for low solar activity and IG12 = 100 for high solar activity, and then performing a linear interpolation between them (see [14] for further details).
Since the IRI-2001 version, the foF2 modeling is also accompanied by a storm option developed by Fuller-Rowell et al. [38] and Araujo-Pradere et al. [39,40] to represent the ionosphere under magnetically disturbed conditions. It is based on the 33 h prior history of the ap magnetic index and gives reliable results, especially at mid latitudes.
Concerning the hmF2 modeling, IRI proposes three options. The first one has been developed by Bilitza et al. [41] with M(3000)F2 ionosonde data and is based on the anticorrelation between hmF2 and the propagation factor M(3000)F2 [35], with M(3000)F2 values from CCIR [35] or URSI [31] mapping procedures (see [14] for further details). The second and the third ones have been recently developed, respectively, by Altadill et al. [42] with hmF2 ionosonde data, and by Shubin et al. [43,44] with hmF2 ionosonde and COSMIC RO data, and they are both based on the spherical harmonic formalism.

Real-Time IRI and the IRTAM Method
In recent years, different data-assimilation techniques have been applied with the intention of improving the IRI output. The aim is to move from the climatological representation provided by the standard IRI model to a description of the ionospheric weather conditions based on the ingestion of real-time measurements.
Bilitza et al. [45] used worldwide ionosonde data from 1986 to 1989 to obtain equivalent IG indices. Komjathy et al. [46], Hernandez-Pajares et al. [47], Ssessanga et al. [48], and Habarulema and Ssessanga [49] used Global Navigation Satellite System data to determine R12 and IG12 equivalent indices. Recently, Pignalberi et al. [14,50,51] proposed a new data-assimilation method, based on ionosonde data, to update the IRI model in the European region through the calculation of an effective IG12; their procedure has been recently updated by assimilating also vTEC values [52]. Pezzopane et al. [11,12] and Pietrella et al. [15] assimilating ionosonde F2-layer peak parameter measurements first determined an effective sunspot number which is used by the Simplified Ionospheric Regional Model [16], and then applied an interpolation technique to assimilate into IRI the full electron density profile recorded by ionosondes.
The assimilation can be carried out as a post-processing activity, without requiring a real-time analysis; however, in recent years, big steps forward have been taken towards a real-time IRI, performing a real-time assimilation. With regard to this, very good results have been achieved by Galkin et al. [17,18], who proposed the IRTAM method (http://giro.uml.edu/RTAM, accessed on 3 August 2021). IRTAM assimilates real-time foF2, hmF2, B0, and B1 measurements from the worldwide network of Digisonde stations (the Global Ionospheric Radio Observatory-GIRO) and uses the URSI procedure [31] to represent the difference between data and model and update the 988 coefficients of the corresponding harmonic expansion for the specific time of the assimilation. As a final step, IRTAM generates near real-time maps of foF2, hmF2, B0, and B1 every 15 min.
IRTAM uses the URSI procedure [31] to describe both foF2 and hmF2; however, while for foF2 this is exactly that of IRI, for hmF2 this was specifically implemented for IRTAM according to the work of Brunini et al. [53], and differs from any IRI hmF2 formulation. The Brunini et al. [53] procedure uses the URSI one to directly map the hmF2 values, that are those modeled by the Bilitza et al. [41] formulation using the M(3000)F2 values provided by the URSI mapping procedure. Practically, this process is nothing else than a simple hmF2 re-mapping, with the purpose to standardize the formalism between foF2 and hmF2. By virtue of this, the IRTAM hmF2 modeled values cannot be considered as a direct updating of the IRI hmF2 values. Moreover, in this paper, IRI hmF2 values are those output by the Shubin et al. [44] default option, whose formulation differs from the IRTAM one and is also based on different datasets. On the contrary, the IRTAM foF2 modeled values are a direct update of the IRI foF2 values, because both models rely on the same formalism.
IRTAM is a four-dimensional data assimilation method because it does not merely assimilate the current data from the GIRO Digisonde network at the time of assimilation but takes into account also the prior 24 h history of the ionosphere at the assimilated station. This approach has the advantage of increasing the robustness of IRTAM by smoothing out data jitter, outliers, and low-confidence values through the diurnal Fourier analysis [17,18]. Moreover, this approach is really suited to IRI because the time and spatial variations are strictly connected through the 988 coefficients [14]. The IRTAM dataassimilation algorithm is named NECTAR (Non-linear Error Compensation Technique for Associative Restoration, [18]). As a first step, for each assimilated station, NECTAR considers the 24 h values recorded by the ionosonde prior to the assimilation time and calculates the differences between observed and modeled (by IRI) values; the procedure is the same for both foF2 and hmF2. The 24 h time series of detrended values is used to describe the corresponding diurnal trend trough a 6th order Fourier analysis (the same as that of IRI) plus a linear term, for a total of 14 time coefficients. The second step concerns the spatial interpolation of the detrended diurnal coefficients at locations different from assimilating sites. This task is accomplished by NECTAR through a recurrent Hopfield neural network optimizer [54], which is used as a spatial interpolator to smoothly spread the information from the assimilating sites to the entire global grid. Once the global grid of detrended diurnal coefficients is obtained, the third step consists in retrieving the 76 spatial coefficients to be used as correction terms in the Jones and Gallet spatial harmonic expansion. The output of the NECTAR method is then 14 (diurnal) × 76 (spatial) correction coefficients to be added to the original ones. The original coefficients are intended to represent the climate behavior of the ionosphere; as a consequence, the correction coefficients calculated by IRTAM describe the departures from the climatological behavior of the ionosphere.

Observations from Ground-Based Ionosondes
Observations of the F2-layer peak ionospheric characteristics, foF2 and hmF2, measured by 40 ground-based ionosonde stations, located at different latitudes in both hemispheres, during the last two solar cycles (from 1 January 2000 to 31 December 2019), are considered as reference. The selected ionosonde stations are listed in Table 1 with the corresponding geographic coordinates, quasi-dipole (QD) magnetic latitude [55], modip [36], and time coverage of the dataset. Figure 1 illustrates the spatial distribution of the 40 ionosonde stations.  Table 1. The yellow circles depict the ionosonde stations location with the corresponding identification number (Table 1, first column). Blue dashed curves depict the QD magnetic parallels, while red ones depict the modip parallels. Ionosonde measured foF2 and hmF2 values were downloaded from the Digital Ionogram DataBASE [56]. For each station, ionograms were recorded by DPS Digisondes [57], and autoscaled by the Automatic Real-Time Ionogram Scaler with True height analysis (ARTIST) software [58]. ARTIST flags the reliability of autoscaled parameters through the Confidence Score (C-Score) parameter [59] ranging from 0 to 100. For this study, only the most reliable values were considered, namely those with C-Score ≥ 75. foF2 and hmF2 time series have a fifteen-minute time sampling (at minutes 0, 15, 30, and 45 of each Universal Time (UT) hour) according to the sounding repetition rate of most of the ionosondes. In Figure 2a, the percentage of available foF2 and hmF2 values per year, also considering the applied filtering based on the C-score value, is graphically represented for each ionosonde station. By virtue of the fifteen-minute time sampling, 100% of available values per year correspond to 35,040 (35,136 for a leap year).
The considered dataset includes the last two solar cycles as depicted by the F10.781 solar index, i.e., the 81-day running mean of the F10.7 solar index, in Figure 2b. F10.7 is the solar radio flux at 10.7 cm wavelength (2800 MHz) [60] and represents one of the most used solar activity proxies for ionospheric modeling. In particular, its 81-day running mean was used in order to smooth its short-time variability. F10. 7  The F10.781 solar index thresholds were selected by considering the solar activity level experienced in the last two solar cycles, and on the basis of the available datasets of ionosonde and COSMIC derived observations.

Observations from Space-Based COSMIC/FORMOSAT-3 Satellites
Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC/FORMOSAT-3, hereafter COSMIC) was a six low-Earth-orbit microsatellites constellation launched on 15 April 2006. The mission was a collaborative project between the National Space Organization in Taiwan and the University Corporation for Atmospheric Research in the United States. COSMIC satellites were deployed into a circular orbit, with 72° of inclination, at about 800 km of altitude and a separation angle of 30° in longitude between neighboring satellites [61]. Each satellite carries a GPS RO receiver capable of measuring the phase delay of radio waves from GPS satellites as they are occulted by the Earth's atmosphere, thus providing an accurate determination of the ionospheric vertical electron density profile up to the COSMIC satellite altitude. COSMIC RO data were downloaded from the COSMIC Data Analysis and Archive Center (CDAAC, https://data.cosmic.ucar.edu/gnss-ro/cosmic1/, accessed on 3 August 2021).
Specifically, COSMIC retrieved foF2 and hmF2 ionospheric characteristics from 22 April 2006 to 31 December 2018 were considered in this study. For this time range, a total of 3,626,729 COSMIC electron density profiles were available, and the most reliable ones were selected by applying the filtering procedure described in the "Methods" section of Pignalberi et al. [62]. As a consequence of such a filtering procedure, the COSMIC dataset that was used in this study for validation purposes was reduced to 1,791,602 profiles.

IRI and IRTAM Models Runs
In this study, we focus on the F2-layer peak characteristics modeled by the IRI-2016, which is the current version of the IRI model. Specifically, IRI was run for the same time periods and locations covered by the ionosondes and COSMIC datasets by using the IRI Fortran code available at the IRI website (http://irimodel.org/, accessed on 3 August 2021). In this way, a one-to-one comparison between measured and modeled foF2 and hmF2 values is guaranteed.
IRI foF2 values were modeled through the URSI coefficients [31], while hmF2 values were modeled through the Shubin et al. [44] option. Moreover, the IRI model was run with the storm option [39,40] "ON" to take into account the magnetic activity disturbance effect on modeled foF2 and hmF2 values.
Likewise, the IRTAM procedure was also run for the same time periods and locations covered by the ionosondes and COSMIC datasets. The IRTAM runs were made by using the IRTAM Fortran package available on the Global Assimilative Model of Bottomside Ionosphere Timeline (GAMBIT) Consortium website (http://giro.uml.edu/GAMBIT/, accessed on 3 August 2021), with the application of the GAMBIT coefficients (https://ulcar.uml.edu/GAMBIT/GambitCoefficients/, accessed on 3 August 2021) for modeling both foF2 and hmF2. GAMBIT coefficients are the URSI ones corrected by IRTAM through the NECTAR method (see Section 2.2) on the basis of assimilated data. Since GAMBIT coefficients are available from the beginning of 2000, the validation analysis here presented is restricted to years 2000-2019. It has to be pointed out that the number of stations assimilated by IRTAM has changed over the years due to the number of available GIRO Digisondes able to stream information in real-time to GAMBIT. This number has changed from about 10 at the beginning of 2000 to about 60 at the end of 2019.
As a consequence, we should expect different performance of IRTAM for different years; specifically, an increase of its performance as time goes by.

Statistical Metrics Adopted in the Validation Process
To evaluate the IRI and IRTAM models' performances, different statistical metrics were calculated between measured and modeled foF2 and hmF2 values. Specifically, the mean of the residuals (Res. Mean) between measured and modeled (by IRI and IRTAM) values, the root mean square error (RMSE), the normalized root mean square error (NRMSE), and the Pearson correlation coefficient (R) are the considered statistical metrics: where xx = fo when considering foF2 and xx = hm when considering hmF2. Concerning the validation with the COSMIC dataset also the root mean square percentage error (RMSPE) was considered: Differently from NRMSE, when calculating RMSPE, the normalization of the residuals between measured and modeled values is performed point-by-point before calculating the root mean square value. The use of RMSPE is more suitable than that of NRMSE for datasets with very different absolute values for different diurnal, seasonal, or solar activity conditions. In principle, this is true for both datasets, the COSMIC and the ionosonde ones. Nevertheless, the COSMIC dataset is global, that is it varies with the geographical coordinates, while the location of ionosondes is fixed. This is why RMSPE is suited especially for the COSMIC dataset.
In order to evaluate the deviations between models and observations, the residuals deviation ratio parameter Rcw [18] was also computed: where εc represents the residuals between "climatological" (i.e., IRI modeled) and measured values, and εw are the corresponding residuals between "weather" (i.e., IRTAM modeled) and measured values. The statistical parameter Rcw allows us to directly compare IRI and IRTAM models and easily evaluate the corresponding performances. In fact, by definition, Rcw = 1 when both models are on par in terms of their accuracy, while Rcw > 1 when IRTAM is more accurate than IRI and Rcw < 1 when IRI is more accurate than IRTAM.

Data Binning
Both foF2 and hmF2 exhibit variations at very different temporal and spatial scales, that are also a consequence of the varying solar and magnetic conditions. To highlight the main climatological diurnal, seasonal, spatial, solar and magnetic activity variations embedded in the data, both modeled and measured values were binned according to different schemes that are here reported in detail.
The diurnal variation was studied in a twofold way: 1. Data were binned as a function of the LT with fifteen minute-wide bins (96 bins in total). 2. Data were collected in three separate diurnal sectors as a function of the solar zenith angle (SZA): The seasonal variation has been studied in a threefold way: The solar activity variation was studied by selecting the three levels of solar activity (low, mid, and high) defined in Section 3.1 on the basis of the F10.781 solar index, and by binning data as a function of F10.781 in bins two s.f.u.-wide.
The magnetic activity variation was studied by selecting three levels of magnetic disturbance on the basis of the global ap index [63]: For what concerns the spatial variation, in the case of the ionosonde dataset this is directly related to the station's geographic position. Conversely, for the COSMIC dataset, this was investigated using geographic coordinates, the magnetic QD latitude and modip. Specifically: 1. Data were binned as a function of the geographic coordinates in bins that are 2.5°wide in latitude and 5°-wide in longitude. 2. Data were binned as a function of the modip latitude in bins that are 2.5°-wide. 3. Data were binned according to three ranges of modip values: High modip: modip < −60° AND modip > 60°.

Graphical Representation of the Statistical Results
Several graphical representations, depending also on the dataset used for validation, were adopted to highlight at best the different climatological variations exhibited by the data.
The following graphical representations were used for both the ionosonde and the COSMIC datasets:


Statistical distribution of residuals between measured and modeled (by IRI and IRTAM) foF2 and hmF2 values.  Density plots between measured and modeled foF2 and hmF2 values, with the corresponding best linear fit line.
The following graphical representations were used only for the ionosondes dataset: Finally, statistical probability distributions of the residuals deviation ratio parameter Rcw were calculated for the entire ionosonde and COSMIC foF2 and hmF2 datasets.

Statistics on the Full Dataset
The full foF2 ionosonde dataset includes 10,133,987 observations obtained cumulating measurements recorded from 2000 to 2019 at the 40 ionospheric stations under consideration. This dataset has been used to calculate the two-dimensional density plots of IRI-foF2 vs. ionosonde-foF2 and IRTAM-foF2 vs. ionosonde-foF2, along with the histograms of corresponding residuals, as shown in Figure 3. The comparison between IRI and IRTAM results highlights that overall IRTAM turns out to be more accurate and precise than IRI. In fact, the statistical distribution of IRTAM residuals around the zero value is less widespread than the IRI one, and Res. MeanIRTAM = 0.005 MHz is smaller than Res. MeanIRI = 0.051 MHz. Furthermore, the comparisons RMSEIRTAM = 0.674 MHz vs. RMSEIRI = 1.002 MHz and NRMSEIRTAM = 12.270% vs. NRMSEIRI = 18.258% show also that the IRI foF2 modeling precision is significantly improved by IRTAM. This point is also reinforced by the IRTAM density plot which appears less scattered than the IRI one (RIRTAM = 0.959 vs. RIRI = 0.908).
Similar conclusions can be made also when IRI and IRTAM performances are evaluated on both the cumulative foF2 time series binned according to Section 4.2, and to the three different level of solar activity defined in Section 3.1, as highlighted by Table 2, and the foF2 time series available for each considered ionospheric station, as highlighted by Table 3. In fact, Table 2 shows values of Res. Mean, RMSE and NRMSE that are always lower for IRTAM, which means that IRTAM models foF2 better than IRI during daytime, nighttime, solar terminator hours, in all seasons, and during different solar and magnetic activity conditions. The improvement made by IRTAM is quite consistent and, more importantly, does not show any particular diurnal, seasonal, or solar activity dependence (as also evidenced by [29]), i.e., IRTAM outperforms IRI in all the conditions. The assimilation of real-time data is particularly efficient during magnetically disturbed conditions where the ionosphere exhibits very steep spatial gradients of the electron density distribution and very fast time variations.   Table 3 instead shows the statistical parameters (1-4) calculated for the foF2 time series recorded in each of the 40 ionosonde stations listed in Table 1. This table highlights how the IRI foF2 modeling accuracy and precision are improved by IRTAM in practically every ionosonde station regardless the location. Since the number and distribution of the stations assimilated by IRTAM has changed over the years, and critically depends on the real-time availability of data, it is difficult to distinguish between assimilated and nonassimilated stations. This is why we extended our analysis also to COSMIC data that are completely independent of ionosondes' ones and are not assimilated by IRTAM. Table 3. Same as Table 2, but for each of the foF2 datasets relative to the 40 ionosonde stations listed in Table 1.

Diurnal, Seasonal, and Solar Activity Statistics for Different Zonal Sectors
In this section, we are going to show the IRI and IRTAM performances for 5 of the 40 ionospheric stations listed in Table 1 In order to simplify the IRI-IRTAM comparison to the reader, all figures have been generated with the same scale. Since the main goal of the paper is the comparison between IRI and IRTAM performances, we will focus on the differences between the statistical metrics provided by IRI and IRTAM, without going into detail about the corresponding physical and theoretical explanation of the observed trends.     6 show that at Jicamarca for LSA, IRTAM models foF2 more accurately than IRI, especially during pre-sunrise and post-sunset hours. An improvement in the foF2 modeling accuracy made by IRTAM is still observed for MSA around sunrise hours for equinoctial and summer months and during daytime hours for some equinoctial months. For both LSA and MSA, at sunrise both models show high percentages of error (at least 25%). Nevertheless, during nighttime, a significant reduction in the percentage error made by IRTAM can be appreciated, which means a good improvement in the foF2 modeling precision. For HSA, around sunrise in equinoctial and summer months, IRTAM performs better than IRI in terms of accuracy. The relatively high absolute errors (at least 2.0 MHz) made by IRI in equinoctial and summer months are strongly smoothed by IRTAM. Except for the hours around sunrise, where the percentage error remains pretty high (at least 25%), a very strong reduction in the percentage error made by IRTAM is observed in all seasons, which means a considerable improvement in the foF2 modeling precision. Figures 7-9 show that at Ascension Island overall, for LSA, during daytime hours and for all seasons, IRI models foF2 less accurately than IRTAM. A very similar behavior of the two models is observed for the absolute error. However, during the central hours of the day, IRTAM provides a greater foF2 modeling precision, being its percentage errors lower than those of IRI. Overall, for both MSA and HSA, the comparison between IRI and IRTAM Res. Mean grids suggests that IRI is less accurate than IRTAM. The IRTAM model reduces both the absolute and the percentage error, during the first hours of the night and daytime hours, which highlights an improvement in the IRI foF2 modeling precision. Both IRI and IRTAM provide bad performance during post-sunset hours, for which extremely high percentage errors, up to 30% and more, are observed. Figures 10-12 show that at Rome for LSA, a better performance in terms of foF2 modeling accuracy is clearly observed in favour of IRTAM, while RMSE and NRMSE grids show very similar patterns, thus indicating that the precision with which IRI and IRTAM model foF2 is comparable. Overall, for MSA, IRI and IRTAM performances are comparable in terms of accuracy. For HSA, a greater accuracy is provided by IRTAM, during presunrise hours and daytime hours for summer months. For both MSA and HSA, the IR-TAM model lowers the values of the absolute and percentage errors, which means that an improvement in the foF2 modeling precision is achieved. Figures 13-18 show that at Sondrestrom and Thule for the three considered solar activity levels a remarkable reduction in the Res. Mean, RMSE, and NRMSE values is observed when passing from IRI to IRTAM, which means that a significant improvement in the foF2 modeling is achieved in terms of both accuracy and precision. The improvement brought by IRTAM is remarkable at equinoctial and winter months.
From a visual inspection of foF2 grids obtained for the three solar activity ranges, some general considerations on how the IRI and IRTAM performances depend on solar activity can be made: no clear Res. Mean trend related to the solar activity for both IRI and IRTAM is found; IRI and IRTAM models show an increasing RMSE trend as the solar activity increases; a clear decreasing NRMSE trend is detected as the solar activity increases at Jicamarca, Ascension Island, and Rome for IRTAM, and at Thule and Sondrestrom for IRI. For the other cases it is difficult to establish a clear solar activity dependence. However, patterns characterized by higher NMRSE are almost always observed for lower solar activity.

Validation Results for foF2 Based on Radio Occultation Observations
The full COSMIC dataset, comprising 1,791,602 selected electron density profiles measured from 2006 to 2018, was used to retrieve reliable foF2 values for calculating the two-dimensional density plots of IRI-foF2 vs. COSMIC-foF2 and IRTAM-foF2 vs. COSMIC-foF2, along with the histograms of corresponding residuals, as shown in Figure 19. Differently from Figure 3, Figure 19 shows that IRTAM does not improve the foF2 modeling accuracy of IRI. In fact, the IRTAM statistical distribution of residuals around the zero value, very similar to the IRI one, shows a value of Res. MeanIRTAM = 0.095 MHz that is a little bit greater than Res. MeanIRI = 0.067 MHz.
Moreover, the comparisons RMSEIRTAM = 1.120 MHz vs. RMSEIRI = 1.053 MHz and NRMSEIRTAM = 17.748% vs. NRMSEIRI = 16.668% show a slight worsening of the foF2 modeling precision made by IRTAM. This is also supported by the IRTAM density plot which appears more scattered than the IRI one (RIRTAM = 0.888 vs. RIRI = 0.901). However, given that the statistical quantities calculated with the two models do not differ significantly, we can claim that the IRI and IRTAM performances based on the COSMIC foF2 data can be considered somewhat comparable. The statistical results reported in the following Table 4 show that only in a few cases is Res. MeanIRTAM < Res. MeanIRI, while in all cases RMSEIRI < RMSEIRTAM, NRMSEIRI < NRM-SEIRTAM, and RIRI > RIRTAM. This means that the same concerns held for the entire COSMIC foF2 dataset (Figure 19) hold also when IRI and IRTAM performances are evaluated on the COSMIC foF2 time series binned according to the procedures described in Section 4.2, and considering the three levels of solar activity defined in Section 3.1.      20 highlights significant differences between the two models. Specifically, the actual ionospheric conditions are better represented by IRTAM over those areas where a high concentration of ionosondes is available for data assimilation. Specifically, a clear improvement in the foF2 modeling is observed over Europe, South Africa, North America, and Brazil, where the RMSPE percentage error is remarkably reduced when IRTAM is run. Vice versa, over the oceans, deserts (e.g., Sahara, Gobi, Kalahari) and Antarctica as well, a better performance from the IRI side is observed. These circumstances can explain why the results of Figure 19 show a global balance between the two models. As expected from an assimilative model, IRTAM improves considerably IRI predictions over those regions where it is possible to assimilate a large number of data. However, a little unexpected is the observed situation that IRTAM does not seem to reproduce the IRI background in areas very far away from the assimilation sites, which is a constraint that it is usually requested by an assimilative model (see, for example, [11,14]). In fact, since IR-TAM works on the residuals between measured values and IRI values through NECTAR (see Section 2.2), it should correct the background model close to assimilation sites and relax to the background far away from assimilation sites. This is not the case represented in Figure 20. A possible explanation of this could lay on the global nature of the correction made by IRTAM to URSI coefficients. In fact, even if the diurnal correction is maximum at assimilation sites (Section 2.2), corrected diurnal coefficients are used to calculate the corrected spatial coefficients through the Jones and Gallet procedure; as a consequence, the correction affects the spatial harmonic functions that have a global nature. Because of this, the effect of the correction at one assimilation site can influence also a very distant location. On the one hand, this is a clear advantage because the correction has a global nature and is robust against sporadic data outliers but, on the other hand, the quality of the procedure critically depends on the spatial distribution of assimilation sites. It is clear from Figure 20 that the current Digisondes' global distribution is an important constraint for IRTAM. This accounts for the IRTAM performances that are better at mid latitudes than at low latitudes, as visible in Figures 20-23. Figure 21 shows that IRI and IRTAM RMSPE maps obtained as a function of modip and local time are very similar: (a) both IRI and IRTAM show the best performance during daytime hours at northern high latitudes; (b) a worsening of both IRI and IRTAM performance is detected during the central hours of the day around 30° N; c) the highest absolute and percentage errors occur in the latitudinal band (30° S, 30° N) around sunrise hours. This last circumstance was somehow expected from the IRI side because, being an empirical climatological model, IRI can have some limitations in describing the variability of the ionospheric characteristics when they present wide variations in a relatively small time window, as is the case at solar terminator hours. Nonetheless, the bottom panel of Figure  21 tells us that between 30° S and 30° N IRI performance was noticeably improved by IRTAM around sunrise hours. This fact constitutes a clear evidence of how data assimilation can play a key role in improving foF2 modeling in certain situations. Both models exhibit quite large errors during nighttime at mid-low latitudes, as already evidenced by Jicamarca and Ascension Islands grids (Figures 4-9). Figure 22 shows that from a seasonal point of view the performances of IRI and IR-TAM are similar. Both models show most of the problems at low/equatorial latitudes throughout the year, especially in correspondence to the crests of the equatorial ionospheric anomaly, mostly the southern one. For what concerns the northern crest the output of IRI seems to be better than the IRTAM one. At mid and high latitudes, the seasonal dependence of the error associated with both models is remarkable, with lowest errors during summer months and highest during winter months. Figure 23 shows the IRI and IRTAM RMSPE maps obtained as a function of modip and the solar activity. This figure highlights a clear improvement of the IRI performance made by IRTAM only for MSA conditions. The poor results observed for LSA from the IRI side were expected, because the conditions characterizing the solar minimum in 2008/2009 [64,65], much lower and prolonged than earlier minima, represented a challenge for IRI [66][67][68]. At the same time, it is somewhat unexpected that IRTAM cannot improve the IRI output for LSA, even though the model assimilates foF2 values that well represent the very low solar activity conditions of that period. As seen before, most of the improvement made by IRTAM is narrowed to mid latitudes.

Statistics on the Full Dataset
The full hmF2 ionosonde dataset comprises 10,133,987 measurements obtained cumulating time series recorded from 2000 to 2019 at the 40 considered ionospheric stations. This dataset has been used to calculate the two-dimensional density plots of IRI-hmF2 vs. ionosonde-hmF2 and IRTAM-hmF2 vs. ionosonde-hmF2, along with the histograms of corresponding residuals, as shown in Figure 24.
The comparison between IRI and IRTAM highlights that the IRTAM model slightly improves the hmF2 modeling accuracy. In fact, the IRTAM and IRI statistical distributions of residuals around the zero value are very similar and, on the other hand, Res. MeanIRTAM = −1.619 km is just a little bit smaller than Res. MeanIRI = 1.804 km. Additionally, the comparisons RMSEIRTAM = 27.965 km vs. RMSEIRI = 32.993 km and NRMSEIRTAM = 10.097% vs. RMSEIRI = 11.912% testify a slight improvement of the hmF2 modeling precision made by IRTAM. This point is also supported by the IRTAM density plot and the correlation coefficient (RIRTAM = 0.845 vs. RIRI = 0.766).
Compared to the foF2 results (ionosonde dataset, Figure 3), in this case the improvement made by IRTAM is lower. However, we have to consider that, for what concerns hmF2, IRTAM cannot be considered as a direct updating of IRI because of the Brunini procedure [53] for mapping hmF2 (see Section 2.2) and, more importantly, because we applied the current default hmF2 IRI option, i.e., that of Shubin et al. [44]. Hence, IRI and IRTAM hmF2 models have to be considered as completely unrelated. Overall, the statistical results reported in Table 5 show that IRTAM models hmF2 better than IRI in terms both of accuracy and precision, although the degree of improvement is lower than that obtained for foF2 (Table 2).  Table 6 shows the statistical parameters (1-4) calculated for the hmF2 time series recorded in each of the 40 ionosonde stations listed in Table 1. Differently from the foF2 results shown in Table 3, the number of cases for which Res.MeanIRI < Res.MeanIRTAM and for which Res.MeanIRTAM < Res.MeanIRI is practically the same. This circumstance is in accordance with the very similar IRI (1.804 km) and IRTAM (−1.619 km) Res.Mean absolute values calculated when the whole hmF2 ionosonde dataset is considered (bottom row of Table 5). This fact points out that the two models are practically equivalent in terms of accuracy. Nonetheless, the IRTAM RMSE and NRMSE values are always lower than the IRI ones, with the exception of the stations of Anyang and I-Cheon, Kwajalein, and Rome. This fact highlights that IRTAM models hmF2 with a precision better than that of IRI. Table 6. Same as Table 5, but for each of the hmF2 datasets relative to the 40 ionosonde stations listed in Table 1. The number of counts on which the statistics were calculated is reported in the rightmost column.

Diurnal, Seasonal, and Solar Activity Statistics Variations for Different Zonal Sectors
The same criteria described in Section 5.2 for the selection of the ionospheric stations and the validation of the results based on ground-based foF2 ionosonde observations, were also adopted to validate the results of IRI and IRTAM models for hmF2. Therefore, as done in Section 5.2, the statistical grids for the same five ionospheric stations previously selected are shown in Figures 25-39 for hmF2.   This improvement is observed mostly around sunrise and post-sunset hours. Moreover, IRTAM reduces the error made by IRI in the hmF2 modeling around noon in summer.  show that at Ascension Island, for the three considered solar activity levels, the IRI hmF2 modeling accuracy is improved by IRTAM especially during nighttime hours. Overall, for LSA, the hmF2 modeling precision of IRI and IRTAM can be considered equivalent during pre-sunrise and daytime hours, since the RMSE and NRMSE grids at these hours show similar patterns. IRI instead shows a better output during post-sunset hours. For MSA and HSA, concerning the modeling precision, a better performance of IRI is observed during pre-sunrise hours.  show that overall at Rome, for the three considered solar activity levels, IRTAM improves the IRI modeling accuracy, in particular for sunrise and daytime hours in summer. Concerning the modeling precision, both models are somewhat equivalent independently of solar activity. Nonetheless, a slight improvement made by IRTAM is found during the central hours of the day for summer months.  show that at Sondrestrom, for LSA and MSA, the IRI hmF2 modeling accuracy during post-sunset hours is remarkably improved by IRTAM, in particular for winter and equinoctial months. Nevertheless, a worsening of IRTAM performance is observed during pre-sunrise hours in winter. For HSA, a drop of the hmF2 modeling accuracy is observed from the IRTAM side during daytime hours in summer. For LSA and MSA, RMSE and NRMSE grids show very similar patterns, thus indicating that the precision with which IRI and IRTAM model hmF2 is on the whole comparable. Nevertheless, it must be pointed out that a clear improvement when passing from IRI to IRTAM is observed during post-sunset hours in winter.  show that at Thule, for the three considered solar activity levels, a considerable improvement in the hmF2 modeling accuracy is achieved by IRTAM during the whole day especially in winter and equinoctial months. Overall, RMSE and NRMSE grids show that IRTAM improves the hmF2 modeling precision, independently of solar activity. In particular, the improvement is appreciated during the whole day in winter and during daytime hours in summer.
From a visual examination of hmF2 grids achieved for the three solar activity ranges, some generic conclusions on how the IRI and IRTAM performances depend on solar activity can be drawn: no clear Res. Mean trend related to the solar activity for both IRI and IRTAM is found; neither IRI nor IRTAM show a clear RMSE trend related to the solar activity; for both models, a decreasing NRMSE trend is observed as the solar activity increases.

Validation Results for hmF2 Based on Radio Occultation Observations
The full COSMIC dataset, comprising 1,791,602 electron density profiles measured from 2006 to 2018, was used to retrieve reliable COSMIC hmF2 values for calculating the two-dimensional density plots of IRI-hmF2 vs. COSMIC-hmF2 and IRTAM-hmF2 vs. COS-MIC-hmF2, along with the histograms of corresponding residuals, as shown in Figure 40. The figure shows that IRTAM slightly improves the hmF2 modeling accuracy because Res. MeanIRTAM = 0.140 km is smaller than Res. MeanIRI = 3.728 km.
The same cannot be said for the hmF2 modeling precision. In fact, the comparisons RMSEIRTAM = 32.223 km vs. RMSEIRI = 22.446 km and NRMSEIRTAM = 11.609% vs. NRMSEIRI = 8.087% show that the IRTAM absolute and percentage errors are by far higher than those of IRI. Moreover, the IRTAM density plot appears much more scattered than that of IRI (RIRTAM = 0.733 vs. RIRI = 0.881). The statistical results reported in Table 7 show that in some cases Res. MeanIRTAM < Res. MeanIRI, while in many cases RMSEIRTAM > RMSEIRI, NRMSEIRTAM > NRMSEIRI, and RIR-TAM < RIRI, which means that the same considerations made for the full COSMIC hmF2 dataset hold also when IRI and IRTAM performances are evaluated according to the procedures described in Section 4.2, considering the three levels of solar activity defined in Section 3.1.      Figure 40 shows that, differently from foF2, the values of RMSEIRI and NRMSEIRI calculated for hmF2 are relatively smaller than those of RMSEIRTAM and NRMSEIRTAM. This is reflected also in the distribution of residuals that in the case of IRTAM appears more widespread around the zero than that of IRI, and this has a strong impact on the maps shown in Figures 41-44. In fact, these maps clearly highlight how IRTAM, independently of the analysis (spatial, diurnal, seasonal, and dependent on solar activity), does not improve the hmF2 modeling made by IRI. This result might be affected to a some extent by the fact that IRI hmF2 values have been calculated through the Shubin et al. [44] model, which is partly based on electron density profiles collected by COSMIC between 2006 and 2012.

Final Analyses and Comparisons between IRI and IRTAM
The results shown in the previous sections are here summarized through the residual deviation ratio Rcw defined in Section 4. Rcw is a statistical parameter that in general is very suitable to assess definitively the performance of one model over another. To this end, the IRI and IRTAM performances are evaluated analyzing the distributions (in a logarithmic scale) of the residuals' deviation ratio calculated on both the full foF2 ionosonde/COSMIC dataset ( Figure 45) and the full hmF2 ionosonde/COSMIC dataset ( Figure 46).  The log10(Rcw) distribution shown in Figure 45a is clearly "shifted" towards the positive values. Specifically, the mean value of the distribution equal to +0.212 highlights that overall IRTAM performs better than IRI by a factor of about 1.6 when considering the foF2 ionosonde dataset. Instead, the log10(Rcw) distribution shown in Figure 45b is quite symmetric with respect to the zero value. This fact is strongly supported by the mean value of the distribution, which is equal to −0.001. Therefore, the IRI and IRTAM performances can be considered equivalent when considering the COSMIC foF2 dataset.
The log10(Rcw) distribution shown in Figure 46a is quite symmetric around the zero. Specifically, the mean value of the distribution equal to +0.057 (which corresponds to an improvement factor of about 1.14) highlights that IRTAM and IRI provide quite comparable outputs when considering the hmF2 ionosonde dataset. The log10(Rcw) distribution shown in Figure 46b is clearly "shifted" towards the negative values. Specifically, the mean value of the distribution equal to −0.188 points out that IRI performs better than IRTAM when considering the hmF2 COSMIC dataset by a factor of about 1.5.
The same analysis based on log10(Rcw) distribution was applied by Galkin et al. [18] on a dataset of foF2 values recorded by 59 ionosondes during May-June 2019. They found results very similar to the ones shown in Figure 45a, with IRTAM improving IRI by a factor of about two (about 0.3 in the logarithmic scale of Figure 45). The slight differences are due on the one hand to the fact that to test IRTAM Galkin et al. [18] used only data from assimilated stations, while in this study we used also non-assimilated stations, and on the other hand to the larger extension of our dataset covering different seasons and solar activity levels. Results similar to those of Galkin et al. [18] were obtained also by Vesnin [29] by considering a larger dataset covering one solar cycle, but again using only data from assimilated stations to test the model. Vesnin [29] investigated the IRTAM performance also for hmF2 and found that IRTAM improved IRI by a factor of about 1.8. However, in that analysis the oldest Bilitza et al. [41] hmF2 IRI option was used as comparison. When using the newest Shubin et al. [44] default IRI hmF2 option, we find much lower differences between IRTAM and IRI, thus confirming the very important step forward made by IRI about the hmF2 modeling, as on the other hand recently outlined by different authors [69][70][71]. Figures 45 and 46 confirm the general picture outlined by the analyses described in Sections 5-8, i.e., the comparison with ionosonde data highlights how IRTAM significantly improves the foF2 prediction made by IRI, while for hmF2 the performances are quite similar between the two models. Since IRTAM assimilates both foF2 and hmF2 from the GIRO network, we would have expected a similar improvement also in the hmF2 prediction. Besides the obvious differences due to the application of the newest Shubin et al. [44] IRI hmF2 default option, two important points need to be highlighted. First, the IR-TAM hmF2 description is based on the mapping procedure introduced by Brunini et al. [53], which introduces residuals in the range from −10 to 10 km when compared to the original hmF2 values obtained from the Bilitza et al. [41] formulation. The second point is inherent to the hmF2 derivation from ionograms. In fact, while foF2 is a parameter that is directly obtained from ionograms as the maximum ordinary frequency reflected by the ionosphere, hmF2 has to be derived trough a mathematical inversion procedure that, starting from critical frequencies measured at different virtual heights, allows obtaining the electron density values at real heights [72]. This inversion procedure is sensitive to different possible error sources due to the E-valley presence, the interpretation of the F2-layer cusp made by ARTIST, and in general the quality of the ionogram echo traces. All of these matters may represent possible sources of error, that are estimated in the order of 10 km [73]. From the above considerations, it clearly emerges that to obtain reliable hmF2 values is more difficult than to get foF2 ones, even with data assimilation. This is a point that requires further improvements and refining of both the data assimilation procedure and the measurement technique itself.
The comparison with the F2-layer peak characteristics derived from COSMIC RO showed a general worsening of the IRTAM performances in comparison with the IRI ones. Specifically, the comparison with the COSMIC dataset (Figures 20 and 41) highlighted that IRTAM improves the IRI foF2 and hmF2 predictions mainly in regions characterized by a dense ionosonde network. This suggests the extent to which IRTAM is tied to data assimilated by ionosondes and to the corresponding spatial distribution. Since assimilated data are used by IRTAM to update the coefficients of the spherical harmonic analysis underlying the IRI description, we would have expected that the improvements were not restricted to the ionosondes' locations but would embrace at least the whole latitudinal sector where assimilated ionosondes are located. Moreover, the IRTAM description should fade towards that of the IRI in regions where the effect of the assimilated data can be considered negligible. These two points are very important, impacting on IRTAM global performances, and need to be deeply investigated for future versions of IRTAM. Currently, IRTAM assimilates data from about 60 GIRO Digisondes. With the ever-increasing number of available Digisondes, able to provide real-time data, and as a consequence of their more homogeneous spatial distribution, a continuous improvement of the IRTAM performance is expected. However, even if in the future the availability of ionosonde data should increase for both the time resolution and the spatial coverage, the fact that IRTAM is so tied to the underlying IRI model (i.e., to the URSI formalism) represents a limit for the improvement and development of IRTAM itself. In fact, IRTAM through NECTAR, to minimize at assimilation sites the mismatch between measured and IRI modeled values, calculates the corrections to be added to the URSI coefficients, but the order of the diurnal and spatial harmonics is left unchanged. This means that steep spatial gradients and fast time variations that are below the limits that can be resolved with the current spatial and temporal resolution of the URSI formalism would not be represented by IRTAM, even with an increased availability of assimilated data. Since steep spatial gradients and fast time variations are customary under specific Space Weather conditions, and the aim of data-assimilation methods is the reliable representation of such conditions, this poses serious limitations for the IRTAM model that its developers should bear in mind in the future.

Conclusions
In the present paper, we compared the IRI and IRTAM models; the latter, being a real-time version of IRI, is based on the assimilation of ionosonde measurements. In order to assess the performance of the two models, two different datasets have been considered: (1) foF2 and hmF2 from ground-based ionosonde observations; (2) foF2 and hmF2 from space-based COSMIC RO observations. Through different analyses and comparison methodologies, we highlighted the main performances exhibited by both IRI and IRTAM for different locations and under different diurnal, seasonal, solar and magnetic activity conditions.
The main results of the study are:  When ionosonde observations are considered for validation, IRTAM improves significantly the IRI foF2 modeling while it slightly improves the IRI hmF2 modeling.  When COSMIC observations are considered for validation, IRTAM improves neither the IRI foF2 modeling nor the IRI hmF2 modeling.
These results highlight that IRTAM, in contrast to most of assimilation models, has ample room for improvement. The points that in our opinion deserve specific attention are: the bad performance of IRTAM when modeling foF2 at low latitudes; the global hmF2 modeling made by IRTAM which is often unreliable, especially in areas far away from the assimilating sites, where the representation made by IRTAM is at times really different from that of the IRI background; the fact that IRTAM performances are too dependent on the assimilated ionosondes location.
The improvement in the near real-time specification of the ionospheric F2-layer peak characteristics is becoming more and more important nowadays for telecommunication purposes and for Space Weather applications in general. For example, Hartman et al. [74] have recently applied IRTAM as the Floating Potential Measurement Unit (FPMU) backup system that will be used to support the International Space Station (ISS) program. IR-TAM foF2 maps were used by Froń et al. [75] to provide global maps of the ionospheric equivalent slab thickness (τ) parameter that are delivered through the GAMBIT Explorer software (http://giro.uml.edu/GAMBIT, accessed on 3 August 2021). An improved realtime specification of τ on a global basis is very important because this parameter describes the shape of the ionospheric electron density profile; thus, an improved specification of τ can help empirical models such as IRI in the description of the profile shape, especially the topside part [76][77][78][79].
Since in the incoming years the applications based on a near real-time specification of the ionospheric conditions will increase in number, an ever more reliable and robust representation of the ionosphere will become of outstanding importance. This is why near real-time data-assimilation models such as IRTAM need continuous improvement and refining, on the one hand to improve the climatological description of the ionosphere made by IRI, and on the other hand to pave the way for a reliable ionospheric weather description. Acknowledgments: This publication uses data from ionospheric observatories made available via the public access portal of the Digital Ionogram Database (http://ulcar.uml.edu/DIDBase/, accessed on 3 August 2021) of the Global Ionosphere Radio Observatory in Lowell, MA. The authors are indebted to the observatory directors and ionosonde operators for the significant investments of their time, effort, expertise, and funds needed to acquire and provide measurement data to academic research. The authors thank the COSMIC/FORMOSAT-3 team for making freely available Radio Occultation data through the COSMIC Data Analysis and Archive Center (CDAAC, https://data.cosmic.ucar.edu/gnss-ro/cosmic1/, accessed on 3 August 2021). Ivan Galkin and the entire GAMBIT Consortium (http://giro.uml.edu/GAMBIT/, accessed on 3 August 2021) are acknowledged for providing access to IRTAM computations (http://giro.uml.edu/IRTAM, accessed on 3 August 2021/), IRTAM Fortran package, and to the GAMBIT coefficients (https://ulcar.uml.edu/GAM-BIT/GambitCoefficients/, accessed on 3 August 2021). The IRI team is acknowledged for developing and maintaining the IRI model and for giving access to the corresponding Fortran code via the IRI website (http://irimodel.org/, accessed on 3 August 2021). The ap magnetic activity index was downloaded from NASA's Space Physics Data Facility of the Goddard Space Flight Center (https://spdf.gsfc.nasa.gov/pub/data/omni/high_res_omni/, accessed on 3 August 2021). The F10.7 solar index was downloaded through the OMNIWeb Data Explorer website (https://omniweb.gsfc.nasa.gov/form/dx1.html, accessed on 3 August 2021) maintained by the NASA.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are frequently used in this manuscript: