Precise Orbit Determination for Climate Applications of GNSS Radio Occultation including Uncertainty Estimation

Global Navigation Satellite System (GNSS) Radio Occultation (RO) is a highly valuable remote sensing technique for probing the Earth’s atmosphere, due to its global coverage, high accuracy, long-term stability, and essentially all-weather capability. In order to ensure the highest quality of essential climate variables (ECVs), derived from GNSS signal tracking by RO satellites in low Earth orbit (LEO), the orbit positions and velocities of the GNSS transmitter and LEO receiver satellites need to be determined with high and proven accuracy and reliability. Wegener Center’s new Reference Occultation Processing System (rOPS) hence aims to integrate uncertainty estimation at all stages of the processing. Here we present a novel setup for precise orbit determination (POD) within the rOPS, which routinely and in parallel performs the LEO POD with the two independent software packages Bernese GNSS software (v5.2) and NAPEOS (v3.3.1), employing two different GNSS orbit data products. This POD setup enables mutual consistency checks of the calculated orbit solutions and is used for position and velocity uncertainty estimation, including estimated systematic and random uncertainties. For LEOs enabling laser tracking we involve position uncertainty estimates from satellite laser ranging. Furthermore, we intercompare the LEO orbit solutions with solutions from other leading orbit processing centers for cross-validation. We carefully analyze multi-month, multi-satellite POD result statistics and find a strong overall consistency of estimates within LEO orbit uncertainty target specifications of 5 cm in position and 0.05 mm/s in velocity for the CHAMP, GRACE-A, and Metop-A/B missions. In 92% of the days investigated over two representative 3-month periods (July to September in 2008 and 2013) these POD uncertainty targets, which enable highly accurate climate-quality RO processing, are satisfied. The moderately higher uncertainty estimates found for the remaining 8% of days (∼5–15 cm) result in increased uncertainties of RO-retrieved ECVs. This allows identification of RO profiles of somewhat reduced quality, a potential benefit for adequate further use in climate monitoring and research.


Introduction
Precise orbit determination (POD) has become indispensable for many space-borne applications that monitor the Earth's climate system, one of which is Global Navigation Satellite System (GNSS) [1,2] radio occultation (RO) [3][4][5]. The RO remote sensing technique is considered highly valuable for atmosphere and climate sciences since it provides high-vertical-resolution measurements over the troposphere and stratosphere with global coverage, high accuracy, long-term stability, and virtually all-weather capability [6,7]. During an occultation measurement GNSS signals scan the atmosphere in limb sounding geometry and arrive with a time delay at the receiving RO satellite in low Earth orbit (LEO), which is due to the signal's refraction in the Earth's atmosphere. A vertically rising or setting occultation event is observed depending on whether the GNSS transmitter satellite rises or sets behind the Earth's horizon from the viewpoint of the rapidly moving LEO receiver satellite.
The atmospheric excess phase path, which denotes the difference between the geometric straight-line distance between transmitter and receiver satellite and the signal's bended propagation path, can be derived from the measured time delays and reliable POD results. This traceability to fundamental time standards ensures long-term stability and no need for calibration [8][9][10]. The short-term stability of an individual RO event of about 1 to 2 min duration is satisfied by highly stable clock oscillators. Finally, essential climate variables (ECVs) [11], such as temperature, pressure, and tropospheric water vapor, can be derived from the raw measurements utilizing a dedicated RO retrieval. Resulting, the GNSS RO space-geodetic observing system delivers an abundance of geographically distributed vertical profiles of ECVs all over the globe, which can be averaged and used for climate studies.
Due to the unique properties of occultation soundings introduced above and global availability since 2001, records of basic RO measurements (i.e., atmospheric excess phase or derived Doppler shift) have the potential to serve as Fundamental Climate Data Record (FCDR): A well-characterized, global, and long-term stable data record for the derivation of accurate and stable ECVs, globally and covering timespans from days to decades [12,13]. However, in order to fully exploit this potential and provide climate benchmark data for meteorology, climate research, climate monitoring as well as for calibration and validation, the accuracy of RO data needs to be assessed. That includes quantification of remaining uncertainties throughout the entire retrieval, starting from raw measurement data with the modeling of the observation geometry, to derivation of atmospheric excess phase data, to the final thermodynamic ECVs (Figure 1).
While the quality and high accuracy of RO in the upper-troposphere and lower-stratosphere regions is well-acknowledged (e.g., [14][15][16]) and high consistency has been assessed between different RO retrievals of leading international processing centers [10,17,18], a rigorous uncertainty estimation and propagation throughout the entire RO retrieval remains an important but incomplete task. In addition, the former occultation processing system at WEGC [16] lacks this utility since the retrieval starts from external excess phase data and thus does not tie to the raw measurements and the physical unit of time. The new Reference Occultation Processing System (rOPS) [19,20] developed at WEGC, on the contrary, aims to establish such a fully traceable processing comprising all retrieval steps [19,21,22]. For this reason, the rOPS features the integration of RO low-level data processing starting from the raw satellite observation data. This low-level processing comprises the RO observation geometry modeling within the daily system modeling (DSM) and the level 1a (L1a) excess phase processing [23] as part of the occultation data processing (ODP) chain (cf. Figure 1). The DSM provides an advanced setup for POD of the RO receiver satellites in LEO, including the capability to routinely assess the uncertainties of the computed orbit data for climate-quality processing, which is presented in this paper.
In order to ensure highly consistent and accurate RO-derived ECVs tied to the raw measurements, precise orbit positions, velocities, and clock estimates of the GNSS transmitter satellites and LEO RO receiver satellites need to be determined and their attributed uncertainty assessed. In this study we focus on the assessment of orbit position and velocity, whereas the satellites clock characteristics and stability will be discussed in a separate publication focusing on RO excess phase processing.
Numerous past studies have evaluated the relationship between LEO orbit accuracy and atmospheric parameters derived from RO [3,[24][25][26][27][28][29]. Based on the orbit accuracy specifications deduced from these studies, the POD processing presented in this study aims for a daily LEO orbit accuracy within 5 cm in position and 0.05 mm/s in velocity respectively, enabling highest quality RO retrieval results. In order to test these orbit requirements we routinely calculate different orbit solutions for mutual consistency check employing two independent POD software packages Bernese GNSS software v5.2 [30] (in short "Bernese" hereafter) and NAPEOS v3.3.1 [31], and use GNSS orbit and clock data from different orbit data archives, the Center for Orbit Determination in Europe (CODE) [32] and the International GNSS Service (IGS) [33].  Reference Occultation Processing System (rOPS) relevant to this study, comprising the input data for the precise orbit determination (Input Data; orange box), the daily system modeling (DSM; green box), and the occultation data processing (ODP; red box). For spelling of institution acronyms (in Input Data and DSM boxes) and further description see Section 3. The ODP box outlines the main level 1 and level 2 retrieval steps (L1a to L2b) from excess phase to essential climate variables (ECVs).
The computation of estimated random and systematic orbit uncertainties is based on orbit inter-comparison between the different solutions, the analysis of satellite laser ranging (SLR) residuals, and random error propagation within Bernese. Finally, days exceeding the orbit accuracy target specifications are associated with increased uncertainty estimates and flagged for the subsequent occultation data processing. Building (also) on these POD uncertainty estimates, the uncertainty propagation along the RO level 1 and level 2 data processing chain (ODP in Figure 1) ultimately delivers thermodynamic ECVs jointly with uncertainty estimates. This processing leads to reference data products enabling added-value for climate monitoring and applications from the co-estimated uncertainties.
Following this introduction, Section 2 starts with an overview of the processed missions. Section 3 introduces the data and software used followed by the description of POD and SLR technique, and the uncertainty estimation approach. Results are presented and discussed in Section 4. In Section 5 a summary and conclusion is provided.

Missions and Spacecraft Payload
The life-time of the RO satellite missions selected for the present study of POD for climate applications span the entire RO era, starting with first continuous long-term RO measurements obtained by CHAllenging Minisatellite Payload (CHAMP) in February 2001. The Gravity And Climate Experiment (GRACE) mission provided RO measurements since 2006, although it has been in space since 2002 and previously activated the RO receiver for testing [34]. Meteorological operational satellite A (Metop-A), the first of a series of three meteorological satellites developed in cooperation by the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) and the European Space Agency (ESA), switched on its RO instrument 8 days after its launch on 19 October, 2006, followed by Metop-B on 17 December, 2012, before the series was completed with the launch of Metop-C from Guiana Space Centre, Kourou on 7 November, 2018 [35].
The CHAMP satellite was launched on July 15, 2000 into a near-polar but non sun-synchronous orbit of 87.3 • inclination and a mean initial orbit altitude of 454 km which decreased over time despite intermediate orbit raise maneuvers, to 325 km approaching the decay. The German Research Centre for Geosciences (GFZ) managed mission provided measurements for atmospheric, magnetic and gravity field research before its end-of-lifetime after over 10 years in space [36]. The onboard BlackJack Global Positioning System (GPS) receiver [37], manufactured by the National Aeronautics and Space Administration's (NASA) Jet Propulsion Laboratory (JPL), is connected to a multi-antenna system including a zenith-looking choke ring antenna for navigation tracking of inter-satellite links at a rate of 0.1 Hz, an anti-velocity facing antenna array for occultation tracking at 50 Hz, and a nadir looking antenna for altimetry measurements.
The receiver front-end and electronics support dual-frequency tracking in 16 × 3 channels for C/A, P1, and P2 code respectively. Out of all channels, 12 are allocated for POD (effectively only a maximum of 10 GPS satellites were tracked simultaneously [37,38]) and the remaining 4 channels are reserved for operation in occultation mode for atmospheric sounding [27,39] or alternatively for altimetry measurements. The satellite's orientation relative to the inertial reference frame is recorded in quaternions (1 Hz) by two dual-headed Advanced Star Camera (ASC) assemblies developed by the Technical University of Denmark (DTU) [40,41]. Also onboard, a Laser Retroreflector (LRR), consisting of four cube corner prisms designed by GFZ [42,43], is mounted on the platform. The LRR passively reflects short laser pulses from an emitting ground station on Earth, thus allowing for the calculation of satellite laser ranges [44].
The GRACE twin satellite mission was a joint mission between NASA and the German Aerospace Center (DLR) and launched into similar near-circular orbits of 89.5 • inclination and about 500 km altitude on 17 March, 2002. The two identical satellites were separated approximately 220 km in nominal orbit for the primary mission goal of measuring climate-relevant variations of the Earth's gravity field [45,46]. As a follow-on mission to CHAMP, GRACE features a similar design and payload, with an identical ASC and LRR mounted on the platform. Supplementary, the K-Band ranging system enables measuring the separation changes between the two satellites to micrometer precision in order to map the Earth's gravity field [47]. Equipped with a similar BlackJack GPS receiver as CHAMP [48], signals received with the zenith antenna are used for POD and the aft-looking antenna provides RO measurements for sensing physical properties of the atmosphere [34,49].
As secondary mission goal RO instrumentation was tested on both GRACE satellites for shorter periods before continuous activation of RO measurements on 22 May, 2006 on GRACE-A [27]. The majority of measurements were obtained by GRACE-A intermitted by shorter periods of occultations by GRACE-B (Jul-Dec 2014, Jun-Oct 2015, Apr-Sep 2016) when swapping maneuvers took place, making GRACE-B the trailing satellite [16]. Therefore, we limit the evaluations in this paper to the quality assessment of GRACE-A, seen as representative for both flight models. After decommissioning and atmospheric reentry of GRACE on 24 December, 2017 (GRACE-B) and 10 March, 2018 (GRACE-A) the success of the mission was continued with the launch of GRACE follow-on in May 2018 [50,51]. However, as of early 2020, RO data from this follow-on mission are not yet available to the community.
As part of the EUMETSAT Polar System the Meteorological Operational (Metop) satellite series consists of three flight models Metop-A/B/C, which were placed sequentially in time (see above) in a sun-synchronous polar orbit of about 98.7 • inclination and an altitude between 796 and 884 km [52]. The centerpiece for GPS inter-satellite tracking is the onboard GNSS Receiver for Atmospheric Sounding (GRAS) developed by Saab Ericsson Space [53,54]. The GRAS unit provides dual-frequency navigation and occultation tracking at 12 × 3 channels for L1 C/A and L1/L2 P(Y). A zenith-looking antenna observes GPS satellite signals at 8 channels with a sampling rate of 1 Hz for POD [28,29]. Furthermore, 4 channels are shared by two high-gain beam forming antennas looking in flight velocity and anti-velocity direction for rising and setting occultations respectively, supporting closed loop (50 Hz) and open loop (1 kHz) measurements [55].
For orientation of the satellite a nominal alignment is applied for the POD. In contrast to CHAMP and GRACE, the Metop satellites are not equipped with an LRR for orbit verification from the terrestrial laser tracking network. Up to this day, the Metop satellite series serves with an almost constant number of occultations per day and thereby underscores its importance as a reliable long-term backbone mission for numerical weather prediction and climate research [56].

Data and Tools
The processing setup presented in this study depends on various data sources (see also Figure 1): The GNSS transmitter satellite orbits are acquired from analysis centers (ACs) of the well-established International GNSS Service IGS, whereas data for the GPS-based POD of the receiver satellites in LEO are provided by the respective mission operators and scientific RO groups. SLR data are obtained from the International Laser Ranging Service (ILRS, https://ilrs.cddis.eosdis.nasa.gov/data_and_products/ data_centers) [57,58], as a service of the International Association of Geodesy.
The GNSS orbit, clock, and Earth orientation parameter (EOP) are taken from the CODE reprocessing product series 2015 (repro2015, https://dx.doi.org/10.7892/boris.75876.3) [59] and for mutual cross-check of the derived LEO orbits from the IGS second reprocessing (repro2, http://acc.igs. org/reprocess2.html) [60]. The IGS multi-AC product combines weighted solutions from 9 contributing ACs, one of which is CODE. Although the combined IGS solutions comply with an orbit accuracy of 2.5 cm in satellite position and 75 ps in satellite clock (http://www.igs.org/products) [61,62], their implementation for RO-oriented POD is limited through the lack of provision of clock products at 30 s sampling within repro2 [63].
The CODE reprocessing series repro2015, on the contrary, offers clock estimates with a sampling of 5 s starting January 1, 2003, in addition to 30 s clocks [64]. Another substantial difference to previous versions of the CODE reprocessing is the first use of the extended Empirical CODE Orbit Model (ECOM2) within repro2015 [65]. The present 5 min clock sampling of the combined IGS solution induces a deterioration of the clock interpolation to high rate occultation data at 50 Hz sampling rate and thus also in the subsequent RO processing. Therefore the JPL solution (another highly weighted orbit product within IGS repro2 besides CODE) is currently selected, for cross-check, as a second GNSS orbit input dataset. This will be revised with the availability of the upcoming third reprocessing of the IGS (http://acc.igs.org/repro3/repro3.html).
Both data centers provide long-term stable data products until February 2015 with consistent alignment to ITRF2008. With regard to the impending POD reprocessing effort of the continuous RO era at WEGC starting from 2001 up to recently, both data sets are prolonged with final products from CODE and IGS, which are aligned to ITRF2014.
GPS receivers onboard of RO satellites record navigation tracking data from the zenith antenna for orbit processing and occultation tracking data from the limb-looking atmospheric sounding antenna(s) (Section 2). For POD with Bernese and NAPEOS the binary data from the receiver are either already provided in, or need to be converted to Receiver Independent Exchange Format (RINEX) [66]. Additionally, information in form of either measured or simulated Euler angles or quaternions provides knowledge on the satellite's attitude. The GPS observation data from the LEO receiver are acquired from COSMIC Data Analysis and Archive Center CDAAC (https://cdaac-www.cosmic.ucar.edu/ cdaac) for CHAMP and GRACE, JPL (https://podaac.jpl.nasa.gov) for GRACE, and EUMETSAT (https://eoportal.eumetsat.int) for Metop. The attitude data are acquired from CDAAC (CHAMP: as measured and provided by GFZ; Metop: nominal type) and JPL (GRACE).
The conversion from binary to RINEX format might include a pre-processing and screening of the data, as it is the case for GRACE. In the current setup two different sources for GRACE observation data are exploited: (1) CDAAC: L1A raw data which show inferred clock biases reaching up to 0.5 µs before reset (native to the GRACE clock oscillator [34,49]); (2) JPL: L1B screened data, where the clock trend is eliminated in the screening process but at the same time the phase observations are altered [67]. Primarily, the latter data are used among highly precise geodetic applications. For the application of RO, however, it is important to note that the native clock offset is required for the correction of the occultation time stamps and range calculation (for removal of GPS receiver oscillator drift in zero differencing). Optimally handled, both data (L1A and L1B) should achieve closely similar results in POD. However, the processing of the large clock trends poses challenges to the POD software packages that might not be optimally handled in the software versions used for this study and somewhat decreases the quality of the resulting orbits.
The internal quality assessment and uncertainty estimation is based on two independent software packages: Bernese GNSS software v5.2 [30], a widely used multi-GNSS data processing software capable of LEO-POD developed at AIUB, and NAPEOS v3.3.1 [31], another state-of-the-art GNSS data processing software developed and maintained by the European Space Operations Centre (ESOC) of ESA. Both software packages are capable of multi-GNSS and support the SLR tracking system for orbit validation. The SLR results presented in this paper are calculated with the Bernese software.
A similar overall processing strategy for LEO POD is inherent to the individual software packages but they vary in detail (Appendix A illustrates the setups used). Furthermore, Bernese was used to propagate random error estimates based on software extensions by AIUB [68,69]. The different models and processing setups of the Bernese and NAPEOS packages as used in this study are summarized in Table 1. The subsequent uncertainty calculations are carried out by separate routines of the rOPS.

Precise Orbit Determination
Overall, there are three different approaches to the orbit determination of satellites in LEO: (1) kinematic; (2) dynamic; and (3) reduced-dynamic. For kinematic POD [70] no force models are applied to the LEO, and the satellite positions are solely derived geometrically from the GNSS observations. The independence from dynamical models makes the kinematic approach attractive for many applications, e.g., gravity field determination [71]. However, kinematic orbit determination is sensitive to data problems in GNSS measurements, bad GNSS constellation (limited visibility), and data gaps [72]. Employing dynamic modeling [73] can overcome the strong dependence on the quality of the GNSS observation. Gravitational and non-gravitational force models are introduced to constrain the estimate of the orbit and ensure a continuous orbit in periods of low visibility and data gaps. The quality of the resulting orbit solution is limited to the quality of the force modeling.
The strong dependence on the force models may be reduced by a so-called reduced-dynamic approach [74][75][76]. By applying additional empirical parameters [68] the satellite trajectory is adjusted to the observations and the residuals are minimized, which results in orbit solutions of highest quality. For the application of RO, reduced-dynamic approaches are favored over a kinematic approach because of the need of RO for the orbit to be consistent, continuous, and of high precision. These requirements ensure the most accurate interpolation of the position, velocity, and clock estimates to the occultation measurements of an average duration of 2 min with a 50 Hz sampling rate. The usage of empirical orbit parameters differs between POD software packages. NAPEOS follows a more dynamical approach with fewer empirical parameters applied compared to the setup of Bernese as used for this study [77].
Introducing CODE and JPL GNSS orbits (15 min) and clocks (30 s sampling, see Section 3.1) the LEO orbit determination is performed using Bernese and NAPEOS. In general the processing setup is similar to both software packages. Bernese employs a discrete kinematic approach for the initial code observations-based orbit and uses orbit integration to obtain a continuous a-priori orbit. In NAPEOS, however, initial values from the implemented Bancroft algorithm [78] serve as a basis for a dynamic orbit fit through the observations. Subsequently, iterative data screening (cycle slip and outlier detection) takes place before the final orbit integration. In Bernese the final orbit is established by means of empirical accelerations to the estimation process [79]. Models and parameters applied in the orbit processing are summarized in Table 1.
For processing raw observation data with large inherent clock estimates, such as for GRACE (cf. Section 3.1), the aforementioned processing scheme of Bernese was extended introducing an additional a-priori orbit calculation and modified data screening options.

Orbit Solutions
In order to enable internal quality assessment and uncertainty estimation for the LEO RO receiver orbit, independent from orbit products of external providers, three POD runs are conducted routinely and in parallel at WEGC (see also Figure 1): • WEGC-BC: The orbit solution selected for the following RO data processing. The calculation is based on Bernese and the GNSS orbit and clock data from the CODE reprocessing in 2015; referred to as primary solution hereafter. • WEGC-BI: Same Bernese setup as for WEGC-BC above but employing GNSS orbit and clocks from JPL (as part of IGS repro2). • WEGC-NC: This solution is calculated using NAPEOS, introducing the GNSS orbit and clock products from CODE as for WEGC-BC.
Among these three orbit solutions the latter two only serve for the quality assessment of the primary solution and the uncertainty estimation through inter-comparison and are not used in the further RO processing.
Additionally, for this study external orbit data from other leading processing centers are used for the validation of the internal primary solution. This comprises orbits from UCAR ( [89,90]), EUMETSAT [28,91], and AIUB [92,93]. First inter-agency comparison at WEGC was carried out within a prototype study of several test days [94]. The extensive UCAR/CDAAC RO data archive holds standard re-processed orbits from their reprocessing campaign 2016 for CHAMP and Metop and post-processed orbits for GRACE. EUMETSAT provided orbit products used for comparison for Metop-A/B and CHAMP. AIUB, as a leading institute for orbit determination provides reliable orbit solutions for the CHAMP and GRACE missions. Table 2 summarizes the different orbit solutions used in this study.

Satellite Laser Ranging
Laser ranging to LEO satellites is a highly valuable complementary technique for validation of GNSS-derived orbits, since it uses completely independent measurements from the global geodetic ground station network of the ILRS [44,95]. SLR observations are obtained by measuring the two-way signal travel time from an SLR ground station to a satellite equipped with a Laser Retroreflector (LRR), which passively reflects the short laser pulses emitted. For most ILRS stations a one-way range precision typically ranging from 0.3-1.8 cm is achieved [58,96]. The signal travel time, corrected for propagation effects and translated into a range R, is then compared to the range from the orbit determination: Due to the measurement geometry the SLR residuals ∆R SLR for high-orbiting GNSS satellites are most sensitive to the radial component of the orbit solution. For satellites in LEO, which are in general tracked at lower elevation angles, however, SLR range is sensitive to all orbital components. In an idealized concept of zero mean error and even distribution of the standard deviation in all directions the 3D position standard error can be approximated by σ pos = √ 3σ SLR [44]. This idealization disregards that LEO passes are more sensitive in north-south direction to the along-track component and in east-west direction to cross-track component [97].
The SLR residuals then serve in our context as a measure for orbit accuracy, supporting the detection and estimation of potential systematic errors. Station coordinates for the observation modeling were taken from the Satellite Laser Ranging Frame 2008 (SLRF2008, https://ilrs.cddis.eosdis. nasa.gov/science/awg/SLRF2008.html) in alignment with the GNSS orbit products used in the POD (Table 1).

Uncertainty Modeling
The estimation of uncertainties attributed to POD poses a challenging task since the accuracy of an orbit in an absolute sense is difficult to determine, but measures for orbit uncertainty can be derived through validation. Hence, as part of the RO processing at WEGC the novel POD setup of the rOPS provides a capability to compute estimated systematic and random position and velocity uncertainties for LEO and GNSS satellites from orbit inter-comparison (internal solutions, Table 2), satellite laser ranging residuals, and propagated random errors in LEO results.
In a conservative approach all systematic uncertainty estimates of an orbit found to be below an a-priori defined conservative bound uncertainty (u CB ) are set to this low-bounding value. Otherwise, if exceeding this bound, they are kept at the derived value (as described below). In addition, the overall smaller estimated random uncertainties contribute in a root-mean-squares (RMS) sense to the total combined uncertainty budget for receiver (Rx) and transmitter (Tx) satellites respectively, as follows: where u r denotes the estimated random uncertainty and u s the total estimated systematic uncertainty. u s is per definition always ≥ u CB as described for the LEO satellites next.

LEO Systematic Uncertainty Estimates
LEO satellite orbits with derived systematic uncertainties within a 3D root-mean-squared (RMS) u Rx CB of 5 cm (position) and 0.05 mm/s (velocity) threshold promote negligible errors into the final atmospheric variables derived from RO. Furthermore, this threshold is regarded as conservative, since it fulfils the RO quality requirement of an error not higher than about 0.05 mm/s in along-occultation-path velocity [3,29], which is the most critical POD variable for RO retrieval quality.
In our approach several sources are exploited for the estimation of the LEO orbit uncertainties. First, the inter-comparison of different orbit products routinely calculated (Section 3.3) allows for an assessment of the generic consistency between the individual orbit solutions. Since the employed POD software packages, processing setup, parametrization, models, and GNSS data used are regarded as sufficiently independent, possible systematic errors are considered plausibly estimated. In addition, independent SLR tracking can also reveal systematic deficiencies in the orbit determination and help identify the most accurate solution. The following three measures are derived on a daily basis for checking the WEGC-BC primary solution: The orbit computations are carried out by Bernese with different GNSS orbit, clock, and EOP data and might reveal deficiencies in the transmitter orbit modeling. (c) WEGC-BC vs. SLR difference (SLR validation cross-check): Analysis of SLR residuals for missions equipped with an LRR. This validation represents snapshots because of the limited number of available SLR measurement sequences per day.
The overall 3D-RMS estimate is calculated from epoch-wise differences in radial, along-track, and cross-track components as well as the RMS from the 1D SLR residuals, resulting in one daily scalar estimate per measure. From these three scalar uncertainty measures for (a), (b), and (c), termed u (a) , u (b) , u (c) hereafter, the final systematic uncertainty estimate u Rx s for the primary orbit solution is derived as follows: Such a variance combination of the uncertainty measures ensures a reasonable RMS estimate even if the individual components are not fully independent (according to the Guide to the Expression of Uncertainty in Measurement by the Joint Committee for Guides in Metrology (JCGM) [98]; see also [99] (pp. 29-44 therein)). Note that laser ranging residuals only contribute their measure u (c) to the calculation if available. The position uncertainties (in units cm) are empirically found (cf. Figure 2) to map with a factor around 1/1000 to the velocity uncertainties (in units mm/s); hence we adopted a conversion factor of 1/1000 for this mapping.

LEO Random Uncertainty Estimates
For the estimation of random uncertainty we consider the error propagation in Bernese based on formal errors of position and velocity as introduced in [68,69]. However, the estimation from the variance-covariance information of the parameter estimation of the POD depends mainly on the observation constellation and hence is considered an optimistic estimate. We therefore conservatively apply a coverage factor of 2 for an increased confidence interval for the estimate. In line with advice in JCGM [98], suggesting coverage factors typically in the range 2 to 3, the role of such a factor is to ensure that non-quantified components of the overall uncertainty are accounted for within the expanded uncertainty range. In our case the factor of 2 serves this purpose. The radial, along-track, and cross-track uncertainty components of position and velocity ( Figure 2) are also combined in a 3D-RMS sense.

GNSS Uncertainty Estimates
The quantification of systematic and random uncertainty estimates for the GNSS satellites is based on the exploitation of existing analyses from the GNSS expert community. For the random uncertainty we introduce a fixed conservative bound estimate of 1 cm, adapted from the analysis of formal errors of GPS [100]. The systematic uncertainty is specified with a 3 cm conservative bound estimate, denoted u Tx CB hereafter, established from SLR investigations [101,102] and long-arc analysis [62,63]. Although only two GPS spacecraft were equipped with an LRR, their position quality is expected to be of similar level as they are less influenced by gravitational forces through their higher orbit. In addition to the conservative threshold values the accuracy codes (AccCode) given in the header of the SP3 orbit files (ftp://igs.org/pub/data/format/sp3c.txt) serve as a daily indicator of formal accuracy and may increase the systematic uncertainty estimate for the individual GNSS satellites as follows: u Tx s = max u Tx CB , 2 AccCode , where the accuracy codes map as an integer base-2 logarithm: AccCode = log 2 (σ) [60]. Analogous to the combined uncertainty of the LEO, the systematic and random uncertainty estimates for the GNSS are combined in a RMS sense (Equation (2)). The GPS satellite orbit uncertainty estimates depicted in Figure 3 suggest a factor around 1/5000 from position to velocity for the random errors. In order to provide also some room for unmodeled systematic effects, the corresponding uncertainty estimates for velocity are conservatively computed with a conversion factor of 1/3000 from GNSS position uncertainty estimates (in units cm) to velocity uncertainty estimates (in units mm/s).

Results and Discussion
Three different orbit solutions are calculated routinely within WEGC's POD processing (listed in Table 2). In order to assess the quality of the orbit solution designated for further RO processing (WEGC-BC), different analyses are performed. The WEGC-BC primary solution is compared to internal control runs (WEGC-BI, WEGC-NC) and orbit solutions from external processing centers (Section 4.1). Satellite laser ranging residuals are analyzed (Section 4.2), and measures for the orbit uncertainties are derived (Section 4.3). The results cover the analysis for CHAMP and GRACE-A in July to September 2008 and the same time period in 2013 for Metop-A/B.

Orbit Comparison
As outlined in Section 3.5, the inter-comparison of the different precise orbit products allows for a realistic assessment of the accuracy and might reveal possible systematic errors. Figure 4 shows daily 3D-RMS differences in position between the WEGC-BC solution and externally provided orbits. The corresponding position and velocity differences in radial, along-track, and cross-track directions are summarized in Table 3. When comparing the different orbit solutions overall best agreement is found for the comparison of WEGC-BC and WEGC-BI. This can be expected from the fact that in this case the differences are simply governed by the impact of different GNSS orbit, clock, and associated EOP data while using the same POD software with consistent processing settings. For all missions the RMS of the daily 3D-RMS series does not exceed 1.4 cm. However, a single day for GRACE-A of about 4 cm difference can be found. The WEGC-NC solution generally compares at a similar level to WEGC-BC as the external solutions, with the exception of decreased differences found between WEGC-BC and AIUB for CHAMP.
CHAMP. For the CHAMP mission best agreement is found with AIUB with an RMS of about 1.8 cm (note: 3 days from AIUB are missing in the provided data). Solutions from WEGC-NC, EUMETSAT, and UCAR compare at a similar level, where the comparison to UCAR features a slightly better comparison with an RMS of about 4.0 cm. Individual days still exhibit differences above the 5 cm threshold.
GRACE. Differences from orbit comparison for GRACE-A presented are higher than results from space-geodetic applications, such as gravity field recovery [46]. However, results from those applications are based on screened observation files (available from JPL) where the clock error was detrended and altered substantially, compared to the raw observation files available at CDAAC. As already outlined in Section 3.1, dealing with such large clock trends poses challenges in handling the data by the POD software packages. This leads to comparison results indicating a slightly degraded quality but still within the conservative bounds for RO. Note that while the WEGC-NC and AIUB solutions are based on these screened L1B products, all other solutions are based on the raw L1A observation data.   In Figure 4b increased values of comparison in September 2008 and missing days end of August for UCAR can be seen for GRACE-A, which is in accordance with periods when the satellite clock error is peaking before it is reset. Also, increased along-track differences (Table 3) compared to solutions based on pre-screened RINEX suggest a shift in time. In particular for the comparison of the primary solution, and the one by AIUB, the use of different low-level input data leads to increased differences for GRACE-A (4.3 cm) compared to the same comparison conducted for CHAMP (1.8 cm). With respect to the AIUB solution as reference, which can be considered a particularly reliable reference due to its K-band validation [69], WEGC orbits compare better than UCAR (WEGC-NC 1.5 cm, WEGC-BC 4.3 cm, WEGC-BI 4.7 cm; UCAR 6.5 cm).
Metop. In principle, a similar characteristics of the comparison results is found for both Metop satellites. Apart from the WEGC-BI solution, best agreement is found against WEGC-NC with an RMS of about 4.0 cm and 4.2 cm for Metop-A and Metop-B, respectively. Individual days of the comparison are slightly exceeding the 5 cm threshold, and furthermore days with increased differences against UCAR can be observed in August and September leading to an overall RMS close to and slightly exceeding 5 cm for Metop-A and Metop-B, respectively. These findings are consistent with earlier results from orbit inter-comparison for Metop [28]. Note that both Metop satellites undergo a maneuver in the considered period. Additionally, two days with missing attitude from CDAAC remain currently unprocessed.
Finally we note that in general all orbit solutions for the presented missions (including those from UCAR, EUMETSAT, and AIUB) satisfy the introduced RO-application-oriented target specification of 5 cm in position and 0.05 mm/s in velocity.

Satellite Laser Ranging Validation
Limited to the CHAMP and GRACE-A missions, which are equipped with an LRR, Table 4 summarizes the laser ranging residuals statistics for the two missions using normal point data available from (1) all stations and (2) a high-quality subset of 12 stations. The station selection was applied in order to discriminate between different performance levels among the stations of the ILRS network [44]. The results comprise the internal and external POD solutions derived from GPS measurements. Figure 5 illustrates SLR residuals that were calculated based on observations from the selected high-quality stations only and the WEGC-BC primary orbit.
In general, the mean and standard deviation for the validated orbit solutions (Table 4) decreases with the limitation to high-quality (HQ) stations. For GRACE-A, however, the WEGC-BC, WEGC-BI, and UCAR solutions, which are all solutions based on unscreened L1A GPS navigation tracking data, show a larger mean for the HQ station selection than for all stations. Furthermore, regardless of the station selection, those solutions exhibit standard deviations more than two times larger than the WEGC-NC and AIUB solutions, which are based on the screened L1B GPS navigation tracking data (see Section 3.1). The SLR residuals in Figure 5 are based on high-quality station subset after applying a 20 cm outlier rejection and an elevation mask of 10 degrees. The high quality stations contribute from 70% (GRACE-A) to 71% (CHAMP) to all normal tracking points, a number generally high since the GRACE and CHAMP missions were top-ranked according to the ILRS mission priority list (https: //ilrs.cddis.eosdis.nasa.gov/missions/mission_operations/priorities/priorities_20080926.html). Overall, the SLR intercomparison results do not exhibit notable systematic variations within the time periods considered. With a similar number of observations (CHAMP: 8271, GRACE-A: 8533) both data series show a small mean. However, the larger scattering which is observed for GRACE-A (Figure 5b) manifests in an increased standard deviation of about 1 cm compared to CHAMP. Furthermore, it was found that using screened L1B RINEX data for the GRACE-A POD reduces the standard deviation (comparable to CHAMP), although this is not applicable in the case of RO (Section 3.1).

Uncertainty Estimation
The results of the uncertainty estimation as introduced in Section 3.5 are illustrated in Figure 6, and summarized in Table 5 (3-month mean and standard deviation range), for the RO LEO satellites under investigation and two representative GPS transmitter satellites possibly involved in an occultation event (one GPS delivering good estimates and one slightly degraded). The component estimates in Table 5 use plausible fractional weighting of the 3D variance (u 2 3D ) [69,89]. In case of a typical day, when the geometry modeling system (Figure 1) delivers results within highest-quality demands for RO processing, the 3D-RMS uncertainties are generally obtained near the conservative bound of 5 cm in position and 0.05 mm/s in velocity for the LEO satellite.
LEO satellites. The estimated random uncertainties for all missions (Figure 6, panels (a) to (d)) show consistent behavior over time, with lowest values found for GRACE. On days where the uncertainty calculation failed due to missing input (i.e., the POD of one of the control orbit runs was not successful) the estimated random uncertainty is set to a conservative estimate of 2 cm, which leads to a slightly increased combined uncertainty estimation on those days for GRACE. For the majority of days, the raw-estimated systematic uncertainty component stays within the conservative bound. In particular this applies to most days for GRACE and the Metop satellites. For CHAMP some more days are found that moderately exceed the conservative-bound threshold.
GPS satellites. With regard to the GPS transmitter satellites ( Figure 6, panels (e) and (f)), the estimated random uncertainties are fixed to the predefined conservative value of 1 cm (cf. Figure 3). The estimated systematic uncertainty, as derived from the accuracy codes, reveals some days of decreased orbit accuracy. Since the estimated uncertainties are empirically mapped with a conversion factor of 1/3000 from position to velocity (Section 4.3), it is well visible that the contribution of GPS velocity uncertainties is quite minor compared to the velocity uncertainties derived for the LEO satellites. Table 5. Random and systematic position uncertainty estimates (3D) and decomposed with 3D-variance weighting factors w into radial (w r = 0.25), along-track (w a = 0.35), and cross-track (w c = 0.4) components.

Conclusions
Climate benchmark data derived from GNSS RO require accurate and robust POD of the GNSS transmitter and LEO receiver satellites taking part in an occultation measurement. In this paper we presented a novel setup for routine quality assessment and uncertainty estimation of the daily LEO receiver satellite orbits independent from external validation sources. The GNSS orbit uncertainty estimates were complemented by building on existing error estimates from the GNSS community. As part of WEGC's rOPS we provide estimates of systematic and random uncertainties associated with the orbit determination, deduced from different LEO POD runs, SLR measurements, and formal uncertainty analysis. We focused on the validation of the WEGC primary orbit solution, which is chosen to deliver the receiver orbit positions, velocities, and clock estimates for subsequent RO processing and the derivation of atmospheric profiles of ECVs.
For the performance assessment of the extended POD setup we investigated 3 representative months of data from July to September in 2008 (CHAMP and GRACE-A) and 2013 (Metop-A/B). The comparison of the WEGC primary solution with orbit solutions from internal control runs showed reliable agreement within 5 cm in position and 0.05 mm/s in velocity, which is the target threshold specification set for high-quality orbit products for RO climate applications. This threshold is also satisfied for most days in an inter-comparison with orbits from the external providers UCAR, EUMETSAT, and AIUB. SLR residuals using selected high-quality ground stations, available and used for CHAMP and GRACE-A, exhibited a mean and standard deviation of −0.15 ± 1.83 cm and 0.17 ± 2.88 cm, respectively, well within the target ranges. Interchange of the GNSS orbit data used with one POD software, without altering the processing setup, shows relatively small differences in the orbit inter-comparison. This suggests that different POD software implementations and configuration settings are more relevant to the uncertainty estimation than potential quality differences in the GNSS orbit data products.
Overall, the uncertainty estimates are found to be within the specified target thresholds for 92% of the days considered. The remaining 8% of days exhibit higher uncertainties of order 5 to 15 cm. This is still adequate for RO processing and at the same time results in somewhat increased uncertainty estimates of the derived ECVs, after propagation through the rOPS retrieval chain. These results suggest a high processing standard and robustness, shared among all processing centers, and indicate rOPS POD system readiness for long-term climate reprocessing of RO data records from CHAMP, GRACE, and the Metop satellite series.
Further modifications of the observation geometry system of rOPS may be applied in future with the inclusion of other RO missions and enhancements of the basic uncertainty estimation. As another important RO mission, we currently integrate POD processing of FORMOSAT-3/COSMIC (Formosa Satellite mission-3/Constellation Observing System for Meteorology, Ionosphere, and Climate) [103]. This will show the retrieval performance implications of the rOPS POD setup for the case of overall degraded POD quality (about 15-25 cm position uncertainties [90]), due to less favorable attitude behavior of the rather small COSMIC spacecraft and restrictions in processing observations from two POD antennas.
For the time being, the further RO processing of FORMOSAT-3/COSMIC data, as well as of the RO data from the more recent FengYun-3 operational satellite series [104] and the commercial CubeSat constellation of Spire [105], is prepared with using the LEO orbits from these data providers together with fixed adopted values of position and velocity uncertainty estimates based on their orbit quality assessments. As part of the future multi-satellites portfolio, we will also process the very recent FORMOSAT-7/COSMIC-2 RO data [106]. Another aspect is that the processing of orbit arcs is currently limited to 24 h arcs, restricting calculation of orbit overlap statistics or long-arc analysis, which in future might serve as an additional measure in the uncertainty estimation process.
As a next step, the upcoming RO reprocessing at WEGC will provide a thorough long-term performance check of the current implementation of this new POD subsystem and may motivate additional modifications and improvements. Based on the encouraging results of this study we expect to obtain high-quality ECV data records for climate monitoring and research.

Acknowledgments:
The authors are grateful to the Center for Orbit Determination in Europe (CODE Bern), the International GNSS Service (IGS), and the Jet Propulsion Laboratory (JPL Pasadena) for their reprocessing effort and the publicly provided access to GNSS orbit, clock, and earth orientation parameter data. We are equally thankful to the International Laser Ranging Service (ILRS) for their effort to collect and publicly provide satellite laser ranging observations. We also appreciate and acknowledge the availability of GPS observations of low Earth orbit satellites shared by UCAR/CDAAC Boulder, JPL Pasadena, and GFZ Potsdam. We thank D. Arnold for the discussion and support in processing with the Bernese GNSS software and H. Peter and T. Springer for discussion and support with the NAPEOS software. We furthermore thank colleagues from WEGC's RO & Climate team, in particular A. Leuprecht, for support on the rOPS framework setup and computations. We also thank F. Ladstädter from WEGC for his help and discussions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:   Figure A1. Schematic illustration of the orbit processing setup of the Bernese and NAPEOS software packages used in this study for the LEO POD. The main processing steps from input data (top) to POD results (bottom) are shown and software-internal routines used for specific steps are indicated by their names (green names near arrows; for details see the software documentation [30,31]). For explanation of acronyms in various boxes see the Abbreviations list above.