Multi-GNSS PPP solutions with different handling of system-specific receiver clock parameters and inter-system biases

In the multi-GNSS solutions integrating GPS, GLONASS, Galileo, and BeiDou, the receiver clock may be treated twofold; the clock parameter may be estimated for each GNSS separately or the common clock can be estimated, e.g., for GPS, with inter-system biases (ISBs) for other systems. The latter strategy reduces the number of estimated independent clock parameters per epoch almost by a factor of four because the clock parameters are estimated epoch-wise, whereas ISBs are estimated as constant values for the entire day or month. Due to the discontinuities in reference satellite clocks, the estimated ISBs and receiver clock parameters have also to be reinitialized at day boundaries. This raises questions about whether only the common clock has to be reset or all ISB values and what is the impact of the reinitialization of clock parameters with covariance values when estimating system-specific clock parameters. We analyze the effects of different types of stochastic modeling applied to the parameters of clocks and ISBs. In this study, we test five different strategies to clock handling in multi-GNSS kinematic Precise Point Positioning derived continuously for one month. We found that two solutions can be considered equivalent: (1) estimating system-specific clocks and (2) estimating the common clock with ISB and resetting at day boundaries the common clock parameter and ISBs. Oppositely, resetting only the common clock parameter or assuming that the ISB keep their stabilities over long periods is insufficient to obtain superior results of station coordinates and reliable time transfer results.


Introduction
The development of Global Navigation Satellite Systems (GNSS) has contributed to the use of multi-GNSS as a common solution employing GPS (G), GLONASS (R), Galileo (E), and BeiDou-3 (C) (Montenbruck et al. 2017;Closas and Gusi-Amigo 2017). The advantage of such a solution is that more observations can be selected and, as a result, the accuracy of the results is increased. Precise Point Positioning (PPP, Kouba et al. 2017) achieves accuracies down to the millimeter range, depending on the measurement mode used (Bisnath and Gao 2009). Multi-GNSS PPP improves the accessibility of the solution by improving the tracking geometry, which is especially important in areas with a highly obscured horizon (Xia et al. 2019).
When processing data from a single system, the hardware delay is assimilated to the receiver clock and does not significantly affect the position estimation. In contrast, inter-system biases (ISBs) or system-specific clock differences must also be considered when using multi-GNSS data. ISB includes the hardware delay caused by different signal decrypting schemes and circuits for different satellite systems (Chen et al. 2015;Zhou et al. 2019;Hong et al. 2019;Mi et al. 2019). In multi-GNSS PPP, the receiver clock parameter can be estimated in two ways. The clock can be estimated separately for each system, or a single common clock parameter can be determined for all GNSS together with ISBs. Receiver clocks estimated as white noise for a single GNSS system show substantial noise (Weinbach and Schön 2011;Ge et al. 2018;Krawinkel and Schön 2021). Furthermore, due to the discontinuity at the day boundary, a receiver clock discontinuity can be observed (Guo and Zhang 2014;Lyu et al. 2019;Qin et al. 2020), which is caused by changing the tracked satellites and selecting a different reference clock. Such adverse effects directly affect the precision of positioning using the PPP technique.
We analyze the stability of GNSS station coordinates based on different handling of the receiver clock parameters and ISBs in multi-GNSS solutions. We show the impact of different clock and ISB handling on the repeatability of estimated station coordinates and time transfer results. Moreover, we analyze the impact of resetting clock parameters and ISBs at the day boundaries. The analysis is conducted for selected International GNSS Service (IGS, Johnston et al. 2017) station coordinates estimated with multi-GNSS PPP solutions based on the following GNSS: G, R, E, and C with considering different types of clocks, such as hydrogen maser (HM), cesium (Cs), rubidium (Rb), and crystal oscillators (XO). All selected stations belong to the multi-GNSS Experiment (MGEX) pilot project (Montenbruck et al. 2017) and ensure the tracking of four GNSS systems to analyze system-specific clock errors. The consistency of the MGEX orbit and clock products has recently been increased (Steigenberger and Montenbruck 2020), however, differences between GNSS systems still have to be properly handled in multi-GNSS solutions.
This study considers (1) the receiver clock parameters estimated independently for each system and (2) the common receiver clock parameters and ISBs estimated. The ISBs are estimated on a daily basis as well as at monthly intervals. The estimation of system-specific clocks per epoch increases the number of estimated independent clock parameters almost by a factor of four with respect to the solution with the estimation of the common clock and ISBs for other GNSSs.
The next two sections describe the methodology and the relationship between determining clock parameters independently for each GNSS system and the common clock with ISB. We have indicated strategies for different types of stochastic modeling superimposed on the parameters of clocks and ISBs. The results and main advantages of both solutions on the station coordinate repeatability and time transfer results are then shown in the following sections.

Methodology
We process multi-GNSS dual-frequency pseudorange (code) P and carrier phase L observations using the undifferenced uncombined model (Schönemann 2014) of the PPP technique: with where s denotes the satellite number and M is the corresponding GNSS ( M ∈ {G, R, E, C});i is the number of the frequency f ; s 0 is the geometric distance between the position of the satellite s antenna phase center (X s , Y s , Z s ) and the a priori position of the receiver r antenna phase center X r , Y r , Z r ; c is the speed of the light; t s and t M r are satellite and system-specific receiver clock offsets, respectively; T is the troposphere delay and m s is its corresponding mapping function; Δ s P and Δ s L are satellite, receiver, and site displacement effect corrections for P and L measurements, respectively; ISB G M is the receiver inter-system bias between M and G (for M = G it is set to 0);I s is the slant ionosphere delay; N s i is the carrier-phase ambiguity that corresponds to its original wavelength i ; e s r is the direction vector and X r is the position correction vector. For system-specific clock estimation, all ISB G M = 0. For the solutions with ISB G M determination, one common t M r parameter is derived, i.e., t G r = t R r = t E r = t C r . The analysis includes clocks connected to receivers at 16 selected IGS MGEX stations. Each selected station tracks at least four systems: G, R, E, C. Figure 1 shows the IGS stations analyzed with an indication of the clock type. Table 1 provides the most relevant information for each station. The stations are sorted by clock type.
Multi-GNSS stations were selected based on analyses of the quality of clock stability for different types of clocks and the distribution of these stations in multiple locations (Mikoś et al. 2023). All PPP solutions are determined in the GNSS-WARP (Wroclaw Algorithms for Real-Time Positioning) software (Hadaś 2015). Undifferenced and uncombined dual-frequency GNSS observations are processed and carrier-phase ambiguities are estimated as float values. The solution is computed based on the sequential least squares method (Hadas et al. 2019). We derive the solution using the final orbits of satellites and clocks from the Center for Orbit Determination in Europe (CODE, Prange et al. 2017) MGEX analysis center. The kinematic (1) (4) X r = X r Y r Z r T PPP mode was used for the analyses, as indicated in the description of the strategy in Table 2 the station coordinates and troposphere parameters are common to all GNSS systems, whereas clocks and ISB differ in particular solution strategy.

Multi-GNSS PPP model based on an independent determination of the clock parameter in each system
Determining the clock parameters independently in each system increases the number of independent unknowns up to four in each estimated epoch when using a quad-system solution. The number of estimated independent clock parameters  per day equals 4 * (24 * 60 * 2) = 11, 520 for each station, assuming a 30 s GNSS data sampling rate. Despite the fact that the GNSS receiver is physically connected to one clock, the system-specific clock parameters differ because they typically absorb some differences between GNSS systems, temperature and electromagnetic variations (Mikoś et al. 2023).
The system-specific difference may also emerge from observation geometry and determinability of clock parameters.

Multi-GNSS PPP model based on ISB estimation for clock parameter
If we estimate the common clock parameter epochwise together with the ISBs once a day, the number of independent parameters will be reduced, equaling 3 + (24 * 60 * 2) = 2, 883 for 30 s sampling and four GNSSs. In this study, we use GPS as the common clock, which is the reference for R, E, and C when estimating ISBs. The most important question concerning the ISBs is their short-and long-term stability over time. In the following sections, the advantages and disadvantages of the solution with the determination of the common clock parameter and the ISBs will be examined in relation to other GNSS. We also tested different strategies for resetting the ISB and clock parameters i.e. (1) resetting them on a daily basis, (2) resetting only the common clock and keeping ISB estimates continuously over the entire month, as well as (3) solution without resetting the ISB and the common clock parameter.

Solution strategies
The clock parameter is analyzed using five different solution strategies to find the best multi-GNSS PPP processing strategy, as shown in Table 3.
In GREC, the clock parameters are estimated independently for each system. The GREC clocks are estimated as white noise without any time correlations considered. In GREC + res, the clock parameters are estimated as in GREC, along with resetting these parameters at the boundary of the day. During the resetting procedure, not only is the parameter reinitialized but also all off-diagonal elements in the covariance matrix. This means that the variances are set to the a priori values and the correlations between the clock parameters and other parameters are set to zero at the day boundaries. Resetting should account for a change of the reference clock in CODE orbit products because for each day, a different reference clock is selected, which should be fully uncorrelated with the clock values from the previous day. Please note that in GREC and GREC + res strategies, the number of estimated independent parameters is the same.
In G + ISB, the common clock parameter is estimated together with the ISBs, which are estimated as constants. This strategy assumes that the differences in ISBs are constant over long spans, a month in this case, and can be absorbed by ISBs. The receiver is connected to one physical clock, whose values are estimated epoch-wise by the common clock parameter.
In G + ISB + res, the common clock parameter is estimated together with the daily ISBs estimated daily parameters and all these parameters and covariances are reset at the boundary of the day. This strategy assumes that the ISB vary over time, and thus, ISBs should be estimated on a daily basis, whereas clock discontinuities should be accounted for by resetting the common clock. The G + ISB + res strategy slightly increases the number of independent estimated parameters by three ISB values per day when compared to G + ISB.
In G + ISB + resG, the common clock parameter is estimated with the ISBs and resetting the common clock parameter at the boundary of the day. Thus, ISBs are estimated as constant parameters, monthly-wise, whereas discontinuities in the reference clock should be absorbed by resetting the common clock parameter.

Results
For the following analyses, the repeatability of station coordinates is evaluated by the interquartile range (IQR) and standard deviation (STD) of the estimated coordinates in each strategy, as well as by the fast Fourier transform (FFT). The quality of the frequency stability of the clock parameters is evaluated using the Modified Allan Deviation (MDEV).

Results of the experiment for AREG and BRUX stations
The results shown in this section are based on two stations with Rb and HM clocks, i.e., AREG and BRUX, respectively. Figure 2 shows the results for GREC and G + ISB solutions. The upper part shows the values for clock estimates, whereas the lower part shows the epoch-wise clock differences between adjacent epochs. These are the differences of the values estimated every 30 s for each system. All values are shown in metric units; for the conversion to time units, one has to divide results by the speed of light.
For the AREG station with Rb clock, the raw clock estimates are very close to each other in the GREC solution. Different results are clearly visible for each system for the BRUX station equipped with HM. Comparing the results for the two solutions shown in Fig. 2, the STD of the G clock for the AREG station reaches 3.550 m and for the common clock reaches 3.591 m for GREC and G + ISB, respectively. For station BRUX, the STD of the G clock reaches 0.359 m and for the common clock 0.420 m for GREC and G + ISB, respectively. Therefore, the G clock estimates from the GREC solution are more stable compared to the common clock estimate from the G + ISB.
The epoch difference in the GREC solution shows the highest short-term stability for the G and E systems. For the AREG station, the STD value reaches 0.031, 0.034, 0.031, and 0.032 m for G, R, E, and C, respectively, while for station BRUX the respective STDs are 0.008, 0.012, 0.008, and 0.013 m. Comparing the results of the two solutions shown in Fig. 2 bottom, the STD of the G clock for station AREG reaches 0.031 m and for the common clock 0.033 m for GREC and G + ISB, respectively, while the STD value for BRUX station is 0.008 and 0.010 m for GREC and G + ISB. Figure 3 shows the ISB results for two solutions. The left column shows the results for the GREC solution with reconstructed epoch-wise ISB values as differences between system-specific R, E, C clocks and G clocks. The scale of the vertical axis is identical in both cases. The values for each clock parameter were reduced by the value of the last epoch to increase the readability of Fig. 3.
The right column shows the results for the G + ISB solution. Due to the significant difference in the scatter of the results in the early and late epochs, a zoom is shown for successive solutions. The vertical axis for the AREG station is identical for R-G and E-G, whereas it is different for C-G due to a large scatter. For the BRUX station results, all solutions have identical axis scales. The labels on the horizontal axis for the zoomed figures are identical to those from the main figures. In the GREC solution, the AREG station obtains larger ISB values for all three systems than BRUX. The STD value is 0.238, 0.197, and 0.579 m for R-G, E-G, and C-G, respectively, for AREG. The BRUX station results do not have large variabilities for R-G and E-G; the STD value is 0.198, 0.192, and 0.579 m for R-G, E-G, and C-G, respectively in GREC. In both cases, the C-G  solution showed a large discontinuity on April 21, which worsened the STD statistics. For the G + ISB solution, the results for the BRUX stations converge much faster than for AREG. For R-G and E-G, the values converge to a certain level already on the first day, around the second hour and keep their stability till the end of the month. For C-G, an asymptotic convergence for neighboring epochs starts from the second epoch and stabilizes after several hours. In the case of the results for the AREG station, a similar relationship can be seen for E-G as for the BRUX station. However, there is a disturbance in the solution for the R-G and C-G results, which may be due to the quality of the observations, changes in satellite clocks and satellite-specific ISBs during the solution. In the further part of the solution, starting from day three, there is almost no variability in the results and the ISB values converge.
The stability of the clocks is shown in Fig. 4, based on MDEV analyses for the GREC and G + ISB solutions. For the GREC solution and for AREG with the Rb clock, the amount of noise is so large that the results shown in each system have almost identical stability. For the BRUX station with an HM clock, the stability of the result in each system shows different results. In particular, the results for system C show worse stability than the other systems. However, the results for G, R and E are similar, except for the first epochs in R, which are of special interest for the time transfer. When comparing the results for G in both solutions, they show very similar results, however, the very first epochs seem to be more stable in GREC. Figure 5 shows the MDEV error values for both solutions. For the AREG station, the convergence time to reach the MDEV error level of 10 −13 takes about 6 × 10 4 s, so more than 16 h. For the BRUX station, the convergence time to the error level of 10 −13 for both solutions is faster, lasting 9 × 10 2 s. The results for both stations show higher accuracy in the initial analysis period for the GREC solution when compared to G + ISB. Figure 6 shows three coordinate components for AREG and BRUX. For both stations, a more accurate result is obtained in the GREC solution, for which there are far fewer outliers compared to G + ISB. Most of the discontinuities in G + ISB coincide with the day boundaries. The STD of station coordinates for the AREG station is 0.014, 0.017, and 0.055 m for the North (N), East (E), and Up (U) components, respectively in the GREC solution. In the G + ISB solution, the STD value is 0.020, 0.027, and 0.060 m for N, E, and U, respectively. For the BRUX station, the STD value in the GREC solution is 0.011, 0.012, and 0.025 m for N, E, and U, respectively. In the G + ISB solution, it is 0.029, 0.023, and 0.042 m for N, E, and U. To summarize this part, we can conclude that estimating a single clock parameter and ISB values is insufficient for multi-day and multi-GNSS PPP, despite the receiver being connected to one physical clock. This fact is mainly caused by the reference clock changes in the satellite products used for positioning. Superior station coordinate and clock stabilities can be obtained from the solution GREC in which, for each system, different epochwise clock parameters are estimated, despite increasing the Fig. 4 Analysis of the stability of clock parameters for GREC and G + ISB using MDEV for AREG and BRUX stations number of the estimated parameters by a factor of four compared to G + ISB.
The following results show analogous analyses for the other three solutions with different types of clock resets, i.e., GREC + res, G + ISB + resG, and G + ISB + res. Figure 7 shows results for the AREG and BRUX clock estimates in three solutions. The characteristics of the results of this analysis are similar to the results obtained in Fig. 2 for the first two solutions without clock resets; however, differences at the several-centimeter level occur. For the BRUX station for the G system for GREC + res solution and the common clock for the other two solutions, Fig. 7 shows slight differences for different computation strategies, even though the vertical axis of Fig. 7 has a level of several meters. For the HM clock in BRUX, there is much less noise in the solution than for the Rb clock at AREG. The system-specific clock estimates in GREC + res for BRUX are not parallel as they tend to converge at the end of the series, especially for the G and E clocks. Figure 8 shows the ISB results for three solutions with resetting ISBs. The left column shows the results for the GREC + res solution for AREG and BRUX with the identical vertical axis, the middle column shows the results for the G + ISB + res solution, whereas the right column shows the results for the G + ISB + resG solution.
For the AREG station, the STD value in the GREC + res solution is 0.236, 0.197, and 0.577 m for the reconstructed ISBs, i.e., R-G, E-G, and C-G, respectively. For the G + ISB + res solution, the STD of estimated ISBs equals 0.226, 0.415, and 0.599 m for R-G, E-G, and C-G, respectively. For the G + ISB + resG solution, the ISBs are not reinitialized, which results in STD values below 1 mm when excluding the first epochs. For the BRUX station, the STD value of the GREC + res solution is 0.195, 0.193, and 0.193 m for R-G, E-G, and C-G, respectively. For the G + ISB + res solution, the STD value is 0.356, 0.261, and 0.550 m for R-G, E-G, and C-G, respectively. For both stations, smaller variabilities in ISB are achieved in the GREC + res solution than in G + ISB + res despite more independent parameters are estimated in GREC + res. For the G + ISB + resG solution, the results retrieved are similar to the G + ISB solution in Fig. 3. Table 4 shows the STD values for the ISB of each of the strategies. GREC and GREC + res are relatively close to each other, similarly to G + ISB and G + ISB + resG strategies. However, for the G + ISB + res strategy, the particular ISB differs from the results obtained with the GREC and GREC + res strategies. This information is crucial as it indicates the variations in observation noise among different strategies. Despite resetting the ISB at the boundary of the day for both the GREC + res and G + ISB + res strategies, STD value of ISB differs for these strategies. Figure 9 shows the MDEV analysis for the solutions GREC + res, G + ISBr + resG, and G + ISB + res. As in the analysis in Fig. 4, there are no noticeable differences between the solutions for the G clocks. The MDEV stability shows similar characteristics in each solution for BRUX and AREG stations.
Analyzing the MDEV errors in Fig. 10, the most accurate results are obtained for the G + ISB + resG and G + ISB + res solutions. In both cases, the solution is more accurate than for G + ISBrec without considering the reset at the boundary of the day, as shown in Fig. 5. The GREC solution in Fig. 5 is more accurate than the GREC + res  Fig. 10 when comparing the two solutions with system-specific clock estimates. Figure 11 shows the series of estimated station coordinates for all three solutions considering a reset at the day boundary. The most accurate solution for the AREG station is G + ISB + res. The STD value for this solution is 0.013, 0.018, and 0.053 m for N, E, and U, respectively. For BRUX stations, the most accurate solution is GREC + res. The STD value for this solution is 0.011, 0.012, and 0.025 m for N, E, and U, respectively. However, for BRUX stations, the G + ISB + res solution gives similar results: 0.012, 0.013, and 0.025 m for N, E, and U, respectively. The G + ISB + res solution for BRUX stations also has fewer single outliers than GREC + res.
For AREG, the most optimal solution is G + ISB + res out of all solutions with and without resetting clocks parameters. For BRUX, the GREC, GREC + res and G + ISB + res solutions achieve similar results. This shows that for the less stable clock types, the choice of solution is much more important for the final solution than for the more stable clock types, where similar results are obtained by different strategies to determine the clock parameters in a given system. G + ISB + res is also characterized by a reduced number of estimated parameters when compared to GREC and GREC + res. Figure 12 shows the differences between the reference solution G + ISB + res and the other solutions for the G clock or the common clock estimates. Similar results are obtained for the differences between the solutions against GREC and GREC + res and between the solutions against G + ISB and G + ISB + resG. The STD value for AREG stations reaches 0.126, 0.126, 0.206 and 0.216 m for the differences  . 9 Analysis of the stability of clock parameters for GREC + res, G + ISB + resG, and G + ISB + res using MDEV for AREG and BRUX stations Fig. 10 MDEV errors for the clock parameter for the G system for AREG and BRUX stations from GREC + res, G + ISB + resG, and G + ISB + res against GREC, GREC + res, G + ISBrec and G + ISB + resG, respectively. The STD value for the BRUX station reaches 0.073, 0.074, 0.162, and 0.164 m for the differences relative to GREC, GREC + res, G + ISB, and G + ISB + resG, respectively. Figure 12 reveals large differences between the common clock and G clock estimates depending on the number of system-specific clock estimates and clock or ISB resets. The differences reach 0.75 m (or 2.5 ns) with large clock discontinuities for different days. The G + ISB + res solutions are more consistent with GREC and GREC + res than with G + ISB and G + ISB + resG, which means that estimating ISB as constant monthly values is insufficient to obtain stable clock results.

Determination of the coordinate stability for each strategy
To identify the most stable solution, Fig. 13 shows the 3D position repeatability for all stations, including all solutions with the STD and IQR values. STD values assume the normal distribution of data, thus, are more sensitive to outliers than IQR values that include only data between the 75th and 25th percentiles.
The least stable solutions are G + ISB and G + ISB + resG. Similar results for the other three solutions are obtained with similar STDs varying within a few mm. The main advantage of the G + ISB + res solution is the smaller number of clock parameters estimated per day compared to GREC and GREC + res solutions. Comparing the STD and IQR results for each station shows how many outliers there are in these solutions. From the results in Fig. 13 alone, it is impossible to distinguish the type of clock used. For XO, Rb, Cs, and HM, a similar relationship between the solutions can be seen, which means that the observation geometry of the satellites in each system has a greater effect on the quality of the results than the clock at the station. The results for the YEL2 station show much worse stability for STD values than the other stations, despite using HM. Accuracy is also inferior for the NKLG station, which uses the XO clock. Figure 14 shows the FFT results for N, E, and U coordinate components for selected stations considering all five solutions. The G + ISB and G + ISB + resG solutions are the weakest in terms of FFT noise. In both cases, the noise and diurnal signals significantly differ from the other solutions. In most cases, for GREC, GREC + res, and G + ISB + res solutions, the characteristic signals can be determined for 8 h, 12 h, and 24 h, especially for the stations PTBB, REDU and ZIM2 for the E and U components. In the case of the other four stations, the characteristic signals are not so clear, but the 8 h, 12 h and 24 h signals for the individual components are still noticeable. Some stations also show a 6 h signal, such as REDU for the U component. When analyzing the PPP station coordinate results, one must be very careful with interpreting the diurnal signals and their harmonics because their amplitudes may significantly differ when using different strategies to clock handling. For ZIM2, the amplitudes of the 24 h signal of the U component equal 18 and 24 mm in G + ISBrec + res and GREC Fig. 12 Differences between the reference G + ISB + res clock and other clock solutions for the G system for AREG and BRUX stations solutions, respectively, despite the fact that both solutions have been recognized as similar in terms of the station coordinate repeatability. Table 5 shows the median STD values of station coordinates determined based on all stations analyzed for the entire period of analysis. For G + ISB and G + ISB + resG, the results are the worst as observed in previous analyses. Superior accuracies are obtained for the other solutions -GREC, GREC + res, and G + ISB + res. For the 3D median values, G + ISB + res provides superior results and is slightly better than GREC by 0.4 mm. GREC + res performs slightly worse than G + ISB + res by 1.2 mm. The choice of the G + ISB + res solution is justified both by the fewer parameters and by the accuracy of the solution. Figure 15 shows the percent values of coordinate residuals smaller than 1 cm, between 1 and 2 cm, and 2-3 cm, etc., for N, E and U components relative to the mean value for stations AREG, BRUX, CRO1, and TID1 and five strategies. For all coordinate components for each station, the G + ISB + res strategy gives the most accurate results because the number of coordinate residuals is the largest in the range from 0 to 1 cm. The largest discrepancies occur at TID1 for the E component, with a 23% difference between the most and least accurate results for residuals below 1 cm. The BRUX station has the most stable results for each coordinate component. For the other stations, the E coordinate gives inferior results when compared to BRUX. Similar results can be obtained from Fig. 13, where the STD value for the 3D shows larger errors for the AREG, CRO1 and TID1 stations than for the BRUX station. In Fig. 15, the AREG station achieves 84, 67, and 37% of the  . 15 Percentage of coordinate residuals smaller than 1 cm, from 1 to 2 cm, from 2 to 3 cm, etc., with respect to the mean value for the N, E, and U components. Results are shown for AREG, BRUX, CRO1, and TID1 stations for all five solutions residuals for N, E, and U, respectively, in the range between 0 and 1 cm for the best strategy G + ISB + res. The BRUX station achieves 78, 75, and 42%, whereas CRO1 achieves 77, 66, and 28% for N, E, and U, respectively, better than 1 cm. By far the largest number of values in the least accurate range above 10 cm are found for the TID1 and CRO1 stations. However, the number of coordinate residuals above 10 cm is similar for the GREC, GREC + res, and G + ISB + res strategies, which showed similar performance in all analyses. The results shown in Fig. 15 directly indicate the G + ISB + res strategy as the most accurate. G + ISB + res is better than other solutions by at least 1% of the values with coordinate residuals between 0 and 1 cm. Figure 16 illustrates the differences in the time transfer between station PTBB and stations BRUX, TID1, and YEL2 in the form of an MDEV diagram. The reference solution is G + ISB + res for the station PTBB. For YEL2-PTBB, large differences of 5 ×10 −10 between GREC + res, G + ISB, and G + ISB + res solutions are obtained for short averaging intervals, which are especially important for the time transfer. The differences between GREC and G + ISB + res solutions are much smaller of 5 ×10 −11 for short intervals, even when clocks from different systems are used, e.g., E clock from the GREC solution and common clock from the G + ISB + res solution. For long averaging intervals of 4 ×10 3 s, the stability for each solution converges to a similar result for YEL2-PTBB. For TID1-PTBB in the GREC and GREC + res solutions, the results for each estimated clock coincide at the level of 4 ×10 −11 for short intervals. Interestingly, larger differences are obtained for two G clocks and the common clock from GREC and G + ISB + res than between Galileo from GREC and the common clock from G + ISB + res. For averaging intervals of 2 × 10 5 s, almost all results coincide, except for the results obtained for the C clocks. In the case of BRUX-PTBB, the most stable solution is for the G and E for solutions with independent clocks and the common clock for solutions with the ISB for short intervals and the G clock or the common clock for long intervals. For the other systems, the clocks show inferior stability. Unlike TID1-PTBB and YEL2-PTBB, for BRUX-PTBB the results for the G clock or the common clock reach the level of 2 × 10 −14 , while for the other two MDEV differences, the stability is worse by about one order of magnitude. Figure 16 shows that the time transfer results may significantly differ when using a different strategy for handling the receiver clock and ISB, whereas the selection of the GNSS system plays a minor role, especially for GPS and Galileo.

Conclusion
The strategy used to determine the clock parameters in multi-GNSS solutions, or a clock parameter along with the ISB, is a very important choice in order to obtain the most accurate results of station coordinates and time transfer. This article has shown five different strategies to evaluate the advantages and disadvantages of each solution. Two solutions can be considered almost equivalent, i.e., (1) GREC, with epoch-wise estimating of system-specific clocks as random variables and (2) G + ISB + res, with the epoch-wise estimation of a common clock together with ISBs as daily constants. Resetting only the common clock parameter, as in G + ISB + resG is insufficient to obtain superior results. For solutions with inferior observation geometry, GREC provides slightly worse results than G + ISB + res. GREC requires a large number of independent clock parameters to be estimated. Resetting clocks and ISB is mandatory when estimating ISB. For the GREC solution, resetting the clock parameters can be skipped. Neglecting the reset of the ISB leads to large coordinate outliers at day boundaries and inferior clock estimates. The G + ISB + res solution leads to the largest number of station coordinate differences smaller than 1 cm. Median STD values of the 3D station coordinates are 39.0, 39.8, 69.0, 66.3, and 38.6 mm for GREC, GREC + res, G + ISB, G + ISB + resG, and G + ISB + res, respectively. The G + ISB + res solution requires significantly fewer parameters and achieves slightly better stability results than solutions with independently determined clocks. When determining the common clock, it is necessary to reset this parameter at the day boundary due to the changing of tracked satellites and the selection of a different reference clock. This issue is not related to the performance of individual GNSS systems, individual receivers, or the PPP algorithm with ISB estimation.
We tested also another strategy for determining the common clock and the ISB-G + ISB + resISB which involves resetting only the ISB at the day boundary. However, upon analysis, we argue that this strategy is not superior in terms of accuracy with respect to G + ISB + res or GREC.
The time transfer between the PTBB reference station and the BRUX, TID1, and YEL2 stations has been analyzed employing the clock and ISB from the G + ISB + res solutions. For the intercontinental time transfer between YEL2 and PTBB, the clocks based on the GREC solution are more consistent with the common clock from the G + ISB + res solution than for the common clock from G + ISB or GPS clock from GREC + res. The difference in the time transfer can reach even 5 ×10 −10 for 30 s averaging intervals in MDEV analysis. For TID1 and PTBB, large differences occur for the BeiDou system. In the case of the BRUX and PTBB time transfer, the best stability is achieved when using the clock parameters of the GPS system or the common clock for the strategy with ISB. For all stations, Galileo performs slightly better than the GPS clock or the common clock for short-term averaging intervals, whereas the GPS clock or the common clock provides better results over long periods. In all cases, the GREC-based clocks are most similar to the G + ISB + res-based common clock, whereas the common clock from G + ISB + resG, G + ISB, and clocks from GRES + res for YEL2-PTBB lead to larger discrepancies. However, the choice of strategy for clock handling leads to the MDEV differences of several 10 −11 of time transfer and dozens of millimeters in terms of the 3D STD values.
Similar results were obtained for GNSS stations connected to high-quality atomic clocks, such as hydrogen masers, and GNSS receivers employing internal oscillators, which means that the stability of the receiver clock has a minor impact on the multi-GNSS positioning results based on solutions with different number of estimated clock parameters. However, selecting a proper strategy for receiver clock handling will substantially influence the repeatability of the multi-GNSS PPP positioning. Moreover, it directly impacts the time series of station coordinates, e.g., on the amplitude of diurnal and semidiurnal signals of the station height component, which was revealed by the FFT analysis.