A simplified approach for generating GNSS radio occultation refractivity climatologies

Abstract. The possibility of simplifying the retrieval scheme required to produce GNSS radio occultation refractivity climatologies is investigated. In a new, simplified retrieval approach, the main statistical analysis is performed in bending angle space and an estimate of the average bending angle profile is then propagated through an Abel transform. The average is composed of means and medians of ionospheric corrected bending angles up to 80 km. Above that, the observed profile is exponentially extrapolated to infinity using a fixed a priori scale height. The new approach circumvents the need to introduce a "statistical optimisation" processing step in which individual bending angle profiles are merged with a priori data, often taken from a climatology. This processing step can be complex, difficult to interpret, and is generally recognised as a potential source of structural uncertainty. The new scheme is compared with the more conventional approach of averaging individual refractivity profiles, produced with the implementation of statistical optimisation used in the EUMETSAT Radio Occultation Meteorology Satellite Application Facility (ROM SAF) operational processing. It is shown that the two GNSS radio occultation climatologies agree to within 0.1% from 5 km up to 35–40 km, for the three months January, February, and March 2011. During this time period, the new approach also produces slightly better agreement with ECMWF analyses between 40–50 km, which is encouraging. The possible limitations of the new approach caused by mean residual ionospheric errors and low observation numbers are discussed briefly, and areas for future work are suggested.


Introduction
GNSS radio occultation (GNSS-RO) measurements are likely to have an increasingly important role in climate monitoring and climate reanalyses in the coming years, as the duration of the time series lengthens.They provide accurate geophysical information with a high vertical resolution at altitudes from around 10 km to 40 km.The GNSS-RO measurement characteristics, and details of the retrieval scheme required to produce geophysical variables from the measurements, have been described by Kursinski et al. (1997).We emphasise that the retrieval of geophysical variables from the raw measurements requires the use of a priori information.For example, the retrieval of refractivity profile information from bending angle profiles requires the extrapolation of the bending angle profile to infinity with an a priori model of the bending angles.The need for extrapolation is evident from the Abel transform, relating refractive index, n, to ray bending angle, α, through, where a is the impact parameter and x = nr, with r being the radius of a point on the ray path.The upper limit of the integral necessitates extrapolation of the observed bending angle profile which is finite.In addition, the signal to noise of the observed bending angle profile falls as the tangent height of the bending angle increases.This is because to first order the bending angle values fall exponentially with tangent height, but the measurement error in the stratosphere is relatively constant with height.For example, the bending angle at Published by Copernicus Publications on behalf of the European Geosciences Union.

122
H. Gleisner and S. B. Healy: Simplified generation of RO refractivity climatologies a tangent height of 60 km is typically ∼ 5 microradians, and the measurement noise is ∼ 1-2 microradians.Therefore, the extrapolation step and merging with the a priori bending angle information is often combined with a smoothing of the retrieved bending angles over an extended vertical interval usually above ∼ 40 km, in a processing step often referred to as "statistical optimisation".The "optimised" bending angle profiles are then used to retrieve refractivity profiles and, subsequently, profiles of temperature, geopotential height and pressure.The climatologies of the geophysical parameters are commonly derived by binning and averaging the individual retrieved profiles.
Most of the GNSS-RO processing centres have implemented some form of statistical optimisation with varying degrees of complexity.The specific details of the implementations are outlined in Ho et al. (2009Ho et al. ( , 2012)).In general, the statistical optimisation processing step can be written in a matrix/vector form as where α, α b and α o are the vectors of optimised, background (a priori) and observed bending angles, respectively, and K is the gain matrix, which can be written in terms of the error covariance matrices assumed for the observed and background bending angle profiles.The ROtrends project, in which GNSS-RO climatologies produced by different processing centres are compared (Ho et al., 2009(Ho et al., , 2012;;Steiner et al., 2012), has shown that that differences in the statistical optimisation are a major source of "structural uncertainty".This has also been noted in previous studies (von Engeln, 2006).Structural uncertainty refers to errors that arise as a result of subjective choices made in the processing of the observations (e.g., Thorne et al., 2005), which may vary in between the processing centres.The implementations of the statistical optimisation at the processing centres differ because of differences in the a priori atmospheric state information used to estimate α b , and in the assumed error covariance matrices used in the gain matrix K.
Although it may be necessary to implement quite complex statistical optimisation routines for retrieving individual profile information from GNSS-RO measurements, we will show that this is not the case when producing climatologies.The statistical optimisation complicates the interpretation of the GNSS-RO climatologies of the geophysical parameters and leads to structural uncertainty.An alternative approach, which circumvents the statistical optimisation processing step, and simplifies the retrieval, is to perform averaging in bending angle space, and then propagate the estimate of the averaged bending angle information through the retrieval to produce the climatologies of the geophysical variables.If the retrieval scheme is linear, this should produce the same results as averaging the individual geophysical profiles, so differences in the climatologies are a useful indication of non-linearity in the retrieval scheme, and may be regarded as a way of examining the influence of the statistical optimisation in the processing.The averaged bending angle information should be considerably smoother than individual profiles, so the smoothing introduced by the statistical optimisation should not be required.Although it must be emphasised that the new approach still requires the use of a priori information, because of the extrapolation of the averaged bending angle profiles, this can be done in a relatively simple manner, and the sensitivity of the refractivity values below 40 km to the a priori is small.Ringer and Healy (2008) suggested monitoring climate signals directly in bending angle space in order to remove the problems associated with statistical optimisation, but some users prefer to work with retrieved geophysical parameters.Therefore, the present study investigates the use of averaged bending angle profiles to produce the GNSS-RO refractivity climatology, and compares this new approach with the profile by profile approach using the climate retrieval scheme used at the DMI, which includes an implementation of statistical optimisation.It is demonstrated that the new approach produces essentially the same results -to within ∼ 0.1 %in the 5-35 km vertical interval, where the information content of GNSS-RO is largest.This demonstrates that statistical optimisation is not required when generating climatologies from GNSS-RO measurements, and that the processing of the observations can be simplified considerably.We note that a similar independent study using the Jet Propulsion Laboratory GNSS-RO retrieval system has been published recently (Ao et al., 2012).While the averaging of the bending angle profiles is treated somewhat differently from the present study, the general conclusions concerning the applicability of the average-profile approach are similar.
The data used in this study are described in Sect. 2. Section 3 describes the processing of the data to refractivity.The main results, comparing two approaches for generating refractivity climatologies, are given in Sect.4, and the discussion and conclusions are given Sect. 5.

Data
The present study is based on occultations observed by the GNSS-RO instruments onboard the FORMOSAT-3/COSMIC satellites (henceforth referred to as COSMIC) during the first 3 months of 2011.The excess-phase and satellite orbit data were obtained from the CDAAC database at UCAR, and further processed into bending angles and refractivities at the RO Meteorology Satellite Application Facility (ROM SAF) which is a decentralised operational RO processing centre under EUMETSAT, hosted by the DMI.
The CDAAC database contains about 125 000 occultations with excess phase data for the studied 3-month time period.After processing and several steps of quality screening, 84 % of these profiles remain: about 30 000 in January, 30 000 in February, and 45 000 in March 2011.The data are distributed into 5-degree latitudinal bins based on a nominal latitude assigned to each occultation.
In the evaluation of the results, we also use refractivity climatologies generated from ECMWF analysis refractivity profiles co-located with the observed profiles.ECMWF analysis data mapped to bending angle space are also used as a reference in Sect.3.3 when discussing upper-level bending angle distributions.The ECMWF bending angles above about 60 km depend not only on the model data itself, but also on extrapolation above the model top.However, this is of less importance since we only use the ECMWF bending angle profiles as a reference to describe the variability of the monthly mean and median profiles.

The Abel transform
The Abel transform (Eq. 1) relates the refractive index, n, to ray bending angle, α.The integral in Eq. ( 1) is solved by assuming that α can be approximated as a linear function of impact parameter, a, between successive observation levels j where α j is the bending angle at impact parameter a j .The sub-integrals, can then be solved analytically, The refractivity (defined as N = 10 6 (n − 1)) at height x is given by the sum of contributions from the successive atmospheric layers, from height x to the top of the atmosphere.Through the approximation ln(n) ≈ 10 −6 N, the Abel transform becomes fully linear in α (see Appendix).
The errors due to the discretisation of the Abel transform contribute to the structural uncertainty.Studies of the consistency of forward and inverse Abel transforms show that the errors associated with the discretisation scheme described above are on the order of 0.1 % in the stratosphere (Lewis, 2008).It should be noted that this is not an inherent limitation of the Abel inversion, but a consequence of the chosen discretisation.

Individual bending angle profile processing
The standard procedure to retrieve refractivity from noisy bending angles is to smooth and merge the observed bending angle profiles individually with a priori data taken from a climatology.This is often a complex processing step that involves methodological choices and determination of various parameters.This section briefly outlines the method used in the ROM SAF operational system in order to highlight the complexity of the approach.The next three paragraphs summarise the more detailed description given in Lauritsen et al. (2011).
The implementation of statistical optimisation (SO) is based on the optimal linear combination (OLC) approach suggested by Gorbunov (2002).The method combines the ionospheric correction and the SO in a single least squares processing step.
The background profile used in the SO is found through a global search of a small library of MSIS-90 profiles (Hedin, 1991) mapped to bending angle space.Each of the MSIS bending angle profiles is individually scaled and shifted in log(bending angle) space, in order to fit -in a least squared sense -a smoothed version of the observed ionospheric corrected bending angle profile between 20 and 70 km.The ionospheric corrected bending angle profile used here is given by the standard linear combination (LC) of the L1 and L2 bending angles (Gorbunov, 2002).The least square fit is only performed at altitudes where the MSIS bending angle deviate by less than 30 % from the LC bending angle.The profile that minimises the magnitude of the two retrieved scaling and shifting parameters is selected.The method differs somewhat from the approaches outlined by Lohmann (2005) and Gobiet and Kirchengast (2004).
The error variance assigned to the selected a priori bending angle profile is estimated from the mean deviation of this profile from the observed LC bending angle in the 12-35 km height interval.The ionospheric noise variance is estimated from the upper part of a strongly smoothed version of the L1-L2 signal.These variance estimates are then used to solve a set of linear equations, corresponding to Eq. ( 2), which gives the statistically optimised bending angles, noting that all vertical error correlations are ignored in this calculation.
The optimised bending angles undergo an Abel transform by computing Eqs. ( 5) and ( 6).The retrieved refractivity profiles, N (x), are converted to functions of mean-sea level altitude, N (H ), and then interpolated to a regular altitude grid.Finally, they are averaged within monthly latitude bins to form the zonal monthly means, N (H ).These monthly means, generated by the ROM SAF operational system, are used in Sect. 4 as a reference in the evaluation of the mean bending angle processing described below.

Mean bending angle profile processing
Rather than handling the noise in the bending angles separately for each profile, the noise may be suppressed by averaging a large number of profiles.This reduces the random errors and raises the altitude where instrumental errors and ionospheric residuals become significant.The extrapolation of the mean bending angle profiles to infinity could be done with statistical optimisation, similar to the single-profile processing.However, in this study we explore a simpler approach that involves using observational data up to 80 km from where the bending angles are exponentially extrapolated.
As a result of neutral-atmosphere variability, instrumental errors, and ionospheric residuals there is substantial profileto-profile variation of the bending angles.Averaging suppresses this variability but above about 50 to 60 km, the mean bending angle profiles exhibit large-scale wiggles, occasionally leading to negative bending angles (Fig. 1a).The median bending angles do not show the corresponding wiggles, although smaller-scale variations are evident (Fig. 1b).The large variability of the mean bending angle profiles, and the significantly smaller variability of the median profiles, are seen more clearly in Fig. 2, where we have subtracted ECMWF analysis data mapped to bending angle space.Here, the ECMWF data are only used as a reference from which to measure variability.
The relative difference between mean and median is an indicator of the skewness of the bending angle distribution.In the troposphere and lower stratosphere, similar tendencies are seen in the ECMWF analyses mapped to bending angle space and in observed data (Fig. 3), indicating that the observed mean-median differences have a real neutralatmosphere origin, rather than being a consequence of measurement errors.In contrast, the observed bending angle distributions above 50 to 60 km (Fig. 3) are strongly influenced by large measurement errors, such as residual ionospheric errors.The mean value can no longer be considered a robust estimate of the central value of the neutral-atmosphere bending angle distribution.In the presence of strongly deviating bending angles, the median value provides a more robust estimate of the central value of the bending angle distribution (Press et al., 1992), which is consistent with the results in Figs.1-3.Ao et al. (2012) have not discussed the robustness of the mean bending angle estimates in their work, and we note that this may account for the wiggles in the observed mean profile shown in their Fig.4a.
Consequently, we use mean bending angles only up to 50 km altitude.In the height interval between 50 and 60 km, the mean values are merged with the median values through a simple linear combination, and from 60 km up to 80 km we use the median values.Above 80 km, we assume an exponential fall-off with a constant scale height of 7.5 km.This approach of extrapolating the observed bending angles with a fixed a priori scale height appears to be sufficiently accurate  when we are primarily concerned with refractivity below 40 km, where the GNSS-RO information content is largest.We have also tested scale heights of 5 km and 10 km, and found that they have negligible impact on the refractivity below values 40 km.A simple analytical estimate shows that the ray bending above 80 km only contributes about 0.1 % to the refractivity at 40 km altitude (see Appendix).Hence, realistic errors made in the exponential extrapolation to infinity only have small effects in the stratosphere.However, we will revisit the value scale height, in the context of retrieving temperature profiles from the the refractivity information, in future work.
The bending angle means and medians are generated from bending angle profiles interpolated to a common impact altitude grid.Before inverting an average bending angle profile with the Abel integral (Eq.1), we must change the height variable from impact altitude, H a , to impact parameter, a. Impact altitude is related to impact parameter through where R c is the local radius of curvature of the reference ellipsoid and u is the undulation, i.e. the height of geoid above the reference ellipsoid.Even though each profile within a bin has its own radius of curvature, we convert the height scale using a single radius of curvature for the whole bin.This is chosen as the average of the local radius of curvature of the geoid, where R c,i + u i is the radius of curvature for occultation i and m is the number of occultations in the bin.The height conversion using a mean radius of curvature, R c , assumes that the mean bending angle, α(H a ), at impact altitude H a is identical to the mean bending angle, α(a) , at impact parameter a.The spread in impact parameters for a fixed impact altitude, in combination with a nonlinear α(a) relation, leads to a bias in the mean bending angle.

Atmos
The output from the Abel transform (Eq. 1) is mean refractivity, N (x), as a function of x = nr.This is converted to mean refractivity as a function of altitude, H , using the radius of curvature given by Eq. ( 8), After interpolation to a regular altitude grid, these mean refractivities can be directly compared to mean refractivities obtained from single-profile processing.Any differences are due either to nonlinearities in the processing, or to the different upper-level handling of the bending angles.

Results
Zonal monthly mean refractivities for the first 3 months of 2011 were computed using the ROM SAF standard singleprofile processing and the new average-profile processing.
Both refractivity climate data sets were generated from the same set of COSMIC bending angle profiles binned into 5degree latitude zones.A corresponding set of monthly mean In the troposphere and lower stratosphere, the observed mean-median differences are very similar to the collocated ECMWF data, indicating a real neutral-atmosphere origin.In contrast, the observed mean-median differences above 50 to 60 km are strongly influenced by large measurement errors.

Fig. 3.
Relative differences between bending angle means and medians for observed COSMIC data (upper panel), and collocated analysis data (lower panel).Data are from January 2011 binned into 5-degree latitude bands.In the troposphere and lower stratosphere, the observed mean-median differences are very similar to the collocated ECMWF data, indicating a real neutralatmosphere origin.In contrast, the observed mean-median differences above 50 to 60 km are strongly influenced by large measurement errors.
refractivities was generated from ECMWF analysis data colocated with the COSMIC occultations.
Figure 4 shows the relative differences between the monthly mean refractivities obtained from inverting the average bending angles, and those obtained from inverting individually optimised bending angle profiles.The relative differences at 40 km altitude are on the order of 0.1 %, and at 50 km the differences are around 1 %.The differences appear to have a latitudinal structure, being positive at low latitudes and predominantly negative at higher latitudes.At northern hemisphere high latitudes, the differences are more variable than at lower latitudes and in the Southern Hemisphere.This may be related to the high bending angle variability in the winter hemisphere.
In the lower troposphere the differences become larger, reaching up to around 1 % nearest to the surface.The latter is most likely a consequence of non-linearities due to the use of a mean radius of curvature in the conversion of height scales within each latitude bin (see Sect. 3.3).
As a comparison, the structural uncertainties related to the discretisation of the Abel integral are on the order of 0.1 % (see Sect. 3.1).The weak nonlinearities of the Abel inversion itself (Eq. 1) cause errors that are an order of magnitude smaller, decreasing exponentially with altitude at the density scale height from a maximum value of around 0.02 % nearest to the surface in the tropics (see Appendix).
Figures 5-7 show the relative differences between monthly mean refractivities from COSMIC data and from co-located  ECMWF data, using average-profile inversion (upper panels) and single-profile inversion (lower panels).Above the lowest few kilometers in the troposphere and below about 35 km in the stratosphere, the differences with respect to ECMWF analysis data are almost identical for the two inversion methods.This indicates that for the studied time period, the singleprofile statistical optimisation has not caused any significant effects of the a priori below 35 km.Another side of this is that there are no obvious benefits of the additional complexity introduced into the processing by the statistical optimisation, compared to the new simplified approach.In the height range 35 to 50 km, there are small but discernable differences of up to 1 % between the two methods.Routine monitoring of bending angle departure statistics at ECMWF generally shows a positive bending angle bias relative to ECMWF analyses, peaking near 40 km for all GNSS-RO instruments.This bending angle bias contributes to the refractivity bias near 40 km, shown in Figs.5-7.During the studied 3-month period, the co-located ECMWF analysis data appear to agree slightly better with the refractivities obtained from the average bending angles, than with the refractivities obtained from the standard single-profile processing scheme.

Discussion and conclusions
One of the motivations for inverting the average bending angles, rather than the individually optimised bending angle profiles, is to remove a source of structural uncertainty from the processing.It has been shown that differences in the bending angle optimisation procedures constitute a source of systematic differences between GNSS-RO processing centres, particularly above 25 km and at high latitudes (Steiner et al., 2012).
While there may be compelling reasons for applying statistical optimisation to retrieve geophysical profile information from individual bending angle profiles, we have shown that this is not the case when generating zonal monthly mean refractivities.We find that the average-profile processing and the single-profile processing produce monthly mean refractivities that are nearly identical in the stratosphere up to around 35-40 km.In the 40-50 km altitude range, we find relative differences between the two inversion methods up to around 1.0 %.For the 3 months included in this study, the average-profile inversions tend to agree slightly better with monthly mean refractivities from ECMWF analysis data than the single-profile inversions, although more extensive studies One possible area of concern is whether the new approach is more sensitive to residual ionospheric biases, because of the use of observed bending angle values up to 80 km.We do not see any evidence of this potential problem below 40 km in the results presented in this study, but further work will be required to confirm this at solar maximum conditions.High solar activity is likely to produce noisier individual bending angle profiles, and consequently less smooth mean bending angle profiles.However, note that SO schemes are designed to reduce random errors and extrapolate the bending angle profiles to infinity, rather than reduce measurement bias.Furthermore, some SO schemes include an explicit adjustment of the background climatology to the observed bending angle profile, meaning that residual ionospheric errors could still cause biases in the SO bending angles and refractivity profile.
Another consideration is whether the average-profile approach will be sufficiently accurate to produce refractivity climatologies prior to 2006, before the launch of the COS-MIC satellites which led to an order of magnitude increase in the GNSS-RO data numbers.We have performed a preliminary investigation to test how robust the new approach is by randomly removing 85 % of the data and repeating the study with the remaining 15 % of the data.The results (figures not shown) indicate that the mean refractivity differences below 40 km are similar to those shown in Figs.4-7.These preliminary results are encouraging, but we have not yet conclusively demonstrated that the technique will work with just CHAMP measurements at solar maximum conditions.This will be done in future work.In this context it is important to note that the uncertainty in any GNSS-RO refractivity climatology during the period prior to the COSMIC mission will be greater than during the COSMIC mission irrespective of the processing method, because of the large change in data numbers.The key issue to be investigated is whether the average profile processing increases the uncertainty when data numbers are reduced, when compared with the standard approach to producing the climatology.This question requires further investigation.
The presented average-profile approach is mainly a tool for generating refractivity climatologies in the stratosphere.In the lowest few kilometers, nearest to the surface, the effects of nonlinearities introduced by the conversion of height scales become evident.The remedy for this may be to compute a correction based on the known profile-to-profile variability of the atmosphere's radius of curvature.The optimal way of doing that needs to be studied.
Inverting average bending angle profiles is a much simpler approach than the often quite complex implementations of SO schemes.Shorter execution times, more robust implementation of processing software, and the simplification of error propagation are definitive advantages.It also affords a test of the standard retrieval schemes for potential impacts from the SO step.The findings by Steiner et al. (2012), showing that systematic differences above 25 km and at high latitudes can be attributed to differences in the SO implementations, suggests that this testing would be useful.
We conclude that it is possible to retrieve monthly mean refractivity profiles directly from average bending angles without statistical optimisation, using a very simple a priori bending angle model above 80 km.These results are very encouraging, and suggest that much simpler retrieval systems than are currently used can be employed when producing refractivity climatologies from GNSS-RO measurements.

Fig. 1 .
Fig. 1.(a) Observed bending angle means, and (b) observed bending angle medians, for COSMIC data from January 2011 within 5-degree latitudinal bins.Above 55 to 60 kilometers an increasing impact of measurement errors shows up as large-scale wiggles in the means, occasionally leading to negative bending angles.The medians do not show the corresponding wiggles, although higher-frequency variations are evident.

Fig. 2 .Fig. 1 .
Fig. 2. Mean bending angle profiles (left panel), and median bending angle profiles (right panel), for January 2011 within 5-degree latitudinal bins, after subtraction of ECMWF analysis data mapped to bending angle space.Above 50 to 60 kilometers, the lower variability of the median profiles indicates that they provide a more robust measure of the central value of the bending angle distribution.15

Fig. 2 .
Fig. 2.Mean bending angle profiles (left panel), and median bending angle profiles (right panel), for January 2011 within 5-degree latitudinal bins, after subtraction of ECMWF analysis data mapped to bending angle space.Above 50 to 60 km, the lower variability of the median profiles indicates that they provide a more robust measure of the central value of the bending angle distribution.

Fig. 3 .
Fig.3.Relative differences between bending angle means and medians for observed COSMIC data (upper panel), and collocated ECMWF analysis data (lower panel).Data are from January 2011 binned into 5-degree latitude bands.In the troposphere and lower stratosphere, the observed mean-median differences are very

Fig. 4 .
Fig. 4. Relative differences between monthly mean refractivities obtained by average-profile inversion and single-profile inversion.Data for the first three months of 2011.

Fig. 5 .Fig. 5 .
Fig. 5. Relative differences between monthly mean refractivity from COSMIC data and from collocated ECMWF data, using average-profile inversion (upper panel) and single-profile inversion including statistical optimization (lower panel).Data for January 2011.

Fig. 6 .Fig. 6 .
Fig. 6.Relative differences between monthly mean refractivity from COSMIC data and from collocated ECMWF data, using average-profile inversion (upper panel) and single-profile inversion including statistical optimization (lower panel).Data for February 2011.

Fig. 7 .Fig. 7 .
Fig. 7. Relative differences between monthly mean refractivity from COSMIC data and from collocated ECMWF data, using average-profile inversion (upper panel) and single-profile inversion including statistical optimization (lower panel).Data for March 2011.