Interactive comment on “ Processing and validation of refractivity from GRAS radio occultation data ” by K . B . Lauritsen

Scientific Quality: 2 The scientific quality of this manuscript is good. It provides a logical presentation of the research, it cites and discusses previous work. It follows the OLC approach by Gorbunov (2002), but it should provide more specific details on how this implementation is different from Gorbunov 2002 and whether these algorithm differences are significant. This work does not mention data gaps that are known to exist in GRAS BA data. The authors should discuss how they process through these gaps and what impact they may have on inversion errors.


Introduction
Radio occultation (RO) uses the radio signals continuously transmitted by the GPS satellites to measure the phase change as the radio signal path skirts the Earth's atmosphere on its way from the transmitting GPS to a receiver on another orbiting satellite.The phase measurements can be processed into vertical profiles of atmospheric parameters, such as refractivity, temperature, pressure and humidity (Kursinski et al., 1997;Healy and Eyre, 2000).RO is an additional application of GPS signals and RO instruments are currently in orbit on both operational and research satellites.Using Correspondence to: K. B. Lauritsen (kbl@dmi.dk)satellite observations more effectively will improve weather forecasting as well as climate change monitoring.RO data are becoming increasingly important for these applications (Anthes, 2011).
The radio occultation method should be regarded as complementary to passive atmospheric sounders; it has a high vertical resolution in atmospheric regions and it operates on completely different measurement principles (Collard and Healy, 2003).The fact that RO measurements do not merely reproduce other measurements is clearly shown within the field of Numerical Weather Prediction (NWP), where assimilation of RO data has a substantial positive impact (Healy et al., 2005;Healy and Thépaut, 2006;Cardinali, 2009).One of the key advantages of RO measurements is that they can be assimilated without bias correction.Therefore, they can potentially improve the assimilation of satellite radiance measurements by correcting model biases and providing "anchor points" to prevent adaptive, variational bias correction schemes drifting towards the NWP model climatology (see e.g.Dee, 2008).Within the field of climate monitoring (and for detection of climate change) the possibilities to accurately observe climate trends and to make bias corrections from independent measurements based on completely different measurement principles are very important, see e.g.Leroy et al. (2006).
The Global Navigation Satellite System Receiver for Atmospheric Sounding (GRAS) is an RO instrument onboard the Metop-A satellite, operated by the European Organisation for the Exploitation of Meteorological Satellites (EU-METSAT) (Luntama et al., 2008).GRAS has occultation antennas looking forward and backward relative to Metop's flight direction, and both are able to track two GPS satellites simultaneously.With the current GPS constellation of around 30 satellites, more than 650 occultations are tracked every day.The profiles are irregularly distributed across the globe, providing good overall spatial coverage.GRAS is capable of tracking in "phase-locked loop" (at 50 Hz) and in "raw sampling mode" (also called "open-loop"; at 1000 Hz) (Bonnedal et al., 2010).The resulting data profiles extend from the lowest part of the atmosphere up to about 80 km.
Metop-A was launched in October 2006 and is part of the EUMETSAT Polar System (EPS).The EPS is designed for operational data provision, which means that observations are rapidly made available to users.The near real time (NRT) timeliness requirement for level 1 data is availability to users within 2 h 15 min after sensing time.Each orbit of data (about 100 min) is down-linked over Svalbard, processed at EUMETSAT, and disseminated to users.Users include e.g.NWP centers worldwide.The data are also processed further into so-called level 2 geophysical and level 3 climate products at the Satellite Application Facilities (SAF), which are specialized development and processing centers in member states of EUMETSAT.The NRT timeliness on operational SAF products is 3 h.For more information and a general introduction of EPS, see Klaes et al. (2007).
Processing of operational GRAS/Metop bending angle (BA) data from EUMETSAT is done within the GRAS SAF at the Danish Meteorological Institute (DMI), with partners: European Centre for Medium-range Weather Forecasts (ECMWF), Institute d'Estudis Espacials de Catalunya (IEEC, Spain), and Met Office (UK).The objectives of the GRAS SAF are to deliver operational RO products from the GRAS instruments onboard the Metop satellites, and to supply the Radio Occultation Processing Package (ROPP) containing modules for pre-processing and assimilation of the RO data into NWP models (Offiler et al., 2008).The near real-time GRAS SAF data products consist of profiles of refractivity, temperature, pressure, and humidity, whereas off-line and reprocessing products also include bending angles and gridded climate data (Lauritsen et al., 2008).The GRAS SAF receives NRT level 1 bending angle data processed by EUMETSAT Central Applications Facility (CAF).These data are further processed to vertical profiles of refractivity (level 2) using state-of-the art inversion algorithms.The products are formatted as BUFR files and disseminated over the Global Telecommunication System network to NWP users worldwide within 1:41 h (average value), 1:48 h (90 % of the profiles), and close to 100 % of the profiles within the NRT timeliness of 3 h from observation time (see, e.g. the GRAS SAF website: http://www.grassaf.org).
In the present paper we focus on the validation of refractivity profiles derived from operational NRT bending angle data from GRAS/Metop.The remaining parts of the paper are organized as follows: in Sect. 2 we give a brief description of the data used in this study.In Sect. 3 we give a short overview of the processing and the GRAS RO data.Section 4 presents the analyses of GRAS refractivity data.Section 5 contains discussions and conclusions.

Description of GRAS data
The present study is based on the occultations observed by the GRAS instrument during the first 10 days of December 2010.The about 6000 occultations observed during this time period have a global distribution typical for observations made from a satellite in a Sun-synchronous, low-Earth orbit: approximately uniform in longitude but with a higher density (per area unit) toward the poles, and with a characteristic local-time structure.Figure 1 shows the geographical distribution of occultation events.
For a Sun-synchronous satellite like Metop-A, with a nodal precession rate of the satellite orbit that precisely matches the Earth's orbit around the Sun, the distribution of observations in space and local time is nearly stationary.The local times for equatorial crossings do not change with time.Figure 2 shows the distributions in longitude, latitude and local times for the studied period.The overall distribution is largely governed by the orbit of the Metop-A satellite whereas the detailed structure also depends on the GRAS instrument observational characteristics and on the transmitting GPS satellite orbits.The gaps around local times 09:30 and 21:30 -which is when the Sun-synchronous Metop-A satellite passes the equator -are due to the fact that the GRAS instrument observes the horizon at a distance of around 3000 kilometers.The GPS satellites are predominantly setting and rising in certain azimuths as viewed from Metop-A, and due to the instrument characteristics observations are done most efficiently in the forward and backward directions.The equatorial region cannot be observed from directly above the equator.The distribution of local times broadens toward higher latitudes, but the diurnal cycle is never fully sampled except near the poles which is relevant for the generation of climate data sets (Pirscher et al., 2007).

Processing of GRAS data
The processing to refractivity is based on operational bending angles processed by EUMETSAT CAF using the GRAS Processing Product Facility (PPF) version 2.16 from 1 July 2010.The bending angles consist of L1 and L2 retrieved from the observed L1 and L2 excess phases and a linear combination (LC) of L1 and L2 in order to reduce the influence of the ionosphere.The bending angles are based on geometric optics and accordingly the BA data below 10 km are not optimal for NWP and refractivity calculations (von Engeln et al., 2009).The operational input data may also contain gaps and missing values for bending angles and latitudes and longitudes.We do not process over data gaps but cut off the BA data at the first instance of a missing value for either the bending angle or the latitude and longitude coordinates.If there are more than one block of data we use the longest block.
For the processing to the refractivity in the upper stratosphere and mesosphere, the GRAS SAF operational system includes an approach referred to as optimal linear combination (OLC) of bending angles (Gorbunov, 2002).This approach gives the statistically optimized neutral atmosphere bending angle calculated from the observed L1 and L2 bending angles as well as a background profile based on a spectral representation of the MSIS90 climatological model (Hedin, 1991) transformed to bending angle space.Thus, it combines ionospheric correction and statistical optimization (SO) in a single least squares framework.The OLC method is an optimal filter based on the estimate of the variances and covariances of L1 and L2 signals and noises neglecting crosscorrelations between different impact parameters (Gorbunov, 2002).The background profile used in the SO for a given occultation is found through a global search in a small library of MSIS bending angles.By using a global search it is possible to find a background profile that fits the observations better than the local MSIS profile, while still being realistic in a climatological sense.In the searching process the MSIS bending angles are scaled and shifted (in bending angle logspace) in a least squares fit to the observed (non-optimized) LC bending angle between 40 and 60 km.However, the fit is not performed at altitudes where the MSIS bending angle deviates more than 30 % from the observed LC bending angle.More than 30 % deviation can be explained by ionospheric residual noise in the LC bending angle and stronger noise in the L2 channel.Currently about 6 % of all profiles are discarded because of such variations.The chosen background profile is the one where the scaling and shifting is the smallest in the sense described in more detail below.This way of choosing the background profile differs from the approach by Lohmann (2005) as well as that of Gobiet and Kirchengast (2004) who suggested to choose the model profile with the best least squares fit.In order to eliminate impact parameter ambiguities, a method consisting of finding the nearest monotonic impact parameter sequence with respect to the L2-norm is employed.
The main steps in the processing from GRAS bending angle to refractivity are summarized below: 1. Elimination of impact parameter ambiguities.

Linear interpolation to a fixed impact parameter grid
with steps of about 100 m.
3. Smoothing of bending angles using a sliding cubic polynomial fit with a window of about 1 km.-Identification of the model bending angle for which (lnα, β − 1) has the smallest L2-norm.
-Estimation of ionospheric signal and noise variance using the highest part (above 50 km) of the occultation.There is an upper height limit that excludes strong variations from the signal-band noise estimates.The limit is estimated dynamically and never exceeds 80 km.
-Calculation of relative mean deviation of neutral bending angle from the scaled and offset model bending angle using the data at heights 12-35 km (giving an estimate of the model variance).
-Solving a set of linear equations taking into account the estimated variances.
6. Inversion to refractivity via the Abel transform using piece-wise analytical integration and asymptotic correction (i.e.assuming exponential decrease of the bending angle profile above ∼100 km).
The refractivity profiles are quality-controlled and flagged as "bad" (and does not appear in the validation results shown in Sect.4) if one of the following is true: (i) refractivity profile does not reach below 20 km (4-5 % of profiles); (ii) one or more points in the refractivity profile below 35 km differ by more than 10 % from the corresponding profile obtained from ECMWF fields (1-1.5 % of profiles); (iii) refractivity profile reaches below the ECMWF model surface (0.8-1.0 % of profiles); (iv) refractivity is negative (0.05-0.1 % of profiles).These criteria flag about 7 % of all profiles (percentages for each criteria are shown in parentheses).

Results
In this section, the retrieved refractivity profiles are compared to the corresponding ECMWF profiles obtained by forward modelling analysis fields (analyses are available every six hour at 00:00, 06:00, 12:00, and 18:00 UTC; we use the analysis closest to the observation time) of temperature, pressure and humidity to refractivity as a function of altitude at the location of the occultations.The location is provided with the bending angle data, and is for each occultation given by the latitude and longitude of the point on the straight line between the occulting GPS satellite and the Metop satellite that grazes the Earth's surface.This corresponds to a tangent point altitude between 10 and 15 km (Foelsche et al., 2011).Statistics are separated into high latitudes (above 60 • N and below 60 • S), mid latitudes (30-60 • S and 30-60 • N), and low latitudes (between 30 • S and 30 • N).
Figure 3 shows the results for the first 10 days of December 2010 (cf.Fig. 1).Globally (Fig. 3a), the mean difference to the model profiles are within 0.2 % in most of the interval below 30 km (except near the surface).Above 30 km, there is an increasing positive bias relative to the model profiles, which mainly is believed to be due to a bias in the ECMWF model around 40 km (S.Healy, personal communication, 2009).For reference the figure also contains the statistical comparison to ECMWF short-range forecasts.It is observed that this yields slightly larger standard deviations above about 10 km and a slightly larger bias at about 30-40 km.Although the approach implemented for statistical optimization at the GRAS SAF seeks to avoid introducing biases (by finding a climatological profile fitting the data between 40 and 60 km), a bias may nevertheless be introduced at this level of processing, possibly because of limitations in the approach and limitations in the MSIS90 climatology.Below about 8 km, data are affected by atmospheric multipath, which results in a bias in retrieved refractivity because the operational GRAS data are processed using geometrical optics.The standard deviation is about 0.5 % between 8 and 25 km.The standard deviation increases significantly above 30 km mainly because of the exponential decrease of refractivity with altitude.The statistics at high (Fig. 3b), mid (Fig. 3c), and low (Fig. 3d) latitudes are somewhat similar to the global statistics, though with a significantly increasing bias in the lower troposphere as one goes towards lower latitudes.This again is due to increased multi-path in the moister tropical and extra-tropical regions, in combination with the geometrical optics processing.At low latitudes (Fig. 3d), the standard deviation is a minimum (about 0.25 %) in the upper troposphere and increases to a local maximum (about 0.7 %) near the tropopause.
The right panels in each of Fig. 3a-d, show the number of profiles as a function of altitude.These indicate the penetration depth, and show that less profiles reach into the lower troposphere at lower latitudes.In general, the cut-off in individual profiles is determined by the tracking of the L2 signal.Without the L2 signal, the standard ionospheric correction cannot be performed and therefore the bending angle is currently not processed for those altitudes where the L2 signal is missing.In particular for rising occultations, the L2 signal is not tracked in the lower troposphere, and profiles do not reach much below 10 km.However, a common approach  is to extrapolate the difference between the L1 and the L2 bending angle to lower altitudes (Kuo et al., 2004), and it is anticipated that the processing of bending angle at EUMET-SAT CAF will include such extrapolation in the near future (A. von Engeln, personal communication, 2010).This allows the generation of an extrapolated L2 signal, and thus, the ionospheric correction can be performed down to altitudes where the L1 CA signal is tracked in phase-locked loop mode.Figure 4 shows the penetration depth (the altitude of the lowest point in the profiles) for rising (left) and setting (right) occultations separately, and without the "L2 extrapolation" (top) and including "L2 extrapolation" (bottom).Thus, including the "L2 extrapolation", most profiles will reach the lower troposphere.The improvement in penetration depth is most noticeable for rising occultations, but also setting occultations will benefit from the "L2 extrapolation".The evaluation of the quality of the extrapolated profiles is currently ongoing in collaboration with EUMET-SAT CAF.

Conclusions
In this paper we have discussed operational bending angle data from December 2010 from the operational GRAS/Metop satellite, and shown statistical comparisons between retrieved refractivity generated at the GRAS SAF and ECMWF analyses.For global averages, the mean differences to ECMWF analyses are smaller than 0.2 % below 30 km, with standard deviations around 0.5 % for altitudes between 8 and 25 km.The penetration depth for rising occultations is generally poor, which is related to the tracking of the L2 signal.An extrapolation of the L2 bending angle to lower altitudes, currently under evaluation by the EUMET-SAT CAF and the GRAS SAF, is anticipated to improve the penetration in the near future.
Current RO bending angle GRAS data are obtained from closed loop sampling and derived from geometrics optics inversion (von Engeln et al., 2009).Therefore the data only have limited information about the atmosphere below about 8 km.Future GRAS data will be based on wave optics inversion and include both closed loop and open loop (raw sampling) data.This will significantly improve the profile penetration statistics especially in tropical regions (see, e.g.Gorbunov et al., 2011).

Fig. 1 .
Fig. 1.Geographical distribution of occultations observed by the GRAS instrument onboard Metop-A during the first 10 days of December 2010.

Fig. 1 .
Fig. 1.Geographical distribution of occultations observed by the GRAS instrument onboard Metop-A during the first 10 days of December 2010.

Fig. 2 .
Fig. 2. Latitude, longitude, and local time (LT) distribution of occultations observed by the GRAS instrument onboard Metop-A during December 2010.

Fig. 2 .
Fig. 2. Latitude, longitude, and local time (LT) distribution of occultations observed by the GRAS instrument onboard Metop-A during December 2010.

Fig. 3 .
Fig. 3. Statistical comparisons to ECMWF analyses (blue) and short-range forecasts (red) of refractivity derived from GRAS occultations for the first 10 days of December 2010.(a): Global, (b): high latitudes (above 60 • ), (c): mid latitudes (between 30 • and 60 • ), (d): low latitudes (below 30 •).The mean difference (solid) and standard deviation around the mean (dashed) are indicated in the panels to the left, whereas the number of profiles included in the statistics is indicated in the panels to the right.

Fig. 3 .Fig. 4 .
Fig. 3. Statistical comparisons to ECMWF analyses (blue) and short-range forecasts (red) of refractivity derived from GRAS occultations for the first 10 days of December 2010.(a) Global, (b) high latitudes (above 60 • ), (c) mid latitudes (between 30 • and 60 • ), (d) low latitudes (below 30 •).The mean difference (solid) and standard deviation around the mean (dashed) are indicated in the panels to the left, whereas the number of profiles included in the statistics is indicated in the panels to the right.
Searching through a bending angle model (for every month, every 10 • latitude, and every 20 • longitude) based on the MSIS90 climatological model.