Wave-optics uncertainty propagation and regression-based bias model in GNSS radio occultation bending angle retrievals
- 1A.M. Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences, Moscow, Russia
- 2Hydrometeorological Research Centre of Russian Federation, 123242, Moscow, B. Predtechensky per., 11-13, Russia
- 3Wegener Center for Climate and Global Change (WEGC), University of Graz, Graz, Austria
- 4Institute for Geophysics, Astrophysics, and Meteorology–Institute of Physics, University of Graz, Graz, Austria
Correspondence: Michael E. Gorbunov (email@example.com)
A new reference occultation processing system (rOPS) will include a Global Navigation Satellite System (GNSS) radio occultation (RO) retrieval chain with integrated uncertainty propagation. In this paper, we focus on wave-optics bending angle (BA) retrieval in the lower troposphere and introduce (1) an empirically estimated boundary layer bias (BLB) model then employed to reduce the systematic uncertainty of excess phases and bending angles in about the lowest 2 km of the troposphere and (2) the estimation of (residual) systematic uncertainties and their propagation together with random uncertainties from excess phase to bending angle profiles. Our BLB model describes the estimated bias of the excess phase transferred from the estimated bias of the bending angle, for which the model is built, informed by analyzing refractivity fluctuation statistics shown to induce such biases. The model is derived from regression analysis using a large ensemble of Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) RO observations and concurrent European Centre for Medium-Range Weather Forecasts (ECMWF) analysis fields. It is formulated in terms of predictors and adaptive functions (powers and cross products of predictors), where we use six main predictors derived from observations: impact altitude, latitude, bending angle and its standard deviation, canonical transform (CT) amplitude, and its fluctuation index. Based on an ensemble of test days, independent of the days of data used for the regression analysis to establish the BLB model, we find the model very effective for bias reduction and capable of reducing bending angle and corresponding refractivity biases by about a factor of 5. The estimated residual systematic uncertainty, after the BLB profile subtraction, is lower bounded by the uncertainty from the (indirect) use of ECMWF analysis fields but is significantly lower than the systematic uncertainty without BLB correction. The systematic and random uncertainties are propagated from excess phase to bending angle profiles, using a perturbation approach and the wave-optical method recently introduced by Gorbunov and Kirchengast (2015), starting with estimated excess phase uncertainties. The results are encouraging and this uncertainty propagation approach combined with BLB correction enables a robust reduction and quantification of the uncertainties of excess phases and bending angles in the lower troposphere.
The bending angle (BA) and atmospheric profiles retrieval chain for Global Navigation Satellite System (GNSS) radio occultation (RO) data includes many steps involving linear and (moderately) nonlinear transformations, starting from excess phase and amplitude measurements (Gorbunov et al., 2006). Error or uncertainty propagation through the geometric optical part of the retrieval chain has been investigated in a series of theoretical and empirical studies (Kursinski et al., 1997; Syndergaard, 1999; Palmer et al., 2000; Rieder and Kirchengast, 2001; Kuo et al., 2004; Steiner and Kirchengast, 2005; Schreiner et al., 2007; Scherllin-Pirscher et al., 2011b; Scherllin-Pirscher et al., 2011a; Scherllin-Pirscher et al., 2017; Innerkofler et al., 2016; Schwarz et al., 2016; Schwarz et al., 2017a; Schwarz et al., 2017b; Li et al., 2016; Li et al., 2017).
The uncertainty propagation through the wave-optical bending angle retrieval block was investigated recently for large-scale (systematic) and small-scale (random) uncertainties by Gorbunov and Kirchengast (2015), including simulation results demonstrating random uncertainty propagation. Such wave-optical retrieval is essential in the lower troposphere (altitudes below 5 km), where the RO observations are subject to several specific uncertainties not present higher up in the atmosphere, including effects from low signal-to-noise ratio, multipath propagation, and super-refraction (Sokolovskiy, 2001, 2003; Xie et al., 2006; Ao, 2007; Xie et al., 2010; Sokolovskiy et al., 2010).
A thorough treatment of systematic uncertainty and its propagation from excess phase to bending angle in the lower troposphere, including the aim to correct for the known boundary layer bias (BLB) in standard lower troposphere RO retrievals, often termed the “negative refractivity bias” (Sokolovskiy et al., 2010; Gorbunov et al., 2015), is lacking so far. Also, the propagation of both estimated systematic and estimated random uncertainties through the wave-optical chain, complementary to the geometric-optical uncertainty propagation work of Schwarz et al. (2017a, b), has not been investigated and demonstrated yet. This study focuses on providing these missing investigations and on demonstrating BLB correction for a representative large ensemble of real RO data from the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) mission as well as introducing a complete uncertainty propagation approach. The findings and algorithms obtained are used in the development of the new reference occultation processing system (rOPS) including an RO retrieval chain with integrated uncertainty propagation (Kirchengast et al., 2015, 2016a, b).
Our starting points for the BLB model construction are the approach based on refractivity fluctuations introduced by Gorbunov et al. (2015) and the recent study of RO systematic errors by Gorbunov (2014). Refractivity fluctuations constitute an external factor that results in a systematic shift in the signal phase due to its physical nature rather than any deficiency of the processing algorithm. Although this model cannot be looked at as a complete explanation of the bias, it serves as a convenient structural model that allows exposing probable candidates for the role of the objective characteristics of the signal received that may correlate with the bias. These characteristics will hereafter be referred to as predictors in the BLB model. In particular, it has already been shown by Gorbunov (2014) that the bending angle can serve as such a predictor. Further predictors and the complete BLB model setup based on a regression-modeling approach are described in this study.
This approach results in the BLB and (residual) systematic uncertainty model formulated in terms of tropospheric bending angles. In order to incorporate this uncertainty modeling into the RO retrieval chain with integrated uncertainty propagation, it needs to be transferred into the equivalent excess phase BLB and (residual) systematic uncertainty estimate. For its propagation, a perturbation approach or the approximation derived by Gorbunov and Kirchengast (2015) can be employed. In that paper we discussed the propagation of excess phase to bending angle uncertainty through the Fourier integral operator (FIO) used for the bending angle retrieval (Gorbunov and Lauritsen, 2004). This uncertainty propagation uses the stationary phase approximation, which allowed for the derivation of simple propagation formulae.
In order to now transform the bending angle uncertainty into the equivalent excess phase uncertainty, we use the inverse FIO, which was recently employed by Gorbunov (2016) for the retrieval of reflected rays from RO data. Specifically, the systematic uncertainty is evaluated for every RO event in the form of estimated profiles of bending angle BLB and (residual) systematic uncertainty. These estimates are then transformed into the equivalent BLB and (residual) systematic uncertainty of the excess phase, where they complement the estimated random and basic systematic uncertainty of the excess phase, which is available separately from the preceding step of excess phase processing (Innerkofler et al., 2016; Schwarz et al., 2017a, b). Both together are used as input to the wave-optical uncertainty propagation.
The paper is organized as follows. In Sect. 2 we describe the empirical BLB model, based on a regression analysis guided by the understanding that refractivity fluctuation statistics induce such biases, as well as a simple (residual) systematic uncertainty model for the BLB-corrected bending angles. Section 3 describes the wave-optical propagation of estimated systematic and random uncertainties from excess phase to bending angle for the methodology also recalling the key results needed from Gorbunov and Kirchengast (2015) and Gorbunov (2016). In Sect. 4 we discuss the results of the application of the BLB correction based on a large ensemble of COSMIC RO data from representative test days throughout the year 2008. Section 5 provides our conclusions.
The BLB model is formulated to be capable of providing bending angle BLB profiles over the lower troposphere up to 5 km impact altitude, corresponding to about 4 km (mean-sea-level) altitude, with the primary bias effects occurring within the atmospheric boundary layer below about 2 km altitude. Here we describe its setup by first introducing the underlying refractivity fluctuations model (Sect. 2.1), which is then followed by the BLB model description (Sect. 2.2). The model is built as a regression model using adaptive functions based on predictors available for each RO event, including impact altitude, latitude, bending angle, BA standard deviation, canonical transform (CT) amplitude, and the CT amplitude fluctuation index as the main ones. The selection of the predictors is explained in Sect. 2.3 and their use in constructing the adaptive functions in Sect. 2.4.
Along with the description we illustrate the performance of the BLB model to quantify the boundary layer biases based on the predictors, underpinning that the BLB profiles obtained for individual RO events can be effectively used for BLB correction and lead to a significant reduction of systematic uncertainty. A simple model for the estimated residual systematic uncertainty after the BLB profile subtraction, which accounts for the residual bias and the uncertainty (indirectly) incurred from the use of ECMWF analysis profiles as regression reference, is described in Sect. 2.5.
2.1 Underlying model of refractivity fluctuations
In order to formulate our approach to the bending angle BLB in terms of the negative refractivity bias (Sokolovskiy et al., 2010), we use the fluctuation-based model introduced by Gorbunov et al. (2015). This model is used as a simple structural model that allows finding good candidates for the objective characteristics of the observed signals that correlate with the bias. Figure 1 shows an example profile of the refractivity structure constant and the corresponding relative difference statistics of an ensemble of bending angle and refractivity profiles. The latter were obtained by comparison of the modeled “truth” based on ECMWF refractivity fields, used as reference, and perturbed data based on the same ECMWF fields but with random refractivity fluctuations according to the profile superimposed. The profile was tuned to realistically represent BLB statistics of RO observations and the wave-optics propagator (WOP) package (Gorbunov, 2011) was used to realistically compute the bending angles.
It is visible in Fig. 1 that the refractivity fluctuations lead to a negative refractivity bias of up to about 2 % in the boundary layer and an associated negative BLB in bending angle of up to about 5 %, which is typical of biases seen in real RO data. Random differences (standard deviation) reach realistic values as well – about 1.5 % in refractivity and about 5 % in bending angle.
To put these simulation results into direct context with real data, Fig. 2 shows another set of difference statistics for bending angles and refractivities, from low latitudes to high latitudes, where we again used the modeled truth from ECMWF fields as reference but now to illustrate the differences of observed profiles from COSMIC. The COSMIC data were processed by the OCC (occultations) package for RO data processing, as described in Gorbunov et al. (2006). These results confirm that the refractivity fluctuations model, with corresponding settings, can reproduce the systematic and random error behavior of RO bending angles and refractivities in the boundary layer. A somewhat higher level of RMS deviations (standard deviation) seen for the COSMIC data, compared to Fig. 1, is likely caused by the fact that ECMWF fields themselves deviate from the real atmospheric state (see, e.g., the error modeling performed by Scherllin-Pirscher et al., 2011b, 2017).
Figure 3 presents a latitude–longitude map of the COSMIC–ECMWF refractivity differences at a height of 0.6 km in terms of systematic relative refractivity deviation. These results illustrate the regional variations in refractivity bias behavior and are similar to those presented in Xie et al. (2006, 2010) and Gorbunov (2014).
Our further strategy of the bias correction consists in the following. We perform the numerical simulation of occultation events with superimposed fluctuations and analyze different objective characteristics of RO signals in order to find those that correlate with the simulated bias. These characteristics will be referred to as predictors. Using this set of predictors, we also compare the simulation results with the processing of real COSMIC observations. We assume that this will allow the formulation of a model for BLB correction, which will also effectively mitigate biases in the retrieved refractivity profiles and further-derived atmospheric profiles. We have to formulate the BLB model with a flexible functional behavior in order to reliably serve its purpose.
2.2 Bending angle BLB model from regression to adaptive functions
We model the BLB by a predictor-based empirical model that is flexible enough to capture the BLB behavior by suitable predictors under widely variable predictor value ranges for individual RO events. Because the dependence of the BLB model profiles from predictors is unknown a priori, we solve for this dependence in the form of the linear combination of a set of linear and nonlinear functions of the predictors. We refer to these functions as adaptive functions. The model estimate of the regression coefficients of the linear combination is based on the comparison of a large set of bending angle observations with a reference data set.
In this study, introducing a first reliable BLB model version, the observations are from the COSMIC mission and the reference data set consists of gridded fields of meteorological variables from ECMWF. The ECMWF data have their own systematic uncertainty, which is taken into account by letting these uncertainties flow into the estimated residual systematic uncertainty of bending angle profiles after BLB correction (Sect. 2.5).
The BLB model is formulated as follows. We used a set of COSMIC bending angle observations, including 24 representative days from year 2008. We adopted the 15th and 16th day of every month, amounting in total to about 54 000 RO events. We used the corresponding ECMWF fields as basis for obtaining the “true” reference bending angles. To this end, we employed the wave-optics propagator (Gorbunov, 2011) to generate the bending angle profiles from the ECMWF refractivity fields. We then performed a regression of the differences of observed and reference bending angles in the lower troposphere with respect to the chosen adaptive functions (Sect. 2.4). The adaptive functions are formulated in terms of predictors, which are evaluated from objective characteristics of every RO event, without using the reference data (Sect. 2.3). These ingredients allow for the derivation of regression coefficients, which upon their estimation complete the BLB model then ready to be applied based on predictors from a given RO event.
Because we need to derive the regression model for widely diverse BLB behavior, we start with very general regression relations. Consider two series of random variables, vector xi and scalar series yi, where the lower index i enumerates the realizations. We will define the components of xi predictors, because we approximate the random variables yi as a linear combination of predefined adaptive functions φj of xi. The number of predictors, and of associated adaptive functions, is much smaller than the number of realizations (difference profiles of observed and reference bending angles in the lower troposphere). We write the overdetermined system of equations as
or in the vector form,
This system has a pseudo-inverse solution, i.e., the vector α that minimizes the discrepancy
is obtained as the least-squares solution of this overdetermined problem in the form
Now consider a numerical estimation of α that allows for an evaluation readily augmentable in terms of the number of realizations and adaptive functions. Preparing the quadratic form
we have available matrix as a square symmetric matrix that can be evaluated by the summation over any existing set of realizations of xi. Similarly, using the transform
we have available vector z as a vector that can also be evaluated by the summation over any existing set of realizations of xi and yi. Finally, it is straightforward in this formulation to obtain the regression coefficients from
For convenience, matrix and vector z can be redefined in terms of averaging over the ensemble of realizations. Denoting N as the number of realizations, this is performed by dividing both and z by N,
Practically, normalization can also be an issue, depending on the number of adaptive functions. If their number is as high as about 200, such as in our study (Sect. 2.4), then even a small change in the normalization factor is raised to the 200th power when evaluating the matrix determinant. This may result in overflow or underflow in the matrix inversion. Therefore, the numerical algorithm requires accurate tuning of the normalization factor in order to ensure a stable and robust inversion of matrix .
After having solved for the regression coefficient vector α, it can be used within Eq. (3), which then serves as the BLB model applicable to any given RO event. It will provide the estimated bending angle BLB profile y for the RO event when its predictors are used to specify the model matrix .
2.3 Predictors for the model's adaptive functions
Here we consider the predictors that we may reasonably choose. Besides predictors depending on RO event altitude and latitude (discussed separately below), we adopt the following four predictors that are derived from observational RO data – all as function of impact parameter p within the lower troposphere (below an impact altitude of 4.5 km): (1) bending angle ϵ(p), (2) bending angle standard deviation δϵ(p), (3) normalized CT amplitude ACT(p), and (4) CT amplitude fluctuation index β(p). Bending angle standard deviation is the bending angle standard error estimate based on radio-holographic analysis (Gorbunov et al., 2006). The CT amplitude (Gorbunov, 2002; Gorbunov and Lauritsen, 2004) is the normalized energy distribution over rays in the impact parameter space. We use the CT amplitude normalized in such a way that it should equal unity in vacuum. The CT amplitude fluctuation index β(p) is defined as
where is a smoothing operator (low-pass filter) for which we use a 2 km smoothing width.
Figure 4 shows the scatter plot of fluctuation-affected bending angles versus reference bending angles for the fluctuation-model simulations (like for Fig. 1) and the COSMIC observations (like for Fig. 2). In both cases the asymmetry with respect to the diagonal is visible (fluctuation-affected bending angles are tentatively smaller than reference ones). This indicates that the bending angle itself can serve as one meaningful predictor of (negative) boundary layer biases.
Figure 5 shows scatter plots of the difference of fluctuation-affected and reference bending angle profiles versus bending angle standard deviation (top), normalized CT amplitude (middle), and CT amplitude fluctuation index (bottom) for simulations (left) and COSMIC observations (right).
Comparing the behavior of these predictors, their correlation with the bending angle difference is clearly more salient in the simulations but some smaller asymmetry can also be noticed for the COSMIC observation differences. We therefore kept all four predictors in this study and left possible further reduction of these predictors (and associated adaptive functions) to future fine-tuning of the BLB model regression. An important conclusion from these comparisons is that the fluctuation model alone does not explain the patterns observed in the real observations. However, the role of this model is to help in finding reasonable predictors. The further bias correction procedure is only based on the predictors that can be readily derived from observations, rather than on the fluctuation model.
In addition to these four predictors we utilize the RO event coordinates (impact altitude z and latitude λ), where , with RLC being the local radius of curvature and Ugeoid the geoid undulation applying to the event location. We use the impact altitude z directly and in the form of the following six trigonometric functions of z:
where zmin and zmax are the limits of impact altitude wherein the BLB profiles are evaluated (equal to 1.5 and 4.5 km). Latitude λ is used in the form of another six trigonometric functions of λ,
Altogether we therefore use Np=17 predictors, including impact altitude, the four observation-derived predictors, six functions of impact altitude, and six functions of latitude. This number of predictors exceeds that in radiation correction schemes, where six to eight predictors are typically used (e.g., Zhu et al., 2014).
2.4 Construction of the model's adaptive functions
General adaptive functions as we use here are constructed in the form of different degrees of the predictors and their cross products from degree zero, which produces unity, up to some maximum degree Dp,
We use a maximum degree of Dp=3 and apply some additional constraints in order to reduce the number of adaptive functions. For the six trigonometric functions of impact altitude (Eq. 14), taking their degrees beyond degree 1 and their cross products is not allowed as these will not be linearly independent from other trigonometric functions of the impact altitude. The same applies to the six trigonometric functions of latitude (Eq. 15) for which we therefore also disregard degrees beyond degree 1 and cross products.
For our choice of Dp=3 we thus obtain Nf=214 adaptive functions. To understand this number, consider different degrees of predictors. There is one zero-degree function (unity). There are 17 functions of degree 1 (the 17 predictors). There are functions of degree 2. There are functions of degree 3. Therefore, we arrive in total at adaptive functions, which provide the needed flexibility for the highly variable BLB profile behavior while still allowing for a robust estimation of the regression coefficients. If future fine-tuning of the regression model were to reduce the number of predictors, the number of adaptive functions would reduce accordingly.
2.5 Simple residual systematic uncertainty model
As described in Sect. 2.2, after obtaining the regression coefficient vector (Eq. 10), we can use it within the regression model (Eq. 3), which then serves as the BLB model applicable to any given RO event. It provides the bending angle BLB model profile for the RO event, δαBLB(z), based on its predictors depending on location (impact altitude, latitude) and bending angle and CT amplitude characteristics (Sect. 2.3).
Given this basis, we define a simple initial systematic uncertainty model for the BLB-corrected bending angle profiles of the lower troposphere, , which consists of two components: (1) an estimated lower-bound ECMWF-reference field-induced systematic uncertainty, , that accounts for the uncertainty from using the ECMWF analysis fields as the regression reference, which have their own (small) systematic deviations from the truth, and (2) an estimated residual bias uncertainty after BLB correction by subtracting the BLB model profile, , since the empirical–statistical BLB regression model can never fully fit the individual bias situation of an RO event.
From experience with estimated biases of ECMWF analysis fields in other studies (e.g., Li et al., 2013, 2015; Scherllin-Pirscher et al., 2017; Li et al., 2017), we formulate the model for the ECMWF-reference field-induced systematic uncertainty profile as a fractional model () with a linear increase downward over the lower troposphere towards the surface,
where αrefEC(z) is the ECMWF reference bending angle profile, zmin and zmax are the limits of impact altitude (set to 1.5 and 5.0 km), and is the fractional uncertainty at zmin empirically set to 0.25 %. For perspective, the bending angle uncertainties obtained this way correspond in terms of temperature to uncertainties from about 0.2 K near 4 km impact altitude to 0.6 K near the surface (for details on uncertainty relations among RO-derived variables see Scherllin-Pirscher et al., 2011b, 2017, and references therein).
The estimated residual bias uncertainty profile after BLB correction is formulated from experience with other bias corrections, such as sampling bias correction (e.g., Scherllin-Pirscher et al., 2011a, 2017), and based on BLB correction performance results with test ensembles during this study in a straightforward fractional form,
where rresBLB is the systematic uncertainty reduction factor empirically set to 0.2, i.e., expressing that due to the BLB correction the bias in the bending angle profile is reduced by a factor of 5.
For the estimated residual systematic uncertainty finally attributed to the BLB-corrected lower-tropospheric bending angle at any impact altitude, we then simply adopt the larger one of the two uncertainties,
implementing the lower-bound uncertainty role of in cases where the estimated residual bias uncertainty of individual RO events according to Eq. (20) is occasionally very small. is a smoothing operator (low-pass filter) with a 0.4 km filter width that we use to ensure adequate smoothness of the resulting profile also over those altitude levels where the two uncertainty components cross in their magnitude.
The propagation of systematic and random uncertainties through the wave-optical retrieval chain was investigated by Gorbunov and Kirchengast (2015), where a simple approximation was derived and verified based on numerical simulations (as summarized in Sect. 1). The approximation considers the excess phase as function of time, Ψ(t), and its systematic (small-scale) and random (large-scale) uncertainties, Σ1(t) and Σ2(t), respectively. The uncertainty in the impact parameter space (Gorbunov and Lauritsen, 2004) is then evaluated as , where t(p) is the time of observation of the ray with impact parameter p.
Practically the application of this approximation was shown by Gorbunov and Kirchengast (2015) to work well for propagating random uncertainties (covariance matrices), while in sensitivity tests and evaluations for this study we found that it does not work sufficiently well for propagating systematic uncertainties, due to the large-scale nature of such (increment) profiles not transforming smoothly under FIO operations (Gorbunov and Lauritsen, 2004). Similarly, given the BLB and residual systematic uncertainty model being formulated in terms of bending angle, the inverse transformation of these into the equivalent excess phase bias and uncertainty proves to be not straightforward either.
The reason and underlying problem is that the perturbation of the excess phase due to the superimposing of the systematic uncertainty of the bending angle is not smooth. The variation in the bending angle profile in each realization results in the different phase perturbation corresponding to a different ray manifold with a different caustic structure. Therefore, the excess phase perturbation has a complicated nonlinear relation with the phase (eikonal) uncertainty in impact parameter space, and this perturbation corresponds to a complicated coherent signal being a superposition of multiple signals corresponding to different rays.
To overcome this difficulty, we do apply the linearized approximation only for the propagation of random uncertainty, i.e., the covariance propagation according to Gorbunov and Kirchengast (2015) – Eqs. (29) and (30) therein. This is applied within the rOPS wave-optical retrieval, for both GNSS frequencies, right after the bending angle profiles themselves have been retrieved by the (forward) FIO in the CT2 implementation (Gorbunov and Lauritsen, 2004; Gorbunov, 2011). The BLB and estimated systematic uncertainty propagation is then computed, in a consistent way for bending angles and excess phases, with a perturbation approach in a three-step sequence as follows.
First, the BLB profile and its estimated systematic uncertainty profile after BLB subtraction are computed according to Sect. 2.5 for the lower-tropospheric bending angle profile at the L1 frequency, for the location and characteristics (i.e., the applicable predictors) of the given RO event. It is not computed for the second (L2) frequency, since the L2 profiles are generally more noisy (making BLB estimation difficult) and not further used at impact altitudes below 5 km. Below this level, where the neutral atmospheric excess phase always exceeds several hundreds of meters, the dual-frequency ionospheric correction instead always uses L1–L2 difference bending angles extrapolated from above (Schwarz et al., 2017b), avoiding noise amplification and mitigating potentially adverse effects on top-of-boundary-layer estimates recently pointed out by Sokolovskiy et al. (2016).
Second, the BLB-corrected L1 bending angle profile and this profile perturbed by the estimated systematic uncertainty profile are each projected back to excess phase by applying the inverse FIO approach recently introduced by Gorbunov (2016). This provides the BLB-corrected L1 excess phase profile and, from the difference of the two back-projected profiles, the estimated systematic excess phase uncertainty profile pertaining to it. The latter BLB-related systematic uncertainty is then added (in a root-mean-square sense) to the basic systematic excess phase uncertainty available from the raw processing towards excess phase (Innerkofler et al., 2016), yielding the total estimated systematic excess phase uncertainty profile.
Third, the BLB-corrected L1 excess phase profile and this profile perturbed by the total estimated systematic uncertainty profile are processed again through the standard (forward) FIO CT2-wave-optics retrieval in order to obtain a BLB-corrected retrieved bending angle profile, for a consistency check with the original BLB-corrected bending angle profile, as well as the total estimated systematic bending angle uncertainty profile, from the difference of the two CT2-retrieved bending angle profiles. The systematic bending angle uncertainty profile at the second (F2) frequency is finally obtained from also processing the L2 excess phase profile perturbed by its associated systematic uncertainty through the wave-optics retrieval and estimating it from the difference of the resulting perturbed bending angle profile to the one originally retrieved from the unperturbed L2 excess phase.
Despite of the complexities from the nonlinearities involved, we obtain in this way a consistent set of excess phase and bending angle profiles together with their estimated systematic and random uncertainties, which are BLB-corrected at the L1 frequency in the lower troposphere. The extra computational expense for the uncertainty propagation due to the nonlinearity is reasonably limited to one additional forward and inverse FIO operation at L1 frequency, which is required for the perturbation approach to systematic uncertainty propagation. This is similar to the uncertainty propagation work of Schwarz et al. (2017a, b), where the perturbation approach is also needed in a small number of steps (during geometric-optics bending angle retrieval and dry-air temperature retrieval) for the systematic uncertainty propagation.
Here we evaluate the consistency of the BLB-corrected bending angles and their associated retrieved refractivities by (i) using the original BLB-corrected bending angles and (ii) back-projecting the original BLB-corrected retrieved bending angles to obtain BLB-corrected excess phases and then retrieving the bending angles again. This provides a basic validation of our procedure as described in Sect. 3; to limit the extent of this paper, the detailed inspection and validation of the uncertainty propagation itself is left for a follow-on study.
We investigated the BLB correction of an independent ensemble of COSMIC-retrieved bending angles employing our BLB model, as in Sect. 2.1 using ECMWF analysis fields as reference. Figure 6 shows the COSMIC–ECMWF difference statistics of bending angles (left) and refractivities (right) after bending angle BLB correction. These statistics were evaluated for a set of 12 days of COSMIC data from year 2008, including the 17th day of every month, amounting in total to about 26 000 RO events. This implies that these COSMIC and ECMWF ensembles of profiles are independent from the ones used in the derivation of the BLB model regression coefficients (the 15th and 16th day of every month; see Sect. 2.1).
Cross-checking these results with results from COSMIC and ECMWF ensembles using the 16th day of every month (not separately shown), we find them practically indistinguishable in terms of their difference statistics. This indicates the statistical homogeneity of the data sets and the robustness of the BLB model. Furthermore, from comparing Fig. 6 with Fig. 2, it is clear that the BLB correction achieves a substantial decrease in the boundary layer biases by about a factor of 5, which is consistent with the systematic uncertainty reduction factor rresBLB=0.2 (Eq. 20). Immediately above the boundary layer, above about 2 km altitude, the BLB-corrected profiles possibly contain slightly increased uncertainty, at a small magnitude, which is accounted for by the reference field-induced lower-bound uncertainty (Eq. 19) included in the systematic uncertainty model up to 5 km impact altitude. This may be improved in the future by a further-refined BLB model design.
Figure 7 shows the COSMIC–ECMWF difference statistics of bending angles (left) and refractivities (right) based on the BLB-corrected bending angles obtained by back-projecting the original BLB-corrected bending angles by the inverse FIO to obtain BLB-corrected excess phases and then retrieving the bending angles again. Except for about the lower half of the kilometer above surface where there is possibly some degradation, the results are found very close to those shown in Fig. 6 for the original BLB-corrected bending angles. This indicates the basic validity and robustness of our approach to transfer the BLB-corrected bending angles to BLB-corrected excess phases. Future more detailed inspection of the full uncertainty propagation approach according to Sect. 3 will consolidate this encouraging initial validation.
Figure 8 presents a latitude–longitude map of the COSMIC–ECMWF difference at a height of 0.6 km in terms of systematic relative refractivity deviation, similar to Fig. 3, but after the BLB correction applied to the underlying bending angles. This plot indicates that, although the overall average bias is minimized, there are some regional maxima and minima. Some of them correspond to the areas with a sharp marine boundary layer (Xie et al., 2006, 2010; Gorbunov, 2014), where the negative bias is reduced but still remains. Other regions with larger deviations are located above northern Africa and Australia, where there is a positive overcorrection. The latter regions correspond to a similar terrain type, i.e., dry desert areas. This indicates the need for refined predictors, taking into account such regional effects, in order to further mitigate in a next step these more specific biases.
In this study we developed a regression-based approach for modeling and propagating the atmospheric boundary layer biases and associated (residual) systematic uncertainties within the wave-optical retrieval chain of the reference occultation processing system, which is a new RO processing system with integrated uncertainty propagation that focuses on calibration–validation and climate applications.
Currently, there is no quantitative physical model describing BLB in RO data, although there was a series of studies discussing different mechanisms resulting in BLB. The starting point encouraging and informing our BLB model design was the fluctuation-based explanatory modeling of the well-known negative refractivity bias problem in the boundary layer. We showed that it is possible to achieve a reasonable agreement with observed bending angle and refractivity biases by modeling fluctuation statistics consistent with reasonable tropospheric profiles of the refractivity structure constant .
Based on this understanding we can robustly assume that reliable modeling of the bending angle BLB, and subsequent use of the model for BLB correction, will also effectively mitigate biases in the retrieved refractivity profiles and further-derived atmospheric profiles. However, given the highly variable refractivity fluctuations affecting individual RO events in reality, which implies a complex dependence of the bending angle BLB on the location and the data characteristics of individual RO profiles, we found it necessary to implement a BLB model with a very flexible functional behavior in order to reliably serve its purpose.
We therefore have chosen a versatile empirical regression-modeling approach and found suitable predictors of the BLB in the lower-tropospheric bending angle, including the bending angle and its standard deviation, CT amplitude and its fluctuation index, impact altitude and its trigonometric functions, and trigonometric functions of latitude. Degrees and cross products of these predictors were used to form a set of flexible adaptive functions that served as the basis for the BLB model, which was then obtained by regression to a large ensemble of COSMIC and ECMWF profile differences. Also, a simple (residual) systematic uncertainty model was formulated, applying to the bending angles after BLB correction. For any given RO event, the BLB model profile can be computed based on the predictors that purely depend on the event location and the characteristics of the bending angle and CT amplitude profiles.
Together with the linearized wave-optics (random) uncertainty propagation approach described by Gorbunov and Kirchengast (2015), we used the new approach to formulate the algorithmic sequence for wave-optical retrieval of bending angles from excess phases including consistent BLB correction and the associated random and systematic uncertainty propagation. Evaluating the consistency of the BLB-corrected bending angles and their associated retrieved refractivities we achieved a successful basic validation of the new procedure: we found that the BLB correction delivers a substantial decrease in the boundary layer biases by about a factor of 5, which is consistent with our initial model of residual systematic uncertainty.
Our bias model uses ECMWF fields as a reference; therefore, it involves the biases that are intrinsic to these ECMWF fields. However, the same approach can be applied together with an independent estimate of the ECMWF biases. In this study, we assumed that ECMWF biases form a small fraction of the observed systematic COSMIC–ECMWF differences.
These results are encouraging for follow-on work in the near future that can provide a refined BLB model design and a detailed inspection and validation of the complete wave-optical retrieval and uncertainty propagation as introduced in this study. In this way, the rOPS geometric-optical bending angle retrievals (Schwarz et al., 2017b), generally available reliably from the middle troposphere upwards, can be complemented and merged, from the upper troposphere downwards, with these wave-optical bending angle retrievals. Jointly this provides high quality of the RO data and their integrated uncertainty estimates from the stratosphere down close to the surface.
The code used in this study does not belong to the public domain and cannot be distributed.
COSMIC radio occultation data are freely available. To get access to them, it is necessary to sign up at the website of the COSMIC Data Analysis and Archiving Center (CDAAC): http://cdaac-www.cosmic.ucar.edu/cdaac/ (follow the “Sign up” link for further details). ECMWF analyses are not free products and can only be obtained subject to licensing conditions depending on country and other factors. Information about ECMWF datasets and availability from the archive is provided at http://www.ecmwf.int/en/forecasts/accessing-forecasts; the commercial catalogue can be found at http://www.ecmwf.int/en/forecasts/datasets/catalogue-ecmwf-real-time-products.
Both authors formulated the initial approach to integrating wave-optical uncertainty propagation into the reference occultation processing system (rOPS) and the overall study design. MG conceived and developed the bias-modeling approach and the boundary layer bias model, performed the computational work and the analysis, prepared the figures, and wrote the first draft of the manuscript. GK provided input on the bias and uncertainty modeling design and feedback during the work and significantly contributed to the writing of the manuscript. Both authors contributed to consolidating the manuscript for submission and publication.
The authors declare that they have no conflicts of interest.
This article is part of the special issue “Observing Atmosphere and Climate with Occultation Techniques – Results from the OPAC-IROWG 2016 Workshop”. It is a result of the International Workshop on Occultations for Probing Atmosphere and Climate, Leibnitz, Austria, 8–14 September 2016.
Work on Sects. 1 and 2 was supported by the Russian Foundation for Basic
Research (grant No. 16-05-00358-a). Work on Sects. 3 and 4 was supported by
the Austrian Research Promotion Agency FFG within the Austrian Space
Applications Programme ASAP (ASAP-9 project OPSCLIMPROP and ASAP-10 project
OPSCLIMTRACE). We acknowledge Taiwan's National Space Organization (NSPO) and
the University Corporation for Atmospheric Research (UCAR) for providing the
COSMIC RO data via the COSMIC Data Analysis and Archiving Center (CDAAC). We
acknowledge the European Centre for Medium-Range Weather Forecasts (ECMWF)
for providing the atmospheric analysis fields.
Edited by: Axel von Engeln
Reviewed by: two anonymous referees
Gorbunov, M.: Statistical analysis of systematic errors in RO measurements, ROM SAF Visiting Scientist Report 20, Danish Meteorological Institute, Copenhagen, available at: http://www.romsaf.org/Publications/reports/romsaf_vs20_rep_v11.pdf (last access: 20 December 2017), (SAF/ROM/DMI/REP/VS20/001), 2014. a, b, c, d
Gorbunov, M.: Development of wave optics code for the retrieval of bending angle profiles for reflected rays, ROM SAF CDOP-2, Visiting Scientist Report 27, Danish Meteorological Institute, European Centre for Medium-Range Weather Forecasts, Institut d'Estudis Espacials de Catalunya, Met Office, available at: http://www.romsaf.org/Publications/reports/romsaf_vs27_rep_v10.pdf (last access: 20 December 2017), (SAF/ROM/DMI/REP/VS27/001), 2016. a, b, c
Gorbunov, M. E. and Kirchengast, G.: Uncertainty propagation through wave optics retrieval of bending angles from GPS radio occultation: Theory and simulation results, Radio Sci., 50, 1086–1096, https://doi.org/10.1002/2015RS005730, 2015. a, b, c, d, e, f, g
Gorbunov, M. E. and Lauritsen, K. B.: Analysis of wave fields by Fourier integral operators and its application for radio occultations, Radio Sci., 39, RS4010, https://doi.org/10.1029/2003RS002971, 2004. a, b, c, d, e
Gorbunov, M. E., Lauritsen, K. B., Rhodin, A., Tomassini, M., and Kornblueh, L.: Radio holographic filtering, error estimation, and quality control of radio occultation data, J. Geophys. Res., 111, D10105, https://doi.org/10.1029/2005JD006427, 2006. a, b, c
Gorbunov, M. E., Vorob'ev, V. V., and Lauritsen, K. B.: Fluctuations of refractivity as a systematic error source in radio occultations, Radio Sci., 50, 656–669, https://doi.org/10.1002/2014RS005639, 2015. a, b, c
Innerkofler, J., Pock, C., Kirchengast, G., Schwärz, M., Jäggi, A., Andres, Y., Marquardt, C., and Schwarz, J.: SI-traceable radio occultation excess phase processing with integrated uncertainty estimation for climate applications, in: OPAC-IROWG 2016 International Workshop, Seggau Castle, Austria, 8–14 September 2016, 2016. a, b, c
Kirchengast, G., Schwärz, M., Schwarz, J., Scherllin-Pirscher, B., Pock, C., Innerkofler, J., Fritzer, J., Steiner, A. K., Plach, A., Danzer, J., Foelsche, U., and Ladstädter, F.: Reference Occultation Processing System (rOPS) for cal/val and climate: a new GNSS RO retrieval chain with integrated uncertainty propagation, in: IROWG-4 International Workshop, Melbourne, Australia, 16–22 April 2015, 2015. a
Kirchengast, G., Schwärz, M., Schwarz, J., Scherllin-Pirscher, B., Pock, C., Innerkofler, J., Proschek, V., Steiner, A. K., Danzer, J., Ladstädter, F., and Foelsche, U.: Employing GNSS radio occultation for solving the global climate monitoring problem for the fundamental state of the atmosphere, in: EGU General Assembly 2016, Austria Center Vienna, Austria, 17–22 April 2016, 2016a. a
Kirchengast, G., Schwärz, M., Schwarz, J., Scherllin-Pirscher, B., Pock, C., Innerkofler, J., Proschek, V., Steiner, A. K., Danzer, J., Ladstädter, F., and Foelsche, U.: The reference occultation processing system approach to interpret GNSS radio occultation as SI-traceable planetary system refractometer, in: OPAC-IROWG 2016 International Workshop, Seggau Castle, Austria, 8–14 September 2016, 2016b. a
Kuo, Y.-H., Wee, T.-K., Sokolovskiy, S., Rocken, C., Schreiner, W., Hunt, D., and Anthes, R. A.: Inversion and error estimation of GPS radio occultation Data, J. Meteorol. Soc. Jpn., 82, 507–531, 2004. a
Kursinski, E. R., Hajj, G. A., Schofield, J. T., Linfield, R. P., and Hardy, K. R.: Observing Earth's Atmosphere with radio occultation measurements using the Global Positioning System, J. Geophys. Res., 102, 23429–23465, 1997. a
Li, Y., Kirchengast, G., Scherllin-Pirscher, B., Wu, S., Schwärz, M., Fritzer, J., Zhang, S., Carter, B. A., and Zhang, K.: A new dynamic approach for statistical optimization of GNSS radio occultation bending angles for optimal climate monitoring utility, J. Geophys. Res.-Atmos., 118, 13022–13040, https://doi.org/10.1002/2013JD020763, 2013. a
Li, Y., Kirchengast, G., Scherllin-Pirscher, B., Norman, R., Yuan, Y., Fritzer, J., Schwärz, M., and Zhang, K.: Dynamic statistical optimization of GNSS radio occultation bending angles: advanced algorithm and performance analysis, Atmos. Meas. Tech., 8, 3447–3465, https://doi.org/10.5194/amt-8-3447-2015, 2015. a
Li, Y., Kirchengast, G., Scherllin-Pirscher, B., Schwärz, M., Nielsen, J. K., and Yuan, Y. B.: A new algorithm for the retrieval of atmospheric profiles from GNSS radio occultation data in moist air conditions, in: OPAC-IROWG 2016 International Workshop, Seggau Castle, Austria, 8–14 September 2016, 2016. a
Li, Y., Kirchengast, G., Scherllin-Pirscher, B., Schwärz, M., Nielsen, J. K., Wee, T.-K., Ho, S.-P., and Yuan, Y. B.: A new algorithm for the retrieval of atmospheric profiles from GNSS radio occultation data in moist air conditions and cross-evaluation among processing centers, J. Geophys. Res.-Atmos., in submission, 2017. a, b
Palmer, P. I., Barnett, J. J., Eyre, J. R., and Healy, S. B.: A nonlinear optimal estimation inverse method for radio occultation measurements of temperature, humidity, and surface pressure, J. Geophys. Res., 105, 17513–17526, 2000. a
Rieder, M. J. and Kirchengast, G.: Error analysis and characterization of atmospheric profiles retrieved from GNSS occultation data, J. Geophys. Res., 106, 31755–31770, 2001. a
Scherllin-Pirscher, B., Kirchengast, G., Steiner, A. K., Kuo, Y.-H., and Foelsche, U.: Quantifying uncertainty in climatological fields from GPS radio occultation: an empirical-analytical error model, Atmos. Meas. Tech., 4, 2019–2034, https://doi.org/10.5194/amt-4-2019-2011, 2011a. a, b
Scherllin-Pirscher, B., Steiner, A. K., Kirchengast, G., Kuo, Y.-H., and Foelsche, U.: Empirical analysis and modeling of errors of atmospheric profiles from GPS radio occultation, Atmos. Meas. Tech., 4, 1875–1890, https://doi.org/10.5194/amt-4-1875-2011, 2011b. a, b, c
Scherllin-Pirscher, B., Steiner, A. K., Kirchengast, G., Schwärz, M., and Leroy, S. S.: The power of vertical geolocation of atmospheric profiles from GNSS radio occultation, J. Geophys. Res.-Atmos., 122, 1595–1616, 2017. a, b, c, d, e
Schreiner, W., Rocken, C., Sokolovskiy, S., Syndergaard, S., and Hunt, D.: Estimates of the precision of GPS radio occultations from the COSMIC/FORMOSAT-3 Mission, Geophys. Res. Lett., 34, L04808, https://doi.org/10.1029/2006GL027557, 2007. a
Schwarz, J., Kirchengast, G., and Schwärz, M.: Integrated uncertainty propagation for GNSS radio occultation: from excess phase profiles to dry-air variables, in: OPAC-IROWG 2016 International Workshop, Seggau Castle, Austria, 8–14 September 2016, 2016. a
Schwarz, J., Kirchengast, G., and Schwärz, M.: Integrating uncertainty propagation in GNSS radio occultation retrieval: From bending angle to dry-air atmospheric profiles, Earth Space Sci., 4, 200–228, 2017a. a, b, c, d
Schwarz, J., Kirchengast, G., and Schwaerz, M.: Integrating uncertainty propagation in GNSS radio occultation retrieval: from excess phase to atmospheric bending angle profiles, Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2017-159, in review, 2017b. a, b, c, d, e, f
Sokolovskiy, S., Rocken, C., Schreiner, W., and Hunt, D.: On the uncertainty of radio occultation inversions in the lower troposphere, J. Geophys. Res.-Atmos., 115, D22111, https://doi.org/10.1029/2010JD014058, 2010. a, b, c
Sokolovskiy, S., Zeng, Z., Schreiner, W., Hunt, D., Weiss, J., Lin, J., Kuo, Y.-H., and Anthes, R.: Specific errors of ionospheric correction in the troposphere induced by horizontal inhomogeneity of electron density in the ionosphere, in: OPAC-IROWG 2016 International Workshop, Seggau Castle, Austria, 8–14 September 2016, 2016. a
Sokolovskiy, S. V.: Modeling and inverting radio occultation signals in the moist troposphere, Radio Sci., 36, 441–458, 2001. a
Steiner, A. K. and Kirchengast, G.: Error analysis for GNSS radio occultation data based on ensembles of profiles from end-to-end simulations, J. Geophys. Res., 110, D15105, https://doi.org/10.1029/2004JD005687, 2005. a
Syndergaard, S.: Retrieval analysis and methodologies in atmospheric limb sounding using the GNSS radio occultation technique (PhD thesis), Sci. Rep. 99-6, Danish Meteorol. Inst., Copenhagen, Denmark, 1999. a
Xie, F., Syndergaard, S., Kursinski, E. R., and Herman, B.: An approach for retrieving marine boundary layer refractivity from GPS occultation data in the presence of superrefraction, J. Atmos. Ocean. Tech., 23, 1629–1644, 2006. a, b, c
Xie, F., Wu, D. L., Ao, C. O., Kursinski, E. R., Mannucci, A. J., and Syndergaard, S.: Super-refraction effects on GPS radio occultation refractivity in marine boundary layers, Geophys. Res. Lett., 37, L11805, https://doi.org/10.1029/2010GL043299, 2010. a, b, c
Zhu, Y., Derber, J., Collard, A., Dee, D., Treadon, R., Gayno, G., and Jung, J. A.: Enhanced radiance bias correction in the National Centers for Environmental Prediction's Gridpoint Statistical Interpolation data assimilation system, Q. J. Roy. Meteor. Soc., 140, 1479–1492, https://doi.org/10.1002/qj.2233, 2014. a