Dynamic statistical optimization of GNSS radio occultation bending angles : an advanced algorithm and its performance analysis

Introduction Conclusions References

Although the RO technique has been rather successful, it still suffers from some weaknesses.For example, the RO observations are affected by higher-order ionospheric effects and observation errors at high altitudes (> 30 km) (Bassiri and Hajj, 1993;Danzer et al., 2013;Liu et al., 2013Liu et al., , 2015)).These errors propagate downward from bending angles to refractivity through the Abel integral and also degrade the Published by Copernicus Publications on behalf of the European Geosciences Union.accuracy of the retrieved temperature and other atmospheric profiles (Healy, 2001;Rieder and Kirchengast, 2001;Gobiet and Kirchengast, 2004;Steiner and Kirchengast, 2005).Therefore, it is very important to have a best-possible initialization of the ionosphere-corrected bending angles at high altitudes for more accurate climate monitoring.
Statistical optimization is a commonly used method to initialize RO bending angles at high altitudes (e.g., Sokolovskiy and Hunt, 1996;Gorbunov et al., 1996;Hocke, 1997;Healy, 2001;Gorbunov, 2002;Gobiet and Kirchengast, 2004;Gobiet et al., 2007).It is a generalized least-squares approach that combines an observed RO bending angle profile with a background bending angle profile (Turchin and Nozik, 1969;Rodgers, 1976Rodgers, , 2000)).The weights of the two types of bending angles are determined by the inverse of their error covariance matrices.The statistical optimization equation used is (Healy, 2001;Gobiet and Kirchengast, 2004) where α s is the statistically optimized bending angle, α b and α o are the respective (unbiased) background and observed bending angle profiles, and C b and C o are the corresponding error covariance matrices.
In statistical optimization, the more accurately the error covariance matrices represent the error characteristics, the more accurate is the optimized bending angle profile.However, it is not straightforward to obtain such suitable error covariance matrices, especially for the background bending angle since they are not supplied together with common climatological models nor is the construction a straightforward task.Therefore, previous approaches have usually simplified the calculation of the error covariance matrices.
A typical approach is to estimate the background error covariance matrix by assuming a constant relative standard error of the background bending angle and a simple error correlation structure like exponential fall-off over an atmospheric scale height (Healy, 2001;Rieder and Kirchengast, 2001;Gobiet and Kirchengast, 2004) or disregarding correlations (Sokolovskiy and Hunt, 1996;Hocke, 1997;Gorbunov, 2002;Lohmann, 2005;Gorbunov et al., 1996Gorbunov et al., , 2005Gorbunov et al., , 2006)).Similarly, the observation error covariance matrix is formulated from estimating the observation error at a defined mesospheric altitude range (where the RO signal is weak) and using simple exponential fall-off error correlations (Healy, 2001, Gobiet andKirchengast, 2004) or again just ignoring the latter.These rough estimations generally result in inaccurate error covariance matrices and therefore result in inaccurate optimized bending angles that degrade the accuracy of subsequently retrieved atmospheric profiles.More details on the various schemes are provided by Li et al. (2013), Sect.2.1 therein.
Improved accuracy in optimized bending angles was obtained when using an improved statistical optimization algorithm to initialize ionosphere-corrected bending angles (Li, 2013;Li et al., 2013).Li et al. (2013) used European Centre for Medium-Range Weather Forecasts (ECMWF) shortrange (24 h) forecast fields as background bending angles.Their background error covariance matrix was accurately and realistically estimated using large ensembles of ECMWF short-range forecast, analysis, and RO observed bending angles.It was constructed using daily global fields of estimated background uncertainty profiles and a daily global-mean correlation matrix.The background uncertainty profile was dynamically estimated taking into account its variations with latitude, longitude, altitude, and day of year.They not only calculated the random errors of background bending angles using large ensembles of ECMWF and RO data but also empirically modeled the potential systematic background uncertainty and finally combined these two uncertainties to formulate the background uncertainty.The global-mean correlation matrix was also calculated using large ensembles of ECMWF analysis and forecast fields.Finally, the biases in the background bending angles were corrected to avoid the potential effects on optimized bending angles.
Since this first-version dynamic statistical optimization algorithm dynamically estimated the background error covariance matrix only, it is hereafter referred to as the b-dynamic algorithm ("b" represents background) in this study.The bdynamic algorithm was evaluated by Li et al. (2013) against the Occultation Processing System version 5.4 (OPSv5.4)algorithm developed by the Wegener Center for Climate and Global Change (WEGC) (Pirscher, 2010;Ho et al., 2012;Steiner et al., 2013).It was found that the b-dynamic algorithm significantly reduced random errors of the optimized bending angles and left less or about equal levels of residual systematic errors.The quality of the subsequently retrieved refractivity and temperature profiles was also improved.In addition, even the dynamically estimated background error correlation matrix alone was able to improve the optimized bending angles.
The aim of this study is to obtain even more accurate and reliable atmospheric profiles for optimal climate monitoring by advancing the b-dynamic algorithm to a complete dynamical estimation of both background and observation uncertainties and correlations.This is accomplished by employing a realistically estimated observation error covariance matrix in addition to the b-dynamic formulation.The observation error covariance matrix is constructed with dynamically estimated observation uncertainties for each occultation event and daily global-mean correlation matrices from large ensembles of data.The observation uncertainty is calculated using bias-corrected difference profiles of observed RO bending angle profiles relative to co-located ECMWF forecast bending angle profiles.The global-mean correlation matrix is calculated using multiple days of RO bending angle profiles and co-located ECMWF analysis bending angle profiles.In addition, the basic b-dynamic algorithm is updated to obtain even more robust background error covariance matrices.Fi-nally, the stability of the new dynamic algorithm is evaluated using full months of RO data from CHAMP and COSMIC.
The structure of this paper is as follows.Section 2 introduces the methodology of the dynamic statistical optimization algorithm.Section 3 evaluates its performance using both simulated and observed RO data.Finally a summary and conclusions are given in Sect. 4.

The new dynamic statistical optimization algorithm
Figure 1 shows the algorithmic steps of the dynamic statistical optimization algorithm.The new algorithm mainly includes two parts: (1) the dynamic estimation of the background error covariance matrix and the bias correction of background bending angles, and (2) the dynamic estimation of the observation error covariance matrix.Information on background/observation uncertainty and on the background/observation correlation matrix is needed for constructing complete background/observation error covariance matrices.The uncertainty at any vertical level is the square root of the diagonal value of the error covariance matrix at that level.The correlation matrix includes correlation functions for all vertical levels.Each such correlation function comprises the correlation coefficients of the error at the vertical level where it peaks to the errors at all other vertical levels.In summary, the background/observation uncertainties and the corresponding correlation matrix together formulate the background/observation error covariance matrices (Gaspari and Cohn, 1999).
Assuming that k = 1,. . ., N occ are sequentially numbered RO events of a day, then for each occultation event k the dynamic algorithm estimates the (unbiased) background bending angle profile α k b and its corresponding error covariance matrix C k b , as well as the observation error covariance matrix C k o .Using these three quantities together with the observed bending angle α k o , the statistically optimized bending angle profile α k s can be determined as The algorithm for the estimation of α k b and C k b has been described in detail by Li et al. (2013) as part of the b-dynamic algorithm.It will be briefly summarized in Sect.2.1, focusing on the key algorithmic steps and the advances in the dynamic algorithm.In Sect.2.2, details on the computation of C k o and other issues that are critical to the capability of the dynamic algorithm are provided.

Dynamic estimation of the background error covariance matrix and bias calibration of background bending angles
The dynamic estimation of the background error covariance matrix includes three algorithmic steps: (1) construction of basic daily background fields (blue boxes in the left part of Fig. 1), (2) preparation of the derived daily background fields (green boxes), and (3) dynamic estimation of the background error covariance matrix (orange boxes).
In the first step, daily fields of the basic background variables are prepared using a 10 • latitude × 20 • longitude grid (with the base cell centered at 5 • N, 10 • E), at 400 levels from 0.2 to 80.0 km with 200 m steps.This construction of background variables allows us to suitably capture large-scale background error dynamics as a function of latitude, longitude, (impact) altitude, and time, and it yields daily fields on a global 18 × 18 × 400 grid.Compared to the b-dynamic algorithm, which used 200 representative impact altitude levels from 0.1 km to 80.0 km with non-equidistant spacing, this new scheme allows direct use of these variables for the next step of calculation, avoiding additional vertical interpolation of all variables and particularly also within correlation matrices.
The data used to calculate these basic background variables include ECMWF analysis fields and corresponding 24 h forecast fields with a T42L91 resolution at 00:00 and 12:00 UTC, and observed RO bending angles.In calculating the mean variables in each grid cell, time averaging over 7 days (from 3 days before to 3 days after the day of interest) and horizontal averaging over geographic domains of at least 1000 km × 3000 km (over 10 • latitude × 60 • longitude cells from 60 • S to 60 • N, poleward over larger longitude ranges of 95 • from 60 to 70 • N/S, 120 • from 70 to 80 • N/S, and 270 • from 80 to 90 • N/S) were used.Compared to the b-dynamic scheme, which used 5 days of data only and smaller geographic regions (1000 km × 1000 km) for averaging, this update allows more data to be used for more reliable statistical estimates (especially important for mean observed bending angles) at each 10 • latitude × 20 • longitude grid point, while still capturing the slow variations of the mean field well.These calculations include mean variables, the construction of error correlation matrices, and empirical modeling of biases.For details see Li et al. (2013).
The second step involves the preparation of the derived daily background fields.These specific statistical quantities include (i) the forecast-minus-analysis standard deviations s f−a , which represent the estimated random uncertainty of the background bending angles; (ii) the estimated uncertainty of the mean background bending angle u b ; and (iii) the difference between the mean forecast bending angle and the mean background bending angle ᾱf−b .
In the third step, the background error covariance matrix is calculated from the fields obtained in step 2. Co-located profiles of ᾱk f−b , u k b , and s k f−a are derived by bi-linear horizontal interpolation to the RO event location at all vertical levels.The combined background standard uncertainty profile u occ b is then calculated as Herein the bias coverage factor f bc is introduced as a userdefined parameter that can be employed to penalize the estimated bias-type uncertainty u k b relative to the estimated random uncertainty s k f−a .This enables minimizing the influence from potential residual background biases, relative to observation uncertainty, on the resulting optimized profile α k s .In this study the bias coverage factor was chosen to linearly decrease with altitude, setting f bc = 15 at 30 km (strong penalty in lower stratosphere) and f bc = 1 at 80 km (no penalty at top boundary).This choice was found to be useful for climate applications (more discussion of f bc , including for comparison also example cases with constant f bc = 1 and f bc = 5, is given in Sect.3.2).Using the background standard uncertainty and the global-mean correlation matrix R f−a = R f−a,ij , the background error covariance matrix In order to effectively reduce the residual bias in the background bending angle α k b (cf.Eq. 2) to within the estimated uncertainty of the background mean u k b , α k b is calculated by subtracting the forecast-minus-background mean difference ᾱk f−b from the co-located forecast bending angle α k f : Figure 2 illustrates the estimated relative uncertainty of the forecast-minus-analysis standard deviation 100•(s f−a / ᾱa ) (top) and the modeled relative systematic bias of the ECMWF analysis bending angle 100 • (b a / ᾱa ) (bottom).b a is the main component of u k b (cf.Li et al., 2013) and ᾱa is the ECMWF analysis mean bending angle; b a therefore illustrates the size of the systematic uncertainty term f bc • u k b in Eq. (3) when f bc = 1, i.e., when no user-defined bias penalty is employed.Any application of a bias coverage factor f bc >1, such as used in this study for representing climatequality retrievals, magnifies this term accordingly.
As a function of latitude (left), forecast-minus-analysis standard deviation reveals largest errors at high altitudes over the Antarctic region.The relative standard deviations are larger than 10 % near 80 km, decreasing to ∼ 5 % at 70 km and to ∼ 3 % at 50 km, and remain at ∼ 1 % below 25 km.In non-polar regions, the standard deviations amount approximately to 3 to 4 % near 80 km, decrease to 2 % at 75 km and remain within 1 to 2 % below.The day-to-day temporal evolution of s f−a at high southern latitudes (60 to 70 • S) (right) reveals relatively smooth uncertainty conditions for the entire month of July, without much temporal dynamics in the uncertainty.In the course of testing, it was found that the estimates are consistent with flow-dependent forecast-minusanalysis error estimates produced by ECMWF's Ensemble of Data Assimilations (EDA) system (Isaksen et al., 2010;Bonavita et al., 2011;M. Bonavita, ECMWF, personal communication, 2012).
Regarding variations of 100 • (b a / ᾱa ) as a function of latitude (bottom left), the bias-type uncertainties are also largest at high altitudes in the Southern Hemisphere.The relative uncertainties are larger than 7 % near 80 km, decreasing to 5 % at 60 km, and remain smaller than 3 % below 40 km.In nonpolar regions, the relative uncertainties amount to 4 % near 80 km, decreasing to below 2 % also at 40 km.The temporal evolution of the systematic uncertainty 100 • (b a / ᾱa ) over a month (bottom right) shows that these relative uncertainties also reveal little sub-monthly variations, due to the way in which they are constructed (Li et al., 2013).
Figure 3 shows exemplary global-mean correlation functions (left) and associated correlation lengths (right) for the 5, 15, and 25 July 2008.The correlation functions, plotted for three representative height levels (30, 50, 70 km), are rather similar over the month.They have a main peak of nearly Gaussian shape with negative side peaks at each side.Further outward, small secondary positive peaks occur, after which the functions essentially approach zero.The correlation lengths increase rather smoothly with altitude from about 0.8 km at 20 km to 6 km at 80 km.They also show little variation over the example month of July 2008.
Overall this behavior indicates that in months without larger atmospheric anomalies (such as for example sudden stratospheric warming at high latitudes; e.g., Manney et al., 2008) a daily update of correlation matrices is not necessarily needed.In a long-term application, however, it is not clear when and where some (transient) anomalies may occur, so a

Y. Li et al.: Dynamic statistical optimization of GNSS radio occultation bending angles
daily update of the background fields was used as a cautious baseline.

Dynamic estimation of the observation error covariance matrix
The error covariance matrix of the observed bending angle C k o is calculated using an observation uncertainty profile u occ o , estimated on a per-event basis, and a global-mean error correlation matrix R occ o : Different to the OPSv5.6 and the b-dynamic algorithm, which estimate the observation uncertainty u occ o between about 65 and 80 km with an MSIS climatology model bending angle profile as a reference (Li et al., 2013), the full dynamic algorithm estimates u occ o as a vertical profile over the stratopause region and mesosphere, using the co-located ECMWF forecast bending angle profile as a reference.
More specifically, the first step is to subtract the co-located forecast bending angle profile α k f from the observed bending angle profile α k o : The difference profile α k o is then smoothed with a 15 km window moving average (from 45 km to the top bound of the profile, usually 80 km).The resulting smoothed difference profile is denoted as α k o .The next step is to subtract the smoothed difference profile α k o from the original difference profile α k o in order to obtain a delta-difference profile α k o that essentially contains only random errors: Finally the observation uncertainty at any impact altitude level i (corresponding to the impact altitude of z i ) of each occultation event, u occ o,i , is calculated as where n is the number of sample points between z i + 7.5 km and z i − 7.5 km and where i z−7.5 and i z+7.5 denote the corresponding impact altitude indices.Equation ( 9) is applied from 45 km to the top of the profile, providing u occ o,i estimates from 52.5 to 72.5 km; below (above) the value at 52.5 km (72.5 km) is extended downward (upward) just as a constant value.This construction ensures that variations over the mesosphere can be accounted for while data below the stratopause, where the estimated delta-difference profile α k o may also contain atmospheric variability noise (e.g., from gravity wave activity), are not allowed to influence the estimate.7) and (8), for six example RO events from 15 July 2008.The simMetOp events are simulated using the End-to-end GNSS Occultation Performance Simulation and Processing System software version v5.6 (EGOPSv5.6) in the same way as the sim-MetOp ensemble was produced by Li et al. (2013); for more details see Sect. 3 below.In each panel, the red line shows the original difference profile α k o ("delta"), the blue line the smoothed difference profile α k o ("deltasmooth"), the green line the delta-difference profile α k o ("deltadelta"), and the magenta line the observation uncertainty u occ o ("Uncert").It can be seen that Eqs. ( 7) and ( 8) are robust in removing systematic errors, leaving a good random signal, from which the uncertainty u occ o is reliably estimated.Estimated uncertainties are smallest (near 0.5 µrad) for the simMetOp events (top) that mimic MetOp/GRAS-type high performance receiver errors without additional noise effects (Li et al., 2013) and are distinctively larger for real data from CHAMP (middle) and COSMIC (bottom).The uncertainties of CHAMP and COSMIC bending angles can be quite variable, as illustrated, and can reach 10 µrad or more for CHAMP events, while it is typically only around 2 µrad or so for COSMIC.Figure 4 also shows that the variation of uncertainty, which often may come from variations in ionospheric small-scale noise being added onto the RO receiver-related noise, can be well captured over the mesosphere.
The global-mean observation error correlation matrix R occ o is estimated using the difference profiles between the observed and co-located ECMWF analysis bending angle profiles.In this estimation, we first construct a global-mean error covariance matrix using all available difference profiles from 3 days before to 3 days after the day of interest, i.e., using the same weekly smoother as for the background estimations.
Then we derive the global-mean correlation matrix by dividing all elements of the covariance matrix by their corresponding square roots of diagonal values.RO observations used for this calculation included data from COSMIC and GRACE for January and July 2008, and MetOp-A for July 2008 (no processed data were available for January 2008).The obtained error correlation matrix R occ o is used based on Eq. ( 6) for the construction of the observation error covariance matrices C k o , for all RO events k of the day.We did not include a bias calibration step similar to the background bias calibration (Eq.5) at the observation side.Such a step may complement the C k o estimation in the future, however, for subtracting residual ionospheric biases in the observed bending angle profiles α k o before they enter the calculation of mean observed bending angles (as part of estimating ᾱk f−b ; Li et al., 2013) and the statistical optimization equation (Eq.2).
Figure 5 illustrates representative observation error correlation functions (extracted from R occ o ) and associated correlation lengths in the identical format as shown in Fig. 3 for the background error correlations, for ease of intercomparison "delta" is the difference profile of the RO ionosphere-corrected bending angle to the co-located ECMWF forecast profile used as a reference, "deltadelta" is the delta-difference profile after subtracting a smoothed profile "deltasmooth" from the difference profile "delta", and "Uncert" is the resulting observation uncertainty estimate; for detailed description see Sect.2.2.For ease of intercomparison the layout is the same as in Fig. 3.
of the different characteristics.By comparing Figs. 3 and 5, it can be seen that the observation error correlation functions are basically similar in shape to the background error correlation functions (main peak, negative side peaks, smaller positive secondary side peaks) albeit with significantly shorter correlation lengths.In addition, the functional shape of side peaks is not as smooth as the one for the background.The correlation length is essentially constant with altitude (above 30 km), amounting to about 0.8 km.
The intra-monthly variation is essentially negligible within the given July 2008 test month (applies also to January 2008, not shown), pointing to room for further improvement of the utility of the estimation for long-term processing, e.g., considering larger ensemble sizes and sub-global regions.These slow dynamics of the observation error covariance matrices, and of the background error covariance matrices as discussed in Sect.2.1, enable reliable use also in near-real-time or fasttrack processing (i.e., processing within 3 h or within followon day of observations).Instead of using 7 days centered about the day being processed (including 3 days before and after the center day) 7-day-history data (from the previous day to 7 days prior) may be used in these cases, with insignificant degradation in performance.

Other improvements of the new algorithm
In the b-dynamic algorithm of Li et al. (2013), the statistical optimization was applied exactly down to 30 km.However, for some noisy RO events especially from CHAMP, ionosphere-corrected bending angles can be still noisy around 30 km impact altitude.For these noisy events there can be a sharp change of bending angle characteristics from the rather smooth statistically optimized profile above 30 km to the rather noisy purely observed profile below 30 km.On the other hand, some (simulated) events may have very small observation uncertainties at the altitudes below 40 km, which may lead to degraded robustness of the matrix inversion (C b + C o ) −1 , due to large differences of observation and background uncertainties.In order to safeguard against these effects, we have improved the algorithm as follows.
First, we gracefully adjust the observation uncertainty u o computed according to Sect.2.2 based on adjusting the observation-to-background uncertainty ratio r obu = u o /u b below 40 km, in order to ensure that it approaches a robust small value at the bottom of the statistical optimization range.We modify the r obu profile to make it linearly transit from the r obu value prevailing at 40 km to r obu = 0.1 at 28 km.Typical values within 28 and 40 km may range from near 0.1 (simMetOp) to around 1 (CHAMP), so the linear transit to 0.1 at 28 km generally implies a decrease of the r obu for real data.In order to avoid any possible sharp change of the r obu profile at 40 km, we smooth the resulting profile by a moving average filter with 2 km width.Using the modified r obu profile, the observation uncertainty is then reconstructed as , enforcing dominance of the observation information in the optimized profile from 40 km to 28 km.Alternatively to this u o (z i ) adjustment used in this study, the background uncertainty profile may be adjusted, u b (z i ) = u o (z i )/r obu (z i ), which keeps the effect on the optimized profile the same while leaving the observation uncertainty unchanged.
Second, we apply the statistical optimization down to 28 km and then apply a half-sine-weighted transition across 32 to 28 km between the statistically optimized bending angles and purely observed bending angles.That is, the weighting function over this transition altitude range, w(z i ), is formulated as where z soT is the statistically optimized-to-observed bending angle transition altitude, set to 30 km, and z soT is the statistically optimized-to-observed bending angle transition half width, set to 2 km.Employing w(z i ) from 32 to 28 km in the form we get a well-defined gradual transition from optimized bending angle to observed bending angle.When these two improvements to the b-dynamic algorithm are applied, the sometimes "spiky" behavior of some (CHAMP) profiles near the bottom boundary of statistical optimization discussed by Li et al. (2013) disappeared.
Another issue requiring caution is the robustness of matrix inversions, especially related to the weighting matrix in Eq. ( 2), C b (C b + C o ) −1 , containing the inversion of the summary matrix of the observation and background error covariance matrices.Since the observation error correlation functions can be insufficiently smooth and since the main peaks are close to Gaussian shape, the summary matrix can become close to ill-conditioned and cannot be accurately inverted in this case, inducing undue noise into the resulting weighting matrix.
This could be overcome by replacing the main peak by a 5th-order polynomial function as described by Gaspari and Cohn (1999), which approximates a Gaussian shape and was also successfully used in the context of matrix inversion by Steiner and Kirchengast (2005), Eq. ( 5) therein (note that a typo leaked into the z 2 term as cited in this equation for z < 1; the correct denominator is "3" instead of "2", as in the z 2 term cited for the 1 < z branch).After using this approximation, the inversion was sufficiently robust.Alternatives for ensuring robustness include the use of even-larger ensemble sizes in matrix construction and more advanced methods of matrix inversion such as truncated singular value decomposition.

Evaluation of the new dynamic statistical optimization algorithm
The dynamic algorithm was implemented in the EGOPSv5.6 software (Fritzer et al., 2013) to enable a complete RO retrieval.The EGOPSv5.6 system was also used to simulate RO observations (simMetOp events) and to retrieve atmospheric profiles.The standard RO data processing chain within this system is the OPSv5.6 retrieval.
We evaluated the new dynamic algorithm against this OPSv5.6algorithm, which is, in terms of statistical optimization formulation, the same as the OPSv5.4algorithm used for comparison by Li et al. (2013).Briefly, the OPSv5.6 algorithm uses ECMWF short-range forecast bending angles as a background and employs exponential fall-off functions to express the correlations of both background and observation uncertainties.The background uncertainty is modeled as amounting to 15 % of background bending angles.The observation uncertainty is estimated as the standard deviation of observed bending angles relative to co-located MSIS model bending angles in the impact altitude range from 65 km to about 80 km.For more detailed information on OPSv5.6/v5.4see Pirscher (2010), Steiner et al. (2013), andSchwaerz et al. (2013).
In addition to the OPSv5.6 intercomparison, the atmospheric profiles retrieved by the dynamic algorithm are compared with those retrieved by the b-dynamic algorithm (Li et al., 2013) and with those by the UCAR/COSMIC Data Analysis and Archive Center (CDAAC) Boulder.
The data sets used for the evaluation include simulated MetOp data (simMetOp) as well as real observed CHAMP and COSMIC data.simMetOp data were simulated in the same way as by Li et al. (2013), using moderate ionosphere conditions in the forward simulations and using observational errors representing MetOp/GRAS-type receiving system errors.
As a basis for the CHAMP and COSMIC retrievals, excess phase and orbit data were downloaded from UCAR/CDAAC Boulder (CDAAC data version 2009.2650 for CHAMP and 2010.2640 for COSMIC).CDAAC atmospheric profiles (atmPrf) used for the comparison of retrieved profiles are mainly from the same CDAAC data version.Recently reprocessed atmospheric profiles provided by CDAAC (version 2014.0140) are also used for comparison, but due to their very recent release only CHAMP has been used so far.We denote the CHAMP data version 2009.2650 as "CDAAC" in figure legends, and version 2014.0140 as "CDAAC new ".In the evaluation, retrieved RO profiles are shown relative to co-located reference profiles.For the CHAMP and COS-MIC data, these co-located reference profiles were extracted from ECMWF analysis fields; for simMetOp data the "true" ECMWF analysis field profiles from the forward modeling were used as a reference.This is the same setup as was used by Li et al. (2013).

Algorithm performance for individual profiles
Figure 6 illustrates the effects of statistical optimization on individual bending angle profiles by a few representative RO events.The left panels show the background and observation uncertainties as well as the observation-to-background (Obsto-Bgr) weighting ratio , which expresses on a percentage scale how the information is weighted between observations and background.r obw is a convenient approximate variable for this purpose, which exactly applies if the covariance matrices in the statistical optimization equation (Eq.2) are diagonal.Observation uncertainty is smallest for the simMetOp event (top), largest for the CHAMP event (middle), and in between for the COSMIC event (bottom).At 60 km these observation uncertainties are roughly 0.4, 3, and 1.5 µrad for simMetOp, CHAMP, and COSMIC, respectively.These differences in observation uncertainty yield the largest r obw for the sim-MetOp event and the smallest for CHAMP.Related to this, the altitude where both observations and background receive equal weight (r obw = 50 %), which may be considered the transition altitude below which the observation information dominates the retrieval, is highest for simMetOp (> 60 km), medium for COSMIC (around 55 km), and smallest for CHAMP (near 45 km).Ongoing follow-on work on large RO data sets of many months analyzes the r obw statistically and confirms that these few example events are typical for the respective data sources.
Since bending angles increase roughly exponentially with decreasing altitude, as seen in the middle column of Fig. 6, differences among the various retrieved profiles seem to be small.The right panels, however, actually show the differences of optimized bending angle profiles relative to their reference.For the simMetOp event, bending angle differences from the dynamic algorithm are smallest over all altitudes, confirming the high utility of the algorithm, since here the "true" profile from forward simulation serves as a reference.The bending angle differences of the b-dynamic algorithm are similar and the values are also small.The dif-ferences from the OPSv5.6 algorithm are largest and significantly noisier.
For the CHAMP event the relative differences from the dynamic algorithm, the b-dynamic algorithm, and CDAAC are smaller than those from the OPSv5.6 algorithm below 50 km.Above about 50 km, the relative differences from the dynamic algorithm increase and are largest.For the COSMIC event, bending angles from the dynamic, b-dynamic, and OPSv5.6 algorithms are rather similar below 50 km.Above 50 km, differences from the dynamic algorithm are tentatively largest.For this event, the differences of CDAAC are generally larger than of the other three algorithms.
Inspecting further individual RO events (not shown) confirmed that the relative differences of simMetOp data from the dynamic algorithm are consistently smaller and smoother than those from the other approaches.This underlines the robust capability of the dynamic algorithm for improving the quality of the ionosphere-corrected bending angles.For CHAMP and COSMIC measurements, the relative differences from the dynamic algorithm are also generally smaller and smoother than those from other algorithms below 50 km.However, above 50 km, the differences from both the dynamic algorithm and from CDAAC are generally larger than those from the OPSv5.6 and b-dynamic approaches.This does not mean that bending angle profiles from the dynamic and CDAAC algorithms are not accurate at high altitudes, however; the result mainly depends on the determination of the weights of the background and observed bending angles in the statistical optimization.
In the new dynamic algorithm, the estimated relative background errors at high altitudes are generally larger than those of the OPSv5.6 algorithm (e.g., in the polar winter regions, the relative background errors of the dynamic algorithm are usually larger than 30 % around 60 km, while the OPSv5.6 algorithm assumes a constant error of 15 % at all altitudes).At the same time the dynamically estimated observation uncertainties are usually smaller than those of OPSv5.6, which often sets a large standard value (22 µrad) for events that are very noisy at high altitudes.Therefore the new dynamic algorithm gives significantly more weight to the observed bending angles at high altitudes.If a user of the dynamic algorithm wanted less observational weighting, this would need to be tuned by the bias coverage factor f bc (cf.Eq. 3), which would accordingly modify the r obw .For example, using a constant setting of f bc = 5, as used for the algorithmic introduction by Li et al. (2013), would significantly change the r obw compared to the linear f bc altitude dependence used in this study (for more discussion on the effects of f bc choice see Sect.3.2 below).Furthermore, similar to what Li et al. (2013) did for the b-dynamic algorithm, we also inspected the effects of using two different correlation matrices, i.e., the estimated global-mean correlation matrices of our dynamic algorithm and the simple analytical correlation matrices constructed by exponential fall-off correlation functions with a correlation length of 10 km for background and 2 km for ob-  servation as used in the OPSv5.6 formulation.The dynamically estimated uncertainties were used for both cases since we are interested only in the differences from the correlations in this particular test.
Figure 7 shows the comparison results for these two types of correlation matrices, again using three exemplary events from simMetOp, CHAMP, and COSMIC, showing absolute bending angles (left) and differences to reference (right).The simMetOp event highlights that bending angle differences from the full correlation case are much smoother and smaller than those from the exponential fall-off correlation.For the real CHAMP and COSMIC events, the magnitudes of the differences from the two cases are similar, but also here it is clearly evident that the use of the full correlation leads to smoother differences than the use of exponential fall-off correlation.We conclude that the use of adequately realistic correlation matrices is preferable.

Statistical performance evaluation results
In this section the performance of the dynamic, b-dynamic, and OPSv5.6 statistical optimization algorithms are evaluated using simulated MetOp data from 15 January and 15 July 2008 and monthly CHAMP and COSMIC observations from January and July 2008.In addition, atmospheric profiles retrieved and provided by UCAR/CDAAC for the same time periods are used for comparison.The mean systematic differences between retrieved and reference profiles and the associated standard deviations are calculated and analyzed in bending angle, refractivity, and temperature profiles, similar to the statistical performance evaluation of the b-dynamic algorithm by Li et al. (2013).
In order to detect and exclude outlier profiles from the statistical profile ensembles, the quality of retrieved profiles is checked as follows.Bending angle profiles are checked from 25 to 80 km, and a profile is flagged as bad if a bending angle exceeds a threshold at any impact altitude level, which was defined based on careful sensitivity tests as the maximum of either 40 µrad absolute or 25 % relative deviation from the co-located ECMWF short-range forecast bending angle.In practice, looking at it from the top downwards, the transition from the absolute to the relative criterion occurs roughly around 35 km, where bending angle values start to exceed 160 µrad.
Refractivity and temperature profiles are checked in the same way as used for OPSv5.6/v5.4 (Schwaerz et al., 2013;Steiner et al., 2013;Pirscher, 2010), i.e., the deviation from co-located ECMWF analysis profiles at any altitude level must not exceed 10 % for refractivity between 5 and 35 km and 20 K for temperature within 8 and 25 km.These quality checks are performed on all profiles retrieved with the EGOPS software.For profiles provided by CDAAC, we check the CDAAC quality flag and use only profiles flagged to be of good quality.
Figure 8 shows the systematic differences and standard deviations of optimized bending angle profiles of the global ensembles of simMetOp from 15 January and 15 July 2008, and of CHAMP and COSMIC events from January and July 2008.For simMetOp (top), it is clear to see that the performance of the dynamic algorithm outperforms the b-dynamic algorithm and the OPSv5.6 algorithm, exhibiting smallest systematic differences and associated standard deviations.Compared to the OPSv5.6 algorithm, the best improvement is found between 40 and 60 km.These results are very en-couraging and confirm the fundamental capabilities of the dynamic algorithm.
Comparison of the CHAMP (middle) and COSMIC (bottom) results from the dynamic algorithm with the OPSv5.6 and CDAAC algorithms shows that the bending angle standard deviation from the dynamic (and b-dynamic) algorithm is again generally smaller than that of the OPSv5.6 and CDAAC results.Above 45 km for CHAMP, and above 55 km for COSMIC, the standard deviations from the dynamic algorithm exceed those from the OPSv5.6 algorithm.This is due to increased weight of noisy RO bending angles in the mesosphere compared to OPSv5.6 as discussed in Sect.3.1 above.Standard deviations from both CDAAC data versions are larger than those from the other methods, and particularly the new data version (shown for CHAMP) exhibits largest standard deviation already from about 35 km upwards.
Systematic differences of CHAMP and COSMIC data (differences are calculated against co-located ECMWF analysis profiles) are rather similar for the dynamic, b-dynamic, and OPSv5.6 algorithms.The systematic differences from CDAAC algorithms are also similar, in particular below about 35 km, but they are larger and feature somewhat different characteristics above about 40 km for July 2008.In particular the new data version (shown for CHAMP) exhibits different (oscillatory) behavior both in January and July.These results indicated that the new dynamic algorithm is robust and competitive in providing optimized profiles with biases minimized in a best-possible manner, as should be expected from its realistic account for both observation and background uncertainties and error correlation structures.
Furthermore it can be seen, in particular from the CHAMP results (CHAMP data have highest observational noise), that the improved treatment of the transition to purely observed data around 30 km has mitigated the sharpness of the change in standard deviation.
In order to discuss the effects of different choices of f bc on the resulting optimized bending angles, we compared the choice of this study for linear altitude dependence (see Sect. 2.1; termed f bc = 1 To 15 here) with the choice in the algorithmic introduction by Li et al. (2013) (f bc = 5) and with a reference case intentionally making no use of the bias penalty option (f bc = 1).Figure 9 shows the comparative re-Figure 9. Systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines) of statistically optimized bending angles, relative to "perfect" simulated bending angles or co-located ECMWF analysis bending angles used as a reference, of the global ensemble of simMetOp events on 15 January and 15 July 2008 (upper two panels), and of CHAMP and COSMIC events from the complete months of January and July 2008 (middle and bottom panels, respectively).Results from three different bias coverage factor choices in the dynamic algorithm -i.e., f bc = 1 (green), f bc = 5 (blue), and f bc = 1 To 15 (red) -as well as from the OPSv5.6 algorithm (black) are shown.
sults for these three f bc choices in the same statistical result format as used for Fig. 8; the results shown for "Dynamic f bc = 1 To 15" and OPSv5.6 replicate the ones of Dynamic and OPSv5.6 of Fig. 8 for context.It can be seen that for sim-MetOp the f bc = 5 and f bc = 1 To 15 results are rather similar, both for systematic differences and standard deviations, while the f bc = 1 choice indicates how a strong un-penalized weight of background profiles forces the optimized solution towards the background at high altitudes (here visible above about 45 km).As is to be expected, any residual biases in the background will therefore best survive in the optimized solution in the case of an f bc = 1 choice.This is the very reason why the background is intentionally penalized if the dynamic scheme is employed for climate-quality retrievals that strive to minimize background influence.For the CHAMP and COSMIC data, standard deviations are largest for the f bc = 1 To 15 case, medium for the f bc = 5 case, and smallest for the f bc = 1 case.This is in line with the expectation of how the r obw changes under these different choices since the observational noise is more and more mitigated the more relative weight the background receives.The effect on the systematic differences (for CHAMP and COS-MIC relative to the co-located ECMWF analysis profiles) is relatively small in these global-scale statistics; but also here the direction is that no bias penalty forces the results towards the background mean state at high altitudes.
Overall Fig. 9 demonstrates that the sensitivity to the detailed quantitative choice of f bc is fairly weak, as can be seen from the moderate differences between these three cases with very different f bc choices.This is favorable since it implies that no detailed quantitative tuning of f bc is needed in practice; the f bc as the only free user-defined variable in the new dynamic algorithm is rather a clear and transparent option to predefine the influence of background information according to what users deem suitable for their application.
In order to evaluate the performance of different statistical optimization algorithms in different latitude regions, the systematic differences and standard deviations of optimized bending angles were calculated for five latitudinal bands in addition to the global case (90 • S to 90 • N), including tropics (TRO, 20 • S to 20 • N), Southern Hemisphere/Northern Hemisphere subtropics and midlatitudes (SHSM/NHSM, 20 • S/N to 60 • S/N), and Southern Hemisphere/Northern Hemisphere polar regions (SHP/NHP, 60 • S/N to 90 • S/N).Figures 10,11,and 12 show the statistical results for the global case (for context, same as in right column of Fig. 8) and for these five regions for simMetOp (Fig. 10), CHAMP (Fig. 11), and COSMIC (Fig. 12).The July 2008 results are shown, which are found to be well representative; the latitude-resolved data characteristics in January 2008 are similar.
Figure 10 shows that the performance of the dynamic, bdynamic, and OPSv5.6 algorithms are rather similar globally and in the Northern Hemisphere (NHSM, NHP).In the Southern Hemisphere, and in particular in the SHP region (Antarctic winter in July), the conditions are evidently more challenging, such that the OPSv5.6 algorithm accrues increased biases in the upper stratosphere above 50 km.The new dynamic algorithm underscores its good and reliable basic performance in all regions, both in terms of biases and standard deviations.
Figures 11 and 12 for the real bending angle data show that the performance of all algorithms in all latitude bands except SHP (Antarctic winter) is consistent and generally similar to the performance visible from the global ensemble, which we discussed along with Fig. 8 above.In SHP, both systematic differences and standard deviations are markedly larger.In particular the CDAAC results, and most so the new CDAAC data version, exhibit relatively large systematic differences  at altitudes above 35 km (up to around 5 %) and also standard deviations closely reaching or exceeding 10 % even at altitudes near 50 km.Error characteristics for January (not shown) in the NHP region (Arctic winter) generally mirror the SHP July (Antarctic winter) error characteristics.
We consider the new dynamic algorithm in this context to confirm its robust performance also for real data in all regions, although the lack of a "true" reference in these cases does not allow for strong conclusions.In evaluating future long-term processing application of the algorithm, we will also include stricter validation against independent colocated data of high quality over the stratosphere and mesosphere from other sources such as the Envisat/MIPAS and TIMED/SABER satellite instruments (Remsberg et al., 2008;García-Comas et al., 2012).
Figures 13 and 14 show the global statistics results for refractivity (Fig. 13) and temperature (Fig. 14) for simMetOp (top), CHAMP (middle), and COSMIC (bottom).These refractivity and temperature results reflect the results for the bending angles in a filtered manner, after having passed through the Abelian integration (refractivity) and in addition the hydrostatic integration (temperature), which lead to smoothing and downward propagation of biases and to reduction of standard deviations (e.g., Gobiet and Kirchengast, 2004;Steiner and Kirchengast, 2005).Due to this downward propagation, the differences from the various algorithms become smaller, and results are closely similar below about 40 km and in most cases even above.The most notable differences from the consistent behavior of the different algorithms are those of the CDAAC processings above about 50 km, which exhibit the largest systematic differences and standard deviations, and the additional deviations of the new CDAAC version (shown for CHAMP) already from about 35 km upwards.
Again we consider the performance of the new dynamic algorithm robust and encouraging for larger-scale applications, which may also include further adjustments of parameters like f bc and of averaging domains for constructing the dynamic uncertainty and correlation information.
Figure 15 extends the view of Fig. 7 on the sensitivity to the choice of correlation modeling to a statistical view.It depicts the results of global statistics for simMetOp (top) and COSMIC (bottom), for bending angle (left), refractivity (middle), and temperature (right), from either operating the Bending angle (left), refractivity (middle), and temperature (right) systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines), relative to their "perfect simulated" or co-located ECMWF analysis data used as a reference, of the global ensemble of simMetOp (top) and COSMIC (bottom) events from 15 July 2008, using either the realistic global-mean correlation matrix of the new dynamic method ("full correlation") or simple exponential fall-off correlation as in the existing OPSv5.6 ("exp.falloffonly").The number of events (NoE) in each statistical ensemble is also indicated in the panels.
full dynamic algorithm or from using simplified correlation modeling with the exponential fall-off approximation.The simMetOp results show, in line with the results of Fig. 7, that the use of the realistically modeled full correlations is a superior choice, though the reduction of systematic difference relative to the "true" reference is small (after the Abel or hydrostatic integrations) in refractivity and temperature.The COSMIC results indicate that the choice of correlation modeling strongly impacts the standard deviation and to a more limited degree also the systematic differences.While this behavior does not itself imply a preference, it is very clear that the choice of the realistic full correlation modeling will be the physically more sound and more adequate approach also for real data.

Summary and conclusions
This study presented a new dynamic statistical optimization algorithm to initialize RO ionosphere-corrected bending angle profiles at high altitudes for optimal climate monitoring throughout the stratosphere.This dynamic algorithm uses multiple days of ECMWF analysis, ECMWF short-range (24 h) forecast, and RO observation data to realistically estimate background and observation error covariance matrices.Both the background and observation error covariance matrices are constructed with geographically varying uncertainty estimation and with a global-mean correlation matrix estimated on a daily basis.The b-dynamic algorithm recently introduced by Li et al. (2013) was used as a starting point and provided for the estimation of the background error covariance matrix and the bias correction of background bending angles.
The main advancements of the new dynamic algorithm compared to this previous algorithm are that it (1) adds a dynamically estimated observation error covariance matrix with altitude-dependent observation uncertainty and a realistically calculated global-mean correlation matrix; (2) updates the algorithm of the calculation of basic statistical mean variables by using ECMWF and RO data from a longer time window and larger geographical regions for more accurate and reliable estimation; and (3) eliminates weaknesses that existed near the lower boundary of statistical optimization (30 km) by improving the uncertainty formulation and transition to purely observed data across this boundary.
We illustrated and discussed key variables of the dynamic background and observation error covariance matrices, including systematic and random uncertainties and correlation functions, in order to provide insight and show the realistic character of their behavior.Both the random and systematic background uncertainties appear to be largest in the polar regions of the winter hemisphere at mesospheric altitudes.The observation uncertainties capture variations with altitude, especially in the mesosphere, and can well represent the error characteristics of RO events, from high-quality simulated data to comparatively noisy CHAMP data.The observation error correlation functions show a similar shape to background data, but with less functional smoothness and with much smaller correlation lengths of about 0.8 km (while background error correlation lengths range from about 1 km near 20 km to about 6 km near 80 km).All uncertainty and correlation estimates were found to exhibit little sub-monthly variations during the test months January and July 2008.In the case of anomalous sub-monthly conditions (e.g., sudden stratospheric warming) that will occasionally happen during long-term processing periods, we expect more variation, however.
The new dynamic algorithm was evaluated mainly against the one currently used in the OPSv5.6 system, using simulated MetOp data on single days (15 January and 15 July 2008) and real observed CHAMP and COSMIC data from two full months (January and July 2008).The following was found for the new dynamic algorithm, in particular compared to OPSv5.6: (1) it can reduce systematic errors (biases) and standard deviations of optimized bending angles, as proven by simMetOp data including "true" reference profiles from end-to-end simulations, and subsequently also benefits the error characteristics of retrieved refractivity and temperature profiles; (2) it can reduce the random errors of optimized bending angles in the stratosphere for real data, as evaluated for CHAMP and COSMIC, still at the same time leaving less or about equal residual systematic error (bias) in the bending angles; (3) it can better account for the observational noise in the mesosphere, leading to larger standard deviations than OPSv5.6 there from greater weight of the observations in the optimized profiles, albeit without applying any artificial observation uncertainty values in the case of high noise levels.
Beyond the evaluation of the new dynamic algorithm against OPSv5.6,atmospheric profiles from UCAR/CDAAC were also intercompared, including use of very recently released CHAMP data from the newest (2014) CDAAC data version.It was found that CDAAC bending angles generally exhibit markedly higher standard deviations above about 35 km and that in particular the new data version shows comparatively large systematic differences and standard deviations.The reasons for this new-version behavior deserve further study.
Overall, compared to previous simplified approaches of statistical optimization, the dynamic algorithm presented here, which realistically estimates both background and observation error covariance matrices, contains high capabilities for future large-scale implementation.The evaluation of the algorithm provided clear evidence that it can deliver reliable and accurate atmospheric profiles for atmosphere and climate applications.The results therefore indicate high suitability for employing the new dynamic approach in the processing of long-term RO data into a climate record, leading to well-characterized and high-quality atmospheric profiles over the entire stratosphere.

Figure 1 .
Figure 1.Schematic illustration of the algorithmic steps of the dynamic statistical optimization approach; for description see Sect.2.1 and 2.2.

Figure 2 .
Figure 2. Variability of relative standard deviations of forecast-minus-analysis bending angle differences (upper two panels) and of the systematic uncertainty of the mean analysis bending angle (bottom two panels) as function of latitude (left) and of day of month (right), respectively.

Figure 3 .
Figure 3. Global-mean error correlation functions from the background error covariance matrix (left), for the 5, 15, and 25 July 2008 at three representative impact altitude levels (30, 50, and 70 km), and estimated correlation lengths of the correlation functions (right) at all impact altitude levels from 20 to 80 km for the same 3 days.

Figure 4
Figure 4 illustrates the observation uncertainty profile u occ o , and intermediate variables from Eqs. (7) and (8), for six example RO events from 15 July 2008.The simMetOp events are simulated using the End-to-end GNSS Occultation Performance Simulation and Processing System software version v5.6 (EGOPSv5.6) in the same way as the sim-MetOp ensemble was produced byLi et al. (2013); for more details see Sect. 3 below.In each panel, the red line shows the original difference profile α k o ("delta"), the blue line the smoothed difference profile α

Figure 5 .
Figure 5. Global-mean error correlation functions from the observation error covariance matrix (left), for the 5, 15, and 25 July 2008 at three representative impact altitude levels (30, 50, and 70 km), and estimated correlation lengths of the correlation functions (right) at all impact altitude levels from 20 to 80 km for the same 3 days.For ease of intercomparison the layout is the same as in Fig.3.

Figure 6 .
Figure 6. Background and observation bending angle uncertainty profiles as well as observation-to-background (Obs-to-Bgr) weighting ratio (left); statistically optimized bending angle profiles from the OPSv5.6,b-dynamic, dynamic, and CDAAC algorithms together with their reference profile (middle); and difference of the optimized profiles to the reference profile (right).Three example events from 15 July 2008 are illustrated, from simMetOp (top), CHAMP (middle), and COSMIC (bottom), respectively.

Figure 7 .
Figure 7. Statistically optimized bending angle profiles together with their reference profile (left) and their difference to the reference profile (right), of three example events from simMetOp (top), CHAMP (middle), and COSMIC (bottom) from 15 July 2008, using either the realistic global-mean correlation matrix of the new dynamic method ("full correlation") or simple exponential fall-off correlation as existing in OPSv5.6 ("exp.falloffonly").

Figure 8 .
Figure 8. Systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines) of statistically optimized bending angles, relative to "perfect" simulated bending angles or co-located ECMWF analysis bending angles used as a reference, of the global ensemble of simMetOp events on 15 January and 15 July 2008 (upper two panels), and of CHAMP and COSMIC events from the complete months of January and July 2008 (middle and bottom panel, respectively).Statistics of the OPSv5.6 (black), b-dynamic (blue), dynamic (red), CDAAC (version 2009.2650for CHAMP and version 2010.2640 for COSMIC, green), and CDAAC new (version 2014.0140 for CHAMP, magenta) statistical optimization methods are shown.The number of events (NoE) used in the ensemble of each statistical calculation is also indicated in each panel.

Figure 10 .
Figure 10.Systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines) of statistically optimized bending angles, relative to "perfect" simulated bending angles used as a reference, of simMetOp events on 15 July 2008.Statistics for the OPSv5.6 (black), b-dynamic (blue), and dynamic (red) statistical optimization algorithms are shown for six different regions: Global (90 • S to 90 • N), TRO (tropics, 20 • S to 20 • N), SHSM (Southern Hemisphere subtropics and midlatitudes, 20 to 60 • S), NHSM (Northern Hemisphere subtropics and midlatitudes, 20 to 60 • N), SHP (Southern Hemisphere polar region, 60 to 90 • S), and NHP (Northern Hemisphere polar region, 60 to 90 • N).The number of events (NoE) used in the ensemble of each region is also indicated in the panels.

Figure 11 .
Figure 11.Systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines) of statistically optimized bending angles, relative to co-located ECMWF analysis bending angles used as a reference, of CHAMP events from July 2008.Statistics for the OPSv5.6 (black), b-dynamic (blue), dynamic (red), CDAAC (version 2009.2650,green), and CDAAC new (version 2014.014,magenta) statistical optimization algorithms are shown for the same six regions as in Fig. 10.The figure layout is the same as for Fig. 10.

Figure 12 .
Figure 12.Systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines) of statistically optimized bending angles, relative to co-located ECMWF analysis bending angles used as a reference, of COSMIC events from July 2008.Statistics for the OPSv5.6 (black), b-dynamic (blue), dynamic (red), and CDAAC (version 2010.2640,green) statistical optimization algorithms are shown for the same six regions as in Fig. 10.The figure layout is the same as for Figs. 10 and 11.

Figure 13 .
Figure 13.Systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines) of retrieved refractivity profiles, relative to "perfect" simulated refractivity or co-located ECMWF analysis refractivity used as a reference, for the global ensemble of simMetOp events on 15 January and 15 July 2008 (top panels) and of CHAMP events (middle panels) and COSMIC events (bottom panels) from the complete months of January and July 2008.Statistics of the OPSv5.6 (black), b-dynamic (blue), dynamic (red), CDAAC (version 2009.2650for CHAMP and version 2010.2640 for COSMIC, green), and CDAACnew (version 2014.0140 for CHAMP, magenta) statistical optimization methods are shown.The figure layout is the same as for Fig. 8.

Figure 14 .
Figure 14.Systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines) of retrieved temperature profiles, relative to "perfect" simulated temperature or co-located ECMWF analysis temperature used as a reference, for the global ensemble of simMetOp events on 15 January and 15 July 2008 (top panels) and of CHAMP events (middle panels) and COSMIC events (bottom panels) from the full months of January and July 2008.The figure layout is the same as for Figs. 8 and 13.

Figure
Figure15.Bending angle (left), refractivity (middle), and temperature (right) systematic differences (SysDiff, light lines) and standard deviations (SD, heavy lines), relative to their "perfect simulated" or co-located ECMWF analysis data used as a reference, of the global ensemble of simMetOp (top) and COSMIC (bottom) events from 15 July 2008, using either the realistic global-mean correlation matrix of the new dynamic method ("full correlation") or simple exponential fall-off correlation as in the existing OPSv5.6 ("exp.falloffonly").The number of events (NoE) in each statistical ensemble is also indicated in the panels.