## Assimilation of Radio Occultation Data in the Global Meteorological Model GME of the German Weather Service

D. Pingel and A. Rhodin

Abstract The assimilation of GPS radio occultations within the three-dimensional variational data assimilation system of the German Weather Service requires GPS radio occultation bending angle forward operators. To optimize the forward operator setup, different one- and three-dimensional bending angle forward operators are evaluated. The innovation statistics for radio occultation data from the CHAMP, GRACE-A, and FORMOSAT-3/COSMIC satellites are compared with estimates based on the background- and observation errors specified in the assimilation scheme. A numerical experiment is performed to assess the impact of assimilated radio occultation data on the weather forecast scores, to be compared to an experiment assimilating conventional observation data only.

### 1 Data Assimilation at the DWD

In the assimilation step of the numerical weather forecasting procedure, preceding the forecast run itself, the state of the meteorological model is updated using observational information. Conventional sources of meteorological observations are in-situ data such as radiosondes, ground station measurements, data from aircrafts, and buoys. They are complemented by satellite remote-sensing measurements, e.g., radiances of a given wavelength at nadir.

Global positioning system radio occultation (GPS RO) data is a relatively new and valuable source of observational information introducing additional information on temperature and humidity into the model (Healy et al., 2007b; Wickert et al., 2009). Benefits of RO data are good vertical resolution, independence of cloudy conditions (that generally afflict radiance measurements, cf. McNally (2002)), the lack of fundamental biases, and a nearly uniform global coverage.

At the German weather service (Deutscher Wetterdienst, DWD), the assimilation of observational data (including GPS radio occultation data) will be performed by

Deutscher Wetterdienst (DWD), Offenbach, Germany e-mail: [email protected]

the three-dimensional variational (3D-Var) assimilation method. The assimilation is done in a statistically optimal way, taking into account the observations y, the forecast xb (model background) from the previous analysis and error characteristics of both model and observations, given by their respective error covariance matrices B and R. Therefore, a thorough specification of observational and forecast errors is of high relevance in data assimilation. Within the assimilation process, a cost function J containing penalty terms for deviances of the control variables x (gridded model field) from the background xb on the one side, and of the observations y from their model equivalents H(x) on the other side is minimized:

J(x) = Jb + Jo = 2[x - xb]TB-1[x - xb] + 1[y - H(x)]TR-1[y - H(x)]. (1)

A forward operator H maps the model background data into observation space. For in situ measurements, H simply performs an interpolation to the location of the observation. In the case of remote sensing measurements, H is more complex and generally nonlinear. The analysis x is given by the zero of the gradient of the cost function J with respect to the control variables x, dJ = B-1 [x - xb] + HtR-1 [y - H(x)] (2)

d x with H being the Jacobi matrix of H. Tangent-linear and adjoint operators calculate the product of H with a vector to the right and to the left, respectively. The analysis increment is x - xb = [HtR-1H + B-1]-1HtR-1 [y - H(xb)] = BHt[HBHt + R]-1[y - H(xb)].

The solution of Eq. (3) is performed in observation space (Physical-space Statistical Assimilation System, PSAS). This implies a reduction of the size of the numerical equations to solve significantly compared to a minimization in model space. It is done by substituting x - xb = BHTz (4)

and solving

for z in observation space. The nonlinear problem is solved in an iterative procedure with a combination of an outer Newton algorithm and an inner preconditioned Conjugate Gradient (CG) algorithm. Nonlinearities are accounted for by iterating Eq. (5) with H linearized at the current estimate x;, denoted by H; (Daley and Barker, 2000):

Subsequently, the analysis increments of the model state variables are determined by post-multiplication of the solution z with the forecast error covariance according to Eq. (4).

Basically, phase delay, bending angle, or refractivity fields can play the role of the observations to be assimilated. At the DWD, bending angles have been chosen as observations, as they allow for a reliable specification of significant background and observational errors and can be derived without any assumptions regarding horizontal homogeneity or isotropy. Hence an appropriate forward operator H is to be implemented to map the model background at the location of the observation to the corresponding first guess bending angle.

### 2 Evaluation of Bending Angle Forward Operators

An essential step in the preparation of the assimilation of radio occultations is to find a reasonable setup for the forward operator that is suitable for data assimilation under the specific terms of operational application. One- and three-dimensional bending angle forward operators are evaluated, with the standard deviation of the innovations (differences of observation and first guess value) being considered as a benchmark criterion. While the three-dimensional ray-tracing operator shows the smallest standard deviation, the result of optimized one-dimensional forward operators come close to the performance of the three-dimensional operator.

The evaluation of a possible bending angle forward operator setup, providing a background profile of bending angles a(p) as function of impact height p, is a prerequisite for the assimilation of RO observations. A three-dimensional ray tracing model and a one-dimensional model have been considered to serve as forward operators:

The ray-tracing operator H3d (Gorbunov and Kornblueh, 2003) is supposed to be the best fit to the physical reality, as it takes into account the along-ray and transverse-ray horizontal gradients of the temperature and humidity atmospheric fields as well as the drift of the ray's tangential points in the course of an occultation. A certain disadvantage of the ray-tracing operator is the high demand of computing resources (time and memory).

The one-dimensional forward operator is based on the assumption of spherically symmetric atmospheric fields. It applies an inverse Abel transform

to the refractivity profile

(with c1 = 77.6 K hPa 1 and c2 = 3.73 ■ 105 K2 hPa 1) derived from the temperature T, pressure p, and water vapor partial pressure pv above a single point in the horizontal, the occultation point xoc. By reducing the horizontally extended occupation event to a single point xoc, the tangential point drift is neglected. Furthermore, the choice of xoc introduces a certain indefiniteness to the forward model. To model the drift of the tangential points at least to some extent, three versions of the one-dimensional forward operator have been implemented for test purposes, each of them considering a differently determined xoc: The first version H1d performs the inverse Abel transform at xoc as given in the corresponding observation data set of bending angles. For the second version H1dmn, xoc is determined as the mean of the individual ray's tangential points in the lowest 20 km of the occultation profile (indicated by the additional abbreviation "mn" in the acronym "1dmn"). The third version H2d applies the inverse Abel transform to each ray of the occultation profile separately, with xoc being the tangential point of the ray. This approach is expected to model the tangential point drift best compared to the two preceding versions.

The data set to test the performance of the different forward operators is a CHAMP (CHAllenging Minisatellite Payload) and GRACE-A (Gravity and Climate Recovery Experiment) phase delay data set provided by the GeoForschungsZen-trum (GFZ) Potsdam for January and July 2005, processed to bending angles at the DWD with the method of canonical transform CT2 (Gorbunov and Lauritsen, 2002, 2004a,b; Gorbunov et al., 2006). DWD global meteorological model GME (Majewski et al., 2002) three-hourly forecast fields served as model background. The statistical properties of the innovations, i.e., of the differences of observation and first guess, are then taken into account as a criterion to evaluate the different forward operators.

Figure 1a-c shows zonal averages of the ratio of the (O - B ) standard deviation of the three versions H1d, H1dmn, and H2d of the one-dimensional forward operator to the respective value of the three-dimensional ray-tracing operator H3d for January 2005 (bin size 1 kmx10°). For all three one-dimensional forward operators, the values of the ratios are larger than unity for most areas and heights (they are strictly larger than unity when considering global means). This indicates that the ray-tracing forward operator performs best in terms of standard deviation.

When comparing the different one-dimensional forward operators, the operator H1d, accepting the occultation point as given by the CT2, performed poorest (Fig. 1a). The corresponding standard deviation exceeds the respective value for H3d by up to 8%. However, the ratio is significantly smaller for the tropics compared to the extra tropics. This is related to the fact that due to differences of the atmospheric dynamics horizontal gradients of temperature and humidity are less pronounced in the tropics than in the extra tropics. As the provision for these gradients is one of the characteristics of the operator H3d compared to H1d, the bending angle values of the one- and three-dimensional operators tend to be comparable in the tropics.

Compared to the operator H1d, the operator H1dmn shows a reduced standard deviation for impact heights above ~8 km (Fig. 1b). This indicates that the averaging process brings xoc closer to the tangential points of the rays in the upper troposphere and lower stratosphere.

The first guess values of the version H2d, applying the inverse Abel transform to each individual ray, are closest to the ray-tracing results (Fig. 1c) as expected.

Fig. 1 Evaluation of bending angle forward operators: Ratio of the standard deviation of the one-dimensional forward operators (a) H1d (top),

(c) H2d (bottom) with respect to the ray tracing forward operator H3d, January 2005

st.dev. 1 dmn/3d 200501

st.dev. 1 dmn/3d 200501

st.dev. 2d/3d 200501

st.dev. 2d/3d 200501

Above an impact height of ~10 km, the difference to the operator H3d is less than 1% for most areas. Additionally, the ratio of the standard deviations is more uniform for H2d than for H1d and H1dmn.

However, for the lowermost rays especially at tropical latitudes, the more refined version H2d performs worse than the more simple versions H1d and H1dmn. This effect might be related to the fact that for most occultation events the positions of the tangent points of the lowest lying rays lie off the more or less straight line described by tangent points of higher rays. These "kinks" might be an inappropriate choice as a location for georeferencing points. A thorough investigation of this effect is necessary, though.

Below an impact height of km, the ray tracing forward operator H3d is significantly better in terms of standard deviation. This may be attributed to the fact that below this height the bending angle signal contains predominantly humidity information (especially in the tropics). Water vapor is, contrary to temperature, a quantity that fluctuates significantly on a horizontal scale comparable with the extent of the occultation. These fluctuations are considered for by the operator H3d, but not by the one-dimensional forward operators.

Figure 2 shows the same quantity as Fig. 1c, but for July 2005. The shift of the maxima of the ratio of the standard deviations in the lower troposphere towards the north is clearly visible, corresponding to the increase of humidity in these regions in the summer of the northern hemisphere.

The results of the evaluation of the bending angle forward operators are confirmed in another way: To expose the differences of the performance of the operators (in terms of standard deviations) more clearly, we consider the innovation statistics for those observations (y) whose respective first guess bending angle st.dev. 2d/3d 200507

st.dev. 2d/3d 200507

0.98 0.99 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08

Fig. 2 As Fig. 1c, but for the summer month

0.98 0.99 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08

Fig. 2 As Fig. 1c, but for the summer month values (H(xb)); differ by at least a given fraction of the ray-tracing first guess value ( Hsd(Xb));

|(Hf(Xb)); - (H3d(Xb)); | > e ■ (H3d(Xb))i (9)

with F e {1d, 1dmn, 2d} denoting the version of one-dimensional forward operator and components of H denoted by (H);. Each observation value (y) selected by Eq. (9) corresponds to background values evaluated with the different versions of the one-dimensional forward operator and with the 3d operator. To measure the closeness of the observations to the background value, we determined the innovation (y - H(xb)) statistics for each of the different combinations of observations y and background values H1d(xb), H1dmn(xb), H2d(xb), H3d(xb). Figure 3 shows the standard deviations (percentage of first guess value, global average, July 2005) for data sets selected this way with operators H1dmn and H3d, and e = 5%. The standard deviation of H3d is significantly smaller than the corresponding value of H1dmn, i.e., the first guess values derived with H3d tend to be closer to the observations than those calculated with H1dmn. Corresponding calculations with other combinations of forward operators and different values of e clearly indicate that the performance

Fig. 3 Evaluation of bending angle forward operators: Standard deviation of H1dml] and H3d (e = 5%) for observations selected according Eq. (9)

0.25

0.25

of the bending angle forward operator in terms of standard deviation increases in the order Hu ^ H^mn ^ Hzd ^ Hid.

Although there are significant differences on the level of innovation statistics, the overall deviances of the optimized one-dimensional forward operators from the ray tracing result are only minor (1-2%) for impact heights above 8 km. As the statistical weight of GPS RO observations is small below this height (large observation errors due to humidity gradients), this suggests the usage of an appropriately optimized one-dimensional bending angle forward operator for the operational assimilation, in agreement with results of other studies (Poli 2006; Healy et al., 2007a).

For the following monitoring and assimilation experiments, the optimized one-dimensional forward operator H1dmn is applied rather than the ray-tracing operator.

### 3 Monitoring

A long-term monitoring of the radio occultation observations is arranged. To this end, the statistical parameters of the differences of observations and background values (O - B) are determined on a monthly basis. Generally, the observed bending angles agree with the background equivalents with a deviation that corresponds to the assumptions made for the background and observation errors in the assimilation system.

The bending angle observations are compared with the corresponding background equivalents, derived from the forecast fields of the GME by application of the one-dimensional inverse Abel transform forward operator (version H1dmn, see Sect. 2). Near real time (NRT) bending angle data sets of CHAMP and GRACE-A from GFZ and FORMOSAT-3/COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) by CDAAC/UCAR (COSMIC Data Analysis and Archive Center/University Corporation for Atmospheric Research) are considered. In addition, a third data set is generated by applying the CT2 method (at the DWD) to near real time phase delay data sets of CHAMP and GRACE-A by GFZ. The CT2 method allows to estimate values of the bending angle observational error, which are not included in the near real time bending angle data set from GFZ. Thus, additional statistical information on the occultation data is gained.

A quality control is applied to the GPS RO data to be monitored: A 3a (O - B) check is applied to remove outliers. Observations with a relative observation error larger than 5% (with respect to observed bending angle value) are discarded, which removes observations mostly in the tropical lower troposphere (influence of humidity fluctuations) and in the stratosphere (residual stratospheric fluctuations induce additional noise). In addition, the value of the bending angle is confined to the interval [0, 0.02] to exclude errors due to ducting processes. These parameters of the quality control are also applied in the assimilation experiment (see Sect. 4).

The statistical properties of the corresponding differences, the innovations, are taken into account to assess the consistency of model and observational bending angle values and corresponding error estimates and are discussed in the following paragraph.

Figure 4a-c shows the global distribution of the number of RO profiles (bin size 1 km x 10°) available after quality control. The different orbit geometry of CHAMP/GRACE-A on the one hand and COSMIC on the other hand leads to the different patterns in the respective diagrams.

Figures 5, 6, 7, 8, and 9 show zonal averages (bin size 1 km x 10°) of statistical parameters of the (O — B) innovations of the three data sets for the time of March 2007, to be discussed in the following.

As the CHAMP and GRACE-A near real time data not yet include estimates for the observational error, a simple, partially linear functional dependence of the observation error on the impact height (Healy and Thepaut, 2006) is assumed for this data set (see Fig. 5a): The relative observation error decreases linearly from 10-1% for 0 km to 10 km impact height. Above 10 km, the observation error is assumed to be 1% of the bending angle value, until it reaches a lower absolute limit of 6 x 10—6 radians. This observation error model is in contrast to the zonally varying observation error for CHAMP and GRACE-A by the CT2 method (Fig. 5b) and COSMIC (Fig. 5c).

Stratospheric observation errors are to the largest part due to ionospheric fluctuations. For both the COSMIC and the CT2 derived CHAMP/GRACE-A bending angle data sets observation errors are estimated from differences between bending angles derived by geometrical optics and a climatology. For the tropospheric sections of the occultation profile, signal tracking in multipath areas (to be resolved by radio-holographic approaches) dominate the observation error. For COSMIC data, observation errors are derived by direct propagation of fluctuation in the excess phase signal to bending angles. For the CT2 method, the spectral width of the running spectra of the transformed signal wave field is used for the estimation. For details of the methods involved in the estimation of observation errors, see Gorbunov (2002); Gorbunov et al. (2006) (CHAMP/GRACE-A with CT2) and CDAAC/TACC (2007a,b) (FORMOSAT-3/COSMIC) and references therein.

The influence of tropical humidity on the observation error is consistent with higher values of the relative observation error in the lower latitudes. Above 20 km impact height, a difference in the observation errors of CHAMP and GRACE-A by CT2 on the one hand and COSMIC on the other hand is noticeable. It is due to the usage of a different climatology for estimation of the stratospheric component of the error. Comparing the near real time bending angles by GFZ with the bending angles obtained by applying CT2 at the DWD, it is apparent that the CT2 generally tends to keep more occultation data than the method of GFZ (cf. Fig. 4a,b in 10-20 km impact height), but assigns a relatively high observation error, especially in the lower and upper troposphere (see Fig. 5b). High values for the observation error cause more observations to be discarded at these height levels, see Fig. 4b. For all data sets, the estimated observational errors (Fig. 5a-c) are generally smaller than the model background errors (Fig. 6, model background errors are the same for the different RO data sets), especially in the upper troposphere of the extra tropics. The background error of the bending angle Bb is calculated from the background error covariances of temperature and relative humidity Bfspecified in the 3D-Var by application of the linearized observation operator:

Fig. 4 Number of GPS RO profiles (bin size 1 km x 1Ö0) for (a) CHAMP/GRACE-A (bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by CDAAC/UCAR)

Fig. 4 Number of GPS RO profiles (bin size 1 km x 1Ö0) for (a) CHAMP/GRACE-A (bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by CDAAC/UCAR)

Fig. 5 Estimated bending angle observation error for (a) CHAMP/GRACE-A (functional approach, bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by CDAAC/UCAR)

Fig. 5 Estimated bending angle observation error for (a) CHAMP/GRACE-A (functional approach, bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by CDAAC/UCAR)

Fig. 6 Estimated bending angle background error in the GME model (percentage of first guess value, derived from the background error covariances at t and rh by means of the tangent linear operator ffidmn)

Fig. 6 Estimated bending angle background error in the GME model (percentage of first guess value, derived from the background error covariances at t and rh by means of the tangent linear operator ffidmn)

The background error distribution of the bending angles are large in the extratropical jet stream regions where atmospheric variations are large. Another maximum is found in the lower tropical troposphere. Here humidity variations give the largest contribution to the bending angle background errors, reaching values of up to 15%. The observation errors estimated by the COSMIC and CT2 algorithm (Fig. 5) show maxima in the tropics as well. If no quality control was applied, these observational errors would be much larger and would reach values similar to the background errors. In the lower tropics only relatively few observations pass the criteria (cf. Fig. 4) so that the observational error is of the same order or less than the background error, which is a prerequisite for using the observations in the assimilation.

Figure 7a-c displays the standard deviation a of the innovation for the three data sets. The standard deviation has generally higher values in the tropics, due to the more pronounced influence of humidity. A maximum of the standard deviation at the impact height of ~18 km is caused probably by the presence of gravity waves. The deviations above 25 km result from ionospheric fluctuations which can only be partly corrected by taking into account the differences of the L1 and L2 phase. Remaining large differences of the retrieved bending angles from climatology are usually suppressed by fitting to the climatology. Reasonable bending angle values at higher altitudes are important if temperature profiles are derived from the bending angles by the Abel inversion. They are not essential if a forward operator is applied to the atmospheric profiles to derive model counterparts to the retrieved bending angles so long as the errors of the retrieval are specified correctly. This is done by the CT2 retrieval algorithm. In this case the model background will be used for

Fig. 7 Standard deviation of innovations (O — B)/B (percentage of first guess value, H1dmn bending angle forward operator) for (a) CHAMP/GRACE-A (bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by UCAR/CDAAC)

st.dev. (0-B)/B [%] GFZNRT all 200703 1dmn st.dev. (0-B)/B [%] GFZNRT all 200703 1dmn

Fig. 8 Ratio of standard deviation of innovations (O — B) with respect to estimated standard deviation for (a) CHAMP/GRACE-A (bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by UCAR/CDAAC)

D. Pingel and A. Rhodin (a) st.dev.(0-B)/est.err. [%] GFZNRT all 200703 1dmn

D. Pingel and A. Rhodin (a) st.dev.(0-B)/est.err. [%] GFZNRT all 200703 1dmn

Fig. 8 Ratio of standard deviation of innovations (O — B) with respect to estimated standard deviation for (a) CHAMP/GRACE-A (bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by UCAR/CDAAC)

(b) st.dev.(0-B)/est.err. [%] GFZNRTCT2 all 200703 1dmn

(b) st.dev.(0-B)/est.err. [%] GFZNRTCT2 all 200703 1dmn

(c) st.dev.(O-B)/est.err. [%] COSMIC all 200703 1dmn

(c) st.dev.(O-B)/est.err. [%] COSMIC all 200703 1dmn

Fig. 9 Bias of innovations (O — B)/B (percentage of first guess value, ff1dmn bending angle forward operator) for (a) CHAMP/GRACE-A (bending angles by GFZ), (b) CHAMP/GRACE-A (phase delays by GFZ, bending angles by CT2 at DWD), (c) COSMIC (bending angles by UCAR/CDAAC)

bias (0-B)/B [%] GFZNRTCT2 all 200703 1dmn bias (0-B)/B [%] GFZNRTCT2 all 200703 1dmn

bias (0-B)/B [%] COSMIC all 200703 1dmn bias (0-B)/B [%] COSMIC all 200703 1dmn

the regularization in the variational assimilation scheme. The innovation's standard deviation is compared to its estimate, the quadratic sum of the observation error and background error eest = (e2bs + ebg)1/2. Figure 8a-c shows the ratio of the standard deviation to the estimated deviation a/eesi. In general, the CHAMP and GRACE-A observed bending angles agree with the corresponding background equivalents reasonably well and within the assumptions made for the background and observation errors in the assimilation system. The ratio is less than unity for most areas and height levels, indicating a better agreement of observations and background values than estimated from the error specification. In the upper tropical troposphere, the relatively high value of the ratio a/eest again indicates the presence of gravity waves not accounted for in the background error model.

The bias of the innovations as depicted in Fig. 9a-c is reasonably small, but shows a remarkable vertically oscillating pattern in the upper troposphere and stratosphere in the tropics and mid-latitudes, which might be related to properties of the GME background model. In addition, the overall bias of the CT2 processed CHAMP and GRACE-A data (Fig. 9b) seems to be slightly more negative than the corresponding NRT result by GFZ and COSMIC (Fig. 9a,c).

### 4 Assimilation Experiments

The assimilation setup, including the optimized bending angle forward operator, is tested in assimilation experiments. The impact of the radio occultation observations on the weather forecast quality is assessed by an experiment assimilating RO data from the six FORMOSAT-3/COSMIC satellites in addition to conventional data only. A positive impact on the forecast scores in the southern hemisphere is observed.

A first assimilation experiment with RO observations has been performed, in order to estimate the maximal impact of the RO data and to apply the necessary optimization to the assimilation system. The assimilation experiment is carried out for the time of May 2007 by the DWD 3D-Var with a three-hourly assimilation time window. The assimilated observations are RO bending angle data from the six FORMOSAT-3/COSMIC satellites and conventional in-situ-data. In the corresponding control experiment, conventional data is the only observational source.

Figure 10a shows a significant reduction of the standard deviation of the differences of the experiment's analysis to the IFS analysis (ECMWF) for the geopoten-tial at 500 hPa in the southern hemisphere. For the northern hemisphere, the average change in the analyzed fields is neutral (not shown).

Subsequent forecast runs are performed to assess the impact on the forecast results. The forecast quality can be quantified by the anomaly correlation coefficient (ANOC) of an atmospheric field, measuring the correlation of deviances of forecasts with those of coinciding analysis with respect to the model climatology. Figure 10b shows the anomaly correlation coefficient of the geopotential field at 500 hPa in the southern hemisphere for the control experiment, the experiment

Deviation of 500hPa geopotential to IFS analysis control SH GPS RO SH

Deviation of 500hPa geopotential to IFS analysis control SH GPS RO SH

05/10 05/12 05/14 05/16 05/18 05/20 05/22 05/24 05/26 05/28 05/30 06/01 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00

ANOC geopotential 500 hPa SH 2007051400 till 2007053100 (18 forecasts)

05/10 05/12 05/14 05/16 05/18 05/20 05/22 05/24 05/26 05/28 05/30 06/01 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00

ANOC geopotential 500 hPa SH 2007051400 till 2007053100 (18 forecasts)

72 96

forecast time (h)

Fig. 10 Assimilation of GPS RO observation data in experiment, southern hemisphere: (a) standard deviation of analysis differences of geopotential at 500 hPa to IFS, (b) Anomaly correlation (ANOC) coefficient

72 96

forecast time (h)

Fig. 10 Assimilation of GPS RO observation data in experiment, southern hemisphere: (a) standard deviation of analysis differences of geopotential at 500 hPa to IFS, (b) Anomaly correlation (ANOC) coefficient

additionally including RO observations, and an experiment that assimilates conventional observations and AMSU-A data (using 18 forecasts). The assimilation of RO results in a significant increase of the ANOC in the southern hemisphere, with a similar result for the temperature field. It is remarkable that the improvement of the forecast scores due to RO is already half of the improvement to be expected when assimilating AMSU-A radiance observations in addition to conventional data. In the northern hemisphere, still a slight degradation of the ANOC is visible for geopotential and temperature, the reason for this is still to be identified.

The verification of temperature and humidity against corresponding measurements of radio soundings and aircraft data indicates a significant reduction of the standard deviation for the experiment with assimilation of RO observations in the southern hemisphere (not shown).

### 5 Conclusions

The performance of three different versions of a one-dimensional and a three-dimensional ray-tracing bending angle forward operator is evaluated for use in a three-dimensional variational data assimilation system. The ray-tracing operator takes horizontal gradients and the ray's tangential point drift into account. The different versions of the one-dimensional operator all apply the inverse Abel transform, assuming spherically symmetric atmospheric fields, but reflect the tangential point drift each to a different degree. Considering the standard deviation of the differences of first guess to observed occultation data above the impact height of 8 km, the three-dimensional and the optimized one-dimensional forward operator differ by less than 2%. Therefore, a choice of a properly optimized, numerically less expensive one-dimensional forward operator for the use in an operational data assimilation system seems legitimate. Monitoring of near real time radio occultation bending angle data sets from CHAMP, GRACE-A, and FORMOSAT-3/COSMIC has been carried out for the time of several months. The standard deviation of the observations-first guess differences are well within the error bound estimation of the assimilation system. An experiment assimilating bending angle observations shows a significant improvement of the forecast quality in the southern hemisphere when compared to a control experiment assimilating conventional data only, proving radio occultation data to be a valuable source of meteorological information.

Acknowledgements We thank the GeoForschungsZentrum (GFZ) Potsdam for reliably providing CHAMP and GRACE-A data sets of radio occultation data, both offline and near real time processed. The FORMOSAT-3/COSMIC project is gratefully mentioned for free provision of radio occultation data. We thank Michael Gorbunov for making the CT2 processing method and both bending angle forward operators available to us. The German Ministry for Education and Research supported the research project NRT-RO related to the near real time provision and usage of radio occultation data within the GEOTECHNOLOGIEN research program.

References

CDAAC/TACC (2007a) Atmospheric data inversion (ROAM): Algorithms for inverting radio occultation signals in the neutral atmosphere. http://cosmic-io.cosmic.ucar.edu/cdaac/doc/ documents/roam05.doc, http://tacc.cwb.gov.tw/cdaac/doc/documents/roam05.doc CDAAC/TACC (2007b) Ionospheric data inversion (GMRION): Algorithms for inverting radio occultation signals in the ionosphere. http://cosmic-io.cosmic.ucar.edu/cdaac/doc/ documents/gmrion.pdf, http://tacc.cwb.gov.tw/cdaac/doc/documents/gmrion.pdf Daley R, Barker E (2000) NAVDAS Source Book 2000: NRL Atmospheric Variational Data

Assimilation System. Tech. rep., Naval Research Laboratory, Monterey, Canada Gorbunov ME (2002) Ionospheric correction and statistical optimization of radio occultation data.

Radio Sci 37(5):1084, doi:10.1029/2000RS002370 Gorbunov ME, Kornblueh L (2003) Principles of variational assimilation of GNSS radio occultation data. Report 350, Max-Planck-Institut fur Meteorologie, Bundesstrasse 55, D-20146 Hamburg

Gorbunov ME, Lauritsen KB (2002) Canonical Transform Methods for Radio Occultation Data. Sci Rep 02-10, Danish Meteorological Institute Copenhagen, available at http://www.dmi.dk/ dmi/Sr02-10.pdf

Gorbunov ME, Lauritsen KB (2004a) Analysis of wave fields by Fourier integral operators and their application for radio occultations. Radio Sci 39(RS4010), doi:10.1029/2003RS002971 Gorbunov ME, Lauritsen KB (2004b) Canonical transform method for radio occultation data. In: Kirchengast G, Foelsche U, Steiner AK (eds) Occultations for Probing Atmosphere and Climate, Springer, Berlin Heidelberg, pp 61-68 Gorbunov ME, Lauritsen KB, Rhodin A, Tomassini M, Kornblueh L (2006) Radio holographic filtering, error estimation, and quality control of radio occultation data. J Geophys Res 111(D10105), doi:10.1029/2005JD006427 Healy SB, Thepaut JN (2006) Assimilation experiments with CHAMP GPS radio occultation measurements. Q J R Meteorol Soc 132:605-623, doi:10.1256/qj.04.182 Healy SB, Eyre JR, Hamrud M, Thepaut JN (2007a) Assimilating GPS radio occultation measurements with two-dimensional bending angle observation operators. Q J R Meteorol Soc 133 (626):1213-1227

Healy SB, Wickert J, Michalak G, Schmidt T, Beyerle G (2007b) Combined forecast impact of GRACE-A and CHAMP GPS radio occultation bending angle profiles. Atmos Sci Let 8(2): 43-50

Majewski D, Liermann D, Prohl P, Ritter B, Buchhold M, Hanisch T, Paul G, Wergen W (2002) The operational global icosahedral-hexagonal gridpoint model GME: Description and highresolution tests. Mon Wea Rev 130:319-338 McNally AP (2002) A note on the occurence of cloud in meteorologically sensitive areas and the implications for advanced infrared sounders. Q J R Meteorol Soc 128:2551-2556, doi:10.1256/qj.01.206

Poli P (2006) Effects of horizontal gradients on GPS radio occultation observation operators. II: A Fast Atmospheric Refractivity Gradient Operator (FARGO). Q J R Meteorol Soc 130(603):2807-2825, doi:10.1256/qj.03.229 Wickert J, Michalak G, Schmidt T, Beyerle G, Cheng C, Healy SB, Heise S, Huang CY, Jakowski N, Köhler W, Mayer C, Offilier D, Ozawa E, Pavelyev A, Rothacher M, Tapley B, Viehweg C (2009) GPS radio occultation: Results from CHAMP, GRACE and FORMOSAT-3/COSMIC. Terr Atmos Oceanic Sci 20(1):35-50, doi:10.3319/TAO.2007.12.26.01(F3C)

## Post a comment