Characteristics of the IGS receiver clock performance from multi-GNSS PPP solutions

Global navigation satellite system (GNSS) receivers belonging to the International GNSS Service (IGS) are equipped with different types of clocks, such as internal crystal quartz clocks, rubidium and cesium atomic clocks, as well as hydrogen masers. These clocks are characterized by different phase and frequency accuracies and stabilities, resulting in different systematic clock time series patterns. We analyze the clock offsets between different GNSS systems, provide noise characteristics of the undifferenced and differenced clock parameters, and detect systematic patterns of the clocks. The time series of the receiver clocks are dominated by the diurnal, semidiurnal, and sometimes terdiurnal signals with amplitudes up to several meters. Hydrogen masers provide the highest clock stability, and the lowest is by internal clocks. However, there are also groups of very stable internal clocks that perform similarly to low-performing hydrogen masers and rubidium clocks. The interquartile ranges for epoch-differenced clock parameters fall between 3 and 250 mm for the best hydrogen masers and the worst internal clocks, respectively.


Introduction
The main tasks of global navigation satellite systems (GNSS) include the determination of the receiver position, velocity, and timing. The International GNSS Service (IGS, Johnston et al. 2017) was established with the aim of providing high-quality GNSS orbit, clock, Earth rotation, troposphere, and ionosphere parameters (Beutler et al. 1999) and realizing the GNSS-based international terrestrial reference frames (Rebischung et al. 2016). Most of the IGS stations employ internal temperature compensated crystal clock oscillators (XO). However, some stations were equipped with external clocks, such as rubidium (Rb), cesium (Cs) atomic, or hydrogen maser (HM) clocks, which guarantee high-frequency stability. Receiver atomic clocks can be used for time transfer, and some are used for this purpose (Petit 2021) or realization and distribution of the Coordinated Universal Time (UTC, Wu et al. 2021). Modeling the receiver clock parameters provided by a stable clock improves the GNSS solutions, particularly the up component in the kinematic precise point positioning (PPP) solution (Krawinkel and Schön 2016).
Many IGS analysis centers employ a two-step procedure for deriving final solutions. First, the double-differenced GNSS observations are used for deriving orbits, Earth rotation parameters, and station coordinates, which are then reused to determine satellite and receiver clocks based on undifferenced GNSS observations (Bock et al. 2009;Prange et al. 2020). This pragmatic procedure was established in the 1990s to eliminate the issue of a large number of receiver and satellite clock parameters to be estimated for each observation epoch in the final IGS network solutions. Forming double differences eliminates all issues related to the clock parameters and most issues related to biases. However, it deteriorates the observational geometry, especially the up coordinate component and troposphere delays which are strongly correlated with the receiver clock parameters (Gao and Chen 2004;Weinbach and Schön 2011). The current trends in the IGS rely on common estimation of all parameters with proper handling of satellite and receiver clock and all biases, taking into account correlations between clock and all other parameters (Strasser et al. 2019). The stable atomic clocks connected to the IGS receivers will allow for future clock modeling between observational epochs, further stabilizing the GNSS-derived parameters in the network and PPP solutions. Wang and Rothacher (2013) demonstrated the possibility of the up component stabilization in GPS-based solutions with a stable oscillator and proper clock modeling. The stability and performance of the receiver clocks must thus be first analyzed in terms of their applicability for stochastic and deterministic modeling.
We analyze the receiver clock parameters for selected IGS stations (IGS receiver clocks) derived from GPS, GLO-NASS, and Galileo PPP solutions. The receiver clocks are estimated as the white noise parameter for each GNSS. All calculations are performed using the GNSS-WARP software (Wroclaw Algorithms for Real-time Positioning). The GNSS-WARP software is implemented in Matlab and was originally designed for real-time GNSS processing and deriving receiver coordinates, troposphere delays, and clock corrections. In addition, the software allows data processing in post-processing mode and operating in simulated realtime mode using real-time products recorded with the BKG Ntrip Client (BNC) in ASCII files. The undifferenced and uncombined PPP model is used in the processing with ambiguities estimated as float values. Observations are processed using modified least squares adjustment with the propagation of the covariance matrix, which is similar to a Kalman filter approach (Hadaś 2015;Hadas et al. 2019). In this study, we derive a solution using final satellite orbits and clocks from the Center for Orbit Determination in Europe (CODE) and the European Space Agency (ESA) analysis centers.
IGS stations equipped with different receiver types, different antennas, and connected to different atomic clocks are used for the clock analysis. All selected stations belong to the Multi-GNSS Experiment (MGEX) pilot project (Kouba 2009;Montenbruck et al. 2017) and ensure tracking of three GNSS systems to analyze system-specific clock errors. Zajdel et al. (2020Zajdel et al. ( , 2022 identified system-specific errors emerging from different aliasing, revolution periods of satellites, and multipath effects in multi-GNSS (GPS, GLONASS, Galileo) solutions. This raises the question of whether similar systematic errors can be identified in the multi-GNSS clock parameters. There are very few studies identifying systematic patterns and periodicities in GNSS clock time series; Specifically, those derived from MGEX receiver clock parameters and biases empirically derived from GNSS data. The empirical receiver GNSS clock performance depends not only on the quality of the clock and its environment but also on the quality of GNSS orbit and clock products, geometry of observations, the number of tracked satellites, quality of phase and code GNSS observations, local multipath, and receiver performance. Therefore, the stability of the estimated receiver clock parameters may remarkably differ from the stability solely provided by the external clock connected to the GNSS receiver.
This study aims to identify and analyze systematic errors and provide noise characteristics of the system-specific clock parameters for IGS stations connected to different types of internal and external clocks. The various atomic clocks are compared with the identification of misbehaving atomic clocks at selected IGS stations. We provide the inter-epoch stability characteristics of GNSS-derived clock corrections. These analyses will make it possible to introduce stochastic modeling superimposed on the clock parameter, which will allow for GNSS coordinates and tropospheric delays of higher accuracy than from a solution in which the clock is estimated as white noise, and the inter-epoch clock correlations are neglected.

Methodology
Clock stability analysis for receivers is conducted based on 27 selected MGEX stations that employ different receiver and antenna types and are connected to different clocks. Each selected station tracks at least three systems: GPS (G), GLONASS (R), and Galileo (E). Figure 1 shows the selected, globally distributed, multi-GNSS stations connected to different types of internal and external atomic clocks denoted using different colors. Figure 1 also includes information about the manufacturer of the receiver. Additional information about the analyzed stations is given in Table 1. Technical data provide the status for the period for which the analyses are performed.
Throughout the selection of multi-GNSS stations, different manufacturers were considered: Septentrio, Leica, Trimble, Javad, and stations equipped with different clock types: XO, Rb, Cs, and HM. Observations from 4.05.2020 to 11.05.2020 (DoY 125-131, 2020) were used for the analysis of a representative network of stations. The seven-day data allow us to verify the clock jumps at the 24-h boundary for solutions based on CODE and ESA analysis centers. Analyzed stability, periodicities, and systematic patterns for individual clocks are similar over a longer period. Most of the expected signals in clocks are related to the diurnal signals, satellite revolution periods, and harmonics of these periods. Also, the multipath effect repeats every ground track repetition period. The data were processed in the GNSS-WARP software in a static PPP mode employing the strategy described in Table 2. The clock parameters are estimated as epoch-specific parameters independent for each observation epoch and each GNSS system. An uncombined, undifferenced PPP model is used in the computation. In contrast, the wet tropospheric delay is estimated with a constraint between epochs of 4mm √ h , where h stands for one hour.
The station coordinates are common for all GNSS systems and are derived as static parameters. This study aims to evaluate the clock parameter stability and identify systematic periods for potential stochastic clock modeling to eventually improve GNSS positioning. To be comparable with the GNSS positioning results, clock results are expressed in meters, which can easily be expressed in time units by dividing by the speed of light.

Results
Various statistical analyses are conducted for each constellation to characterize clock stabilities statistically. These include standard deviations (STD) of estimated clock values and the differentiation between neighboring epochs, fast Fourier transform (FFT), deriving the amplitudes of diurnal, semidiurnal, and terdiurnal signals, as well as modified Allan deviation (MDEV) and Hadamard deviation (HDEV) analysis. Inter-system biases based on ESA and CODE products are also designated to verify the stability of delays between systems. Most of the results have been referred to the PTBB stations to remove the effect of the reference clock changes in the reference CODE and ESA products. Thus, PTBB equipped with a stable HM was adopted as the reference station.  are referred to the PTBB clock. The internal clock from WROC shows the worst performance. Its FFT reveals characteristic signals with dominating diurnal and semidiurnal signals with amplitudes of about 2.4 and 3.1 m. For DAE2, the FFT includes diurnal and terdiurnal characteristic signals with amplitudes reaching 0.8 m. For BRUX, the FFT shows no characteristic periodic signals. All GNSS systems provide similar clock error repeatabilities for each station  shown, which means that the clock errors must be associated with local environmental conditions, such as temperature, electromagnetic stabilities, GNSS cabling, and other common GNSS errors rather than the limitations and repeatabilities of the GNSS constellations because GPS, GLO-NASS, and Galileo are characterized by different revolution periods. GLONASS clocks show sporadic jumps that are related to insufficient observation geometry. The clock offsets for GLONASS are larger than those between GPS and Galileo, corresponding to substantial inter-system biases for GLONASS.

Raw clock analysis and FFT
We compare solutions based on raw CODE and ESA products to separate product-specific errors from solutionspecific errors using BRUX clock parameter. Figure 5 at the top shows larger discontinuities at daily boundaries for ESA than for CODE products, which can be related to the fact that CODE generates 3-day overlapping solutions (Lutz et al. 2016). For the middle day of the analysis, the CODE products led to a positive drift of the clock parameter, whereas ESA products led to a negative drift (Fig. 5), which can be related to the different selection of the reference clocks in both analysis centers' products affected by different linear   drifts. GPS, GLONASS, and Galileo have similar clock offsets corresponding to the inter-system biases. The lower part of Fig. 5 shows values that are referred to the same reference clock at the PTBB station. In this case, the daytime boundary jumps are almost entirely reduced. However, some jumps are introduced by the reference clock into the GLONASS clock parameter at BRUX. The clock drift and offsets are similar when based on ESA and CODE products, which means that the clock differentiation efficiently removes the differences in the reference products provided by different IGS analysis centers, which are mainly caused by various selections of a reference clock. Again, GPS, GLONASS, and Galileo have similar clock offsets corresponding to the inter-system biases.
The mean clock offset at the 24-h boundary for the raw ESA products is 171, 169, 181 mm for GPS, GLONASS, and Galileo, respectively (see Fig. 5, top, left). In contrast, the mean offset is equal to 23, 95, 23 mm for the raw CODE products for GPS, GLONASS, and Galileo. The results show better stabilities for the solutions based on CODE products, mainly for GPS and Galileo. However, after deriving clock parameter differences with respect to the PTBB clock, no more jumps at day boundaries can be identified in either CODE and ESA products.
The different noise characteristics for different clocks can be found when analyzing differenced clocks against the clock at the PTBB station from different stations in Fig. 6. An interesting phenomenon is a distinctive positive offset for GLONASS in relation to the other systems for the majority of stations, especially for stations equipped with Trimble NETR9 receivers.
The best high-frequency clock stability is obtained for the HM clocks when comparing the clocks shown in Fig. 6, which is shown by less noise and data spikes. However, the initial offsets for Rb clocks are larger than for XO and Cs clocks. ZIM2 and BOR1 are located in Switzerland (Zimmerwald) and Poland (Borowiec), respectively, but are equipped with different types of clocks and the same type of receiver. For ZIM2 and BOR1, the clock estimates with diurnal and semidiurnal variations are almost the same, which local electromagnetic conditions cannot explain. The similarity of both estimated clocks suggests that BOR1 has been incorrectly described in a log file and relies on the internal clock instead of the external HM clock. ZIM2 and BOR1 share the same receiver type and a similar observation geometry above the stations, as they both are located in Central Europe. However, the local multipath conditions are different for both stations. Therefore, we infer that the observation geometry and receiver type may dominate clock estimates for the stations equipped with internal clocks.
Although WROC, ULAB, and SEYG are equipped with internal oscillators, the clock quality for the latter station is much better than that for the first two. None of the analyzed stations shows signals related to the satellite revolution periods of different GNSS constellations. Instead, the observed signals are always harmonics of a synodic day for all satellite systems. WROC and ULAB are equipped with Leica and Javad receivers, respectively, whereas SEYG is equipped with the same Trimble receiver as employed in ZIM2. Thus, distinct differences may occur for clock parameters even when using the same receiver and clock pair (Fig. 6).
The Fourier analysis of GNSS clocks shows distinct signals of 12 h and 24 h for most of the stations. The amplitudes of diurnal and semidiurnal signals are shown in Fig. 7 after referring the values to PTBB. For each station, similar values for each system, including GPS, GLONASS, and Galileo, are recorded. From the comparison of the clock types, the HM clocks seem to provide the smallest amplitude values for the 12 h and 24 h signals, except for BOR1 and HOB2. For the remaining stations with HM clocks, the mean Fig. 6 Values for individual clocks at selected stations using ESA products for May 2020, as the day of the month. The results are referred to the PTBB clock amplitudes of the diurnal signals with reference to PTBB are 114, 140, and 94 mm for GPS, GLONASS, and Galileo, respectively. For semidiurnal signals, the mean amplitudes equal 53, 63, and 23 mm for GPS, GLONASS, and Galileo, respectively. Consequently, for HM clocks, due to their high stability, temperature causes less destabilization of the clock than for XO, Rb, and Cs.

Inter-epoch clock differences for adjacent epochs of observations
The clock inter-epoch differences are analyzed between the adjacent epochs. It means that for the clock time series, we computed the difference between the one epoch and the previous one with the observation interval of 30 s. Figure 8 shows results of the clock differences of adjacent epochs for two stations to verify the short-term clock stability. The short-term stability is important in the context of the proper stochastic modeling of the clock parameter in static and kinematic PPP.
The time series in Fig. 8 shows a large variability of the difference clock values between the two stations. Additionally, the results are similar for GPS, GLONASS, and Galileo systems. The STD of raw and epoch-differenced clock readings are given in Table 3 as G, R, E and dG, dR, and dE, respectively. Table 3 shows that the IGS stations with time laboratories, such as BRUX and PTBB, provide the highest time stability standards. For station BRUX, the best receiver clock parameter obtained for Galileo provides an STD of 466 and 4 mm for raw and differenced clocks, respectively, whereas the STD of the ZIM2 clock parameter equals 1676 and 32 mm, respectively. Table 3 also shows STD values for inter-system biases for several other stations. The clock at station PTBB provides similar values as the clock at station BRUX, which implies similar clock stability. Figure 9 shows the stability of the inter-system bias for clocks at six stations. The results are referred to the PTBB clock. Interestingly, it is impossible to distinguish which station is connected to the atomic clock and which one relies on the internal clock solely based on the analysis of inter-system clock differences. The inter-system bias is stable in the analyzed period between GLONASS and Galileo clocks with respect to GPS clocks for almost all stations. Despite the fact that the raw ZIM2 clocks are affected by high-frequency noise, the inter-system biases are not affected by high-frequency scatter, which implies that the high-frequency variations affect the system-specific clocks in the same manner. We typically see only the inter-system delays (biases) due to antenna and receiver electronics, e.g., for ZIM2 and BOR1 equipped with the same receiver type, similar biases are obtained. Thus, inter-system analysis alone makes it more difficult to verify the clock type. However, Defraigne and Verhasselt (2018) identified an improved time transfer connection based on a multi-GNSS solution with the inter-system bias handling.
Interquartile range (IQR, the spread difference between the 75th and 25th percentiles of the data) statistics are performed to verify the quality of the estimated clock parameters for the particular system. IQR has been selected due to its robustness to gross errors as opposed to STD. Two statistical analyses are shown in Fig. 10; the results for the raw clocks obtained for each epoch are shown at the top, while the bottom panel depicts the results of the differencing between neighboring epochs. The stations are sorted by clock type. The epoch-undifferenced data with respect to the clock parameter of the PTBB station in the top part of  Although the stability of raw clock readings for some HM is worse even than for XO clocks due to substantial drifts, i.e., the frequency offsets, the difference of the epochs for HM shows a substantial advantage of HM in terms of frequency stability. In contrast to HM, XO clocks are less stable in consecutive epochs. However, they are relatively stable in the long term, which is reflected in minor drift and small IQR values for raw clock values. Figure 11 shows the MDEV and HDEV analysis for seven selected stations referred to the PTBB clock based on the Galileo clock. HDEV is insensitive to clock drifts, as opposed to MDEV. The line style distinguishes the clock type for each station. The best short-term stability is obtained for the BRUX clock, for which the MDEV value is equal to 1.9 × 10 −11 for 30 s integration intervals. MDEV values show similarities between ZIM2 and BOR1 despite the different clock types included. The clock at station WROC is characterized by significantly inferior shortterm stability, which is two orders of magnitude worse than for the best-evaluated stations. For long-term stability, the best results are achieved by BRUX. The MDEV of this station for the 12 h integration time is equal to 0.4 × 10 −12 . This shows that the best HM clocks are characterized by both high short-term and long-term stabilities. Cs and XO clocks provide a similar level of long-term stability in the range from 1.0 × 10 −12 to 7.0 × 10 −12 for 12 h integration time. However, for short-term stability, which is much more important in the context of receiver clock modeling in PPP, Cs and XO clocks are inhomogeneous and indicate significant stability differences. HDEV shows a similar performance of different clocks to the MDEV analysis. Additional analysis of MDEV and HDEV is conducted for each GNSS for three stations equipped with stable HM clocks. Figure 12 illustrates MDEV for the clocks at stations BRUX and YEL2 grouped by GNSS relative to PTBB clocks. The most stable receiver clocks in the short-term domain are obtained for Galileo and also for GPS, for which MDEV are equal to about 1.5 × 10 −11 and 2.0 × 10 −11 for 30 s integration time. Interestingly, for the PTBB-YEL2 short-term stability, the Galileo results are slightly better than those for GPS clocks. GLONASS's initial stability is inhomogeneous and slightly worse than that of GPS and Galileo. MDEV for the integration time equal to 12 h is similar for GPS, GLONASS, and Galileo and is close to 4 × 10 −12 .

MDEV and HDEV analysis
Interestingly, the short-term stability between YEL2 and PTBB is slightly better than that between BRUX and PTBB, even though YEL2 and PTBB are located on different continents, as opposed to BRUX and PTBB. BRUX typically tracks the same set of GNSS satellites due to the proximity of PTBB; therefore, most of the satellite orbit and clock errors cancel out when forming differences of PPP-derived clock values. The high accordance between YEL2 and PTBB results shows that the intercontinental time transfer is currently of superior quality thanks to the high quality of satellite orbits and clocks. The accuracy of time transfer depends on the type and quality of the receiver clock, whereas the orbital products do not limit the time transfer accuracy when  using post-processed satellite orbits and clocks. HDEV in Fig. 12 shows a better performance for Galileo in low-Tau values and BRUX-based GPS in high-Tau values. The low-Tau HDEV points characterize the GNSS time transfer noise, whereas the large Tau shows the intrinsic qualities of the clocks as steered. This suggests that the time transfer BRUX-PTBB and YEL2-PTBB should be more effective when based on Galileo than when based on GPS. However, the high-Tau values show a superior solution for GPS BRUX-PTBB, which may correspond to the long-term clock synchronization between two observatories.

Classification of GNSS receiver clocks
After previous analyses, a classification of GNSS stations into groups is conducted based on the characteristics and stability of the clocks. Due to the similar clock values obtained for the three systems, GPS, GLONASS, and Galileo, the further division relies on GPS. Figure 13 provides a classification into five groups with information about the type of clock and the receiver manufacturer.
The division is based on epoch-differenced clock values. Group 1 contains the stations with the least stable clocks, all of which employ XO clocks. The IQR value for WIND is equal to 206 mm, which constitutes a good representative of Group 1. Group 2 includes stations employing Cs clocks. The IQR value for DAE2 is equal to 67 mm. Group 3 comprises all four types of clocks, which, interestingly, are characterized by a similar clock performance. Group 3 contains the four most stable XO clocks, one Rb, two inferior Cs clocks, and two least stable HM clocks. Group 4 includes the most stable clocks from the whole set, namely HM. For Group 4, the IQR values are at the level of a few mm; e.g., for BRUX, IQR equals 2 mm. It should also be emphasized that the stations in Group 4 contain Septentrio receivers, except for GOLD equipped with Javad. Two HM clocks at stations BOR1 and HOB2 could not be classified into Group 4, as they indicate much worse clock stability when compared to other HM clocks at stations. Again, ZIM2 and BOR1 show similar performance; thus, they are both classified into Group 3. Group 5 contains the most stable Cs clocks. The IQR value for KOUG is equal to 9 mm, which is a typical representative of Group 5. Figure 14 shows the cumulative absolute clock stability values for epoch-differenced clocks, i.e., with stability better than 1 cm, from 1 to 2 cm, from 2 to 3 cm, etc. The stations are arranged according to their group assignment in Fig. 13. For Group 1 with Int, the typical stability of the clock is worse than 100 mm. The typical values for Groups 3 and 4 fall between 30-40 and 20-30 mm, respectively. By far, the most favorable results are for stations in Groups 4 and 5 equipped with the best-performing HM and Cs clocks, respectively. In Group 4, almost all results are better than 10 mm. One has to bear in mind that five Cs clocks are classified into Groups 2 and 5, whereas two HM clocks are in Group 3, which means that the performance of Cs and HM is inhomogeneous and only the proper handling of atomic clocks leads to the highest clock stability. Figure 15 depicts the MDEV and HDEV clock stabilities for stations representing each of the groups with respect to the PTBB clock. Group 1 for short-term stability provides the worst results at the level of 1 × 10 −9 in the 30 s integration intervals. The most stable Group 4 provides the shortterm stability of the order of 1-2 × 10 −11 . For the integration Fig. 13 GNSS stations are classified into five groups by analyzing the GPS-based receiver clocks. The different colors indicate a particular group, type of clock, and manufacturer of the receiver employed at the station intervals of about 1 h, the clocks belonging to the Groups 1, 2, and 3 cannot be clearly distinguished, which indicates a low performance of Cs clocks for these intervals. Moreover, Group 5, with the best Cs clocks, seems to be inhomogeneous and affected by color noise, despite very high stability for the short integration intervals. Long-term stability is highest for clocks belonging to Group 4. An MDEV value as low as 3-30 × 10 −12 is reached for the 12 h integration interval for Group 4 stations, which infers a superior quality of the most stable HM when compared to other types of clocks whose stability is one order of magnitude worse. Despite the insensitivity of HDEV to the clock drifts, the overall picture of short-term and long-term clock stabilities is the same for all groups based on MDEV and HDEV analyses.

Effect of temperature stability
The observed clock stability is affected by the ambient temperature and the attempt to maintain a constant temperature in the receiver environment. Unfortunately, this information is incomplete or not specified in the log files for most of the GNSS stations. For example, according to the log files, the temperature is kept as follows: ASCG (25 °C ± 5 °C), BRUX (21 °C ± 0.2 °C), DLF1 (20 °C ± 3 °C), KOUG (21 °C ± 3 °C), SEYG (25 °C ± 5 °C), WROC (22 °C ± 2 °C). This information directly shows that the greatest temperature stability is for BRUX, where it only varies by 0.2 °C. At the other stations, the temperature may vary by several degrees. The proper temperature stability may also be relevant for clocks.

Conclusion
We evaluated the stability of the clock parameters from PPP solutions for GNSS stations connected to internal clocks and different types of external clocks. No evident differences in the clock stability are found between GPS-and Galileobased clock parameters. GLONASS receiver clock parameters are affected by sporadic jumps. No system-specific periodic signals related to the constellation revolution period or aliasing effects of the repeatability of the satellite and station positions could be found in the time series of clock parameters, which implies that the selected GNSS system has a minor impact on the derived clock values. However, most of the analyzed clocks show distinct diurnal, semidiurnal, and terdiurnal signals with amplitudes even up to 5 m.  HM clocks provide the smallest periodic variations in clock readings; however, some of the HM show substantial secular drifts which is related to high short-term time stability but low long-term time stability or frequency offsets. Secular drifts of HM do not constitute a problem when differencing clocks between adjacent epochs or when introducing relative weighting for the clock parameter in PPP solutions.
Clocks can be classified into different groups based on their characteristics. However, this classification cannot rely solely on the information provided in the log station files because even clocks of the same type connected to the same receiver may lead to different clock characteristics. Therefore, each IGS station requires an individual analysis of the clock stability. In our classification, all possible clocks, HM, Cs, Rb, and XO could be put to Group 3 when using the epoch-differenced clock values, which means that the best-performing internal clocks may have comparable performance to the low-performing HMs. This situation directly shows that the quality of clock stability is affected by the type of clock itself but also by other factors, such as environmental conditions. Thus, the characteristics of the clock are composed of more variables that condition the clock to maintain a given stability. The most stable HM receiver clock parameters were obtained for receivers produced by Septentrio and one by Javad. The stations BRUX, PTBB, and YEL2 used clock models CH1-75A, VCH 1003 M, and VCH-1008, respectively, during the data acquisition period.
HMs are the best option for GNSS receivers from the MDEV analysis and inter-epoch differencing. The IQR values for epoch-differenced clock parameters fall between 3 and 250 mm for the best HM and the worst internal clocks, respectively, which shows the large discrepancy between the stability of clocks that equals two orders of magnitude.
We show that not only the type of the station clock matters but also the GNSS constellation and geometry of observations, which is sometimes inferior for GLONASS. When using different reference products from CODE and ESA, distinct jumps are detectable in the clock time series. However, most of the reference product differences disappear when the clock difference between clocks is formed. The quality of intercontinental time transfer is the same as the intra-continental, which means that the quality of reference satellite orbits and clocks does not limit the time transfer accuracy. We conducted a series of studies of clock stability at ground stations based on various quality indicators-FFT, MDEV, HDEV, inter-epoch, inter-system, and inter-clock differencing, which is new in the literature. Results are presented as a prelude to further research of stochastic clock modeling for improving the PPP position products, which requires a comprehensive characteristic of the receiver clock stability. There are significant differences in estimated receiver clocks when based on raw ESA and CODE products; however, these differences disappear as soon as the reference clock is unified and changed to PTBB. We provided a comparison between GNSS systems to find any system-specific errors in clock parameters related, e.g., to satellite revolution periods or multipath effect and found that the satellite ground track repeatabilities for different GNSS cannot be identified in the clock parameters. Instead, the quality of receiver clock estimates is dictated by the quality of the receiver clock and the GNSS system with the superior performance of Galileo and slightly worse of GPS. Therefore, the analysis of the stability of clocks is essential for proper stochastic clock modeling to improve the station coordinate stability and tropospheric delays in future studies.

Data availability
The GNSS observations used for this study emerge from the MGEX campaign (Montenbruck et al. 2017) of the IGS (Johnston et al. 2017) and can be freely downloaded from the Crustal Dynamics Data Information System, https:// cddis. nasa. gov Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.