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 from 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 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 basic 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 CHAllenging Minisatellite Payload (CHAMP); Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC); and Meteorological Operational Satellite A (MetOp). 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 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.

1 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-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), and Steiner 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.

Figure 1 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.

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 sequential 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).

Figure 2 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.

Table 1 Principal variables for the rOPS L1b uncertainty propagation.

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.

Table 2 Vertical grids, coordinate variables, and specific settings for the rOPS L1b processing system.

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 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 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. (2013, 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; and Danzer 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 (COSMIC) (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; simMetOp 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

2.1 Methods

We follow the Guide to the Expression of Uncertainty in Measurement (JCGM, 2008a, 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 noise-resolving 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-range-correlated 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

$$\begin{array}{}\text{(1)}& Y={\mathbf{A}}^{XY}\cdot X,\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }{\mathbf{C}}_Y& {\displaystyle }=E\left[Y{Y}^T\right]=E\left[\left({\mathbf{A}}^{XY}X\right){\left({\mathbf{A}}^{XY}X\right)}^T\right]\\ \text{(2)}& {\displaystyle }& {\displaystyle }={\mathbf{A}}^{XY}E\left[X{X}^T\right]{\left({\mathbf{A}}^{XY}\right)}^T={\mathbf{A}}^{XY}{\mathbf{C}}_X{\left({\mathbf{A}}^{XY}\right)}^T\end{array}$$

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_X^r$ and the correlation matrix R_X:

$$\begin{array}{}\text{(3)}& C_{X,ij}=u_{X,i}^r\cdot u_{X,j}^r\cdot R_{X,ij}.\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }{u}_Y^s& {\displaystyle }={\mathbf{A}}^{XY}\cdot \left(X+u_X^s\right)-Y={\mathbf{A}}^{XY}\cdot X+{\mathbf{A}}^{XY}\cdot u_X^s-Y\\ \text{(4)}& {\displaystyle }& {\displaystyle }=Y+{\mathbf{A}}^{XY}\cdot u_X^s-Y={\mathbf{A}}^{XY}\cdot u_X^s,\end{array}$$

where $u_X^s$ and $u_Y^s$ 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).

2.2 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_{L_{\mathrm{r},k}}^s\left(t\right)$, random uncertainty profiles $u_{L_{\mathrm{r},k}}^r\left(t\right)$, and correlation matrices ${\mathbf{R}}_{L_{\mathrm{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_{\mathit{r}\mathrm{R}}^s\left(t\right)$, $u_{\mathit{v}\mathrm{R}}^s\left(t\right)$, $u_{\mathit{r}\mathrm{T}}^s\left(t\right)$, and $u_{\mathit{v}\mathrm{T}}^s\left(t\right)$. 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_{L_{\mathrm{r},k}}^s\left(t\right)$. 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\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$, f_T1=1.57542 GHz, and f_T2=1.22760 GHz.

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_{L\mathrm{r}}^s$ in panel (b), estimated random uncertainty profiles $u_{L\mathrm{r}}^r$ 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).

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 L_r,k(t) in panel (a), $u_{L_{\mathrm{r},k}}^s\left(t\right)$ in panel (b), $u_{L_{\mathrm{r},k}}^r\left(t\right)$ in panel (c), and ${\mathbf{R}}_{L_{\mathrm{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).

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_{L_{\mathrm{r},k}}^r\left(t\right)$ 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_{\stackrel{\mathrm{̃}}{\mathrm{m}}}$. Next, the delta profile $\mathit{\delta }{L}_{\mathrm{r}\stackrel{\mathrm{̃}}{\mathrm{m}},k}=L_{\mathrm{r},k}-L_{\stackrel{\mathrm{̃}}{\mathrm{m}}}$ was calculated. After smoothing $\mathit{\delta }{L}_{\mathrm{r}\stackrel{\mathrm{̃}}{\mathrm{m}},k}$ with a 10 km moving-average boxcar (BC) filter, the smoothed profile was subtracted from $\mathit{\delta }{L}_{\mathrm{r}\stackrel{\mathrm{̃}}{\mathrm{m}},k}$ again, to get $\mathit{\delta }\mathit{\delta }{L}_{\mathrm{r}\stackrel{\mathrm{̃}}{\mathrm{m}},k}$, the random noise profile component of L_r,k isolated in this way. Finally, the estimated random uncertainty was determined as

$$\begin{array}{}\text{(5)}& u_{L\mathrm{r},ik}^r=\sqrt{\sum _{j=i-M/\mathrm{2}}^{i+M/\mathrm{2}}\mathit{\delta }\mathit{\delta }{L}_{\mathrm{r}\stackrel{\mathrm{̃}}{\mathrm{m}},jk}^{\mathrm{2}}},\end{array}$$

where M is the number of grid points equivalent to a window width of 10 km. To avoid boundary effects of the filter, $u_{L_{\mathrm{r},k}}^r$ 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])

$$\begin{array}{}\text{(6)}& u_{L\mathrm{r},ik}^r=u_{L_{\mathrm{r},k}}^r\left(z_{\mathrm{aGradr}}\right)+{\displaystyle \frac{z_{\mathrm{aGradr}}-z_{\mathrm{a},i}}{\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}}},\end{array}$$

for all elements of $u_{L_{\mathrm{r},k}}^r\left(t\right)$ 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)

$$\begin{array}{}\text{(7)}& C_{L\mathrm{r},ijk}=u_{L\mathrm{r},ik}^r\cdot u_{L\mathrm{r},jk}^r\cdot R_{L\mathrm{r},ijk}=u_{L\mathrm{r},ik}^r\cdot u_{L\mathrm{r},jk}^r\cdot {\mathit{\delta }}_{ij}.\end{array}$$

For the MC validation of the CP, error profile realizations ${\mathit{ϵ}}_{L\mathrm{r}}^r$ were superimposed onto simulated “true” excess phase profiles $L_{\mathrm{r},k}^{\mathrm{T}}\left(t\right)$. As a source for $L_{\mathrm{r},k}^{\mathrm{T}}$, 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_{L\mathrm{r}}^{r,\mathrm{Std}}$ 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_{L\mathrm{r},ij}={\mathit{\delta }}_{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_{L\mathrm{r}}^{r,\mathrm{Std}}$ 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_{L\mathrm{r},k}^s$ 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_{L\mathrm{r},k}^s$ 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_{L\mathrm{r},\mathrm{1}}^s=\mathrm{0.2}$ mm and $u_{L\mathrm{r},\mathrm{2}}^s=\mathrm{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_{L\mathrm{r},k}^s$ (in units [m]) increases by

$$\begin{array}{}\text{(8)}& u_{L\mathrm{r},ik}^s=u_{L_{\mathrm{r},k}}^s\left(z_{\mathrm{aGrads}}\right)+{\displaystyle \frac{z_{\mathrm{aGrads}}-z_{\mathrm{a},i}}{\mathrm{3}\times {\mathrm{10}}^{\mathrm{7}}}}.\end{array}$$

In order to avoid a sharp kink in the $u_{L_{\mathrm{r},k}}^r$ profiles at z_aGradr, and in the $u_{L_{\mathrm{r},k}}^s$ profiles at z_aGrads, a 2 km width moving-average boxcar filter was applied to smooth these simple uncertainty models around these transition altitudes (for the example profile $u_{L\mathrm{r}}^s$ 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_{\mathit{r}\mathrm{T}}^s=\mathrm{3}$ cm and $u_{\mathit{v}\mathrm{T}}^s=\mathrm{0.01}$ mm s⁻¹, consistent with accuracies for GPS orbits available from GNSS orbit providers like the International GNSS Service (IGS). The receiver position and velocity uncertainties, $u_{\mathit{r}\mathrm{R}}^s=\mathrm{5}$ cm and $u_{\mathit{v}\mathrm{R}}^s=\mathrm{0.05}$ mm s⁻¹ for CHAMP and MetOp, are adopted 4 times smaller than those for COSMIC with $u_{\mathit{r}\mathrm{R}}^s=\mathrm{20}$ cm and $u_{\mathit{v}\mathrm{R}}^s=\mathrm{0.2}$ mm s⁻¹, 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).

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_{L\mathrm{F}}^s$ in panel (b), estimated random uncertainty profiles $u_{L\mathrm{F}}^r$ 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).

3 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 COSMIC 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 $X_{\mathrm{r}}\in \mathit{\left\{}{L}_{\mathrm{F},k}\right(t)$, D_r,k(t), α_G,k(t), α_G,k(z_a), α_M,k(z_a), α_F,k(z_a), and α_r(z_a)}), the estimated systematic uncertainty profiles $u_{X\mathrm{r}}^s$, the estimated random uncertainty profiles $u_{X\mathrm{r}}^r$, representative correlation functions R_Xr,i (with i such that $z_{\mathrm{a},i}\in \mathit{\{}\mathrm{10},\mathrm{30},\mathrm{50},\mathrm{70}\mathrm{km}\mathit{\}}$), 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 $X_{\mathrm{m}}\in \mathit{\{}{L}_{\mathrm{m}},D_{\mathrm{m}},{\mathit{\alpha }}_{\mathrm{m}}\mathit{\}}$ for comparison.

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.

3.1 Doppler shift retrieval

3.1.1 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 high-frequency 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_{L\mathrm{F}}^r$ 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_{L\mathrm{F}}^r$ and R_LF, we use Eqs. (A7) and (A8).

To propagate the estimated systematic uncertainty $u_{L\mathrm{r}}^s$, 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_{L\mathrm{r}}^s$ is chosen to be constant down to z_aGrads, the filter has little effect, and $u_{L\mathrm{F}}^s$, shown in Fig. 4b, is essentially equal to $u_{L\mathrm{r}}^s$, shown in Fig. 3b.

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_{D\mathrm{r}}^s$ in panel (b), estimated random uncertainty profiles $u_{D\mathrm{r}}^r$ 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).

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).

3.1.2 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_{D\mathrm{r}}^r$ (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 and 5e. The correlation 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_{L\mathrm{F}}^s$ and get $u_{D\mathrm{r}}^s$ (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⁻¹ occurs.

Figure 6 Results for geometric optics bending angle profiles α_G (with model profile α_m for comparison) in panel (a), estimated systematic uncertainty profiles $u_{\mathit{\alpha }\mathrm{G}}^s$ in panel (b), estimated random uncertainty profiles $u_{\mathit{\alpha }\mathrm{G}}^r$ 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).

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.

3.2 Bending angle retrieval

3.2.1 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_{\mathit{\alpha }\mathrm{G}}^r$, 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_{\mathit{\alpha }\mathrm{G}}^s$ 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.

3.2.2 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_{\mathit{\alpha }\mathrm{W}}^s$, the random uncertainty profile $u_{\mathit{\alpha }\mathrm{W}}^r$, 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 and Kirchengast (2015, 2018). The WO retrieval and uncertainty propagation results are supplied up to 20 km impact altitude by the WO algorithms.

3.2.3 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_{\mathrm{a}}^{\mathrm{GW}}$ in a transition range $z_{\mathrm{a}}^{\mathrm{GW}}\pm \mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}$, 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).

Figure 7 Results for filtered bending angle profiles α_F (with model profile α_m for comparison) in panel (a), systematic uncertainty profiles $u_{\mathit{\alpha }\mathrm{F}}^s$ in panel (b), random uncertainty profiles $u_{\mathit{\alpha }\mathrm{F}}^r$ 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).

3.3 Atmospheric bending angle derivation

3.3.1 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\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$.

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_{\mathrm{c}\mathrm{2}}=\mathrm{10}/\mathrm{7}$ Hz in this example case.

The relevant covariance-propagated random uncertainties $u_{\mathit{\alpha }F,k}^r$ 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).

Figure 8 Results for atmospheric bending angle profile α_r (with model profile α_m for comparison) in panel (a), systematic uncertainty profile $u_{\mathit{\alpha }\mathrm{r}}^s$ in panel (b), random uncertainty profile $u_{\mathit{\alpha }\mathrm{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).

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_{\mathrm{c}\mathrm{2}}=\mathrm{10}/\mathrm{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⁻¹, is near 1.1 km (Fig. 7f).

3.3.2 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_{\mathit{\alpha }\mathrm{r}}^r$ is therefore considerably larger than $u_{\mathit{\alpha }\mathrm{F}\mathrm{1}}^r$, and $u_{\mathit{\alpha }\mathrm{F}\mathrm{2}}^r$ (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_{\mathit{\alpha }\mathrm{r}}^s$, 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.

4 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

$$\begin{array}{}\text{(9)}& C_{L\mathrm{r},ij}^{\mathrm{MC}}=u_{L\mathrm{r},i}^{r,\mathrm{Std}}\cdot u_{L\mathrm{r},j}^{r,\mathrm{Std}}\cdot {\mathit{\delta }}_{ij},\end{array}$$

by a large number M of draws $L_{\mathrm{r}}^{\mathrm{T}}+{\mathit{ϵ}}_{L\mathrm{r},j}^r$ (with $j\in \mathit{\{}\mathrm{1},\mathrm{\dots },M\mathit{\}}$ 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\in \mathit{\left\{}{L}_{\mathrm{F},kj}\right(t),D_{\mathrm{r},kj}(t),{\mathit{\alpha }}_{\mathrm{G},kj}(z_{\mathrm{a}}),{\mathit{\alpha }}_{\mathrm{F},kj}(z_{\mathrm{a}}),{\mathit{\alpha }}_{\mathrm{r},j}(z_{\mathrm{a}}\left)\mathit{\right\}}$ and $k\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$). From these individual realizations the mean profiles X^MC and the covariance matrices ${\mathbf{C}}_X^{\mathrm{MC}}$,

$$\begin{array}{ll}{\displaystyle }{\mathbf{C}}_X^{\mathrm{MC}}=& {\displaystyle \frac{\mathrm{1}}{M-\mathrm{1}}}\left[\left(X_{\mathrm{1}}-X^{\mathrm{MC}}\right){\left(X_{\mathrm{1}}-X^{\mathrm{MC}}\right)}^T+\mathrm{\dots }\right.\\ \text{(10)}& {\displaystyle }& {\displaystyle }\left.+\left(X_M-X^{\mathrm{MC}}\right){\left(X_M-X^{\mathrm{MC}}\right)}^T\right],\end{array}$$

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 ${\mathbf{C}}_X^{\mathrm{CP}}$ are compared to the MC-derived matrices ${\mathbf{C}}_X^{\mathrm{MC}}$. In order to be able to attribute potential changes between CP and MC covariance matrices better, we decompose C_X into $u_X^r$ and R_X (Eqs. A7 and A8), and compare them separately.

Figure 9 Results from the validation of CP covariance matrices ${\mathbf{C}}_{L\mathrm{r}}^{\mathrm{CP}}$ (“CP”) by MC covariance matrices ${\mathbf{C}}_{L\mathrm{r}}^{\mathrm{MC}}$ (“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_{L\mathrm{r}\mathrm{1}}^r$ in light blue); the middle column for f_T2 (CP in red, MC in black, in panel b $u_{L\mathrm{r}\mathrm{2}}^r$ 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 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.

Figure 9a (light blue) and b (orange) show the random uncertainties $u_{L\mathrm{r},\mathrm{1}}^r$ and $u_{L\mathrm{r},\mathrm{2}}^r$, respectively, which characterize the input distribution and from which the random error profiles ${\mathit{ϵ}}_{L\mathrm{r},j}^r$ are drawn. They also show the CP results for the random uncertainty $u_{L\mathrm{F}\mathrm{1}}^r$ (dark blue in Fig. 9a) and $u_{L\mathrm{F}\mathrm{2}}^r$ (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_{D\mathrm{r}\mathrm{1}}^r$ (Fig. 9d), $u_{D\mathrm{r}\mathrm{2}}^r$ (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_{\mathrm{r}}^{\mathrm{T}}$ 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_{\mathit{\alpha }\mathrm{F}\mathrm{2}}^r$ is smaller than $u_{\mathit{\alpha }\mathrm{F}\mathrm{1}}^r$, even though $u_{\mathit{\alpha }\mathrm{G}\mathrm{2}}^r$ was larger than $u_{\mathit{\alpha }\mathrm{G}\mathrm{1}}^r$. 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.

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_{L\mathrm{F}\mathrm{1}}^r$ (heavy) and systematic $u_{L\mathrm{F}\mathrm{1}}^s$ (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).

In order to demonstrate that a full CP is necessary to propagate random uncertainties correctly, we also calculated random uncertainties $u_{\mathit{\alpha }\mathrm{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.

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_{X\mathrm{r}}^r$ (heavy) and mean systematic uncertainty $u_{X\mathrm{r}}^s$ (light) profiles (a, d, g, j); the latter is shown as $\mathrm{10}\times u_{L\mathrm{F}\mathrm{1}}^s$ (in a) and $\mathrm{100}\times u_{X\mathrm{r}}^s$ (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).

5 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 random 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_{L\mathrm{F}\mathrm{1}}^r$ and $u_{L\mathrm{F}\mathrm{1}}^s$ 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_{L\mathrm{F}\mathrm{1}}^r$ 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_{L\mathrm{r}\mathrm{1}}^r$ 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_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}$ (cf. Table 2).

Compared to CHAMP, the mean random uncertainty $u_{L\mathrm{F}\mathrm{1}}^r$ (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 COSMIC 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_{L\mathrm{F}\mathrm{1}}^r$ 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_{L\mathrm{F}\mathrm{1}}^s$ – 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_{L\mathrm{F}\mathrm{1}}^s$, $u_{L\mathrm{F}\mathrm{1}}^r$, 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. 11g–i).

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_{\mathit{\alpha }\mathrm{r}}^s$, 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.

6 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.

Appendix A: Algorithm description

In this appendix the rOPS L1b uncertainty propagation algorithm is introduced, following the L1b retrieval chain (Fig. 2; Sect. 3) step by step, starting with excess phase profile L_r as input and proceeding to L_F, D_r, α_G, α_M, α_F and finally α_r. The relevant variable definitions and symbol explanations are summarized in Tables 1 and 2. A fully detailed algorithmic description is provided by Kirchengast et al. (2017b).

If not stated otherwise, elements of the vector-type vertical profiles are addressed using subscript i (with $i\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N\mathit{\}}$), and optionally j (with $j\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N\mathit{\}}$), running from top downward towards the bottom of the profile, where N is the number of vertical grid levels. Until the interpolation of all quantities to the common monotonic impact altitude grid z_a, all quantities are provided on an equidistant 50 Hz time grid t with grid points t_i.

All steps in Sects. A1 and A2 are applied to each of the GNSS transmitter channels' carrier frequencies f_Tk, as also indicated by the index k in Fig. 2. In the notation of these sections we therefore suppress the index k for the convenience of simplified readability. Also for conciseness we write the estimated systematic uncertainty equations only for the total systematic uncertainties u^s and briefly address the type of the relevant components (whether basic systematic uncertainty u^b or apparent systematic uncertainty u^a) in the surrounding text.

A1 Doppler shift retrieval

A1.1 Basic low-pass filtering

The Doppler differentiation (item 1.4 in Fig. 2) would potentially amplify high-frequency noise in the excess phase profiles. To avoid this amplification, a BWS low-pass filter (e.g., Smith, 1999) is applied onto the excess phase profiles first (item 1.2 in Fig. 2).

For this basic filtering the relative cutoff frequency f_c∕f_s is set to 0.05, equivalent to f_c=2.5 Hz, 21 grid points, or a cutoff period ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{1}/f_{\mathrm{c}}=\mathrm{0.4}$ s, for the standard sampling rate f_s of 50 Hz used for all RO profiles in the L1b processor of the rOPS. The corresponding sample width of the Blackman window $\stackrel{\mathrm{̃}}{M}$ (with samples $m\in \mathit{\{}\mathrm{0},\mathrm{\dots },M\mathit{\}}$) is set to $\stackrel{\mathrm{̃}}{M}=\mathrm{2}\cdot f_{\mathrm{s}}/f_{\mathrm{c}}$, yielding 41 grid points. This ensures a reliable filter performance, also allowing the vertical resolution of the filtered data to be robustly quantified.

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 (“$\mathit{\delta }{L}_{\mathrm{Fm}}^{\mathrm{BWS}}$”, red), and the boxcar filtered profiles with M=21 points (“$\mathit{\delta }{L}_{\mathrm{Fm}}^{\mathrm{BC}\mathrm{21}}$”, green) and M=11 points (“$\mathit{\delta }{L}_{\mathrm{Fm}}^{\mathrm{BC}\mathrm{11}}$”, blue).

With such a design, the BWS low-pass filter combines efficient removal of high-frequency noise with a narrow smoothing window. The BWS filter thus achieves a better smoothing effect, while keeping a w_LF of higher resolution than a simple moving-average BC filter. Based on a time segment of a few seconds of the excess phase delta profile of the COSMIC example event (also used for Figs. 3 to 8), Fig. A1 illustrates how the BWS filter compares to boxcar filters of 11 and 21 grid points. The corresponding filter functions are 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 “boxcar-equivalent width”, is therefore its full width at half maximum (see Fig. A1a), corresponding to $\stackrel{\mathrm{̃}}{M}/\mathrm{4}+\mathrm{1}$ samples with our design.

The actually used sample width M of the BWS filter is equal to $\stackrel{\mathrm{̃}}{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\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N\mathit{\}}$, and N being the profile length in grid points), the filter width M is thus

$$\begin{array}{}\text{(A1)}& {\displaystyle }M=\left\{\begin{array}{ll}\stackrel{\mathrm{̃}}{M}& \text{ for }\stackrel{\mathrm{̃}}{M}/\mathrm{2}<i<N-\stackrel{\mathrm{̃}}{M}/\mathrm{2}\\ \mathrm{2}i-\mathrm{1}& \text{ for }\mathrm{1}\le i\le \stackrel{\mathrm{̃}}{M}/\mathrm{2}\\ \mathrm{2}(N-i)+\mathrm{1}& \text{ for }N-\stackrel{\mathrm{̃}}{M}/\mathrm{2}<i<N\end{array}\right..\end{array}$$

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 $\mathit{\delta }{L}_{\mathrm{rm}}=L_{\mathrm{r}}-L_{\mathrm{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)

$$\begin{array}{}\text{(A2)}& L_{\mathrm{F}i}=L_{\mathrm{m}i}+\sum _{j=\mathrm{0}}^N{A}_{ij}^{\mathrm{BWS}}\cdot \mathit{\delta }{L}_{\mathrm{rm}j}\end{array}$$

for the filtered excess phase, where $j\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N\mathit{\}}$. The band matrix operator A^BWS has elements

$$\begin{array}{}\text{(A3)}& {\displaystyle }{A}_{ij}^{\mathrm{BWS}}=\left\{\begin{array}{ll}\mathrm{0}& \text{ for }j<i-M/\mathrm{2}\text{ and}\\ & \text{ for }j>i+M/\mathrm{2}\\ w_{j-i+M/\mathrm{2}}& \text{ for }i-M/\mathrm{2}<j<i+M/\mathrm{2}\end{array}\right..\end{array}$$

The central filter weight $w_{\mathrm{0}+M/\mathrm{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/\mathrm{2}$ (and therefore $\mathrm{0}\le m\le M$), each single BWS weight is calculated using

$$\begin{array}{}\text{(A4)}& w_m={\displaystyle \frac{w_{\mathrm{raw},m}}{{\sum }_{m=\mathrm{0}}^M{w}_{\mathrm{raw},m}}}\end{array}$$

and

$$\begin{array}{}\text{(A5)}& {\displaystyle }{w}_{\mathrm{raw},m}=\left\{\begin{array}{l}\frac{\sin \left(\mathrm{2}\mathit{\pi }{f}_{\mathrm{c}}/f_{\mathrm{s}}\left(m-M/\mathrm{2}\right)\right)}{m-M/\mathrm{2}}[\mathrm{0.42}-\mathrm{0.5}\cos \\ \left(\mathrm{2}\mathit{\pi }\frac{m}{M}\right)+\mathrm{0.08}\cos \left(\mathrm{4}\mathit{\pi }\frac{m}{M}\right)]\\ \text{ for }m\ne M/\mathrm{2}\\ \mathrm{2}\mathit{\pi }{f}_{\mathrm{c}}/f_{\mathrm{s}}\\ \text{ for }m=M/\mathrm{2}\end{array}\right..\end{array}$$

The estimated random uncertainty is then propagated by covariance propagation (item 1.3 in Fig. 2),

$$\begin{array}{}\text{(A6)}& {\mathbf{C}}_{L\mathrm{F}}={\mathbf{A}}^{\mathrm{BWS}}\cdot {\mathbf{C}}_{L\mathrm{r}}\cdot {\left({\mathbf{A}}^{\mathrm{BWS}}\right)}^{\mathrm{T}}.\end{array}$$

The random uncertainty profile $u_{L\mathrm{F}}^r$ 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

$$\begin{array}{}\text{(A7)}& u_{L\mathrm{F},i}^r=\sqrt{C_{L\mathrm{F},ii}}\end{array}$$

and

$$\begin{array}{}\text{(A8)}& R_{L\mathrm{F},ij}={\displaystyle \frac{C_{L\mathrm{F},ij}}{u_{L\mathrm{F},i}^r{u}_{L\mathrm{F},j}^r}}.\end{array}$$

The correlation length profile l_Lr has elements

$$\begin{array}{}\text{(A9)}& l_{L\mathrm{r},i}={\left.{\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}\right|}_i\cdot |t_i-t(R_{L\mathrm{F},ij}=\mathrm{1}/e\left)\right|\end{array}$$

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.,

$$\begin{array}{}\text{(A10)}& u_{L\mathrm{F}i}^s=\sum _{j=\mathrm{0}}^N{A}_{ij}^{\mathrm{BWS}}\cdot u_{L\mathrm{r}j}^s.\end{array}$$

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,

$$\begin{array}{}\text{(A11)}& {\mathit{\tau }}_{L\mathrm{F}}\approx {\displaystyle \frac{\mathrm{1}}{f_{\mathrm{c}}+\mathrm{\Delta }{f}_{\mathrm{c}}/\mathrm{2}}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{f}_{\mathrm{c}}}},\end{array}$$

with our design choice $\stackrel{\mathrm{̃}}{M}=\mathrm{2}(f_{\mathrm{s}}/f_{\mathrm{c}})$ and with the BWS filter stopband-to-passband transition width being (Smith, 1999)

$$\begin{array}{}\text{(A12)}& \mathrm{\Delta }{f}_{\mathrm{c}}\approx {\displaystyle \frac{\mathrm{4}{f}_{\mathrm{s}}}{\stackrel{\mathrm{̃}}{M}}}.\end{array}$$

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:

$$\begin{array}{}\text{(A13)}& w_{L\mathrm{F},i}={\left.{\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}\right|}_i\cdot {\mathit{\tau }}_{L\mathrm{F}},\end{array}$$

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,

$$\begin{array}{}\text{(A14)}& \mathit{\delta }{L}_{\mathrm{Fm}}=L_{\mathrm{F}}-L_{\mathrm{m}},\end{array}$$

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 five-point derivative δD_rm,i is given by

$$\begin{array}{ll}\text{(A15)}& {\displaystyle }& {\displaystyle }\mathit{\delta }{D}_{\mathrm{rm},i}={\left.{\displaystyle \frac{\mathrm{d}\mathit{\delta }{L}_{\mathrm{Fm}}\left(t\right)}{\mathrm{d}t}}\right|}_i{\displaystyle }& {\displaystyle }={\displaystyle \frac{-\mathit{\delta }{L}_{\mathrm{Fm},i-\mathrm{2}}+\mathrm{8}\mathit{\delta }{L}_{\mathrm{Fm},i-\mathrm{1}}-\mathrm{8}\mathit{\delta }{L}_{\mathrm{Fm},i+\mathrm{1}}+\mathit{\delta }{L}_{\mathrm{Fm},i+\mathrm{2}}}{-t_{i-\mathrm{2}}+\mathrm{8}{t}_{i-\mathrm{1}}-\mathrm{8}{t}_{i+\mathrm{1}}+t_{i+\mathrm{2}}}}\end{array}$$

for each of the frequencies (e.g., Syndergaard, 1999). This can be expressed in matrix form as

$$\begin{array}{}\text{(A16)}& D_{\mathrm{r},i}=D_{\mathrm{m},i}+\mathit{\delta }{D}_{\mathrm{rm},i}=D_{\mathrm{m},i}+\sum _{j=\mathrm{1}}^N{A}_{ij}^{\mathrm{L}\mathrm{2}\mathrm{D}}\cdot \mathit{\delta }{L}_{\mathrm{Fm},j},\end{array}$$

using matrix operator A^L2D with

$$\begin{array}{ll}\text{(A17)}& {\displaystyle }& {\displaystyle }{\mathbf{A}}^{\mathrm{L}\mathrm{2}\mathrm{D}}={\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{12}\mathrm{\Delta }t}}\left[\begin{array}{rrrrrrrrr}-\mathrm{18}& \mathrm{24}& -\mathrm{6}& \mathrm{0}& \mathrm{0}& \mathrm{0}& & \mathrm{0}& \mathrm{0}\\ -\mathrm{6}& \mathrm{0}& -\mathrm{6}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\ddots }& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{8}& \mathrm{0}& -\mathrm{8}& \mathrm{1}& \mathrm{0}& & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{8}& \mathrm{0}& -\mathrm{8}& \mathrm{1}& \mathrm{\ddots }& \mathrm{0}& \mathrm{0}\\ & \mathrm{\ddots }& & \mathrm{\ddots }& & \mathrm{\ddots }& & \mathrm{\ddots }& \\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\ddots }& -\mathrm{8}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& & \mathrm{0}& -\mathrm{6}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& & \mathrm{24}& -\mathrm{18}\end{array}\right],\end{array}$$

where $\mathrm{\Delta }t=t_{i+\mathrm{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

$$\begin{array}{}\text{(A18)}& {\mathbf{C}}_D={\mathbf{A}}^{\mathrm{L}\mathrm{2}\mathrm{D}}\cdot {\mathbf{C}}_{L\mathrm{F}}\cdot {\left({\mathbf{A}}^{\mathrm{L}\mathrm{2}\mathrm{D}}\right)}^{\mathrm{T}}.\end{array}$$

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) using

$$\begin{array}{}\text{(A19)}& u_{D\mathrm{r},i}^r=\sqrt{C_{D\mathrm{r},ii}}\end{array}$$

and

$$\begin{array}{}\text{(A20)}& R_{D\mathrm{r},ij}={\displaystyle \frac{C_{D\mathrm{r},ij}}{u_{D\mathrm{r},i}^r{u}_{D\mathrm{r},j}^r}}.\end{array}$$

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.,

$$\begin{array}{}\text{(A21)}& u_{D\mathrm{r},i}^s=\sum _{j=\mathrm{1}}^N{A}_{ij}^{\mathrm{L}\mathrm{2}\mathrm{D}}\cdot u_{L\mathrm{F},j}^s.\end{array}$$

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

$$\begin{array}{}\text{(A22)}& D_{\mathrm{r},i}=\left[v_{\mathrm{R},i}\cos \left({\mathit{\varphi }}_{\mathrm{R},i}\right)-v_{\mathrm{T},i}\cos \left({\mathit{\varphi }}_{\mathrm{T},i}\right)\right]-{\dot{r}}_{\mathrm{RT},i},\end{array}$$

where

$$\begin{array}{}\text{(A23)}& {\mathit{\varphi }}_{\mathrm{R},i}={\mathit{\eta }}_{\mathrm{R},i}-\arcsin \left({\displaystyle \frac{a_{t,i}}{r_{\mathrm{R},i}}}\right)\end{array}$$

and

$$\begin{array}{}\text{(A24)}& {\mathit{\varphi }}_{\mathrm{T},i}=\left(\mathit{\pi }-{\mathit{\eta }}_{\mathrm{T},i}\right)-\arcsin \left({\displaystyle \frac{a_{t,i}}{r_{\mathrm{T},i}}}\right)\end{array}$$

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_{\mathrm{R},i}=\left|{\mathit{v}}_{\mathrm{R},i}\right|$ is the receiver velocity; $r_{\mathrm{R},i}=\left|{\mathit{r}}_{\mathrm{R},i}\right|$ 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); and ${\dot{r}}_{\mathrm{RT},i}={\left|\frac{d({\mathit{r}}_{\mathrm{T}}-{\mathit{r}}_{\mathrm{R}})}{dt}\right|}_i$ 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:

$$\begin{array}{}\text{(A25)}& {\mathit{\alpha }}_{\mathrm{G},i}={\mathit{\theta }}_{\mathrm{RT},i}-\arccos \left({\displaystyle \frac{a_{t,i}}{r_{\mathrm{R},i}}}\right)-\arccos \left({\displaystyle \frac{a_{t,i}}{r_{\mathrm{T},i}}}\right),\end{array}$$

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)

$$\begin{array}{}\text{(A26)}& z_{\mathrm{a},i}=a_{t,j\mathrm{1}}-h_{\mathrm{G}}-R_{\mathrm{C}},\end{array}$$

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_{D\mathrm{r}}^r$ and the uncertainties of the bending angle $u_{\mathit{\alpha }\mathrm{G}}^r$, using

$$\begin{array}{}\text{(A27)}& u_{\mathit{\alpha }\mathrm{G}\left(t\right),i}^r\approx {\left.-{\left({\displaystyle \frac{\mathrm{d}{a}_{\mathrm{SL}}}{\mathrm{d}t}}\right)}^{-\mathrm{1}}\right|}_i{u}_{D\mathrm{r}\left(t\right),i}^r,\end{array}$$

where a_SL is the straight-line impact parameter. These bending angle uncertainties $u_{\mathit{\alpha }\mathrm{G}}^r$ 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

$$\begin{array}{}\text{(A28)}& u_{\mathit{\alpha }\mathrm{G}\left(z_{\mathrm{a}}\right),i}^r\approx {\left.-{\left({\displaystyle \frac{\mathrm{d}{a}_t^T}{\mathrm{d}t}}\right)}^{-\mathrm{1}}\right|}_i{u}_{D\mathrm{r}\left(t\right),i}^r,\end{array}$$

where $a_t^T$ is the “true” impact parameter. We use the forward-modeled impact parameter a_tm instead (i.e., adopt $a_{t\mathrm{m}}=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

$$\begin{array}{}\text{(A29)}& u_{\mathit{\alpha }\mathrm{G},i}^r=f_{u\mathit{\alpha }\mathrm{lin}}\cdot u_{\mathit{\alpha }\mathrm{G}\left(z_{\mathrm{a}}\right),i}^r.\end{array}$$

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_{\mathit{\alpha }\mathrm{G}}^r$ 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 uncertainties (item 2.8 in Fig. 2),

$$\begin{array}{}\text{(A30)}& C_{\mathit{\alpha }\mathrm{G},ij}=u_{\mathit{\alpha }\mathrm{G},i}^r\cdot u_{\mathit{\alpha }\mathrm{G},j}^r\cdot R_{D\mathrm{r},ij}.\end{array}$$

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_{at}^s$, and consequentially the bending angle uncertainty $u_{\mathit{\alpha }\mathrm{G}}^s$, are taken into account: systematic errors in the Doppler shift, i.e., $u_{D\mathrm{r}}^s$; systematic errors in the velocities of the satellites, i.e., $u_{\mathit{v}R}^s$ and $u_{\mathit{v}T}^s$; and systematic errors in the positions of the satellites, i.e., $u_{\mathit{r}R}^s$ and $u_{\mathit{r}T}^s$. The latter two orbit-borne types are interpreted as apparent systematic uncertainties (Sect. 2.2), while the excess phase-borne uncertainty $u_{D\mathrm{r}}^s$ is a basic systematic uncertainty.

For the propagation of these estimated systematic uncertainties to $u_{\mathit{\alpha }\mathrm{G}}^s$, 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_{\mathrm{r},i},r_{\mathrm{R},i},r_{\mathrm{T},i},v_{\mathrm{R},i},v_{\mathrm{T},i})$, with $i=\mathrm{1},\mathrm{2},\mathrm{\dots },N$. Adopting the first-order deviations to represent the estimated systematic uncertainties, we obtain

$$\begin{array}{ll}\text{(A31)}& {\displaystyle }& {\displaystyle }{u}_{at,i}^s={\displaystyle \frac{\mathrm{1}}{k_{at,i}}}{\displaystyle }& {\displaystyle }\sqrt{{\left(u_{D\mathrm{r},i}^s\right)}^{\mathrm{2}}+{\left(k_{v\mathrm{R},i}\cdot u_{v\mathrm{R},i}^s\right)}^{\mathrm{2}}+{\left(k_{r\mathrm{R},i}\cdot u_{r\mathrm{R},i}^s\right)}^{\mathrm{2}}+{\left(k_{v\mathrm{T},i}\cdot u_{v\mathrm{T},i}^s\right)}^{\mathrm{2}}+{\left(k_{r\mathrm{T},i}\cdot u_{r\mathrm{T},i}^s\right)}^{\mathrm{2}}},\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{k}_{at,i}={\left.{\displaystyle \frac{\partial D_{\mathrm{r}}}{\partial {\mathit{\varphi }}_{\mathrm{R}}}}\right|}_i\cdot {\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{R}}}{\partial a_t}}\right|}_i+{\left.{\displaystyle \frac{\partial D_{\mathrm{r}}}{\partial {\mathit{\varphi }}_{\mathrm{T}}}}\right|}_i\cdot {\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{T}}}{\partial a_t}}\right|}_i=-v_{\mathrm{R},i}\\ {\displaystyle }& {\displaystyle }\cdot \sin {\mathit{\varphi }}_{\mathrm{R}i}\cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{r_{\mathrm{R},i}^{\mathrm{2}}-a_{t,i}^{\mathrm{2}}}}}-v_{\mathrm{T},i}\cdot \sin {\mathit{\varphi }}_{\mathrm{T}i}{\displaystyle \frac{\mathrm{1}}{\sqrt{r_{\mathrm{T},i}^{\mathrm{2}}-a_{t,i}^{\mathrm{2}}}}},\\ {\displaystyle }& {\displaystyle }{k}_{v\mathrm{R},i}={\left.-{\displaystyle \frac{\partial D_{\mathrm{r}}}{\partial v_{\mathrm{R}}}}\right|}_i=-\cos {\mathit{\varphi }}_{\mathrm{R}i},\\ {\displaystyle }& {\displaystyle }{k}_{v\mathrm{T},i}={\left.-{\displaystyle \frac{\partial D_{\mathrm{r}}}{\partial v_{\mathrm{T}}}}\right|}_i=-\cos {\mathit{\varphi }}_{\mathrm{T}i},\\ {\displaystyle }& {\displaystyle }{k}_{r\mathrm{R},i}={\left.{\displaystyle \frac{\partial D_{\mathrm{r}}}{\partial {\mathit{\varphi }}_{\mathrm{R}}}}\right|}_i\cdot {\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{R}}}{\partial r_{\mathrm{R}}}}\right|}_i={\displaystyle \frac{v_{\mathrm{R},i}\cdot \sin {\mathit{\varphi }}_{\mathrm{R}i}\cdot a_{t,i}}{r_{\mathrm{R},i}\sqrt{r_{\mathrm{R},i}^{\mathrm{2}}-a_{t,i}^{\mathrm{2}}}}},\\ \text{(A32)}& {\displaystyle }& {\displaystyle }{k}_{r\mathrm{T},i}={\left.{\displaystyle \frac{\partial D_{\mathrm{r}}}{\partial {\mathit{\varphi }}_{\mathrm{T}}}}\right|}_i\cdot {\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{T}}}{\partial r_{\mathrm{T}}}}\right|}_i={\displaystyle \frac{v_{\mathrm{T},i}\cdot \sin {\mathit{\varphi }}_{\mathrm{T}i}\cdot a_{t,i}}{r_{\mathrm{T},i}\sqrt{r_{\mathrm{T},i}^{\mathrm{2}}-a_{t,i}^{\mathrm{2}}}}}.\end{array}$$

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 $\mathrm{\Delta }{\dot{r}}_{\mathrm{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_{\mathit{\eta }}^s\approx \mathrm{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

$$\begin{array}{ll}\text{(A33)}& {\displaystyle }& {\displaystyle }{u}_{\mathit{\alpha }\mathrm{G},i}^s={\displaystyle }& {\displaystyle }\sqrt{{\left(u_{\mathit{\theta }\mathrm{RT},i}^s\right)}^{\mathrm{2}}+{\left(k_{at,i}\cdot u_{at,i}^s\right)}^{\mathrm{2}}+{\left(k_{r\mathrm{R},i}\cdot u_{r\mathrm{R},i}^s\right)}^{\mathrm{2}}+{\left(k_{r\mathrm{T},i}\cdot u_{r\mathrm{T},i}^s\right)}^{\mathrm{2}}},\end{array}$$

where

$$\begin{array}{ll}\text{(A34)}& {\displaystyle }& {\displaystyle }{k}_{at,i}={\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial a_t}}\right|}_i={\displaystyle \frac{\mathrm{1}}{\sqrt{r_{\mathrm{R},i}^{\mathrm{2}}+a_{t,i}^{\mathrm{2}}}}}+{\displaystyle \frac{\mathrm{1}}{\sqrt{r_{\mathrm{T},i}^{\mathrm{2}}+a_{t,i}^{\mathrm{2}}}}},{\displaystyle }& {\displaystyle }{k}_{r\mathrm{R},i}={\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial r_{\mathrm{R}}}}\right|}_i=-{\displaystyle \frac{a_{t,i}}{r_{\mathrm{R},i}\sqrt{r_{\mathrm{R},i}^{\mathrm{2}}-a_{t,i}^{\mathrm{2}}}}},\\ {\displaystyle }& {\displaystyle }{k}_{r\mathrm{T},i}={\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial r_{\mathrm{T}}}}\right|}_i=-{\displaystyle \frac{a_{t,i}}{r_{\mathrm{T},i}\sqrt{r_{\mathrm{T},i}^{\mathrm{2}}-a_{t,i}^{\mathrm{2}}}}}.\end{array}$$

In practice we separately calculate the basic and apparent systematic uncertainty estimates ($u_{\mathit{\alpha }\mathrm{G}}^{\mathrm{b}}$ from the first RHS terms in Eqs. A31 and A33, $u_{\mathit{\alpha }\mathrm{G}}^{\mathrm{a}}$ from the orbit-borne terms) and afterwards obtain $u_{\mathit{\alpha }\mathrm{G}}^s$ 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 (2015, 2018). Along with the state profile, the systematic uncertainty profile $u_{\mathit{\alpha }\mathrm{W}}^s$, 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_{\mathit{\alpha }\mathrm{W}}^r$ 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_{\mathit{\alpha }\mathrm{W}}^s$ is composed of a basic systematic uncertainty $u_{\mathit{\alpha }\mathrm{W}}^{\mathrm{b}}$, propagated through the wave-optical retrieval from the excess phase uncertainty $u_{L\mathrm{r}}^s$, and an apparent systematic uncertainty $u_{\mathit{\alpha }\mathrm{W}}^{\mathrm{a}}$, 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 $\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}=\mathrm{2}$ km around transition altitude $z_{\mathrm{a}}^{\mathrm{GW}}$, allowed within 9 km to 14 km, estimated from α_G data quality. The resulting merged bending angle profile α_M is

$$\begin{array}{}\text{(A35)}& {\mathit{\alpha }}_{\mathrm{M},i}={\mathit{\gamma }}_i\cdot {\mathit{\alpha }}_{\mathrm{G},i}+(\mathrm{1}-{\mathit{\gamma }}_i)\cdot {\mathit{\alpha }}_{\mathrm{W},i},\end{array}$$

where the weighting profile γ is formulated as

$$\begin{array}{}\text{(A36)}& {\mathit{\gamma }}_i=\left\{\begin{array}{l}\mathrm{1}\text{for }{z}_{\mathrm{a},i}\ge z_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}\\ \mathrm{0.5}\cdot \left[\sin \left(\frac{\mathit{\pi }}{\mathrm{2}}\cdot \frac{z_{\mathrm{a},i}-z_{\mathrm{a}}^{\mathrm{GW}}}{\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}}\right)+\mathrm{1}\right]\\ \text{ for }|z_{\mathrm{a},i}-z_{\mathrm{a}}^{\mathrm{GW}}|<\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}\\ \mathrm{0}\text{for }{z}_{\mathrm{a},i}\le z_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}.\end{array}\right.\end{array}$$

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_{\mathrm{a}}^{\mathrm{GW}}\pm \mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}$) we can assign

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{C}_{\mathit{\alpha }\mathrm{M},ij}=\\ \text{(A37)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{C}_{\mathit{\alpha }\mathrm{G},ij}& \text{ for }{z}_{\mathrm{aTop}}>z_{\mathrm{a},i}>z_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}\\ & \text{ and }{z}_{\mathrm{aTop}}>z_{\mathrm{a},j}>z_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}\\ C_{\mathit{\alpha }\mathrm{W},ij}& \text{ for }{z}_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}>z_{\mathrm{a},i}>z_{\mathrm{aBot}}\\ & \text{ and }{z}_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}>z_{\mathrm{a},j}>z_{\mathrm{aBot}}\\ \mathrm{0}& \text{ for }{z}_{\mathrm{aTop}}>z_{\mathrm{a},i}>z_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}\\ & \text{ and }{z}_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}>z_{\mathrm{a},j}>z_{\mathrm{aBot}}\\ \mathrm{0}& \text{ for }{z}_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}>z_{\mathrm{a},i}>z_{\mathrm{aBot}}\\ & \text{ and }{z}_{\mathrm{aTop}}>z_{\mathrm{a},j}>z_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}\end{array}\right.,\end{array}$$

while within the merging zone we can assign

$$\begin{array}{}\text{(A38)}& C_{\mathit{\alpha }\mathrm{M},ij}={\mathit{\gamma }}_i{\mathit{\gamma }}_j{C}_{\mathit{\alpha }\mathrm{G},ij}+(\mathrm{1}-{\mathit{\gamma }}_i)(\mathrm{1}-{\mathit{\gamma }}_j)C_{\mathit{\alpha }\mathrm{W},ij},\end{array}$$

wherein i is understood such that $z_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}>z_{\mathrm{a},i}>z_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}$ and j such that $z_{\mathrm{aTop}}>z_{\mathrm{a},j}>z_{\mathrm{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_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}>z_{\mathrm{a},j}>z_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}$ and $z_{\mathrm{aTop}}>z_{\mathrm{a},i}>z_{\mathrm{aBot}}$, are calculated according to the same formula.

Due to the linear relation between α_M, α_G, and α_W, expressed by Eq. (A35), a bias $u_{\mathit{\alpha }\mathrm{G}}^s$ in the GO bending angle and a bias $u_{\mathit{\alpha }\mathrm{W}}^s$ 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_{\mathit{\alpha }\mathrm{M}}^s$ according to (item 2.9 in Fig. 2)

$$\begin{array}{}\text{(A39)}& u_{\mathit{\alpha }\mathrm{M},i}^s={\mathit{\gamma }}_i\cdot u_{\mathit{\alpha }\mathrm{G},i}^s+(\mathrm{1}-{\mathit{\gamma }}_i)\cdot u_{\mathit{\alpha }\mathrm{W},i}^s.\end{array}$$

In practice this formulation is again applied separately for the basic and apparent systematic uncertainty estimates, afterwards obtaining the $u_{\mathit{\alpha }\mathrm{M}}^s$ 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_{\mathrm{a}}^{\mathrm{GW}}+\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}$, equal to the WO resolution w_αW below $z_{\mathrm{a}}^{\mathrm{GW}}-\mathrm{\Delta }{z}_{\mathrm{a}}^{\mathrm{GW}}$ 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 GPS case (with $k\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$, f_T1=1.57542  GHz, and f_T2=1.22760  GHz) are in the prime focus of this section.

This concerns in particular special provisions for the minor (L2) channel noise filtering and its tropospheric extrapolation. In general the algorithms are applicable for any of the available GNSS systems, however; if the minor channel (f_T2) delivers similar data quality to the major one (f_T1), the special provisions for the former will practically have no effect.

A3.1 Adaptive low-pass filtering and minor-channel extrapolation

Before applying the dual-frequency ionospheric correction, the merged bending angle state profiles α_M,k(z_a) at the common z_a grid are filtered with further BWS low-pass filter operations, and the minor channel is extrapolated.

For α_M1 the filter is set to the same cutoff frequency as the basic BWS filter preceding the Doppler derivation (i.e., f_c1=2.5 Hz), ensuring a reliable reference resolution and basic smoothness of the whole merged profile. For filtering of α_M2 a (GPS L2) noise-minimization algorithm is used, following the approach of Sokolovskiy et al. (2009) for optimal filtering for ionospheric correction. We search for minimized noise employing a flexible cutoff frequency f_c2∈{2.5, 2, 10∕7, 1, 5∕7, 0.5 Hz}, corresponding to using cutoff periods τ_c2 from 0.4 to 2 s and sample widths of M=40 to M=200 (for BWS filter design details see Sect. A1.1).

We adopt the cutoff frequency f_c2 for α_M2 filtering that minimizes the noise fluctuations of the ionosphere-corrected atmospheric bending angle delta profile $\mathit{\delta }{\mathit{\alpha }}_{\mathrm{rm}}^{fc\mathrm{2}}\left(z_{\mathrm{a}}\right)={\mathit{\alpha }}_{\mathrm{r}}^{fc\mathrm{2}}\left(z_{\mathrm{a}}\right)-{\mathit{\alpha }}_{\mathrm{m}}\left(z_{\mathrm{a}}\right)$ when evaluated over the mesospheric altitude range between 50 and 70 km (similar to the functional minimization of Sokolovskiy et al., 2009; Eq. 4 therein). At these high altitudes the residual atmospheric mean signal after subtraction of the forward-modeled signal α_m(z_a) is very small (< 0.03–0.3 µrad), and therefore the noise level representative for the given RO event is well quantifiable.

The weight matrix of the BWS filter, ${\mathbf{A}}_k^{\mathrm{BWS}}$, is determined for both frequencies analogously to Eqs. (A3) to (A5). When using the baseband approach with model profile α_m to create the delta profile δα_Mm with elements

$$\begin{array}{}\text{(A40)}& \mathit{\delta }{\mathit{\alpha }}_{\mathrm{Mm}i,k}={\mathit{\alpha }}_{\mathrm{M}i,k}-{\mathit{\alpha }}_{\mathrm{m}i},\end{array}$$

the filtered bending angle is then (item 3.1 in Fig. 2)

$$\begin{array}{}\text{(A41)}& {\mathit{\alpha }}_{\mathrm{F}i,k}={\mathit{\alpha }}_{\mathrm{m}i}+\sum _{j=\mathrm{0}}^N{A}_{ij,k}^{\mathrm{BWS}}\cdot \mathit{\delta }{\mathit{\alpha }}_{\mathrm{Mm}j,k},\end{array}$$

where $i,j\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N\mathit{\}}$ and $k\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$.

Due to the stronger power of the L1 signal for (most of) the GPS satellites, the GPS signals of both frequencies are not of the same quality, and the L2 data (for those satellites where encrypted and hence power-degraded L2 signals are transmitted) do not reach down as far as the L1 data (i.e., z_aBot2>z_aBot1). If due to this reason α_F2 does not reach down as far as α_F1 and z_aBot2≤z_aBot2Max (with z_aBot2Max currently set to 15 km), a tropospheric bending angle extrapolation (TBAE) is applied in order to artificially extend α_M2 to also reach down to z_aBot1 (item 3.3 in Fig. 2).

Briefly summarized, this TBAE is currently implemented as follows. A linear gradient profile for the difference profile between the two bending angles, ${\mathit{\alpha }}_{\mathrm{F}\mathrm{12}}=({\mathit{\alpha }}_{\mathrm{F}\mathrm{1}}-{\mathit{\alpha }}_{\mathrm{F}\mathrm{2}})$, is estimated by a least squares fit over a sufficiently wide impact altitude range from z_aBot2 upward (as wide as the extrapolation range, at least 10 km). This gradient profile is then linearly extended down to z_aBot1 and subtracted from α_F1, to obtain the extrapolated part of α_F2 from z_aBot2 to z_aBot1. If z_aBot2>z_aBot2Max, then no TBAE is performed since the extrapolation range is considered too large. Details are provided by Kirchengast et al. (2017b), where the most recent version of the atmospheric bending angle derivation is described that includes this α_F12 extrapolation in a further advanced form.

For the propagation of the estimated random uncertainty we get (item 3.2 in Fig. 2)

$$\begin{array}{}\text{(A42)}& {\mathbf{C}}_{\mathit{\alpha }\mathrm{F},k}={\mathbf{A}}_k^{\mathrm{BWS}}\cdot {\mathbf{C}}_{\mathit{\alpha }\mathrm{M},k}\cdot {\left({\mathbf{A}}_k^{\mathrm{BWS}}\right)}^{\mathrm{T}}\end{array}$$

for the bending angle error covariance matrices of the leading (k=1) and minor (k=2) channel.

If a TBAE is applied to α_F2, the random uncertainty of α_F2 below z_aBot2 is equal to that of α_F1, because the noise is “copied” from α_F1 since the linear gradient profile from fitting α_F12 is noise-free. As a consequence, in these cases, we set the matrix elements of C_αF2 to (item 3.4 in Fig. 2)

$$\begin{array}{}\text{(A43)}& C_{\mathit{\alpha }\mathrm{F}\mathrm{2},ij}=\left\{\begin{array}{ll}{\mathbf{C}}_{\mathit{\alpha }\mathrm{F}\mathrm{2},ij}& \text{ for }{z}_{\mathrm{aTop}}>z_{\mathrm{a},i}>z_{\mathrm{aBot}\mathrm{2}}\\ & \text{ and }{z}_{\mathrm{aTop}}>z_{\mathrm{a},j}>z_{\mathrm{aBot}\mathrm{2}}\\ {\mathbf{C}}_{\mathit{\alpha }\mathrm{F}\mathrm{1},ij}& \text{ for }{z}_{\mathrm{aBot}\mathrm{2}}>z_{\mathrm{a},i}>z_{\mathrm{aBot}\mathrm{1}}\\ & \text{ and }{z}_{\mathrm{aBot}\mathrm{2}}>z_{\mathrm{a},j}>z_{\mathrm{aBot}\mathrm{1}}\\ \mathrm{0}& \text{ for }{z}_{\mathrm{aBot}\mathrm{2}}>z_{\mathrm{a},i}>z_{\mathrm{aBot}\mathrm{1}}\\ & \text{ and }{z}_{\mathrm{aTop}}>z_{\mathrm{a},j}>z_{\mathrm{aBot}\mathrm{2}}\\ \mathrm{0}& \text{ and }{z}_{\mathrm{aTop}}>z_{\mathrm{a},i}>z_{\mathrm{aBot}\mathrm{2}}\\ & \text{ for }{z}_{\mathrm{aBot}\mathrm{2}}>z_{\mathrm{a},j}>z_{\mathrm{aBot}\mathrm{1}}\end{array}\right..\end{array}$$

C_αF1 and C_αF2 can then be decomposed as needed into $u_{\mathit{\alpha }\mathrm{F}\mathrm{1}}^r$, R_αF1, and $u_{\mathit{\alpha }\mathrm{F}\mathrm{2}}^r$, R_αF2, respectively. Kirchengast et al. (2017b) describe the most recent version consistent with a further advanced form of the TBAE, where the separate assignments according to Eq. (A43) are no longer needed.

The estimated systematic uncertainties $u_{\mathit{\alpha }\mathrm{M},k}^s$ (in practice the basic and apparent systematic uncertainty estimates separately) are filtered with the same filter settings as for the state profiles (item 3.1 in Fig. 2) and are thus obtained in the form

$$\begin{array}{}\text{(A44)}& u_{\mathit{\alpha }\mathrm{F}i,k}^s=\sum _{j=\mathrm{0}}^N{A}_{ij,k}^{\mathrm{BWS}}\cdot u_{\mathit{\alpha }\mathrm{M}j,k}^s.\end{array}$$

Since these are smooth profiles, they are marginally changed by this low-pass filtering. The systematic uncertainty component contributed by the TBAE to the estimated systematic uncertainty is added after the ionospheric correction (see next subsection).

As for the basic low-pass filtering of excess phases (Sect. A1.1), the resolution profiles of the filtered bending angles w_αF1 and w_αF2 are determined by the cutoff frequencies f_c1 and f_c2 of the BWS filters, following Eqs. (A11) and (A13).

A3.2 Ionospheric correction

Based on the filtered and sometimes extrapolated state profiles α_F1 and α_F2, the ionospheric refractive effects are corrected for by the standard dual-frequency correction of bending angles (Vorob'ev and Krasil'nikova, 1994) used in the f_T1–f_T2 difference profile form (Sokolovskiy et al., 2009) (item 3.5 in Fig. 2). For the elements of the retrieved atmospheric bending angle profile α_r we thus get

$$\begin{array}{}\text{(A45)}& {\mathit{\alpha }}_{\mathrm{r},i}={\mathit{\alpha }}_{\mathrm{F}\mathrm{1},i}+{\mathit{\gamma }}_{f\mathrm{T}\mathrm{12}}\cdot \mathit{\delta }{\mathit{\alpha }}_{\mathrm{F}\mathrm{12},i},\end{array}$$

where

$$\begin{array}{}\text{(A46)}& \mathit{\delta }{\mathit{\alpha }}_{\mathrm{F}\mathrm{12},i}={\mathit{\alpha }}_{\mathrm{F}\mathrm{1},i}-{\mathit{\alpha }}_{\mathrm{F}\mathrm{2},i}\end{array}$$

and

$$\begin{array}{}\text{(A47)}& {\mathit{\gamma }}_{f\mathrm{T}\mathrm{12}}={\displaystyle \frac{f_{\mathrm{T}\mathrm{2}}^{\mathrm{2}}}{f_{\mathrm{T}\mathrm{1}}^{\mathrm{2}}-f_{\mathrm{T}\mathrm{2}}^{\mathrm{2}}}}.\end{array}$$

Propagated through the operator of the ionospheric correction (Eq. A45, currently used here in the classical form with f_T1 and f_T2 terms), the estimated random uncertainty of the resulting atmospheric bending angle, expressed by the error covariance matrix C_αr (item 3.6 in Fig. 2), is obtained as

$$\begin{array}{}\text{(A48)}& {\mathbf{C}}_{\mathit{\alpha }\mathrm{r}}={\left(\mathrm{1}+{\mathit{\gamma }}_{f\mathrm{T}\mathrm{12}}\right)}^{\mathrm{2}}{\mathbf{C}}_{\mathit{\alpha }\mathrm{F}\mathrm{1}}+{\mathit{\gamma }}_{f\mathrm{T}\mathrm{12}}^{\mathrm{2}}{\mathbf{C}}_{\mathit{\alpha }\mathrm{F}\mathrm{2}}.\end{array}$$

C_αr can then also be decomposed into $u_{\mathit{\alpha }r}^r$ and R_αr with the usual equations (cf., e.g., Eqs. A19 and A20).

Equation (A45) is as well applied to propagate the estimated systematic uncertainty (in practice the basic and apparent systematic uncertainty estimates separately) through the ionospheric correction using (item 3.5 in Fig. 2)

$$\begin{array}{}\text{(A49)}& u_{\mathit{\alpha }\mathrm{r},i}^s=u_{\mathit{\alpha }\mathrm{F}\mathrm{1},i}^s+{\mathit{\gamma }}_{f\mathrm{T}\mathrm{12}}\cdot \left(u_{\mathit{\alpha }\mathrm{F}\mathrm{1},i}^s-u_{\mathit{\alpha }\mathrm{F}\mathrm{2},i}^s\right),\end{array}$$

where it is assumed that the systematic errors in α_F1 and α_F2 are positively correlated, i.e., have the same sign, and the associated uncertainty estimates are hence subtracted from one another (as the bending angles are in Eq. A46). This assumption is reasonable, since the same sources of non-ionospheric systematic effects apply to both frequency channels (Doppler shift, orbit velocity, and orbit position uncertainties).

In the case of TBAE, Eq. (A49) needs to be supplemented below z_aBot2, since additional uncertainties $u_{\mathit{\alpha }\mathrm{2}\mathrm{TE}}^s$ arise from the errors made in the fitting parameters and in the extrapolation model (linear extrapolation) of the TBAE. Hence, for the range $z_{\mathrm{aBot}\mathrm{2}}>z_{\mathrm{a},i}\ge z_{\mathrm{aBot}\mathrm{1}}$,

$$\begin{array}{}\text{(A50)}& u_{\mathit{\alpha }\mathrm{r},i}^s=u_{\mathit{\alpha }\mathrm{r}}^s\left(z_{\mathrm{aBot}\mathrm{2}}\right)+u_{\mathit{\alpha }\mathrm{2}\mathrm{TE},i}^s,\end{array}$$

with $u_{\mathit{\alpha }\mathrm{2}\mathrm{TE}}^s$ being the conservative estimate for additional (apparent) systematic uncertainty within the extrapolated impact altitude range. We set $u_{\mathit{\alpha }\mathrm{2}\mathrm{TE}}^s$ to zero at z_aBot2 and linearly increase it from there downwards with a gradient of 1 µrad per 10 km (an experience-based best guess; cf. Scherllin-Pirscher et al., 2011b, a, who also address aspects of such tropospheric extrapolation in their discussions of error sources). It is interpreted as an apparent systematic uncertainty estimate, since due to the linear fit-based TBAE construction its event-to-event bias character will be essentially random (Scherllin-Pirscher et al., 2011b).

Also, the ionospheric correction currently applied in the rOPS is just a first-order correction, which will leave higher-order residual ionospheric errors in α_r (e.g., Syndergaard, 2000; Danzer et al., 2013; Liu et al., 2013, 2015; Healy and Culverwell, 2015). The uncertainty from higher-order residual ionospheric biases (RIBs), $u_{\mathrm{RIB}}^s$, is therefore added to the propagated (basic) systematic uncertainty. $u_{\mathrm{RIB}}^s$ is interpreted as basic systematic uncertainty, since the higher-order ionospheric residuals may not vanish in ensemble-of-events averaging. The other non-ionospheric sources of systematic errors and the RIBs can be reasonably assumed to be uncorrelated. The total estimated systematic uncertainty of the retrieved atmospheric bending angle α_r hence is

$$\begin{array}{}\text{(A51)}& u_{\mathit{\alpha }\mathrm{r},i}^s=\sqrt{{\left(u_{\mathit{\alpha }\mathrm{r},i}^s\right)}^{\mathrm{2}}+{\left(u_{\mathrm{RIB}}^s\right)}^{\mathrm{2}}}.\end{array}$$

Based on previous studies (e.g., Liu et al., 2013, 2015; Danzer et al., 2013, 2015), $u_{\mathrm{RIB}}^s$ is taken to be constant along the entire profile and is estimated to amount to 0.05 µrad. These last two components, $u_{\mathit{\alpha }\mathrm{2}\mathrm{TE}}^s$ and $u_{\mathrm{RIB}}^s$, are indicated as item 3.7 in Fig. 2. It is clear that this initial systematic uncertainty estimation can be significantly improved by future dedicated work on better quantifying and (if suitable) correcting for the systematic uncertainty components.

The resolution of the retrieved bending angle, w_αr, essentially corresponds to the higher resolution of the two bending angle profiles α_F1 and α_F2, and thus generally closely matches w_αF1 in most cases. As a simple but robust and suitable estimate, assuming that the resolutions of α_r and α_F1 scale in the same way as the correlation lengths l_αr and l_αF1 (derived from R_αr and R_αF1 as described for R_LF in Eq. A9), we compute w_αr as

$$\begin{array}{}\text{(A52)}& w_{\mathit{\alpha }\mathrm{r},i}={\displaystyle \frac{l_{\mathit{\alpha }\mathrm{r},i}}{l_{\mathit{\alpha }\mathrm{F}\mathrm{1},i}}}\cdot w_{\mathit{\alpha }\mathrm{F}\mathrm{1},i}.\end{array}$$

In concluding, we note that the atmospheric bending angle derivation algorithms used in this study – i.e., the adaptive filtering, TBAE, and ionospheric correction parts as described in this section – have recently received further advancement towards a form fully based on the combination of α_F1 and the difference profile α_F12 (rather than of α_F1 and α_F2), more aligned with the concept of Sokolovskiy et al. (2009). A detailed description of this most recent version is found in Kirchengast et al. (2017b).

Appendix B: Variance propagation for comparison

The full covariance propagation applied to propagate random uncertainties requires numerically “expensive” matrix operations, and therefore considerable efforts were made to seize opportunities for reducing the number of numerical operations (e.g., by only calculating with those elements of the band matrix A^BWS which lie within the width of the filter window).

However, as demonstrated in Sect. 4, simplification to a mere variance propagation (i.e., only considering the diagonal elements of the covariance matrices) is not reasonably possible because it leads to an unacceptable overestimation of random uncertainties. This overestimation occurs since the influence of the covariance elements – and thus for example the partially compensating impact of the negative side peaks in the correlation functions – is disregarded.

Here we state the two equations used to obtain the variances-only propagation results shown for comparison purposes in Fig. 9: the estimated random uncertainty was propagated through the BWS filter using

$$\begin{array}{}\text{(B1)}& {\left(u_{\mathit{\alpha }\mathrm{F}i,k}^r\right)}^{\mathrm{2}}=\sum _{j=\mathrm{0}}^N{\left(A_{ij,k}^{\mathrm{BWS}}\right)}^{\mathrm{2}}\cdot {\left(u_{\mathit{\alpha }\mathrm{M}j,k}^r\right)}^{\mathrm{2}},\end{array}$$

and subsequently through the ionospheric correction using

$$\begin{array}{}\text{(B2)}& {\left(u_{\mathit{\alpha }\mathrm{r},i}^r\right)}^{\mathrm{2}}={\left(\mathrm{1}+{\mathit{\gamma }}_{f\mathrm{T}\mathrm{12}}\right)}^{\mathrm{2}}\cdot {\left(u_{\mathit{\alpha }\mathrm{F}\mathrm{1},i}^r\right)}^{\mathrm{2}}+{\mathit{\gamma }}_{f\mathrm{T}\mathrm{12}}^{\mathrm{2}}\cdot {\left(u_{\mathit{\alpha }\mathrm{F}\mathrm{2},i}^r\right)}^{\mathrm{2}}.\end{array}$$