Introduction

Biases (also called hardware or instrument delays) exist for each satellite and receiver for pseudorange and carrier-phase measurements of different frequencies and signals. It is important to note that the true absolute bias of a single signal is unsolvable because bias parameters are one-to-one correlated with clock offsets. The IGS (International GNSS Service) clock products follow the convention that the ionosphere-free linear combination of two pseudorange signals for each constellation (e.g., C1W and C2W for GPS) is equal to zero. Users only need to correct pseudorange biases if they are using different signals (i.e., C1C) by using the corresponding DCB (differential code bias) corrections (Johnston et al. 2017; Kouba and Springer 2001; Montenbruck et al. 2014). The DCB products were later evolved to the OSB (observable-specific signal bias) products by combining the DCB observation model with the IGS clock datum condition (Villiger et al. 2019; Wang et al. 2020).

Satellite phase biases are biases with respect to the IGS clock products, which are in turn aligned to pseudorange measurements. Prior to the development of single-receiver ambiguity fixing techniques, satellite phase biases were not relevant because in dual-frequency processing the constant part of a satellite phase bias can be absorbed by the corresponding float phase ambiguities while the varying part is assumed to go into satellite clock offsets. However, as centimeter- or even millimeter-level absolute positioning results were demanded in various fields, the classical precise point positioning (PPP) (Zumberge et al. 1997) was combined with ambiguity resolution and evolved to PPP-AR (ambiguity resolution) (Laurichesse et al. 2010; Lin et al. 2023; Mervart et al. 2008) and PPP-RTK (real-time kinematic) (Li et al. 2017; Zhang et al. 2019).

To support PPP-AR applications, dedicated satellite phase bias products are needed. Ge et al. (2005) provided users with satellite wide-lane and narrow-lane biases, while Loyer et al. (2012) and Laurichesse et al. (2009) incorporated satellite narrow-lane biases into satellite clock offsets. For the later, the datum of satellite clock products is defined by two phase signals (e.g., L1W and L2W for GPS), which is slightly inconsistent with the legacy IGS clock products. Schaer et al. (2021) combined satellite DCBs, wide-lane biases, narrow-lane biases together with the IGS clock datum condition and determined OSBs on each signal. A dedicated exchange format for OSBs was later agreed as SINEX_BIAS version 1.00 (Schaer 2016). Concurrent with the introduction of the IGS20 frame (starting from GPS week 2238), a new bias convention was introduced, which requires consistent use of an antenna phase center (APC) model in the Hatch–Melbourne–Wübbena (HMW) (Hatch 1983; Melbourne 1985; Wubbena 1985) and geometry-free linear combinations, and also identifies the applicable model as part of the DCB/OSB product (Schaer 2022).

More and more analysis centers are providing satellite phase bias products in the frame of IGS or as a local service (Duan et al. 2023, 2021; Geng et al. 2019; Li et al. 2018; Strasser et al. 2019; Uhlemann et al. 2015). For the users, it is crucial to know if the available satellite phase bias products are sufficiently precise and consistent to conduct PPP-AR applications. The commonly used method is to check the performance of ambiguity fixing rates and positioning results in PPP-AR applications, which, however, are not sensitive to phase bias errors of a specific satellite. In order to obtain more robust satellite products for PPP-AR applications, Banville et al. (2020) made a preliminary combination of satellite clock/bias products from individual analysis centers over one week. However, potential outliers in satellite phase biases (which are almost inevitable in long-time periods) need further discussions in the combination task.

With the above background, this work aims at proposing a method completely independent of measurements to detect outliers in satellite phase biases. The methodology developed for this purpose is introduced in the following section and subsequently applied to different sets of orbit, clock, and bias products for GPS and Galileo satellites available from IGS analysis centers. It offers a lean alternative to PPP-based validation of biases and can be advantageously applied for quality control and combination of GNSS products. Following the analysis and discussion of the achieved results, summary and conclusions are presented.

Methodology

Considering dual-frequency GPS observations from semi-codeless P(Y)-code tracking without loss of generality, the observation model for the ionosphere-free combination of carrier-phase ranges can be expressed as:

$$\begin{aligned} {\text{IF}}\left( {L_{{{\text{L1W}}}} ,L_{{{\text{L2W}}}} } \right) & = r^{s} - r_{r} + c\left( {{\text{d}}t_{r} - {\text{d}}t^{s} } \right) + {\text{IF}}\left( {\zeta_{r,1}^{s} ,\zeta_{r,2}^{s} } \right) + c \cdot {\text{IF}}\left( {b_{{{\text{L1W}}}}^{s} ,b_{{{\text{L2W}}}}^{s} } \right) + c \cdot {\text{IF}}\left( {b_{{r,{\text{L1W}}}} ,b_{{r,{\text{L2W}}}} } \right) \\ & \quad + \frac{{cf_{2} }}{{f_{1}^{2} - f_{2}^{2} }}N_{{{\text{WL}}}} + \frac{c}{{f_{1} + f_{2} }}N_{1} + T_{r}^{s} + {\text{IF}}\left( {\lambda_{1} , \lambda_{2} } \right) \cdot \omega_{r}^{s} + \varepsilon + \smallint \\ \end{aligned}$$
(1)

(Hauschild 2017). Here IF denotes the ionosphere-free linear combination for frequencies \(f_{1}\), \(f_{2}\) (wavelength\({\lambda }_{1}, {\lambda }_{2}\)) of two phase measurements\({L}_{{\text{L}}1{\text{W}}},{ L}_{{\text{L}}2{\text{W}}}\), two satellite and receiver phase biases \({b}_{{\text{L}}1{\text{W}}}^{s},{ b}_{{\text{L}}2{\text{W}}}^{s},{b}_{r,{\text{L}}1{\text{W}}},{b}_{r,{\text{L}}2{\text{W}}}\), and the related phase center offset/variation corrections \({\zeta }_{r,1}^{s},{ \zeta }_{r,2}^{s}\). \({r}^{s}\) and \({r}_{r}\) denote satellite and receiver positions, \({\text{d}}{t}_{r}\) and \({\text{d}}{t}^{s}\) represent receiver and satellite clock offsets, \({T}_{r}^{s}\) denotes the tropospheric range delay, \({\upomega }_{r}^{s}\) represents the phase windup correction (Wu et al. 1992), and c denotes the speed of light in vacuum. \({N}_{{\text{WL}}}\) denotes the integer wide-lane ambiguity and \({N}_{1}\) represents the narrow-lane ambiguity. Finally, \(\varepsilon\) and \(\epsilon\) describe contributions of modeling errors and measurement errors, respectively.

Modeled contributions include satellite orbits, clock offsets and biases, which are commonly provided by the respective products of individual IGS analysis centers. The consistency of different products may be assessed by considering differences of (1) for the respective products at common epochs. As a special case, the consistency of products from a single analysis center may be evaluated at the day boundary by considering the difference (\({\Delta }_{AB}\)) of (1) between the 24 h epoch of day \(n\) and the 0 h epoch of day \(n+1\).

$$\begin{aligned} 0 & = \Delta_{AB} r^{s} - r_{r} + c\left( {\Delta_{AB} {\text{d}}t_{r} - \Delta_{AB} {\text{d}}t^{s} } \right) + c \cdot \Delta_{AB} {\text{IF}}\left( {b_{{{\text{L1W}}}}^{s} ,b_{{{\text{L2W}}}}^{s} } \right) \\ & \quad + c \cdot \Delta_{AB} {\text{IF}}\left( {b_{{r,{\text{L1W}}}} ,b_{{r,{\text{L2W}}}} } \right) + \frac{{cf_{2} }}{{f_{1}^{2} - f_{2}^{2} }}\Delta_{AB} N_{{{\text{WL}}}} + \frac{c}{{f_{1} + f_{2} }}\Delta_{AB} N_{1} + \Delta_{AB} \varepsilon \\ \end{aligned}$$
(2)

Linearizing the contribution of the satellite position, the geometric range may be expanded as

$$\begin{aligned} \Delta_{AB} r^{s} - r_{r} \approx \Delta_{AB} \left( {e^{T} r^{s} } \right) & = e^{T} e_{r} \Delta_{AB} r + e^{T} e_{a} \Delta_{AB} a + e^{T} e_{c} \Delta_{AB} c \\ & = \Delta_{AB} r + \left[ {\left( {e^{T} e_{r} - 1} \right)\Delta_{AB} r + e^{T} e_{a} \Delta_{AB} a + e^{T} e_{c} \Delta_{AB} c} \right] = \Delta_{AB} r + \Delta_{AB} \varepsilon_{{{\text{LOS}}}} \\ \end{aligned}$$

with \(e=\frac{{{\varvec{r}}}^{s}-{{\varvec{r}}}_{r}}{\Vert {{\varvec{r}}}^{s}-{{\varvec{r}}}_{r}\Vert }\) denoting the line of sight unit vector, \({e}_{r}\), \({e}_{a}\), \({e}_{c}\) denoting unit vectors in the radial, along-track, and cross-track direction, and \(r, a, c\) denoting the corresponding components of the satellite position vector. The primary contribution to the range difference stems from the difference of the radial orbit component, while smaller, line-of-sight dependent contributions result from the along-track and cross-track differences. Considering a global average over the station in the field of view of the given satellite, roughly 14% of the along-track and cross-track orbit errors contribute to the RMS (root mean square) of the associated line-of-sight range error \({{\Delta }_{AB}\varepsilon }_{{\text{LOS}}}\) for GNSS satellites in medium Earth orbit (Montenbruck et al. 2018). Representative RMS values of the non-line-of sight position errors \(\sqrt{{\left({\Delta }_{AB}a\right)}^{2}+{\left({\Delta }_{AB}c\right)}^{2}}\) over daily intervals for different IGS products amount to 9 mm for GPS and 12 mm for Galileo satellites yielding \({{\Delta }_{AB}\varepsilon }_{{\text{LOS}}}\) values of about 1.3 and 1.7 mm, respectively. Orbit day-boundary discontinuities of a single analysis center can be two times larger. For practical purposes, we lump \({{\Delta }_{AB}\varepsilon }_{{\text{LOS}}}\) into \({\Delta }_{AB}\varepsilon\) and consider only the combined error in the subsequent analysis.

Equation (2) still contains a variety of receiver specific contributions that hamper a direct analysis of the orbit/clock/bias product consistency. To eliminate most of these contributions, we form a further difference of (2) between satellite \({\text{i}}\) and a selected reference (\({\text{ref}}\)) satellite.

$$\begin{aligned} - \Delta_{AB} r^{{i - {\text{ref}}}} - c\Delta_{AB} {\text{d}}t^{{i - {\text{ref}}}} - c\Delta_{AB} {\text{IF}}\left( {b_{{{\text{L1W}}}} ,b_{{{\text{L2W}}}} } \right)^{{i - {\text{ref}}}} - \frac{{cf_{2} }}{{f_{1}^{2} - f_{2}^{2} }}\Delta_{AB} N_{{{\text{WL}}}}^{{i - {\text{ref}}}} & = \Delta_{AB} \frac{c}{{f_{1} + f_{2} }}N_{1}^{{i - {\text{ref}}}} \\ & \quad + \Delta_{AB}^{{i - {\text{ref}}}} \varepsilon \\ \end{aligned}$$
(3)

Within this expression, \({\Delta }_{AB}{N}_{\text{WL}}^{i-{\text{ref}}}\) denotes the product-minus-product and satellite-minus-satellite double difference of the wide-lane ambiguities that would be determined by a user of the two products. Its value may be obtained by forming, in analogy with (1), the observation model for the HMW combination of dual-frequency pseudorange and phase observations. The wide-lane ambiguity double difference

$$\Delta_{AB} N_{{{\text{WL}}}}^{{i - {\text{ref}}}} = {\text{round}}\left[ {\frac{{f_{1}^{2} - f_{2}^{2} }}{{cf_{2} }}\Delta_{AB} {\text{HMW}}\left( { b_{{{\text{C1W}}}}^{s} , b_{{{\text{C2W}}}}^{s} , b_{{{\text{L1W}}}}^{s} , b_{{{\text{L2W}}}}^{s} , b_{{{\text{C1W}}}} ,b_{{{\text{C2W}}}} ,b_{{{\text{L1W}}}} ,b_{{{\text{L2W}}}} } \right)^{{i - {\text{ref}}}} } \right]$$
(4)

can then be expressed as the double difference of the HMW combination of satellite code and phase biases \({b}_{{\text{C}}1{\text{W}}}^{s}, {b}_{{\text{C}}2{\text{W}}}^{s}, {b}_{{\text{L}}1{\text{W}}}^{s}, {b}_{{\text{L}}1{\text{W}}}^{s}\) in the two products.

Overall, the left-hand part of (3) represents a quantity that can be computed directly from satellite orbit, clock and code/phase bias products. Its value should be equal to an integer number of narrow-lane wavelength plus noise, and the fractional part of a narrow-lane wavelength represents a measure of the inconsistency of the two satellite products A and B. We name this fractional term as Differences of satellite Orbits, Clock offsets and narrow-lane Biases (DOCBs). In order to minimize the impact from the selected reference satellite, we employ a zero-mean condition as a further step for the DOCBs per constellation. DOCBs can either be considered to assess the consistency of products from a single analysis center across multiple days or to assess the daily consistency of products from two different analysis centers. Both applications are independently discussed in the subsequent sections.

Consistency of satellite phase bias products between two consecutive days

DOCB values at the midnight epoch for satellite products from a single analysis center are computed to check the consistency of satellite products between two consecutive days. Figure 1 illustrates the computation strategy same as also given in equation (3). It is important to note that although only the difference at the midnight epoch is computed DOCB values in this section represent the consistency of satellite phase biases between two days as satellite phase bias products available in IGS are given as daily bias products.

Fig. 1
figure 1

Schematic display of the DOCB computation at the midnight epoch between two consecutive days

Full size image

CODE (Center for Orbit Determination in Europe) final, CODE rapid, CNES/CLS (Center National d’Etudes Spatiales/Collecte Localisation Satellites) final and WUHN (Wuhan University) MGEX rapid GPS and Galileo satellite orbit, clock offset and bias products are analyzed. Time periods range from day 001 to day 180 of the year 2023, which covers a wide range of Sun elevation angels above the orbital plane (\(\beta\)-angle) for each satellite. Satellite product information is shown in Table 1. Biases of GPS C1W, C2W, L1W, L2W and Galileo C1C, C5Q. L1C, L5Q signals are assessed. APC (antenna phase center) models for bias products are indicated by each analysis center in the header of the BIAS file.

Table 1 Satellite product information
Full size table

CODE (Schaer et al. 2021) and WUHN (Geng et al. 2022) analysis centers provide GNSS satellite orbits and clock offsets at the midnight epoch in their daily final (rapid for WUHN) products, which allows the computation of day boundary discontinuities for both satellite orbits and clock offsets between two consecutive days without extrapolations. Satellite code and phase biases are provided as daily constant values, assuming that potentially varying components are, at least partly, absorbed by satellite clock offsets. On the other hand, CODE rapid satellite clock products, and CNES/CLS final satellite orbit and clock products both do not contain the midnight epoch. In order to compute DOCBs for these two products, we first extrapolate CNES/CLS satellite orbits to the midnight epoch using the dynamic orbit fitting method (Dach et al. 2015). Then, we compute satellite clock differences between epoch 24:00:00 and 23:59:30 of the CODE final products for each satellite. Finally, we add the computed clock differences to CODE rapid and CNES/CLS final clock products at epoch 23:59:30 and deduce satellite clock offset for each satellite at epoch 24:00:00. We have proven that this orbit/clock extrapolation method introduces additional noise of 0.01–0.02 narrow-lane cycles.

The left figures of Figs. 2 and 3 show DOCB values as a function of time for GPS and Galileo satellites, while the right figures show DOCB values as a function of the \(\beta\)-angle. Results of days 17 and 18 in the CNES/CLS products have larger outliers for almost all the satellites due to unknown reasons and are therefore removed from the analysis. About 0.2, 1.1, 2.0 and 0.1% of the total DOCB values for the GPS satellites exceed 0.15 narrow-lane cycles for CODE final, CODE rapid, CNES/CLS final and WUHN MGEX rapid satellite products, respectively. Galileo satellite phase bias products show fewer outliers, about 0.1, 0.9, 0.4 and 0.1% of the total DOCB values exceed 0.15 narrow-lane cycles for CODE final, CODE rapid, CNES/CLS final and WUHN MGEX rapid satellite products, respectively. DOCB values of satellite products from all the analysis centers show a slight dependency on the beta angle due to deficiencies in the used radiation pressure models.

Fig. 2
figure 2

GPS day-boundary DOCB values as a function of time (left) and Sun beta angle (right) for different analysis center products, each shifted by 0.5 narrow-lane cycles, std denotes the standard deviation of the individual time series

Full size image
Fig. 3
figure 3

Galileo day-boundary DOCB values as a function of time (left) and Sun beta angle (right) for different analysis center products, each shifted by 0.5 narrow-lane cycles, std denotes the standard deviation of the individual time series

Full size image

We would like to mention that the computed DOCB values rely on satellite clock offsets of a single epoch, which could enlarge the noise of the computed DOCB values. The STDs (standard deviation) of DOCBs over 180 days for CODE final and WUHN products are about 0.02–0.03 narrow-lane cycles and are about 0.04–0.06 narrow-lane cycles for CODE rapid and CNES/CLS products. The main reason could be that CODE final and WUHN solutions employ a ground network of more than 300 stations, while CODE rapid and CNES/CLS use a network of about 100 stations. CODE rapid products employ similar orbital models as those in CODE final products but are based on a sparser ground network. Also, the clock extrapolation to the midnight epoch necessary for CODE rapid and CNES/CLS products contributes to the additional noise.

CNES/CLS satellite phase bias products are validated by their provider by checking the consistency of the resolved integer ambiguities in the overlapping time periods between two consecutive days based on measurements (Katsigianni 2019; Katsigianni et al. 2019). All the detected outliers in satellite phase biases are listed in the GRG_ELIMSAT_all.dat file,Footnote 1 which can be taken as an external reference to check the performance of our method. Figure 4 illustrates the computed DOCBs of CNES/CLS final satellite products before (GPS, GAL) and after (GPSn, GALn) removing all the satellites that are indicated in the GRG_ELIMSAT_all.dat file. We find that the outliers that exceed 0.15 narrow-lane cycles drop from 2.0 to 0.2% for GPS and from 0.4 to 0.1% for Galileo, which demonstrates that our method is able to independently identify almost all the outliers indicated in the GRG_ELIMSAT_all.dat file.

Fig. 4
figure 4

DOCB results of CNES/CLS satellite products before (GPS, GAL) and after (GPSn, GALn) removing outliers in GRG_ELIMSAT_all.dat, each shifted by 0.5 narrow-lane cycles

Full size image

Comparing satellite phase bias products from two analysis centers

The IGS started combining GPS satellite orbits from analysis centers in 1993 (Beutler et al. 1995; Johnston et al. 2017; Kouba and Springer 2001). Since then, a lot of efforts have been invested to make the combined satellite products more precise and robust (Banville et al. 2020; Chen et al. 2023; Griffiths 2019; Mansur et al. 2020; Sośnica et al. 2020; Zajdel et al. 2023).

Same as for satellite orbit and clock offset combinations, outlier detection is crucial for the combination of satellite biases as well. In this section, we show the performance of satellite orbit, clock offset and bias differences between two analysis centers before and after removing the detected outliers. Equation (3) is also valid in the comparison of daily satellite phase bias products between two analysis centers or campaigns. The difference is that \({\Delta }_{AB}{r}^{i-{\text{ref}}}\) and \({\Delta }_{AB}{\text{d}}{t}^{i-{\text{ref}}}\) in (3) here denote the mean differences of satellite orbits and clock offsets of two products over one day (Duan and Hugentobler 2021).

Figure 5 shows daily DOCB results as a function of \(\upbeta\)-angle by comparing cor, grg and wum GPS (left figure) and Galileo (right figure) satellite products to cod products. Outliers are not negligible in all the solutions, and GPS satellites exhibit clearly more outliers than Galileo satellites. Now, we exclude all the detected outlier-satellites identified with the aforementioned method in Figs. 2 and 3 from Fig. 5 and show the remaining DOCB results in Fig. 6. We find that outliers do not vanish but their number is greatly reduced. It is worth to notice that cod-grg and cod-wum results show large and dense outliers for satellites with the \(\upbeta\)-angle close to zero, which could be attributed to the differences in the employed satellite attitude and the radiation pressure models.

Fig. 5
figure 5

Daily DOCB results by comparing cor, grg and wum GPS (left figure) and Galileo (right figure) satellite products to those of cod

Full size image
Fig.6
figure 6

Daily DOCB results by comparing cor, grg and wum GPS (left figure) and Galileo (right figure) satellite products to those of cod after removing the detected outliers in Figs. 1 and 2

Full size image

Both the comparison of the detected outliers with the GRG_ELIMSAT_all.dat file from CNES/CLS and the reduction of the number of large differences in satellite bias comparison between two analysis centers confirm that our method is efficient. This is helpful for the PPP-AR users and for the combination tasks because it makes it possible to decide which satellite phase biases should be used.

Noise Characteristics of DOCBs for one analysis center over time

As studied by Montenbruck et al. (2012), the L1/L5-minus-L1/L2 clock differences for the GPS Block IIF-1 satellite show clear time variations related to the change of the onboard temperature. This indicates that phase biases of GPS Block IIF satellites for frequency L1, L2 and L5 might vary along with time. Similar work was done in Duan et al. (2023) to assess the E1/E5a-minus-E1/E5b, E1/E5a-minus-E1/E5 and E1/E5a-minus-E1/E6 clock differences for all the Galileo satellites. Different than for the GPS Block IIF satellites, clock differences between any dual-frequency for Galileo satellites are constant over one day, which implies that satellite phase biases of each Galileo frequency are constant over time. In the dual-frequency case, satellite phase biases essentially represent the difference between code and integer phase clock. However, this difference is canceled in the DOCB.

DOCBs (\({\Delta }_{AB}^{i-{\text{ref}}}\)) determined in (3) contains the non-line-of sight orbit errors, modeling differences and the phase measurement noise in the determination of satellite products from individual analysis centers. For satellite products from a single analysis center, the mean of the DOCBs at the day-boundary per satellite over time should thus be zero. In this section, we check the noise characteristic of DOCBs over time. DOCB results of cod and wum are employed because a number of GPS and Galileo satellite DOCBs are free of gaps and outliers, which is the precondition to conduct the Allan Deviation analysis (Riley and Howe 2008).

Figure 7 shows the behavior of the Modified Allan Deviation of the DOCBs. Up to integration time intervals of 10 days, DOCBs from both analysis centers show a white noise behavior. However, for time intervals longer than 10 days, DCOBs show different performance between cod and wum products, which could be due to the use of different ground networks (wum uses more stations than cod) or different radiation pressure models. In general, DOCBs for Galileo satellites are more stable than those of GPS satellites for both solutions. One of the reasons could be due to different impact of a varying ground network and the receiver-dependency on satellite biases. The observed white noise characteristic is not relevant for PPP-AR applications but could be helpful for the product providers to determine satellite phase bias products with appropriate constraints to the solutions of the previous day.

Fig. 7
figure 7

The modified Allan deviation of day-boundary DOCBs over 180 days

Full size image

Summary and conclusions

GNSS satellite phase bias products are available from several IGS analysis centers to support PPP-AR applications. However, it is difficult to assess directly the precision of satellite phase bias products. The commonly used approach is to assess the positioning performance in PPP-AR applications. Aside from the computational load and complexity of a PPP-based validation, errors or outliers in phase bias of a specific satellite are not directly visible in this process but lumped into the overall observation residuals. This contribution presents an alternative method which is computationally lean and independent of ground measurements to detect outliers and inconsistencies in satellite orbit, clock and phase bias products. The goal is to provide PPP-AR users with quality-controlled satellite phase bias products and to make the future combination of satellite phase bias products less affected by outliers. In accord with the current practice of daily products assuming constancy of satellite phase biases over 24 h, DOCBs are computed at midnight epoch between two consecutive days to assess the consistency of consecutive data products. The method would likewise be applicable, though, for future ultra-rapid products with sliding shifts of the processing intervals or phase bias products with sub-daily validity intervals.

Results of 180 days show that about 0.2, 1.1, 2.0 and 0.1% of the total DOCB values for GPS satellites exceed 0.15 narrow-lane cycles for CODE final, CODE rapid, CNES/CLS final and WUHN MGEX rapid satellite products, respectively. Galileo satellite phase bias products show fewer outliers, about 0.1, 0.9, 0.4 and 0.1% of the total DOCB values exceed 0.15 narrow-lane cycles for CODE final, CODE rapid, CNES/CLS final and WUHN MGEX rapid satellite products, respectively. Outliers detected in the CNES/CLS products are in line with those found by their own assessment as indicated in the GRG_ELIMSAT_all.dat file. By removing all the detected outliers in individual products, comparisons of daily satellite phase bias products are less affected by large differences. However, notable inconsistencies in satellite orbit, clock and bias values remain during the eclipse season especially when the \(\upbeta\)-angle is close to zero due to the deficiencies of the employed radiation pressure models. The developed method may be used by product providers to identify and eliminate (or mark) outliers in the estimated phase biases and generate a list similar to the CNES/CLS GRG_ELIMSAT_all.dat file. Eventually, the product combination provider and the PPP-AR users may also verify, which satellite phase biases should be used.