Observed impact of the GNSS clock data rate on radio occultation bending angles for Sentinel-6A and COSMIC-2

Radio occultation (RO) measurements require the tracking of signals from a global navigation satellite system (GNSS) by a low-Earth-orbit (LEO) satellite as the signals travel through different layers of the atmosphere. The orbit and clock solutions for the GNSS constellations affect these measurements in two ways: they are needed to obtain a zero-differencing GNSS-based orbit and clock solution for the LEO, and they directly enter the processing of each single RO profile, where the orbit and clock information for the transmitter (GNSS) and receiver (LEO) satellites is required. In this work, we investigate how different GLONASS and GPS orbit and clock solutions affect the statistical properties of RO profiles by comparing our results with forward-modeled bending-angle profiles based on data from the European Centre for Medium-Range Weather Forecasts (ECMWF) short-range forecasts. Given that GNSS orbits are relatively smooth, this study focused on the effect of different transmitter clock data rates, and we tested the range from 1 to 30 s (specifically, the rates of 1, 2, 5, 10 and 30 s). The analysis is based on the reprocessing of Sentinel-6A data (4 months in 2021, September to December, or about 110 000 occultations) and of a smaller sample of recent COSMIC-2/FORMOSAT-7 data (5–7 August 2023, or about 9000 occultations). We find that at impact heights above about 35 km, GLONASS bending-angle statistics markedly improve with the use of higher-rate clock information. For GPS, the statistics are better for more recent GPS blocks, and a rate of 5 s provides a marginal improvement over the 30 s rate for all blocks. In the same impact height range, higher-rate GLONASS clocks also significantly reduce the vertical error correlation. These results are likely the manifestation of the different kinds of short-timescale behavior of the atomic clocks on board the GPS and GLONASS constellations.

1 Introduction

The calculation of bending-angle (BA) profiles based on signals from a global navigation satellite system (GNSS) traveling between a GNSS satellite and a low-Earth-orbit (LEO) spacecraft through the Earth's atmosphere requires accurate knowledge of the position, velocity and clock of the two space vehicles. In the radio occultation (RO) processing described in this work, the GNSS positions and clocks are provided as auxiliary data. They are used explicitly in the computation of each BA (through, e.g., the Doppler equation; Kursinski et al., 1997) and implicitly, given that they are required to obtain the LEO orbit and receiver clock solutions.

The GNSS space vehicles have revolution periods in excess of 10 h, and their orbits, which are largely controlled by gravitational dynamics, are relatively smooth and are typically provided at a rate of few minutes. The GNSS clocks are needed to synchronize the receiver clock when performing the LEO precise orbit determination (POD). Due to the random stochastic noises that affect GNSS atomic clocks, a smaller sampling interval is required to obtain accurate interpolations, as needed by the bending-angle-retrieval process. As a reference, the Center for Orbit Determination Europe (CODE) currently provides orbits and clocks with a 5 min and a 30 s sampling rate, respectively, for rapid products (latency of 18 h; Dach et al., 2023b) and 5 min and both 5 and 30 s sampling rate for orbits and clocks for final products (latency of 2 weeks; Dach et al., 2023a). The Jet Propulsion Laboratory (JPL) Global Differential GPS (GDGPS) real-time (RT) GNSS products are available with a 1 min orbit sampling and a 1 s clock sampling. Table 1 provides the relevant information on these and other sets of GNSS products that have been used in this study.

RO represents a high-rate application of GNSS products, since a typical radio occultation event has a duration of few minutes at most. At these time intervals, different GNSS constellations and different clock hardware within a given constellation show a range of clock stabilities (e.g., Hauschild et al., 2013; Griggs et al., 2015). The implication is that a proper description of the behavior of a clock (i.e., to have a model of its evolution) will require a product with a rate that depends on the clock properties. In other words, more unstable clock hardware benefits from the use of a higher-rate clock correction.

The analysis of Sentinel-6A (S6A) RO data performed at the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) using different GNSS data streams provided the motivation for this study, as described in Sect. 2. The experiment setup is described in Sect. 3. The bulk of the analysis and the results in Sect. 4 are based on a 4-month-long batch of data from the S6A RO instrument (Sect. 4.1) and are complemented by a smaller batch of COSMIC-2/FORMOSAT-7 data (Sect. 4.2). The discussion and conclusions are presented in Sects. 5 and 6, respectively.

2 Motivation

The S6A spacecraft is the latest member of the family of altimetry reference missions (e.g., Donlon et al., 2021). For the first time, such a mission is equipped with an RO instrument, which was built by JPL (see also Paolella et al., 2024, this special issue). The receiver is connected to a POD antenna, which tracks only GPS signals, and to two occultation antennas, which track both GPS and GLONASS signals, located in the velocity and anti-velocity direction. The entire mission is operated by EUMETSAT. For the RO instrument, JPL is responsible for the provision of near-real-time (NRT) RO data, while EUMETSAT is responsible for non-time-critical (NTC) RO products. These two operational streams are based on the same level-zero data, received by EUMETSAT from the spacecraft through its ground segment and redistributed to JPL. JPL also hosts one of the analysis centers of the international GNSS service (IGS) and, through an operational agreement, provides the auxiliary GNSS data for both RO data streams. The JPL GPS auxiliary products provided to EUMETSAT for the NTC processing come from the JPL final stream, with orbits and clocks at 15 m and 30 s, respectively. The JPL GLONASS products come from the JPL RT stream but with orbits at 15 m and clocks at 1 s (we refer to this set, which mixes RT and final products, as JPL-EUM; see Table 1). The NTC data stream (from now on, we refer to this stream as OPE) has a formal timeliness of 60 d, but it is kick-started as soon as the GNSS auxiliary data are available. Typically, EUMETSAT receives these data within 2 weeks of the sensing time.

Figure 1 Sentinel-6A robust statistics against ECMWF-based short-range forecasts for two data sets for the period 24 February–31 March 2024, corresponding to about 38 000 single bending-angle profiles for each set. (a) Operational data stream (OPE), where the POD is based on Bernese 5.2 (Sect. 3.1) and the GNSS auxiliary data are from JPL, with a 30 and 1 s clock rate for GPS and GLONASS, respectively. (b) STC offline data stream, where the POD is based on Bernese 5.4 (Sect. 3.1) and the GNSS auxiliary data are CODE rapid products, with clocks at 30 s for both GPS and GLONASS.

Table 1 Summary of the different sets of GNSS products that have been used for this study.

¹ The rate of the CODE orbit products was increased from 15 m to 5 m after the adoption of the new IGS20 reference frame (starting on 27 November 2022). ² For the year YYYY, products are found at http://ftp.aiub.unibe.ch/CODE/YYYY/ (last access: 9 July 2025); ³ http://ftp.aiub.unibe.ch/CODE/ (last access: 9 July 2025). ⁴ https://www.gdgps.net (last access: 9 July 2025). ⁵ https://cddis.nasa.gov/archive/gnss/products/realtime (last access: 9 July 2025), https://sideshow.jpl.nasa.gov/pub/JPL_GNSS_Products/RealTime (last access: 9 July 2025). ⁶ Since these are a downsampled version of the JPL RT stream, no latency can be defined for them. ⁷ These products are not publicly available but are part of an operational agreement between JPL and EUMETSAT.

In addition to the OPE processor for monitoring purposes at EUMETSAT, we deployed an offline short-time-critical processor (referred to as the STC processor), which uses the CODE rapid products (COR) as GNSS auxiliary data (Dach et al., 2023b) and has a timeliness of about 24 h. As part of this monitoring activity, for both the STC and the OPE processors, robust statistics are computed for the bias and standard deviation (Maronna et al., 2019) with respect to the forward-modeled bending-angle profiles obtained from the temperature, pressure and humidity profiles extracted from the European Centre for Medium-Range Weather Forecasts (ECMWF) short-range forecasts. The forward modeling is done with the Radio Occultation Processing Package (ROPP; Culverwell et al., 2015).

A comparison of the statistics between the STC and OPE data streams provides a first indication of the importance of the GNSS clock data rate for the quality of the derived RO products. In Fig. 1 we plot the bias and the standard deviation for both the OPE and the STC processor, using recent S6A data from the period between 24 February and 31 March 2024, which corresponds to about 38 500 individual bending-angle profiles per processor. For this period, the total number of products obtained with the two processors is slightly different at about the 3 % level. The differences can be ascribed to the effect of the different POD software, the different GNSS auxiliary data used, and the possibility of a different number of level-zero data due to occasional time delays in the satellite's data transmission (such data are received through so-called archive dumps), which could be missed by the STC processor. The bias profiles between the two processors and for the two constellations are self-consistent and hardly distinguishable. For the standard deviation, the situation is different. Both for STC and for OPE the GPS clock products have a 30 s sampling rate. The different data source (COR for the STC and JPL-EUM for the OPE; see Table 1) does not affect the overall standard deviation for GPS. For GLONASS, the STC processor uses the COR 30 s clocks, while the OPE processor is fed with 1 s JPL-EUM clock products. Above an impact height of about 35 km, the GLONASS standard-deviation curves (green) are markedly different.

Previous work has pointed out the lesser short-term stability of GLONASS clocks with respect to GPS ones, at least starting with Block II-F satellites (Griggs et al., 2015; Hauschild et al., 2013). The decreased stability impacts the BA uncertainty, as shown in a preliminary analysis of a small batch of COSMIC-2 data by Yao et al. (2023). Less stable clocks would require a higher-rate correction to compensate for their noise, and the results of Fig. 1 point in this direction.

Figure 2 Modified Allan deviation for the GPS (a) and GLONASS (b) constellations for 31 December 2021, calculated using JPL RT clock data (Table 1). Please note the different ranges of the y axes.

Given that the work of Hauschild et al. (2013) and Griggs et al. (2015) is several years old, and both the GPS and the GLONASS constellations have been evolving in the meantime, we show in Fig. 2 the modified Allan deviation (e.g., Griggs et al., 2015) for the 1 s JPL RT clock products for 31 December 2021, a date which is part of the main S6A processing campaign below (Sect. 4.1). At that time in the GPS constellation there were 15, 12 and 5 Block II-R, Block II-F and Block III satellites, respectively. The results for the GPS satellites are similar to those of Fig. 7 of Griggs et al. (2015), with the addition of the group of Block III satellites, which show a better performance on the short timescales (<10 s) and are similar to results of Block II-F between 10¹ and 10² s. The GLONASS constellation was almost completely composed of GLONASS-M satellites, with a single GLONASS-K1 satellite, which displays a slightly better performance at the very short timescales.

3 Setup of the experiments

Unless otherwise noted, in what follows all the experiments described are based on the reprocessing of radio occultation data from the S6A and COSMIC-2/FORMOSAT-7 missions, from level zero up to BA. These two missions feature the same TriG instrument, a hardware component that builds on the former BlackJack/IGOR receivers and is a POD and RO receiver able to track signals from multiple GNSS constellations (see Esterhuizen et al., 2009, for additional details), equipped with an ultra-stable oscillator (Tien et al., 2012). The level-zero data are decoded using JPL-provided software compliant with the International Traffic in Arms Regulations (ITAR). The RO processing is performed with Yet Another Radio Occultation Software (YAROS) developed and maintained by EUMETSAT, and detailed information on the processing approach can be found in Paolella et al. (2024, this special issue). Here, we provide some details on the POD processing setup and on the GNSS auxiliary data employed.

3.1 POD

The LEO orbit and clock solutions are obtained with the latest version of the Bernese GNSS software (version 5.4) using zero-differenced GNSS data (Jäggi et al., 2007; Bock et al., 2011). Bernese 5.4 allows for the modeling of non-gravitational forces on the LEO, namely air drag and direct and reflected radiation pressure (Mao et al., 2019). To make use of these models, a representation of the shape of the spacecraft and of the material properties of its exterior layers is required. In our simulations we use a macro-model for S6A based on Montenbruck et al. (2021), but no shape model is employed for COSMIC-2. Bernese 5.4 has the capability of performing zero-differencing integer-ambiguity resolution using the Bias-SINEX products (Villiger et al., 2019; Schaer et al., 2021). However, these products are not part of the operational agreement with JPL, and our solutions are based on the float-ambiguity resolution approach. Integer-ambiguity fixing has been shown to perform better than the float-ambiguity approach (e.g., Montenbruck et al., 2018). However, the use of the float-ambiguity solution in the BA processor is more than sufficient for generating high-quality RO products (e.g., Kursinski et al., 1997; Montenbruck et al., 2008; Innerkofler et al., 2020). The POD software embedded in the operational processor in 2021 was Bernese version 5.2. The OPE plot in Fig. 1 is the only figure based on data processed with Bernese version 5.2. This version has neither the non-gravitational force models nor the zero-differencing integer-ambiguity-resolution algorithms available. However, as we argue below in this section, differences in the POD solution (due to, e.g., the different models employed by the POD software) do not affect the BA statistics.

Figure 3 Comparison along the radial (a), cross-track (b) and along-track (c) directions of two S6A orbit solutions for the month of November 2021. In each panel the mean ± the root mean square (RMS) is indicated. The two solutions are obtained with the same software (Bernese 5.4) but with two different GNSS GPS auxiliary data products (Table 1): JPL 30 s (15 min orbits and 30 s clocks) and CODE 5 s (15 min orbits and 5 s clocks). The date format is year/month/day.

We assess the quality of a given POD solution in several complementary ways. There is an internal quality assessment within Bernese, where the final so-called reduced-dynamic orbit is compared to the kinematic orbit solution. A kinematic orbit is obtained from the geometrical solution for the position of the LEO at the epochs in which observations are available (e.g., Švehla and Rothacher, 2005). A reduced-dynamic orbit indicates an orbit solution where the dynamical forces originating from the extended gravitational field of the Earth, from the Sun, etc. are taken into account. The model cannot realistically take into account all possible forces, hence “reduced”, and some empirical parameters compensate for this (see, e.g., Montenbruck et al., 2005, and references therein.). Bernese handles the data used to generate these two orbits in the same exact way (same data rejection parameters, same data weighting approach, etc.). As a result, the comparison of the reduced-dynamic orbit to kinematic orbit cannot detect any problem originating from, e.g., an error in the attitude modeling or in the definition of the position of the POD antenna, since these kinds of errors would affect the two orbits in the same exact way. Thus, we also use different POD software, the NAvigation Package for Earth Observation Satellites (NAPEOS; Springer et al., 2011), to cross-check the consistency of the Bernese-based reduced-dynamic solution with the NAPEOS-based solution. As a result of being an altimetry reference mission, in addition to the RO POD antenna (RO-POD), S6A features another POD antenna (GNSS-POD), which tracks both GPS and Galileo and is connected to a completely independent receiver with its own ultra-stable oscillator (e.g., Donlon et al., 2021). The comparison of the GNSS-POD solution with the RO-POD solution provides a receiver-independent way to cross-check the solution (of course the clock solution in this case cannot be compared, given that the two receivers have different oscillators). Furthermore, the GNSS-POD-based solution itself is routinely compared with solutions obtained by the members of the Copernicus POD quality working group (CPOD QWG) (Fernández et al., 2024), which EUMETSAT is part of. The above checks allow us to perform several complementary software-, receiver- and processing-center-independent cross-comparisons to validate our solution.

Figure 4 Comparison of different S6A POD solutions in the 3D-RMS space. The reference solution, against which the 3D RMS is computed, is the CPOD QWG combined solution (see main text for details). The pink and cyan solutions are the same as shown in Fig. 3 and are based on data of the RO-POD receiver, which tracks only GPS. The other three solutions are obtained by processing data of the GNSS-POD receiver, which tracks both GPS and Galileo. AIUB: Astronomical Institute of the University of Bern (their solution for 19 November was an outlier and has been removed); DLR: Deutsches Zentrum für Luft- und Raumfahrt; JPL: Jet Propulsion Laboratory. The date format is year/month/day.

The focus of this work is to investigate the effect of different GNSS clock data rates on the RO BAs (see Fig. 1), and the effect is both direct, since they enter explicitly in the processing of each single RO profile, and indirect, through their influence on the POD solution. In Fig. 3 we compare two S6A solutions for the month of November 2021 (which falls within our main test campaign; see Sect. 4.1). One is obtained using JPL RT GPS products, downsampling clocks to 30 s and orbits to 15 min, and the other using CODE final products with clocks at 5 s and orbits at 15 min (sets JPL 30 s and COD 05 s in Table 1). The RINEX observation files, decoded from level-zero data provided by the RO receiver at the rate of 1 s, are downsampled to 30 s when using JPL products and to 10 s when using CODE products. The use of 10 s RINEX data is the customary choice in POD reprocessing campaigns (see the discussion related to the regular service review in Fernández et al., 2024). However, in the POD processing, using a RINEX data rate higher than the clock data rate of the GNSS products would not be beneficial but would rather increase the noise of the solution (Dach et al., 2015). Thus, for JPL 30 s, we use RINEX data downsampled to 30 s. It is clear from Fig. 3 that there is good agreement between the two solutions, with biases below the millimeter level and RMS below 2 cm in any given direction. The comparison for the velocity solution (not shown) has both biases and RMS below 0.01 mm s⁻¹. These results are well within the standard target specifications of 5 cm and 0.05 mm s⁻¹ for RO missions (e.g., Innerkofler et al., 2020).

Figure 5 Robust statistics of S6A GLONASS occultations against ECMWF-based short-range forecasts, for the period September to December 2021. Each set corresponds to bending-angle profiles processed with different GLONASS clock data rates, as indicated in the legend (see also Table 1). Above an impact height of about 35 km, there is a clear advantage in using high-rate GLONASS clocks.

As an additional example of the consistency of the POD results, we plot in Fig. 4 the three-dimensional root mean square (3D RMS) for a set of solutions, which includes the two cases of Fig. 3, based on data recorded by the RO-POD receiver, along with three solutions, based on data recorded by the GNSS-POD receiver, obtained by three members of the CPOD QWG. In Fig. 4, we use the CPOD QWG combined solution as the reference solution (against which the 3D RMS is computed), which is a weighted average of the solutions provided by all the members of the CPOD QWG (see Sect. 3.6.2 in  Fernández et al., 2024, for details). Typically, the spread of the different orbits in the 3D-RMS space is well below 3 cm. In Fig. 4 the only exception is the RO-POD solution obtained using JPL 30 s GNSS products (Table 1), which is expected given that it is the only set of data based on a real-time (and thus less accurate) GNSS data stream.

Figure 6 Robust statistics of S6A GPS occultations against ECMWF-based short-range forecasts, for the period September to December 2021. Each set corresponds to bending-angle profiles processed with different GPS clock data rates, as indicated in the legend (see also Table 1). Unlike the case of S6A GLONASS occultations (Fig. 5), the standard deviation is largely unaffected by the GPS clock data rate used, as seen in (b), which is a zoomed-in view of the 50–60 km impact height range.

Figure 7 Above 50 km robust standard deviation of S6A GPS occultations against ECMWF-based short-range forecasts (as in Fig. 6, right panel), for different GPS clock data rates (legend) and GPS transmitter types. (a) Block II-R; (b) Block II-F; (c) Block III. From left to right the standard deviation improves at any given clock data rate. However, the effect of increasing the GPS clocks' data rate is not the same for each block (i.e., panel).

3.2 GNSS auxiliary data

For the results presented in the following sections, the GNSS auxiliary data that we employ are JPL RT data (Table 1). The original orbit products are downsampled to 15 min. Different sets of GNSS auxiliary data are created by modifying the clock products. We use a straightforward downsampling (i.e., decimation) of the clocks to create four additional sets of clock products with rates of 2, 5, 10 and 30 s (JPL dd s products, Table 1). Another obvious approach would be to fit the data of the 1 s product to obtain a smoothed version before downsampling it. We verified that these products outliers are not a concern. Furthermore, given that the choices that need to be made when fitting (function to use, degree of the polynomial in a polynomial fit, etc.) may have an effect on the results, we proceeded with the straightforward downsampling.

Figure 8 Robust statistics of COSMIC-2 GLONASS occultations against ECMWF-based short-range forecasts, for the period 5–7 August 2023. Each set corresponds to bending-angle profiles processed with different GLONASS clock data rates, as indicated in the legend. Above an impact height of about 35 km, there is a clear advantage in using high-rate GLONASS clocks, as was the case for S6A (Fig. 5).

Figure 9 Robust statistics of COSMIC-2 GPS occultations against ECMWF-based short-range forecasts, for the period 5–7 August 2023. The curves have been obtained using 30 s GPS clock products. Each set corresponds to a different GPS block, as indicated in the legend. As was the case for S6A for a given clock data rate (Fig. 7), Block III is better than Block II-F, which is in turn better than Block II-R.

Figure 10 Variation in the L1 phase obtained with different clock rates P_𝚌𝚛 (legend) with respect to the phase obtained using 1 s clocks, P_𝟶𝟷𝚜. The panels refer to two GLONASS (a, b) and three GPS (c, d, e) setting occultations. (Note that the vertical axis range is different between the top and the bottom row.)

4 Results

Given the self-consistency of POD solutions obtained using different GNSS products and approaches (Figs. 3 and 4), the focus of this section is on the effects of using GNSS clock products with different data rates in the BA data processing. The GNSS transmitter clock bias plays an important role when processing the reconstructed signal phases (i.e., the ones obtained by the level-zero decoder) with a zero-differencing algorithm (Beyerle et al., 2005). The main test uses S6A data, and a validation test is performed using COSMIC-2 data.

For all the experiments presented in this section, the POD solution for the LEO has been obtained using the GPS products of the JPL 30 s data set (see Table 1) with clocks at 30 s and, correspondingly, RINEX data at 30 s (Sect. 3.1). For the RO processing, each GNSS auxiliary data set has the same orbits and Earth rotation parameters. Five different sets are used, corresponding to differences in the data rate of the clock products, namely 1, 2, 5, 10 and 30 s (JPL dd s products, Table 1), as discussed in Sect. 3.2.

4.1 Sentinel 6A

We processed 4 months of S6A data, from 1 September to 31 December 2021, which corresponds to around 110 000 occultation profiles. The results of the bias and standard deviation against ECMWF-data-based forward-modeled profiles are presented per constellation.

For GLONASS, Fig. 5 shows that while the bias is unaffected by the use of different GLONASS clock data rates (in the 1–30 s range), the standard deviation markedly improves for impact heights above about 35 km. At 40 km the standard deviation decreases by about 1 % going from a 30 s to a 1 s clock product. Using 30 s GLONASS clocks, a 5 % standard deviation is reached at an impact height of about 45 km, while with a 1 s clock, this occurs at about 52 km.

Figure 6 shows that, as is the case for the GLONASS occultations, the robust bias for the GPS occultations is practically unaffected by the different GPS clock data rates employed. Unlike GLONASS, though, the standard deviation only shows a minor dependence on the GPS clock data rate. However, as is visible in the standard deviation at high impact heights between 50 and 60 km (Fig. 6, right panel), there is not a clear trend in the standard deviation with the increasing clock data rate, contrary to the case of GLONASS (Fig. 5). Given that the clock behavior of the GPS constellation shows a block dependence (Fig. 2), we also plot the bias and standard-deviation curves at high impact heights for the different GPS blocks separately (Fig. 7). Overall, for a given clock data rate, the standard deviation improves going from Block II-R to Block II-F to Block III.

However, for a given block, the effect of increasing the GPS clock data rate does not affect the standard-deviation curve consistently. For Block III (Fig. 7, right), going from 30 to 1 s clocks brings a steady increase in the standard deviation at a given height, with the effect being really minor between 30 and 5 s but becoming noticeable with 2 and 1 s clocks. For Block II-F (Fig. 7, middle), all sets are mostly coincident, with a slight degradation of the standard deviation when using the 1 s clock. For Block II-R (Fig. 7, left), increasing the clock data rate improves the standard deviation. The improvement is steady in going from 30 to 5 s. The 2 and 5 s curves coincide, while going to a rate of 1 s increases the standard deviation.

4.2 COSMIC-2

To verify the results observed using S6A data, we used data from the entire COSMIC-2 fleet (e.g., Ho et al., 2020) in a smaller-scale test using 3 d in 2023 (5–7 August), which corresponds to about 9200 occultations. These results are based on processing of level-zero data available on the UCAR server (UCAR COSMIC Program, 2019).

As was the case for S6A, the increase in the clock data rate for GLONASS occultations improves the standard deviation with minor effects on the bias (Fig. 8). For GPS occultations we only present results for GPS clocks at 30 s (Fig. 9), and these results for the standard deviation are consistent with what was found with S6A, where there is an improvement in going from Block II-R to Block II-F to Block III (Fig. 9). Results for higher GPS clock rates are consistent with those obtained with S6A. The standard deviation for COSMIC-2 is worse than for S6A, despite having the same receiver. We attribute this difference to a POD solution of lesser quality: no macro-model has been used for COSMIC-2, the satellites orbit at a lower altitude (about 550 km versus 1350 km for S6A) where atmospheric perturbations are larger, and the solar activity increased in 2023 with respect to 2021. Future work will further investigate possible improvements to the COSMIC-2 statistics through improved POD configurations.

5 Discussion

The analysis performed in this work points to the importance of using higher-rate GLONASS clock data to obtain the best bending-angle statistical performances in terms of standard deviation at impact heights in excess of 35 km, with improvements observed up to 1 s, which is the highest rate tested (Figs. 1, 5 and 8). There seems to be no simple dependence of the standard-deviation curves for the GPS constellation as a function of the clock data rate (Fig. 7), and the 30 s clock products provide comparable performance to the GLONASS 1 s products (Figs. 1, 6 and 9). These results confirm expectations based on the analysis of GNSS clock noise at RO-relevant timescales (Hauschild et al., 2013; Griggs et al., 2015).

Figure 11 For S6A GLONASS (a) and GPS (b) occultations, the vertical error correlation, shown here as the full width at half maximum (FWHM) evaluated from the spread around the diagonal of the vertical-correlation matrix at each level. The data set used for this figure is the same S6A data set used in Figs. 5 and 6.

The above S6A results are based on global statistics of about 110 000 occultations. However, it is informative to also look at single occultations to understand how the GNSS clock sampling rate affects the low-level processing of the measurements. We focus here on the zero-difference phase P_𝚌𝚛 (linear units), which we define as the NCO (numerically controlled oscillator) phase P_𝙽𝙲𝙾 (cycles) corrected for the receiver (τ_𝚛𝚎𝚌) and transmitter (τ_{𝚝𝚛𝚊𝚗}) clock bias:

$$\begin{array}{}\text{(1)}& P_{\mathtt{\text{cr}}}=c\left[{\displaystyle \frac{P_{\mathtt{\text{NCO}}}}{f}}-\left({\mathit{\tau }}_{\mathtt{\text{rec}}}-{\mathit{\tau }}_{\mathtt{\text{tran}}}\left(\mathtt{\text{cr}}\right)\right)\right].\end{array}$$

In the equation, f and c are the frequency of the signal and the speed of light, respectively. Given that in our experiments the LEO POD solution is always the same (and thus τ_𝚛𝚎𝚌 does not change), P_𝚌𝚛 can be used to assess the effect of the transmitter clock bias error in the measurements, since the bias error does depend on the clock rate cr as indicated with the notation τ_{𝚝𝚛𝚊𝚗}(cr). The quantity τ_{𝚝𝚛𝚊𝚗}(cr) is the value of the transmitter clock bias linearly interpolated to the measurement times.

Figure 12 Modified Allan deviation for 31 December 2023 for (a) the Galileo constellation, calculated using JPL RT Galileo clock data (available from the websites listed in note 5 of Table 1), and (b) the BeiDou constellation, calculated using CODE MGEX data (available at http://ftp.aiub.unibe.ch/CODE_MGEX/CODE/, last access: 9 July 2025).

We calculate the difference in P_𝚌𝚛, for cr equal to 2, 5, 10 and 30 s, with respect to P_01 𝚜, which we regard as the reference. For the frequency L1, in Fig. 10 we plot the variation in the corrected phase P_𝚌𝚛 for each constellation and for each piece of hardware. We note the following:

• 1.

  GLONASS P_𝚌𝚛 variations (top row) are about a factor-of-3 larger than the GPS P_𝚌𝚛 variations (bottom row).

• 2.

  For GLONASS, increasing the clock rate significantly reduces the range of these differences.

• 3.

  For GPS, the reduction in the variations with increasing clock rate is minor, with the possible exception of the 30 s curve for Block II-R-A.

Similar features characterize the frequency L2 (not shown). Since the bending angle is proportional to the time variation in the phase (e.g., Kursinski et al., 1997), it is reasonable to conclude that what we see in this single-occultation-level analysis maps to the trend (or lack thereof) of the standard-deviation curves in Figs. 5, 6, 7, 8 and 9. It is also consistent with the better Allan deviation shown by the GPS satellites (Fig. 2).

We also investigated the impact of the different GNSS clock rates on the vertical error correlation. The higher the vertical correlation, the lower the information content of a given profile (Bowler, 2020), with implications for the assimilation of these data in numerical weather prediction (NWP) models (e.g., Gilpin et al., 2019; Bowler, 2020). The results are shown in Fig. 11. With behavior that is by now familiar, increasing the clock rate for GLONASS leads to a marked improvement in the vertical error correlation, while for GPS the variation is minimal. Once again, it makes sense that sampling the noisier GLONASS clocks with a higher rate leads to a less-correlated set of bending-angle profiles. Finally, it is interesting to note that while there is no obvious trend in the standard deviation across the GPS blocks with respect to different clock rates (Fig. 7), in the vertical-correlation plot, the higher the clock rate, the better the vertical correlation, even though the improvement is limited.

The observed connection of the Allan deviation (i.e., the stability of the given clock hardware) with the standard deviation associated with the robust statistics of the BA profiles (Figs. 5 and 6), the phase variations in single occultations (Fig. 10) and the vertical error correlation (Fig. 11) points at the possibility of using the Allan deviation to guide our expectation for future analysis of Galileo and BeiDou data sets. In Fig. 12 (left) we show the modified Allan deviation for the 1 s JPL RT clock products for the Galileo constellation on 31 December 2023. As with GPS (Fig. 2), the older satellites (in-orbit validation, IOV) have a higher Allan deviation than the newer satellites (full operational capability, FOC). Over the whole interval range shown in Fig. 12 (left), the values are comparable to those of the GPS constellation in Fig. 2. The JPL RT stream only provides GPS, GLONASS and Galileo products. For BeiDou, we used the CODE data from the Multi-GNSS EXtension (MGEX; Prange et al., 2020), whose clock products have a rate of 30 s, and the Allan deviation can only be computed for averaging intervals of 1 min or longer (Fig. 12, right). In this limited region the values for the newer BeiDou-3 satellites are below those of GPS and Galileo. The curves in Figs. 2 and 12 are in line with recent estimates of the Allan deviation of GNSS clock hardware (Jaduszliwer and Camparo, 2021). Based on these results, and under the assumption that the BeiDou curves extend smoothly at short averaging intervals, we expect that using 30 s Galileo and BeiDou products will provide results similar to those we reported for GPS products in this study.

6 Conclusions

This work focused on the impact of using different GPS and GLONASS clock rates on the occultations recorded by the S6A and COSMIC-2 RO receivers. We evaluated the statistical properties of BA profiles obtained using GNSS clock products with different rates by comparisons with forward-modeled BA profiles based on ECMWF short-range forecasts. The use of higher-rate clock products has almost no impact on the bias statistics, but it has an effect on the standard deviation above about 35 km. For GPS occultations, there is a marginal improvement in going from 30 to 5 s, and overall, more recent blocks show better standard deviation. For GLONASS, the improvement is very pronounced and continues for clock rates up to 1 s, which is the highest rate we tested (Figs. 5 and 6). Similarly to the standard deviation, the vertical error correlation for GLONASS occultations markedly improves with the use of high-rate clock products, while the improvement is limited for GPS occultations (Fig. 11). These results are also reflected in the refractivity statistics derived from the Radio Occultation Meteorology Satellite Application Facility (Fig. 3.17 in  Syndergaard and Lauritsen, 2021).

The results we obtained are a manifestation of the stability of the GNSS clock hardware that can be characterized through the Allan deviation, which is overall better for the GPS constellation with respect to the GLONASS constellation (Fig. 2). In the near future, the RO instrument on board S6A and the one in the forthcoming follow-up mission, S6B (launch expected from late 2025 onwards) will start collecting Galileo signals on the RO antennas, in addition to GPS and GLONASS (EUMETSAT Public Document EUM/LEO-JASCS/DOC/21/1253471, 2024). This upgrade will provide a new data set that can be used to extend the analysis performed here to include the Galileo constellation and test our expectations based on the Allan deviation (Fig. 12). In addition, modern RO missions also exploit Galileo and BeiDou signals, as exemplified by the large data set of the Radio Occultation Modeling Experiment (ROMEX; McHugh et al., 2023; Anthes et al., 2024).

For RO observations, the best combination of accuracy and small uncertainty is expected in the 5–30 km range (Kursinski et al., 1997), also referred to as the “RO sweet spot”. Our results are consistent with this expectation in the accuracy (i.e., bias), but in terms of uncertainty (i.e., standard deviation), they show an increase below 10 km, as also reported in other studies (e.g., Anthes et al., 2022, and references therein). Recently, RO observations have been included as key observations in the assessment of temperature trends in the upper troposphere and lower stratosphere (Masson-Delmotte et al., 2021, Figure TS. 10, p. 70). Given that GLONASS higher-rate clock data clearly improve the uncertainty and reduce the vertical error correlation, RO processing centers should aim at using this type of product, given its potential impact on both NWPs (e.g., Lonitz et al., 2021) and climate studies (e.g., Gleisner et al., 2022).