Errors in GNSS radio occultation data : relevance of the measurement geometry and obliquity of profiles

Introduction Conclusions References Tables Figures


Introduction
Global Navigation Satellite System (GNSS) radio occultation (RO) measurements are increasingly used to validate other measurement techniques (e.g., Kuo et al., 2005), due to their high accuracy and precision (e.g., Schreiner et al., 2007;Steiner et al., 2009).For this purpose, it is important to be aware of the characteristics of RO measurements.A comparison of RO data with point measurements, like those from radiosondes, will always show some differences, due to the different nature of the measurement techniques.Since RO is a limb sounding technique, the data contain an inherent along-ray horizontal averaging over 200 km to 300 km (Kursinski et al., 1997).This gives rise to errors when retrieved profiles (derived under the assumption of spherical symmetry) are interpreted as being local vertical profiles.In the stratosphere and above these along-ray representativeness errors are small, but in the troposphere they are significant.Ahmad and Tyler (1999) analyzed these errors with an analytical approach, Healy (2001) with the aid of a simulation study.Later works by Sokolovskiy et al. (2005) and Syndergaard et al. (2005) developed so-called non-local observation operators and showed that these representativeness errors can be significantly reduced in future data assimilation systems by interpreting the retrieved profiles as complicated weighted averages of the three-dimensional atmospheric field.Since the contribution is focused around the tangent point of each ray (the point of closest approach), RO measurements can alternatively be regarded as profiles that follow the tangent point trajectory (TPT), but only in case of spherical symmetry this profile would exactly equal a profile of point measurements along the TPT.RO profiles are vertical in a sense that they contain height-resolved information about the state of the atmosphere, but they are far from being vertical in a strict geometric sense (cf. von Engeln et al., 2004).The departure from verticality depends on the measurement geometry, specifically on the orientation of the orbit planes of the transmitting and the receiving satellite.The work by Foelsche and Kirchengast (2004a, b) has indicated that horizontal variability errors (in all atmospheric parameters) increase with increasing obliquity of the RO profiles.We have built on this study and performed an endto-end simulation study, employing ray tracing through atmospheric fields in high-resolution, to determine the dependence of RO representativeness errors on the measurements geometry, and to investigate if the errors can be reduced, when validation is not performed with vertical reference profiles but with profiles that follow the TPT of the RO event.Section 2 describes the experimental setup, Sect. 3 the method to estimate the TPT from observed data.Results on the error analysis are presented in Sect.4, followed by conclusions in Sect. 5.

Experimental setup
The EGOPS 5.5 (End-to-end Generic Occultation Performance Simulator, Version 5.5) software tool (Fritzer et al., 2009) was used to generate simulated phase measurements and retrievals of the observables bending angle, radio refractivity, total air pressure, geopotential height, and (dry) temperature.Here we focus on the results on (dry) temperature, the parameter, which is most likely to be compared with measurements from other techniques.
"Dry temperature", T dry , means that temperature, T , is calculated from the observed refractivity with the assumption that water vapor is zero.At altitudes above 8 km (polar winter) and 14 km (tropics) the difference between T dry and T is always well below 0.1 K and T dry can be considered equivalent to T (Foelsche et al., 2008).At lower altitudes there is increasing water vapor influence, but T dry can always be calculated from reference data, if information on temperature and water vapor is available (e.g.Foelsche et al., 2008;Eqs. 2-3 therein).

Measurement geometry
We assumed a nominal constellation of 24 Global Positioning System (GPS) satellites as transmitters and a GRAS (GNSS Receiver for Atmospheric Sounding) sensor onboard the MetOp-A satellite with a nominal low Earth orbit (LEO) altitude of ∼820 km and an inclination of 98.7 • as RO receiver.With this constellation, more than 600 rising and setting occultations per day can be obtained (Luntama et al., 2008).We simulated measurements over a 24 h period on 22 June 2007, the date of the atmospheric fields used for the forward modeling (the specific day being arbitrary).Using this setup (summer in the Northern Hemisphere and winter in the Southern Hemisphere) a representative sample of meteorological situations can be expected.
We separated occultation events into seven 10 • azimuth sectors relative to the boresight direction of the receiving antenna aligned with the LEO orbit plane, to reflect the increasing deviation of the RO profiles from verticality with increasing angle-of-incidence.In contrast to the actual GRAS receiver on MetOp-A, which measures RO data up to 55 • azimuth, we used data up to an azimuth angle of 70 • , which has, e.g., also been used as azimuth boundary by the COSMIC Data Analysis and Archive Center (CDAAC, http://cdaac-www.cosmic.ucar.edu/cdaac/index.html).A schematic illustration of this division into azimuth angle sectors is given in Fig. 1, left panel: Sector 1, for example, with four symmetric sub-sectors (red) contains the RO events that are closest to the orbit plane of the receiving satellite.Table 1 summarizes the simulation design in terms of numbers of events per sector and the percentage distribution of events in (2x) 30 • latitude bands.We obtained a total of 636 occultation events during the selected 24 h period.This number is slightly smaller than the number Table 1.Distribution of the occultation events in the seven azimuth sectors and mean characteristics of the tangent point trajectories: azimuth angles of the sectors; number of profiles in each sector; percentage of profiles at low, mid, and high latitudes; mean elevation angle of the tangent point trajectories in four altitude levels; mean tangent point movement between 25 km and 2 km altitude; mean tangent point movement between 80 km and the minimum altitude; mean duration of the occultation events between 80 km and the minimum altitude; mean latitudinal tangent point movement between 80 km and the minimum altitude (positive for South-North movement).

Sector 1
Sector 2 Sector 3 Sector 4 Sector 5 Sector 6 Sector 7 Azimuth angle 0 of occultation events that can be obtained with the actual GRAS receiver (about 650, von Engeln et al., 2009), since the number of available GPS satellites is currently larger than the nominal constellation of 24.
The geographic distribution of the RO events is shown in the right panel of Fig. 1.Overall, there is a quite uniform distribution in latitude, but if we look at the distribution of RO events in individual sectors, we can see some interesting and characteristic features.From sectors 1 and 2 there are no RO events in the tropics, i.e., equatorward of 15 • latitude: The high inclination of the MetOp satellite (98.7 • ) together with the comparatively low inclination of the GPS satellites (55 • ) means that there are no GPS satellites close to the orbit plane in limb viewing geometry, when the MetOp satellite is near the equator.On the other hand, there are no RO events beyond 68 • and 64 • latitude in sector 4 and 5, respectively.More than two thirds of the RO profiles in sector 5 are concentrated in the latitude band between 30 • S and 30 • N, but almost half of the profiles in sector 1 are confined to latitudes beyond 60 • .This behavior inhibits a (meaningful) subsequent presentation of the results as function of sector and latitude so that we will focus to show the results as function of sector and the relation to latitude is to be kept in mind.

Forward modeling
High resolution (T799L91) analysis fields from the European Centre for Medium-Range Weather Forecasts (ECMWF) for 22 June 2007, were used to generate simulated quasi-realistic atmospheric phase delays.The horizontal resolution (T799) corresponds to 800 × 1600 points in latitude and longitude, respectively, and thus furnishes about twelve grid points within the typical horizontal resolution of an occultation event of ∼300 km (e.g., Kursinski et al., 1997).This dense sampling is useful and important to have a sufficient representation of horizontal variability errors in the occultation measurements.The 91 vertical levels (L91, hybrid pressure coordinates) extend from the surface to 0.01 hPa, being most closely spaced in the troposphere.Above the vertical domain of the ECMWF analysis field, the MSISE-90 climatological model (Hedin, 1991) was used.
As this study is focused on altitudes below 25 km, we made the assumption that ionospheric residual errors can be neglected (Kursinski et al., 1997).Forward modeling was thus employed without the ionosphere, resulting in considerable savings in computational time.We performed highprecision 3-D ray tracing with sub-millimeter accuracy (Syndergaard, 1999) and a sampling rate of 10 Hz, for forward modeling the events through the ECMWF analysis (refractivity) fields.The ray-tracing approach was employed, since it provides the actual TPT in three dimensions.Furthermore, this approach is consistent with the subsequently performed geometric optics retrieval.Since the ray tracing implemented in EGOPS 5.5 stops in case of superrefraction or multipath in the lower troposphere, when sharp vertical refractivity gradients are encountered (e.g., Sokolovskiy, 2003), the simulations do not account for these effects.
U. Foelsche et al.: Errors in GNSS radio occultation data

Geometry of RO profiles
Figure 2 shows the horizontal movement of the tangent points of the ray-traced RO profiles in the seven azimuth sectors (shifted by multiples of 50 km, to allow for a better visual separation of the results).RO events in sector 7, as most extreme example, start at 300 km distance at an altitude of 80 km, showing a tangent point movement of up to about 1000 km.Table 1 provides quantitative measures for the "non-verticality" of the RO profiles at different altitudes.The increase of the horizontal movement with increasing azimuth angle is obvious; it is caused by the GNSS and LEO satellite orbit planes being less and less co-planar.Furthermore, there is a clear increase of obliquity with decreasing altitude, since the ray bending results in a cumulative slowing of the downward movement of the tangent point, while the satellites move (and the Earth rotates) with unchanged velocity (resulting in the horizontal displacement of the tangent point).
RO profiles in sector 1 at altitudes between 70 km and 75 km (with negligible atmospheric influence) are still close to vertical, with an average elevation angle of 85.1 • .But this value reduces to 59.6 • at altitudes between 20 km and 25 km and even to 7.0 • near the ground (computed between 2 km and 3 km).For RO profiles in sector 7, the corresponding values are 15.0 • , 10.2 • , and 1.7 • , respectively.Elevations between 2 km and 3 km, and also the mean tangent point movements between 25 km and 2 km (last line of Table 1), are only computed for those profiles that actually reach an altitude of 2 km.This results in a tentatively stronger weighting of the more vertical profiles in each sector, and therefore a more representative estimate, since the profiles with extreme "horizontal smearing" (at low latitudes with strong ray bending and sharp vertical refractivity gradients) tend to stop at higher altitudes also in practice.
For context we also include the mean tangent point movement from 80 km altitude to the minimum altitude reached, which increases from 81.3 km in sector 1 to 555.7 km in sector 7. It is worth noting, that most of this lateral movement happens below 25 km altitude, where ray bending is strongest.Also the mean duration of the occultation event (between 80 km and the minimum altitude) increases clearly from 43.3 s in sector 1 to 99.6 s in sector 7. The South-North component of the TP movement can be expected to have the largest impact for climate applications.For individual profiles we can anticipate a systematic error, when, e.g., the upper part of the RO profile is further north, but the lower part further south than a vertical reference profile.When large ensembles of RO profiles are averaged for climate applications the latitudinal TP movements of the individual profiles will cancel to a high degree.This can already be seen in our ensemble of events, which has only been collected during one day: The South-North movement of the mean profile in sector 1 (last row of Table 1) is only 0.64 km.In sec-   tors 6 and 7, however, there are still uncompensated mean TP movements in South-North direction of 16.3 km and 53.2 km, respectively.

Observation system modeling and retrieval processing
Realistic observation system errors (including error sources like orbit uncertainties, receiver noise, local multipath errors and clock errors) have been superimposed on the simulated phase measurements.For the receiving system simulation, we used the specifications and error characteristics of the GRAS instrument (Luntama et al., 2008), similar to Steiner and Kirchengast (2005).The precise orbit determination (POD) error model contains satellite positioning and velocity errors (modeled randomly based on specified standard errors), where the along-ray velocity standard error of 0.05 mm s −1 , typical for modern POD systems, is the dominant error source.The radial position standard errors of the LEO and the GNSS satellites were set to 0.4 m and 0.2 m, respectively, a conservative bound for modern POD performance.Receiver noise was modeled as white Gaussian noise with a LEO antennae noise temperature of 150 K and a (single-side) tracking loop bandwidth of 10 Hz.
Regarding atmospheric profile retrieval, we applied a geometric optics bending angle retrieval scheme.The core of this algorithm, transforming phase delays to bending angles, is the algorithm described by Syndergaard (1999), which was enhanced to include inverse covariance weighted statistical optimization (with prior best-fit a priori profile search) as described by Gobiet and Kirchengast (2004).Since the forward modeling has been performed without an ionosphere, ionospheric correction was omitted.Refractivity profiles have been computed using an Abel transform and atmospheric parameters using a dry air retrieval, respectively, both steps employing algorithms of Syndergaard (1999).We did not undertake to separately analyze temperature and humidity as for the present analysis we preferred to inspect the parameter dry temperature, which does not require a priori information in the troposphere.

Reference profiles
All retrieved profiles have been differenced against corresponding "true" profiles, computed from the ECMWF analysis fields.We have used three types of "true" reference profiles: (1) at the mean RO event location, (2) along the "true" TPT (as computed by the ray tracing algorithm), (3) along the retrieved TPT.
We define the mean location of an RO profile as the latitude and longitude of the point, where the straight-line connection between the transmitting and the receiving satellite during the occultation event touches the Earth's ellipsoidal surface (corresponding to the tangent point location of real RO profiles between 10 km and 15 km altitude).Figure 2 indicates that this definition of the mean tangent point provides quite a good trade-off for a reasonable fit in the troposphere and stratosphere (in all sectors), especially in the upper troposphere and lower stratosphere, where the RO data quality is best.Since this definition is purely geometric the mean TP can be determined without an actual retrieval.
The exact TPT cannot be determined without performing ray tracing through the "true" atmosphere.A very good approximation of the TPT can be estimated from the GNSS and LEO satellite positions and retrieved impact parameter profiles, however, as described in the following section.

Estimation of the tangent point trajectory
Given the derived impact parameter, a, for a given ray path between the GNSS and the LEO satellite, a unit vector, rT , with origin at the center of refraction and pointing in the direction of the tangent point for that ray, can be estimated by considering the geometry in Fig. 3.Under the assumption of spherical symmetry it is clear that rT can be determined from the two vectors āL and āG (both having the length a) as (1) Each of āL and āG are obtained by considering the dashed triangles in Fig. 3, e.g., for the right side, where âG and rG denote unit vectors in the directions of āG and rG , respectively, n is the unit vector normal to the occultation plane (obtained from rG × rL ), and D G = r2 G − a 2 .Similarly we find The vectors āG = a âG and āL = a âL now readily follow from ( 2) and (3), respectively, and the unit vector rT is obtained from (1).The TPT is now found by repeating the calculations for each pair of satellite positions during the occultation and multiplying rT with the distance of the tangent point from the origin, r T .This, in turn, is found from Bouguer's law, a = n(a)r T , where n(a) is the refractive index for a given impact parameter, derived via the Abel transform.It is worth noting that the latitude of the tangent point, in practice obtained via the Cartesian coordinates of rT , is an estimate of the geographic latitude (not the geocentric one), since the calculations are done with respect to the center of refraction after a correction for the ellipsoid (Syndergaard, 1998).

Results
Based on the differences between retrieved dry temperature profiles and the different types of corresponding reference profiles (Sect.2.5) we have computed difference error statistics.Figure 4 shows the dry temperature errors, with the "true" vertical profiles at mean tangent point locations as reference.Between about 7 km and 25 km altitude, errors are small, reflecting the high accuracy and precision in this "core region" of the RO method (cf.Schreiner et al., 2007).Already in this representation a tendency of errors to increase with increasing angle of ray incidence (azimuth sector) is clearly visible.However, the bias (the mean of the difference profiles) does not exceed 0.5 K and standard deviations of the differences are smaller than 1 K, even in sector 7. Errors are particularly small in the altitude range between 10 km and    15 km, where the estimation of the mean tangent point location assures a good match, even with the vertical reference profile.
Below about 7 km altitude there is a clear increase in systematic and random errors, due to the growing influence of horizontal atmospheric variability.This confirms earlier studies with simulations based on atmospheric fields with lower horizontal resolution (Kursinski et al., 1997;Foelsche et al., 2004a, b).The decrease in the number of ensemble members with decreasing altitude (left side-panels in Fig. 4) is due to the fact that the ray tracer stops when severe superrefractive or multipath structures are encountered (e.g., Sokolovskiy, 2003;Beyerle et al., 2006).This happens frequently in the lower tropical troposphere, where pronounced humidity variations cause sharp vertical refractivity gradi-ents (e.g., Kursinski et al., 1997).Error statistics are only shown as long as there are at least ten ensemble members in each sector (corresponding to an altitude of about 1 km and above).
The setup of Fig. 5 is the same as in Fig. 4, but this time the "true" reference profiles have been extracted along the retrieved tangent point trajectories of the RO events.The overall behavior is very similar, but the magnitude of errors is markedly smaller (with very few exceptions).In order to allow for a better quantitative inspection of the errors in the different azimuth sectors, we computed representative mean values for the lower stratosphere (LS, calculated between 20 km and 25 km altitude), the upper troposphere (UT, 5 km-10 km) and the lower troposphere (LT, below 5 km). Figure 6 (left panels) shows mean dry temperature standard deviations for the three selected altitude intervals.In case of vertical reference profiles (red), there is a clear error increase with increasing azimuth angle (especially beyond sector 3), while standard deviations remain about constant up to sector 6, if the reference profiles are extracted along the "true" (blue) or along the retrieved (green) tangent point trajectories (0.3 K in the LS, 0.7 K in the UT, and about 2 K in the LT).If the vertical reference profile is used, the corresponding values in sector 6 are 0.6 K, 1.0 K, and 3.9 K, respectively (i.e., larger by up to a factor of 2).Furthermore, it can be seen that the retrieved tangent point is almost a perfect fit to the "true" one in the LS and UT, and a very good approximation in the LT.
It should also be noted that the comparatively small errors in the UT and LS (in case of vertical reference profiles) are partly due to the fact that the mean tangent point is estimated at altitudes between 10 km and 15 km.(Sect.2.5).Alternative definitions of the TP, e.g. at the minimum altitude reached by the event, would provide a better fit in the LT, but only for the high price of larger errors in the altitude range, where RO data achieve their highest quality and their largest impact in Numerical Weather Prediction (NWP).For this context we also show errors in refractivity (as relative quantities, due to exponential decrease of refractivity with altitude), since refractivities (and bending angles) are usually the RO parameters that are assimilated into NWP models (right panels of Fig. 6).The behavior is very similar, confirming  (Rieder and Kirchengast, 2001) and previous simulation studies (Foelsche and Kirchengast, 2004a, Fig. 4 and 7 therein), which showed that relative errors in refractivity closely mirrors those in dry temperature.
Figure 7 shows the mean absolute bias, since positive and negative deviations would partially cancel, when computing just the mean value.In the LS, the use of the different reference profiles results in negligible dry temperature differences (left panels) up to sector 3. Biases are generally small, but there is a marked increase with increasing azimuth angle: from 0.05 K to 0.35 K for vertical reference profiles) and from 0.05 K to 0.20 K for reference profiles along the 3-D trajectories.A candidate for such a behavior is an increasing misfit of the local radius of curvature, which is used in the ellipsoid correction, estimated at the mean tangent point location, and assumed to be constant for the RO event (cf.Syndergaard, 1998).
The UT bias is, except for sector 5, very similar for all types of reference profiles, since the estimated mean tangent point location is still quite a good approximation in this altitude interval.The mean absolute UT bias is about 0.1 K, and there is only a slight increase with increasing azimuth angle.There is no obvious reason for the deviation in sector 5 (it is, e.g., not caused by few profiles with large biases), but the respective panel of Fig. 4 shows that the small biases towards the lower troposphere set in slightly higher in the UT in this sector.The most remarkable property of sector 5 is that 67% of the RO profiles are concentrated between 30 • S and 30 • N.
The LT bias does not show a clear azimuth sector dependence if the reference profiles are taken along the "true" and retrieved TPT.Mean absolute bias values range between 0.2 K and 0.6 K.The interpretation is complicated by the complex behavior of the bias profiles in this altitude range (cf.Figs. 4 and 5).Up to sector 5 there is no big difference if the vertical reference profile is used, but the bias in sector 6 and 7 is significantly larger (0.8 K and 1.2 K, respectively).The LT is the only altitude interval, where reference profiles along the retrieved TPT show a noticeable difference from along the "true" TPT (with bias differences of up to 0.2 K), but their use still leads to considerably smaller errors than the use of vertical reference profiles in sectors 6 and 7. Relative refractivity errors (right panels on Fig. 7) are again very similar to those in dry temperature.

Summary and conclusions
We have performed an end-to-end simulation study for investigating (dry) temperature errors in GNSS radio occultation (RO) data in the troposphere and lower stratosphere.Thereby we focused on the dependence of systematic and random errors on the measurement geometry, specifically on the azimuth angle of the RO measurements with respect to the orbit plane of the receiving satellite.Furthermore, we determined whether the (apparent) errors are reduced when the reference profile is not taken vertical at the mean tangent point, but along the retrieved tangent point trajectory (TPT) of the RO profile.
The exact TPT can only be determined by performing ray tracing, but our results confirm that the estimated TPT -calculated from observed impact parameters and bending angles -is a very good approximation to the "true" one.Errors in retrieved atmospheric profiles decrease significantly, by up to a factor of 2, if the RO data are exploited along this retrieved TPT.The TPT of RO profiles should therefore be exploited whenever this is possible, especially when individual profiles are considered, e.g., in NWP applications.
Systematic and random errors in RO data tend to increase with increasing ray incidence angle (less pronounced, if the TPT is properly taken in to account), since the increasing obliquity of the RO profiles leads to an increasing sensitivity to departures from local spherical symmetry.Up to an azimuth angle of 30 • this effect is small, even if the RO profiles are assumed to be vertical.For applications requiring highest accuracy and precision it is advisable to exclude RO profiles with ray incidence angles beyond an azimuth of 50 • .

Figure 1 .Fig. 1 .
Figure 1.Ray incidence azimuth sectors (left) and geographic distribution of RO events in each sector (cylinder projection; right).Upright open triangles denote rising occultations while upside-down filled triangles denote setting occultations.

Figure 2 .
Figure 2. Horizontal movement of the tangent points with decreasing height for the seven azimuth sectors (shifted by multiples of 50 km for better separation). 17

Fig. 2 .
Fig. 2. Horizontal movement of the tangent points with decreasing height for the seven azimuth sectors (shifted by multiples of 50 km for better separation).

Figure 3 .Fig. 3 .
Figure 3. Illustration of the occultation geometry showing a ray path between the GNSS and LEO satellites with indication of the tangent point (TP).The ray asymptotes and a line from the origin through the TP are indicated by thin dotted lines.The two dash-lined right-angled triangles are equivalent (but mirrored and rotated) to the two right-angled triangles defined by the ray asymptotes and the satellite positions.The partial circle denotes a radius of unit length and origin at the center of refraction.
Figure 4. Difference error statistics for dry temperature in sectors 1 -7, with the "true" vertical profile at mean tangent point location as reference.

Fig. 4 .
Fig.4.Difference error statistics for dry temperature in sectors 1-7, with the "true" vertical profile at mean tangent point location as reference.

Figure 7 .Fig. 7 .
Figure 7. Mean absolute bias of dry temperature (left) and refractivity (right) in the lower stratosphere between 20 km and 25 km altitude (top), in the upper troposphere between 5 km and 10 km (middle), and in the lower troposphere below 5 km (bottom; note the different yaxis range).