1. Introduction

Radio occultation(RO)is a limb-soundingremote-sensing technique used for measuring the physical properties of a planetary atmosphere. When aradio signal travels from a global navigation satellitesystem(GNSS)transmitter to a receiver in the low-earth orbit satellite, it is refracted. The delay or thechange of bending due to the atmospheric refractioncan be precisely measured, and then used to infer theelectron density in the ionosphere and the refractivitycontaining temperature and water vapor informationin the stratosphere and troposphere. Compared withother satellite remote-sensing observations, RO observations have the advantage of high accuracy, high vertical resolution, all-weather global coverage, no systembias, etc. They therefore have great potential value inweather and climate applications(Eyre, 1994; Ware et al., 1996; Kursinski et al., 1997; Anthes et al., 2000;Wickert et al., 2001).

Many advanced numerical weather prediction(NWP)centers, such as the European Centre forMedium-Range Weather Forecasts, have shown thatoccultation data are one of the most important setsof observations in the global observation system, asthey can significantly improve the analysis and forecast quality for both the Northern and Southern Hemispheres, particularly the Southern Hemisphere and the stratosphere where conventional observations arescarce(e. g., Healy et al., 2005; Cucurull et al., 2006;Aparicio and Deblonde, 2008; Buontempo et al., 2008;Cucurull and Derber, 2008; Healy, 2008; Poli et al., 2009; Rennie, 2010). In addition, RO observations donot require bias correction; they can "anchor" the biascorrection applied to satellite radiances and help identify NWP model biases(Collard and Healy, 2003).

Although early occultation experiments(Healy and Thepaut, 2006)proved RO observations' potentialcontribution to NWP, it was not until the launch ofFORMOSAT-3/COSMIC(Formosa Satellite mission3/Constellation Observation System for Meteorology, Ionosphere, and Climate; simply referred to as COS-MIC hereafter)in April 2006 that the operational assimilation of RO data became a reality. The mainreason is that the COSMIC mission consists of six loworbit satellites, which can provide more RO observations(approximately 2000 profiles daily), and the dataare available in near real time, with less than 3 h of latency(Anthes et al., 2008). Besides, two operationalRO missions, i. e., the METOP-A launched in October 2006(Loiselet et al., 2000) and the METOP-Blaunched in September 2012, provide more than 1300profiles daily. Other research RO missions, such asCNOFS and TerraSAR-X(Wickert et al., 2005, 2008), also provide a considerable number of RO observations.

GRAPES(Global/Regional Assimilation and PrEdiction System)is the new generation operationalNWP system of the China Meteorological Administration(CMA). It was developed in response to anincreasingly urgent need for more accurate numerical weather forecasts in recent years(Xue, 2006; Xue and Chen, 2008). The regional and global GRAPESsystems have been run in routine operation or pre-operation in the National Meteorological Center ofCMA. As a new type of satellite observation, the ROdata and their assimilation into the GRAPES are discussed in detail in this paper, including the design ofthe observation operator, the strategy of operationalimplementation, and analysis of the assimilation impact on forecast experiments.

The remainder of this paper is organized as follows. Section 2 provides a simple description of theGRAPES variational assimilation system(GRAPES-Var). Section 3 describes the RO observation preprocessing system. Section 4 elaborates on the RO refractivity observation operator and observation error. Sections 5 and 6 present the results of short- and long-term RO observation assimilation experiments, respectively. A summary of the key findings and further discussion are presented in the final section. 2. GRAPES-Var

The early GRAPES-Var(Zhuang et al., 2005; Xue et al., 2008)did not match the GRAPES forecastmodel in either the type of coordinate used or the analyzed atmospheric state vectors. This caused manyerrors in coordinate and variable transformations. Recently, we have updated the GRAPES-Var from thepressure coordinate to the terrain-following height coordinate(Xue et al., 2012). Namely, the current definition of the coordinate and atmospheric state variables in GRAPES-Var is completely consistent withthe GRAPES forecast model; the initial field is obtained directly after data assimilation and without anyadditional transformation. Therefore, the state vectors are

where u and v are the components of wind, q is specific humidity, and m denotes the mass field, whichis either Exner pressure π or potential temperatureθ. Although GRAPES is a non-hydrostatic model, we do not consider the non-hydrostatic problem ofGRAPES-Var at present, and we choose π as the statevariable and θ as the derivative of π. We also do notanalyze hydrometeors. Figure 1 presents a schematicdiagram of the atmospheric state and analysis variables of new GRAPES-Var in the Charney-Phillipsvertical grid. The expression of the RO observationoperator is required to suit such a definition of thecoordinate and variables. 3. RO data preprocessing

The RO measurements under different data processing stages, from the phase and amplitude, bendingangle, and refractivity, to the retrieval of the temperature and water vapor pressure, can be assimilated. The favorite RO measurements assimilated in an operational NWP system are ionosphere-corrected bending angles and refractivity(Eyre, 1994; Kuo et al., 2000). The calculation of bending angle is affected bythe height of the model top(Aparicio and Deblonde, 2008; Rennie, 2010); therefore, refractivity is the bestchoice for current GRAPES with the relatively lowmodel top of 32. 5 km. 3. 1 Conversion of the coordinate

WGS84(World Geodetic System 1984)is a coordinate system designed for GNSS satellites. Occultation observations are reported based on the WGS84geoid geometric altitude, whereas GRAPES employsthe geopotential height vertical coordinate. Therefore, we first need to convert the geometric altitudeto geopotential height for RO observations in the datapreprocessing. The conversion relation follows that of Rennie(2010). 3. 2 Quality control and data thinning

Before data assimilation, it is essential to applyquality control(QC)procedures. The most basic QCprocess is to reject the outliers that are away from themean status of most data. Three QC procedures areapplied in GRAPES. The dual weights mean and st and ard deviation method proposed by Lanzante(1996)is applied first for spatial consistency checks. This removes the refractivity observations that deviate fromthe bi-weight mean by more than five times the biweight st and ard deviation. It is performed over high(60°-90°N/S), medium(30°-60°N/S), and low(30°S-30°N)latitudinal belts for each vertical level. Thesecond QC procedure is to reject the data that couldpossibly contain super-refractions. Suspicious observations are below regions wheredN/dz< -50 N-unitkm^-1 and d²N/d²z> 100(N-unit)² km^-2 are sequentially flagged. Here, N is refractivity and z is height(Poli et al., 2009).

The background check is not performed in the pre-processing but in the data assimilation system. Thecriterion used in the background field check is

where the subscripts o and b represent the observation and background, and σ_o is the observation error st and ard deviation.

Since the vertical resolution of the occultationdata is much higher than that of the model, theymust be thinned. Our approach is to assimilate oneobservation nearest each vertical analysis level, whichnot only avoids the high correlation of observation errors, but also retains enough fine observation information.

In summary, the preprocessing procedure of ROdata is as follows. First, the data are read; then, the coordinates are converted from geometric altitudeto geopotential height, and quality control is per-formed; finally, the data are thinned in the verticaldirection. 4. Observation operator and observation error4. 1 Design of the observation operator

The observation operator of refractivity is expressed as

where N is the refractivity in N-unit, p is the pressurein hPa, T is the temperature in Kelvin, and e is thewater vapor pressure in hPa. Pressure, temperature, and water vapor pressure are not the atmospheric statevariables in GRAPES; instead, we use Exner pressure, specific humidity, and the wind components of u and v. Therefore, the calculation of Eq. (3)uses the formulaewhere c_p = 1004. 5 J kg^-1 K^-1 is the specific heat capacity of dry air at constant pressure, R_d = 287. 05J kg^-1 K^-1 is the specific gas constant for dry air, and p₀ = 1000 hPa is the reference surface pressure. Corresponding tangent linear operators areFinally, the tangent linear operator of refractivity isexpressed as

As a result of using the Charney-Phillips staggering grid, the pressure, temperature, and specific humidity are not on the same vertical level, the calculation of model refractivity is defined at the θ level. Thismeans that the error induced by the variable transformation and vertical interpolation is less, because weonly interpolate π; otherwise, we need to interpolateθ and q. The parameter π at the θ level is calculatedby assuming a linear variation of π with geopotentialheight:

Corresponding pressure at the θ level is then obtainedbased on Eq. (5).

The mean virtual temperature Tv at the θ level iscalculated by using the hydrostatic equation

where g = 9. 80665 m s^-2 is the mean gravitationalacceleration.

As GRAPES uses the terrain-following height-based coordinate, the slope of the model level mightlead to a large horizontal interpolation error. Therefore, we first calculate the background refractivity profile of the four grid points nearest the lowest tangential point of occultation observation, and then use bilinear interpolation to obtain the background refractivity profile at this tangential point. Finally, we linearly interpolate lnN to the occultation observationlocation:

where , and the subscript at-obsdenotes the observation location. We do not take thedrift of the tangential point into account. 4. 2 Test of the observation operator

To test the correctness of the refractivity observation operator, we design an idealized experiment inwhich the observation is at the analysis grid point. The variational assimilation problem is attributed tosolving the following cost function:

where vector x represents the atmospheric state, xbis the background field, vector y represents the observations, B and O are the covariance matrices of thebackground and observations, H(x) represents the observation operator, and the superscript T represents amatrix transpose.

For a single observation, the analysis increment is

where H' is the tangent linear operator. Assumingthat the observation is a logarithm of refractivity and the analysis variable is π, Eq. (15)is expressed asWhen there is only one observation, the product ofthe vectors in the two brackets of the right-h and termis a definite value, and the profile of dπ_a depends onB^πH^T =(HB^π)^T.

If the impact of moisture is ignored, that is, onlythe first part of Eq. (3)is considered, the tangentlinear operator is

Note the staggering relationship between N and π;thus, Eq. (17)should be expressed on the θ level:where h_{k, k+0+5} and h_{k, k-0+5} are the pressure coe±-cients at levels k+0. 5 and k-0. 5(observations are located at level k)Here, in Eqs. (19) and (20), the first subscript represents the location of the observation, and c_v = c_p-R_dis the specific heat capacity of dry air at constant volume. The pressure coe±cients at other levels are zeroexcept for these two levels.

Based on the above equation, we test the correctness of the observation operator in GRAPES-1DVar, which has a consistent configuration with GRAPES-3DVar. Given the following parameters: N_o = 82. 339N-unit at level k = 21(equivalent to 10357 m), σ_o =0. 415 N-unit, θ_k = 333. 9 K, π_{k+0. 5} = 0. 649031, π_{k-0. 5}= 0. 6732736, σπ = 0. 0004, z_{k+0. 5} = 10780 m, and z_{k-0. 5} = 9951 m, the numerical solution and modelsolution are found consistent(see Table 1).

4. 3 Observation error

The observation error of occultation refractivity isobtained according to the statistical result of Rennie(2010), which is a function of latitude and altitude. It is divided into three areas: low latitude(30°S-30°N), midlatitude(30°-60°N/S), and high latitude(60°-90°N/S). Figure 2 shows the st and ard deviation(σ_o)of normalized observation error of refractivity. Itis obtained by using the equation

where σ is the st and ard deviation of the normalizedobservation error, and is a function of latitude and height. 5. Impact of the RO assimilation: Experimental results5. 1 Experimental setup

Two experiments with and without RO data assimilation are performed in the global GRAPES system. They are called the control experiment and theRO experiment. In the control experiment, only theconventional observations are assimilated, which include radiosonde, synop, ships, aircraft, and clouddrift wind observations. The RO experiment thenadds the COSMIC and METOP-A RO observationsto the basic observations of the control experiment. The horizontal resolution of the outer loop is 1°, and that of the inner loop is 1. 125°, correspondingto the spherical filter's cutoff wavenumber of T106. The length of the assimilation time window is 6 h. The incremental digital filter is used for the initialization. We employed the GRAPES forecast modelversion available by the end of 2011. The modeltop is 32. 5 km with 36 layers. The horizontal resolution is 1±. The time integration step is 600 s. The physical parameterizations include the WSM6(WRF Single-Moment 6-class)microphysics parameterization(Hong and Lim, 2006), RRTM(Rapid Radiative Transfer Model)radiation parameterization(Mlawer et al., 1997), CoLM(Common L and Model)l and surface parameterization(Dai et al., 2003), MRF(Medium-Range Forecast)planetary boundary layerparameterization(Troen and Mahrt, 1986), SAS(Simplified Arakawa-Schubert)mass flux parameterizationscheme(Arakawa and Schubert, 1974), and so on. Theexperiments run from 0000 UTC 1 July 2011 to 1800UTC 31 July 2011.

Figure 3a shows an example of the horizontaldistribution of RO profile locations. There are 404COSMIC occultation profiles(dots) and 118 GRAS(Global navigation satellite system Receiver for Atmospheric Sounding)occultation profiles(triangle)fromMETOP-A satellite within each assimilation window. The occultation events mainly occur in the mid-highlatitude areas because the COSMIC and METOP-Asatellites have a high inclination. Figure 3b shows thenumber of occultation data before and after the background quality check in each assimilation time windowin the RO experiment, and also the number of RO profiles(each profile contains several data). The averagenumber of RO profiles is about 464, with 14379 ROdata at each assimilation time window, and the average number of RO data is 12828 after the backgroundcheck. About 10. 8% of the RO data are rejected.

In the one-month RO assimilation experiment, the percentage of RO observations is only about 9%, which is quite low, compared with 47% of atmosphericmotion vectors, 24% of radiosonde data, 16% of aircraft observations, etc., because the number of ROdata available in the global observing system is relatively small at present. However, the 9% occultation data surprised us in subsequent experiments, asit plays a very important role. 5. 2 Comparison with the observations

The impact of RO data on the analysis can beidentified by comparing the departure of the background and analysis to the observation(called BMO and AMO, respectively). Figure 4 shows the statistical results of the AMO and BMO profiles. Since themagnitude of refractivity from the surface to 30 km isabout 2, all the statistical computations are normalized by the observation error st and ard deviation. Thespecified observation error is reasonable, because theAMO profile is closer to zero than the BMO profile is, and the root-mean-square error(RMSE)of analysisis significantly reduced, indicating that the assimilation of occultation data has clearly a positive effect. The difference between BMO and AMO in the uppertroposphere and the stratosphere is larger, indicatingthat the model error is larger there. However, it canbe effectively controlled after the RO data are assim-ilated. The strongest effect is between 10 and 25 km, where the RO data have the highest accuracy; whilethe weaker effect is above 25 km and below 10 km. Besides the data accuracy; the error might be relatedto the vertical structure of GRAPES, especially nearthe model top.

5. 3 Comparison with NCEP analysis

The NCEP analysis data are used to evaluatethe impact of RO observation assimilation. GRAPESis a developing NWP system, in which the observation types assimilated are still limited; whereas theNCEP analysis data accommodate assimilation of avariety of satellite-remote sensing data, making it avaluable dataset to compare and examine the problems in GRAPES. Figure 5 presents comparisons between GRAPES and NCEP analysis. The geopotentialheight difference between the control experiment and the NCEP analysis is mainly located in the SouthernHemisphere, the stratosphere of the tropics, and thehigh latitudes in the Northern Hemisphere. These differences are clearly reduced after the RO observationsare assimilated. The difference between the RO experiment and the control experiment is similar to thatbetween the control experiment and the NCEP analysis, especially over 40°-70°S, where the number ofradiosonde and aircraft observations is only 5% of thedata in the same latitude zone in the Northern Hemisphere. As a result, once the RO observations areassimilated, a significant improvement in the analysis can be seen. The situation with respect to temperature analysis is similar to that of geopotentialheight. However, the gradient difference above 20 hPais larger, which might be related to the vertical structure of GRAPES, but this needs a further investigation.

For the wind analysis, the differences with and without the assimilated RO data are concentrated inthe tropics as well as in the Southern Hemisphere(Fig. 5d). The improvement in the wind analysis with ROobservation assimilation in the equatorial tropics isnot so obvious as that in the Southern Hemisphere(Fig. 5e). This illustrates that the wind analysis ofGRAPES in the tropical region is not good. Giventhe fact that the coupling between wind and pressurefields is weak in the tropics, and the RO data are inessence a type of thermal temperature observation, itis acceptable that the addition of RO observations exerted less of an impact on the wind analysis, and theimpact might have been offset by some systematic error if there is any.

Figure 6 shows the time evolution of RMSE at 10hPa for the control and the RO experiment. The difference between the control and the RO experiment inthe Northern Hemisphere is not evident, with a magnitude of 30 gpm, due to assimilation of considerableamount of conventional observations. Since there areonly a few conventional observations in the SouthernHemisphere, the RMSE of the control experiment is30 gpm on the first day, and grows to 320 gpm onthe 23rd day. This phenomenon is common above 100hPa; and the higher the model level, the more obviousthe phenomenon, reflecting the fact that the model error cannot be controlled in the Southern Hemispherewithout satellite data(such as RO)being assimilated. The same situation is observed for the wind and temperature fields(figures omitted).

5. 4 Impact of RO assimilation on forecast

Figure 7 shows the anomaly correlation coefficients(ACCs)for the 500-hPa geopotential height calculated for the GRAPES 8-day forecast as a functionof time. Regardless of geographical location, a steadyimprovement of the RO experiment in ACCs is evidentcompared to the control run. The largest improvementis in the Southern Hemisphere, where the largest improvement rate is observed on the fourth day and theforecast validity range is increased by approximately0. 8 day. The contrast in ACCs between the RO and the control experiment in the Northern Hemisphere isnot so obvious as that in the Southern Hemisphere, because of the high density of radiosonde observationsin the Northern Hemisphere. The ACCs of the ROexperiment are even slightly lower than those of thecontrol experiment in East Asia beyond 96 h of forecast.

6. Contribution of RO assimilation to the stability of GRAPES assimilation and forecast cycling system

After completing the one-month RO data assimilation experiment, rather than stop the experimentwe let it continue for another 8. 5 months. Therefore, the actual period of the RO experiment is from 0000UTC 1 July 2011 to 1800 UTC 16 April 2012*. Ourpurpose is not only to evaluate the impact of the ROobservation assimilation, but also to further evaluatethe performance of the newly developed GRAPES-Varsystem.

As an example, Fig. 8 presents the evolution ofthe mean difference and RMSE between the GRAPES and NCEP geopotential height analysis in the RO experiment at 500 hPa for both hemispheres. Figure 8illustrates that the number and distribution of the observations are important factors that affect the analysis quality of the data assimilation system. The number of radiosonde stations is greater in the NorthernHemisphere, and more observations are in the troposphere so the RMSE and bias are general smallerin the Northern than in the Southern Hemisphere. Furthermore, the RO observations gradually dominate in the upper troposphere and stratosphere forboth hemispheres; therefore, the difference with theNCEP analysis for the Northern and Southern Hemispheres gradually reduces with height(figure omitted). The analysis skills in both the Northern and Southern Hemisphers are still relatively comparable evenafter nine-month integration, especially in the stratosphere where the number of conventional observationsdecreases and RO observations dominate with a uniform global coverage.

Figures 9 and 10 show the monthly geopotentialheight and wind pattern of the GRAPES and NCEPanalysis at 500 and 20 hPa in the 9th month. In arelatively simple data assimilation system without alarge amount of satellite remote-sensing observationdata assimilated, the system can still run stably formore than nine months, and the difference with NCEPanalysis is fixed at a stable magnitude all the time. The similarities between the GRAPES and NCEPanalysis illustrate that the RO data have become oneof the basic components of GRAPES-VAR. The ROdata not only improve the short-term analysis and forecast, but also guarantee the running stability ofthe GRAPES data assimilation and forecast cyclingsystem.

7. Conclusions and discussion

The GNSS RO has become an important component of the global observing system, and possessesgreat value in global data assimilation. The RO dataassimilation experiment with GRAPES has shownthat, although the number of occultation data is small, no more than 10% of the total observations in the assimilation experiment, they still can greatly improvethe analysis quality and forecast ability. Furthermore, assimilation of the RO data plays a particularly important role in supporting the stable running of theGRAPES assimilation and forecast cycling system.

Although the experiment system is relatively simple, in which only the conventional and RO observations are assimilated, it is easy to identify the system'sproblems. For example, it is clear that the analysis isworse in the lower troposphere and near the modeltop. On the one h and , this might be related to thespecification of observation error; on the other h and , it is possible that the problem lies in the model's vertical resolution. The tuning of the observation error and other prescribed parameters such as the criterionof background field check, is still needed in a real operational assimilation system with many more satelliteremote-sensing data assimilated, in order not to amplify the impact of the RO data.

Currently, the RO observation operator inGRAPES is a local refractivity operator, which is simple and has a low computational cost. Bending angleis a rawer product than refractivity. The advantageof bending angle assimilation is that it considers eachRO observation value that represents the atmosphericcumulative effect in the path of the radio wave propagation, and errors introduced by the spherical symmetry assumption are smaller and the observation erroris simpler. Therefore, we need to develop the bendingangle operator with the rise of the GRAPES modeltop. Additionally, the spherical symmetry assumptionis not valid where the horizontal gradient of refractivity is large; therefore, a non-local refractivity ora ray-tracing forward operator should be taken intoaccount.