Assimilation of GPS radio occultation data at DWD

We describe the status of the assimilation of bending angles from GPS radio occultations in the 3D-Var for DWD’s operational global forecast model GME (“Global Model for Europe”). Experiments show that the assimilation of GPSRO data leads to a significant reduction of biases in the analyses of temperature, humidity and wind in the upper troposphere and the stratosphere, as well as a better r. m. s. fit in the comparison to radiosondes. The impact on forecasts is most prominent in the data sparse Southern Hemisphere, but is also quite notable in the Northern Hemisphere extratropics. The positive results found in the impact experiments lead to the implementation of the assimilation of GPS radio occultations from GRACE-A, FORMOSAT-3/COSMIC and GRAS/MetOp-A into the operational suite on 3 August 2010. We also show some initial results from assimilation experiments using radio occultation data from the German research satellite TerraSAR-X.


Introduction
Data assimilation, i.e., the determination of the state of the atmosphere from observations plays a central role in numerical weather forecasting.Common observational data sources include in-situ data such as radiosondes, weather station measurements, and data from aircrafts and buoys.Nowadays, satellite remote sensing provides an increasing amount of information.Profiles of the thermal structure of the atmosphere can be retrieved from microwave radiance measurements, although with vertical resolution less than what is achievable by balloon ascents.
Global Positioning System (GPS) Radio Occultations (RO) (Kursinski et al., 2000) are a relatively new source of observational data with potentially high vertical Correspondence to: H. Anlauf (harald.anlauf@dwd.de)resolution (of the order of 100 m).They provide information on temperature and humidity from the atmospheric refraction of the GPS radio signals as measured by a satellite on a low-earth orbit.In contrast to radiance data GPS radio occultations are unaffected by clouds or precipitation and are bias-free (at least at the level of the raw observations).These properties make them interesting in particular for inter-instrument calibration studies, e.g. in climate research.Finally, the horizontal coverage of the radio occultation observations is close to uniform, which is beneficial for their use in data assimilation.Thus, GPSRO data have proven to be a valuable source of information in atmospheric research and weather prediction, and are now being used operationally at many NWP centers, e.g.MetOffice, ECMWF, NCEP, Météo-France, Environment Canada, and JMA (for an overview see the talks and proceedings of the GRAS SAF workshop (2008) and OPAC workshop (2010) and references cited therein).
In the past, the DWD was involved in two projects with focus on the use of GPSRO observations in data assimilation.Within the "Radio Occultation Processing Intercomparison Campaign" (ROPIC), the processing of initial phase delay observations to bending angles with different numerical algorithms, implemented at the participating centres, were compared.At DWD the canonical transform CT2 (see Gorbunov et al., 2006, and references cited therein) was evaluated.The results were validated with bending angles calculated from the ECMWF model background by a (presumably rather accurate) ray-tracing operator based on geometrical optics (Gorbunov and Kornblueh, 2003).
The second project, funded by the German Ministry for Education and Research, was carried out in collaboration with the GeoForschungsZentrum Potsdam (GFZ) and focused on fast processing and provision of GPSRO data and their use in numerical weather prediction (NWP).The GFZ provided GPSRO data of the CHAMP and GRACE research satellites and set up a near real-time processing chain to derive bending angles from phase delay measurements.At DWD, different bending angle forward operators were evaluated (3-D ray-tracing, 1-D Abel transform, multiple Abel transforms at tangential points of the respective rays) for the use in the 3D-Var data assimilation system of the DWD global model.First impact studies with assimilation of GP-SRO data at DWD (at that time on the IBM RS/6000 Power5 system running under AIX) have shown that these new observations were beneficial for NWP (see Pingel and Rhodin, 2009).
This paper reviews the developments at DWD that lead to the operational use of GPSRO observations in the global numerical weather prediction suite.

Assimilation of GPS radio occultations
At the German weather service, the operational assimilation of observational data into the global forecast model GME (Majewski et al., 2002) is performed by a three-dimensional variational assimilation system with a physical space analysis scheme (3D-Var PSAS) (Daley and Barker, 2000).
A variational data assimilation scheme determines the analysis by a combination of a preceeding short-range forecast (x b ) and the observations (y o ) with the help of a cost function J (x) that contains penalty terms for deviations of the analysis state (x) from both the forecast and the observations, weighted by their respective statistical uncertainties.
Here H denotes the non-linear observation operator, and R and B model the covariances of observation and background error, respectively.The 3D-Var PSAS performs the numerical minimization of the cost function in observation space.The assimilation of radio occultations is performed using bending angles as observations.The model counterpart is calculated from the refractivity field which depends on the atmospheric fields of pressure, temperature and humidity.The ray-tracing operator (Gorbunov and Kornblueh, 2003) leads to best results in terms of standard deviation of the difference of observation to background (Pingel and Rhodin, 2009) by a small margin, but it has a high computational cost.For the use of GPS radio occultations in operations we currently use a one-dimensional forward operator based on the Abel integral.To take into account the drift of the tangential points relative to the nominal occultation point at least to some extent, an effective occultation point is determined as the mean of the individual rays' tangential points in the lowest 15 km of the occultation profile. 1 The refractivity calcula-tion uses the expression given in (Rüeger, 2002) and recommended in GRAS SAF Report 05 (GRAS SAF project, 2011) for the neutral atmosphere (see also Healy, 2011 for a critical review and a discussion of non-ideal gas effects which we neglect throughout), N ≡ 10 6 (n − 1) = 77.6890p T − 6.3938 where T denotes the temperature in Kelvin, and p and p w denote pressure and partial water vapour pressure in hPa.
The bending angle, α, for a (locally) spherically symmetric atmosphere is a function of impact parameter, a, and can be calculated from the refractive index, n, where x = nr for vertical coordinate r.Above the model top (5 hPa, corresponding to approximately 36 km) the refractivity profile is extended using the MSIS-90 climatology.
The original implementation of Eq. ( 3) by Gorbunov (Gorbunov and Kornblueh, 2003) performs the integration using a spline interpolated version of the numerator to a grid regularly spaced in x with a rather small spacing (δx ∼50 m) even for high altitudes, which is inefficient.For the operational code which has to run on NEC SX-9 vector processors we performed several adaptations and optimizations, including a splitting-up of the integral into two ranges.The lower part still uses constant integration steps, whereas the upper part uses steps linearly increasing towards the upper integration limit.The numerical integration formula used is exact for a linearly varying numerator in the integrand of Eq. (3) for the lower part.For the upper part we assume an exponential variation of the numerator between levels similar to Healy and Thépaut (2006).In order to keep systematic differences between the original and modified implementation small in the troposphere and lower stratosphere, the height for the range splitting was set to 30 km.
An error model which accounts for observational errors of each occultation individually, based on the actual situation as provided with COSMIC data or derived by the canonical transform (Gorbunov et al., 2006) should be advantageous.However, in order to use the same model for all data we always applied the error model described in Healy and Thépaut (2006), where it is assumed that the combination of forward model error and observation error varies only with impact height, h, which is defined as h = a − R c , where R c is the radius of curvature provided with the measurement.The standard deviation of the error is assumed to decrease linearly from 10% of the observed value to 1 % for impact heights from 0 to 10 km, and to stay at 1 % above 10 km, with an absolute lower limit of 6 × 10 −6 radians.We do not use the sentative for the rays passing the troposphere where largest horizontal gradients are to be expected.observation error provided with the COSMIC data, except when it exceeds the value from the error model.Since we use bending angles for assimilation, we assume that observation errors are vertically uncorrelated.

Quality control of GPSRO observations
Before assimilation, several checks are applied to the profiles of bending angles: -As a first step we examine the product quality flag provided by the processing center.Profiles with nonnominal quality, non-nominal processing of excess phase or bending angle will only be monitored.
-For each ray in a profile, we check the consistency of provided tangent point and georeferencing point, by allowing for a maximum great circle distance of currently 3 • , corresponding to a maximum spatial separation of 330 km.In experiments using data before mid-2010 this lead to a rejection of a subset of rays within some COSMIC occultations passing over the 180 • meridian.These occultations typically contained several rays with no reported longitude, and some rays in the profile having geographical locations inconsistent with the occultation geometry.This finding was reported to the CDAAC team.The problem was not observed for COS-MIC near real-time data disseminated after 9 July 2010 and appears to have been fixed.
-We exclude occultation events with profiles starting above 20 km impact height.This mainly affects rising occultations featuring abnormally large first-guess departures.(See also Poli et al., 2008, who use a limit of 10 km.) -An observation-minus-first guess check of 4σ is applied to each ray to remove outliers.This check is chosen less tight than for conventional observations, as a compensation for a possibly poor representation of the first-guess uncertainty in the upper troposphere and in the stratosphere, as well as for model bias in the stratosphere.Note that the 3D-Var also employs a variational quality control during minimization that reduces the weight of observations with large departures from the first guess.
-Observations with a relative observation error larger than 10 % or an absolute observation error larger than 0.01 rad are discarded.This removes observations mostly in the tropical lower troposphere (influence of humidity fluctuations) and in the stratosphere (residual stratospheric fluctuations induce additional noise).
-The value of the bending angle is confined to the interval [0,0.02]rad to exclude errors due to excessive bending angles (caused by e.g.ducting processes).
-Adjustment of observation errors: To achieve a better balance in the assimilation process, individual observation errors are adjusted such that e o /e bg > 0.5.
-Clipping of the lowest section of GPSRO profiles: resolution of nonlinearities by radio-optical processing can cause a non-monotonous bending angle profile.Beneath a given impact height parameter (presently 8000 m) the GPSRO profile section below such points of non-monotonicity is discarded when the bending angle decreases by more than 3 times the assigned observational standard deviation.There is currently no check for the occurrence of super-refraction based on the refractivity lapse rate (Poli et al., 2009).
For the GPSRO observations passing quality control we perform a vertical thinning, so that no more than one ray will be retained per model level.For impact heights up to 10 km H. Anlauf et al.: Assimilation of GPS radio occultation data at DWD we thin to a vertical spacing of 400 m, between 10 km and 20 km to a vertical spacing of 500 m to 700 m, and between 20 km and 30 km to a vertical spacing of 2 km.The thinning process includes an exponential smoothing of the bending angle profiles in order to remove small-scale fluctuations and to improve the representativity of the observational data.For the assimilation we select impact heights in the range from 3 km to 30 km independent of latitude.

Monitoring of data and impact experiments
An initial monitoring of GPSRO measurements was performed at DWD with near real-time data sets of the CHAMP and GRACE-A missions provided by GFZ, and COSMIC by UCAR (University Corporation for Atmospheric Research) (Pingel and Rhodin, 2009).Impact experiments showed clear benefits of the use of these data.However, the computational cost of the initial implementation was very high, partly due to slow convergence of the 3D-Var when GPSRO were assimilated in addition to other observations, and partly due to the code for the bending angle forward operator needing optimization.The latter issue became even more important when DWD moved the operational system to the NEC SX-9, but it was eventually solved.
The investigation of the slow convergence turned up an inconsistency between the non-linear and tangentlinear/adjoint operator which was finally resolved.The inconsistency was discovered by comparing the gradient of the observational cost function for the contribution of a single ray evaluated in the tangent-linear approximation to finite differencing of the non-linear operator as a function of step size.Figure 1 shows the expected behaviour of the normalized difference of the gradients when the increment for finite differencing is chosen as a random vector and the scale of the increment is given by the background error.For large step sizes, the discrepancy is caused by non-linearity of the observation operator, while for very small step sizes rounding error dominates.The level of the plateau (10 −7 ...10 −8 ) shown in Fig. 1 for moderate step sizes is consistent with the expected accuracy for numerical differentiation with IEEE double precision.For the inaccurate implementation of the tangent-linear operator there was a broad plateau at the level of 10 −1 ...10 −2 .Furthermore, we found that convergence was much faster when surface pressure observations were removed from the assimilation.This problem was finally resolved by changing the Variational Quality Control (VQC) for surface pressure.VQC is implemented by use of a non-quadratic formulation (i.e. a non-Gaussian probability density function) for the observational cost function.In general a Huber-function (Huber, 1973) like expression is chosen which represents a Gaussian p. d. f. with fat tails. 2 This formulation was used for all observed quantities except surface pressure, because 2 Our formulation consists of a two times differentiable function which is very similar to the piecewise quadratic or linear Huber surface pressure observations had a large number of outliers which were not well represented by the fat tails of the Huber p. d. f. but better by a superposition of a Gaussian and a flat distribution.Therefore the latter formulation (see e.g.Andersson and Järvinen, 1999, and references cited therein) was previously used for surface pressure.The convergence of the minimization was considerably improved when the Huberfunction was applied to surface pressure as well. 3We explain the interference of surface pressure observations and radio occultations by the fact that RO basically observe refractivity as a function of geometrical height.Refractivity is closely related to atmospheric pressure and therefore RO provide information on pressure as a function of height which otherwise is only constrained by surface pressure observations.All other observation types only provide information on pressure gradients (temperature, thickness).
After the above issues had been addressed, we performed impact experiments for an extended period.These experiments exposed several problems in the Antarctic: when comparing short-range forecasts of both temperature and wind speed to radiosonde observations, we found an increase of the mean departures in the Antarctic stratosphere.In an experiment where radio occultations were blacklisted south of 65 • S, these biases clearly decreased to the level of the operational system.On the other hand, this increase in bias appeared to have no significant negative impact on other regions, so no blacklisting was done when the GPSRO were operationally implemented.Later investigations revealed that part of the increased bias resulted from an inappropriate choice of the vertical correlations in the background error covariance model which had been taken over from the previous forecast model resolution.
Furthermore, there is the well-known problem of a negative bias of the bending angles derived from GRAS closedloop data in the lower troposphere up to about 10 km.In a conservative approach one would blacklist this part of the profiles (Healy, 2008;Poli et al., 2008;Rennie, 2010).On the other hand, according to the error model described above, the assigned observation error for GPSRO observations also increases with penetration into the troposphere, leading to a reduced observation impact.
We performed a comparison of experiments which assimilated GRAS data down to 3 km and 8 km impact height, respectively.The negative bias in the bending angles appears to lead to a systematic but small reduction of moisture when comparing the respective analyses, but we did not see a significant detrimental effect in the comparison of short-range forecasts to radiosondes, nor did we see a significant effect function.The ensuing cost function has a positive definite Hessian and -in absence of strong non-linearities of the observation operator -no multiple minima, leading to enhanced convergence of the minimization algorithm.
3 The large number of outliers for surface pressure stems from persistent biases of some of these observations.This problem may be cured in the future by applying a bias correction. in the global precipitation.For this reason we currently apply the same selection and quality control criteria for GRAS as for other GPSRO, but we will reassess these parameters for GRAS when we find evidence for a negative impact. 4 4 As a side note, we would like to mention that we also assimilate other observations where the presence of a significant bias is established, but their omission would lead to worse forecasts.In the case of aircraft temperature, a bias correction scheme is under development.

Assimilation results
The main impact experiments were performed during the period from 19 May through 2 August 2010, using the operational model configuration with a horizontal resolution of 30 km, 60 vertical layers up to 5 hPa, and 3-hourly intermittent data assimilation.The control experiment consisted of the operational assimilation and forecast system, with GP-SRO measurements being passively monitored but not used.The GPSRO experiment employed exactly the same setup as the operational system except for the GPSRO measurements www.atmos-meas-tech.net/4/1105/2011/being actively used, with selection criteria as described above and applying the exact data cutoff time of the routine to the GPSRO observations as well.
The set of observations assimilated in the GPSRO experiment therefore consisted of in-situ observations from SYNOP stations, BUOYs, SHIPs, aircraft, radio-and dropsondes, PILOTs, and space-borne remote sensing data including sea surface winds from ASCAT, atmospheric motion vectors from geostationary and polar orbiting satellites, AMSU-A brightness temperatures from the MetOp-A and NOAA polar orbiting satellites, and GPSRO bending angles from GRACE-A, FORMOSAT-3/COSMIC 1-6, and GRAS on MetOp-A.
At DWD, the operational global data assimilation system uses a static bias correction scheme for the satellite radiances.The coefficients used in the bias correction are reevaluated off-line from time to time, to counteract drifts of satellite brightness temperature sensors and imperfect model climate.The bias correction is also susceptible to changes that influence the model state in the stratosphere.For the long-term impact experiment using GPS radio occultations, the coefficients of the bias correction therefore had to be adjusted independently from the operational system.
Figure 2 shows the comparison of the mean departures (observation minus first-guess) and RMS differences of radiosonde temperatures in three latitude bands (20 • N-60 • N, 20 • S-20 • N, and 60 • S-20 • S) for the reference run, which is the operational system, and the GPSRO experiment during July 2010.Evidently the bias and RMS error of the shortrange forecasts is reduced, with the impact being largest in the Southern Hemisphere (where less observational data is available than in the Northern Hemisphere) and in the strato-sphere.A comparison of the count of temperature observations used after variational quality control (displayed in the central columns in Fig. 2) shows that more data are actually used, which can be explained by the better consistency of the short-range forecasts of the GPSRO experiment with radiosonde measurements.
Since the refractivity of air depends on the water vapour pressure (cf.Eq. 2), we should expect some improvement in the comparison to radiosonde humidity observations.Disappointingly, this effect is rather small in the Northern Hemisphere and in the Tropics.We do find a significant positive effect in the Southern Hemisphere, see Fig. 3. Note, however, that it is restricted mostly to the upper troposphere where the humidity contribution to refractivity is expected to be rather small.We therefore think that this improvement is an indirect effect of the general improvement of the short-range forecast with respect to vertical structure rather than a direct effect of the analysis.
In order to assess the persistence of the information gained by the assimilation of radio occultations, we performed comparisons of forecasts to each experiment's own analyses.Figure 4 shows the average RMS difference of the 50 hPa temperature as a function of forecast time for the Northern Hemisphere (20 • N-90 • N), Tropics, and Southern Hemisphere (90 • S-20 • S) for July 2010.Evidently there is a large reduction in RMS error in the GPSRO experiment which significantly increases with forecast lead time.However, sometimes there is a minor degradation at shorter forecast lead times, as can be seen from the 12 h forecast time data points for the Northern and Southern Hemispheres.We find this behaviour for temperature forecast also at other stratospheric levels.As a possible explanation for this phenomenon we suggest that the analysis state is either not well balanced or not consistent with the forecast model's climate.
The impact on forecast quality may be demonstrated by showing the improvement of the anomaly correlation score of the 500 hPa geopotential.For the Northern Hemisphere this impact was rather small, but we found a very nice increase of this score for the Southern Hemisphere, cf.Fig. 5.
As the overall assessment of the impact of GPS radio occultations on forecasts was positive, it was decided to switch on their assimilation into the operational system on 3 August 2010.

Radio occultation data from TerraSAR-X
The German research satellite mission TerraSAR-X carries as a secondary payload an IGOR GPS receiver (Beyerle et al., 2010), which is of the same type as those installed on the COSMIC micro-satellites.TerraSAR-X is capable of providing roughly 200 setting occultations per day, which are processed by GFZ Potsdam and disseminated via GTS through DWD since 21 October 2010.
We performed an experiment that monitored and assimilated the radio occultation data from TerraSAR-X in addition Mean-Sea-Level Pressure, NH GME Routine TSR-X Experiment Fig. 7.The RMS difference of mean-sea-level pressure as a function of forecast time in the Northern Hemisphere for the operational routine (dashed blue line) and the experiment including GPSRO data from TerraSAR-X (thick red line).The statistics is based on 56 forecasts (00Z and 12Z) corresponding to 28 days of data (22 October-18 November 2010).The error bars give a rough indication of the uncertainty of the mean as estimated from the time series of the RMS differences.The comparisons are performed against each experiment's own analyses.
to those used in the operational suite.Figure 6 compares the global mean and standard deviation of the first-guess departures and of the analysis departures of the bending angles from TerraSAR-X and COSMIC FM-1 for the period from 22 October to 18 November 2010.Taking into account that the TerraSAR-X occultations generally have a smaller penetration into the lower troposphere than COSMIC, we con-clude that the quality of the TerraSAR-X data is quite similar to that of COSMIC.The comparison of forecasts including these additional data to the operational forecasts showed either marginal (see Fig. 7 for the RMS differences of surface pressure forecast in the Northern Hemisphere for the first 96 h) or no significant impact.This was to be expected, as TerraSAR-X amounts to roughly 10% of the currently available radio occultation observations.Nevertheless, since we found no negative impact, and since the number of occultations observed by the COSMIC mission is expected to decrease due to its end-of-life approaching, it was decided to switch on operational assimilation of TerraSAR-X data into the 3D-Var on 9 December 2010.At the same time, modified vertical background error correlations were activated to mitigate the problems encountered in the Antarctic stratosphere (see Sect. 2.2).

Conclusions
The above results demonstrate that GPS radio occultations provide valuable information for numerical weather prediction, and that their assimilation leads to improved atmospheric analyses and forecasts, confirming the findings of other operational weather services.Being a source of information independent of and complementary to previously used data, they may help to detect problems and identify deficits in the data assimilation system.
Experiments using radio occultation data from GRACE-A, COSMIC, and GRAS showed a significant positive impact on DWD's global forecast system, so that assimilation of these data was turned on 3 August 2010.Evaluation of the quality of GPSRO data from TerraSAR-X data was positive, although the additional impact was very small.In light of the already diminishing number of data from the COSMIC mission, and in order to be able to gain the most benefit from this kind of observational data, the assimilation of TerraSAR-X at DWD was activated on 9 December 2010.

Fig. 1 .
Fig. 1.Difference of the gradient of the observational cost function for the contribution of a single ray evaluated in the tangent-linear approximation to finite differencing of the non-linear operator as a function of step size.

Fig. 2 .
Fig. 2. The mean departure (left) and RMS difference (right) of radiosonde temperature observations to 3-h forecasts for the reference run (dash-dotted blue line) and the GPSRO experiment (continuous red line).The rows show the comparisons for observations in the latitude bands 20 • N-60 • N, 20 • S-20 • N, and 60 • S-20 • S (from top to bottom).The central columns display the number of samples used in the experiment (red) and the difference of experiment to reference (black).

Fig. 3 .Fig. 4 .
Fig. 3.The mean departure (left) and RMS difference (right) of radiosonde humidity observations to 3-h forecasts for the reference run (dash-dotted blue line) and the GPSRO experiment (continuous red line) in the Southern Hemisphere (60 • S-20 • S).The central columns refer to the number of used observations, cf.caption of Fig. 2.

Fig. 5 .
Fig. 5. Anomaly correlation of the 500 hPa geopotential height in the Southern Hemisphere (90 • S-20 • S) of the reference (GME routine, blue line) and of the GPSRO experiment (red line), based on the statistics from 31 forecasts during July 2010.The shaded blue area and the red error bars indicate the respective 90 % confidence intervals derived from the bootstrap method.

Fig. 6 .
Fig.6.The global mean (left) and standard deviation (right) of observation minus first-guess (continuous lines) and observation minus analysis (dashed lines) of the TerraSAR-X (blue) bending angle departures for impact heights from 3 km to 30 km, and compared to COSMIC FM-1 (red).The departures are normalized with respect to the observed bending angles.The statistics is based on 28 days of data (22 October to 18 November 2010).The central columns display the number of samples passing first-guess checks and being used after the variational quality control within a 1-km interval (left: TerraSAR-X, right: COSMIC).Impact heights between 21-22 km, 23-24 km, and 25-26 km were sparsely populated after vertical thinning.