Ionospheric corrections tailored to the Galileo High Accuracy Service

The Galileo High Accuracy Service (HAS) is a new capability of the European Global Navigation Satellite System that is currently under development. The Galileo HAS will start providing satellite orbit and clock corrections (i.e. non-dispersive effects) and soon it will also correct dispersive effects such as inter-frequency biases and, in its full capability, ionospheric delay. We analyse here an ionospheric correction system based on the fast precise point positioning (Fast-PPP) and its potential application to the Galileo HAS. The aim of this contribution is to present some recent upgrades to the Fast-PPP model, with the emphasis on the model geometry and the data used. The results show the benefits of integer ambiguity resolution to obtain unambiguous carrier phase measurements as input to compute the Fast-PPP model. Seven permanent stations are used to assess the errors of the Fast-PPP ionospheric corrections, with baseline distances ranging from 100 to 1000 km from the reference receivers used to compute the Fast-PPP corrections. The 99% of the GPS and Galileo errors in well-sounded areas and in mid-latitude stations are below one total electron content unit. In addition, large errors are bounded by the error prediction of the Fast-PPP model, in the form of the variance of the estimation of the ionospheric corrections. Therefore, we conclude that Fast-PPP is able to provide ionospheric corrections with the required ionospheric accuracy, and realistic confidence bounds, for the Galileo HAS.


Introduction
The Galileo High Accuracy Service (HAS) is a new capability of the European Global Navigation Satellite System (GNSS). The Galileo HAS will transmit differential corrections to the navigation message already being broadcast by Galileo and GPS satellites, to enable precise point positioning (PPP) on a global scale (Malys and Jensen 1990;Zumberge et al. 1997). The Galileo HAS will be free of charge, following the European Commission Implementing Decision 2018/321 (EC 2018). The Galileo program has defined and is currently testing the message standards on the E6 frequency band (Fernández-Hernández et al. 2020;Borio et al. 2020). The approach selected for disseminating the HAS corrections is the state space representation (SSR), which separates dispersive and non-dispersive errors affecting the ranging signals. The Radio Technical Commission for Maritime Services (RTCM) has standardized the SSR correction format in its standard 10403.3 (RTCM 2016). The early testing phase ("phase 0") of the Galileo HAS only includes satellite orbit and clock corrections (i.e. non-dispersive effects), but its initial service ("phase 1") will include also inter-frequency biases (IFBs) and its full capability ("phase 2") will also include ionospheric delay corrections, at least in Europe. When in its full capability, the user performance requirements are 20 cm 95% horizontal accuracy, 40 cm 95% vertical accuracy, and a convergence time of 100 s (GSA 2020). This short convergence time will be achieved thanks to an accurate ionospheric model computed in real time, in such a way that any user could derive its ionospheric corrections within the service area (Rovira-Garcia et al. 2015).
Regarding the ionospheric corrections, two alternatives exist for broadcasting the total electron content (TEC) for GNSS applications. On the one hand, the Slant TEC (STEC) corrections mechanism selected, for instance, by the Centimetre Level Augmentation Service (CLAS) that is broadcast by the Quasi Zenith Satellite System (QZSS) (MEC 2015). This STEC approach is closer to Network Real Time Kinematic (NRTK) Remondi (1985) than to PPP, and it is suitable to service a region such as the Japanese archipelago. Indeed, in a regional/local context it is affordable to transmit the STECs associated with the measurements of a dense network of receivers without requiring a huge bandwidth for the message. However, this approach is not possible in a continental/global context, which would require an unaffordable bandwidth for the message, especially in a multi-constellation scenario.
On the other hand, ionospheric corrections can be broadcast using vertical TEC (VTEC) values at a set of ionospheric grid points (IGPs). This VTEC approach is currently used by the Satellite-Based Augmentation System (SBAS) (RTCA 2016), as it is more suitable for covering a wide area such as an entire continent with only few tens of stations, regardless the number of operating satellites. The Galileo HAS foresees to use a VTEC correction system, and the fast precise point positioning (Fast-PPP) model described in (Rovira-Garcia et al. 2015;Rovira-Garcia et al. 2016a) can provide such corrections.
Provided that satellite orbits and clocks can be computed in real time with accuracies at the level of few centimetres in the so-called geodetic filter (Zhang et al. 2018), one should demand similar accuracies to the ionospheric corrections. In this sense, a value of one total electron content unit (TECU), where 1 TECU corresponds to 16.24 cm in the L1/E1 GNSS frequency band (1575.42 MHz), could be taken as a reference for the required accuracy for the ionospheric corrections in order to provide a HAS (Rovira-Garcia et al. 2015,2016b. Different users can benefit of having ionosphere corrections with an accuracy not much worse than the errors in the satellite orbits and clocks. Firstly, multi-frequency receivers shorten the time required to obtain a decimetre level of accuracy (i.e. the convergence time), as it is required in the phase 2 of Galileo HAS. This is because ionospheric corrections help the navigation filter to separate (i.e. decorrelate) the carrier phase ambiguities from the other parameters estimated: coordinates, time offsets and troposphere. Secondly, as it was shown in Rovira-Garcia et al. (2020), single-frequency receivers improve the navigation accuracy to the sub-metre level, because the ionospheric contribution to the navigation error, the largest error contribution besides receiver multipath and noise, is largely mitigated.
As it was shown in Rovira-Garcia et al. (2020), the accuracy of ionospheric models is not homogenous. In regions well-sounded by multiple reference stations, the ionospheric corrections present rather low variances, whereas the accuracy of ionospheric models degrades in poorly sounded regions as those far from the continents in open sea. In such cases, the ionospheric model should not degrade the accuracy or convergence time of the classical solution based on the ionospheric-free (IF) combination, which is usually taken as a reference. For that reason, the VTEC corrections of Fast-PPP are estimated together with its expected error, similarly to the grid ionospheric vertical error indicator (GIVEI) used in SBAS (RTCA 2016). That is, at each IGP, the value of the VTEC is distributed together with an indicator of the quality of the ionospheric correction, extracted from the co-variance matrix of the estimation of the Fast-PPP model. Far away from the reference stations, the GIVEIs of the ionospheric model should be high enough that the ionospheric correction does not help nor bias the navigation solution.
The aim of this contribution is to present some recent upgrades to the Fast-PPP ionospheric model tailored to the Galileo HAS (phase 2), which include the multi-constellation capability and to confirm the performance obtained in the past with GPS data [see Rovira-Garcia et al. (2016a) and Rovira-Garcia et al. (2020)]. In those works, the advantages of the Fast-PPP ionospheric model with respect to other models were shown focusing on the signal-in-space domain and in the position domain, respectively.
The manuscript is organized as follows: the next section details the data sets used for the present study. Then, we describe the geometry of the model that computes the VTECs, in particular the two-layer ionospheric grid, addressing some of the requirements to broadcast such ionospheric corrections. The methodology section presents the strategy used to compute the model, including the combinations of measurements and different strategies to solve the carrierphase ambiguities. The results section show some examples obtained at seven permanent stations with different baseline distances from the reference receivers involved in the Fast-PPP corrections computation. While results use a global monitoring network and grid, three stations located in Europe may be representative of the future Galileo HAS ionospheric correction performance. The manuscript finalizes presenting a summary of the main findings and conclusions.

Data set
The Fast-PPP corrections are computed by a central processing facility (CPF), which retrieves data from permanent stations belonging to the networks of the International GNSS Service (IGS) (Beutler et al. 1999) and the Rede Brasileira de Monitoramento Continuo (RBMC), specifically code and carrier phase measurements at multiple frequencies and from multiple constellations. In addition to the so-called observables, Fast-PPP also requires accurate coordinates of the permanent receivers in the network and satellite orbits with enough accuracy. We have run the Fast-PPP CPF for one week starting on day 315 in 2019 and for a week starting on day 131 in 2021. In order to evaluate the ionospheric activity, Fig. 1 depicts the values of the Along Arc TEC Rate (AATR) index (see Juan et al. 2018) every 300 s for two receivers at high latitudes, YELL (Canada) and TRO1 (Norway), and one receiver at a low latitude, CHPI (Brazil). According to Juan et al. (2018), AATRs larger than one TECU/min result from high ionospheric activity. Hence, there are only some hours of moderate activity during the week in 2019, whereas during the week in 2021, there are several periods with high ionospheric activity. Figure 2 depicts around 180 permanent stations from the IGS that were collecting GPS, Galileo and GLONASS data to compute the Fast-PPP model. With this number of permanent receivers, it is possible to build a global ionospheric map (GIM) and, in some areas like Europe or Brazil, one has an enough density of permanent receivers which allows to provide accurate ionospheric corrections. In addition to these permanent stations, we have used also the data from other seven permanent receivers implemented as rover receivers (blue boxes). The data from these receivers are used just in the geodetic filter in order to fix their carrier phase ambiguities and, consequently, to have confident values that can be compared with the model predictions using the test defined in Rovira-Garcia et al. (2016a). Notice that three of these receivers are located in mid-latitude regions, including two in Europe and one in the USA, while the other three are located in Brazil, i.e. in a region where it is expected large ionospheric gradients and rates. Table 1 depicts the rover coordinates and the distances to the nearest station used to compute the Fast-PPP ionospheric corrections. As it can be seen, these distances from the nearest reference receiver are around 100, 200 or 300 km, at mid, high and low latitudes, which are distances where it is expected Fig. 2 Distribution of permanent stations from the IGS network. The red empty squares depict receivers used to compute the Fast-PPP ionospheric model, whereas the filled blue squares depict rover receivers that have not participated in the model calculus that the ionospheric model would have the enough quality for providing HAS, i.e. an error of about one TECU (see, for instance, Rovira-Garcia et al. (2016b). In particular, one of the Brazilian receivers (named CEFT) is located more than 1,000 km to the East from the nearest permanent station. Therefore, this rover receiver can be considered as an extreme case where the ionospheric corrections would not help to improve the navigation solution.
In order to contextualize the performance of the Fast-PPP ionospheric model, we have applied the same test to the rapid combined GIM of the IGS, termed IGRG. The IGRG GIM is Table 1 RMS of the post-fit residuals to Eq. (11) for seven rovers located at different distances from the nearest reference station used to compute the Fast-PPP ionospheric model. a single-layer model that assumes the ionospheric delay to be concentrated in a thin shell at 450 km of altitude. Note that the IGRG GIM is obtained in post-process mode using forward and backward process, whereas Fast-PPP is computed in real time using only the forward mode. Finally, the IGRG GIM contains TEC maps every 2 h, whereas Fast-PPP every 5 min.
In spite of the large refresh time, IGRG yields similar results to other IGS GIMs with shorter time updates and therefore is representative of the IGS GIMs.

Requirements for a HAS ionospheric model
Besides the accuracy of the ionospheric model, one of the key points of model tailored for HAS is the bandwidth of the message to be broadcast to the users, that is, the number of parameters to be computed and transmitted in real time and the required refresh time, in such a way that the optimum ionospheric model should be a trade-off between accuracy and message bandwidth. Galileo satellites will transmit an ionospheric correction message over Europe when Galileo HAS reaches its full capability (i.e. the aforementioned "phase 2"). In order to transmit the ionospheric correction message swiftly, and also due to the possibly limited monitoring capability, the IGP grid can be reduced to the European continent. According to the current HAS SIS message structure and coding scheme (Fernández-Hernández et al. 2020), such ionospheric message could fit into 26 pages to cover the main European land. Therefore, a user could receive the ionospheric values in well under 100 s, most probably around 30 to 60 s depending on the final message specifications (e.g. elevation mask, service area and quantization).
The following subsections provide preparatory discussions of key aspects of the HAS ionospheric message, exploring possible solutions. In particular, we deal with some considerations on the geometric used to compute the Fast-PPP ionospheric model and how a user must take those into account in order to compute the ionospheric corrections for any receiver-satellite ray.

Geometry of the ionospheric model
Regarding the accuracy requirement, Fast-PPP uses a duallayer model (see, for instance, Orús et al. 2021) to correct the ionospheric delay. The first layer contains 7176 IGPs at a height of 270 km and accounts for the main ionospheric content, whereas the second layer contains 1792 IGPs, i.e. about one-fourth of the IGPs in the first layer, at a height of 1600 km and accounts for the plasmaspheric content. In this manner, Fast-PPP overcomes the main simplification of the thin-shell models at a unique height, allowing some degree of vertical 3D representation of the TEC distribution. Notice that, for any ground observation, these very different heights for the ionospheric layers will have a clear different obliquity factor on each layer. This will allow to distinguish an ionospheric delay occurring at the top layer from one occurring at the bottom layer. This capability represents a reduction of the mismodelling with respect to a single-layer assumption, which is a key factor to meet the accuracy requirement of one TECU for HAS (Rovira-Garcia et al. 2015,2016b. Because the main drivers of the ionospheric state are the geomagnetic field and the Sun, the Fast-PPP ionospheric model defines the IGPs based on local time (LT). This allows enlarging the refresh time of the model. Moreover, in order to reduce the number of parameters to be transmitted, the IGPs are defined on a grid on LT and the MOdified DIP (MODIP) latitude (Rawer 1963). The VTEC grid is equally spaced over each MODIP band, which implies an irregular distribution of the IGPs in terms of longitude and latitude. Specifically, the IGPs of the bottom layer are equally spaced every 2.5 degrees in MODIP, whereas the IGPs of the top layer are equally spaced every 5 degrees of MODIP.
Regarding the LT spacing between IGPs, it decreases approximately with the cosine of MODIP: LT layer MODIP layer / cos(MODIP) (1) Figure 3 illustrates the Fast-PPP grid at each layer. The plots in the left column use the LT-MODIP coordinates, whereas the plots in the right column depict the irregular distribution of its geographic coordinates (longitude and latitude at 0 h of Universal Time). As it can be seen, this choice provides a high density of IGPs at low latitudes (i.e. low MODIP). Provided that enough coverage of permanent stations is available at such regions, this IGP distribution improves the accuracy of the Fast-PPP ionospheric model in regions around the geomagnetic equator, where great temporal and spatial gradients are experienced at the ionosphere. The irregular grid presents an additional benefit, which is the reduction in the number of IGPs to be estimated in comparison with a regular grid. As a numerical example considering the first layer with a regular resolution of 2.5°× 2.5°in latitude and longitude, would require 10,500 IGPs instead of Fig. 4 Bi-linear interpolation of the VTEC at four IGPs surrounding the IPP the 7179 IGPs of the irregular grid using Eq. (1). This corresponds to a reduction of 30% not only in the computational load in the CPF, but also in the message bandwidth.

Usage of ionospheric corrections
The proposed grid based on an irregular distribution of the IGPs requires a different interpolation scheme to the one used in SBAS (RTCA 2016) or in the IONosphere map EXchange format (IONEX) standard (Schaer et al. 1998), when Fast-PPP users need to compute the ionospheric corrections. The procedure is illustrated in Fig. 4, where the ray between the user receiver and the satellite intersects each ionospheric layer at the so-called ionospheric pierce point (IPP). Then, assuming that the IPP is surrounded by four IGPs located in the same layer at a given time, the user interpolates the VTEC at the four IGPs using the linear distances x and y: where V a and V b are computed as By simple algebraic manipulation, it follows that the algorithm is equivalent to the linear relationship: From Eqs.
(2)-(4), it follows that in the case of a regular grid (i.e. x 1 x 2 ) the algorithm corresponds to the interpolation method adopted in SBAS.

Methodology
During the last years, a trend has been consolidated in the GNSS parameter estimation consisting on the processing of the raw measurements in an undifferenced and uncombined manner (see, for instance, Odijk et al. 2016). This general approach includes also integer ambiguity resolution (IAR) and even the estimation of the parameters of an ionospheric model. Thanks to these estimation strategies, the different parameters can be determined with a high accuracy. However, in a worldwide context, this way of processing involves tens of thousands of parameters that should be estimated in real time in order to provide a HAS. For instance, for a global network with 200 receivers each one tracking the measurements of, let us assume, 30 satellites at three frequencies, this would require to estimate around 2·10 4 carrier phase ambiguities together with close to 10 4 parameters of the ionospheric model.
One of the main advantages of the data used to compute the Fast-PPP ionospheric model is that it does not have to directly deal with (i.e. estimate) the ambiguities of the carrier phase measurements at the same time as the ionospheric delays. In order to do so, a previous independent module called geodetic filter performs IAR, handling nondispersive combinations of GNSS signals. Indeed, firstly, the ionosphere-free (IF) combination of carrier phases is used to determine the ambiguity in such combination, B IF , using PPP models, the precise knowledge of station coordinates, and satellite orbits and clocks. The B IF is estimated as a real number in what is known as a floated ambiguity. Notice that thanks to the knowledge of the receiver coordinates, the confidence (formal error) of the B IF ambiguity, i.e. the square root of the diagonal elements in the covariance matrix of the estimates, is at the level of 1 cm since the beginning of the arc of data.
In parallel, the carrier phase ambiguity of the wide lane (WL) combination, B WL , is calculated through the Hatch-Melbourne-Wübenna (HMW) combination (Hatch 1982), that is, the difference of the wide lane (WL) combination of carrier phases and the narrow lane (NL) combination of code pseudoranges. For a satellite j and a receiver i, B W L (expressed in metres) can be written as: Being λ WL the wavelength of the WL combination λ WL c f m − f n , f m and f n are two different frequencies selected by the Fast-PPP operator, c is the speed of light, N WL the integer part of the carrier phase ambiguity (expressed in cycles) and δ WL the instrumental delay of the satellite or the receiver at the WL combination.
After few minutes of taking data, the double difference (DD) of N WL between pairs of satellites and receivers can be fixed to an integer number. This can be done thanks to the wavelength of the WL combination (e.g. 86 cm using GPS frequencies L1 and L2 and 75 cm using Galileo frequencies E1 and E5a), which is large enough to allow DD IAR, despite the noise of the B WL introduced by the NL combination of pseudoranges. Once DD IAR is performed for B WL , a link appears between the four WL ambiguities involved in the DDs, which reduces dramatically the formal error of each individual estimate of B WL .
The final step of the geodetic filter consists in obtaining the ambiguity in the geometry-free (GF) combination, B GF , for a receiver i and a satellite j, following (Sanz et al. 2013): Figure 5 illustrates the IAR process for three stations and one GPS satellite. The top plot depicts the values of the confidence level of the estimates, whereas the bottom plot illustrates the estimated B GF values. Once the confidence level is reduced, carrier phase ambiguities are fixed, which for this example occurs around 11 h. As a consequence of the IAR, the different B GF values for different receivers or satellites remain linked in the estimation process. That is, they evolve in time in the same manner. This occurs because a constraint between them has been added to the Kalman filter, which increases the robustness of the overall ambiguity estimation.
The IAR performed before the Fast-PPP ionospheric model estimation represents a clear advantage with respect to other state-of-the-art approaches that solve the B GF arc by arc in an independent (i.e. disconnected) manner. This is the case of the carrier phase to code levelling (CCL) process (Mannucci et al. 1998), which estimates B GF averaging the difference between the GF combination of code pseudoranges P GF and carrier phase measurements L GF for each continuous arc of the samples: Figure 6 illustrates the advantage of estimating B GF with the geodetic filter over the CCL. For that purpose, we depict, for all arcs of data within day 315 of 2019, the difference of the two instantaneous estimations of B GF with respect to their final estimate, B GF (t final ), which is assumed the best estimate of the ambiguity because it involves all the data in the arc. Therefore, the figure compares the drift of the estimates ignoring possible final biases. As expected, during the first epochs in the instantaneous estimations, the errors of both methods are large. However, in almost all the arcs of the geodetic filter approach, the differences between the instantaneous and final estimates is below one TECU since the beginning. In contrast, the B GF estimates with the CCL approach depict variations of several TECU, i.e. larger than the required accuracy for a HAS. Because the input data of the ionospheric model are the unambiguous STEC (STEC L GF − B GF ), the more accurate the B GF are, the more accurate the ionospheric estimates will be.
Besides the improvement of the unambiguous STEC precision, performing IAR at the CPF facilitates the IAR process at the user side. Indeed, once IAR is done inside the geodetic filter, it is easy to compute receiver and satellite phase biases (Rovira-Garcia et al. 2021). These satellite biases can be broadcast every few minutes by a HAS in order to enable IAR on the user side. Figure 7 depicts, for two GPS and two Galileo satellites, the real-time estimates of the phase biases of the WL combination (left) and the L1 (right). As it can be seen, the noise in the phase bias estimates is usually well below a tenth of one cycle, which guarantees the IAR. Only the phase bias of L1 can present some problems in the epochs when the satellite is poorly tracked (e.g. orbiting over the oceans). In any case, the noise remains well-below half a cycle.

Testing the ionospheric model
The assessment of the accuracy of any ionospheric model is a key, and difficult, point that has been addressed using several methods. The difficulty of such assessments relies on the procedure to obtain confident ionospheric values with  enough accuracy to be used as a reference. Indeed, direct ionospheric measurements, using different techniques, are affected by unknown biases that shall be removed before using such measurements as the reference. For instance, Orús et al. (2021) used simulated data in order to discard the problem with the biases, focusing the assessment on the different geometry of the ionospheric model. However, unlike simulated data, actual ionosphere has irregularities at different spatial and temporal scales that can affect the estimation of the ionospheric parameters such as hardware biases.
In order to assess ionospheric models with actual data, Orús et al. (2005) introduced the self-consistency test (SCT), which nowadays is widely used. SCT is based on the variation of L GF with respect to a L GF measurement done at a reference time (t ref ), that is, L GF (t) − L GF (t ref ). Note that such differentiation cancels the carrier phase ambiguity B GF present in the carrier phase measurements, because all the measurements belong to the same continuous arc of data that share the same ambiguity value. Then, the SCT consists on comparing the variation of L GF , in metres, with the variation of the ionospheric model predictions for these measurements, , in TECU, defining the variation of the STEC prediction error as: where i corresponds to a receiver and j to a satellite, the factor α m 40.3 · 10 16 / f 2 m converts an ionospheric delay from TECU to metres of delay at the frequency f m .
Once the carrier phase ambiguity B GF is removed, the variation of L GF along the arc has an error at the level of 1 cm (at the level of a tenth of one TECU); this error level should be enough for assessing ionospheric models tailored for HAS.
In this way, the metric to quantify the SCT results could be the RMS of STEC j i or any other statistic like percentiles. There are different choices to select t ref for performing the comparison in Eq. (8). In Orús et al. (2005), the SCT was defined using the comparison between two epochs with the same satellite elevation, whereas Hernández-Pajares et al.
(2017) selected t ref as the epoch with the maximum satellite elevation in the dSTEC tests, which is based on the same idea. It is easy to see that both the SCT and the dSTEC tests, based on the STEC variation along the arc, are equivalent to estimate, per each continuous arc, a carrier phase ambiguity (B GF ) as: where the metric for the test is the residual of this fitting, e.g. the RMS of the post-fit residual. Note that for a worldwide network of about one hundred stations and one hundred of satellites (in a multi-constellation scenario), Eq. (9) involves estimating tens of thousands ambiguities (whose exact number depends on how many arcs occur per receiver-satellite pair). This huge number of ambiguities to be estimated implies that the assessment is done just in a local/regional scale (what the IPP sweeps during the length of the arc of data), in such a way that any regional bias of the ionospheric model could be absorbed by the B GF estimates. Therefore, because the bias values are arbitrary (or the selection of t ref ), only the standard deviation of the fitting of Eq. (9), arc by arc, makes sense for the assessment.
In contrast, we can take advantage from IAR, which allows us to determine the integer values of the carrier phase ambiguities at each frequency, N j mi and N j ni . Then, Eq. (9) becomes: where the IAR values of the ambiguities have been subtracted in the left-hand side. The δ ki and δ j k are instrumental delays (i.e. the phase biases) for the receiver i and the satellite j at each frequency k. Note that the right-hand side of Eq. (8) are constant values, which can be estimated for instance once per day, as it was proposed in Rovira-Garcia et al. (2016a): where the constant values k include the instrumental delays of each receiver and each satellite. Note that in this last test, the fitting involves only the sum of the number of receivers and the number of the satellites considered (some hundreds). In this way, the post-fit residual of the global fit over the whole network becomes a confident measure of the self-consistency of the ionospheric delay estimates. That is, for the same previous case of a network of one hundred receivers and one hundred satellites, the test using Eq. (11) would estimate less than two hundred parameters, two orders of magnitude less than in the SCT or dSTEC test involving Eq. (9). This is the reason why the test proposed in Rovira-Garcia et al.
(2016a) is one of the assessments routinely used by the Fast-PPP ionospheric model. Figure 8 illustrates, for two rover receivers, the residuals of the estimation of the daily constant k values in the right-hand side of Eq. (11), after subtracting two different ionospheric models I to the unambiguous GF combination of carrier phase measurements. The black dots depict the performance of the Fast-PPP ionospheric model obtained in real time, whereas the red squares depict the results obtained by the IGRG GIM. The results, extended in Table 1 for the seven rovers and for every week in each period, confirm those presented in Rovira-Garcia et al. (2016a). That is, in wellsounded areas, the root mean square (RMS) of the residuals of Fast-PPP is two or three times lower than the RMS of the IGS GIM. As expected, there is less improvement when the station is located in poorly sounded areas.
In Table 1, we have split the performance into Galileo and GPS measurements, in order to see differences between the two constellations (notice that IGRG is computed using only GPS data). However, the resulting residuals are quite similar for the two constellations. Note that stations REDU (Redu, Belgium), MARS (Marseille, France) or TRO1 (Tromso, Norway) may be representative of future ionospheric corrections of Galileo HAS over central Europe, assuming a monitoring network with a similar density over Europe and surroundings as that used in our test.
As it can be seen in Table 1, except for the farthest receiver CEFT (at more than 1000 km from the nearest reference receiver), the RMS of the error of the ionospheric corrections using the Fast-PPP ionospheric model is clearly smaller than one TECU and is several times smaller than those using the IGRG ionospheric model. This confirms the results presented, for instance, in Rovira-Garcia et al. (2020) using just GPS data. In the case of CEFT, these large ionospheric errors would not help to improve the navigation solution. However, these errors are still clearly smaller than those using the IGRG GIM.
We now turn our attention to the assessment of the actual errors (Fig. 9) and the predicted errors (Fig. 10) for both the IGRG GIM (top rows) and the Fast-PPP model (bottom rows) for the two analysed periods in 2019 (left columns) and 2021 (right columns). The chosen metric that is depicted in every plot is the complementary cumulative distribution function (CCDF), to allow assessing the percentiles of actual and predicted errors of both ionospheric models. In every plot of Fig. 9, the 5% value in the CCDF is indicated with a dashed line in order to identify the 95th percentile of each error distribution.
As it can be seen in the bottom plots of Fig. 9, 99% of the Fast-PPP ionospheric corrections have actual errors below one TECU in well-sounded areas and in mid-latitude stations. We can observe that the 95% error percentiles, which can be identified thanks to the dashed line, do not degrade between the right and left column plots, despite the fact that some rovers experience an ionospheric activity substantially higher in the 2021 period. Therefore, one can conclude that Fast-PPP is able to provide the required accuracy for HAS in such conditions. The performance is degraded for stations at low latitude or at the highest latitude, where in these regions the Fast-PPP model uses a higher process noise for computing the VTECs, in order to account for the larger ionospheric activity expected in this area. However, for these receivers and during the two periods, the 95% error remains under 2 TECU. In contrast, the top row plots depict the errors for the IGRG GIM, where, in the best conditions (REDU station in mid-latitude), more than the 10% of the ionospheric corrections have an error larger than one TECU.
Looking the RMS of Table 1 or the percentiles inferred from Fig. 9, it is worth noting that the ionospheric errors in the present work are significantly lower than those in Rovira-Garcia et al. (2020). This is because the results on that paper were obtained during year 2014, close to Solar Cycle maximum. Consequently, using the ionospheric corrections of this work for navigation purposes would produce even better positioning results than those presented in Rovira-Garcia et al. (2020). Figure 10 depicts the predicted errors of the ionospheric corrections, which are derived from the covariance matrix of the estimates. It is interesting to see that in the case of the Fast-PPP model, these formal errors provide some kind of guarantee to the quality of the ionospheric corrections, i.e. larger predicted errors are linked to larger actual errors in the ionospheric corrections depicted in Fig. 9. This is not the case for the IGRG GIM, where the larger actual errors in the

Conclusions
Ionospheric models with an accuracy at the level of one TECU can be useful for high-accuracy applications, such as the Galileo HAS. The present work has described some novel characteristics of the Fast-PPP ionospheric model that contributes to provide such accuracy. The first one is related to the geometry of the model, which, as it was shown in previous works, uses a dual-layer description of the ionospheric and plasmaspheric delays and the use of the MODIP latitude and LT to distribute the IGPs. The regular distribution based on MODIP increases the resolution at equatorial latitudes. Unlike other ionospheric models, the IGP distribution in LT is not regular. This choice reduces 30% of the IGPs needed to be computed and broadcast. Based on this model, some aspects concerning the requirements and usage of ionospheric corrections for Galileo HAS have been discussed. A second novel characteristic is related to the quality of the measurement data used as input to compute the Fast-PPP model, which uses IAR instead of CCL to obtain unambiguous carrier phase measurements feeding the ionospheric model.
The accuracy of the Fast-PPP ionospheric model has been examined using rover receivers that have not participated into the computation of the model, confirming accuracies at the level of one TECU in well-sounded regions. In addition, the Fast-PPP model provides a realistic prediction of the error through the formal errors, so that HAS users can decide whether to apply such ionospheric corrections. CPF used in the study. GGC performed the data gathering and processing. ARG and CCT analysed the results and wrote the manuscript. IFH, RO and DB contributed to the technical discussion and organization of the manuscript. All authors provided advice and critically reviewed the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. The present work was supported in part by the European Space Agency contract IONO4HAS 4000128823/19/NL/AS, by the project RTI2018-094295-B-I00 funded by the MCIN/AEI 10.13039/501100011033 which is co-founded by the FEDER programme and by the Horizon 2020 Marie Skłodowska-Curie Individual Global Fellowship 797461 NAVSCIN.

Data availability
The authors acknowledge the use of data and products provided by the International GNSS Service and the Rede Brasileira de Monitoramento Continuo, as detailed in the section Data set. All the data and products are publicly available through the respective organizations' websites.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.