Assimilation of GPS radio occultation data for Numerical Weather

With the availability of approximately 4,000 radio occultation soundings per day within three hours of observation, COSMIC has the potential to contribute significantly to global and regional weather analysis and prediction. However, the basic radio occultation measurements (phase delays) are very different from traditional meteorological measurements (i.e., temperature, water vapor), and to effectively assimilate them into weather prediction models is a challenging task. Over the past five years, considerable progress has been made in the development of an effective strategy for the assimilation of GPS radio occultation data. In this paper, we (1) review the measurement and data reduction procedures, (2) discuss the error characteristics of the GPS radio occultation data, (3) discuss the various strategies for data assimilation, (4) review results from recent data assimilation research, and (5) provide suggestions for future research. Results from recent studies have led to the conclusion that the best strategy to assimilate GPS radio occultation data is a mixture of bending angles below 10 km and refractivity above 10 km using a variational approach. The assimilation of GPS radio occultation data is likely to have a significant positive impact on global and regional weather prediction through improved definition of water vapor, temperature and wind fields. Although refractivity and bending angles are not directly related to the winds, the assimilation of GPS data leads to improvements in the wind analysis through internal model dynamic adjustments. In order to make optimal use of GPS radio occultation data in weather analysis and prediction, considerable research is needed in: (l) better characterization of GPS measurement errors, particularly in the lower troposphere, (2) improving the computational efficiency and optimizing the strategy of bending angle and refractivity assimilation, and (3) performing a set of observing system simulation experiments with realistic simulation of GPS radio occultation data. These research tasks should be conducted prior to the launch of the COSMIC satellites, so that we can fully realize the potential of COSMIC data in global and regional weather prediction.


INTRODUCTION
The lack of data over the oceans and other remote regions contributes greatly to the uncertainties in the initial state of global weather-prediction models, which in turn limits their forecast skill (Rabier et al. 1996).Nutter et al. (1998) calculated the standard deviation of the 500-mb geopotential height difference between NCEP (National Centers for Environmental Prediction) and ECMWF (European Centre for Medium Range Forecast) global analyses for five winters.They found that the major differences between these two global analyses occur primarily over the Pacific Ocean, Atlantic Ocean, and the Polar regions in the northern hemisphere (Fig. 1), where radiosonde observations are sparse (Fig. 2b).Calculations made using more recent data from January 6 through February 27, 1997 yielded essentially the same results (not shown).
The first experimental evidence that radio occultation data may have a positive impact on weather prediction systems was obtained from statistical comparisons of the GPS/MET retrieved refractivity and temperature profiles to global Numerical Weather Prediction (NWP) model analysis (Rocken, et al., 1997).It appeared, that the agreement of the GPS/MET data with NWP was noticeably better over data-dense (U.S., Europe) versus data-sparse (Pacific Ocean) regions.Since there are no physical reasons to expect different errors of the GPS/MET data at those regions, this difference must be attributed to the quality of the NWP analysis which obviously has larger errors over data sparse regions.Similarly, Leroy (1997) compared time-averaged analyses of geopotential heights computed from GPS/MET (Ware et al. 1996) data with operational analyses over the same period from ECMWF and found large differences over datasparse regions such as the South Pacific Ocean.These studies suggest that the GPS radio occultation data is likely to have a significant positive impact on global analyses and global weather prediction (Anthes et al. 2000).
Given the fact that microwave sounding data from geostationary and polar orbiting satellites are routinely incorporated into global analyses, one might ask why there should be significant differences in analyses over remote areas.Differences in the global analyses can be attributed to a number of factors, including differences in physical parameterizations and numerics of the assimilation models and details of the assimilation procedures.The fact that the differences between NCEP and ECMWF analyses are small over areas rich in radiosonde data suggests that these differences are most likely data related.A possible explanation is that the current satellite observations over the ocean have insufficient accuracy and/or vertical resolution and therefore do not provide sufficiently strong constraints to global data assimilation systems, allowing the systems to have a higher degree of freedom in analyzing the data.Consequently, different data assimilation systems produce different analyses over data sparse regions.1985-1986 and 1989-1990.Contours are every 5 m, with values exceeding 15 m shaded (from Nutter et al. 1998).
With a total of 4,000 soundings per day, uniformly distributed around the globe, COSMIC will provide independent and accurate data over the Earth.Figure 2a shows the projected COSMIC sounding distribution over a 24-h period.The radio occultation soundings from COSMIC represent a factor of 3~4 increase in the number of upper-air observations over the existing global radiosonde network (which has a total of ~900 stations, Fig. 2b).The high vertical resolution of radio occultation soundings provides important information on the vertical structure of temperature and pressure fields in the stratosphere and the troposphere.Through continuous assimilation, COSMIC radio occultation data will transfer information to the wind fields, as shown by Kuo et al. (1997) and Leroy (1997).With improved descriptions of temperature, pressure and wind fields, important synoptic-scale atmospheric circulation systems (which are the drivers of surface cyclones and fronts, and their related weather) will be better described in the model initial conditions over the oceans and other data sparse regions.As reviewed by Anthes et al. (2000), GPS/MET radio occultation data have the following characteristics: (1) no instrument drift, (2) minimally affected by aerosols, clouds and precipitation, (3) high vertical resolution, and (4) relatively uniform global coverage.The COSMIC system with a constellation of eight low-Earth orbiting (LEO) satellites will have a global "refresh rate" of about 100 min.This means that the same region of the Earth will be sampled once every 100 min.Moreover, COSMIC data are expected to be available to the operational centers within three hours of observations.These features provide some significant advantages over radiosondes and space-based sounders.When used together with these complementary systems, COSMIC data offer an opportunity to significantly improve the skill of weather prediction models.However, despite these advantages, the use of GPS radio occultation data in a global analysis is not a trivial matter.First of all, the raw measurements of GPS radio limb soundings are the phase delays and amplitudes of the GPS radio signals.It takes a number of steps to reduce the data to the traditional meteorological variables of temperature, pressure and water vapor.Because of the ray traversing geometry, the GPS radio occultation data have unique characteristics, which is very different from the point measurement of a radiosonde or an "area-average" measurement of a microwave sounder.A single GPS radio occultation measurement represents a weighted average over a "pencil-like" volume (the Fresnel volume).The effective horizontal scale of the measurement along the ray is approximately several hundred km, while the cross-ray scale is quite small (~1 km).Another important issue is related to the fact that over regions with significant vertical refractivity gradients, more than one ray may arrive at a given receiver at the same time, a problem known as multipath propagation (Gorbunov and Gurvich 1998).Under such conditions, the retrieved refractivity profiles may contain significant errors.To assimilate GPS radio occultation data effectively into a weather prediction model, one needs to correctly process them and to properly account for the measurement characteristics and measurement errors.
Over the past five years, a considerable amount of effort has been devoted to developing an effective strategy for the assimilation of GPS radio occultation data into weather prediction models.Eyre (1994) concluded that the optimal strategy is to directly assimilate the bending angle data into weather prediction models.Zou et al. (1995) tested the impact of idealized refractivity data using a dry version of the MM5 four-dimensional variational (4DVAR) system.Kuo et al. (1997) performed a set of observing system simulation experiments (OSSEs) to assess the impact of GPS refractivity data on the prediction of an extratropical cyclone.These papers concluded that the GPS refractivity data have the potential to significantly improve the quality of the model initial conditions and the subsequent forecasts.Zou et al. (1999) incorporated forward ray-tracing code into a three-dimensional variational (3DVAR) data assimilation system and assessed the advantages and disadvantages of assimilating refractivities versus bending angles.Matsumura et al. (1999) and Zou et al. (2000) assimilated bending angles from 30 GPS/MET occultation soundings using the NCEP 3DVAR system.Although the amount of the real GPS/MET data was limited, Zou et al. (2000) demonstrated a small positive impact from the assimilation of bending angles.
In the next section, we will discuss the measurement characteristics of GPS radio occultation data and the strategies for the assimilation of these data.Section 3 discusses the error characteristics of the GPS data.Key results from the recent data assimilation studies are reviewed in Section 4. Section 5 discusses remaining scientific issues on the assimilation of GPS radio occultation data, and provides suggestions for future research.

GPS/MET radio occultation measurements
Before discussing various strategies for the assimilation of GPS radio occultation data, it is important to review the basic measurements and the various steps of data processing.A detailed discussion of the GPS radio occultation observations and data processing may be found in a number of papers, see, for example, Hardy et al. (1992), Gorbunov and Sokolovskiy (1993), Gorbunov et al. (1996a), Gorbunov et al. (1996b), Hocke (1997), Kursinski et al (1997), Steiner et al. (1999), andFeng andHerman (1999).Therefore, we only provide a brief summary here.Each of the 24 GPS satellites continuously transmits signals at two L-band frequencies, L1 at 1.57542 GHz (~19 cm) and L2 at 1.2276 GHz (~24.4 cm).In a GPS occultation, the GPS receiver on a low Earth orbit (LEO) satellite tracks the occulted signals from a GPS satellite whenever a GPS satellite rises or sets behind the Earth (Fig. 3).Atmospheric refractivity affects the phase and amplitude of the radio waves, and those phases and amplitudes along with the positions and velocities of satellites are the fundamental observables.If we consider wave propagation in terms of rays, then the effect of the atmosphere can be characterized by the total bending angle α, as a function of the impact parameter, a, (Fig. 3) which are derived from measurements of the phase and amplitude of the received signal.The bending angles for the L1 and L2 signals are calculated from the slope of phase front under the assumptions of single ray propagation and spherical symmetry of refractivity.The first assumption allows for relation of the observed Doppler frequency shift to satellite velocities and directions of a ray at satellite positions.The second assumption puts an additional constraint on the directions of a ray.This permits a solution for both directions and thus the calculation of the bending angle and impact parameter of the ray (Vorob'ev and Krasil'nikova 1994).The center of sphericity is put into the center of local curvature of the reference ellipsoid under estimated tangent points of rays (Syndegaard 1998).In the lower troposphere, especially in the tropics, strong gradients of refractivity induced by gradients in water vapor may cause multipath propagation.At present, two radioholographic techniques are used to account for multipath propagation (Karayel and Hinson 1997;Gorbunov and Gurvich 1998;Mortensen and Hoeg 1998;Gorbunov et al. 1999;Hocke et al. 1999;Sokolovskiy 2000).These techniques use both phase and amplitude of the received signal.One technique is based on back propagation of the received complex electromagnetic field onto an auxiliary trajectory closer to the limb where multipath effects are expected to be smaller.Bending angles are then reconstructed through the slope of the phase front by assuming single ray propagation.The main problem with this technique is positioning of the auxiliary plane in a single-ray domain.The second technique uses sliding-aperture spectral analysis of the received complex electromagnetic signal, and it associates local spectral maxima in each aperture with multiple rays arriving at that aperture.Frequencies of those spectral maxima, along with position of the aperture, define directions of propagation and thus bending angles and impact parameters of the corresponding rays.The main problem with this technique is detection of the local maxima in the spectrum of a radio occultation signal, which may be rather complicated.Both techniques have been tested using signals simulated with the use of global atmospheric models and were found to satisfactorily solve for multipath propagation (Gorbunov et al. 2000).However, their performance has yet to be tested under more realistic conditions that include very complicated small-scale structures of water vapor in the tropical lower troposphere.
We note that to calculate bending angles from phase, precise positions of the GPS and LEO satellites are necessary.Once ) (a α is calculated, it does not depend on the satellite positions and may be specified with the geographical coordinates of tangent point and azimuthal direction of ray for each given a .This simplifies the following processing (inversion), and, especially, the use of ) (a α for variational assimilation.
The effect of the ionosphere is removed by utilizing the 1 st order dependency of the ionospheric refractive index on the frequency of a signal 2 ) 1 ( , and the approximate account for L1 and L2 ray separation (Vorob'ev and Krasil'nikova 1994) Equation (1) provides sufficient accuracy in the troposphere and lower stratosphere.However, the residual errors after (1) introduce the main error source in the upper stratosphere, even under moderate ionospheric conditions.This is due to ionospheric irregularities and the impossibility of accurate account for the ray separation.Syndegaard (2000) suggests a more complicated technique to better solve for the ray separation.

Initialization (optimization) of the bending angles
At high altitudes the observational noise in ) (a α , mostly the residual ionospheric noise, overshadows the signal in ) (a α produced by density variations in the neutral atmosphere.To minimize the effect of that noise on the reconstructed refractivity at lower altitudes after integral Abel inversion (see below), the derived observational ) (a α are routinely replaced by climatological values of ) (a α above some altitude.Another way is to solve for an optimal vector ) (a opt α based on the observational vector and climate vector, each accompanied by its error characterization (Sokolovskiy and Hunt, 1996).This technique results in a smooth transition from observations at low altitudes to climatology at high altitudes and it minimizes the error propagation downward in a statistical sense.The optimization is not necessary if ) (a α is used for variational assimilation.

Calculation of refractivity (Abel inversion)
In a spherically symmetric atmosphere, the refractive index ) (r n is related to ) (a α through the integral transform pair (Fjeldbo et al., 1971); the second one is called Abel inversion, where x = rn(r ) is the refractional radius.After n(x) is reconstructed from (2b), n is calculated as a function of altitude over reference ellipsoid, z = x / n(x) -r e , where r e is the local curvature radius of reference ellipsoid.

Reconstruction of pressure, temperature and water vapor
With the neglect of small terms related to liquid or frozen water, the refractivity N = (n-1) × 10 6 in the neutral atmosphere is related to atmospheric pressure, P, temperature, T, water vapor partial pressure, P w through (3) Equation ( 3) and the hydrostatic equation are sufficient to derive P and T from N in a dry atmosphere.In a moist atmosphere additional information is necessary.Because the refractivity is more sensitive to w P than to P or T , it is more feasible to reconstruct w P by using independent estimates of P or T from models or other sources, rather than vice versa (Ware et al., 1996).A better approach is not to use exact P or T from a model, but to solve for the optimal one-dimensional atmospheric state , each accompanied by its own error characterization (1DVAR; Kuo et al. 1998).
In summary, temperature and/or water vapor estimates may be obtained from GPS radio occultation observations according to the following steps: a. Measurements of phase and amplitude of the two GPS radio signals (L1 and L2).b.Calculation of bending angles for L1 and L2, respectively, based on the precise position and velocities of the GPS and LEO satellites.c.Removal of ionospheric effects and reconstruction of the neutral atmospheric bending angles as a function of impact parameter.d.Optimization of the observational bending angles with the use of climatology at high altitudes.e. Reconstruction of the vertical profile of refractivity using the Abel inversion technique.f.Inference of water vapor (or temperature) profiles from the refractivity profiles with ancillary data (e.g. from independent observations or a global analysis/prediction).

Assimilation strategy
Before attempting to assimilate GPS radio occultation data, we must determine what variable to assimilate.In principle data at the end of each of the above-mentioned steps can be used in the assimilation.With different assumptions used in the various stages of data reduction, and with different amounts of computation time required to assimilate the GPS radio occultation data in different forms, this is not a trivial question.To a large degree, which variable to assimilate also depends on the "capability" of an assimilation system.
The GPS radio occultation soundings are time-continuous measurements.In order to properly incorporate the GPS data into a numerical weather prediction model, it is desirable to perform fourdimensional data assimilation (FDDA).Currently, there are two approaches for time-continuous assimilation: one is Newtonian relaxation (sometimes called nudging), and the other is four-dimensional variational data assimilation.In the relaxation method (Anthes et al. 1974;Kuo and Guo 1989;Stauffer and Seaman 1990), an extra forcing term is added to the model predictive equation to "nudge" or "relax" the model solution toward the observations.For a given variable, this can be written as the following: where F represents the usual forcing terms in the model equations, G is a positive relaxation coefficient, and W a positive weighting factor that expresses the confidence level with respect to the accuracy of the observation, α o is the observations.Since the data can only be assimilated through the predictive equations, the relaxation method can only assimilate "direct" observations -measurements in the form of the predictive variables (i.e., wind, temperature, water vapor, pressure).Therefore, to assimilate the GPS radio occultation data using the relaxation method, it is necessary to reduce the measurements to water vapor and temperature.
As discussed by Eyre (1994), the variational data assimilation (Lorenc 1986; Le Dimet and Talagrand 1986) is a more promising technique for the assimilation of GPS radio occultation data.In variational assimilation the analysis of the atmospheric state x has to be found by minimizing a cost function of the form: where b x is the background estimate of x (e.g. from a short-term forecast), B is the covariance of expected error in b x , y is a vector of observations, } {x H is the observation operator, which calculates the observation that would be made given the state x , O is the covariance of expected error in the observations, F is the covariance of expected error in the observation operator, and ) (x J c represents an additional "penalty function", through which other dynamical or physical constraints can be imposed.For a weather prediction model, x is a vector of predictive variables at all the model grid points, y contains all the observations made within the time window of the assimilation.In order to obtain the optimal initial condition ( 5), the gradient of J with respect to the initial condition ( x ): x ) that provides the "optimal" fit between observations, the background, and the penalty terms.
An important advantage of the variational data assimilation is that the comparison of measurement with analysis (model) takes place in the space of the observed variables, not in that of the analysis (model) variables.Through the use of the observation operator, one can assimilate any observations that can be expressed as a function of the analysis (model) variables.For example, we can assimilate refractivities or bending angles instead of water vapor or temperature.As indicated in (5), to properly assimilate an observation we must have estimates of the expected errors in the measurements and the observation operator.To answer the question of "in what form should the GPS radio occultation observations be assimilated", we consider several factors: (a) complexity of the observation operator, (b) assumptions used in the data reduction, (c) additional data required for assimilation, and (d) estimates of the measurement errors.In principle, the GPS radio occultation data can be assimilated in a form ranging from the raw frequency and amplitude measurements to the derived temperature and water vapor profiles.The advantages and disadvantages of assimilation GPS radio occultation data in various forms are summarized in Table 1.

Raw phase and amplitude data
The advantage of assimilating raw L1 and L2 phases and amplitudes is that there is no pre-processing, and thus no assumptions are used.However, to model these data in the presence of multipath propagation, the use of diffractional observation operator and its tangent linear and adjoint models in a data assimilation system is necessary.This is a very complicated task that would require tremendous amounts of computer time and it is not realistic.It is more feasible to first account for multipath propagation by using radioholographic techniques, and to obtain either phase on the auxiliary trajectory (back propagation), or bending angles (spectral technique) and then to assimilate them.If multipath propagation is not significant (which is the case normally above the troposphere), then the raw phase may be assimilated (amplitude does not contribute significant complementary information).

Raw phase data
Here the assumption of single path propagation is used, and thus the diffractional observation operator is not necessary, instead ray tracing may be applied.However, in this case the observation operator is still too complicated.This is because phase (and excess phase as well) depends on the positions of the satellites.This requires the assimilation model to keep track of the information on the satellite positions, which complicates the processing.Thus, the observation operator must find a ray connecting two given points, which requires the use of an iterative (shooting) technique due to the nonlinearity of the ray equations.Shooting methods are computationally expensive.Moreover, when using L1 or L2 phases separately, an ionospheric model is necessary to account for the ionospheric effect on phase.None of the existing ionospheric models can predict the ionospheric effect as well as it can be removed by dual frequency correction.This correction may be performed with L1 and L2 excess phases similarly to bending angles (1).But then before finding the ray, the positions of satellites must be corrected for the ionospheric bending, which again complicates the problem.

L1 & L2 bending angles
An advantage of assimilating the bending angles is that they may be specified by the impact parameter, geographic coordinates of tangent point and azimuth of a ray, instead of satellite positions.This is due to the fact that bending is accumulated only on the section of a ray inside a refractive medium (atmosphere and ionosphere), while phase is accumulated along the entire ray path, including space.Thus, instead of finding the ray which connects two given points we have to find the ray with given initial conditions at the tangent point, and therefore shooting methods are not necessary.Ray tracing may be started at tangent point in both azimuth directions, then continued up to a certain altitude and stopped.Although bending angles are calculated from the observational data with the assumption of spherical symmetry, this assumption does not introduce problems because the same assumption is made in the observation operator.Instead of using true directions at both ends of a ray, bending angle has to be calculated from one true direction, and another one reconstructed with the use of the Snell's law in the assumption of spherical symmetry.Therefore, the "comparison" between modeled and observed bending angles (through the calculation of the cost function) should not be affected by this assumption.

Ionospheric free bending angle
The assimilation of the neutral atmospheric bending angle is superior to the assimilation of the bending angles for each frequency.It has all the advantages of bending angle assimilation (with less complicated observation operator), and none of the disadvantages associated with the modeling of the ionospheric effects in the assimilation model.By calculating the linear combination of L1 and L2 bending angles, one can remove the ionospheric effects better than using the data from a global ionospheric model.Ray tracing in the observation operator is necessary up to much lower altitudes (up to the top of the neutral atmosphere only) than when assimilating L1 or L2 bending angles separately.

Refractivity
The calculation of refractivity requires the use of the Abel inversion.As an approximation, the refractivity that is derived from the bending angles can be regarded as the "averaged" refractivity along the ray path; this results in heavier weight given to the region near the ray perigee point (because of density variation with height).It does not represent a "point value" of refractivity at the ray perigee (except in the special case of spherical symmetry).For the assimilation of refractivities, one has two choices.The first is to develop an observation operator that follows exactly the same steps of how the observed refractivity is derived.This means that the observation operator would include the bending angle calculation and the Abel inversion.For the sake of convenience, we will call this observation operator N GPS H .The second approach is to assume that the observed refractivity represents an acceptable "approximation" to the local refractivity at the ray perigee point, can potentially lead to significant errors over regions where the horizontal gradient of refractivity is large (and therefore the local spherical symmetry assumption is not valid).This point will be discussed later.

Water vapor or temperature
The use of water vapor or temperatures derived from GPS radio occultation observations is advantageous in that the data are in a form that is the same as the model variable, and can be readily used by a wide variety of objective analysis and data assimilation systems.However, the disadvantages associated with the assimilation of water vapor or temperature data are numerous.First, we have to assume GPS N is the same as Local N . We then have to use ancillary data (temperature or water vapor) from an independent source to derive the water vapor or temperature profile.The accuracy of the derived variable is strongly affected by the accuracy of the ancillary data used.

Discussion
Based on the above analysis, only two approaches appear to be feasible for the assimilation of GPS neutral atmospheric observations: (1) the ionospheric free atmospheric bending angles or (2) the refractivity.Indeed, the assimilation research over the past five years has concentrated on the evaluation of these two approaches.These studies will be reviewed in the following section.Before that several notes are necessary.
It follows from (2a) that the bending angle depends on the refractivity above the ray tangent point.Thus bending angles calculated at the top of an atmospheric model do not contain significant information about the atmospheric state reproduced by that model (except that introduced by extrapolation of the model fields upward).On the other hand, refractivity retrieved from the observational data does contain information about true atmospheric state (unless observational bending angles are overshadowed by noise right above the top model altitude and it results in heavy weight of climatology after their initialization prior to Abel inversion).Thus it appears feasible (and actually desirable) to assimilate Abel-retrieved refractivity (treating it as local) at high altitudes regardless of what variable was assimilated at lower altitudes.
Atmosphere-induced variations of ) (a α or ) (r N as reconstructed from radio occultation data may be characterized by a spectrum.When sampling these observations, the sampling rate must be not less than the spectral bandwidth to avoid aliasing (internal radio occultation sampling rate is normally large enough).In the troposphere, where radioholographic techniques with increased vertical resolution are applied, the required sampling rate (discretization) may be rather large, much larger than the spectral bandwidth of variations of ) (a α outside the model reproducible bandwidth is not assimilated, and the use of high-resolution observations is not warranted, especially in the case of assimilation of ) (a α where calculation of a large number of rays is very computationally expensive.In this case, sampling of the observation data and discretization in the observation operator (i.e., the number of rays) should be reduced to be consistent with the resolution of the model, and anti-aliasing filtering has to be applied prior to sampling of the observational data.

ERROR CHARACTERISTICS OF RADIO OCCULTATION OBSERVATIONS
Error characterization is necessary to assign weights for observational data (introduced by observational error covariance matrix O in ( 5)).Without suitable weights one cannot obtain a final analysis that provides an "optimal" fit to the background and the observations.Unfortunately, all GPS radio occultation data assimilation studies performed so far use a highly simplified form of observational error covariance matrix (or weighting), due to lack of information on the error characterization of actual radio occultation, especially, in the lower troposphere.There are different error sources of radio occultation data that are dominant at different altitudes.To make optimal use of GPS radio occultation data for numerical weather prediction, we need to obtain realistic error statistics and incorporate this information into a data assimilation system.
One group of error sources includes thermal noise (mostly from receiving antenna), orbit and clock errors and local multipath.The magnitudes and the spectral content of these errors are well known (Hardy et al. 1994;Kursinski et al. 1997).These errors are not very important for weather forecast models since at the altitudes covered by the models (normally below 35km), the corresponding errors of bending angles are significantly smaller than their variations due to the model atmosphere (the signal).
The most important error source in the stratosphere appears to be the residual noise after the ionsopheric correction.Figure 4 shows the deviations of the observed bending angles as a function of ray tangent point from those calculated with the use of climatology (CIRA -Committee on Space Research International Reference Atmosphere Model) for 3 GPS/MET occultations.Dashed lines correspond to 20% of the climatology and are estimates of the span of the expected neutral atmosphere-induced variations of bending angles.The magnitude of these variations should exponentially decrease upward.However, as seen in Fig. 4 (especially 4c), above some altitude the magnitude of bending angle variations remains almost constant and it is significantly larger than expected.This is mostly due to the ionospheric noise which overshadows the neutral atmosphere-induced variations of bending angles.GPS/MET data processing includes low-pass filtering which limits the spectral domain of the signal.Thus, the residual noise which is visible in Fig. 4 occupies the same spectral domain as the signal.As seen in Fig. 4, its magnitude may be very different for different occultations, depending on the ionospheric conditions.A good thing about this noise is that its magnitude is expected to be almost the same at all altitudes for a given occultation.This is because normally the ionospheric conditions are changing at much larger scales than in the neutral atmosphere.This allows evaluation of the magnitude of this noise for each occultation individually at those altitudes where it certainly overshadows the neutral atmospheric variations of bending angles, normally above 60 km.Then this magnitude may be used as a parameter for the observation error covariance matrix for the given occultation.For GPS/MET, the magnitude of bending angle noise (rms deviation of bending angles) was estimated for each occultation between 60 and 80 km (because above 80-90 km the sporadic ionospheric layers may significantly contribute to the noise).Figure 5 shows the distribution of the magnitude of noise for one day of observations, October 21, 1995.It is seen that the mean magnitude of the bending angle noise is about (1-2) 10 -6 rad, however, for certain occultations it is almost 10 times larger.
Horizontal inhomogeneity of refractivity does not introduce an error when assimilating bending angles, however, it does when assimilating Abel-retrieved refractivity as a local variable.A fractional RMS magnitude of this error, as estimated from numerical simulations with global atmospheric models (Gorbunov et al., 1996, Kursinski et al., 1997) is about ~0.2% at altitudes 10-30km.This indicates that direct assimilation of the Abel-retrieved refractivity as a local may be feasible at those altitudes.
The lower troposphere, especially in the tropics, introduces the main problem for error characterization.Kursinski et al., (1997) estimate fractional RMS refractivity errors in the lower troposphere as ~1%.At present it may be considered as a coarse estimate of the lower tropospheric errors.However, this estimate is based on numerical simulations with ray tracing, which strictly is not applicable for the complicated small-scale refractivity structures in the tropical lower troposphere.More realistic simulations based on modeling of wave propagation through high resolution 2D (3D) tropospheric refractivity structures have yet to be done.Coarse error characterization based on comparisons of the retrieved refractivities from GPS/MET to radiosondes would result in significant error overestimation (see discussion of the problem of the GPS/MET negative N-bias by Rocken et al., 1997).The large errors in the lower troposphere may be mostly attributed to the signal tracking technique applied in the GPS/MET receiver.For the COSMIC mission, an alternative tracking technique which is more stable to complicated signal dynamics has to be applied in the lower troposphere (Sokolovskiy, 1999).Testing of the COSMIC receiver by simulated signals prior to launch is important for minimization and better evaluation of the lower tropospheric errors.

Refractivity assimilation
As discussed earlier, the assimilation of refractivity represents a good compromise between the assimilation of bending angles and the assimilation of the derived temperature or water vapor profiles.Zou et al. (1995) were the first to assimilate refractivities in a set of OSSEs using an adiabatic version of the MM5 model (Dudhia 1993;Grell et al. 1994).They first conducted a 20-km simulation of a winter cyclone in March 1992 over the continental United States and used that as the control experiment (the so-called "natural run").Vertical profiles of atmospheric refractivity were extracted from the control experiment at selected temporal and spatial resolutions to simulate atmospheric refractivity observations in an idealized setting.The simulated refractivities were then assimilated into a 60-km MM5 model using the 4DVAR approach with a two-hour assimilation window.The temporal resolution of the simulated refractivity data ranged from 5 min to 30 min, and the spatial resolution varied from 60 km to 180 km.The assimilation was performed with and without the use of auxiliary temperature and wind fields.Zou et al (1995) concluded that the assimilation of atmospheric refractivity is effective in recovering the vertical profiles of water vapor.To a lesser degree, the assimilation of refractivity also improved the temperature field.In order to provide a more realistic simulation of the refractivity observation, in one experiment they averaged the simulated refractivity data along a distance of about 240 km.They showed that the assimilation of "averaged" refractivity still had a positive impact.
Although the results of Zou et al (1995) are interesting, there are several limitations associated with this study: (1) a model without moist physics is used, (2) the simulated refractivity observations are highly idealized both in terms of spatial and temporal resolution and in terms of measurement characteristics (no ray-tracing was used to simulate the bending or the use of Abel transform to derive refractivity), and (3) the relatively short assimilation window of 2 hours.Nevertheless, this paper was the first to show that the assimilation of refractivity can lead to significant improvement in moisture analysis.Kuo et al. (1997) continued this line of research by conducting OSSEs using a moist version of MM5 and its 4DVAR system.They first conducted a 90-km simulation of an extratropical cyclone, known as the ERICA (Experiment on Rapidly Intensifying Cyclones over the Atlantic) IOP4 storm with a domain that covered most of the northern hemisphere.The results from the control simulation were then used to simulate GPS/MET refractivity observations with varying spatial resolution and distribution.The simulated observations were assimilated into a 180-km MM5 model using the MM5 4DVAR system that includes the adjoint of moist physics.The assimilation was performed during a 6-h window, from 1200UTC to 1800 UTC 3 January 1989, which was 12 to 6 hours prior to the start of rapid cyclogenesis.The simulated GPS/MET refractivity soundings had a spatial resolution of either 180 km or 360 km, and they were either located on regular grids or randomly distributed with varying orientation.The refractivity data were assumed to be available only at the beginning and ending times of the assimilation window.In an experiment (E1) in which refractivity soundings were assumed to be available on a regular 180km grid, Kuo et al (1997) showed that the assimilation of refractivity produced significant improvement, not only in the temperature and moisture fields, but also in the winds.In particular, the assimilation of refractivity induced sizeable analysis increments along the upper-level jet streams (Fig. 6).With improved definition of upper-tropospheric temperature and wind fields, the assimilation of refractivity produced an improved description of the potential vorticity (PV) field.Although the refractivity is not directly related to the winds, through internal model dynamics, improvements in the model temperature fields due to refractivity assimilation produced improvements in the winds.Positive PV anomalies in the upper troposphere and thermal anomalies in the lower troposphere are important ingredients for surface cyclogenesis (Davis and Emanuel 1991).Although the impacts on upper-level PV and low-level thermal anomalies are small and subtle initially, they grow rapidly as the baroclinic cyclone develops.By the end of the 48-h forecast, the experiment with refractivity assimilation produced a surface cyclone that was 13 mb deeper than that of the no-assimilation experiment (Fig. 7).
In the GPS/MET experiment, only a small fraction of the radio occultations penetrated to the lowest 1 km of the atmosphere (Rocken et al. 1997).Moreover, the receiver tended to lose track in the lower part of the troposphere, giving corrupted data.In order to assess the relative importance of low-troposphere GPS/MET radio occultation data, Kuo et al. (1997) performed an experiment (named E2.upper) in which the simulated refractivity data were assumed to be unavailable below 3 km.They found that restricting the refractivity data to altitude 3 km and above considerably reduced the impact of refractivity assimilation on cyclone prediction.The effect on cyclone deepening was reduced from 13 mb in 36 h to 6 mb in 36 h.This again emphasizes the importance of obtaining high quality radio occultation data in the lower troposphere.
Another interesting result from the Kuo et al (1997) study is the different behavior of the improvements in temperature and moisture fields as a result of refractivity assimilation.The assimilation of refractivity data reduced errors in the initial moisture analysis by 50% during the assimilation period (panels (j) and (k) of Fig. 8).However, six hours into the forecast (following the assimilation period), this impact was reduced to about 15%, and it was further reduced to 10% by 12 hours.By 24 hours, the improvement in the moisture field is largely lost.In contrast, the impact on the temperature field, although small initially, continued to increase with time.This suggests that the maximum value of GPS/MET assimilation in the moisture analysis and precipitation forecast is likely to be confined to relatively short periods of time (i.e. 6 to 12 hours).In order to sustain the impact, we need to continuously assimilate GPS radio occultation data as they become available.The traditional, 12-h update cycle (in which analysis is updated once every 12-h) will not be able to derive the maximum benefits from the GPS/MET radio occultation data.Although Kuo's et al (1997) study provides useful information on the potential impacts of GPS data assimilation, this study also suffers several shortcomings.Aside from the usual "identical-twin" problem associated with OSSEs (Arnold and Dey 1986), the major shortcoming is the unrealistic simulation of the GPS radio occultation soundings.Kuo et al (1997) did not perform ray tracing to simulate bending angles, nor did they invoke Abel transform in the calculation of refractivity.The simulated GPS refractivity was simply the "local" refractivity calculated from the model grid point data.
Moreover, it is not realistic to expect all the data are available at the beginning and ending time of the assimilation window and nothing in between.In reality, data are continuously received, and should be continuously assimilated.In their study, the density of "simulated GPS refractivity sounding" would have required a 100-satellite constellation for GPS soundings at 180-km resolution, and a 24-satellite constellation for 360-km resolution.Nevertheless, they showed that refractivity data could be used to improve global and regional analysis and prediction, and provided some physical explanation on how this was accomplished.Zou et al. (1999) developed a methodology for assimilating atmospheric bending angles.They first developed a ray-tracing observation operator that links the atmospheric state to the GPS bending angle measurements.They then developed the tangent linear and adjoint models of the ray-tracing observation operator.They showed that the assimilation of GPS bending angles in a variational analysis induced changes in the temperature and moisture fields that were not limited to the ray perigee point, but occurred in an elongated band of ± 300 km in its occultation plane.On the model σ -levels, the changes in the analysis as a result of the assimilation of one occultation sounding occurred in an area of about 600 km x 600 km, centered around the ray perigee point.Zou et al. (1999) also discussed the advantages and disadvantages of assimilating refractivities versus bending angles.By using a ray tracing observation operator, the occultation process can be better simulated in the model atmosphere, and the GPS bending angle and its derived refractivity can be calculated based on the model data.The GPS refractivity ( GPS N ) derived from the Abel transform can be compared with the local refractivity (

Local N
), which is calculated simply from the model grid point data at the ray perigee point.This comparison offers insights on the potential errors associated with the assimilation of refractivity (treating the GPS derived refractivity as local refractivity at the ray perigee point).Zou et al (1999) compared 62 vertical profiles of GPS N and Local N calculated from the NCEP global analysis at 1200 UTC October 11, 1995 (Fig. 9).They found that, in general, the difference between local and GPS refractivity is not significant above 20 km.The differences become progressively greater in the lower part of the atmosphere.Below 10 km, differences typically range from 3 to 5 refractivity units.In some extreme cases, it can be as large as 30 units (approximately 10% error).Zou et al. (1999) found that for most soundings, the local atmospheric refractivity (

Local N
) is close to the GPS-derived refractivity (

GPS N ).
However, for soundings that are located over regions with significant horizontal refractivity gradients, errors introduced by local spherical symmetry assumption cannot be ignored.In such cases, treating GPS derived refractivity as local refractivity may introduce large errors, and a direct assimilation of bending angles should be used.It is important to recognize that this comparison was made based on NCEP global analysis (T62L28) that has relatively low horizontal and vertical resolution.It is likely that the differences between

GPS N and Local N
at the ray perigee point can be larger in reality (which will have much more small scale variations).It would be desirable to repeat such calculations using the results from a highresolution mesoscale model.
The first successful assimilation of real GPS radio occultation bending angles was carried out by Matsumura et al. (1999) and Zou et al. (2000).In these studies, they incorporated 30 GPS/MET soundings of bending angle into the NCEP 3DVAR data assimilation system.Zou et al. (2000) found that on average, the adjustments in temperature and specific humidity were approximately 0.4 o C and 0.6 g kg -1 , respectively.Individual changes due to bending angle assimilation were as large as 4 o C and 4 g kg -1 (Fig. 10).With only 30 GPS/MET soundings, it was difficult to demonstrate significant impact of GPS bending angle assimilation on the entire global analysis, as the 30 GPS/MET soundings were outnumbered by the conventional upper-air observations (682 temperature soundings, 539 specific humidity soundings, and 1166 wind reports) that were used together with the GPS radio occultation soundings.Nevertheless, Zou et al (2000) found that even over data-rich regions (where there are many radiosonde observations) the assimilation of GPS/MET bending angles provided a slight improvement to the analysis of wind and temperature.In an effort to demonstrate the impact of assimilation of bending angles, Zou et al. (2000) considered three occultation soundings (GPS-1, GPS-2, and GPS-12) that were relatively close to each other (separation between occultation soundings ranges from 334 km to 593 km).The effect of assimilation of bending angles was demonstrated by comparing the analysis resulting from the assimilation of a single sounding with independent occultation observations that were withheld from the analysis.They found that the assimilation of GPS-1 bending angles reduced the errors of the first guess field at the location of GPS-2 and GPS-12 by 8% and 55%, respectively.The improvements due to bending angle assimilation were also reflected in the analyzed temperature and specific humidity fields.

4.3
Computational efficiency Zou et al. (2000) showed that the assimilation of 30 GPS/MET soundings (each with 76 rays) together with the use of the radiosonde data cost about 7990 CPU seconds on a CRAY3 at NCEP, with 7220 CPU seconds spent on the assimilation of bending angles alone.This was more than 10 times the cost for the assimilation of all the other available upper-air sounding data.In contrast, the assimilation of refractivity from the same 30 soundings only cost 657 CPU seconds.Clearly, direct assimilation of bending angles is expensive.It is highly desirable to find ways to reduce the cost.In particular, it would be desirable to find ways to optimally "mix" the assimilation of bending angles with refractivity.One strategy would be to assimilate refractivity in the upper part of the atmosphere, and to assimilate bending angles in the lower part of the atmosphere.As discussed in Section 2, one needs to consider refractivity assimilation near the top of the model simply because of the fact that the bending angles cannot provide useful information about the model.As discussed earlier, the assimilation of refractivity can lead to significant errors over regions with significant horizontal refractivity gradients.This is particularly true in the lower troposphere.In view of the fact that the differences between

GPS N and Local N
decreases with height (see Fig. 9), it would be desirable to use refractivity assimilation above a level where local spherical symmetry assumption does not produce significant errors.Zou et al. (2000) experimented with a several different strategies to use a mixture of refractivities and bending angles.These experiments result in computational savings ranging from 31% to 46%, while producing no visible degradation to the results.
A typical radio occultation takes about one minute.With a receiver taking 50 Hz data, there are a total of ~3,000 samples per sounding.Because of storage space limitations and due to the fact that the "effective vertical resolution," largely determined by the vertical size of the Fresnel zone, is far less than that defined by the 3,000 samples, the number of rays was reduced to ~300 rays through smoothing and interpolation in the GPS/MET data base.
An important question has to be asked for data assimilation: "For a global model with about 30 vertical levels, how many rays should be used?"This is not a trivial question.Two factors need to be considered.First, we want to extract the maximum useful information from the data that can be supported by the model's horizontal and vertical resolution.The other competing factor is computing cost.As discussed earlier in Section 2, an appropriate approach would be to filter the observational data to make information content consistent with that can be reproduced by the model in the observational space.This would reduce the computational cost and avoid aliasing.
Aside from the fact that one needs to find more efficient ways to calculate ray tracing observation operator and its adjoint, a possible way to reduce the computational cost is through "preconditioning."For example, one can attempt to assimilate refractivity data at all altitudes first as a mean of "preconditioning."By assimilating the refractivity data from all the radio occultation soundings, this would accomplish basically 90% of the work at about 10% of the cost.We can then assimilate the bending angles below 10 km, using the results of refractivity assimilation at those altitudes as the "guess" field.Since the "guess" field already contains much of the information from the radio occultation soundings, the number of iterations for bending angle assimilation can be significantly reduced.

SUMMARY AND DISCUSSION
Significant progress has been made over the past five years in developing an effective strategy for the assimilation of GPS radio occultation data into weather prediction models.The basic conclusions from the studies conducted so far can be summarized as the following: • The use of GPS radio occultation data has the potential to significantly improve the accuracy of global and regional analysis and weather prediction.• Neutral atmospheric bending angles and refractivity are the two forms of the GPS radio occultation data that appear most suitable for data assimilation, and the variational data assimilation is likely the best way to assimilate GPS radio occultation data.• A reasonable strategy is to assimilate bending angles in the lower atmosphere (below 10 km) and to assimilate refractivities in the upper atmosphere (above 10 km), due to the following considerations: o Over regions with significant horizontal refractivity gradients (normally, in the lower troposphere), the assumption of local spherical symmetry is not valid.Treating GPS refractivity as local refractivity can induce significant errors.o The assimilation of bending angles is an order of magnitude more computationally expensive than the assimilation of refractivities.o Near the top of a model, bending angle data do not provide useful information about atmospheric refractivity for the model.o Between 10 km to 30 km, treating GPS refractivity as local refractivity does not introduce significant errors (statistically on the order of rms 0.2%).
Although significant progress has been made in developing an effective strategy for the assimilation of GPS radio occultation soundings, many issues remain to be addressed.In particular, we have yet to convincingly demonstrate the "value-added" of GPS radio occultation data to the existing global observing system.The main problem is, with a single receiver and with sub-optimal configuration typically associated with a "piggy-back" experiment, the GPS/MET experiment does not provide sufficient number of soundings for a given day.On a "good" day, GPS/MET provided approximately 125 soundings.This is less than 10% of the number of radiosonde observations.It is exceedingly difficult to demonstrate the impact of GPS radio occultation sounding through real-data study until a constellation mission such as COSMIC is launched.Until then, OSSEs are powerful tools that allow us to assess the impact of radio occultation observations.With the ability to simulate potential observations from a hypothetical constellation with any number of satellites, OSSE can provide very useful insights on the design of an "optimal" GPS/MET constellation.For example, OSSEs can be used to address the question of "what is the optimal number of LEOs needed to improve global weather prediction?"One can also assess the relative importance of different observing systems (i.e., GPS radio occultation versus lidar wind satellite).Unfortunately, the OSSEs conducted so far are limited in scope and are subject to several shortcomings.To perform a more realistic assessment of the potential impact of GPS radio occultation soundings through OSSEs, one needs to consider the following factors: (i) Realistic simulation of radio occultation.GPS radio occultation processes should be simulated through ray tracing, and the GPS refractivity should be derived through Abel transform.It is not sufficient to simulate the GPS refractivity simply with "local" refractivity calculated from the model grid point data.(ii) Realistic simulation of potential soundings from a hypothetical network.One should realistically simulate the distribution of radio occultation soundings in time and space for a given design of a constellation.Distributing the data randomly or on regular grids is not appropriate.This should be easily accomplished through the use of orbit calculation software.(iii) Adding realistic measurement errors to the simulated observations.For OSSEs to be useful, the simulated observations need to contain realistic measurement and data processing errors.(iv) Avoid the identical twin problem.It is well known that when the same model is used for both creating and assimilating simulated observations, the results tend to be overly optimistic.It would be desirable to use one model to generate the simulated observations, and use another independent model to assimilate the observations.(v) Use of high-resolution models.The simulation of GPS observations is limited by the horizontal and vertical resolution of the "natural run."For example, radio occultation signals and inverted data simulated with the use of global atmospheric models have less complicated structure than real observational data because the small-scale structures which are present in the real atmosphere are not reproduced by the global models (Matsumura et al. 1999).Therefore, it is highly desirable to simulate GPS radio occultation soundings using results from a high-resolution mesoscale model.This will allow a more realistic simulation of the GPS bending angle and refractivity profiles.(vi) The value added by GPS radio occultation data.To demonstrate the value added by GPS radio occultation data, the assimilation of GPS radio occultation data needs to be conducted with and without data from other current and future observing systems (i.e., lidar wind satellite, microwave sounders, etc).This would allow a fair assessment of the value of a GPS radio occultation constellation mission, such as COSMIC, to global and regional weather prediction in the presence of the current and future global observing system.
In order to derive the maximum benefits of COSMIC in global and regional scale weather prediction, the research tasks suggested above should be conducted before the launch of COSMIC.

Fig. 2 .
Fig. 2. (a) Projected COSMIC radio occultation soundings over the northern hemisphere every 24 h and (b) the present radiosonde sites.

Fig. 3 .
Fig. 3. Instantaneous occultation geometry for the Global Positioning System (GPS) and a low Earth orbiter(LEO) satellite.The bending angle α and impact parameter a are determined from the Doppler frequency shift, spacecraft position, and spacecraft velocity measurements.2.1.1Reconstruction of bending angles for both L1 and L2 signals.
6) must be calculated.The minimization process starts with the use of forward observation operator ( } {x H ) that transforms the analysis (model) variables ( x ) to the observation space and simulates the observed quantity in the same form, } {x H . Observations (y) are then compared with the simulated quantity ( back to the space of the analysis (model) variables by making use of the adjoint ( projected residual for all observations, together with the background term ) the gradient of the cost function.The error covariance matrices, B , O , and F control weighting for the related terms, respectively.Through the use of nonlinear optimization algorithm, analysis increments are determined that allow the cost function ) (x J to be reduced.This procedure is repeated until a local minimum is obtained.The final analysis is the initial condition ( * o

.
This assumption greatly simplifies the development of the observation operator, as we can easily develop N Local H according to (3).It is obvious that there are no advantages of using N GPS H over bending angle assimilation, as more computation is required, and we do not derive any new information.The use of N Local H greatly reduces the computation of observation operator and its adjoint, and represents considerable computation saving.However, treating the GPS N as Local N

Fig. 4 :
Fig. 4: Deviation of the observed bending angles as function of ray tangent point from those calculated with the use of climatology (CIRA) for 3 GPS/MET occultations.

Fig. 5 .
Fig. 5. Magnitude of residual ionospheric bending angle noise for all occultation soundings taken on 21 October 1995 during GPS/MET.The three occultations used in Fig. 4 are marked with "+".

Fig. 6 .
Fig. 6.(a) 300 hPa wind speed (contour interval = 5 m s -1 ) at the end of the assimilation cycle (1800 UTC 3 January 1989) in an experiment (E1) in which the simulated GPS refractivity soundings are located at 180-km regular intervals, (b) the differences in 300 hPa wind fields between E1 and the experiment without data assimilation valid at the same time (1800 UTC 3 January 1989) (from Kuo et al 1997).

Fig. 8 .
Fig. 8.The rms errors of temperature ((a) -(i), in o C) and specific humidity ((j) -(r ), in g kg -1 ) at 6 hours intervals from 1200 UTC 3 January to 1200 UTC 5 January 1989 for experiments with (solid line) and without data assimilation.The data assimilation was performed during the period of 1200 UTC to 1800 UTC 3 January (from Kuo et al 1997).

Fig. 9 .
Fig. 9.A "spaghetti" map showing the 62 vertical profiles of the differences of

Fig. 10 .
Fig. 10."Spaghetti" maps showing the 30 vertical profiles of analysis increments of (a) refractivity, (b) temperature, and (c) specific humidity as a result of bending angle assimilation.(From Zou et al. 2000).

Table 1 .
Comparison of the advantages and disadvantages of assimilation of GPS data in different forms.