Integrating uncertainty propagation in GNSS radio occultation retrieval : from excess phase to atmospheric bending angle profiles

Global Navigation Satellite System (GNSS) radio occultation (RO) observations are highly accurate, long-term stable data sets, and are globally available as a continuous record since 2001. Essential climate variables for the thermodynamic state of the free atmosphere, such as pressure, temperature and tropospheric water vapor profiles (involving background information), can be derived from these records, which therefore have the potential to serve as climate benchmark data. However, to exploit this potential, atmospheric profile retrievals need to be very accurate and the remaining uncertainties quantified and 5 traced throughout the retrieval chain from raw observations to essential climate variables. The new Reference Occultation Processing System (rOPS) at the Wegener Center aims to deliver such an accurate RO retrieval chain with integrated uncertainty propagation. Here we introduce and demonstrate the algorithms implemented in the rOPS for uncertainty propagation from excess phase to atmospheric bending angle profiles, for estimated systematic and random uncertainties, including vertical error correlations and resolution estimates. We estimated systematic uncertainty profiles with the same operators as used for the ba10 sic state profiles retrieval. The random uncertainty is traced through covariance propagation and validated using Monte-Carlo ensemble methods. The algorithm performance is demonstrated using test-day ensembles of simulated data as well as real RO event data from the satellite missions CHAMP and COSMIC. The results of the Monte-Carlo validation show that our covariance propagation delivers correct uncertainty quantification from excess phase to bending angle profiles. The results from the real RO event ensembles demonstrate that the new uncertainty estimation chain performs robustly. Together with the other 15 parts of the rOPS processing chain this part is thus ready to provide integrated uncertainty propagation through the whole RO retrieval chain for the benefit of climate monitoring and other applications.


Introduction
Observation systems of the free atmosphere, focusing on the range from the top of the atmospheric boundary layer upwards, were historically designed for weather research and forecasting purposes.They have considerable shortcomings from a climate monitoring perspective (Karl et al., 1995), and so the related global climate monitoring infrastructure remains fragile and incomplete even today (Bojinski et al., 2014).The Global Climate Observing System (GCOS) aims to improve the observational foundation for the climate sciences (GCOS, 2015).For this purpose the establishment of climate benchmark data records is essential.To qualify as climate benchmark, records need to be (1) of global coverage, (2) of high accuracy, (3) long-term stable, (4) tested for systematic errors on orbit, (5) and tied to irrefutable standards, and they need to (6) measure Essential Climate Variables (ECVs) (NRC, 2007;GCOS, 2015).
Based on the quality and abundance of Global Navigation Satellite System (GNSS) signal sources, in particular from the Global Positioning System (GPS) so far, the GNSS radio occultation (RO) observation record is globally available (continuously since 2001), long-term stable (due to the so-J.Schwarz et al.: Uncertainty propagation in RO retrieval called self-calibration and high signal stability during the event), and highly accurate (accuracy traceable to the SI second).Due to the self-calibrating property, the accuracy is also ensured on orbit; i.e., there is no need for calibration or bias correction in post-processing on ground (Leroy et al., 2006).The basic RO excess phase data can therefore serve as a fundamental climate data record (FCDR) as defined by GCOS (2010).From this FCDR with its unique properties, ECVs -in particular the thermodynamic ECVs pressure, temperature, and humidity in the free atmosphere -can be derived using an RO retrieval chain.
In order to reliably serve as climate benchmark data record, however, the retrieved ECV profiles and their claimed accuracy -expressed by the uncertainties provided -need to be traceable back to the (small) uncertainties of the FCDR and in turn to the raw data.This requires that (1) the RO retrieval is highly accurate and avoids any undue amplification of uncertainties associated with the quantities in the FCDR and that (2) the uncertainties are propagated through the entire retrieval chain, from the raw data to the ECV profiles, duly accounting for relevant side influences such as from background information.Developed at the Wegener Center of the University of Graz (WEGC), together with international partners, the Reference Occultation Processing System (rOPS) (Kirchengast et al., 2015) aims to establish such a fully traceable RO processing for the first time (Kirchengast et al., 2016a, b).
In Fig. 1 the basic steps of the RO retrieval chain in the rOPS, i.e., the precise orbit determination (POD) and excess phase processing (labeled "L1a" in Fig. 1), the subsequent atmospheric bending angle retrieval ("L1b"), the refractivity and dry-air retrieval ("L2a"), and the moist-air retrieval ("L2b") are sketched.Kursinski et al. (1997) -and more recently Hajj et al. (2002), Anthes (2011), andSteiner et al. (2011) -provided detailed introductions and reviews of the RO technique and its applications in meteorology and climate.Ho et al. (2012) and Steiner et al. (2013) included comparative current RO retrieval chain descriptions of leading international RO processing centers, none of which yet include uncertainty propagation.Empirical error (uncertainty) estimates computed statistically from retrieved RO atmospheric profiles and climatologies have been derived by Kuo et al. (2004), Steiner and Kirchengast (2005), and Scherllin-Pirscher et al. (2011a, b, 2017), the last of which with a focus on climate uses also providing simple analytical error models.These studies and many others have described the RO retrieval chain in detail and have shown the high accuracy and quality of RO data, particularly in the upper-troposphere and lower-stratosphere region.
The aim of the integrated uncertainty propagation in the rOPS is to eventually propagate uncertainties along this entire retrieval chain from the raw measurement data to the ECVs (Kirchengast et al., 2016a, b), whereby the implementation of the rOPS uncertainty propagation occurs in the se- .Schematic view of the main processors of the retrieval chain in the rOPS (L1a, L1b highlighted, L2a, L2b) and the main operators of the L1b processor (1, 2, 3), which are in the focus of this study.quential blocks illustrated in Fig. 1.The L2a processing and uncertainty propagation from atmospheric bending angle to dry-air profiles has already been introduced by Schwarz et al. (2017, SKS2017 hereafter).
This study is a direct complement to the work in SKS2017.Using the same propagation and validation methods as applied in SKS2017, it focuses on the uncertainty propagation from excess phase to atmospheric bending angle profiles, i.e., the L1b processing.As in SKS2017, random uncertainties are propagated using covariance propagation (CP) and validated using Monte Carlo (MC) ensemble methods.As in the L2a processor, we also propagate (conservative bound) estimates for systematic uncertainties along the retrieval chain of the L1b processor.Additionally, correlation length profiles and resolution profiles are provided.
Uncertainty propagation as covariance propagation from excess phase to bending angle profiles has been outlined and demonstrated in a basic form, by Syndergaard (1999) and Rieder and Kirchengast (2001), but not been implemented yet in processing center retrieval chains or applied to real RO data.As visible in Fig. 1, the L1b processor consists of three major retrieval parts, which are expanded into detailed substructure in Fig. 2. We propagate estimated random uncertainties from excess phase profiles to Doppler shift profiles (Sect.(1) in Fig. 2); further to geometric-optics (GO) bending angle profiles, merged with wave-optics (WO) bending L1a observed excess phase processing output file

Random uncertainty estimation
(1) Doppler shift retrieval (2) Bending angle retrieval . Detailed workflow for state retrieval and uncertainty propagation of the main L1b operators from excess phase to atmospheric bending angle profiles (1)-(3) and of the subroutines used in the MC testing framework (a)-(g).The mathematical notation, including all symbols, is introduced in Tables 1 and 2.
angle profiles (2); and finally to atmospheric bending angle profiles (3), using a full-CP approach.In combination with the definitions of the main operators and variables in Table 1, and of the vertical grid and coordinate variables in Table 2, Fig. 2 provides a concise overview on the detailed workflow of the L1b processor.Uncertainty propagation for the WO bending angle retrieval has been implemented and demonstrated for simulated events by Gorbunov and Kirchengast (2015); estimation of , and unit U ∈ {m, m, m s −1 , rad, rad, rad, rad, rad}, comprising elements X r,i .u s X U estimated systematic uncertainty profile of X (with X and U as defined above), comprising elements u s X,i (including estimated basic and estimated apparent systematic uncertainties, u b X,i and u a X,i ).u r X U estimated random uncertainty profile of X (with X and U as defined above), comprising elements u r X,i .R X 1 error correlation matrix of X (with X as defined above), comprising elements R X,ij .C X U 2 error covariance matrix of X (with X and U as defined above), comprising elements m correlation length profile of X (with X as defined above), comprising elements l X,i .τ X s resolution profile of X (with X as defined above) in time domain, comprising elements τ X,i .w X m resolution profile of X (with X as defined above) in altitude domain (along impact altitude), comprising elements w X,i .
X m U model excess phase/Doppler shift/bending angle profiles based on forward modeling of collocated refractivity profiles from ECMWF short-range forecast fields, with X m ∈ {L m (t), D m (t), α m (z a )}, and U ∈ {m, m s −1 , rad}, comprising elements X m,i .x S U profiles of Cartesian position/velocity vectors of the receiving/transmitting satellite relative to the center of refraction, with x S ∈ {r estimated (systematic) uncertainty profiles of x S (with x S and U as defined above), comprising elements u s xS,i .
A BWS 1 BWS filter matrix operator, comprising the Blackman windowed sinc (BWS) low-pass filter weights (normalized filter functions) in the form of a band matrix.A L2D s −1 Doppler differentiation matrix operator, transforming the filtered excess phase profile to the Doppler shift profile.
random and systematic uncertainties for real events including boundary layer bias correction is introduced by Gorbunov and Kirchengast (2018).
Other ongoing rOPS retrieval advancements relevant to this study are the inclusion of the high-altitude initialization algorithm, introduced by Li et al. (2013Li et al. ( , 2015)), in the L2a processor and the reduction of remaining higher-order ionospheric effects in the retrieved bending angle profiles of the L1b processor (based on work by Syndergaard, 2000;Liu et al., 2015;Healy and Culverwell, 2015;andDanzer et al., 2013, 2015).Furthermore, the precise orbit determination (POD) of the RO receiver satellite and the excess phase processing, also including the associated uncertainty propagation, are part of ongoing work (Innerkofler et al., 2017).
Finally, related work and manuscript preparation on a new moist-air retrieval algorithm (L2b) and corresponding L2b uncertainty propagation are ongoing (Li et al., 2018;Kirchengast et al., 2017a).
The paper is structured as follows.In Sect. 2 we introduce the uncertainty estimation, propagation, and validation methods and the data sources and preparation.In Sect.3, with the help of an example RO event, the uncertainty propagation sequence is introduced.In Sect. 4 we present the results from the MC validation of the CP uncertainty estimates.In Sect. 5 the performance of the algorithm is then evaluated using test day ensembles with real data from the RO missions CHAllenging Minisatellite Payload (CHAMP) (Wickert et al., 2001); FORMOSAT-3 Constellation Observing System for Meteorology, Ionosphere, and Climate (COS-MIC) (Anthes et al., 2008); and Meteorological Operational Satellite A (MetOp) (Luntama et al., 2008); as well as with simulated data approximating characteristics of the Meteorological Operational Satellite A (Luntama et al., 2008;simMe-tOp data hereafter).We close with conclusions and outlook in Sect.6.A detailed description of the implemented uncertainty propagation algorithms can be found in Appendix A. 2 Methods and data

Methods
We follow the Guide to the Expression of Uncertainty in Measurement (JCGM, 2008a(JCGM, , b, 2011;;GUM hereafter) and aim to follow terminology as provided by the International Vocabulary of Metrology (JCGM, 2012), a terminology also adopted by the GUM.SKS2017 provides a more thorough introduction, including the motivation for using the respective uncertainty estimation, propagation, and validation methods; we refer the particularly interested reader to this companion (open-access) work and provide the essential methods needed more in a summarized form below.
We categorize uncertainties into estimated random uncertainties and estimated systematic uncertainties.Effects of unpredictable or stochastic temporal and spatial variations in repeated observations, like effects from fluctuations in the atmosphere or the thermal noise of the receiver system, could in principle be estimated by ensemble statistics from multiple RO events.However, since such effects are essentially stationary in a statistical sense, we can estimate their statistics also from individual RO event data, given their high noiseresolving sampling rate.These effects are included in the estimated random uncertainties.
Systematic effects (biases), which can not be quantified using statistical data analysis based on just one individual RO profile, are estimated and corrected for when known, as recommended by the GUM.The remaining residual biases are assumed to stay within a (conservative) bound estimate, which we refer to as estimated systematic uncertainty and by which we aim to provide at least 90 % likelihood coverage (confidence) that residual biases stay within the plus/minus envelope range of this uncertainty.
Depending on their nature, components of the systematic uncertainty that we need to estimate can be fundamentally systematic across different RO events, a subtype we term estimated basic systematic uncertainties, or appear systematic just for individual RO events, a second subtype that we term estimated apparent systematic uncertainties.It is important to distinguish these two subtypes, since the apparent systematic uncertainties will essentially behave as random uncertainties in ensemble averaging over many RO events, such as when generating climatologies, while the basic systematic ones will not average out and therefore fundamentally limit the (absolute) accuracy of ensemble averages such as climatologies.
Since the noise-type effects giving rise to short-rangecorrelated random uncertainties can be considered uncorrelated to the bias-type effects inducing long-range-correlated apparent systematic uncertainties, and since both are uncorrelated to basic systematic uncertainties, it is insightful and possible with due care to estimate and propagate each of these uncertainties independently.
As for the L2a processor (SKS2017), the operators of the L1b processor (i.e., the boldfaced items 1.2, 1.4, 2.1, 2.7, 2.9, 3.1, and 3.5 in Fig. 2) qualify as explicit, multivariate, and linear measurement models, as defined in the GUM, with correlated input quantities.They can therefore be formulated as where the input quantity X and output quantity Y are rank 1 vectors (profiles) of random variables, which we call state profiles.According to the GUM, their random uncertainties can be propagated using when the uncertainties are normally distributed.This assumption is reasonably justified, since the receiving system noise (i.e., thermal noise and residual clock estimation noise) and the ionospheric noise (from scintillations induced by ionospheric irregularities) are essentially normally distributed overall (Kursinski et al., 1997;Syndergaard, 1999;Gorbunov, 2002b;Sokolovskiy et al., 2009).These noise sources are the main contribution to the random uncertainties in the excess phase profiles feeding into the L1b processor.C X and C Y are the covariance matrices of the input and output variables, respectively, and A XY is the linear (or linearized) operator connecting X and Y .Equation (3) formulates how the covariance matrix C X is calculated from random uncertainty estimates u r X and the correlation matrix R X : As a key variable characterizing R X , correlation length profiles l X are estimated from the correlation functions assembled in R X .The algorithm used estimates l X by searching for the distances downward and upward of the correlation functions' main peak at which the correlation function has dropped to a value of 1/e (≈ 0.378).The adopted correlation length estimate is the arithmetic mean of these two upward and downward estimates (as the peak may be somewhat asymmetric).Additionally the correlation length is constrained by the data domain; i.e., the correlation length can never be larger than the profiles' vertical range.
Since the covariance propagation of random uncertainties requires extensive matrix multiplications for each measurement model along the entire retrieval chain, we also tested simpler variance propagation (VP), for which correlations are ignored; Appendix B summarizes the relevant algorithms.However, as shown in Sect.4, variance propagation unduly overestimates random uncertainties, so that covariance propagation is required.
When the operator is linear, as is the case for the applicable L1b operators, estimated systematic uncertainties can be propagated by application of the state retrieval operator on the estimated systematic input uncertainty: where u s X and u s Y are the rank 1 systematic uncertainty profiles of the input and output variables.
In addition to random uncertainties, systematic uncertainties, and the correlation length, we also estimate resolution profiles w X as context information along with the provided random uncertainty profiles (necessary, e.g., because smoothing can decrease random uncertainties, while making resolution coarser).This is enabled by careful selection and formulation of low-pass filter operations, in particular explicit filter cutoff frequency specification as the main driver of the resolution remaining after low-pass filtering.
We note that the (half-)Fresnel scale physical resolution often ascribed to RO bending angle profiles retrieved by geometric-optics methods (e.g., Kursinski et al., 1997;Gorbunov et al., 2004) will generally be somewhat coarser than the filter-limited resolution estimated here.This is intentional to maximize available information in the bending angle profiles provided by the L1b processor.In the rOPS, on input to the L2a processor and before high-altitude initialization by statistical optimization, the resolution of all profiles is brought to a common altitude-dependent resolution, which reflects the half-Fresnel scale (SKS2017).

Data sources and preparation
The input variables needed for the L1b uncertainty propagation, visible in Fig. 2 and defined in Table 1, are the retrieved excess phase profiles L r,k (t) and the associated systematic uncertainty profiles u s L r,k (t), random uncertainty profiles u r L r,k (t), and correlation matrices R L r,k , as well as the orbit positions and velocities of receiver and transmitter satellite -r R (t), v R (t), r T (t), and v T (t) -and their (systematic) uncertainties, u s rR (t), u s vR (t), u s rT (t), and u s vT (t).For due limitation of depth of workflow detail in Fig. 2 we do not separately show the propagation of the basic and apparent systematic uncertainties as they are both identically propagated through the operator chain shown for u s L r,k (t).All variables are provided on the time grid t with elements t i , at f s = 50 Hz sampling rate, and for the two GPS carrier frequencies: f Tk , with k ∈ {1, 2}, f T1 = 1.57542GHz, and f T2 = 1.22760GHz.
We used excess phase state profiles L r,k (t) and the orbit state profiles r R (t), v R (t), r T (t), and v T (t) from 15 July 2008 as a test day ensemble.For CHAMP, COSMIC, and MetOp, orbit state and excess phase profiles were provided by the COSMIC Data Analysis and Archiving Center (CDAAC) of the University Corporation for Atmospheric Research (UCAR), Boulder, Colorado.The End-to-End GNSS Occultation Performance Simulation and Processing System (EGOPS) (Fritzer et al., 2009) was used for generating the simulated MetOp orbit state and excess phase profiles with realistic receiver noise (simMetOp).Figure 3 shows Exploiting the linearity of the (linearized) retrieval operators, the so-called baseband approach (Kirchengast et al., 2016a) is applied throughout the rOPS.Hereby a zero-order model profile is subtracted from the input state profile, and only the remaining delta profile is processed through the operator.After application of the operator, the zero-order model profile of the output state profile is added back to the resulting delta profile.This approach effectively avoids biases from numerical operations on (near-)exponentially varying RO profiles, since the model profiles that we derive from short-range (24 h) forecasts of the European Centre for Medium-Range Weather Forecasts (ECMWF) skillfully subtract the (near-)exponential variation.The remaining increment profiles that we need to treat numerically then appear to be very linear and with low dynamical range, which leads to very low residual numerical errors of operators such as filters and derivatives.
The model profiles used as zero-order states in the retrieval -i.e., L m , D m , and α m (cf.Table 1) -were created from ECMWF short-range (24 h) forecast refractivity fields, accurately forward-modeled to bending angle (α m ), Doppler shift (D m ), and excess phase (L m ) profiles, collocated to the latitude, longitude, and time of the respective RO event processed in the rOPS.The ECMWF fields used have a horizontal resolution of about 300 km (triangular truncation T42)which corresponds to the approximate horizontal resolution of RO profiles (e.g., Kursinski et al., 1997) -and are available at 91 vertical levels (L91).
ECMWF fields were chosen for their proven leading quality (Untch et al., 2006;Bauer et al., 2015) and thus high suitability for serving as zero-order state profiles; any other reasonable model profiles could be used as well since the retrieval results negligibly depend on the exactly chosen zero-order model profiles.For comparison we plotted L m (t) for the COSMIC example case into Fig.3a, which demonstrates that the ECMWF short-range forecast lies very close to L r1 (t) and L r2 (t) and thus is well suited as the model profile.
While in the future the excess phase random and systematic uncertainty profiles will be more rigorously estimated by the rOPS L1a processor (Innerkofler et al., 2017) and provided as input to the L1b processor, they had to be estimated for this study from existing excess phase profiles with realistic noise and simplified modeling.To this end, each estimated random uncertainty profile u r L r,k (t) was estimated based on the noise of the respective retrieved excess phase profile L r,k (t).The noise was determined following the estimation scheme for bending angle observation errors described by Li et al. (2015, Sect.2.2 therein), so we just briefly summarize how we used it here.
First, for both, the retrieved profile L r,k and for the model profile L m , the mean over all grid points between 60 and 70 km was determined.Then L m was offset-corrected towards L r,k by subtracting the difference of these two means from L m , giving the offset-corrected model profile L m .Next, the delta profile δL r m,k = L r,k − L m was calculated.After smoothing δL r m,k with a 10 km moving-average boxcar (BC) filter, the smoothed profile was subtracted from δL r m,k again, to get δδL r m,k , the random noise profile component of L r,k isolated in this way.Finally, the estimated random uncertainty was determined as where M is the number of grid points equivalent to a window width of 10 km.To avoid boundary effects of the filter, u r L r,k was only determined up to z aTop − 5 km and down to z aGradr at 30 km.It was constantly extended at the upper end and extended by a linear gradient below z aGradr , using (in units [m]) for all elements of u r L r,k (t) below z aGradr , roughly following estimates of ESA/EUMETSAT (1998) and the overall behavior of estimates from real excess phase profiles (the latter became too vulnerable to biases and fluctuations to continue using them below 30 km).
Since the noise components responsible for the random uncertainty at excess phase level are essentially uncorrelated at a sampling rate of 50 Hz (Syndergaard, 1999;Hajj et al., 2002), the correlation matrix R Lr is set to unity in the diagonal and to zero outside (i.e., a Kronecker-δ assignment) for both channels (Fig. 3d-f).In case the future excess phase data from the rOPS L1a processor exhibit non-negligible correlations for some data from some of the RO missions, we will account for these correlations in R Lr , since our L1b algorithm (Sect.3) is prepared for full covariance propagation.The elements of the covariance matrix C Lr are hence (item 1.1 in Fig. 2) For the MC validation of the CP, error profile realizations r Lr were superimposed onto simulated "true" excess phase profiles L T r,k (t).As a source for L T r,k , we used an EGOPS-simulated "error-free" CHAMP event from 8 August 2008 (i.e., no receiver system errors superimposed).Using an "error-free" profile as the basis, with the particular simulated profile just serving as a representative RO profile to illustrate the MC validation, allows us to strictly ensure the consistency of the random uncertainty of the input profile with the ensemble of superimposed error profile realizations.
To create the error profiles, a representative u r,Std Lr uncertainty profile was selected from a COSMIC ensemble of uncertainty profiles, created according to Eqs. ( 5) and ( 6).The error profile realizations are random draws from a distribution characterized by these uncertainties, again assuming that R Lr,ij = δ ij , i.e., that there are no correlations between the individual grid levels (item (a) in Fig. 2; Fig. 3f).The same standard profile u r,Std Lr was used as input for the CP to which the MC results are then compared.This MC validation method applied to test the rOPS L1b uncertainty propagation steps is essentially the same as in SKS2017, and it is described therein in more detail.
The estimated systematic uncertainty u s Lr,k was determined based on a simple model roughly following error estimates from ESA/EUMETSAT (1998), with constant uncertainty from 80 km down to z aGrads at 8 km and a linear uncertainty gradient in the troposphere; as noted above, this simplified modeling will be replaced in the future by realistic uncertainty estimates received as L1b retrieval input from the L1a processor (Innerkofler et al., 2017).
The constant u s Lr,k above z aGrads is 0.1 mm for k = 1 and 0.2 mm for k = 2 for MetOp and simMetOp.This uncertainty is interpreted as an estimated basic systematic uncertainty, i.e., as a lower-bound estimate of available accuracy.
For CHAMP and COSMIC we set u s Lr,1 = 0.2 mm and u s Lr,2 = 0.4 mm, to roughly reflect the fact that these RO receivers are lower-cost instruments with lower gain, and thus somewhat lower tracking performance, than the RO receiver on MetOp (e.g., Luntama et al., 2008;Angerer et al., 2017).From z aGrads downwards, u s Lr,k (in units [m]) increases by In order to avoid a sharp kink in the u r L r,k profiles at z aGradr , and in the u s L r,k profiles at z aGrads , a 2 km width movingaverage boxcar filter was applied to smooth these simple un-certainty models around these transition altitudes (for the example profile u s Lr is visible in Fig. 3b).The orbit position and velocity uncertainties of the transmitter and the receiver satellites show little variation within the short duration of an individual RO event of about 45 s to 2 min (Innerkofler et al., 2017) and can be assumed to be constant biases.They are thus counted to the systematic uncertainties, more precisely the apparent systematic uncertainties, since the actual values of the orbit-borne biases will generally change in a pseudo-random manner from event to event.
We set the transmitter position and velocity uncertainties to u s rT = 3 cm and u s vT = 0.01 mm s −1 , consistent with accuracies for GPS orbits available from GNSS orbit providers like the International GNSS Service (IGS).The receiver position and velocity uncertainties, u s rR = 5 cm and u s vR = 0.05 mm s −1 for CHAMP and MetOp, are adopted 4 times smaller than those for COSMIC with u s rR = 20 cm and u s vR = 0.2 mm s −1 , as found by ongoing rOPS-related POD studies (Innerkofler et al., 2017), consistent with previous literature (e.g., Montenbruck et al., 2009;Schreiner et al., 2010).

Algorithm sequence and example results
In this section the L1b uncertainty propagation algorithm sequence is introduced.We illustrate the effects of the algorithm on the main uncertainty variables by way of the COS-MIC example case already used for Fig. 3.
For each L1b retrieval step -i.e., segments (1), (2), and (3) in Fig. 2 -the results for the principal variables are shown in Figs. 4 to 8.These variables are the state profiles X r (with , and α r (z a )}), the estimated systematic uncertainty profiles u s Xr , the estimated random uncertainty profiles u r Xr , representative correlation functions R Xr,i (with i such that z a,i ∈ {10, 30, 50, 70 km}), and the correlation length profiles l Xr and resolution profiles w Xr .Along with the dual-frequency state profiles, we also show the collocated forward-modeled short-range forecast profiles, i.e., model profiles X m with A concise definition of the variables involved is provided in Table 1, as introduced above.The summary description in this section is complemented by a complete step-by-step description of the algorithm along the entire L1b retrieval chain in Appendix A, which is organized for convenience into the same sequence of subsections.
To simplify the notation in the description, we suppress index k whenever steps are applied in an identical way to the data of both GNSS L-band channels with frequencies f T1 and f T2 .Only if the two channels are treated differently, such as in Sect.3.3, is the index considered again.For conciseness we also do not illustrate both the estimated basic and estimated apparent systematic uncertainty but rather the total estimated systematic uncertainty as the overall result.

Basic low-pass filtering
A Blackman windowed sinc (BWS) low-pass filter with a filter cutoff frequency f c = 2.5 Hz (boxcar-equivalent filter width of 0.2 s) (item 1.2 in Fig. 2) is applied onto the excess phase profile L r (t), before the Doppler differentiation (item 1.4 in Fig. 2), to avoid an amplification of highfrequency noise in the phase profile by the derivative operation.This filter suppresses the noise; consequentially the filtered excess phase profile L F (t), shown in Fig. 4a, is expected to have random uncertainties u r LF of smaller magnitude, but correlated over the length of the filter window.The uncertainties obtained through the implemented algorithm confirm these expectations, i.e., random uncertainty profiles in Fig. 4c are less than a third in magnitude of those in Fig. 3c, and Fig. 4d-e show how the correlation functions widened and the correlation length and vertical resolution reached ∼ 0.5 and ∼ 0.6 km, respectively, above about 30 km impact altitude (Fig. 4f).
The random uncertainty propagation algorithm, i.e., the covariance propagation from C Lr to C LF , is described by Eq. ( A6) and item 1.3 in Fig. 2 and is justified by Eq. (2).To obtain u r LF and R LF , we use Eqs.(A7) and (A8).To propagate the estimated systematic uncertainty u s Lr , which characterizes long-range-correlated offsets or biases, we use the same BWS filter as for the state profile, i.e., making use of Eq. ( 4).Because the input uncertainty profile u s Lr is chosen to be constant down to z aGrads , the filter has little effect, and u s LF , shown in Fig. 4b, is essentially equal to u s Lr , shown in Fig. 3b.
The resolution profile w Lr is determined by the filter width according to Eqs. (A11) and (A13).After the BWS filtering, the resolution is roughly equal to the correlation length l LF , amounting to ∼ 0.6 km above about 30 km impact altitude and becoming finer downwards due to the increasing refraction (Fig. 4f).

Doppler shift derivation
The next step is a five-point differentiation operation (item 1.4 in Fig. 2) used to calculate the Doppler shift profile D r (t) from the filtered excess phase profile L F (t).The resulting dual-frequency Doppler shift profiles are plotted along with the model profile D m (t) in Fig. 5a for the example case.
As for the filtered excess phase, we apply CP (Eq.A18, item 1.5 in Fig. 2) to first calculate the covariance matrix C Dr and then extract u r Dr (shown in Fig. 5c) and R Dr (Fig. 5d  and 5e) from it.The choice of the x-axis range shows the random uncertainties increased, but the differentiation actually does increase relative random uncertainties (relative to the state profile).It also causes anti-correlation with neighboring elements, as visualized by the negative side peaks of the correlation functions in Fig. 5d   length l Dr (of the main correlation function peak) decreases accordingly (now smaller than 0.3 km throughout), because the correlation functions fall off more steeply on both sides of the main peak (Fig. 5f).
For calculating the estimated systematic uncertainty, we use the state operator; i.e., we just differentiate u s LF and get u s Dr (shown in Fig. 5b).With the current illustrative choice of input uncertainties the systematic uncertainty of the Doppler shift profile is zero above the transition to the troposphere, where the estimated systematic uncertainty of the excess phase is assumed constant; in the troposphere a Doppler shift offset of ∼ 0.02 mm s −1 occurs.
The resolution profile w Dr shows that the vertical resolution stays unaffected by this operator (cf.Figs.5f and 4f), because the BWS filter width of the preceding low-pass filtering (intentionally) stretched beyond the five neighboring points involved in the differentiation.

GO bending angle retrieval
The next operator is the GO bending angle retrieval in which retrieved GO bending angle profiles α G (t) are calculated from Doppler shift profiles D r (t) and the orbit position and velocity vectors r T (t), r R (t), v T (t), and v R (t) (item 2.1 in Fig. 2) and then interpolated to the (common monotonic) impact altitude grid z a (item 2.6 in Fig. 2).
Figure 6a shows retrieved α G profiles together with the model profile α m .The mildly nonlinear implicit-type bending angle retrieval operator needs to be solved iteratively, and it requires linearization for both random and systematic uncertainty propagation, as described in detail in Appendix A (Sect.A2).Because this retrieval step is performed level by level, keeping levels independent, the GO bending angle retrieval leaves correlation functions and resolution unchanged (cf.Figs.6d-f and 5d-f).
The estimated random uncertainties u r αG , as shown in Fig. 6c, now increase more strongly in the lower stratosphere and troposphere (to about 40 to 50 µrad near 10 km), because they are dependent on the vertical gradient of the impact parameter a t , which is increasingly larger towards lower altitudes from the increasing refraction.
The main contributions to the estimated systematic uncertainty u s αG are induced by systematic uncertainties in orbit velocity and position of the transmitter and the receiver satellite (details in Sect.A2), which in total amount to about 0.05 µrad (Fig. 6b).Compared to this magnitude, the systematic uncertainty contributed by the Doppler shift uncertainty is very small.

WO bending angle retrieval
Due to strong refractivity gradients and multipath effects, the GO bending angle retrieval can be inadequate in the troposphere, and therefore WO algorithms are applied to reconstruct the geometric optical ray structure of the wave field (e.g., Gorbunov, 2002a;Gorbunov and Lauritsen, 2004).
In the rOPS, along with the WO bending angle profile α W (z a ), the systematic uncertainty profile u s αW , the random uncertainty profile u r αW , the correlation matrix R αW , and the resolution profile w αW are retrieved (item 2.7 in Fig. 2).
The WO bending angle retrieval algorithm used is a canonical transform (CT2) algorithm (Gorbunov et al., 2004), and the associated uncertainty propagation algorithm is not described here, but separately by Gorbunov andKirchengast (2015, 2018).The WO retrieval and uncertainty propagation results are supplied up to 20 km impact altitude by the WO algorithms.

Merging of GO and WO bending angle profiles
In the rOPS bending angle retrieval the results from the WO retrieval, α W , are merged with GO retrieval results, α G , at around a transition altitude z GW a in a transition range z GW a ± z GW a , to get merged profiles α M (item 2.9 in Fig. 2).
The determination of the transition altitude and the merging algorithm are described in Appendix A2.3.We use a specialized covariance propagation to propagate the GO and WO uncertainties, expressed by the covariance matrices C αG and C αW , to properly obtain the covariance matrix of the merged bending angle C αM (Eqs.A37 and A38, item 2.10 in Fig. 2).
Because the rOPS implementation of the WO uncertainty propagation (Gorbunov and Kirchengast, 2018) was still in test phase and not yet available for integration into the simulations here, all examples in this study are GO-only; i.e., only the GO retrieval is performed.Results for α M are thus unchanged from those shown in Fig. 6 and not separately illustrated.
In order to nevertheless test and validate the uncertainty propagation of the merging algorithm, WO retrieval results were artificially substituted by the GO results for the MC validation (Sect.4); i.e., GO was used as a proxy for WO since reasonably capturing expected WO variability as indicated by tests of Gorbunov and Kirchengast (2018).

Adaptive low-pass filtering and minor-channel extrapolation
To prepare the merged bending angle profiles α M,k for the ionospheric correction, they are first filtered by another BWS filter operation (item 3.1 in Fig. 2) in order to ensure adequately smoothed bending angle profiles α F,k , with k ∈ {1, 2}.
The chosen filter cutoff frequency for k = 1 is f c1 = 2.5 Hz, same as the basic filtering (Sect.3.1.1),just to ensure clearness of any higher-frequency effects from operators after the initial excess phase filtering (e.g., from Doppler shift derivation that induces short-range anti-correlation effects).For k = 2, the cutoff frequency f c2 is set noise-dependent, between 2.5 and 0.5 Hz (boxcar-equivalent width of 0.2 to 1.0 s).In events in which the α F2 profile does not reach down as far as α F1 , it is extrapolated down to the bottom of α F1 , z aBot .The results for the filtered bending angle state profiles α F,k are displayed in Fig. 7a, together with the associated model bending angle profile α m .The filter has considerably reduced the noise of the profile, particularly for α F2 , where the algorithm selected a cutoff frequency f c2 = 10/7 Hz in this example case.
The relevant covariance-propagated random uncertainties u r αF,k are shown in Fig. 7c (blue and red), illustrating the reduced noise, especially for α F2 .In return, the peaks of the correlation functions broaden (cf.Figs.7d-e and 6d-e), with correlation lengths l αF,k at near 0.4 km for α F1 and above 0.5 km for α F2 (Fig. 7f).
The estimated systematic uncertainty remains largely unchanged (Fig. 7b) due to its smooth character.
The resolution of the filtered bending angle profiles (according to Eqs.A11 and A13) is determined by the cutoff frequencies f c,k of the BWS filters.In the example case it is therefore essentially unchanged for α F1 , while it is significantly decreased for α F2 (cf.Figs.7f and 6f) since f c2 = 10/7 Hz.That is, the resolution w αF2 in the upper stratosphere for example, where the vertical scanning velocity of this RO event is about 3.2 km s −1 , is near 1.1 km (Fig. 7f).

Ionospheric correction
The final step of the L1b processor is the ionospheric correction (item 3.5 in Fig. 2).The atmospheric bending angle α r is obtained by applying a linear dual-frequency combination of α F1 and α F2 , such that ionospheric effects are largely removed (details are described in Sect.A3).The final retrieved atmospheric bending angle α r of the example case is shown in Fig. 8a.The propagation results for the estimated random uncertainty are shown in Fig. 8c.The linear combination of the ionospheric correction amplifies noise and u r αr is therefore considerably larger than u r αF1 , and u r αF2 (cf.Figs.8c and  7c).
Figure 8d shows how the correlation functions -as obtained through covariance propagation -combine the characteristics of the correlation functions from the two matrices R αF1 and R αF2 , while essentially inheriting the α F1 behavior, since the α F2 influence on the ionospheric correction is comparatively minor (see Sect.A3).
The residual higher-order ionospheric effects are accounted for by a "conservative best-guess" value (0.05 µrad, reflecting results of Liu et al., 2015, and Danzer et al., 2013, 2015) and added (in root-mean-square form) to the systematic uncertainty profile u s αr , leading to a total estimated systematic uncertainty in this example case of ∼ 0.07 µrad (Fig. 8b).Within this uncertainty, the one dominating component from orbit uncertainties (∼ 0.05 µrad; cf.Fig. 6b) can be considered an apparent systematic uncertainty that will essentially average out in ensemble averaging (e.g., climatologies), while the other dominating component from residual higher-order ionospheric biases (also estimated ∼ 0.05 µrad as noted above) can be considered a basic systematic uncertainty.For the latter it is therefore useful and prepared for in the rOPS -in line with GUM recommendations and as discussed in the introductory Sect. 1 -to correct for the quantifiable part of it in the future so that the total basic systematic uncertainty may be mitigated down to the ∼ 0.01 µrad level.
The resolution profile w αr of the retrieved bending angle (Fig. 8f) is dominated by the contribution of α F1 that strongly dominates (intentionally by construction) the ionospheric correction results in terms of the small-scale bending angle variability.Similar to the correlation length profile l αr , it is therefore very close to w αF1 and only slightly larger.

Algorithm validation
The GUM advises to use a MC method for uncertainty propagation if the retrieval operators do not fulfill the criteria for a GUM-type CP.In our case the MC method is put to another beneficial use, to validate the results of the CP, as recommended by JCGM (2011).
For the validation of the covariance propagation by the MC method, we sampled the input excess phase profile random error distribution, statistically described by

C MC
Lr,ij = u r,Std Lr,i • u r,Std Lr,j • δ ij , by a large number M of draws L T r + r Lr,j (with j ∈ {1, . .., M} and M = 1000).For each of these M profile realizations, the state retrieval is run through the L1b retrieval chain, to give M realizations of the output variable X j (with X j ∈ {L F,kj (t), D r,kj (t), α G,kj (z a ), α F,kj (z a ), α r,j (z a )} and k ∈ {1, 2}).From these individual realizations the mean profiles X MC and the covariance matrices C MC X , are calculated (items b-g in Fig. 2).Using the same input profile and uncertainty information as used to specify the MC runs (described in Sect.2.2), the retrieval is then also run with covariance-based uncertainty propagation, and the resulting CP-propagated covariance matrices C CP X are compared to the MC-derived matrices C MC X .In order to be able to attribute potential changes between CP and MC covariance matrices better, we decompose C X into u r X and R X (Eqs.A7 and A8), and compare them separately.Figure 9a (light blue) and b (orange) show the random uncertainties u r Lr,1 and u r Lr,2 , respectively, which characterize the input distribution and from which the random error profiles r Lr,j are drawn.They also show the CP results for the random uncertainty u r LF1 (dark blue in Fig. 9a) and u r LF2 (red in Fig. 9b), compared to the MC propagated random uncertainties (black).
The CP and MC lines match very well and show that the implemented CP algorithm delivers correct results for the basic filtering step.For f T2 , the MC uncertainties do not reach down as far as the CP uncertainties, because the shortest of all draws of the large ensemble of size M determines how far down the recombined MC covariance matrix (Eq.10) reaches.Figure 9c compares CP correlation functions R LF,i1 (blue) and R LF,i2 (red) to the corresponding MC correlation functions (black dashed).
The CP and MC correlation functions also agree well.Both capture the narrow peak, broadened by the BWS filter.Again the MC correlation functions fluctuate around zero left and right of the peak, from the finite ensemble size, but it is obvious that the CP delivers the correct off-peak results (i.e., zero; the off-peak elements outside the BWS filter window must nominally be zero).The MC validation (black) of u r Dr1 (Fig. 9d), u r Dr2 (Fig. 9e), and R Dr,i (Fig. 9f, blue and red) demonstrates that the CP through the Doppler shift derivation performs correctly as well.
The next row, Fig. 9g to i, shows the results for the GO bending angle α G (z a ), i.e., after the interpolation of all quantities to the (common monotonic) impact altitude grid z a .For comparison, in Fig. 9a to f, all quantities have been computed on the common time grid ("setting time" relative to time zero at 80 km altitude) with 50 Hz sampling rate; the corresponding impact altitude of the "true" profile L T r is shown for additional convenience on the right-hand-side (RHS) axis.In Fig. 9g to o, these bending angle quantities have been computed on the impact altitude grid; in these cases therefore the corresponding setting time of the "true" profile is shown for additional convenience on the RHS axis.
The results for the filtered bending angle α F follow in Fig. 9j to l.Also here the MC results match the CP result well.Due to the lower BWS cutoff frequency for α F2 , now u r αF2 is smaller than u r αF1 , even though u r αG2 was larger than u r αG1 .Correspondingly the peak of the correlation functions R αF,i2 widened more than those of R αF,i1 (cf.Fig. 9l and i).
Finally, Fig. 9m to o show the CP results for retrieved atmospheric bending angle α r , where Fig. 9n is included as a special cross-comparison in the case of only variance propagation being used instead of CP. Figure 9m and o confirm that CP results are also correct for this final L1b variable, in terms both of random uncertainty and correlation functions.
In order to demonstrate that a full CP is necessary to propagate random uncertainties correctly, we also calculated random uncertainties u r αr based on mere VP from α G to α r for comparison.A description of this VP algorithm (i.e., only diagonal elements of the covariance matrices are considered) is provided in Appendix B. Figure 9n clearly shows that VP would overestimate random uncertainties in α r considerably, pointing to the importance of the complete CP implementation in the L1b retrieval chain, even though the correlation lengths involved in the processing steps are rather small.

Performance demonstration
To statistically evaluate the performance of the new L1b uncertainty propagation algorithm, we also processed a complete test day of real (CHAMP, COSMIC, MetOp) and simulated (simMetOp) data from GNSS RO satellite missions.Figure 10 shows the results for estimated systematic and ran-dom uncertainty profiles, as well as correlation length and resolution profiles for filtered excess phase profiles L F,1 .Figure 11 subsequently illustrates the ensemble mean of the same variables for L F1 , D r1 , α G1 , and α r for the test day ensemble.In Fig. 10 we also co-illustrate the number of events processed for each of the RO missions (middle column).
About 5 % of the total number of processed profiles for each mission have been discarded, because they were detected as outliers based on the magnitude of their random uncertainty profiles (these outliers are not included in the number of profiles shown).All profiles are shown as function of impact altitude, because each of the profiles in the ensembles needed to be interpolated to the same (standard) impact altitude grid, to orderly calculate their mean profiles.
Figure 10a shows u r LF1 and u s LF1 for all ∼ 100 CHAMP events.It is visible (also in Fig. 10d and g) that the random uncertainty is estimated based on excess phase noise between 30 and 75 km and synthetically extended above and below, as described in Sect.2.2.For the large majority of events, u r LF1 lies between about 0.5 and 3 mm in the range between 30 and 75 km.Note that these results show the random uncertainties after the application of the basic BWS filter (Sect.3.1.1),but the input uncertainties u r Lr1 are of similar shape (though larger in magnitude).
Figure 10b shows that the correlation length profiles of the CHAMP ensemble (gray) and its ensemble mean (yellow) are of relatively constant magnitude from 35 to 80 km but then get smaller downward, because the RO event's scan velocity decreases (see Eq. A13).Since the BWS filter determines the vertical resolution and the correlation length at the same time, the resolution profiles w LF1 (Fig. 10c) are quite similar to the correlation length profiles l LF1 (Fig. 10b).
The number-of-events profile shows that most CHAMP events end between 5 and 12 km (Fig. 10b, black).This is because the GO profiles illustrated here are cut off right at the lower end of the GO-WO transition range at z GW a − z GW a (cf.Table 2).
Compared to CHAMP, the mean random uncertainty u r LF1 (Fig. 10d) for the ∼ 1500 events of the COSMIC ensemble is smaller, particularly above 30 km, indicating the improved data quality of this later mission.The mean of the correlation length profiles l LF1 (Fig. 10e) is higher than for CHAMP (Fig. 10b), and correspondingly the resolution of the COS-MIC profiles is also somewhat coarser (Fig. 10f and c).The cutoff frequency and sampling rate -and thus the resolution in time -are set to be the same in the rOPS, irrespective of the missions; these differences hence are due to the different vertical scan velocities of the missions induced by the differences in orbit altitudes (CHAMP ∼ 400 km, COSMIC ∼ 700 km).
For the real MetOp data (available here as a data set from UCAR/CDAAC, as for CHAMP and COSMIC), u r LF1 appears similar to COSMIC (cf.Fig. 10d, g), while for simMetOp (with best possible simulated MetOp-type receiver noise) it is clearly smaller than for COSMIC.From 35  to 80 km the mean random uncertainty profile for simMetOp stays below 1 mm (Fig. 10j).Three individual profiles exhibit comparatively high uncertainties of larger than 2 mm within about 40 to 55 km, however, reflecting that the simMetOp error simulations are capable of partly generating higher-noise profiles of the type more frequently seen in the real MetOp data (Fig. 10g).
On the other hand, the average correlation length/resolution profile of the ∼ 500 real MetOp and ∼ 700 simMetOp ensemble members is very similar, driven by the orbit being essentially the same for the real data and the simulations (Fig. 10h, i, k, l).Compared to COSMIC (Fig. 10e, f), the correlation length and resolution are again somewhat larger/coarser, due to an even somewhat higher scan velocity of the MetOp satellite (∼ 820 km orbit altitude).The systematic uncertainty u s LF1 -just co-illustrated for completeness in Fig. 10a, d, g, and j -is almost left unchanged by the BWS filter and is essentially equal to the preset input uncertainty for all three missions (Sect.2.2).
Figure 11 shows how the u s LF1 , u r LF1 , l LF1 , and w LF1 profiles are on average affected by the uncertainty propagation.The color code for the different satellite missions is the same as in Fig. 10.The propagation effects visible are similar to those already seen in Figs. 3 to 8. The Doppler shift derivation increases the relative uncertainties and reduces correlation length (of the main peak), while the resolution stays the same (Fig. 11d-f).The GO bending angle retrieval leaves correlation length and resolution unchanged, while random uncertainties increase strongly in the lower stratosphere and troposphere due to the increasing refractive effects (Fig. 11gi).
Finally, the BWS filtering before the ionospheric correction decreases random uncertainties and increases correlation length, and resolution somewhat.However, the linear combination of the two bending angle profiles α F1 and α F2 then increases the random uncertainty again (cf.Fig. 11j and g).The adaptive minor-channel cutoff frequency f c2 for the relatively noisy CHAMP profiles is generally lower than for the other two missions, and the filter effect is therefore stronger for CHAMP (indicated by the larger l αr in Fig. 11k) The estimated systematic uncertainty of the atmospheric bending angle u s αr , indicated for completeness in Fig. 11 (left column, inflated by a factor of 10 in panel a and 100 in panels d, g, and j for somewhat better visibility), stays below 0.1 µrad for all three missions.

Conclusions
In order to deliver climate benchmark data sets, it is essential to integrate uncertainty propagation in RO retrievals.In this study we presented the uncertainty propagation algorithm chain from excess phase profiles to atmospheric bending angle profiles (L1b processing), as newly implemented in the rOPS at the WEGC.Along with the basic profile retrieval, we provide estimates for systematic and random uncertainties, error correlation matrices, and vertical resolution profiles, which is unique amongst all existing RO processing systems so far (Ho et al., 2012;Steiner et al., 2013).
We validated the implemented algorithm via comparison to Monte Carlo sample propagation results and demonstrated the performance of the algorithm using real-data ensembles.We find close agreement between the implemented covariance propagation of random uncertainties and the Monte Carlo validation runs, verifying the correctness of the implemented algorithm.The test day ensembles for three different missions (CHAMP, COSMIC, MetOp) show reliable, robust, and consistent results that provide valuable insight and understanding of retrieval chain details.
Together with the integration of the uncertainty propagation algorithm from atmospheric bending angle profiles to dry-air profiles (L2a processing) presented by Schwarz et al. (2017), the rOPS can now provide estimates of systematic and random uncertainty profiles, of error correlation matrices and resolution, and of observation-to-background weighting ratio profiles from excess phase to dry-air profiles.
The next step towards the final atmospheric profiles, currently ongoing, is the introduction of integrated uncertainty propagation for the moist-air retrieval (L2b processing).Implementation of uncertainty propagation for the wave-optics bending angle retrieval and for the orbit determination and excess phase processing (L1a processing) is ongoing as well.
Once completed, the full rOPS retrieval chain will run with integrated uncertainty estimation, a major step towards climate benchmark data provision, and beneficial for the wide diversity of uses in atmospheric and climate science and applications.
Data availability.The RO excess phase and orbit data used in the study are available from UCAR/CDAAC Boulder, CO, USA, at http://cdaac-www.cosmic.ucar.edu/(last access: 26 March 2018).The ECMWF analysis and forecast data are accessible from ECMWF Reading, UK, at http://www.ecmwf.int/en/forecasts/datasets (last access: 26 March 2018).To access the relevant result files of the uncertainty propagation please contact the corresponding author.The developed algorithm is provided in the Appendix.Fm ", green) and M = 11 points ("δL BC11 Fm ", blue).
displayed in Fig. A1a, while Fig. A1b compares the filter results.
It is clearly seen that the smoothing window width of the BWS filter best corresponds to an 11-point boxcar filter (confirmed numerically by minimization of the sum of squared differences between boxcar and BWS filter result), while giving considerably better filtering results (as for example visible between 31.5 and 32.0 s, where the 11-point boxcar filter zigzags around the BWS result).The effective filter width of the BWS filter, which we also term "boxcarequivalent width", is therefore its full width at half maximum (see Fig. A1a), corresponding to M/4 + 1 samples with our design.
The actually used sample width M of the BWS filter is equal to M, except that it decreases at the top and bottom of the profile such that it does not reach beyond the first/last element of the vector to be filtered.At the ith grid point (with i ∈ {1, 2, . .., N }, and N being the profile length in grid points), the filter width M is thus The state profile of the filtered phase L F is obtained using the "baseband approach" (Kirchengast et al., 2016a), i.e., by first subtracting a zero-order model profile L m and applying the filter only to the delta profile δL rm = L r − L m (with the model profile being adequately smooth over the scale of the filter window width).This approach efficiently mitigates residual numerical biases.After the application of the BWS filter, the model profile is added back again.We express the BWS filter as a linear matrix operator A BWS and get (item 1.2 in Fig. 2) for the filtered excess phase, where j ∈ {1, 2, . .., N }.The band matrix operator A BWS has elements The central filter weight w 0+M/2 at j = i is the (M/2)th filter element (according to the definition of the BWS weights below); therefore its index is M/2.With m = j − i + M/2 (and therefore 0 ≤ m ≤ M), each single BWS weight is calculated using and The estimated random uncertainty is then propagated by covariance propagation (item 1.3 in Fig. 2), The random uncertainty profile u r LF and the error correlation matrix R LF are not needed for the subsequent random uncertainty propagation but are calculated from C LF for being available for the L1b output, using and The correlation length profile l Lr has elements computed upward and downward from the main peak of the correlation function and then averaged.Here dz/dt is the scan velocity profile, obtained from using the MSL altitude grid z t calculated as part of the forward modeling towards L m at the corresponding time grid t (cf.Table 2).We note that after the L2a refractivity retrieval also the MSL altitude grid consistent with the retrieved refractivity profile could be used (as described by SKS2017, Appendix A therein), from a repeated forward modeling.The difference for the scan velocity estimate is found to be very small, however, since the forward-modeled z t based on collocated refractivity profiles from ECMWF short-range forecast fields is already sufficiently reliable, and this also keeps the L1b processor as a decoupled predecessor of the L2a processor.
For the estimated systematic uncertainty, interpreted as a basic systematic uncertainty (Sect.2.2), we apply the same low-pass filter as used for the state profile (item 1.2 in Fig. 2), but with no zero-order profile subtracted, i.e., The resolution in time of L F and its uncertainties, τ BW , is the boxcar-equivalent width (cf.Fig. A1a) determined by the cutoff frequency f c of the BWS filter, with our design choice M = 2(f s /f c ) and with the BWS filter stopband-to-passband transition width being (Smith, 1999) Given f c = 2.5 Hz, this results in an effective resolution τ LF = 0.2 s and corresponds to the resolution obtained when applying a 11-point boxcar filter as explained at the beginning of this section above.The filter window intercomparison in Fig. A1a also illustrates this, because the full width at half maximum of the 2.5 Hz 41-point BWS filter is 11 points.
This resolution in time can finally be converted to the vertical (MSL altitude) resolution in space: where, as for the correlation length estimation (Eq.A9), the scan velocity profile is employed to convert from the time domain to MSL altitude domain.

A1.2 Doppler shift derivation
After the application of the BWS filter to the excess phase profiles L r (for both carrier frequencies of the given GNSS system), the state profile of the Doppler is derived from the filtered phase profiles L F (item 1.4 in Fig. 2).To minimize systematic errors from the numerical differentiation to negligible magnitude, the model profile L m is again subtracted from the filtered phase profile, and the resulting delta profile δL Fm is then differentiated.After the derivative, the zero-order Doppler shift model profile D m is added (the latter also available from the forward modeling, in a form strictly consistent with the excess phase model profile L m ).
Based on careful tests of different formulations, we use a five-point derivative scheme.The discretization of this fivepoint derivative δD rm,i is given by for each of the frequencies (e.g., Syndergaard, 1999).This can be expressed in matrix form as using matrix operator A L2D with , where t = t i+1 − t i , being 0.02 s in our case of f s = 50 Hz.
The estimated random uncertainty can then be propagated (item 1.5 in Fig. 2) using The covariance matrix is again (cf.Eqs.A7 and A8) decomposed into estimated random uncertainties and error correlation matrix (item 2.2 in Fig. 2 For the estimated systematic uncertainty, further on interpreted as basic systematic uncertainty (cf.Eq.A10), we apply the derivative operator (item 1.4 in Fig. 2) to the systematic uncertainties, with no zero-order profile subtracted, i.e., The resolution remains unaffected by the Doppler shift derivation, since the five-point sample width of the derivative operator is fully within the 11-point effective filter width (stopband) of the BWS filter applied before, so that τ Dr = τ LF and w Dr = w LF .

A2 Bending angle retrieval A2.1 GO bending angle retrieval
From the Doppler shift state profile D r (again for both frequencies of the given GNSS system) we can derive the impact parameter profile a t and GO bending angle profile α G (item 2.1 in Fig. 2) using first the geometric relation where and for each individual level of the time grid t i , in order to determine a t from sequential application to all levels (Kursinski et al., 1997;Syndergaard, 1999).Here v R,i = |v R,i | is the receiver velocity; r R,i = |r R,i | the receiver radial position; η R,i the angle between the receiver velocity and position vectors; φ R,i then the angle between the receiver velocity and ray path vectors (and all these equivalently for the transmitter); is the time derivative of the distance between the transmitter and the receiver at time t i , i.e., the "kinematic straight-line Doppler shift" to be subtracted in Eq. (A22) to match the (excess) Doppler shift D r,i induced by the atmosphere (and ionosphere).
Based on a t , the elements of the GO bending angle profile α G are subsequently calculated using another geometrical relation: where θ RT,i is the opening angle between the transmitter and receiver position vectors.Syndergaard (1999, Fig. 1.5 therein) provides an illustration of the relevant geometry.All the variables in Eqs.(A22)-(A25) are defined in the occultation plane spanned by the receiver and transmitter position vectors after oblateness correction (Syndergaard, 1998), i.e., after they have been transformed to originate in the Earth ellipsoid's center of local curvature in the occultation plane at the mean tangent point (MTP) location of the RO event.
The MTP location is defined as the geodetic (geographic) location on the WGS84 ellipsoid, where the straight-line path between transmitter and receiver touches this ellipsoid, i.e., where the straight-line tangent height is zero.This can be computed with very high accuracy at the sub-meter level (see Scherllin-Pirscher et al., 2017, for more details on the geolocation accuracy of RO).Using the MTP location's center of local curvature rather than the Earth's center of mass as the origin is essential to ensure that the assumption of spherical symmetry, implicit in Eqs.(A22) to (A25), is accurately valid geometrically.
The impact parameter retrieval is solved iteratively, because it is impossible to rearrange Eqs.(A22) to (A24) into an explicit expression for the retrieval of the impact parameter, but it is mildly nonlinear and converges fast, in particular if the initial guess for a t,i is estimated from the previous level (starting at the top level with the straight-line impact parameter).
After the GO bending angle retrieval, the bending angles of all GNSS frequencies are interpolated to a common monotonic impact altitude grid z a (item 2.6 in Fig. 2), based on the monotonically sorted impact parameter grid of the leading channel, a t1 (i.e., k = 1).
For each element of z a we get (item 2.3 in Fig. 2) where j is the index of the elements of the sorted impact parameter grid a t1 .h G is the geoid undulation (see Scherllin-Pirscher et al., 2017, for a detailed discussion of its use in RO analysis), and R C is the local radius of curvature of the RO event.
Because the impact parameter is only implicitly expressed in Eqs.(A22)-( A24), but GUM-type uncertainty propagation along Eqs. ( 2) and ( 4) requires an explicit measurement model, we make use of a linearization of the bending angle retrieval.We use the approach described by Melbourne et al. (1994), and applied to uncertainty propagation by Syndergaard (1999), for the propagation of the estimated random uncertainty from Doppler shift D r to GO bending angle α G (item 2.5 in Fig. 2).
This linearization establishes a direct relation between random uncertainties of the Doppler shift u r Dr and the uncertainties of the bending angle u r αG , using where a SL is the straight-line impact parameter.These bending angle uncertainties u r αG are relative to the time grid as independent coordinate.To get the desired uncertainties with respect to the impact altitude grid z a (introduced in Eq.A26), the uncertainties of the impact altitude z a need to be transferred to the bending angle, so that the z a grid can subsequently be considered free of error.Syndergaard (1999) showed that this can be done by replacing Eq. (A27) with where a T t is the "true" impact parameter.We use the forwardmodeled impact parameter a tm instead (i.e., adopt a tm = a T t ) and accept the additional error thus incurred, assuming it is smaller than the 2 % relative error due to the linearization estimated by Melbourne et al. (1994).This is a reasonable assumption given the high quality of our forward-modeled profiles derived from ECMWF short-range forecast refractivity fields.
As a consequence we have to accept that the overall inaccuracy of our random uncertainty estimate cannot be brought below 2 %.Therefore, to ensure that our simplified estimate does not underestimate the real uncertainty, we account for the linearization error by multiplying a factor f uαlin = 1.02 to the uncertainty of the retrieved GO bending angle In this way we acknowledge that although the calculation of the state of the bending angle does not make use of the linearization, and therefore the linearization does not increase the uncertainty of the state profile, it may increase the error in the uncertainty estimate itself.Finally, the u r αG profile is also interpolated to the common monotonic impact altitude grid z a .
In the GO approximation, the bending angle values at each grid point only depend on the Doppler shift values of the same grid points; i.e., the existing correlations between the errors at different levels are left unchanged: i.e., R αG = R Dr .The covariance matrix can hence be calculated by recombining the Doppler shift correlation matrix with the propagated www.atmos-meas-tech.net/11/2601/2018/uncertainties (item 2.8 in Fig. 2), For the propagation of the estimated systematic uncertainty (item 2.4 in Fig. 2) three types of potential systematic errors adding to the impact parameter uncertainty u s at , and consequentially the bending angle uncertainty u s αG , are taken into account: systematic errors in the Doppler shift, i.e., u s Dr ; systematic errors in the velocities of the satellites, i.e., u s vR and u s vT ; and systematic errors in the positions of the satellites, i.e., u s rR and u s rT .The latter two orbit-borne types are interpreted as apparent systematic uncertainties (Sect.2.2), while the excess phase-borne uncertainty u s Dr is a basic systematic uncertainty.
For the propagation of these estimated systematic uncertainties to u s αG , Eqs. (A22)-(A24) are linearized around the retrieved state quantities (serving as zero-order state), and no terms higher than first-order are kept.Then φ R and φ T in Eq. (A22) are substituted by the linearized versions of Eqs.(A23) and (A24), and the resulting equation is solved (level by level) for the impact parameter a t,i = f (D r,i , r R,i , r T,i , v R,i , v T,i ), with i = 1, 2, . .., N .Adopting the first-order deviations to represent the estimated systematic uncertainties, we obtain where A number of simplifications have been made to arrive at this result.First, the last term in Eq.A22 is disregarded since errors in the positions are assumed to be constant with respect to the short time duration of an RO event; remaining errors ṙRT after taking the derivative are therefore of higher order.Next, orbit position and velocity uncertainties are both assumed to be constant within the short duration of an event, and the velocity uncertainties obtained are interpreted as uncertainties along the direction of the velocity vector.Consequentially, the uncertainty is also projected along with the vector into the ray path direction.A more conservative estimation (which we consider overly conservative in context) would interpret the uncertainties as ellipsoids at the velocity vectors' heads and would hence take the full magnitude of the uncertainties along the ray path direction (without projection).
Furthermore, since all error sources (the processing of the occultation tracking data and the POD for transmitter and receiver) are essentially independent from each other, the different input uncertainties are assumed to be uncorrelated.Finally, we reasonably assumed the errors of the angle between the position and velocity vectors (η) to be negligible (u s η ≈ 0) for the purpose here, for both the transmitter and receiver.
In order to finally derive the systematic uncertainty of the bending angle from the impact parameter's uncertainty, we continue with a linearization of Eq. ( A25) and arrive at u s αG,i = (A33) .
In practice we separately calculate the basic and apparent systematic uncertainty estimates (u b αG from the first RHS terms in Eqs.A31 and A33, u a αG from the orbit-borne terms) and afterwards obtain u s αG as a combined result, in order to enable separate propagation in subsequent processing steps.
The resolution profile remains unaffected by the bending angle retrieval, since the level-by-level approach of the algorithm does not create extra correlation and further vertical smoothing, so that τ αG = τ Dr and w αG = w Dr .

A2.2 WO bending angle retrieval
After the GO bending angle, the WO bending angle state profile α W (z a ) is retrieved (item 2.7 in Fig. 2) from excess phase profile L r (t) (and its uncertainties) and the amplitude profile A r (t) (and uncertainties) in a WO retrieval following Gorbunov and Kirchengast (2015Kirchengast ( , 2018)).Along with the state profile, the systematic uncertainty profile u s αW , the covariance matrix C αW , and the resolution profile w αW are derived.
The covariance matrix C αW is then decomposed to random uncertainty profile u r αW and correlation matrix R αW in the same form as done above for C Dr (Eqs.A19 and A20) and C LF (Eqs.A7 and A8).The estimated systematic uncertainty u s αW is composed of a basic systematic uncertainty u b αW , propagated through the wave-optical retrieval from the excess phase uncertainty u s Lr , and an apparent systematic uncertainty u a αW , estimated in the lower troposphere as residual bias uncertainty of a regression-based boundary layer bias correction (Gorbunov and Kirchengast, 2018).
The WO bending angle retrieval algorithm and the associated uncertainty propagation algorithm are not explicitly described here; the reader is referred to Gorbunov and Kirchengast (2015) and Gorbunov and Kirchengast (2018).However, we have prepared the merging with the WO bending angle variables (they will be actually merged in when the WO tests within the rOPS are complete), which is described next.

A2.3 Merging of GO and WO bending angle profiles
The α W profile, prepared on the common grid z a , and the α G profile are merged over an upper-tropospheric transition range (item 2.9 in Fig. 2).The gradual transition, weighted by a symmetric half-sine function, has a defined impact altitude transition of half-width z GW a = 2 km around transition altitude z GW a , allowed within 9 km to 14 km, estimated from α G data quality.The resulting merged bending angle profile α M is where the weighting profile γ is formulated as (A36) To determine the random uncertainties for the merged GO-WO input bending angle, we need to merge the covariance matrices of both bending angles.
We can assume both incoming covariance matrices C αG and C αW are provided on the common monotonic target grid z a (i.e., also the WO uncertainties and correlations are interpolated to this common grid before the merger).We further can reasonably assume that there are no cross-correlations between GO and WO errors, given the very different retrieval schemes.Based on this we can compose the covariance matrix of the merged bending angle profile, C αM (item 2.10 in Fig. 2) as follows.Outside the merging zone (i.e., outside of z GW a ± z GW a ) we can assign Due to the linear relation between α M , α G , and α W , expressed by Eq. (A35), a bias u s αG in the GO bending angle and a bias u s αW in the WO bending angle can be combined linearily as well, and we can compute the estimated systematic uncertainty of the merged bending angle u s αM according to (item 2.9 in Fig. 2) In practice this formulation is again applied separately for the basic and apparent systematic uncertainty estimates, afterwards obtaining the u s αM profile as a combined result, in order to allow separate propagation in subsequent processing steps.
The resolution profile of the bending angle, w αM , is equal to the GO resolution w αG above z GW a + z GW a , equal to the WO resolution w αW below z GW a − z GW a and has a transition with transition weight γ i in between, again following the linear formulation such as in Eqs.(A35) and (A39).
Because the integration and testing of the uncertainty propagation through the rOPS WO bending angle retrieval are currently still ongoing, as noted in Sect.A2.2 above, the examples shown in this study are all GO-only; i.e., only the GO retrieval is performed.The merging algorithm as described is ready to include the WO bending angles, however.

A3 Atmospheric bending angle derivation
In order to retrieve the atmospheric bending angle profile α r , ionospheric effects need to be corrected for, using the retrieved bending angles from each transmitter frequency channel.Since the only GNSS constellation currently used for RO is the GPS -except for recent initial data from the Chinese GNOS instrument using BeiDou signals (Liao et al., 2016;Bai et al., 2018) -the data characteristics of the Figure1.Schematic view of the main processors of the retrieval chain in the rOPS (L1a, L1b highlighted, L2a, L2b) and the main operators of the L1b processor (1, 2, 3), which are in the focus of this study.
Interpolation of all quantities to common monotonic impact altitude grid (2.3) Grid generation (2.4) Syst.unc.propagation (2.5) Rand.unc.propagation and R L r,k in terms of representative correlation functions in panels (d) and (e), for a typical COSMIC RO event of the test day ensemble from 15 July 2008 (example case).

Figure 3 .
Figure 3. Input profiles of retrieved excess phase L r (with model profile L m for comparison) in panel (a), estimated systematic uncertainty profiles u s Lr in panel (b), estimated random uncertainty profiles u r Lr in panel (c), representative correlation functions R Lr,i (at 10, 30, 50, and 70 km) in panels (d) and (e), and correlation length l Lr (solid) and resolution profiles w Lr (dotted) in panel (f), which are set zero for these initial essentially uncorrelated input data.All profiles are shown for both GPS carrier frequencies: f T1 (blue) and f T2 (red).

Figure 4 .
Figure 4. Results for filtered excess phase profiles L F (with model profile L m for comparison) in panel (a), estimated systematic uncertainty profiles u s LF in panel (b), estimated random uncertainty profiles u r LF in panel (c), representative correlation functions R LF,i (at 10, 30, 50, and 70 km) in panels (d) and (e), and correlation length l LF (solid) and resolution profiles w LF (dotted) in panel (f).All profiles are shown for both GPS carrier frequencies: f T1 (blue) and f T2 (red).

Figure 5 .
Figure 5. Results for retrieved Doppler shift profiles D r (with model profile D m for comparison) in panel (a), estimated systematic uncertainty profiles u s Dr in panel (b), estimated random uncertainty profiles u r Dr in panel (c), representative correlation functions R Dr,i (at 10, 30, 50, and 70 km) in panels (d) and (e), and correlation length l Dr (solid) and resolution profiles w Dr (dotted, estimated for main peak) in panel (f).All profiles are shown for both GPS carrier frequencies: f T1 (blue) and f T2 (red).

Figure 6 .
Figure 6.Results for geometric optics bending angle profiles α G (with model profile α m for comparison) in panel (a), estimated systematic uncertainty profiles u s αG in panel (b), estimated random uncertainty profiles u r αG in panel (c), representative correlation functions R αG,i (at 10, 30, 50, and 70 km) in panels (d) and (e), and correlation length l αG (solid) and resolution profile w αG (dotted, estimated for main peak) in panel (f).All profiles are shown for both GPS carrier frequencies: f T1 (blue) and f T2 (red); in panels (b) and (f) both profiles are essentially identical (so that blue shadows the red color).

Figure 7 .
Figure 7. Results for filtered bending angle profiles α F (with model profile α m for comparison) in panel (a), systematic uncertainty profiles u s αF in panel (b), random uncertainty profiles u r αF in panel (c), representative correlation functions R αF,i (at 10, 30, 50, and 70 km) in panels (d) and (e), and correlation length l αF (solid) and resolution profiles w αF (dotted, estimated for main peak) in panel (f).All profiles are shown for both GPS carrier frequencies: f T1 (blue) and f T2 (red).

Figure 8 .
Figure 8. Results for atmospheric bending angle profile α r (with model profile α m for comparison) in panel (a), systematic uncertainty profile u s αr in panel (b), random uncertainty profile u r αr in panel (c), representative correlation functions R αr,i (at 10, 30, 50, and 70 km) in panel (d), and correlation length l αr (solid) and resolution profile w αr (dotted, estimated for main peak) in panel (f).

Figure 9
Figure 9 shows the different steps along the retrieval chain from L F,k (t) to D r,k (t), α G,k (z a ), α F,k (z a ), and α r (z a ) in the rows, for k = 1 (GPS f T1 frequency) in the left column and for k = 2 (GPS f T2 frequency) in the middle column.The right column shows multiple representative correlation functions, from near 10 to near 70 km.Due to the limited number of MC draws, the MC results (black lines) show some jitter both in the estimated random uncertainty and in the correlation functions.Since the purpose of the MC results is only to demonstrate the correctness of the CP result, we can disregard this behavior.Figure9a(light blue) and b (orange) show the random uncertainties u rLr,1 and u r Lr,2 , respectively, which characterize the input distribution and from which the random error profiles r Lr,j are drawn.They also show the CP results for the random uncertainty u r LF1 (dark blue in Fig.9a) and u r LF2 (red in Fig.9b), compared to the MC propagated random uncertainties (black).The CP and MC lines match very well and show that the implemented CP algorithm delivers correct results for the basic filtering step.For f T2 , the MC uncertainties do not reach down as far as the CP uncertainties, because the shortest of all draws of the large ensemble of size M determines how far down the recombined MC covariance matrix (Eq.10) reaches.Figure9ccompares CP correlation functions R LF,i1

Figure 9 .
Figure 9. Results from the validation of CP covariance matrices C CP Lr ("CP") by MC covariance matrices C MC Lr ("MC") (first four rows): the consecutive retrieval steps are shown for L F (a-c; in panels a-b also for L r ) and D r (d-f) relative to setting time t, and for α G (g-i) and α F (j-l) relative to impact altitude z a .The left column shows estimated random uncertainties for f T1 (CP in blue, MC in black, in panel a u r Lr1 in light blue); the middle column for f T2 (CP in red, MC in black, in panel b u r Lr2 in orange); the right column representative correlation functions at 60, 36, and 12 km for f T1 (CP in blue, MC in black) and 72, 48, and 24 km for f T2 (CP in red, MC in black).The last row (m-o) shows CP (blue) and MC (black) results for estimated α r random uncertainties (m) and representative correlation functions at 72, 60, 48, 36, 24, and 13 km (o), as well as variance propagation ("VP") results (light blue) for α r in panel (n).

Figure 10 .
Figure 10.Uncertainty propagation results for real-data ensembles from 15 July 2008, for the filtered excess phase profile L F1 of the leading channel (f T1 , GPS L1 frequency).Left column: estimated random u r LF1 (heavy) and systematic u s LF1 (light) uncertainty profiles of each ensemble member (gray), and the ensemble mean (color) for CHAMP (a), COSMIC (d), MetOp (g), and simMetOp (j).Middle column: correlation length profiles l LF1 of each ensemble member (gray), the ensemble mean (color), and the ensemble size profile (black, scale at upper axis) for CHAMP (b), COSMIC (e), MetOp (h), and simMetOp (k).Right column: estimated resolution profile w LF1 of each ensemble member (gray) and the ensemble mean (color) for CHAMP (c), COSMIC (f), MetOp (i), and simMetOp (l).

Figure 11 .
Figure 11.Uncertainty propagation results for real-data ensembles from 15 July 2008 for output profiles of the leading channel (f T1 , GPS L1 frequency).The first row shows results for L F1 (a-c), the second for D r1 (d-f), the third for α G1 (g-i), and the fourth for α r (j-l).The different ensemble mean profiles are shown in colors (CHAMP, yellow; COSMIC, orange; MetOp, red; and simMetOp, violet).Left column: mean random uncertainty u r Xr (heavy) and mean systematic uncertainty u s Xr (light) profiles (a, d, g, j); the latter is shown as 10 × u s LF1 (in a) and 100 × u s Xr (d, g, j) for enabling visibility of these small quantities.Middle column: correlation length profiles l Xr (b, e, h, k).Right column: vertical resolution profiles w Xr (c, f, i, l).

Figure A1 .
Figure A1.Comparison of the Blackman windowed sinc (BWS) low-pass filter and boxcar (BC) filters based on a representative segment (between 30.3 and 32.7 s) of the excess phase profile L r1 of the COSMIC example event.(a) Filter functions for the BWS filter with f c = 2.5 Hz and M = 41 points ("BWS", red) and boxcar filters with M = 21 points ("BC21", green) and with M = 11 points ("BC11", blue), around the central value of the segment (31.5 s).(b) Filter effects on the excess phase profile L r1 from running the filters over the segment.Shown are the unfiltered excess phase delta profile ("δL rm ", light gray), the BWS filtered profile with f c = 2.5 Hz and M = 41 points ("δL BWS Fm ", red), and the boxcar filtered profiles with M = 21 points ("δL BC21Fm ", green) and M = 11 points ("δL BC11Fm ", blue).
,i > z aBot and z aTop > z a,j > z GW a merging zone we can assignC αM,ij = γ i γ j C αG,ij + (1 − γ i )(1 − γ j )C αW,ij , that z aTop > z a,j > z aBot .Because of the symmetry of the covariance matrix, the covariance elements in the merging zone orthogonal to the one above, i.e., for z GW a + z GW a > z a,j > z GW a − z GW a and z aTop > z a,i > z aBot , are calculated according to the same formula.

Table 1 .
Principal variables for the rOPS L1b uncertainty propagation.retrieved excess phase/filtered excess phase/retrieved Doppler shift/retrieved geometric optics bending angle/merged GO WO bending angle/filtered bending angle/retrieved bending angle, with X r

Table 2 .
Vertical grids, coordinate variables, and specific settings for the rOPS L1b processing system.Tk (for GPS transmitters k ∈ {1, 2} denoting the L-band frequencies f T1 = 1.57542GHz and f T2 = 1.22760GHz).altitudegrid,calculated from sorted impact parameters a t,i of the leading channel (f T1 ) bending angle, via z a,i = a t,i − h G − R C .Used as standard vertical grid after interpolation of all dependent quantities to z a .ztm MSL altitude grid corresponding to time grid t, obtained as part of the forward modeling towards α m , D m , and L m (cf.Table1).
m impact altitude of the bottom of the RO profile; it can lie between 25 km and the Earth's surface.Its value can be different for the different GNSS frequencies (i.e., z aBot,k , for k ∈ {1, 2}).