A multi-frequency method to improve the long-term estimation of GNSS clock corrections and phase biases

The space segment of the Global Navigation Satellite System (GNSS) is equipped with highly stable atomic clocks. In order to use these clocks as references, their time offsets must be estimated from ground measurements as accurately as possible. This work presents a multi-frequency and multi-constellation method for estimating satellite and receiver clock corrections, starting from unambiguous, uncombined, and undifferenced carrier-phase measurements. A byproduct of the estimation process is phase biases (i.e., the hardware delays of the carrier-phase measurements occurring at receivers and satellites). The stability and predictability of our clock estimates for receivers and satellites (GPS and Galileo) are compared with those obtained by the International GNSS Service (IGS), whereas the phase biases are assessed against two independent determinations involving combinations of carrier-phase measurements. We conclude that the method reduces day boundary discontinuities in the clock corrections, and that the estimated phase biases reproduce variabilities already observed by other authors.


INTRODUCTION
Time and frequency transfer is one of the main uses of the Global Navigation Satellite System (GNSS), besides geodetic or navigation applications. Stable ground clocks equipped with a passive hydrogen maser (PHM) or rubidium atomic frequency standard (RAFS) can be used as global time references. In such a way, the clocks associated with other ground receivers or satellites can be estimated and referred to any of these reference PHMs.
usually done in a two-step process (Bock et al., 2009). In the first step, all parameters are computed with a low sampling rate (e.g., five minutes). Then, in a second step starting from the solutions for the low-varying parameters, satellite and receiver clock corrections are estimated with a higher sampling rate (e.g., 30 seconds, even 1 second), as seen in Bock et al. (2009) and Rovira-Garcia et al. (2016).
In such an estimation procedure, GNSS measurements are collected by a worldwide network of receivers located at permanent stations with coordinates accurate to the centimeter. This allows subtraction from the measurements' well-known geodetic models, such as those used in the precise point positioning (PPP) technique (Malys & Jensen, 1990;Zumberge et al., 1997). In this manner, PPP accurately corrects antenna offsets, solid and ocean tide displacements, relativistic and gravitational effects, tropospheric propagation delays, and windup effect among other terms.
However, the ionospheric delay experienced by measurements cannot be modeled with enough accuracy. To eliminate first-order ionospheric delay, the standard procedure for estimating satellite and receiver clock corrections uses the dual-frequency ionospheric free (IF) combination of code pseudoranges and carrier-phase measurements (Subirana et al., 2013). The obtained residuals Δ and Δ (i.e., measurements corrected from the PPP models) for any receiversatellite pair can be used to estimate the satellite and receiver clock corrections as white noise processes and the receiver hardware biases between different constellations as a random walk process. As commented before, this usually occurs in the second step at a high throughput, starting from the orbits and clock corrections previously estimated in the first step.
Alternative to the use of IF combinations for estimating clock corrections, some authors like Schöenemann et al. (2011), Odijk et al. (2016), and Strasser et al. (2019) have obtained similar results using uncombined measurements, the so-called raw observation processing, in which ionospheric delays are estimated as white noise processes or through the estimation of additional parameters of an ionospheric model.
Since the early days of GPS, the importance of integer ambiguity resolution (IAR) has been recognized to enhance the accuracy of geodesy applications including clock determination (Blewitt, 1993). IAR reduces the number of parameters to be estimated and avoids the correlations between ambiguities and clocks. Thus, IAR is a standard procedure performed by different analysis centers (ACs) contributing to the IGS. In particular, the Double-Differenced (DD) method (Beutler et al., 1986) is applied to compute daily batches of precise orbits and clock corrections. Moreover, different authors showed the advantages of performing undifferenced IAR in PPP (Ge et al., 2008;Laurichese & Mercier, 2007).
However, as the processing is done in daily batches, IAR does not guarantee that the same integer values for ambiguities can be obtained every day. That is, in the daily estimates, arcs of continuous carrier-phase measurements spanning two consecutive days can present different values each day on their ambiguities and, consequently, in the clock estimates.
This problem is known as day boundary discontinuity (DBD), such as seen for instance in Ray and Senior (2003). If continuity is not ensured after performing the daily IAR process, DBDs in clock estimates can reach the decimeter level (i.e., several tenths of a ns).
The reason behind DBDs is that the noise and the instrumental delays present in the code pseudoranges can vary differently from those in the carrier-phase measurements. A solution to the DBDs proposed in Collins et al. (2010) was to define a different clock for pseudorange and carrierphases (i.e., the decoupled clocks technique). DBDs on satellite clock solutions introduce artificial discontinuities that limit the accuracy of clock analysis (e.g., frequency stability) for periods longer than one day (Collins et al., 2010;Dach et al., 2006;Petit et al., 2015).
Recently, Rovira-Garcia et al. (2021) presented a methodology using only unambiguous carrier-phase measurements [as in Collins et al. (2010)]. The method guarantees the continuity of the unambiguous carrier-phase measurements over days, even when IAR is performed in daily batches and thus is able to reduce the DBDs between the clock estimates of such adjacent daily batches.
The aim of the present work is, on one hand, to extend the study of Rovira-Garcia et al. (2021) to hundreds of days and, on the other hand, to show that the methodology can be also applied when more than two frequencies are involved. Notice that handling more than two frequencies implies addressing different temporal evolution of the phase biases at different frequencies, as pointed out in Montenbruck et al. (2012) for the GPS satellites of Block IIF.

DATA SET
We have applied our methodology to more than 340 days in the year 2017, extending the proof of concept of 16 days of data presented in Juan et al. (2020) and Rovira-Garcia et al. (2021). In order to estimate the satellite and receiver clock corrections, together with the phase biases, we downloaded RINEX observation files containing code pseudoranges and carrier-phase measurements every 30 s. The frequencies used by this method are L1, L2, and L5 for GPS, and E1 (L1), E5a (L5), and E5b (L7) for Galileo. This choice applies to each of the approximately 150 receivers located at permanent stations within the IGS network that were used in the study. The redundant number of receivers allows robust estimation of the phase biases, despite the fact that some regions (e.g., oceans) are not as well covered as continental areas.
In order to compare our estimates, we used RINEX clock files with precise clock solutions computed by different IGS ACs contributing to the Multi-GNSS Experiment (MGEX) (Montenbruck et al., 2017). Specifically, the selected clock products for the study are: the GRM product of the Centre National d'Études Spatiales (CNES), the GBM product of the GeoForschungsZentrum (GFZ), the COM product from the Center for Orbit Determination in Europe (CODE), and the repro4 product from the European Space Operations Center (ESOC). The first three products are publicly available at the IGS MGEX product archive, whereas repro4 is an adhoc process of ESOC in the context of the projects named the Gravitational Redshift Experiment with Eccentric Satellites (GREAT) funded by ESA (Delva et al., 2018;Herrmann et al., 2018;Rovira-Garcia et al., 2021). All ACs provided satellite clock corrections at an interval of 30 s, although COM initially contained estimates every 300 s. Eventually, after day 226 in 2017, COM provided determinations every 30 s as well.
We focused our study on GPS and Galileo clock estimates, which were available from all ACs except ESOC, which only provided estimates for Galileo satellites in its repro4. In this ESOC solution, the reference receiver for every day during the entire 2017 was the station BRUX, equipped with a very stable PHM atomic clock.
In contrast to ESOC, the other ACs typically re-aligned the clock estimates to the broadcast GPS time for the satellites as an underlying timescale (Ray et al., 2017). That is, first the ACs estimated all clock corrections by setting the value of one reference station (that could vary from day to day) to zero, and then, in a second step, the ACs subtracted from the clock estimates the mean value of all satellite clock determinations.
This approach is robust to variations of the reference station from one day to another, thus should produce continuous satellite clocks. However, when the satellite set changes, the mean value of the satellite clocks changes accordingly, in a similar way as the jumps in the differential code biases (DCBs) reported in Sanz et al. (2017). This issue, among others, can produce DBDs as shown in the results section.
The present study focuses on the aforementioned second step, in which satellite and receiver clock corrections are estimated starting from a previous determination of slow varying parameters. In this regard, our starting product is the repro3 solution; the penultimate version of the prod-ucts computed by ESOC for the GREAT study. The main differences between repro3 and repro4 are that repro3 includes satellite orbits and clock corrections for GPS and Galileo and the reference station changes from day to day. The clock corrections examined in the results section correspond to the repro4 product.
Our study sought to evaluate the nominal performance of the clock corrections determined by our approach and those computed by the aforementioned four IGS ACs. Then, in order to perform a fair comparison, we scanned the products computed during the entirety of 2017. We looked for periods in which none of the five products presented a large discontinuity associated, for instance, to the chosen clock reference or any other potential problem associated with the clock determination methodology.
After the examination, we selected two common periods: days 191 through 230 and days 281 through 310. In these two periods, none of the five solutions presents an anomalous behavior and thus are representative of a nominal performance.

METHODOLOGY
We propose a method to estimate clock corrections and its associated inter-frequency biases using only carrier-phase measurements. As in Strasser et al. (2019), our approach uses uncombined carrier-phase measurements (instead of the standard IF combination of carrier phases) in which the integer part of each carrier-phase ambiguity has been fixed through a previous estimate of the phase biases. Moreover, our approach supports more than two frequencies in processing multiple constellations and guarantees the continuity of phase biases over days. In our methodology, the carrier-phase and pseudorange code measurements at a given frequency m between a receiver i and a satellite j can be written as: where is the geometric distance between the receiver and satellite antennae, c is the speed of light, and are the receiver and satellite clock offsets, the factor = 40.3 ⋅ 10 16 ∕ 2 converts from the total electron content unit (TECU) to meters of delay at the frequency , and account for the slant tropospheric and ionospheric delays between the receiver and the satellite. The and are the DCB for the receiver and the satellite, whereas and are the phase instrumental delays, and is the integer carrier-phase ambiguities expressed in cycles and converted to meters by the corresponding wavelength . Finally, and represent noise terms.
The key point of the methodology is to resolve the carrier-phase ambiguities in an undifferenced way. In order to guarantee the continuity of the unambiguous carrier-phase measurements over days, we have to do a rough estimation for the phase biases of both satellites and receivers at all frequencies and constellations. This previous process is described in detail in Rovira-Garcia et al. (2021) and Juan et al. (2020), allowing the estimation of the integer part of the carrier-phase ambiguities using standard PPP models.
Once the integer part of the ambiguities is removed to the raw undifferenced carrier-phase measurements in Equation (1), we then apply again geodetic models in order to obtain the phase residuals Δ (i.e., undifferenced and ambiguity-resolved carrier-phase measurements corrected from effects modeled with PPP). Δ at a given frequency m between a receiver i and a satellite j can be written as: where the instrumental delays + have been redefined in Equation (3) as: Note that Equation (3) is only applied for carrier-phase measurements for which the ambiguities are fixed, otherwise, those measurements are discarded. This corresponds to an IAR success rate around a 90%, similar to Rovira-Garcia et al. (2021) in its 16-day study. The modeling error and other errors such as thermal noise and multipath of Δ in Equation (3) are typically below one centimeter.
The method to solve Equation (3) uses the Kalman filter to estimate the clock corrections of satellites and receivers, the slant ionospheric delay between each receiver-satellite pair, and the phase biases of satellites and receivers. The estimates are produced at regular intervals of 30 s. Note that the value of the phase biases is strongly correlated with the ionospheric term in the right-hand side of Equation (3).
The process noise used in the filter to estimate the clock corrections and ionospheric delays is white noise, whereas the phase biases are estimated as a random walk.
The reason not to estimate them later as a constant offset is to allow studying the evolution of the phases biases with time that can be related to variations in the satellite temperature or in the orbital position (Montenbruck et al., 2012).
In order to solve the rank deficiency in Equation (3), we have to define one reference frequency for each device (satellite or receiver) and one reference clock for referring the rest of the clock corrections. We have selected the receiver CEBR as the reference clock, as this station is equipped with a very stable PHM atomic clock.
For the satellites, the reference frequency can be imposed by assuming that the phase biases are the same for two of the frequencies involved in the processing. That is, for each GPS or Galileo satellite, we impose the constraint: From Equation (3), this is equivalent to take the IF combinations built with frequencies 1 and 2 for each GPS satellite, and 1 and 5 for each Galileo satellite, as the satellite clock reference, respectively.
In the same way, for each receiver, we have to define one reference frequency. The reference frequency can be defined linked to the IF combination of GPS frequencies 1 and 2 , by assuming that the hardware biases are equal in the constraint: which is equivalent to consider the IF combination of these two carrier-phases as the reference for the receiver clock corrections. In this sense, when considering a triplefrequency processing per constellation, the method has to estimate for each receiver two different phase biases for GPS ( 1 , and 5 , ) and three phase biases for Galileo ( 1 , and 5 , and 7 , ) because, for the latter, we have not imposed any constraint relation between them.
Extending the same reasoning to the satellite phase biases, the method estimates two different phase biases for each GPS satellite ( 1 , and 5 , ) and two different phase biases for each Galileo satellite ( 1 , and 7 , ).

RESULTS
This section analyzes the estimates performed by the proposed method. First, we turn our attention to the input data used by the estimation. Second, we analyze the satellite and receiver clock results, specifically their predictability and stability. Finally, we analyze the satellite phase biases observing its temporal evolution. Figure 1 depicts an example of the carrier-phase residuals computed with Equation (3) for up to three frequencies and different Galileo and GPS satellites using their pseudorandom numbers (PRNs). In order to magnify the continuity of the residuals at the day change, we have applied a linear detrending. That is, for each receiver-satellite pair, we have subtracted from the residuals the value of a linear model adjusted by the whole data.

Carrier-phase residuals
In this sense, we eliminate large trends (bias and drift) in the receiver and satellite clock corrections and , but the evolution of all other terms in Equation (3) can still be observed in Figure 1, that is, mainly the clock corrections and ionosphere. The magnitude of the ionosphere variation is at the level of a few tens of decimeters and depends on the frequency as indicated at the top (f 1 ), middle (f 5 ), and bottom (f 2 ) and (f 7 ) panels.
Apart from the different stability of the depicted satellite clock corrections, the detrended residuals confirm that the daily IAR pre-process does not introduce any discontinuity at the midnight epoch. Because these residuals are the input data in Equation (3) for estimating the receiver and satellite clock corrections, and the remaining terms (ionospheric delay and phase biases) should exhibit a continuous behavior, the elimination of DBDs in the input data guarantees the inter-day continuity of the clock solutions.
Notice that in Figure 1 it is also possible to see the different stability of different satellite clocks. Indeed, the noisier clock occurs in the GPS satellite G22 (Block IIR-B), while the clocks of the Block IIF GPS satellites are more stable. It is even possible to distinguish between Galileo satellites which, according to the Notice Advisory to Galileo Users (NAGUS) published by the European GNSS Service Centre (EUSPA, 2017), during these days, E09 was operating with a PHM, whereas E22 had activated its RAFS.

Satellite clocks estimates
For two Galileo satellites operating with PHMs, Figure 2 illustrates the clock estimation performed by the proposed approach together with the products from the IGS ACs over days 191 through 230 and days 281 through 310. As an example, we selected two Galileo satellites operating with PHMs, but other Galileo satellites have similar performances. As in Figure 1, in order to better depict the continuity of the clock estimates over the day changes, we have detrended the clock solutions by subtracting a linear model computed for each AC and for each period of more than one week. Moreover, in order to enhance the plot visibility, the time series of each AC are shifted with one arbitrary value. However, errors in the radial component are not the only source of error affecting satellite clock estimates. Indeed, as it is seen in previous references, other effects related to satellite position, such as geomagnetic field, satellite temperature, or even relativistic correction, can affect the satellite clock solution. In any case, these types of errors related with satellite position are expected to be continuous over time, so they are not the main cause of the DBDs.
Regarding DBDs, it can be seen that the clock corrections estimated by the proposed approach (labeled as gAGE), which relies on a careful IAR pre-process, are continuous regardless of the day change. The remaining discontinuities observed in our clock corrections (e.g., in day 284 for E01 or day 192 for E30) are caused by the starting products (precise orbit and clock corrections) used to model the carrier-phase measurements in Equation (3). In contrast, we can observe DBDs in the clock estimates of other ACs, which reach the decimeter level (i.e., the ns level).
A standard metric to evaluate the stability of clock corrections is through the Allan deviation (ADEV) as defined in Allan (1966). Figure 3 illustrates the ADEVs for time intervals up to 533,000 seconds, which correspond to 6.17 days, for the two previous Galileo satellites E01 and E30. Note that the solution computed by CODE is excluded from the first period, as the clock estimates were produced at a sampling interval of 300 seconds and therefore does not reflect a fair comparison with clock solutions every 30 seconds.
Moreover, we have excluded days 293, 301, and 310 of the CNES solution due to an anomalous behavior affecting all the satellites simultaneously. We can see that our approach produces clock estimates that are similarly as stable as those computed by other IGS ACs. In this sense, the ADEV does not evidence the advantage of the DBD removal observed in Figure 2. This is because the DBD magnitude is approximately a few decimeters and only occurs once per day (every midnight in GPS time), whereas the ADEV computation includes many other points that belong to the same day (i.e., unaffected by DBDs). In some solutions, the fluctuations of the clock estimates are comparable to the DBDs. These two facts dilute the effect of the presence of DBDs in the ADEVs depicted in Figure 3.
In order to emphasize the advantage of the proposed method in reducing DBDs, we used one hour of clock solutions during the last hour of one day to fit a linear model using least squares. Then, this linear model was propagated one hour forward into the following day. Finally, we computed the differences between the linear model predictions and the actual clock estimates over the first hour of the following day. This procedure assessed the extrapolation errors, which is of primary interest for navigation applications, specifically, for high-accuracy service providers questing for an optimal refresh rate of satellite clock corrections. Figure 4 illustrates the extrapolation errors, aggregating all GPS and all Galileo satellites during each period considered in 2017. The histogram of prediction errors obtained from the clock estimates determined by the proposed approach reached higher peaks centered around zero than those obtained by other ACs, evidencing the benefits of the DBD removal.
The peaks of our approach for Galileo are higher than those of GPS, confirming the better predictability (i.e., stability) of Galileo clocks. The satellite clock estimates computed with three frequencies are similar to those previously estimated with two frequencies in Rovira-Garcia et al. (2021). The necessary estimation of the additional phase bias mitigated the redundancy of having an additional frequency.

F I G U R E 4
Histogram depicting the prediction errors after linear extrapolation for all GPS and Galileo satellites in the top and bottom panels, respectively; the left column depicts the first period considered in 2017, whereas the right column depicts the second period; the bin size of the histograms was 5 mm [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] The effect of the lower stability of GPS clocks can also be seen in the histograms of GPS clock estimates of other ACs: they look more similar between different ACs than those of Galileo, suggesting that the lower predictability of GPS clocks dominates the differences in estimation methodologies.
We can see that ESOC and GFZ improve the Galileo determinations from the first to the second period, but the reason is unclear. In contrast, our estimates do not substantially vary from the first to the second period maintaining a smaller distribution of prediction errors in comparison to other ACs, whose distributions are not centered in zero because the DBDs are not necessarily symmetrical.

4.3
Receiver clock estimates Figure 5 depicts the results obtained for receiver clock corrections during eight days in each of the aforementioned periods. In order to better illustrate the continuity of the clock estimates over the midnight epoch, we selected two IGS stations, BRUX and MGUE. These two stations are equipped with stable PHM clocks and hence a high predictability can be expected. Moreover, all the ACs provide clock solutions for these stations during the selected periods. As in the case of the satellites in Figure 2, we have detrended each time series of the clock estimates with a linear model. We can observe that the receiver clock corrections determined by the proposed methodology present a high continuity, because of the mitigation of the DBDs in input carrier-phase measurements.
The most stable solution is obtained by ESOC in the station BRUX, but as noted Section 2, BRUX has been chosen as reference in the daily solutions of ESOC. Apart from this case, the receiver clock estimates using the proposed approach seem to be more stable than the other clock solutions.
We can observe that our solution is slightly noisier than the one obtained by the ACs. The reasons that can explain F I G U R E 5 Receiver clock offsets for permanent stations equipped with a PHM clock: BRUX (top row) and MGUE (bottom row), estimated by the proposed approach (black), ESOC (grey), CODE (red), GFZ (blue), and CNES (green); the left column depicts days 191-199 within the first period in 2017, whereas the right column depicts days 281-289 within the second period in 2017 [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] this fact are: first, our receiver estimates were given every 30 s, whereas the ACs gave their receiver clock estimates every 300 s; and second, our approach used orbits and clock corrections externally computed by ESOC as starting products. Therefore, any difference on the models for computing these products with respect to our model could impact our clock solutions.
The solution of CODE presents a nearly pure linear drift for MGUE during the days 195, 197, 282, 284, and 287, and for BRUX during the day 285, which corresponds to these stations being selected as reference stations in their processing. Figure 6 depicts the ADEV analysis for the two receiver clocks previously shown in Figure 5. In order to compute meaningful ADEVs, we selected days within the first and second period in which the receivers did not experience any reset. The clock results for the stations present lower ADEVs than those obtained for the satellites in Figure 3.
The highest stability (i.e., lowest values of ADEVs) was obtained by ESOC in station BRUX, but it is recalled that it corresponds to the reference station. Moreover, during certain days, other solutions used BRUX or MGUE as a reference station. We can observe that our clock estimates reach 10 −15 in both stations, for long averaging times. However, as in the case of satellites, our ADEVs do not evidence the advantage of the DBD removal observed in Figure 5. This is a consequence of the large number of points that intervened in every plot, whereas decimeter-level DBDs only occur once per day.
In order to quantify the improvement of the receiver clock estimates using the new approach, we have applied the same procedure as the satellites to evaluate the DBD mitigation. That is, to compute the residuals of the first hour of the day with respect to a linear prediction fitted with the clock estimates of the last hour of the preceding day. To enhance the visibility of the comparison, we have calculated the root mean square (RMS) of the residuals during the first hour of each day. We use those RMS values to evaluate the error of the linear extrapolation during the entirety of the two periods. Figure 7 shows the cumulative distribution function (CDF) of the daily RMS of the residuals for the two stations previously depicted in Figure 5. We can observe that the prediction error of gAGE clock estimates is the lowest, with values at the level of 3 cm (approximately 100 ps).
As in the case of satellites, these errors are not far from the nominal accuracy of the IGS final products (approximately 75 ps) and are lower than the errors reported in Ray and Senior (2003) with just GPS satellites.
In contrast, we can see that the clock estimates of the other ACs presented median errors at the level of one decimeter, except ESOC estimates for BRUX station, which was at the same level or better than the clocks determined by our approach. However, BRUX is the reference clock in its daily solutions.

Satellite phase biases
As commented in Section 3, the method estimates the satellite and receiver phase biases simultaneously with clock offsets and ionospheric delays. Because of the one-to-one correlation between the biases and the slant ionospheric delay, their absolute value is not well determined. However, the difference of phase biases (e.g., 2 − 1 , 5 − 1 , or 7 − 5 ) can be accurately determined thanks to the elimination of the ionosphere. In what follows, we evaluate our phase biases with respect to two independent determinations.
First, we use the combination ionosphere-and geometry-free (IGF), , which can be built as the difference of two IF combinations as seen in Li et al. (2020). For instance, using the GPS measurements and frequencies 1 , 2 , and 5 , the IGF for a receiver i, and a satellite j is built as: where Δ = 5 − 1 stands for the phase bias difference between frequencies 5 and 1 . Recall that because of the constraints imposed by Equations (5) and (6), the biases for 1 and 2 cancel for both receiver and satellites in the IF combination of reference. We can proceed in the same manner for Galileo satellites to obtain the differences of phase biases, for instance, for 7 − 1 . Because the unambiguous values are accurately known, Equation (7) allows estimating the phase bias differences of satellites and receivers. In order to solve the rank deficiency of the problem, we have to constrain to zero one of the Δ values in Equation (7), such as the value for the reference receiver CEBR.
In this way, the rest of Δ can be determined in a similar way as the DCB estimation through an ionospheric model [see Sanz et al. (2017)]. During some epochs, it can sometimes occur that there are no observations in common throughout the satellites and stations with the reference. This can generate discontinuities in the clock correction estimates, and in some cases, can impede the estimation. An example of this can be seen in the bottom-left plot of Figure 8, where the red crosses toward the end of day 284.
A second independent estimation of the phase biases can be performed using 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.
The HMW cancels out the terms related to in Equation (1) and Equation (2) and therefore, once having the F I G U R E 8 Time series depicting differences of phase biases for GPS satellites of the Block IIF G08, G09, and G10 in the top-left, top-right, and bottom-left panel, respectively, whereas Galileo satellite E30 (GSAT206 on PHM) is depicted in the bottom-right panel [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] ambiguities resolved, it can be written as: where and correspond to the frequencies 1 and 5 for GPS and 5 and 7 for Galileo. Note the participation of the code pseudorange in Equation (8), unlike Equations (3) and (7) that are based on phase-only measurements.
These comparisons provide us an additional assessment of the suitability of our methodology to handle more than two frequencies. Figure 8 depicts the differences of the estimated phase biases from Equation (3), together with the two independent determinations computed from Equations (7) and (8). In order to magnify this variability, the comparison is done for one Galileo satellite and three GPS satellites of Block IIF, which depict larger fluctuations in comparison with the biases of the Galileo satellite.
We can observe the agreement in the patterns of the differences of estimated phase biases, with those computed directly from measurements, for both the IGF combination and the HMW combination, the latter only biased by the DCB value as expressed by Equation (4). This agreement confirms that the estimation process is fully compatible for three frequencies.
Moreover, it also shows that estimating the hardware delays with a random walk, instead of a bias, is an adequate choice to accommodate periodic variations of Block IIF satellites. Specifically, the time variations of the biases in the third GPS frequency, 5 , which were already pointed out in Montenbruck et al. (2012) to be linked to temperature changes during the eclipsing periods.
It is apparent that the Galileo satellite presents a higher stability with respect to time than the GPS satellites. This is thanks to the fact that Galileo satellites maintain the temperature of the onboard atomic clocks, leaving them highly stable. The drawback of estimating hardware delays such as random walk, instead of biases, is that the use of three frequencies does not strengthen the clock estimates with respect to the estimates using just two frequencies.

CONCLUSIONS
The present contribution shows that the methodology defined in Juan et al. (2020) and Rovira-Garcia et al. (2021) to estimate receiver and satellite clock corrections can be applied using three frequencies. The method requires unambiguous carrier-phases as input measurements. Avoiding pseudorange measurements in the estimation process largely mitigates discontinuities in the clock estimates in the day change.
Unlike other clock solutions, DBDs in our clock estimates are reduced to the level of few centimeters (approximately 100 ps), which is at the same level as the IGS final products (75 ps). This significant DBD reduction of our proposed methodology may benefit the statistical characterization of long-term phenomena correlated with the onboard satellite clocks. The stability of our satellite clock estimates has been measured with the ADEV computed over periods of tens of days, presenting similar stabilities as the clock solutions of other ACs, despite the latter containing DBDs of a few decimeters.
The approach considered here uses several carrierphase measurements without the need of building IF combinations, which is the standard way for estimating clock corrections. Because the estimation process uses three frequencies, clocks offsets are estimated jointly with instrumental phase biases. These hardware delays are estimated as random walk processes, instead of computing them as constant parameters.
This treatment of the process noise allows us to take into account the dependency of the temperature of some of these biases already pointed out by some authors. Differences of our phase bias estimates have been compared with independent determinations built with ionosphere-and geometry-free combinations of GNSS observables. From these comparisons, we had an additional confirmation of the suitability of the methodology.

A C K N O W L E D G M E N T S
The present work was supported in part by the European Space Agency contract (REL-GAL) N.4000122402/ 17/NL/IB, 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. The authors acknowledge the use of data and products provided by the International GNSS Service.