Positioning performance of the Neustrelitz total electron content model driven by Galileo Az coefficients

To aid single-frequency GNSS users, the Neustrelitz Total Electron Content Model (NTCM) has been proposed as an adequate solution to mitigate propagation errors besides the GPS Klobuchar and Galileo NeQuick models. Using the three effective ionization coefficients broadcast in the Galileo navigation message as driver parameters, the version NTCM-GlAzpar, in general, performs equal to or better than NeQuickG when compared to the total electron content domain. In this work, we performed a global statistical validation of the NTCM-GlAzpar model in the position domain by comparing its results with the results of the Klobuchar and NeQuickG models for the first time. For this purpose, we used the GNSS analysis tool gLAB and its customization capabilities in the Standard Point Positioning mode. The data used for model validation correspond to a one-month period of perturbed solar and geomagnetic activity (December 2014) and another one-month period of quiet conditions (December 2019). The data have a worldwide coverage with up to 73 IGS stations. The statistical analysis of the hourly average 3D position error shows that the root mean squared (RMS) values of the Klobuchar model are 6.71 and 2.75 m for the perturbed and quiet conditions, respectively, whereas the NeQuickG model has RMS values of 4.61 and 2.35 m. In comparison, the corresponding RMS values of 4.36 and 2.32 m of the NTCM-GlAzpar model confirm its better positioning performance for both periods. However, we identify also that the performance of NTCM-GlAzpar slightly worsens toward higher latitudes and at night local time. Simple software adaptations and a low computational cost make NTCM-GlAzpar an alternative practicable algorithm to improve the accuracy of GNSS single-frequency applications.


Introduction
The interaction of Global Navigation Satellite System (GNSS) signals with free electrons of the ionosphere causes signal delays, a major error source for single-frequency positioning applications. Indeed, this dispersive effect may cause link related range errors of up to 100 m for frequencies operating between 1 and 2 GHz. In order to mitigate ionospheric propagation errors, the Ionospheric Correction Algorithm (ICA), also known as the Klobuchar model (Klobuchar 1987), and the NeQuickG model have been successfully implemented by the Global Positioning System (GPS, IS-GPS-200H 2013) and Galileo satellite navigation system (OS SIS ICD 2010), respectively. The Klobuchar model is a vertical delay model based on a thin-shell ionospheric approximation. It defines the mean vertical delay at the GPS L1 frequency as a cosine function with varying amplitude and period, depending on the geomagnetic latitude of the ionospheric piercing point (IPP). The maximum of the cosine is fixed at 14 h local time (LT). For nighttime hours, the vertical ionospheric delay is set at a constant value of 5 ns or 9.24 total electron content units or TECU ( 1 TECU = 10 16 electrons m −2 ). In order to model the amplitude and period of the cosine function, the Klobuchar model uses the eight third-order polynomial coefficients a 0 … a 3 ; 0 … 3 broadcast in the GPS navigation message. Since the Klobuchar model is a two-dimensional approach, an obliquity factor or mapping function is required to transform vertical delay into slant delay at the user level.
Differently, NeQuickG is a comprehensive threedimensional model used by Europe's Galileo system to numerically integrate electron density along the GNSS links. The NeQuickG model used for real-time Galileo single-frequency ionospheric correction is driven by a single input parameter called the effective ionization level Az (European Commission 2016). It is expressed as where Az is measured in solar flux units or sfu 1 sfu = 10 −22 Wm −2 Hz −1 and the modified dip latitude , or Modip, is defined as tan = I(cos ) −1∕2 . I is the geomagnetic inclination at 300 km height and is the geographic latitude of the location. The three constants in the equation, a i0 , a i1 and a i2 , correspond to the effective ionization coefficients broadcast in the Galileo navigation message. Hence, the Az magnitude is a proxy of solar activity equivalent to the index F10.7 and varies as a function of Modip of the user location.
In addition to these two prominent ionospheric models, also the Neustrelitz Total Electron Content Model (NTCM) has been developed by the German Aerospace Center (DLR) as an alternative approach to mitigate transionospheric time delay and range errors. The NTCM originally presented by Jakowski et al. (2011) is a two-dimensional empirical ionospheric model of the vertical total electron content (VTEC) that defines it as a linear combination of five subfunctions, as follows: where each subfunction, respectively, represents the dependencies on local time, season, geomagnetic field, latitude anomaly and solar activity. In particular, F 5 is driven by a proxy of the solar extreme ultraviolet (EUV) radiation level represented by the solar radio flux index F10.7. Further modifications of this method include the NTCM-BC model (Hoque and Jakowski 2015) that simplifies the original variant by ignoring long-term TEC dependencies (seasonal and solar cycle variability), and the NTCM-Klobpar model (Hoque et al. 2017(Hoque et al. , 2018 which replaces the driving parameter F10.7 by a so-called Klobpar parameter derived from GPS broadcast coefficients. Hoque et al. (2017) have demonstrated that the NTCM-Klobpar model outperforms the classical Klobuchar model by approximately 40% and 10% during high and low solar activity conditions, respectively. Similarly, the NTCM-GlAzpar version has been proposed by Hoque et al. (2019), taking advantage of the applicability of the Galileo broadcast parameters as a solar radiation proxy.
A detailed overview of the equations and variables describing NTCM-based models is provided in their literature, respectively, referred. In summary, the distinction between the different NTCM versions lies in the linear function that defines the solar activity dependence. In the case of NTCM-GlAzpar, this function is expressed as where k 1 = 1.41808 and k 2 = 0.13985 . The quantities k 1 and k 2 are model coefficients computed from an empirical approach fitted to the GlAzpar and worldwide ionospheric behavior derived from TEC data from the Center for Orbit Determination in Europe (CODE, https:// www. aiub. unibe. ch/ resea rch/ code___ analy sis_ center). For details, we refer to Hoque et al. (2019). Since NTCM considers the latitude dependency as a separate function, the solar activity represented by GlAzpar is Modip independent and therefore simplifies the ionospheric correction method for the user. The magnitude of GlAzpar is derived by taking the root mean squared (RMS) of Az values for Modip values between 70° North and South latitude. The square of the Modip Eq. (1) is integrated analytically between the two limits and divided by the Modip range 140, obtaining the following expression in sfu (Hoque et al. 2019) A validation of the NTCM-GlAzpar model in comparison with the NeQuickG model has been given by Hoque et al. (2019) considering VTEC measurements provided by Global Ionospheric Maps (GIMs) of the International GNSS Service (IGS) and altimeter satellites and slant TEC measurements obtained from over 300 IGS stations globally distributed. The model comparison in the VTEC domain shows that the NTCM-GlAzpar model performs, on a global average, at least equal to or better than the NeQuickG model. Also, its computation time has been found to be 65 times faster than NeQuickG, making it a fast running solution for positioning applications.
In this work, for the first time, we present a global validation of the NTCM-GlAzpar model in the position domain over a selected period of perturbed and quiet ionospheric conditions. Moreover, we compare the 3D position performance of the NTCM-GlAzpar, Klobuchar and NeQuickG models by adapting the software capabilities of the GNSS-Lab Tool suite (gLAB) in the SPP mode. For comparison, also the uncorrected SPP solutions are presented. Therefore, the aim of this investigation is to demonstrate that, by implementing small technical changes, the NTCM-GlAzpar algorithm provides a state-of-the-art correction in GNSS-based single-frequency applications.

Data sources
We selected a period of one month over two years with different solar and geomagnetic activity in order to compare the performance of the different ionospheric models. The (3) Aiming at a statistically significant validation of the ionospheric modeling, we achieved a worldwide coverage with up to 73 IGS stations per day in 2014. The number of stations used daily depended on the availability of stationspecific Receiver Independent Exchange Format (RINEX) navigation files and Solution Independent Exchange Format (SINEX) information. For this work, we have used the SINEX products as a priori reference because, by combining independent estimates from a number of IGS Analysis Centers (ACs) and Global Network Associate Analysis Centers (GNAACs), they provide weekly combined solutions consistent at the 1-2 and 3-4 mm levels for horizontal and vertical station coordinates (Ferland and Piraszewski 2009). For simplicity, the station coordinates of the SINEX files were assumed to be constant during the given GPS week. Particularly, the data used in this work correspond to the GPS weeks 1821-1825 and 2082-2086. Referenced by the stations of 2014 and under the same criterion of data availability, the number of stations used for the period of December 2019 resulted in a sample of fewer stations than those used for 2014. The number of stations for December 2019 reached up to 56 stations per day. Figure 2 shows our achieved global coverage (top) and the daily amount of stations with applicable observation, navigation and reference files for this study (bottom).
The GPS observation and navigation files required for SPP processing were downloaded from the NASA's Archive of Space Geodesy Data (https:// cddis. nasa. gov/ archi ve/). The SINEX information used by gLAB as a priori receiver position was also acquired from the data products repository on the same website. In addition, the effective ionization coefficients required for the NTCM-GlAzpar and NeQuickG models were collected from the respective Galileo navigation files for each specific station. During a day, the Galileo space segment broadcasts at least one, but up to four, sets of ionization coefficients that are driving parameters of NeQuickG. In principle, a set of coefficients is stored in the header of the IGS processed broadcast ephemeris product and in the header of the stationspecific navigation data files. However, the Galileo navigation message does not provide an issue-of-data for the ionospheric correction parameters nor does it ensure that identical coefficients are broadcast by all satellites. Moreover, such values are missing for some stations considered in this investigation, and therefore, the mean values retrieved from the rest of available stations were adopted. Figure 1 depicts the resulting mean values of the ionization coefficients a 0 , a 1 and a 2 per day of data and their corresponding standard deviation. The period of December 2014 traced in blue shows a larger deviation between sets of broadcast ionization coefficients, whereas the red profiles of December 2019 show more stable sets of coefficients. In turn, these calculated intervals of a 0 , a 1 and a 2 were used to compute the ionization parameter GlAzpar according to Eq. (4). In the bottom panels of Fig. 2, a comparison of The error bars denote the daily standard deviation calculated from the overall sample of IGS stations. The two bottom panels trace the corresponding estimation of GlAzpar and a comparison to the F10.7 index traced in black. All these magnitudes are expressed in sfu this solar activity driver to its analogous F10.7 proxy (black line) exhibits a good agreement in magnitude for the perturbed period of 2014, with the exception of the days 350-356 that have flux values larger than 160 sfu. The quiet period of 2019 displays a similar profile between the parameters GlAzpar and F10.7, but with absolute values strongly deviated by around 40 sfu. Other authors have already reported discrepancies between these two drivers. For instance, Hoque et al. (2019) pointed out strong deviations between GlAzpar and F10.7 for flux values larger than 160 sfu when comparing them over the period 2013-2017. This discrepancy is nevertheless not expected to degrade the performance of the NTCM-GlAzpar model since its coefficients are generated based on the direct TEC dependency of GlAzpar , and not F10.7.

Data processing with gLAB
For the processing of GNSS data, we used the gLAB software package version 5.4.4 that can be downloaded from https:// gage. upc. edu/ gLAB/ (Ibañez et al. 2018). gLAB has been developed as a multipurpose process and analysis tool under a European Space Agency (ESA) contract with the research group of Astronomy and Geomatics (gAGE) from the Universitat Politecnica de Catalunya (UPC). Our motivation to use gLAB as a positioning tool for the validation of the ionospheric model NTCM-GlAzpar adheres to its extended area of applicability within ESA's educational and professional context. Validation studies of gLAB have confirmed that this tool keeps at least the same level of accuracy as the GPS Toolkit (GPSTk, Salazar et al. 2010) maintained by the Applied Research Laboratories of the University of Texas, meaning that it is capable of achieving a 3D position accuracy better than 5 cm for Precise Point Positioning (PPP) experiments ). gLAB's customization capabilities offer an interactive Graphical User Interface (GUI) and a Data Processing Core (DPC) that can be adapted to the user's requirements. Indeed, although current versions of this software support ionospheric correction with the Klobuchar model for GPS data and NeQuickG model for Galileo data, we implemented further changes into its DPC to include the NTCM-GlAzpar algorithm. Our modifications to gLAB also made possible to use NeQuickG and the three broadcast ionization coefficients in Galileo data to apply them as driving parameters to GPS data. The NeQuickG version 1.2 used for this investigation was obtained from the European GNSS Service Center (GSC) and is publicly available at https:// www. gsc-europa. eu/ suppo rt-to-devel opers/ nequi ck-g-source-code.
To start a standard processing routine, first, the stationspecific RINEX observation and broadcast navigation files, as well as the SINEX a priori information, were uploaded to the gLAB execution command as input data. For the data modeling, we set a 30 s resolution decimation and an elevation mask for satellites below 10˚ since they may contain increased errors due to low signal-to-noise ratio and multipath. For the rest of the processing parameters, we attached to the default template offered by gLAB, meaning that, between others, the tropospheric mapping was computed as the obliquity factor described by Black and Eisner (1984). This mapping depends only on satellite elevation and it is common for wet and dry components. The default SPP modeling settings also corrected the Differential Code Biases (DCB) due to electronics, antennas and cables of receiver and transmitter devices by using the P2-P1 bias given as the Total Group Delay (TGD) in the RINEX navigation file. The processing was performed by supposing a kinematic receiver and by uniquely using the pseudorange measurement of the GPS L1 C/A (coarse acquisition) code signal for the filtering.
Of particular significance for our scientific purpose was the configuration of the ionospheric correction model to be used. This parameter allowed us to evaluate the positioning performance of four different modeling cases: As a result of the processing with gLAB, a file containing the receiver daily solution message in the Earth Centered Earth Fixed (ECEF) coordinate system was created. As a rule, our decimation of 30 s meant that the output file contained 2880 modeled values per day (86,400 s divided by 30). We used the fields containing the receiver North ( Δ iN ), East ( Δ iE ) and Up ( Δ iU ) difference in relation to the SINEX nominal a priori position to determine the 3D position error according to the expression where i denotes each of the 2880 temporal points over the day. Subsequently, we calculated the hourly mean Δ 3D hour over each of the periods DOY 335-365 of 2014 and 2019 by grouping the processing results of each available IGS station. Figures 3, 4 and 5 display the modeled hourly mean position errors of the 3D, North, East and Up coordinates for stations grouped by high-, mid-and low-latitude ranges, respectively. These position errors as a function of local time are given for both, the perturbed conditions of 2014 (top panels) and the quiet period of 2019 (bottom panels). Each subplot traces the results obtained with no correction model

Evaluation of models and discussion
The resulting mean time series presented in Figs. 3, 4 and 5 represent well the tendency observed among the whole sample of up to 73 stations. We can describe these position error profiles in terms of temporal and spatial effects. Independently of the latitude region and the period of solar activity, in general, the North and East position errors are hardly affected by the application of an ionospheric correction model, reaching positioning fluctuations of a couple of meters along the diurnal LT and being such errors slightly greater for the perturbed period than for the quiet conditions. On the other hand, the error estimation of the Up component is strongly diminished by the application of the ionospheric correction models and is also accountable for the large deviations of the 3D position errors. The perturbed solar activity of December 2014 causes uncorrected 3D range errors of about 8 m at day-time for high-latitude stations and they may increase to more than 20 m for stations at low-latitudes. The quiet period of December 2019 creates little 3D error position variability for high-latitude stations with flat values of around 2.5 m, but for equatorial regions, an error peak at day-time may reach some 7 m. The 3D position errors as a function of local time for every latitude range display a typical Gaussian-shaped profile, except for the 3D curve of the equatorial stations under perturbed solar activity (top panel of Fig. 5). We see a secondary peak at post-sunset LT that causes error deviations of more than 5 m for modeled signals. Such a feature can be explained in terms of post-sunset irregularities that characterize the equatorial ionosphere (Merid et al. 2021). At first glance, the results of Figs. 3, 4 and 5 allow us to resolve that the NTCM-GlAzpar (black solid lines) and NeQuickG (blue dashed lines) models outperform the Klobuchar algorithm (green dash-dotted lines) and the case with non-corrected ionospheric delay (red dotted lines). This is particularly evident from the 3D positioning solutions for both perturbed and quiet solar activity. Compared to the case without the ionospheric correction model, the Klobuchar model also greatly improves the 3D error position deviation. However, we notice that the rest of modeling approaches reach their peak at around 14 h LT for high-and mid-latitude ranges, whereas the peak of the Klobuchar model is shifted to about 18 h LT (Figs. 3  and 4). The reason for this deviation has been attributed to strong perturbations, and therefore, this is not seen for the analysis of the quiet period (Hoque et al. 2018).
For a better perception of the model-specific performance along the worldwide distribution of stations and diurnal variation, we followed a comparative analysis similar to the evaluation of the NTCM-Klobpar ionospheric model presented by Hoque et al. (2018). Namely, with the a priori station coordinates of the SINEX files as a reference, we computed the RMS, mean, standard deviation, 65 and 95 percentiles of the 3D position errors obtained with the Klobuchar, NTCM-GlAzpar and NeQuickG models. The statistical estimates were derived separately for perturbed and quiet periods. To obtain statistically significative results, the data were arranged based on geodetic latitudinal position of the stations ( ) and LT range into the following groups: In this manner, considering the overall number of 73 IGS stations (groups i, ii and iii), these newly established sectors enclosed 34 stations at low-latitudes (group iv), 30 at midlatitudes (group v) and 9 at high-latitudes (group vi). The statistical analysis for each group and ionospheric model is presented in Table 1, except for the case with no application of an ionospheric correction model since it shows the poorest performance. As follows, we discuss and evaluate the performance of the ionospheric correction models for each group of this statistical experiment: i. Complete geographical and diurnal dataset: The overall ionospheric delay mitigation of the NeQuickG and NTCM-GlAzpar models in comparison with an approach without correction or the Klobuchar algorithm is evident, especially for the lapse of 2014. By comparing the two models driven by Galileo ionization coefficients, NTCM-GlAzpar performs better than NeQuickG with RMS values of 4.36 m versus 4.61 m during perturbed activity. Their performance is nevertheless on the same level for quiet conditions, with, respectively, 2.32 m versus 2.35 m RMS values. Table 1 compiles the rest of statistical metrics for this sector, and the first panel of Fig. 6 also graphically compares the RMS 3D position error achieved by the Out of this statistical experiment considering six different subsets of data for perturbed and quiet solar activity, we can clearly state that, in the 3D position domain, the ionospheric correction model NTCM-GlAzpar performs better than the NeQuickG, Klobuchar and non-corrected cases on a global scale (group i). Nevertheless, the results are ambiguous when comparing the two best performing models NTCM-GlAzpar and NeQuickG for horizontal and vertical directions independently. For the perturbed ionospheric conditions of 2014, NeQuickG performs better than NTCM-GlAzpar for correcting the horizontal components (RMS values of 2.53 m versus 2.72 m, respectively), but worse for vertical corrections (RMS values of 3.86 m versus 3.41 m). We see the opposite case for the quiet period of 2019. While NTCM-GlAzpar performs better than NeQuickG for the correction of horizontal position errors (respective RMS values of 1.03 m versus 1.21 m), its performance is slightly poorer for correcting vertical deviations (RMS values of 2.08 m versus 2.02 m). Also, the performance of the Neustrelitz model deteriorates toward higher latitudes. This is evident from the 3D position error histograms for each latitude Fig. 6 Graphical comparison between RMS values of the 3D position error achieved with the four modeling approaches for perturbed and quiet periods. The panels display the results for subsets organized by local time and latitude range range in Fig. 7. In comparison with the rest of modeling approaches, the probability density function of the NTCM-GlAzpar model has a better mode for stations in the latitude sector |0 • − 30 • | , but it is slightly surpassed by the mode of NeQuickG toward stations with latitudes between |60 • and90 • |.
Our analysis of datasets grouped by local time (groups ii and iii) shows that the 3D position estimation through the NTCM-GlAzpar model also deteriorates during nighttime. Therefore, the fact that the accuracy of NTCM-GlAzpar is poorer at night and higher latitudes is linked to the reduced solar activity conditions. With the decreasing of solar activity and irradiation intensity, and during perturbed conditions, the Az parameters which are optimized for NeQuickG pronounce more and more the characteristics of NeQuickG, thus reducing the influence of the solar activity level which is the external driver for NTCM. This is seen in Fig. 1, where F10.7 and Az derived parameters differ significantly under low solar activity conditions. In other words, if the Az parameters were optimized for NTCM-GlAzpar in the same way as for NeQuickG, further improvements in the positioning performance using the Neustrelitz model would be expected.
A poorer positioning performance during periods of low solar activity is true for most of ionospheric empirical models, and also, the correlation between TEC and the F10.7 proxy is reduced. Due to this reason, the quality of the ionization coefficients a 0 , a 1 and a 2 estimated in the Galileo control segment by comparing monitoring stations TEC to model outputs is degraded, resulting in disparities between the GlAzpar and F10.7 parameters as observed for the period of December 2019 (bottom panel of Fig. 2 in red). Indeed, this is also the reason why the benefit of NTCM-GlAzpar, on a global scale, is imperceptible for the quiet period of December 2019 when compared to NeQuickG.

Conclusions
We have performed a comprehensive evaluation of the positioning accuracy obtained by four ionospheric correction modeling cases: no ionospheric correction model, Klobuchar model, NTCM-GlAzpar model and NeQuickG model. We used GNSS data of up to 73 stations worldwide corresponding to a period of perturbed (December 2014) and quiet (December 2019) solar and geomagnetic activity. The data processing and determination of 3D positioning solutions were achieved by utilizing the analysis tool gLAB and its customization capabilities in the SPP template.
As summary of this novel validation of the NTCM ionospheric model driven by Galileo ionization coefficients, we can state that, in the position domain, the NTCM-GlAzpar clearly outperforms the Klobuchar model and slightly surpasses the accuracy of NeQuickG model on a global scale. When comparing the 3D position error results for the whole data sample, the RMS values achieved with NTCM-GlAzpar are of 4.36 m and 2.32 m versus 4.61 m and 2.35 m obtained with NeQuickG for perturbed and quiet conditions, Fig. 7 Histograms of the probability density functions of the four ionospheric modeling approaches for the periods of perturbed (top panels) and quiet (bottom panels) solar activity. The comparison is displayed for the entire sample of stations, as well as for stations grouped by latitude ranges. The bin width of the histograms is set to respectively. It means an improvement in the order of 25 cm and 3 cm for the respective investigated periods. Nevertheless, a statistical analysis per latitude and LT subsets allows us to identify that the accuracy of NTCM-GlAzpar degrades toward high latitudes, at nighttime and for periods of low solar activity. Although the 3D position estimation of NTCM-GlAzpar can be established on a comparable level of accuracy with NeQuickG, its faster computational benefit of around 65 times and simple technical adaptations make the Neustrelitz model an alternative practicable algorithm to improve the accuracy of GNSS single-frequency applications.