Interactive comment on “ Generation of a Bending Angle Radio Occultation Climatology ( BAROCLIM ) and its use in radio occultation retrievals ”

Abstract. In this paper, we introduce a bending angle radio occultation climatology (BAROCLIM) based on Formosat-3/COSMIC (F3C) data. This climatology represents the monthly-mean atmospheric state from 2006 to 2012. Bending angles from radio occultation (RO) measurements are obtained from the accumulation of the change in the raypath direction of Global Positioning System (GPS) signals. Best quality of these near-vertical profiles is found from the middle troposphere up to the mesosphere. Beside RO bending angles we also use data from the Mass Spectrometer and Incoherent Scatter Radar (MSIS) model (modified for RO purposes) to expand BAROCLIM in a spectral model, which (theoretically) reaches from the surface up to infinity. Due to the very high quality of BAROCLIM up to the mesosphere, it can be used to detect deficiencies in current state-of-the-art analysis and reanalysis products from numerical weather prediction (NWP) centers. For bending angles derived from European Centre for Medium-Range Weather Forecasts (ECMWF) analysis fields from 2006 to 2012, e.g., we find a positive bias of 0.5 to 1% at 40 km, which increases to more than 2% at 50 km. BAROCLIM can also be used as a priori information in RO profile retrievals. In contrast to other a priori information (i.e., MSIS) we find that the use of BAROCLIM better preserves the mean of raw RO measurements. Global statistics of statistically optimized bending angle and refractivity profiles also confirm that BAROCLIM outperforms MSIS. These results clearly demonstrate the utility of BAROCLIM.


Introduction
Global data sets of the lower and middle atmosphere (troposphere to upper mesosphere) provide important information to understand atmospheric dynamics of the Earth's climate system.Observational data as well as analysis/reanalysis data and atmospheric models are used, e.g., to study specific atmospheric phenomena such as MJO (Madden Julian Oscillation), ENSO (El Niño-Southern Oscillation), or the QBO (Quasi Biennial Oscillation).Long-term observational records, reanalysis data sets, and atmospheric models can also be used to investigate atmospheric climate change.
Simple empirical atmospheric models are of high utility if a quick but reasonable estimate of the atmospheric state is of main interest.This is of importance, e.g., for simulation studies in the field of atmospheric remote sensing or within the retrieval of atmospheric parameters from remote sensing measurements.For this purpose several research communities use early empirical models like CIRA (Committee on Space Research (COSPAR) International Reference Atmosphere) (Fleming et al., 1990) or MSIS (Mass Spectrometer and Incoherent Scatter Radar) (Hedin, 1991;Picone et al., 2002).
Published in the early 1960s, the CIRA model was the earliest comprehensive climatological model of the atmosphere to contain information up to the thermosphere.This model is based on observational data such as radiosondes, rocket data, and satellite observations.Its fourth version, CIRA-86, includes thermosphere models as well as tables of monthlymean zonal-mean temperature, pressure, geopotential height, and zonal wind from the surface to an altitude of 120 km.The most recent CIRA version (CIRA-2012)  dated versions of empirical models of the Earth's upper atmosphere (above 120 km) (CIRA, 2012).
Recent MSIS model versions (MSIS-90 and NRLMSIS-00) provide information on atmospheric composition, total mass density, and temperature from the ground up to the exosphere also using observations from ground, rockets, and satellites.The MSIS model output depends on time (day of year, universal time, local solar time), location (altitude, latitude, longitude), geomagnetic activity (represented by the magnetic index A p ), and solar activity (represented by the 10.7 cm solar radiation flux, F 10.7 ).
An overview on principal features of a number of global and regional models of the middle atmosphere and thermosphere is given by AIAA ( 2004), summarizing their model content, uncertainties, and limitations.Randel et al. (2004) compared middle-atmosphere climatologies from approximately 10 to 80 km -including historical measurements from rocketsonde winds and temperatures (1970 to 1989), lidar temperature data (1990s), global meteorological analyses, and satellite data -and found notable differences between these data sets.Some biases found in atmospheric analyses were caused by the low vertical resolution of these data as well as the low vertical resolution of some assimilated satellite measurements.
Since 2004 research initiatives have tried to understand and eliminate errors of previous middle-atmosphere models, building new, state-of-the-art analysis and reanalysis products.However, there are still uncertainties and differences in current reanalysis products, and the SPARC (Stratospheretroposphere Processes And their Role in Climate) Reanalysis/analysis Intercomparison Project (S-RIP; Fujiwara et al., 2012) aims at understanding and reasonably interpreting these differences and contributing to future reanalysis improvements in the middle atmosphere.
In this study, we aim at compiling and investigating a global climatological model from recent high-resolution radio occultation (RO) bending angles.
The RO method (Melbourne et al., 1994;Kursinski et al., 1997) is an active limb sounding technique that utilizes radio signals continuously broadcast by Global Positioning System (GPS) satellites.The measured quantity is the phase change of the GPS signal as a function of time, which is a measure for physical atmospheric parameters, in particular bending angle and radio refractivity, from which density, pressure, geopotential height, temperature, and humidity profiles can be derived with very high accuracy (though in the lower moist troposphere auxiliary information is necessary to separate dry-air and moist-air contributions to the refractivity).Precise and stable oscillators aboard the GPS satellites ensure measurement stability and consistency between various RO missions, without the need for instrument-dependent calibrations (Hajj et al., 2004;Schreiner et al., 2007;Foelsche et al., 2011).
The RO bending angle can be used as a climate variable (Ringer and Healy, 2008).In the upper troposphere and above (where moisture is negligible), it mainly depends on atmospheric density, which again depends on pressure and temperature.This means that a change in bending angle can be caused by changes of all these parameters.This has to be kept in mind when analyzing atmospheric variability in terms of bending angle.However, the use of bending angle as a climate variable is superior to other RO-derived atmospheric quantities that are sensitive to the additional use of a priori information used in the retrieval (see below).
First RO profiles of the Earth's atmosphere were provided by the GPS/MET experiment in intermittent periods from 1995 and 1997 (Rocken et al., 1997) 2013), and KOMPSAT-5 (2013) as well as from the first six-satellite RO constellation, Formosat-3/COSMIC (F3C;launched in 2006).
Atmospheric profiles from RO are used for data assimilation in numerical weather prediction (NWP) (e.g., Healy and Thépaut, 2006;Cucurull and Derber, 2008) and in atmospheric and climate research (see Anthes, 2011;Steiner et al., 2011, for reviews).Most studies utilize RO profiles in the upper troposphere and lower stratosphere (UTLS) region, i.e., the altitude range from approximately 5 to 35 km, where RO profiles of pressure, geopotential height, and temperature are known to be of best quality (Scherllin-Pirscher et al., 2011b).
The top altitude of RO measurements depends on the instrument settings but is at least 80 km (e.g., Metop-A and Metop-B).For F3C it is usually about 120 km or higher.However, at these altitudes, individual RO profiles are dominated by measurement noise and ionospheric disturbances, because neutral atmospheric density gradients are small.To derive refractivity profiles via the Abel integral transform (Fjeldbo et al., 1971), and since noisy and erroneous information from high altitudes propagates downwards in the retrieval process, bending angles are combined with a priori information, e.g., from a climatology, usually applying statistical optimization (Rodgers, 2000;Healy, 2001;Gorbunov, 2002;Gobiet and Kirchengast, 2004;Lohmann, 2005).
Formally, the upper limit in the Abel integral transform is infinity (Fjeldbo et al., 1971).If a climatology is used for statistical optimization, it is therefore necessary that it is able to provide a good measure of the bending angle to very high altitudes.Usually, the climatology is used as an extrapolation above the altitude range where the statistical optimization is performed (or some other form of extrapolation is necessary for high-accuracy purposes).However, biases in the a priori information are inevitable.At high enough altitudes biases are non-negligible (in a fractional sense) and affect the quality of derived atmospheric profiles (Ho et al., 2012;Steiner et al., 2013).Therefore high-quality a priori information is of particular importance.Ao et al. (2012) and Gleisner and Healy (2013) introduced a new approach to derive GPS RO climatological products.Instead of averaging individual refractivity profiles, they calculated monthly-mean bending angles and computed mean refractivity as the Abel inversion of the mean bending angle.Ao et al. (2012) showed that the maximum altitude through which the F3C measurements are useful increased substantially, which leads to a reduced bias in climatological averages of refractivity.The main drawback of this approach, however, is that it can only be used to obtain mean atmospheric fields and is not applicable to individual RO profiles.
The generation of the bending angle radio occultation climatology (BAROCLIM) described in this paper is mainly based on the work by Foelsche and Scherllin-Pirscher (2012) and Scherllin-Pirscher (2013).BAROCLIM is based on measurements from 2006 to 2012.Since this period is much shorter than the standard period for climate normals (30 years), BAROCLIM does not represent a long-term mean.However, it is of very high utility to obtain a quick and reasonable estimate of the atmospheric mean state as required, e.g., within the retrieval of atmospheric parameters from remote sensing measurements.
Section 2 gives a detailed description of the RO data set as well as of reference data used in this study.In Sect. 3 we present the method used to construct BAROCLIM, and in Sect. 4 we discuss potential systematic errors and evaluate the model by comparing it to reference data from NWP analysis fields provided by the European Centre for Medium-Range Weather Forecasts (ECMWF).In Sect. 5 we show that BAROCLIM can further be used as a priori information for statistical optimization in RO profile retrievals.Conclusions are drawn in Sect.6.

Data
For the generation of BAROCLIM, we used ionospherecorrected bending angles as a function of impact altitude (impact parameter with the local radius of curvature and geoid undulation subtracted).These profiles were retrieved with the Wegener Center for Climate and Global Change (WEGC) Occultation Processing System version 5.6 (OPSv5.6)(Schwärz et al., 2013) and interpolated on a regular 200 m grid.
Input data to the WEGC processing system are profiles of excess phase and amplitude as well as precise orbit information of GPS and low Earth orbit (LEO) satellites (level 1a data) provided by other data centers.Currently WEGC uses level 1a data provided by University Corporation for Atmospheric Research (UCAR)/COSMIC Data Analysis and Archive Center (CDAAC) for all RO satellite missions.Since recent UCAR/CDAAC processing versions vary for different missions and GPS receivers used for RO measurements are not of the same quality (Foelsche et al., 2011, e.g., found different bending angle noise for different missions), we decided to use only data from the F3C constellation.The UCAR/CDAAC processing version of F3C data used in this study was 2010.2640 for the entire time period starting in 2006.
The processing system at UCAR/CDAAC relies on the Bernese software package to obtain precise orbits of F3C satellites and applies single differencing to remove F3C clock offsets (Schreiner et al., 2009).In the inversion retrieval, WEGC first corrects orbits for the Earth's oblateness (Syndergaard, 1998), removes outliers from L1 and L2 excess phase profiles, and smooths data with a regularization filter (Syndergaard, 1999).Doppler shift profiles are then calculated from excess phase profiles via numerical differentiation.In the upper troposphere and above, bending angles are computed from Doppler shift profiles based on geometric optics (Melbourne et al., 1994).Wave optics retrieval is applied in the lower and middle troposphere (Gorbunov, 2002;Gorbunov et al., 2004;Gorbunov and Lauritsen, 2004).To remove the ionospheric contribution to the bending angle, WEGC applies ionospheric correction on bending angle level following Vorob'ev and Krasil'nikova (1994) and Hocke et al. (2003).Except for the use of F3C open-loop data (Sokolovskiy et al., 2006) and the wave optics retrieval in the lower and middle troposphere, the WEGC OPSv5.6 bending angle retrieval is very similar to the OPSv5.4processing described by Ho et al. (2012) and Steiner et al. (2013).
We used F3C data from August 2006 to December 2012 (more than 6 years of data).The number of F3C measurements increased from 2006 to 2007 until all F3C spacecraft were in their final orbits (F3C/flight model no. 3 (FM-3) did not reach its final orbit altitude because solar panels were stuck, which limited the power and payload operation of this spacecraft).From 2007 to 2009 F3C always tracked more than 60 000 RO events per month (exception: June 2009, UCAR/CDAAC provided approximately 55 000 level 1a profiles).Since 2010 the number of measurements decreased again due to battery degradation of all spacecraft.Furthermore, F3C/FM-3 has been out of contact since 1 August 2010.However, the minimum number of level 1a data provided by UCAR/CDAAC per month in the period between 2006 and 2012 was always larger than 30 000.
Since RO profiles usually do not reach higher than about 120 km (only about 80 km for Metop) and because the true neutral atmospheric bending angle decreases nearly exponentially with height (and therefore measurement noise dominates in a fractional sense in the mesosphere and above), there is an upper limit to which RO data are useful for the generation of a climatology.To extend BAROCLIM to higher altitudes, we used bending angles based on a modified version of the MSIS-90 empirical model of the neutral atmosphere (Hedin, 1991).
The MSIS-90 model is based on data from several satellites, incoherent scatter radar stations, and rocket probes

B. Scherllin-Pirscher et al.: Generation of BAROCLIM and its use in RO retrievals
as well as from tabulations from the Middle Atmosphere Program (MAP) Handbook 16 (Hedin, 1983(Hedin, , 1987(Hedin, , 1991)).Høeg et al. (1995) modified the original MSIS-90 model by smoothing out a discontinuity at 72.5 km and fixing the local apparent solar time at 0 h, the solar radio flux at F 10.7 = 150×10 −22 W m −1 Hz −1 , and the magnetic index at A p = 4.
Subsequently, this modified version of the MSIS model was used to generate spectral models of refractivity and bending angle for the use of fast and efficient modeling and inversion of RO data.The spectral model of refractivity, N, was obtained from the MSIS pressure, p, and temperature, T (N = 77.6 p/T ), for the 15th of every month.The spectral model of bending angle was obtained from the spectral refractivity model via the Abel integral transform.Thus, the MSIS-based bending angles and refractivities contain no information on atmospheric humidity and are given as a function of month, impact height, latitude, and longitude.Both spectral models are based on Chebychev polynomials in the vertical and spherical harmonics in the horizontal and were constructed using an approach similar to the one we now use for BAROCLIM described in Sect.3.4.The MSIS spectral models were later incorporated into the End-to-end Generic Occultation Performance Simulator (EGOPS) (Fritzer et al., 2012) and the Radio Occultation Processing Package (ROPP) (Culverwell, 2013).
To validate BAROCLIM, we used operational ECMWF analysis fields at T42 horizontal resolution (comparable to RO horizontal resolution) and 91 vertical levels.Profiles were extracted at mean RO event location of those F3C measurements, which were used to construct BAROCLIM.We applied a forward model to derive ECMWF bending angles from refractivity (above the ECMWF model top, refractivities were extended with MSIS profiles scaled to fit the ECMWF model at high altitudes).
We will show in Sect. 5 that BAROCLIM can be used as a priori information for the statistical optimization in the processing of RO measurements.To investigate the performance of BAROCLIM as a priori information, we used level 1a data from CHAMP, GRACE-A, SAC-C, and F3C from January 2008 and July 2008 (level 1a data were again provided by UCAR/CDAAC; CHAMP data version was 2009.2650;other satellites had consistent data version 2010.2640) and applied the WEGC OPSv5.6 retrieval using BAROCLIM and MSIS as a priori information for bending angle initialization.Retrieved atmospheric profiles were then compared to operationally retrieved OPSv5.6 profiles that used ECMWF short-range forecasts as a priori in the statistical optimization (Schwärz et al., 2013).

BAROCLIM generation
Figure 1 summarizes the steps of the BAROCLIM generation.Using ionosphere-corrected bending angles as a function of impact altitude on a regular 200 m grid, we first ap-plied a quality control (QC) to identify and exclude bad profiles before averaging bending angles in latitude and altitude and calculating monthly-mean 10 • zonal-mean bending angles.These monthly means include data from 6 or 7 years.Since we aimed at generating a BAROCLIM spectral model, which (theoretically) reaches from the surface to infinity, we extended mean RO profiles with MSIS.We refer to these bending angles as the BAROCLIM discrete model.The BAROCLIM spectral model was then constructed by expansion into Chebychev polynomials and zonal harmonics.Below we describe these steps in detail.

Quality control of individual profiles
Some bending angles are very noisy and/or contain unphysical values.They could strongly affect the quality of the bending angle climatology if they were not properly excluded.To avoid entering these profiles in BAROCLIM, we only used OPSv5.6 profiles that passed the external QC performed at WEGC (Schwärz et al., 2013).This external QC comprises bending angle, refractivity, and temperature profiles, which are compared to co-located profiles from ECMWF analysis fields.The profile is flagged as bad if the difference between the RO and the ECMWF profile is larger than 20 % in bending angle, 10 % in refractivity, and/or 25 K in temperature.However, since these quality checks are only performed in the upper troposphere and lower stratosphere region up to 35 km, profiles passing the external QC can still be very noisy in the upper stratosphere and above.
The inspection of individual bending angle profiles indeed revealed that it is imperative to perform an additional QC.We therefore introduced a threefold approach for an additional outlier rejection.In a first step we rejected all profiles with bending angles smaller than −40 µrad or larger than +40 µrad in the altitude range from 50 to 80 km.In this altitude range, neutral atmospheric bending angles are usually smaller than 20 µrad, and bending angles smaller/larger than ±40 µrad result from very strong (unphysical) noise, possibly sometimes related to ionospheric scintillations (Zeng and Sokolovskiy, 2010).
To also detect very bad profiles below 50 km, we rejected all profiles with bending angles smaller than −20 µrad below 50 km in a second step.Profiles that were flagged as bad by these QC mechanisms were outlier profiles, which could degrade the quality of BAROCLIM.Inspecting remaining bending angle profiles after application of these first checks, we still found some very noisy bending angle profiles (top panel of Fig. 2) and therefore applied an additional QC.
We used all remaining profiles and calculated mean bending angles and standard deviations for 10 • latitude bands for multiyear monthly means and rejected all profiles with bending angles outside of 4 standard deviations (4σ ) from the mean in the altitude range from the surface to 100 km.Fig application of this final QC results in a considerable decrease in standard deviation.
Table 1 gives an overview on the number of profiles provided by UCAR/CDAAC, retrieved at WEGC, and passing BAROCLIM QC.The number of bending angle profiles for a given month used to generate BAROCLIM is always larger than 200 000, except for June, when it is slightly below.In some months, however, the number is even larger than 260 000.These numbers are large enough and sufficient to obtain a smooth BAROCLIM up to about 60 km.

Average over high-quality profiles
The optimal horizontal extent of the regions to calculate a typical climatological mean from high-quality measurements is a trade-off between a sufficiently large number of profiles and atmospheric variability.Our experience of building atmospheric climatologies utilizing RO data (e.g., Foelsche et al., 2008;Scherllin-Pirscher et al., 2011a) showed that 10 • zonal bands were a reasonable choice for calculating mean atmospheric profiles from RO data.These bands range from 90 • S to 90 • N, resulting in 18 zonal bands.
The mean number of profiles per 10 • latitudinal band varied between 11 000 (June) and 15 000 (October).However, the latitudinal distribution of RO events is not uniform, which is due to the orbit characteristics of the GPS and F3C satellites (orbit inclinations of 55 • and 72 • , respectively).The largest number of RO measurements per latitude band is obtained at mid-latitudes, where more than 20 000 profiles enter BAROCLIM.A smaller number of measurements per latitude band is obtained at low (approximately 8000 profiles) and high latitudes (between 2000 and 4000 profiles).We judged that 8000 tropical profiles are sufficient to calculate mean profiles because atmospheric variability at high altitudes is small at low latitudes.We also considered 2000 high-latitude profiles to be sufficient due to the decrease of zonal surface area with increasing latitude.Thus, 2000 profiles are enough to represent a small area.
Multiyear monthly-mean bending angles for 10 • zonal bands for January and July are shown in Fig. 3 together with their standard error of the mean σ mean = σ/ √ N (σ is the standard deviation, and N is here the number of observations used to estimate the mean -not to be confused with refrac- tivity).Bending angles are negative (white areas) somewhat above an altitude of 80 km, and the standard error of the mean is larger than 10 % above 70 km in the winter hemisphere at high latitudes.While negative mean bending angles might be caused by residual systematic ionospheric errors (see Danzer et al., 2013), the high standard error of the mean (in a fractional sense) is a result of the decreasing bending angle with altitude.Below 60 to 70 km, however, mean bending angles are rather smooth and the standard error of the mean generally does not exceed 2 %.

BAROCLIM discrete model
Because of the generally decreasing bending angle with altitude, the mean bending angle (Fig. 3) is error-dominated above 80 km.Therefore we combined the mean bending angles with a priori information to generate a model that is useful also above the mesosphere.A priori information profiles can be obtained from already-existing climatological models or profile data sets.Current state-of-the-art analysis, reanaly- sis, or forecast products from NWP centers do not reach high enough in the atmosphere (the ECMWF model top, e.g., is at 0.01 hPa, corresponding approximately to 80 km) and can therefore not be used for extension to higher altitudes.
Since it is readily available, we decided to use the modified MSIS climatology as a priori information.In order to make maximum use of the information content of the RO data, and since the MSIS-90 climatology might be biased at high altitudes, we did not necessarily take the MSIS profile as representative for a given latitude band and month.Instead we performed a search in the spectral model of the MSIS bending angle climatology on a regular 5 • × 10 • latitude-longitude grid through all months from January to December, and we found the best match to the RO data using a least-squares fit to the RO mean bending angle profile in the altitude range from 60 to 80 km, where RO data quality is still high.To correct remaining background biases, the best-fitting MSIS profile was then multiplied with a fit factor obtained from regression with respect to the mean RO bending angle profile at high altitudes (least-squares adjustment from 60 to 80 km).We found fit coefficients being close to unity (0.99 to 1.01), with exceptions only in Southern Hemisphere winter at high latitudes (80 to 90 • S), where fit coefficients were as small as 0.96, 0.88, and 0.91 in May, June, and July, respectively.
To combine the mean RO bending angle profile with the corresponding MSIS profile, we applied statistical optimization by inverse covariance weighting (Gobiet and Kirchengast, 2004) between 60 and 80 km using an error correlation length of 2 km for the RO profile and an error correlation length of 15 km for the adjusted MSIS model profile.Furthermore, we assumed the MSIS background error to increase linearly from 0 % at 80 km to 15 % at 78 km, kept it constant at 15 % between 78 and 62 km, and then increased linearly again from 15 % at 62 km to 100 % at 60 km.All these percent values refer to the absolute values of the MSIS bending angle at the respective impact altitude level.While the linear increase of the relative error at the top and the bottom end of statistical optimization avoids too sharp transitions, the constant fraction of 15 % was determined empirically by Gobiet and Kirchengast (2004).
The observational error was set to the mean background error between 62 and 78 km and was constant with height (in absolute value, not percentagewise).Using these settings, we obtained smooth statistically optimized bending angles, for which the height where the retrieval to a priori error ratio (Gobiet et al., 2007) equals 50 % is 67.2 km for all profiles.Outside of the transition region from 60 to 80 km, the statistically optimized bending angle equals that of MSIS (above 80 km) or that of the mean RO profiles (below 60 km).
Even though Fig. 3 indicates that the mean bending angles reach down to the surface (2 km impact altitude approximately corresponds to 0 km altitude), mean bending angle quality also decreases in the lower troposphere.This is mainly because of the strongly decreasing number of measurements in the lower troposphere (see also Anthes et al., 2008) but is also due to reduced quality of the measurements.Besides this, averaging bending angles at low impact altitudes is tricky, since the lowest impact altitude in individual profiles (even if profiles are tracked all the way to the sur-face) depends on the bending angle.For the ray grazing the surface in one occultation event with a large bending angle, the impact altitude is larger than for the ray grazing the surface in another occultation event with smaller bending angle.It thus becomes dubious to talk about bending angle at the lowest impact altitude for the mean profile.
Being aware that MSIS is a dry-air climatology (no humidity is included in this model) and accepting that BAROCLIM will not reflect real atmosphere conditions at the lowest altitudes, we decided to use this model to extend BAROCLIM down to the surface.BAROCLIM is therefore, like MSIS, a dry-air model, being clearly wrong in regions were moisture is usually abundant, but for technical reasons smooth bending angles in the lower troposphere close to the surface are necessary when generating the BAROCLIM spectral model.
To extend mean RO bending angles down to the surface, we first extracted the MSIS profile for the given month and latitude and searched for the best fit in longitude.We then applied a gradual transition using a cosine weighting function from the mean RO bending angle α RO to the MSIS bending angle α MSIS .This weighting function was defined as w(z) = 1/2 1 + cos π(z top − z)/ z , and the tropospheric bending angle α trop was obtained from α trop = w(z)α RO (z)+ (1 − w(z))α MSIS .

B. Scherllin-Pirscher et al.: Generation of BAROCLIM and its use in RO retrievals
Since the amount of water vapor in the lower troposphere depends on latitude, we performed RO-MSIS transition between 5 and 10 km from 60 • S/N to 90 • S/N, between 7 and 12 km from 30 • S/N to 60 • S/N, and between 10 and 15 km from 30 • S to 30 • N.
To sum up, our BAROCLIM discrete model is available for every month (January to December) and has a horizontal resolution of 10 • zonal bands and a vertical gridding of 200 m.It relies 100 % on RO data from the upper troposphere up to 60 km.Above 80 km and below 5, 7, or 10 km (depending on latitude), it consists of data-driven adjusted MSIS profiles.

BAROCLIM spectral model
For fast and easy access to BAROCLIM at any latitude and impact altitude, and to make the functionality similar to the MSIS bending angle and refractivity spectral models in EGOPS and ROPP, we expanded the BAROCLIM discrete model in Chebychev polynomials and zonal harmonics.Since the bending angle scale height is more finely structured than the smooth, almost exponentially decreasing bending angle, we expanded a function into Chebychev polynomials, which depends on the bending angle scale height.
First we introduced the variable z = h − h surf (z ≥ 0), where h is impact altitude and h surf is the lowest possible impact altitude corresponding to a hypothetical ray grazing the surface.The lowest impact altitude was estimated from the MSIS climatology using h surf = N MSIS surf ×10 −6 R E , where R E = 6371 km is the mean radius of the Earth and N MSIS surf is MSIS refractivity at the surface extracted at the specific month and latitude and at longitude λ = 0 • .The bending angle α(z) was then extracted from the BAROCLIM discrete model by interpolation to a number of discrete impact heights evaluating z(x) = 100(ln 2 − ln(1 − x)) at k max values of x given by x k = cos(π(k − 1 2 )/k max ) (k = 1, . .., k max ).This mapping yields a finer vertical spacing at low altitudes and coarser vertical spacing at higher altitudes.Having α(z) at these discrete impact heights, the bending angle scale height H S (z) was calculated as where α surf is the bending angle at z = 0 (also extracted from the BAROCLIM discrete model).Chebychev coefficients, c j , were obtained from where and m and b are slope and intercept, respectively, of a straight line fit to the scale height at high altitudes.Finally, j = 1, . .., k max , and k max is the number of extracted Chebychev coefficients for k max −1 Chebychev polynomials (Press et al., 1986).The Chebychev coefficients were then expanded into zonal harmonics.Besides the Chebychev coefficients, also h surf , α surf , m, and b were expanded into zonal harmonics.
In general, zonal harmonics coefficients A n are obtained from a given function f (y), where normally y = cos θ (θ is co-latitude (polar distance)) as (see, e.g., Spiegel, 1979) where P n are Legendre polynomials, n = 1, . .., n max , and n max is the number of extracted zonal harmonics coefficients.
In our case f (y) was c j , h surf , α surf , m, or b.Thus, the final output of the BAROCLIM spectral model was n max zonal harmonics coefficients of k max Chebychev coefficients, surface impact altitude, surface bending angle, and slope and intercept of the straight line.
To reconstruct the bending angle from the BAROCLIM spectral model for a given impact altitude and latitude, we first applied Clenshaw's recurrence formula (Press et al., 1986) for zonal harmonics to obtain the Chebychev coefficients c j as well as h surf , α surf , m, and b.We also applied Clenshaw's recurrence formula to reconstruct G(x) (where x = 1 − 2 exp(−z/100)) before reconstructing bending angles using . (5) More details on the expansion of BAROCLIM into Chebychev polynomials and zonal harmonics as well as their reconstruction can be found in Scherllin-Pirscher (2013).
To settle on the order of the Chebychev polynomials and the degree of the zonal harmonics, we calculated differences between the bending angles from the BAROCLIM discrete model and the BAROCLIM spectral model for different choices of k max and n max , aiming at minimizing these differences.Since the BAROCLIM discrete model has a horizontal resolution of 10 • zonal bands (18 zonal bands), we found minimum differences for 18 zonal harmonics coefficients.Concerning the order of the Chebychev polynomials, we found reasonably good agreement between the BARO-CLIM discrete model and the BAROCLIM spectral model for 64 Chebychev coefficients.When using 128 Chebychev coefficients, the spectral model even reproduces the sharp tropical tropopause.For this reason we decided to use 128 Chebychev coefficients when constructing the spectral model and reconstructing the bending angle, but lower vertical resolution bending angles can be reconstructed using a smaller number of Chebychev coefficients.This could be useful for applications where computational speed is more important than high vertical resolution.The BAROCLIM spectral coefficients (stored in a NetCDF-file) as well as the Fortran-90 code needed to reconstruct bending angles from these coefficients are designed to be included in a future release of the ROPP software.ROPP is free of charge after registration at http://www.romsaf.org.
Figure 4 shows (as an example for the month of January) the BAROCLIM spectral model for 128 Chebychev coefficients and 18 zonal harmonics coefficients.The right panel of Fig. 4 shows that differences between the BAROCLIM discrete model and the BAROCLIM spectral model are, in general, within 0.5 % up to 60 km (a closer inspection of the differences reveals that it is even within 0.3 % in most places).Larger differences are found close to the 60 km altitude level (transition height of RO-only data and statistically optimized RO data) and above 80 km where the absolute amount of the bending angle is so small that even very small differences yield a noticeable percentage value.

Error sources
Atmospheric climatological fields of RO data are affected by (i) random statistical errors, (ii) systematic errors, and (iii) sampling errors (Scherllin-Pirscher et al., 2011a).Random statistical errors include, e.g., receiver thermal noise, clock stability/differencing errors, ionospheric noise, and statistical velocity errors (see e.g., Ramsauer and Kirchengast, 2001).Random statistical errors diminish by averaging over a large number of profiles.Since BAROCLIM is based on a very large number of quality-controlled RO soundings, all contributions from statistical errors are negligible, except at the highest altitudes (cf.Fig. 3).
Systematic errors are more important for BAROCLIM.From the RO measurement and retrieval perspective, these errors include systematic errors in orbit determination, local multipath, residual ionospheric errors, and errors due to as-sumptions in the RO retrieval.Systematic errors of BARO-CLIM also include contributions due to the additional use of MSIS at high and low altitudes.Schreiner et al. (2009) investigated the uncertainty of UCAR/CDAAC precise orbit determination (POD) of F3C satellites and found a velocity error of 0.17 mm s −1 , which approximately corresponds to an F3C bending angle error of 0.05 µrad (1 % at 60 km).Due to the lack of an alternative measurement system onboard F3C, Schreiner et al. (2009) could not give an estimate of the orbit accuracy.If we assume that half of the F3C velocity error is attributable to a systematic error component, the corresponding BAROCLIM error will be 0.5 % at 60 km.
Errors due to local multipath depend on the spacecraft size and on the reflection coefficient (Rocken et al., 2008).For F3C these local multipath errors are estimated to be smaller than 0.05 mm s −1 (Rocken et al., 2008), which corresponds to 0.015 µrad in bending angle when using velocity-errorto-bending-angle-error conversion given by Schreiner et al. (2009).This error corresponds to 0.3 % at 60 km.
Systematic residual ionospheric errors are important for BAROCLIM.In general, ionospheric residual errors depend on the level of ionization at high altitudes, which again depends on local time (i.e., day-versus nighttime conditions due to solar insolation) as well as on solar activity (Kursinski et al., 1997;Schreiner et al., 2011;Danzer et al., 2013).Danzer et al. (2013) showed that this error rarely exceeds 0.1 µrad from mid-2006 to the end of 2011.When averaging over this time period, it is even smaller than 0.05 µrad.Giving a conservative estimate and assuming a 0.1 µrad residual ionospheric error from mid-2006 to 2012 (solar activity increased in 2012), which corresponds to 2 % at 60 km, this error is the most important systematic error source of BARO-CLIM.
Systematic errors due to assumptions in the retrieval process (such as spherical symmetry) are assumed to be small at high altitudes (Rocken et al., 2008).Another BAROCLIM systematic error component results from the additional use of the MSIS model at low (tropospheric) and high (mesospheric and above) altitudes.Large systematic BAROCLIM errors in the troposphere are due to the absence of atmospheric water vapor in the MSIS model.For this reason BAROCLIM is not generally useful for tropospheric studies.Systematic errors from MSIS a priori information used at high altitudes (below 70 km) are assumed to be small due to the way MSIS was used (finding a bending angle profile that fits the mean RO data at high altitudes).
Finally, errors in BAROCLIM are caused by discrete sampling times and locations of RO measurements (sampling error; see, e.g., Foelsche et al., 2008;Scherllin-Pirscher et al., 2011a).The sampling error depends on the number of profiles and atmospheric variability and can be estimated from reference data that reflect true spatial and temporal variability.Using 6 to 7 years of RO data for BAROCLIM with more than 200 000 profiles per month (exception: June with 198 177 profiles), the BAROCLIM sampling error is negligible.However, since BAROCLIM is only based on measurements from 6 or 7 years, BAROCLIM might be biased relative to the long-term mean atmospheric state over 30 years.During the BAROCLIM time period, e.g., several major sudden stratospheric warming (SSW) events occurred in Northern Hemisphere winter (e.g., in January 2009 and 2010), yielding an RO climatology biased towards too high temperatures and too high atmospheric densities (i.e., too large bending angles) at northern high latitudes in these months.Other potential biases associated with the short RO record (e.g., the influence of QBO and ENSO) were estimated to be small.

Comparison to ECMWF
In December 2006 ECMWF started assimilating RO data in its operational assimilation system (Healy, 2007), which implies that ECMWF analyses and RO measurement data are not independent anymore.Healy et al. (2005) and Healy and Thépaut (2006) showed that assimilation of RO data significantly improved the forecast skill of the ECMWF operational system in the upper troposphere and lower stratosphere.Foelsche et al. (2009) found reduced biases in mean ECMWF analysis fields at least up to 30 km after ECMWF started assimilating RO.However, even though the assimilation is performed up to 50 km (von Engeln et al., 2009), a large RO observational error assumed in the assimilation at high altitudes limits the impact of RO data in this region.Systematic differences between the BAROCLIM spectral model and mean ECMWF analysis profiles therefore provide valuable information not only about the quality of BAROCLIM but also on the quality of ECMWF analyses, especially at high altitudes.
Figure 5 shows that the differences are small (< 0.5 %) in the upper troposphere and lower stratosphere region between B. Scherllin-Pirscher et al.: Generation of BAROCLIM and its use in RO retrievals 119 approximately 10 and 35 km.Besides this good agreement, striking features, which are very similar in all months, are (i) large negative differences in the troposphere, (ii) a band of positive differences (mostly < 1 %) above approximately 35 km, and (iii) large positive differences (> 2 %) above 50 km.
Large negative tropospheric differences are caused by BAROCLIM being a dry-air model.Neglecting atmospheric humidity yields unrealistically small bending angles in regions where humidity is high.Positive BAROCLIM minus ECMWF analysis differences above 35 km (< 1 %) and above 50 km (> 2 %) are mainly attributable to biases in ECMWF analyses rather than to BAROCLIM.In general, these biases above 35 km might be related to the bias correction of assimilated radiances from satellite measurements, whereas the specific bias above 50 km is most likely attributable to data from AMSU-A channel 14, which are assimilated without bias correction (S.Healy, ECMWF, personal communication, October 2014).
Similar ECMWF biases at high altitudes have been found with satellite measurements from the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) instrument on the European environmental satellite ENVISAT and from the Microwave Limb Sounder (MLS) instrument on the U.S. Aura satellite (Chauhan et al., 2009).Further comparisons of ECMWF analyses and temperature data from MIPAS and Sounding of the Atmosphere using Broadband Emission Radiometry (SABER, an instrument on the U.S. TIMED (Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics) satellite) also revealed a similar bias structure in ECMWF analyses above 40 km (J.Innerkofler, WEGC, personal communication, October 2014).All these comparisons consistently showed that ECMWF temperatures are too high at high altitudes, yielding a negative measurement minus ECMWF temperature bias.When pressure biases are small, a negative temperature bias corresponds to a positive refractivity bias, which in turn corresponds to a positive bending angle bias (see, e.g., Scherllin-Pirscher et al., 2011b) as shown in Fig. 5.
This comparison clearly shows that BAROCLIM is of very high quality at least up to 60 km and has the potential to validate middle-atmosphere data.

Use of BAROCLIM in RO profile retrievals
The intended aim of BAROCLIM was its use as a priori information in RO profile retrievals.We therefore evaluated its performance by processing occultation data from different RO missions and comparing retrieved atmospheric profiles obtained with different a priori information.
As mentioned in Sect.2, we used level 1a RO data provided by UCAR/CDAAC for all missions and applied the WEGC OPSv5.6 retrieval to obtain ionosphere-corrected bending angles.As bending angle a priori information for statistical optimization we used BAROCLIM and MSIS profiles co-located to RO events (termed "BAROCLIM-Col" and "MSIS-Col", respectively) and BAROCLIM and MSIS profiles that best fit the ionosphere-corrected bending angle (termed "BAROCLIM-SF" and "MSIS-SF", respectively, where SF means "searched" and "fit").This best-fit algorithm was similar to the "enhanced IGAM high-altitude retrieval scheme" described by Gobiet and Kirchengast (2004), searching for the best-fitting MSIS/BAROCLIM profile between 35 and 55 km and performing linear regression to find a multiplication factor to refine the fit to the data from 45 to 65 km.
With this approach we do not necessarily take a profile from MSIS/BAROCLIM corresponding to the latitude and season of the retrieval, but one that fits the data the best at high altitudes.Thus, with the SF approach we use MSIS/BAROCLIM as a library of different profiles representing different (average) atmospheric conditions on Earth.The approach should reduce sensitivity to biases in the climatology, although it does not guarantee that biases in the retrieved profiles are absent.
To assess the performance of BAROCLIM in RO profile retrievals we calculated monthly statistics of raw ionospherecorrected bending angles minus optimized bending angles.Since the purpose of statistical optimization is to reduce random errors, while preserving the mean, the mean difference between raw and optimized bending angles is an indicator of the quality of the background climatology.
Figure 6 shows mean results for January and July 2008 for CHAMP, GRACE-A, SAC-C, and F3C from 30 to 60 km impact altitude.Since noise of individual ionosphere-corrected raw bending angle profiles is high in the middle and upper stratosphere, difference profiles have a very large standard deviation.We did not include this information in the plot but note that it reaches approximately 15 % between 45 km (CHAMP) and 50 km (F3C).
Figure 6 shows that -while bending angle systematic differences of BAROCLIM-SF, BAROCLIM-Col, and OPSv5.6 are very close to 0 for all satellites and both months -MSIS-Col and MSIS-SF data are slightly negatively biased above 40 and 50 km, respectively.Largest differences (> 5 % above 55 km) can be found for MSIS-Col.
Figure 7 shows how differences propagate from statistically optimized bending angle to refractivity using MSIS-SF, BAROCLIM-SF, or OPSv5.6.Since bending angle and refractivity differences shown in Fig. 7 are calculated against ECMWF analyses, zero difference does not necessarily mean that it is close to the truth (cf.Fig. 5 and the discussion there).Bending angle and refractivity systematic differences are very similar for all satellites, with differences amongst the satellites being, in general, smaller than 1.5 % up to Contrary to systematic differences, the magnitude of the standard deviation features distinct satellite-dependent characteristics.Since CHAMP data noise is larger compared to the other satellites, OPS uses more weight of the a priori information, which results in smoother profiles and smaller standard deviation above 40 km.When comparing standard deviations of the three methods, larger standard deviations are found for the two search and fit algorithms than for OPSv5.6.
We conclude that the results using BAROCLIM seem promising, in particular when used in combination with the SF approach.As mentioned, such an approach should reduce the sensitivity to possible biases in BAROCLIM because it is then merely used as a library of different profiles representative of different (average) atmospheric conditions.The fact that BAROCLIM is based on data from only one mission (F3C) and from a limited period of time (2006 to 2012) is therefore not so important in this context; BAROCLIM can be used in this way for other RO missions in the past and in the future as long as the climate in the upper stratosphere does not change drastically in terms of global variations of bending angle.

Summary, conclusions, and outlook
In this study, we used radio occultation data from the F3C mission from August 2006 to December 2012 (more than 6 years of data) to compile a bending angle radio occultation climatology (BAROCLIM).After careful quality control we calculated multiyear monthly means for 10 • zonal bands from the surface up to 100 km.Since mean RO profiles become noisier above 60 km (in a fractional sense) and are error dominated above 80 km, we used a priori information from the MSIS climatology modified for RO purposes and applied statistical optimization from 60 to 80 km to extend BAROCLIM into the thermosphere.We also used the MSIS model to obtain smooth bending angles in the troposphere and performed a cosine transition to mean RO profiles in the middle troposphere.This implies that BAROCLIM is a dry-air model in the troposphere.However, BAROCLIM relies on RO data from the upper troposphere up to 60 km and also contains RO-derived information higher up.BAROCLIM spectral coefficients and the reconstruction code, which is needed to obtain bending angles, are designed to be included in a future release of the ROPP software.This RO package can be downloaded for free after registration at http://www.romsaf.org.
We showed that BAROCLIM is of very high quality in the stratosphere and lower mesosphere, where systematic biases are small.In this altitude range differences between BARO-CLIM and ECMWF analyses (forward-modeled to bending angle) rather show deficiencies in ECMWF analyses than in BAROCLIM.At 40 km, e.g., we find BAROCLIM minus ECMWF analysis bending angle differences of about 0.5 to 1.0 % and sometimes more.Above 50 km, this difference even exceeds 2 %.This is generally consistent with findings by Chauhan et al. (2009), using other types of satellite measurements.
BAROCLIM was originally developed to be used as a priori information for bending angle initialization in RO data processing.We thus evaluated BAROCLIM by comparing retrieved RO profiles initialized with different a priori information provided by BAROCLIM, MSIS, and ECMWF.These comparisons showed that RO bending angles initialized with BAROCLIM are close to raw (unoptimized) bending angles.This means that BAROCLIM-initialized bending angles preserve the mean of the raw measurements, while MSIS-initialized bending angles are slightly negatively biased.Comparison of retrieved RO profiles to ECMWF analyses also indicated that BAROCLIM outperforms MSIS.These results confirmed the capability of BAROCLIM to be used in RO profile retrievals.
The main advantage of BAROCLIM compared to the average bending angle approach proposed by Ao et al. (2012) and Gleisner and Healy (2013) is that utilization of BARO-CLIM yields individual RO profiles rather than climatological fields and individual RO profiles are known to provide accurate information of, e.g., tropopause characteristics, occurrence of multiple tropopauses, or the height of the atmospheric boundary layer.
Our current BAROCLIM spectral model only includes profiles of bending angle.An important BAROCLIM update could comprise its inversion to refractivity, density, pressure, and temperature so that these parameters could be used for other applications as well.

Figure 2 .
Figure 2. F3C bending angles (in µrad) in January in the latitude band from 80 to 90 • S before (top) and after (bottom) the application of the 4σ -outlier rejection criterion.The thick green line shows the mean of all profiles, the thick yellow lines 1 standard deviation, and the thin yellow lines 4 standard deviations.The dashed lines in the lower panel indicate the thin yellow lines in the upper panel.

Figure 3 .
Figure 3. Mean F3C bending angles (in µrad) (top) and their standard deviation (in percent) of the mean (bottom) in January (left) and July (right) as a function of latitude and impact altitude up to 100 km.

Figure 4 .
Figure 4. BAROCLIM spectral model (calculated with 128 Chebychev coefficients and 18 zonal harmonics coefficients) as a function of latitude and impact altitude for January (left) and difference between the BAROCLIM spectral model and the BAROCLIM discrete model (right).

Figure 5 .
Figure 5. Systematic difference between the BAROCLIM spectral model and ECMWF analyses forward-modeled to bending angle as a function of latitude and impact altitude up to 60 km for January (top left), April (top right), July (bottom left), and October (bottom right).

Figure 6 .
Figure 6.Global statistics of systematic difference between raw and statistically optimized RO bending angles for CHAMP, GRACE-A, SAC-C, and F3C for January 2008 and July 2008.Different colors indicate different a priori information used in the retrieval.

Figure 7 .
Figure 7. Global statistics of systematic difference (left) and standard deviation (right) between RO retrievals and ECMWF analyses for January 2008.The top panels show (statistically optimized) bending angle, and bottom panels refractivity.Different line types indicate retrievals for SAC-C, GRACE-A, CHAMP, and F3C.Different colors different indicate a priori information used in the retrieval.
ure 2 shows F3C profiles for the month of January before and after application of this 4σ criterion.This figure reveals that Figure 1.Flowchart of the BAROCLIM generation.

Table 1 .
Number of level 1a F3C data provided by UCAR/CDAAC (first column), number of bending angles profiles retrieved at WEGC (second column), and number of profiles that passed BAROCLIM quality control (last column).Numbers are shown for every month, adding up the data from August 2006 to December 2012.