An Update of the NeQuick-Corr Topside Ionosphere Modeling Based on New Datasets

Abstract

A new analytical formula for H₀, one of the three parameters (H₀, g, and r) on which the NeQuick model is based to describe the altitude profile of the electron density above the F2-layer peak height hmF2, has recently been proposed. This new analytical representation of H₀, called H_0,corr, relies on numerical grids based on two different types of datasets. On one side, electron density observations by the Swarm satellites over Europe from December 2013 to September 2018, and on the other side, IRI UP (International Reference Ionosphere UPdate) maps over Europe of the critical frequency of the ordinary mode of propagation associated with the F2 layer, foF2, and hmF2, at 15 min cadence for the same period. The new NeQuick topside representation based on H_0,corr, hereafter referred to as NeQuick-corr, improved the original NeQuick topside representation. This work updates the numerical grids of H_0,corr by extending the underlying Swarm and IRI UP datasets until December 2021, thus allowing coverage of low solar activity levels, as well. Moreover, concerning Swarm, besides the original dataset, the calibrated one is considered, and corresponding grids of H_0,corr calculated. At the same time, the role of g is investigated, by considering values different from the reference one, equal to 0.125, currently adopted. To understand what are the best H_0,corr grids to be considered for the NeQuick-corr topside representation, vertical total electron content data for low, middle, and high latitudes, recorded from five low-Earth-orbit satellite missions (COSMIC/FORMOSAT-3, GRACE, METOP, TerraSAR-X, and Swarm) have been analyzed. The updated H_0,corr grids based on the original Swarm dataset with a value for g = 0.15, and the updated H_0,corr grids based on the calibrated Swarm dataset with a value for g = 0.14, are those for which the best results are obtained. The results show that the performance of the different NeQuick-corr models is reliable also for low latitudes, even though these are outside the spatial domain for which the H_0,corr grids were obtained, and are dependent on solar activity.

1. Introduction

The topside ionosphere mainly refers to the part of the ionosphere extending from the F2-layer peak height (hmF2) to the overlying plasmasphere (e.g., [1]). In order to obtain reliable total electron content (TEC) values for global navigation satellite system (GNSS) signals correction, its modeling has recently become increasingly important (e.g., [2,3,4,5,6,7,8,9,10]). This is why the most widespread empirical ionospheric models extend the topside ionosphere modeling up to the altitude of the GNSS satellites at about 20,000 km. It is worth highlighting that alternative points of view, showing that a significant signal amplification and delay occur in the D and E ionospheric layers, have been recently proposed (e.g., [11]).

The modeling of the topside ionosphere represents a challenge since it cannot be directly monitored by the widely spread ground-based ionosonde network, and consequently information on the topside electron density can be obtained only through costly techniques like topside sounders satellites, incoherent scatter radars (ISRs), in-situ measurements by low-Earth-orbit (LEO) satellites, and radio occultation (RO) of GNSS signals (e.g., [12]). The topside modeling challenge is testified by the difficulty of the most known empirical ionospheric model, the International Reference Ionosphere (IRI) model [13], to provide a reliable representation of this portion of the ionosphere [14,15,16,17,18]. At the moment, IRI offers four options to model the topside ionosphere [13], including the NeQuick one [19,20,21,22,23], which, before the IRI-2020 version, was the default option.

Several analytical functions were studied to represent the topside ionosphere. The most widely considered are the Chapman [24] and exponential functions, or a linear combination thereof, and the Epstein family functions [25]. In order to obtain a reliable description of the topside shape at different altitude ranges, the most important requirement is a monotonic increase trend of the scale height as the altitude increases (e.g., [26]). This is adopted by both the vary-Chapman [4,7,27,28,29,30,31] and the NeQuick approach [21,26,32] and by the Empirical-Canadian High Arctic Ionospheric Model (E-CHAIM; [33]).

It is worth highlighting that when modeling the topside ionosphere, the empirical scale heights often deviate from the theoretical scale heights [26,34]. Essentially, empirical scale heights are directly obtained from electron density observations and proved to be effective in accurately depicting the observed electron density topside profile through assumed analytical parameterizations. In contrast, theoretical scale heights are deduced from our understanding of the ionospheric plasma’s physical state and chemical composition within the framework of the plasma ambipolar diffusion theory. Pignalberi et al. [26] demonstrated that empirical effective scale heights and theoretical ones vary significantly at the F2-layer peak and in the lower topside and converge in the upper topside, approaching similar values within a few hundred kilometers above the F2-layer peak.

In a recent study, Pezzopane and Pignalberi [35] introduced a new analytical formulation for H₀, one of the three parameters (H₀, g, and r) on which the NeQuick model is based to describe the topside vertical electron density profile. H₀ is the value of the scale height at the height hmF2, g is the gradient of the scale height near the F2-layer peak, and r controls the asymptotic behavior of the scale height at infinity. Obtaining this new analytical formulation for H₀, referred to as H_0,corr, was possible by constraining the NeQuick topside formulation to the characteristics of the F2-layer peak (hmF2 and the absolute maximum of electron density of the F2 layer, NmF2), as outlined by the International Reference Ionosphere UPdate (IRI UP) method across the European region [36,37], and to the electron density measurements obtained by Swarm satellites [38,39]. Since the original NeQuick models H₀ as a function of the bottomside thickness parameter [40,41], it exhibits a strong coupling between the bottomside and the topside part of the ionosphere, which in some cases causes inaccurate results [33,42,43]. The new formulation H_0,corr, also based on in-situ measurements, improved the performance of the NeQuick topside representation at middle latitudes [35,44] and low latitudes [45,46].

To obtain the H_0,corr formulation, Pezzopane and Pignalberi [35] used in situ electron density measurements recorded by Langmuir Probes (LPs) onboard Swarm satellites and the characteristics of the F2-layer peak modeled by the IRI UP method from December 2013 to September 2018. The present work, on the one hand, extends these two datasets until December 2021, and on the other hand, it considers the Swarm-calibrated dataset proposed by Lomidze et al. [47]. Since September 2022, Lomidze-calibrated electron density observations have been stored in the L1B EFI-LP data products [48] and considered as the reference dataset. The extension of the Swarm dataset allowed for covering low solar activity conditions that were previously missing by Pezzopane and Pignalberi [35]. Moreover, the IRI UP method has been refined in terms of the quality check routine used to choose the best variogram model. This improved the functionality of the Universal Kriging method in providing more reliable effective index maps for ingestion in the IRI model to obtain updated maps of the critical frequency of the F2 layer, foF2, and of the F2-layer peak height hmF2 [37]. At the same time, driven by Singh et al. [44] and Singh et al. [46] studies, the role of g is investigated by considering values different from 0.125, which is the value usually considered by the NeQuick topside representation.

To understand which are the best H₀ grids underlying the H_0,corr formulation, among those obtained with different datasets, vertical TEC (vTEC) data for low, middle, and high latitudes, recorded from Precise Orbit Determination (POD) antennas from five LEO satellites (COSMIC/FORMOSAT-3, GRACE, METOP, TerraSAR-X and Swarm A) have been considered. Additionally, a validation based on vTEC data, obtained by integrating the RO electron density profiles measured by COSMIC/FORMOSAT-3 from hmF2 to 600 km above hmF2, is shown. The global extension of the validation dataset allowed assessing the performance of NeQuick-corr also for high and low latitudes, which was not considered by Pezzopane and Pignalberi [35].

Section 2 will describe the datasets used in the study, on one side to update the formulation for H_0,corr, and on the other side to test what is the best formulation for H_0,corr in terms of the underlying numerical grids. Section 3 presents a recall of both the NeQuick topside modeling and the formulation for H_0,corr proposed by Pezzopane and Pignalberi [35]. The updated grids obtained by extending the underlying Swarm and IRI UP datasets until December 2021, and those obtained considering the calibrated dataset proposed by Lomidze et al. [47], are shown in Section 4, along with the description and discussion of the validation analyses based on LEO satellites. The summary and conclusions will be the subject of Section 5.

2. Data and Methods

2.1. Swarm Observations

Swarm is a constellation of three LEO satellites in circular near-polar orbits launched at the end of November 2013 by the European Space Agency (ESA) and still in operation [38,39]. After the commissioning phase ended in April 2014, the satellites reached their final orbital configuration. Two of them (Swarm A and C) fly side-by-side with an inclination of 87.35°, an initial altitude of about 460 km, and an initial East–West separation of about 1–1.5° in longitude; the third (Swarm B) has a different orbital configuration with an inclination of 87.75° and an initial altitude of about 520 km. Each satellite takes about 130–140 days to cover all local times. Swarm satellites provide electron density (N_e) observations through a pair of LPs included in the electric field instrument (EFI) payload, with a sampling frequency of 2 Hz [48,49].

For this study, we consider N_e measurements at 2 Hz recorded between 5 December 2013 and 31 December 2021 by LPs on board the three Swarm satellites. Swarm Level 1b LP data (EFIx_LP_1B) contain a time series of N_e measurements along the satellites’ orbit. The EFIx_LP_1B dataset regularly undergoes upgrades leading to substantial improvements in the data quality [50], which are classified through different product baselines identified by two digits in the EFIx_LP_1B filename. The N_e dataset used in this work concerns the baseline 05, for which no calibration has been applied. The N_e-calibrated dataset has been obtained through the coefficients found by Lomidze et al. [47], who proposed a correction of N_e values based on measurements from ISRs, COSMIC-1 radio occultation, and ionosondes. Swarm N_e data are provided with two flags, one indicating the source of measurements (Flags_LP), the other characterizing the plasma density measurements (Flags_Ne), as reported in the Swarm L1b Product Definition [50]. In this study, only the most reliable measurements, those with Flags_LP = 1 and Flag_Ne ≤ 20, were considered.

This study also relies on POD antenna data [51]. POD antennas are Global Positioning System (GPS) receivers through which TEC values are estimated [52]. Level 2 TEC data contain time series of vertical and slant (both absolute and relative) TEC for each GPS satellite in view (at most eight due to the instrumentation scheme). TEC time series are at 0.1 Hz until 15 July 2014, and since then at 1 Hz. This study considers vTEC data for the year 2014 as measured from Swarm A.

2.2. COSMIC/FORMOSAT-3 Observations

The Constellation Observing System for Meteorology Ionosphere and Climate (COSMIC/FORMOSAT-3, hereafter COSMIC-1) mission was a satellite fleet, launched in 2006 and dismissed in 2020, aimed at studying the atmosphere and the ionosphere through the RO technique [53]. It was a constellation of six microsatellites in a circular orbit of 72° of inclination, with a separation angle of 30° in longitude between adjacent satellites, at about 800 km of altitude. The satellites were all equipped with GPS RO receivers that can be used for the determination of the ionospheric electron density profile by measuring the phase delay of radio waves transmitted from GPS satellites as they are occulted by the Earth’s atmosphere. The satellites are also equipped with POD antennas useful for obtaining vTEC measurements.

In this work, the most reliable COSMIC-1 profiles recorded from April 2006 to December 2018 [54] were selected according to the procedure proposed by Pignalberi et al. [26]. TEC data processed by the University Corporation for Atmospheric Research (UCAR) were included for performance evaluation. Each individual slant TEC observation was transformed into a vTEC value through the mapping function proposed by Foelsche and Kirchengast [55], employing a single-layer height of 1300 km. The goodness of this choice was shown by Zhong et al. [56] and Prol and Hoque [57]. vTEC data recorded by COSMIC-1 in 2008, 2011, and 2014 were considered.

2.3. GRACE Observations

The Gravity Recovery and Climate Experiment (GRACE) satellite is a pioneering scientific mission designed to improve our understanding of the Earth’s gravity field. Launched in March 2002 as a joint project between the National Aeronautics and Space Administration (NASA) and the German Aerospace Center (DLR), GRACE consists of two identical satellites flying around Earth in polar orbits, approximately 220 km apart, at an altitude of about 500 km. GRACE operated successfully for over 15 years, far exceeding its planned mission duration. In addition to its contributions to gravity science, the GRACE mission had practical applications in ionospheric sensing. It was equipped with GPS receivers designed for POD, allowing the extraction of TEC data. vTEC data processed by UCAR for the years 2008, 2011, and 2014 were considered.

2.4. METOP Observations

The first Meteorological Operational (METOP) satellite was launched in October 2006 in a polar orbit at an altitude of about 800 km. It was initially designed for global weather forecasting and climate monitoring, providing temperature, humidity, wind speed, and ozone parameters. As part of the METOP payload, the onboard GNSS Receiver for Atmospheric Sounding (GRAS) is dedicated to acquiring GNSS limb-sounding measurements for atmospheric profiling and precise orbit determination [58]. Recently, Prol et al. [59] conducted a comprehensive evaluation of the METOP-A satellite for ionospheric sensing, concluding that the ionospheric data are as accurate as those from COSMIC-1 and FengYun-3D satellites. Given the great capabilities of the mission, vTEC data from METOP satellites are also included in the current study for evaluation purposes. vTEC data processed by UCAR for the years 2008, 2011, and 2014 were considered.

2.5. TerraSAR-X Observations

TerraSAR-X is a radar imaging satellite developed and operated by DLR in collaboration with the Airbus Defence and Space division. Launched on 15 June 2007, TerraSAR-X can acquire high-resolution radar images of the Earth’s surface in a polar orbit at around 515 km of altitude. The satellite payload has no instrument dedicated to the ionospheric research, but it is equipped with a dual-frequency GPS receiver. As shown by Zakharenkova and Cherniak [60], POD data can be effectively used to determine TEC values to properly describe the variabilities of the ionosphere and plasmasphere. vTEC data processed by UCAR for the years 2008, 2011, and 2014 were considered.

2.6. The IRI UP Method

IRI UP [36] is an empirical data assimilation method based on the IRI model. Its aim is the update of the F2-layer description made by the IRI model through the assimilation of foF2 and of the F2-layer propagation factor for a distance of 3000 km, M(3000)F2, routinely recorded by a network of European ionosondes. foF2 is directly related to NmF2 by means of the relation foF2 = (NmF2/(1.24 × 10⁴))^1/2, where foF2 and NmF2 are respectively expressed in the units of MHz and el/cm³. These assimilated characteristics are exploited by IRI UP to compute, at each ionosonde location, effective values (IG_12eff and R_12eff) of IG₁₂ and R₁₂, which are, respectively, the 12-month running-mean of the ionospheric index IG [61] and of the sunspot number R. Starting from these discrete datasets, corresponding effective index maps are generated through the Universal Kriging method [62]. These maps are then used as input for the IRI model to obtain updated values of foF2 and hmF2 over the European region. The IRI UP method proved to be very powerful in obtaining more reliable F2-layer characteristics for both quiet and magnetically disturbed conditions [36,37]. The IRI UP method, assimilating values at a fifteen-minute resolution by twelve European ionosondes (see Table 1 of Pignalberi et al. [36]), outputs updated maps of foF2 and hmF2 over the European region (from 15° W to 45° E in longitude and from 30° N to 60° N in latitude, with a 1° × 1° spatial resolution) every 15 minutes. IRI UP foF2 and hmF2 maps have been calculated for each Swarm satellite’s passage over the European region from 5 December 2013 to 31 December 2021. Ionosonde data were downloaded from the Digital Ionogram DataBASE [63] by means of the SAO Explorer software (https://ulcar.uml.edu/SAO-X/, accessed on 1 April 2024) developed by the University of Massachusetts, Lowell, and are the output of two different autoscaled systems: ARTIST (Automatic Real-Time Ionogram Scaler with True height analysis) [64,65] and Autoscala [66,67]. Specifically, the ARTIST system was applied to autoscale foF2 and M(3000)F2 data from the ionograms recorded by digisondes [68], while Autoscala was applied to autoscale foF2 and M(3000)F2 data from the ionograms recorded by AIS-INGV ionosondes [69] and the VISRC2 ionosonde installed in Warsaw [70].

Pignalberi et al. [37] have improved the first version of the IRI UP method. In addition to considering for the stations equipped with the ARTIST system only ionograms with a Confidence Score [71] greater than 75 [72], they introduced a filter, valid for both ARTIST and Autoscala, to identify wrongly autoscaled ionosonde data and exclude them from the assimilation process. Moreover, they replaced the preliminary quality check used to choose the variogram model in the Universal Kriging method with a new quality check routine (NQCR) based on the variogram’s residuals. NQCR quantitatively identifies the best variogram model used to generate more reliable effective index maps for ingestion in the IRI model to obtain updated foF2 and hmF2 maps. Unlike Pezzopane and Pignalberi [35], whose analyses were based on the first version of IRI UP proposed by Pignalberi et al. [36], in the present study, we applied the IRI UP improvements proposed by Pignalberi et al. [37].

3. The NeQuick Ionospheric Topside Representation and the H_0,corr Formulation

NeQuick models the topside ionosphere [16,21] through a semi-Epstein layer describing the electron density as a function of the height h, starting from NmF2 at the hmF2 height

$$N_{\mathrm{e}}\left(h\right)=4Nm\mathrm{F}2\frac{\exp \left(\frac{h-hm\mathrm{F}2}{H}\right)}{{\left[1+\exp \left(\frac{h-hm\mathrm{F}2}{H}\right)\right]}^2},$$

(1)

where H is a height-dependent empirical scale height

$$H\left(h\right)=H_0\left[1+\frac{r\cdot g\cdot \left(h-hm\mathrm{F}2\right)}{r\cdot H_0+g\cdot \left(h-hm\mathrm{F}2\right)}\right],$$

(2)

with r and g that are constant factors: g = 0.125 is the gradient of H near the F2-layer peak; r = 100 controls the asymptotic behavior of H at infinity. H₀ is instead the value of H at hmF2 and is equal to

$$H_0=k\cdot B_{2\mathrm{bot}},$$

(3)

where B_2bot is the NeQuick bottomside thickness parameter [40,41], which represents the scale height of the F2-layer bottomside region; k is defined as

$$k=3.22-0.0538fo\mathrm{F}2-0.00664hm\mathrm{F}2+0.113\frac{hm\mathrm{F}2}{B_{2\mathrm{bot}}}+0.00257R_{12};$$

(4)

B_2bot is connected to the electron density derivative at the inflection point (dN_e(h)/dh)_max of the semi-Epstein layer, describing the bottomside F2 layer as follows:

$$B_{2\mathrm{bot}}=\frac{0.04774fo\mathrm{F}{2}^2}{{\left(\frac{\mathrm{d}{N}_{\mathrm{e}}\left(h\right)}{\mathrm{d}h}\right)}_{\max }}.$$

(5)

Mosert de Gonzales and Radicella [73], using ionosonde data, proposed the following empirical relationship for (dN_e(h)/dh)_max:

$${\left(\frac{\mathrm{d}{N}_{\mathrm{e}}\left(h\right)}{\mathrm{d}h}\right)}_{\max }=0.01\exp \left[-3.467+1.714\ln \left(fo\mathrm{F}2\right)+2.02\ln \left(M\left(3000\right)\mathrm{F}2\right)\right].$$

(6)

According to Pezzopane and Pignalberi [35], the main idea behind the development of NeQuick-corr is to obtain corrected values of H₀ directly from measurements collected in the topside ionosphere. To accomplish this task, the NeQuick topside Equation (1) is forced to join the F2-layer peak values (hmF2, NmF2), provided by the IRI UP method, and the topside anchor point (h_sat, N_e,sat), provided by Swarm satellites. Since the IRI UP method provides F2-layer peak data over the European region at 1° × 1° spatial resolution, Swarm measurements were gridded over the same region, with the same spatial resolution, for each track of satellites over Europe, and for the whole period under study. In particular, to remove any possible outlier for any Swarm satellite passage over Europe, each N_e measurement falling in a 1° × 1° grid point was subjected to the Median Absolute Deviation (MAD) filtering [47,74]. Then, the mean of N_e measurements in each grid point was calculated. These means of N_e represent the topside anchor points (h_sat, N_e,sat) considered for the calculation of H_0,corr. In the time span considered in this study, the number of tracks cumulated over Europe was 13,293 for Swarm A, 13,092 for Swarm B, and 13,258 for Swarm C. Once gridded at 1° × 1° spatial resolution, the dataset of N_e means consisted of 298,446 values for Swarm A, 295,316 values for Swarm B, and 299,266 values for Swarm C.

A value of H₀ is calculated for every couple of co-located and simultaneously measured (by Swarm) and modeled (by IRI UP) anchor points, which gives rise to a huge dataset of H₀ over the European region. Due to the temporal resolution of IRI UP maps (fifteen minutes according to the ionosondes’ sounding repetition rate), IRI UP maps and Swarm observations are shifted in time by 7 min and 30 s at most. This time difference is negligible when considering the median climatological behavior of these ionospheric parameters. Computed H₀ values are then modeled as a function of foF2 and hmF2, generating for Swarm A and C a joint H₀ dataset, taking into account that these satellites fly at the same altitude very close to each other.

Specifically, H₀ values are arranged as a function of foF2 and hmF2 by applying the two-dimensional binning procedure described in Pignalberi et al. [75], that is:

• Each calculated value of H₀ is associated to a specific pair of values (foF2, hmF2);
• H₀ values are two-dimensionally binned as a function of foF2 and hmF2, with a bin width of 0.25 MHz and 5 km, respectively, within the following ranges: foF2 ∈ [0, 16] MHz; hmF2 ∈ [150, 450] km;
• If the number of H₀ values in the bin is greater than or equal to 10, the corresponding median is calculated; otherwise, the bin is considered statistically insignificant.

At the end of the analysis, two-dimensional binned grids H_0,AC and H_0,B of H₀, based separately on Swarm A and C measurements and on Swarm B measurements, are generated. The selection of the corresponding H₀ value in these two-dimensional grids allows modeling the topside vertical electron density profile through Equations (1) and (2) once a pair of values (foF2, hmF2) is available. It is worth highlighting that this pair of values can be either modeled or measured. Pezzopane and Pignalberi [35] reported a different performance when using the two grids. In particular, they verified that the use of the H_0,AC grid gives better results from hmF2 to approximately the Swarm A and C altitude, while the use of the H_0,B grid gives better results from the Swarm B altitude to higher ones, and they proposed the following corrected version H_0,corr of H₀:

$$\{\begin{array}{l}{H}_{0,\mathrm{corr}}=H_{0,\mathrm{AC}}\begin{array}{ccc}& \mathrm{for}& h=hm\mathrm{F}2,\end{array}\\ H_{0,\mathrm{corr}}=H_{0,\mathrm{AC}}+\left(H_{0,\mathrm{B}}-H_{0,\mathrm{AC}}\right)\frac{h-hm\mathrm{F}2}{600}\begin{array}{ccc}& \mathrm{for}& h\in \left(hm\mathrm{F}2,hm\mathrm{F}2+600\right),\end{array}\\ H_{0,\mathrm{corr}}=H_{0,\mathrm{B}}\begin{array}{ccc}& \mathrm{for}& h\ge hm\mathrm{F}2+600\end{array}\end{array}$$

(7)

All the variables in (7) are in km. As it was highlighted by Pezzopane and Pignalberi [35], for the very few bins for which H_0,B is lower than or equal to H_0,AC, H_0,corr is always assumed constant and equal to H_0,AC for h > hmF2. Moreover, if for a definite pair (foF2, hmF2) only one value of either H_0,AC or H_0,B is available, this is assumed to be equal to H_0,corr for h > hmF2. The original NeQuick H₀ parameter has to be used for those cases for which neither H_0,AC nor H_0,B is available.

Pezzopane and Pignalberi [35] calculated the two-dimensional binned grids H_0,AC and H_0,B considering each Swarm satellite’s passage over the European region from 5 December 2013 to 31 December 2018. In this work, we extended the underlying Swarm and IRI UP datasets until 31 December 2021. Moreover, the Swarm calibrated dataset proposed by Lomidze et al. [47] has also been considered.

4. Results and Discussion

Figure 1 shows a comparison between H_0,AC and H_0,B median values calculated with and without the Lomidze calibration, considering the dataset from 5 December 2013 to 31 December 2018, namely the one used by Pezzopane and Pignalberi [35]. Figure 2 shows the corresponding residuals and percentage residuals. Figure 3 and Figure 4 are the same as Figure 1 and Figure 2 but consider the updated dataset from 5 December 2013 to 31 December 2021.

Figure 1, Figure 2, Figure 3 and Figure 4 show that H_0,AC and H_0,B values obtained through Lomidze calibrated N_e values are sensibly higher than uncalibrated ones by a factor between 10% and 30%. This is because the Lomidze calibration increased Swarm LP N_e values by a factor depending on the considered Swarm satellite and varying between about 9% and 11% in plasma frequency [47]. Consequently, H_{0,AC Lomidze} and H_{0,B Lomidze} would produce higher N_e values in the topside and then higher vTEC values; this agrees with many past studies that pointed to a general underestimation of vTEC modeled by NeQuick (and then by IRI) (e.g., [35,76]). However, it is worth highlighting that an electron density increase in the topside ionosphere can also be obtained by an increase of the g parameter in the NeQuick topside scale height [32]. This has been recently verified by Singh et al. [46], who pointed out that a value of g = 0.15, instead of 0.125 typically adopted in NeQuick, improves the NeQuick-corr topside representation on a global scale.

On the other hand, Figure 1, Figure 2, Figure 3 and Figure 4 also show that the differences between the H_0,AC and H_0,B values obtained with the two N_e datasets (from 5 December 2013 to 31 December 2018 and from 5 December 2013 to 31 December 2021) are negligible for most pairs (foF2, hmF2). The largest differences are scattered and visible above all for very low foF2 values, as is well visible in Figure 5, showing the corresponding residuals and percentage residuals.

Once H_0,AC and H_0,B grids have been calculated, considering as a reference vTEC values recorded by COSMIC-1, GRACE, METOP, TerraSAR-X, and Swarm A satellites, we performed a comparison test of the following five NeQuick topside descriptions:

• The NeQuick original description (the one represented by Equations (1)−(6)), that until the IRI-2016 version was the default topside option of the IRI model;
• The NeQuick-corr topside description based on IRI UP modeled values and Swarm uncalibrated N_e measurements from 5 December 2013 to 31 December 2021, considering the value of g = 0.125 typically adopted in NeQuick;
• The NeQuick-corr topside description based on IRI UP modeled values and Swarm uncalibrated N_e measurements from 5 December 2013 to 31 December 2021, considering g = 0.15 as suggested by Singh et al. [46];
• The NeQuick-corr topside description based on IRI UP modeled values and Swarm calibrated N_e measurements according to Lomidze et al. [47] from 5 December 2013 to 31 December 2021, considering the value of g = 0.125 typically adopted in NeQuick;
• The NeQuick-corr topside description based on IRI UP modeled values and Swarm calibrated N_e measurements according to Lomidze et al. [47] from 5 December 2013 to 31 December 2021, considering g = 0.14.

The last topside description is the result of the same analysis described by Singh et al. [46], but this time performed on the NeQuick-corr topside description based on Swarm calibrated N_e measurements. This analysis highlighted how a value of g equal to 0.14, instead of the value 0.125 usually considered by NeQuick, improves the NeQuick-corr topside representation based on Swarm calibrated N_e measurements. It is worth noting that vTEC values from POD antenna observations are representative of the topside ionosphere and plasmasphere electron content, i.e., from the LEO satellite altitude to GNSS altitude (about 20,000 km from the ground).

The comparison test was made considering vTEC data recorded in 2008 (low solar activity, yearly mean value of R₁₂ = 4.2), in 2011 (medium solar activity, yearly mean value of R₁₂ = 73.5), and in 2014 (high solar activity, yearly mean value of R₁₂ = 107.2). The selection of these three years ensured the best possible coverage in terms of satellite data from different missions.

For each dataset, the Root Mean Square Error (RMSE) between modeled and measured vTEC values

$$\mathrm{RMSE}\left(\mathrm{TECU}\right)=\sqrt{\frac{{\displaystyle \sum _{i=1}^N{\left({\mathrm{vTEC}}_{\mathrm{modeled},i}-{\mathrm{vTEC}}_{\mathrm{measured},i}\right)}^2}}{N}},$$

(8)

where the index i runs on the N available pairs of modeled and measured vTEC, was calculated by selecting (a) the global dataset, (b) data at high latitudes in the Quasi-Dipole (QD) magnetic latitude [77] range (60° N, 90° N] and (60° S, 90° S], (c) data at middle latitudes in the QD magnetic latitude range (30° N, 60° N] and (30° S, 60° S], and (d) data at low latitudes in the QD magnetic latitude range [30° S, 30° N].

Modeled vTEC values were calculated as follows:

$${\mathrm{vTEC}}_{\mathrm{modeled}}={\displaystyle \underset{h_{\mathrm{LEO}}}{\overset{h_{\mathrm{GNSS}}=20,000}{\int }}{N}_{\mathrm{e},\mathrm{modeled}}\left(h\right)\mathrm{d}h},$$

(9)

where h_LEO is the orbit’s altitude of the five LEO satellite missions considered, h_GNSS is the orbit’s altitude of the GNSS satellites, and N_e,modeled is the electron density modeled according to the five aforementioned topside descriptions. For each of these, the F2-layer peak anchor point was modeled through the IRI model default options: URSI coefficients for foF2 [78] and the Shubin model for hmF2 [79]. Since the four different NeQuick-corr implementations require a pair (foF2, hmF2) to select a value of both H_0,AC and H_0,B, only vTEC values associated with (foF2, hmF2) pairs included in the grids shown in Figure 3 have been modeled. Therefore, through the NeQuick-corr formulation, it was impossible to obtain vTEC-modeled values for the whole validation dataset. This issue affects mainly the low latitudes where foF2 and hmF2 values may significantly differ from those recorded over the European region and used to obtain the grids shown in Figure 3.

Table 1. RMSE between modeled and measured vTEC values for the year 2014. The model with the best performance is bolded.

vTEC RMSE (TECU)—2014
Satellite	Dataset	IRI-NeQuick (g = 0.125)	NeQuick-Corr (g = 0.125)	NeQuick-Corr (g = 0.15)	NeQuick-Corr Lomidze (g = 0.125)	NeQuick-Corr Lomidze (g = 0.14)
GRACE	Global	4.803	3.722	2.938	3.119	2.993
	High latitudes	2.867	3.520	2.738	3.087	2.720
	Middle latitudes	2.952	3.432	2.416	2.752	2.437
	Low latitudes	7.133	4.408	3.887	3.731	4.047
Swarm A	Global	6.959	5.741	4.234	5.040	4.259
	High latitudes	4.955	5.409	4.231	4.897	4.216
	Middle latitudes	4.826	5.500	3.894	4.774	3.917
	Low latitudes	9.971	6.610	4.878	5.730	4.953
TerraSAR-X	Global	4.496	3.269	2.097	2.674	2.117
	High latitudes	2.413	2.935	2.040	2.539	2.044
	Middle latitudes	2.558	3.082	1.860	2.447	1.862
	Low latitudes	6.979	3.887	2.529	3.171	2.585
COSMIC-1	Global	4.905	5.014	3.302	4.619	3.552
	High latitudes	4.118	4.403	3.122	4.116	3.315
	Middle latitudes	4.280	4.762	3.052	4.375	3.308
	Low latitudes	5.992	6.017	3.946	5.521	4.236
METOP	Global	2.165	2.315	1.561	1.964	1.477
	High latitudes	1.543	1.710	1.283	1.500	1.203
	Middle latitudes	1.795	2.163	1.429	1.844	1.367
	Low latitudes	2.877	3.352	2.201	2.784	2.064

shows a summary of the statistical analysis in terms of RMSE for the year 2014, where the model with the best performance is bolded.

The NeQuick-corr topside description based on IRI UP modeled values and Swarm uncalibrated N_e measurements, from 5 December 2013 to 31 December 2021, considering g = 0.15, proved to be the best model for the majority of cases. Moreover, even in cases where it is not optimal, the performance of the NeQuick-corr (g = 0.15) topside description is very close to the optimal, as in the case of the METOP dataset, for which the best RMSE is obtained by the NeQuick-corr Lomidze (g = 0.14) option. In general, the performance of NeQuick-corr (g = 0.15) is similar to that of NeQuick-corr Lomidze (g = 0.14), which highlights the key role of the g parameter.

Tables related to the years 2008 and 2011 show similar results and are given in the Supplementary Material. Concerning the RMSE values, it is worth highlighting their overall decrease as the solar activity decreases, according to the overall decrease of both measured and modeled vTEC values. Table 1 and those of Supplementary Material show that the results are pretty independent of the latitude and, taking into account that the NeQuick-corr has been developed over the European region, this shows that its analytical description is robust also outside the spatial domain originally considered.

The results reported in Table 1 are complemented with a more in-depth analysis based on joint probability distributions (JPDs) and distributions of residuals between measured and modeled vTEC values. Figure 6 and Figure 7 show the JPDs related to TerraSAR-X and METOP, while Figure 8 and Figure 9 show the corresponding distributions of residuals. Similar plots for GRACE, Swarm A, and COSMIC-1 can be found in the Supplementary Material. In the Supplementary Material, the reader can also find the JPDs and distributions of residuals related to TerraSAR-X, METOP, GRACE, and COSMIC-1 for the years 2008 and 2011.

Figure 6, Figure 7, Figure 8 and Figure 9 show that without modifying the original value of g = 0.125, both the original NeQuick description (i.e., one of the options of the IRI model) and the NeQuick-corr ones significantly underestimate the measurements. Instead, the use of optimized g values (g = 0.15 when considering uncalibrated Swarm N_e values; g = 0.14 when considering calibrated Swarm N_e values according to Lomidze et al. [47]) is very effective and improves the NeQuick-corr performance, reducing significantly the vTEC underestimation. Figure 6, Figure 7, Figure 8 and Figure 9 highlight that the NeQuick-corr (g = 0.15) option shows the best performance, followed by the NeQuick-corr Lomidze (g = 0.14) option.

Most of the improvements concern the accuracy of the NeQuick-corr model, i.e., the mean of the residuals, while the dispersion of data is less affected by the g optimization, as testified by the standard deviation of the residuals. Consequently, the accuracy improvement corresponds to lower RMSE values, which, on the other hand, are far from being close to zero. This fact suggests that the significance of the parameter r, the one controlling the asymptotic behavior of the scale height at infinity, namely in the plasmaspheric domain, cannot be understated. Indeed, our study concerns the optimization of both H₀ and g parameters based on data recorded in the first hundreds of kilometers above the F2-layer peak, while a non-negligible fraction of vTEC comes from the plasmasphere where an accurate specification of the r parameter becomes relevant [76]. Pignalberi et al. [80] have recently shown that characterizing globally such parameters through RO observations is not an easy task, especially at low latitudes, and understanding what range of values is to be tested represents a challenge nowadays due to a lack of three-dimensional observations in the ionosphere–plasmasphere transition region. Several recent studies have noted differences in the scale height behavior at the plasmaspheric altitudes compared to the topside ionosphere. For instance, Prol et al. [7] performed a statistical analysis of the scale height using topside sounder observations. They found a notable transition from an approximately linear dependence of the scale height on altitude to the parabolic one at around 1000–1500 km heights. Furthermore, Prol et al. [30] found that using Van Allen Probes data, specifically N_e measurements near the plasmapause, improved our understanding of the plasmasphere. These data could be valuable for future studies focused on tracking variations of the r parameter and might likely lead to a reduction of VTEC residuals shown in Figure 8 and Figure 9.

Another feature that stands out from Figure 8 is that at high latitudes, the distribution of residuals of TerraSAR-X shows a double-hump structure when the NeQuick-corr is considered without optimizing the parameter g. The optimization of the parameter g somehow smooths out this feature. The same holds for Swarm A and GRACE, but not for COSMIC-1 and METOP (see also the Supplementary Material). A possible cause of this behavior is that the analyses have been done without performing any local time and seasonal sorting, and this at high latitudes might have triggered such a double hump structure. The investigation of this particular behavior at high latitudes is, however, out of the scope of the present study.

According to Pezzopane and Pignalberi [35], an additional analysis that has been done is based on COSMIC-1 RO electron density profiles from April 2006 to December 2018. For each RO profile, vTEC is calculated by integrating the topside electron density values from hmF2 to 600 km above hmF2, as follows:

$${\mathrm{vTEC}}_{\mathrm{measured}}={\displaystyle \underset{hm\mathrm{F}2}{\overset{hm\mathrm{F}2+600}{\int }}{N}_{\mathrm{e},\mathrm{COS}\mathrm{MIC}-1}\left(h\right)\mathrm{d}h}.$$

(10)

Moreover, each RO electron density profile provides measured values of foF2 and hmF2 that are used to model the IRI-NeQuick and the four NeQuick-corr topside profiles considered in this study. vTEC modeled values are then calculated as follows:

$${\mathrm{vTEC}}_{\mathrm{modeled}}={\displaystyle \underset{hm\mathrm{F}2}{\overset{hm\mathrm{F}2+600}{\int }}{N}_{\mathrm{e},\mathrm{modeled}}\left(h\right)\mathrm{d}h}$$

(11)

In this way, vTEC_modeled values and vTEC_measured values are based on the same F2-layer peak anchor point. With this analysis, the performance of the models is evaluated for the lowest topside section, i.e., where the correction is believed to be most effective, as it is based on Swarm data recorded in the altitude range of 460–520 km. Moreover, this analysis highlights a possible solar activity dependence effect on models’ performance, something that the RMSE values shown in Table 1 and in those of the Supplementary Material seem to testify.

Figure 10, Figure 11 and Figure 12 show the yearly mean of RMSE values calculated between modeled and measured vTEC values respectively for: Northern high latitudes, for QD magnetic latitudes above 60° N; Northern middle latitudes, for QD magnetic latitudes between 30° N and 60° N; low latitudes, for QD magnetic latitudes between 30° S and 30° N. Pezzopane and Pignalberi [35] performed such an analysis only for middle latitudes.

Figures related to Southern high and middle latitudes show similar results as the Northern hemisphere and can be found in the Supplementary Material.

Figure 10, Figure 11 and Figure 12 highlight that the performance of the models shows a rather marked dependence on solar activity, with the RMSE values increasing as the solar activity increases. Even more significant, the performance among different models presents a different dependence on solar activity variations. For example, at middle latitudes, the Lomidze calibration produces worse results for low solar activity and better results for high solar activity when compared to uncalibrated N_e values.

At high latitudes, the NeQuick-corr description, in most cases, cannot improve the IRI-NeQuick performance. This is rather expected, being the NeQuick-corr formulation based on data recorded over the European region. On the other hand, this is something that seems to be restricted to the lowest topside section because Table 1 and those of the Supplementary Material showed that the NeQuick-corr model performs reasonably well also at high latitudes. Concerning the low latitudes, this analysis confirms that, although the NeQuick-corr is based on data recorded over the European region, its performance is generally good, even though, in this case, both the Lomidze calibration and the optimization of the g value do not bring any improvement. This outcome might be ascribed to the restriction of this analysis to vTEC values from hmF2 to 600 km.

5. Summary and Conclusions

This study is a continuation of the work by Pezzopane and Pignalberi [35]. The numerical grids of H_0,corr were updated by extending the underlying Swarm and IRI UP datasets until December 2021. The extension of the dataset facilitated coverage of low solar activity conditions that were not considered before. In addition, regarding Swarm measurements, it also considered the calibration proposed by Lomidze et al. [47]. Correspondingly, two different sets of grids of H_0,corr have been calculated and tested. The test has been made globally for low, middle, and high latitudes, considering both vTEC values measured by POD antennas and those obtained by integrating the RO profiles recorded by LEO satellites. The main outcomes of the study can be summarized as follows:

• Even though the NeQuick-corr formulation is based on datasets recorded over the European region, its performance is deemed satisfactory for both low, middle, and high latitudes when considering the profile up to the GNSS satellites’ altitude. On the other hand, when considering the lowest topside section, from hmF2 to 600 km above hmF2, it was demonstrated that NeQuick-corr provides adequate performance at low and middle latitudes;
• The study highlighted that considering values of the parameter g other than 0.125 (usually adopted) is very effective in mitigating the vTEC underestimation made by the NeQuick model and significantly improves the NeQuick-corr performance, primarily in terms of accuracy;
• The best performance is obtained by the NeQuick-corr topside description corresponding to g = 0.15 and by the NeQuick-corr Lomidze topside description corresponding to g = 0.14;
• RMSE deviates significantly from 0. This fact suggests that the significance of parameter r, which controls the asymptotic behavior in the plasmaspheric domain, is crucial, and its modeling becomes a necessity for an accurate vTEC modeling;
• The performance of different NeQuick-corr options depends on solar activity, with the RMSE increasing as the solar activity increases. This feature is smoothed out when considering optimized values of g. This suggests that the g parameter is most likely dependent on the solar activity level. In fact, especially at middle latitudes, the application of the optimized g parameter is very effective for high but not for low solar activity years. On the other hand, Pignalberi et al. [80] have recently highlighted this dependence.

While Lomidze-calibrated N_e values are currently considered as the nominal Swarm LP N_e dataset, recent investigations by Smirnov et al. [81], Xiong et al. [82], Pignalberi et al. [83] and Haralambous et al. [84] pointed out how Swarm LP N_e values exhibit local time and solar activity variations that were not described by the Lomidze et al. [47] calibration procedure. In particular, it has been recently demonstrated that the Swarm LPs overestimate the plasma densities on the nightside at low and middle latitudes. This is due to an assumption made in the ion densities calculation that the ambient plasma consists exclusively of O⁺ ions, which is frequently violated at nighttime during low solar activity [81,82]. Therefore, in the future, the calculation of H₀ values through Swarm LP observations could be further refined by applying the recommendation resulting from the Swarm LP calibration and validation ongoing process.

The results of the study highlight the need for a revision of the original NeQuick topside formulation (used also by the IRI model). The link of H₀ with the bottomside, included in the NeQuick-corr formulation, could be abandoned, but also the constant values of g and r should be critically evaluated. Anyhow, our results suggest the use of the NeQuick-corr formulation, instead of the original one, when modeling the topside from hmF2 to the GNSS altitude.