Assessing the quality of ionospheric models through GNSS positioning error: methodology and results

Single-frequency users of the global navigation satellite system (GNSS) must correct for the ionospheric delay. These corrections are available from global ionospheric models (GIMs). Therefore, the accuracy of the GIM is important because the unmodeled or incorrectly part of ionospheric delay contributes to the positioning error of GNSS-based positioning. However, the positioning error of receivers located at known coordinates can be used to infer the accuracy of GIMs in a simple manner. This is why assessment of GIMs by means of the position domain is often used as an alternative to assessments in the ionospheric delay domain. The latter method requires accurate reference ionospheric values obtained from a network solution and complex geodetic modeling. However, evaluations using the positioning error method present several difficulties, as evidenced in recent works, that can lead to inconsistent results compared to the tests using the ionospheric delay domain. We analyze the reasons why such inconsistencies occur, applying both methodologies. We have computed the position of 34 permanent stations for the entire year of 2014 within the last Solar Maximum. The positioning tests have been done using code pseudoranges and carrier-phase leveled (CCL) measurements. We identify the error sources that make it difficult to distinguish the part of the positioning error that is attributable to the ionospheric correction: the measurement noise, pseudorange multipath, evaluation metric, and outliers. Once these error sources are considered, we obtain equivalent results to those found in the ionospheric delay domain assessments. Accurate GIMs can provide single-frequency navigation positioning at the decimeter level using CCL measurements and better positions than those obtained using the dual-frequency ionospheric-free combination of pseudoranges. Finally, some recommendations are provided for further studies of ionospheric models using the position domain method.


Introduction
In the case of the radio waves transmitted from the global navigation satellite system (GNSS), the Total Electron Content (TEC) of the ionosphere delays pseudorange measurements and advances carrier phase measurements by the same amount (Parkinson et al. 1996). Stand-alone user receivers operating with more than one frequency can eliminate 99.9% of ionospheric refraction by using the so-called ionosphericfree (IF) combination ( L IF = f 2 1 L 1 −f 2 2 L 2 f 2 1 −f 2 2 ) of carrier phase measurements L, at two frequencies f 1 and f 2 . By contrast, singlefrequency receivers must select one of several ionospheric models (Rush 1986 and references therein) to account for the Slant TEC (STEC) as accurately as possible. Therefore, for this type of user, it is important to assess the error of the ionospheric delay correction model because the uncorrected, i.e., unmodeled, part of ionospheric delay contributes to the absolute positioning error of GNSS-based applications.
Several tests are available to assess the accuracy of global ionospheric models (GIMs). Usually, these assessments rely on a comparison of the ionospheric model prediction, 4 Page 2 of 12 STEC GIM , against an ionospheric truth, STEC ref . Examples of such reference STECs are measurements from dual-frequency space-borne radar altimeters (Imel 1994;Orús et al. 2003). However, errors in such reference values can reach up to several total electron content units (TECUs), where 1 TECU = 10 16 e − /m 2 and corresponds to 16 cm at the f 1 frequency. Therefore, the applicability of these reference values is limited for assessment of GIMs to within several TECUs of error. Rovira-Garcia et al. (2016b) presented an assessment of GIMs based on reference values accurate to the centimeter by means of an unambiguous and unbiased geometry-free (GF) combination of carrier phase measurements L GF = L 1 − L 2 . Such precise reference values were obtained as STEC ref = L GF − B GF , where the carrier phase ambiguities in the GF combination, B GF = N 1 − N 2 , are solved using a network solution (Blewitt 1989). Such integer ambiguity resolution (IAR) of the carrier phase ambiguities N 1 and N 2 (Banville et al. 2012) exploit Precise Point Positioning (PPP) techniques (Malys and Jensen 1990;Zumberge et al. 1997) to obtain the carrier phase measurements as precise and unambiguous pseudoranges. Despite its accuracy, IAR as part of PPP requires a complex data process where satellite orbits, clocks, and tropospheric delays are estimated at the centimeter level, which limits the applicability of PPP IAR.
Instead of using IAR, a more straightforward strategy to estimate the carrier phase ambiguity term is the so-called Carrier phase to Code Levelling (CCL) process (Mannucci et al. 1998). The carrier phase measurements L GF are "leveled" to the GF combination of code pseudoranges, P GF = P 2 − P 1 , by averaging the difference between them for each continuous arc of the samples. Although code multipath produces leveling errors (Brunini and Azpilicueta 2009), the ionospheric references obtained as STEC ref = L GF − ⟨L GF − P GF ⟩ exhibit an intermediate accuracy that ranges from 1.4 to 5.3 TECUs (Ciraolo et al. 2007;Rovira-Garcia et al. 2016a) without the use of a complex data handling process.
The STEC correction can be evaluated with the so-called self-consistency test (Orús et al. 2005), in which one carrier phase ambiguity per arc is estimated using least squares (LS) from the differences between L GF and STEC GIM . The residuals of the LS fit are the metric for GIM assessment. The advantage of the self-consistency test is its simplicity and its precision. However, its drawback is that it only accounts for the standard deviation of the error; regional biases present in the GIM are neglected because they are absorbed in the arbitrary ambiguity value. In this regard, the methods which previously estimate carrier phase ambiguities, such as IAR or CCL, are more robust for the testing.
Another testing methodology is based on computing position, velocity and time (PVT) of permanent fixed stations with known coordinates and comparing single-frequency positioning errors obtained with different GIMs. The PVT approach to study the accuracy of GIMs is straightforward, as it does not require the complexity of computing reference ionospheric values. However, it faces several difficulties that, if not properly considered, can lead to inconsistent results compared to tests on the ionospheric delay domain, based on the difference between STEC GIM and STEC ref . This means that the comparison of the accuracy among different GIMs in the position domain can be contaminated by errors other than the ionosphere. For example, when using the Standard Point Positioning (SPP) method in Wu et al. (2013) or Hoque et al. (2015), two major error contributions are code pseudorange measurements (thermal noise and multipath of approximately 1 m or greater), and broadcast satellite orbits and clocks (errors at meter level). Therefore, SPP is limited to the assessment of GIMs within several TECUs of error.
In Orús (2017) and Nie et al. (2018), the error contribution from satellite orbits and clocks was mitigated by using precise products from the International GNSS Service (IGS) that are accurate to a few centimeters (Beutler et al. 1999;Dow et al. 2009). Despite this improvement, some inconsistent results remain when comparing to the conclusions of the PVT assessment presented in Orús (2017) with the assessment based on STEC ref values presented in Rovira-Garcia et al. (2016b). In particular, such inconsistent results surfaced for the GIM of the Fast Precise Point Positioning (Fast-PPP) technique, computed using about 150 permanent stations distributed worldwide (Rovira-Garcia et al. 2015), which exploits its accuracy to reduce the convergence time of the PPP navigation solution. In Orús (2017), the GIMs computed by IGS performed similarly to the Fast-PPP GIM. In contrast, in Rovira-Garcia et al. (2016a, b) the Fast-PPP GIM outperformed the IGS GIMs by a factor 2 to 3. This inconsistency has motivated us to use the same data set to compare the same GIMs that are used in Orús (2017). Overall, it is of great interest to investigate whether the accuracy of GIMs can be discriminated using the PVT assessment in a similar manner than using STEC ref values that require complex geodetic processing.
The goal of this contribution is to analyze the reason for obtaining different results when applying the ionospheric reference values and position domain methodologies. For this purpose, some of the remaining error contributions, such as measurement noise, multipath, and outlier presence, are identified. Some recommendations are provided to perform further tests because the single-frequency PPP assessment of ionospheric models is receiving increasing attention from the ionospheric research community (Prol et al. 2018).

Data set
To reproduce the results of Orús (2017), we analyzed the Global Positioning System (GPS) measurements of the entire year 2014 collected by a global distribution of 34 permanent stations (Table 1). The observation Receiver INdependent EXchange format (RINEX) files were downloaded from the Multi-GNSS Experiment (MGEX) network (Montenbruck et al. 2017).
To mitigate the errors in the GPS satellite orbits and clocks, we used the Final Combined IGS Products, at a sampling rate of 900 s for orbits and 30 s for clocks. The Antenna Phase Center (APC) corrections for satellite and receivers were obtained from ANTenna EXchange format (ANTEX) files of the corresponding GPS week.
The satellite hardware delay for the P1-P2 Total Group Delay (TGD) was corrected from the GPS broadcast navigation message. This is consistent with the testing methodology presented in Orús (2017), where the GPS broadcast TGDs substituted those computed for Fast-PPP to bring it to the same level as GPS and Galileo regarding TGD broadcasting. Although any GIM should be used in conjunction with their associated TGDs instead of the broadcast values, the agreement among satellite DCBs computed by different IGS Analysis Centers is at the level of a few tenths of nanosecond (Hernández-Pajares et al. 2009). Then, the mismodeling of using the broadcast TGDs is several times lower than the expected errors in the ionospheric delays of the GIMs under test that are in the range of several TECUs (Rovira-Garcia et al. 2016b). Therefore, using broadcast TGDs for all GIMs does not distort the PVT test and all GIMs are affected similarly.
The present comparison includes the same GIMs as in Orús (2017): the Fast-PPP GIM and two GIMs from IGS, namely the Final Combined IGS Product (IGSG, for short) and the Rapid Product from the Universitat Politècnica of Catalunya (UPC) (UQRG, for short). The STEC correction of these GIMs was obtained following the IONosphere map EXchange format (IONEX) standard defined in Schaer et al. (1998). Every IONEX file contains two types of maps for 24 h and with a fixed interval of 15 min (Fast-PPP and UQRG) and 2 h (IGSG). The first set of maps contains the vertical TEC (VTEC) at Ionospheric Grid Points (IGPs), whereas the second set of maps contains the root-meansquare (RMS) error of the VTEC at each IGP. Although the IONEX standard names it as VTEC RMS error, those values are the diagonal elements of the covariance matrix of the VTEC estimation at each IGP. Therefore, the value depends not only on GNSS measurements but also on the constraint equations imposed to fill data holes and to smooth the VTEC values. Assuming a thin-shell model of the ionosphere, we projected the VTEC and its RMS error to slant values according to the secant of the zenith angle at a mean ionospheric height (Parkinson et al. 1996).
The tested GIMs exhibit some differences. First, IGS GIMs use one single layer at an altitude of 450 km to describe the VTEC, whereas the Fast-PPP GIM uses two layers at altitudes of 268 and 1600 km. Second, IGS GIMs were computed with stations covering all longitude ranges, whereas the original Fast-PPP GIMs used in Orús (2017) involved only stations within the longitude range of − 130° to + 130° (Rovira-Garcia et al. 2016a). This was because the original coverage requirement for the Fast-PPP GIM included the longitudes within ± 110°. Since this difference is important, the figures in the manuscript indicate longitudes outside ± 130° with a gray shadow. To overcome this limitation, we have also used recomputed Fast-PPP GIMs, including all longitudes. The third difference resides in the RMS errors of the VTEC of the three GIMs, as depicted in the example shown in Fig. 1. The RMS values of IGSG (top left) and UQRG (top right) significantly differ, despite the use of similar satellite geometries since both GIMs were derived from the GPS constellation. In addition, the RMS errors of UQRG and IGSG are more uniform than those in the Fast-PPP (bottom panels). The Fast-PPP GIM was derived with lowsmoothing constraints because it targeted the retrieval of STEC in a realistic manner under conditions of high ionospheric activity around the Solar Maximum (Rovira-Garcia et al. 2016a). This required the use of a large process noise at low latitudes, resulting in the large VTEC RMS values shown in Fig. 1 for the Fast-PPP GIM for both the original (bottom left) and the recomputed version (bottom right).
In the original computation of the Fast-PPP GIM, the process noise for the IGPs in the southern hemisphere was reduced to compensate for fewer available receivers. In contrast, this constraint was not imposed for the recomputed Fast-PPP GIM, then the process noise does not depend on the receiver density but on the expected ionospheric variability: large process noise is used for IGPs at low latitude and high latitude. In this sense, the recomputed VTEC RMS values in the southern hemisphere are driven mainly by the data because the constraints are down-weighted, which is a more realistic assumption. This results in low values of VTEC RMS located only at IGPs sounded by the reference stations, see the dashed lines in the bottom right panel of Fig. 1.
GIM users can only know how good the ionospheric model has been solved by noticing the magnitude of the RMS error of the VTEC. In general, it is expected that this value is low in well-sounded areas because many receivers are employed to generate the ionospheric correction. Conversely, the RMS error of the VTEC increases with the distance to the reference stations, with the rate of increase being a function of the degree of smoothing that is introduced by means of constraint equations, applied in the generation of the GIMs.

Methodology
We have executed the test defined in Rovira-Garcia et al. (2016b). The method uses LS to fit the differences between STEC GIM and STEC reference for all receivers and satellites during 24 h, to a daily constant per receiver and a daily constant per satellite: where STEC reference is the unambiguous L GF obtained after performing IAR for all measurements in the network. The metric of the test is the postfit residuals of the estimated K values obtained with LS. The test is executed for every day in 2014 and the postfit residuals are accumulated to have one distribution for every GIM.
We used our in-house open-source GNSS-Lab Tool suite (gLAB) (Ibáñez et al. 2018) to compute the PVT solution. The default PPP configuration in gLAB was slightly modified to be consistent with the previous assessment (Orús 2017) as follows: sampling rate of 30 s, disabling of the carrier phase measurements, using elevation mask of 7°, including satellites under eclipse, using nominal tropospheric correction from satellite-based augmentation systems (SBASs) (RTCA 2016), applying a threshold of 30 m to filter outlier values in the pre-fit residuals over the pre-fit residual median, and disabling the use of the VTEC RMS errors present in the GIM as a weight to the pseudorange measurements: where 0 is the standard deviation of the measurements, assumed 1 m for all satellites as in (Orús 2017), and GIM is the slant value of the RMS error of the VTEC extracted from the IONEX files. As in Orús (2017), the study is done assuming GIM = 0 ; however, since the latter is an important part of the GIMs, results using non-zero GIM are presented in parenthesis next to the results using only 0 .
The first step in the methodology is to compute the reference position to be used as ground truth. Note that precise coordinates for MGEX receivers are available starting in 2015 (Montenbruck et al. 2017), whereas data collection started in 2014. Thus, the daily reference coordinates were estimated with PPP using L IF . After 24 h of data collection and coordinates processed in static mode, the typical accuracy is at the centimeter level (Sanz et al. 2013).
The second step is the PVT computation using single-frequency f 1 . At every epoch, the positioning error of every GIM is evaluated as the difference between the reference position computed in the first step and the coordinates estimated in kinematic mode, with the position coordinates modeled as white noise. We have used two types of measurements for the PVT test: the pseudorange code measurements in order to be consistent with Orús 2017), and CCL measurements to reduce the code noise. Indeed, once the noise of the input measurements has been filtered, the most important remaining source of error is the ionospheric mismodeling. Thus, the accuracy of the ionospheric correction can be properly characterized. Finally, to distinguish PVT errors caused by weak satellite geometries as expressed by high Dilution Of Precision (DOP), the dual-frequency IF solution is also computed using both pseudoranges and CCL measurements. Using CCL measurements requires two passes over the same data set. In the first pass, the carrier phase ambiguities are determined as the average of the difference between P GF and L GF , per each continuous arc of all samples. In the second pass, we compute the PVT in the navigation filter using CCL measurements as if they were precise pseudoranges. In this manner, carrier phase ambiguities are not estimated in the navigation filter simultaneously with the PVT parameters, which is unlike standard single-frequency PPP in which case the estimated ambiguities absorb part of the mismodeling present in the GIM predictions. Therefore, the error in the ambiguity estimation is translated into the PVT, thus impeding a thorough calibration of the GIM accuracy by means of the PVT error.
Steps 1 and 2 are repeated for each day and for each station. The third and final step is to determine the RMS, mode, mean and percentiles (50th, 68th, 95th) of the positioning error distribution. As we will show later, the distribution of 3D positioning errors is not homogeneous. Thus, reflecting the error by percentiles is better than by simple RMS. For instance, the 95 th percentile quantity is the value where 95% of the 3D positioning errors are contained, after having sorted all errors in ascending order.

Results
The accuracy of the three GIMs has been assessed on the ionospheric delay domain, i.e., using STEC ref , and on the position domain using pseudorange and CCL measurements. The comparison between both test domains allows detecting several difficulties that, if not properly considered, can lead to inconsistent results. The present analysis, which uses the percentiles of the distribution functions of the postfit residuals in (1) and of the 3D positioning errors, complements previous assessments based on the use of the RMS that have been reproduced. Namely, the results of the test defined in Rovira-Garcia et al. (2016b) that rely on precise ionospheric reference values and the results of the single-frequency pseudorange (P1 in short) positioning test presented in Orús (2017).  Table 2 shows the results of the ionospheric delay domain test based on (1) and using the IGSG, UQRG and Fast-PPP GIMs. It shows that the original two-layer Fast-PPP GIM reduces the RMS of the postfit residuals of the single-layer GIMs IGSG and UQRG by 57% and 43%, respectively, which is in line with the previous assessments by Rovira-Garcia et al. (2016a, b). Moreover, the recomputed Fast-PPP GIMs, which include all longitudes, improved the 95th percentile and the RMS. Figure 2 depicts the 3D positioning errors, with respect to the reference coordinates obtained with static PPP, using pseudorange measurements. We assess first the results of the IF solution because these should not depend on the receiver location. However, the obtained positioning error percentiles are heterogeneous. Indeed, stations equipped with Javad receivers exhibit errors larger than the other two manufacturers. Although multipath and thermal noise depends on the receiver hardware and antenna environment, this dependency interferes with the test objective of assessing the accuracy of the different GIMs through the GNSS positioning. Figure 2 depicts that the original Fast-PPP GIM produced larger 3D positioning errors than the IF at four stations, as previously acknowledged in Orús (2017). Specifically, at the three stations located at the most western longitudes (FTNA, MAO0 and THTG) and at station SGOC. The origin of the large positioning errors is the same in these four receivers located at latitudes close to the geomagnetic equator, where the ionospheric gradients are important: the nearest reference stations used to derive the original Fast-PPP model were located at large distances causing large RMS errors in the VTEC.

Position domain results with pseudoranges
As commented before, we recomputed the Fast-PPP GIMs for 2014. We added eight reference stations to improve coverage of the Pacific ocean, and we increased the process noise of the recomputed Fast-PPP GIMs to provide RMS errors of the VTEC in a more realistic manner. For the receiver SGOC, despite being selected as a reference station, an error parsing the header of its RINEX files prevented its inclusion into the generation of the original Fast-PPP GIM. Thus, the nearest stations to SGOC were located at distances of 723 and 1772 km. This problem has been solved in the recomputed Fast-PPP GIMs. Table 3 presents the 3D positioning errors accumulating all stations in 2014. The table provides a more informative description than the use of only the RMS or mean values, as discussed later. 3D positioning errors are obtained assuming GIM = 0 in (2), as in (Orús, 2017) and the values in parenthesis correspond to 3D positioning errors obtained using GIM to weight the GIM corrections in (2). The 3D positioning error decreases for all GIMs when its GIM is used, confirming that the RMS errors of the VTEC depicted in Fig. 1 are an essential part of the GIM that should be used.
Numerically, the RMS values for the single-and dualfrequency 3D positioning errors are similar to those reported in Orús (2017) and Nie et al. (2018). Therefore, the outcome of single-frequency P1 positioning using precise ionospheric models is similar to Orús (2017); accurate GIMs can produce navigation errors smaller than those obtained with the dual-frequency IF code combination. The improvement occurs if the GIM is accurate enough, i.e., the mismodeling present in the GIM is smaller than the amplification factor of the pseudorange noise produced by the IF combination.
The two rightmost columns of Table 3 summarize the improvement of the recomputed Fast-PPP GIMs. The recomputed Fast-PPP GIMs improved the 3D positioning errors at the four stations mentioned above, whereas the results were similar for the remaining 30 stations. Because it concerns errors in the range of several tenths of meters, only the largest percentiles of the positioning error distribution and the RMS are clearly distinct in comparison with the original Fast-PPP GIMs used in Orús (2017). On the contrary, the mode does not vary, and the median is similar. As indicated in Orús (2017), outlier values affect the results of the RMS assessment.

Assessing the use of RMS as a metric
This subsection addresses the effect of accumulating the positioning errors into a single indicator. In particular, we analyze the suitability of using the RMS of the 3D navigation errors as a widely used metric to compare the accuracy of GIMs. Figure 3 shows the distribution of the 3D positioning errors for the IF solution and for the three GIMs. The top panel shows the entire distribution, whereas the bottom panel shows errors smaller than 5 m. In the top plot, the analysis of the tails of the error distribution reveals that for all four processing modes, large errors play an important role   Table 1 for 2014. The bin size is 3 cm. The top panel corresponds to the entire distribution, and the bottom panel corresponds to errors smaller than 5 m. The labels correspond to the same products as those in Fig. 2 4 Page 8 of 12 in the calculation of the RMS values of Table 3, as positioning errors of one order of magnitude difference are mixed after being squared.
In the case of the dual-frequency IF solution, large 3D positioning errors are caused by multipath and satellite geometries with large DOP. In addition to these errors, the single-frequency P1 solutions are affected by errors in each GIM, mainly determined by the topology of the network of stations used to derive the ionospheric correction that, as shown in Fig. 1, is neither dense enough nor homogeneous, due to practical reasons; e.g., oceans.
The bottom panel of Fig. 3 shows where the individual distribution of 3D errors is smaller than 5 m. These are cases referring to small DOPs and well-sounded areas, where the GIM corrections are reliable. In such a condition, the GIMs can be thoroughly tested. The peaks of the positioning error distributions, which correspond to the mode values in Table 3, are clearly observed. An important result arises from the comparison of the modes, that is, using the same data and strategy as used in Orús (2017): in well-sounded regions, the original Fast-PPP GIM reduced the mode values of the 3D errors obtained with UQRG and IGSG by 30% and 40%, respectively. In line with the results of Table 2, the mode value of the 3D positioning errors obtained with the recomputed version of the Fast-PPP GIMs is the same as the original Fast-PPP GIMs In terms of the RMS values of 3D positioning errors, the original Fast-PPP GIM values are 8% lower than those of IGSG and 13% larger than the values of UQRG. On the contrary, the recomputed Fast-PPP GIM improves the RMS of the 3D positioning errors obtained with IGSG and UQRG by 37% and 23%, respectively. This inconsistency is related to the small number of available test receivers in poorly sounded areas outside longitudes ± 130°, highlighted in gray in Fig. 2, that produced large positioning errors, as already observed in Orús (2017).

Position domain results with CCL
The positioning performance using pseudorange measurements observed for the recomputed Fast-PPP GIM in comparison to the IGS GIMs is still far from the results found in Table 2, where Fast-PPP GIM reduces the RMS errors of UQRG and IGSG by 53% and 65%, respectively. In this subsection, we analyze the positioning errors obtained with CCL measurements which are of lower noise. Figure 4 illustrates how the 68th and 95th percentiles of the 3D positioning errors decrease in stations where the code noise and multipath are the dominant errors. CCL positioning errors are more homogenous than those obtained with code pseudoranges, thus reducing the receiver-type dependency. We observe that for all stations, the 95th percentile of the 3D positioning error for the dual-frequency IF is at the level of a few decimeters, which is only slightly degraded with respect to the standard kinematic PPP floating the carrier phase ambiguities. Thus, the CCL procedure functions as intended, and only a small error in the CCL alignment is absorbed by the positioning. Obviously, ambiguities estimated with CCL are worse than those of IAR, but they exhibit enough quality to calibrate the accuracy of the GIMs.  Figure 4 depicts how in poor-sounded areas, e.g., stations JOG2 or WUH2, the CCL errors increase compared to Fig. 2. This worsening is attributable to a smaller number of satellites used to compute the PVT solution, after the necessary cycle-slip detection for the CCL process.
The effect of the pseudorange noise in position-based test methodologies can be inferred by comparing the 3D positioning errors obtained with raw pseudorange measurements and with the CCL measurements: graphically from the Cumulative Distribution Functions (CDFs) shown in Fig. 5 and numerically from the mode values of Tables 3 and 4. The Fast-PPP peak value is reduced by 50% (from 54 to 24 cm), whereas the reduction is more limited for the other GIMs: 20% for IGSG (from 90 to 75 cm) and 12% for UQRG (from 78 to 69 cm). Therefore, the pseudorange noise contributed to the positioning errors with Fast-PPP GIM in an appreciable manner, whereas, for the IGS GIMs, the most important source of error was the ionospheric mismodeling.
The reduction of the noise of the input measurements allows inferring the resolution of the position domain test: pseudorange measurements could be sufficient to assess IGS GIMs, whereas testing accurate GIMs as Fast-PPP requires precise CCL measurements.
The top panel of Fig. 6 depicts the entire distribution of positioning errors using CCL measurements, including large errors obtained in poorly sounded areas. Fast-PPP GIM provides smaller positioning errors than the other two GIMs, despite UQRG and IGSG present smaller RMS errors of the VTEC in Fig. 1 than those of Fast-PPP GIM. The RMS errors of the VTEC present in GIMs tailored for navigation should bound as confidently as possible the actual errors. This can be achieved by a careful selection of the process noise and the constraint equations used in the VTEC estimation process by means of Kalman filter, for instance, to fill data holes.
In contrast, the bottom panel of Fig. 6 depicts 3D positioning errors obtained in nominal conditions were all GIMs can be homogeneously tested: well-sounded areas, good DOPs and low measurement noise. The peak comparison, i.e., the mode, reveals that Fast-PPP reduces the mode values of IGSG and UQRG by 68% and 65%, respectively. This improvement is in line with the results based on reference on ionospheric reference values presented in Table 2 and previous assessments (Rovira-Garcia et al. 2016a, b).

Summary and conclusions
Evaluation changes in estimated receiver coordinates (position domain approach) is a straightforward method for comparing the accuracy of GIMs. However, this method presents several difficulties that, if not properly considered, can lead to inconsistent results compared to the traditional ionospheric delay domain tests. We have identified the following  Table 1 for 2014. The top panel corresponds to errors obtained with code pseudorange measurements and the bottom panel corresponds to those obtained with CCL measurements. The labels correspond to the same products as those in Fig. 2  issues that should be considered when performing a thorough test by means of the GNSS positioning error: 1. Sparsity of data The various GIMs are computed using heterogeneous receivers and different distributions of permanent receivers. Therefore, accuracies of VTEC present in the GIMs are not homogeneous. This is accounted for in the RMS errors included in the IONEX files, which are an essential part of the ionospheric correction that should be taken into account. In this sense, comparative assessments should use receivers situated in well-sounded regions where the RMS errors of the VTEC are rather low. Thus, the accuracy of the PVT is expected to be high. On the contrary, ionospheric corrections are misleading when they are associated with large RMS errors. In such cases, GIM corrections should not be used for testing purposes.
2. Metric Although the RMS of the positioning error is a widely used metric to assess the accuracy of GIMs, the mode or percentiles of the distribution of 3D positioning errors provide a more robust statistical comparison. The more heterogeneous the actual positioning errors are, the less meaningful the RMS becomes because it is driven by the presence of outlier values in the tail of the error distribution function. 3. Measurement selection and processing The data used for calibrating the accuracy of GIMs through the PVT error should be as clean as possible from error sources (e.g., multipath, satellite orbits and clocks). In this case, the use of carrier phase measurements is preferable to pseudoranges. Computing the ambiguities offline with the CCL procedure and using the CCL measurements as precise pseudoranges, avoids estimating the carrier phase ambiguities simultaneously with the PVT as in standard single-frequency PPP that absorbs part of the GIM error.
In conclusion, tests based on the position domain are useful assessments that can be performed without a complex processing facility. Once the above-mentioned issues are considered, we are capable of thoroughly assess the accuracy of any GIM and obtain equivalent results to those tests using the ionospheric delay domain. In particular, executing the PVT test with CCL measurements, we showed that, in well-sounded areas, the two-layer Fast-PPP GIM provides mode errors of 0.24 m in 3D, which is several times lower than other GIMs. This result agrees with the tests based on the ionospheric domain presented in previous works.  Table 1 for 2014. The bin size is 3 cm. The top panel corresponds to the entire distribution and the bottom panel corresponds to errors smaller than 5 m. The labels correspond to the same products as those in Fig. 2