Analysis of vertical wave number spectrum of atmospheric gravity waves in the stratosphere using COSMIC GPS radio occultation data( discussion paper )

GPS radio occultation (RO) is characterized by high accuracy and excellent height resolution, which has great advantages in analyzing atmospheric structures including small-scale vertical ﬂuctuations. The vertical resolution of the geometrical optics (GO) method in the stratosphere is about 1.5 km due to Fresnel radius limitations, but full 5 spectrum inversion (FSI) can provide superior resolutions. We applied FSI to COS-MIC GPS-RO proﬁles from ground level up to 30 km altitude, although basic retrieval at UCAR/CDAAC sets the sewing height from GO to FSI below the tropopause. We validated FSI temperature proﬁles with routine high-resolution radiosonde data in Malaysia and North America collected within 400 km and about 30 min of the GPS RO events. 10 The average discrepancy at 10–30 km altitude was less than 0.5 K, and the bias was equivalent with the GO results. Using the FSI results, we analyzed the vertical wave number spectrum of normalized temperature ﬂuctuations in the stratosphere at 20–30 km altitude, which exhibits good consistency with the model spectra of saturated gravity waves. We investigated


Introduction
The Global Positioning System (GPS) radio occultation (GPS-RO) technique is an active limb sounding observation of the Earth's atmosphere and ionosphere using a GPS Introduction  (Ware et al., 1996), was successfully conducted by the University Corporation for Atmospheric Research (UCAR) in collaboration with the Jet Propulsion Laboratory (JPL) in 1995 using a micro-satellite, MicroLab-1. Since this pioneering experiment a number of GPS-RO missions have been carried out, such as the Challenging Minisatellite Payload (CHAMP) (Wickert et al., 2001) (Anthes et al., 2008) was launched as a joint project between UCAR and the National Space Organization (NSPO) of Taiwan. The GPS-RO has become an epoch making measurement technique, greatly 10 impacting diverse applications in the atmospheric sciences, meteorology, boundary layer studies, numerical weather prediction, global climate change, and ionospheric disturbances.
With a single LEO satellite, about 250 GPS-RO events can be obtained per day using a single occultation antenna. The COSMIC mission employs a constellation of six LEO 15 satellites having two forward and backward antennas on each satellite, for a total of about 2500 GPS-RO data points. Because the GPS RO is a limb sounding technique, it has an inherent advantage in height resolution. Further, by applying advanced radio holographic (RH) methods, vertical resolution was as good as a few hundred meters in the lower atmosphere. Such a good height resolution, comparable to a ground-based 20 radiosonde, has not been attained with other conventional satellite techniques.
Most current retrieval techniques for GPS RO assume spherical symmetry of the atmosphere in their inversion process, so the horizontal resolution of a retrieved profile centered at a tangent point is roughly 300 km. Note also that the tangent point of the ray drifts horizontally. Knowledge on spatial variations of the atmospheric field can be down to the upper troposphere, and then the retrieval procedure is switched to FSI from there to the ground.
In this study, we primarily follow the software system at CDAAC, but with modification of the height range for FSI analysis so as to obtain a high-resolution temperature profile in the stratosphere. We validate the accuracy and vertical resolution of FSI retrieval in the stratosphere through comparison with radiosondes. Taking advantage of the increased height resolution, we study the vertical wave number spectra of temperature perturbations at 20-30 km, and discuss consistency with theoretical predictions of saturated gravity waves. We also analyze the valid range of the FSI retrievals. 15

Geometrical Optics (GO)
Here we briefly introduce the fundamentals of the GPS RO retrieval method. The necessary information for retrieval is the amplitude and phase of the occulted GPS radio signals. Based on the precise orbit information of the GPS and LEO satellite, the atmospheric excess phase (delay) is calculated. Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | this hypothesis is a good assumption. But this also implies that the GO method is not valid at shadow boundaries and caustics due to their non-monotonic refractive gradients (e.g., Melbourne, 2004). The first retrieved result is a bending angle profile of the L1 and L2 signals. The received signal is smoothed, with the normal filter for the L1 signal corresponding to 5 a height smoothing of about 0.8-1 km. Because of its poorer quality, however, the L2 signal may not provide a bending angle profile with the same precision and resolution as the L1 signal. A stronger smoothing window is therefore sometimes applied to reduce the effect of L2 noise, giving a height smearing over 3-4 km. Ionospheric calibration then eliminates the ionospheric effect using the L1 and L2 bending angles.
Note that small-scale ionosphere perturbations may still remain as retrieval errors after ionospheric calibration.
In altitudes higher than about 40 km, the magnitude of the neutral atmospheric bending angle becomes very small due to its exponential decrease, sometimes making it impossible to extract a useful signal because the magnitude of receiver system noise 15 remains constant. The uncertainty of the signal at high altitude could introduce a bias that propagates toward lower altitudes through Abel inversion. To reduce this effect, the observed bending angle is subjected to optimal estimation by referring to a bending angle from climate models.
An Abel inversion transforms the optimized bending angle profile to a refractive in-20 dex profile under the basic hypothesis that the atmosphere is spherically symmetric. Pressure and temperature profiles can be derived from the refractive index. Note that in a moist atmosphere the refractive index cannot be uniquely converted to a temperature, so the contribution of water vapor to the refractive index must be separated out by using, for example, a variational method with a numerical weather forecast as its initial 25 value.
The available height range of the GPS-RO profiles is normally below about 40 km. Above 40 km, ionospheric effects and other measurement errors could exceed the atmospheric phase delay, which decreases exponentially with height. On the other hand, 2075 at the lowest part of the troposphere multipath effects become a serious problem by a sharp height gradient of the refractive index due to water vapor. Moreover, below the tropopause the horizontally irregular distribution of water vapor may invalidate the hypothesis of a spherically symmetric atmosphere, as is assumed at higher altitudes. This will also introduce errors to retrieved results. 5 In the interval between 5-30 km, accuracy in the temperature inversion can be expected to be 0. 2-0.4 K (Kursinski et al., 1997). Because the GO method describes electromagnetic waves as optical rays, vertical resolution is limited by Fresnel diffraction. The vertical width of the first Fresnel zone is on average about 1.5 km in the stratosphere, and roughly 0.5 km in the lower troposphere (e.g., Melbourne, 2004).

Full Spectral Inversion (FSI)
Although GO has reasonably good vertical resolution, it is limited by the Fresnel zone. Moreover, GO cannot overcome the multipath problem in the complex atmospheric structures frequently caused by water vapor in the lower troposphere. To solve the multipath problem and obtain a better vertical resolution, radio holographic (RH) meth- 15 ods are introduced to the GPS RO retrieval. The RH methods are types of wave optical methods that invert the occulted radio signal in a complex form. At CDAAC several RH methods are used as an inversion algorithm, including back propagation (BP) Karayel et al., 1997), sliding spectrum (SS) (Sokolovskiy, 2001), canonical transform (CT) (Gorbunov, 2002), and full spectrum inversion (FSI) (Jensen 20 et al., 2003). Among these RH methods FSI provides the same accuracy and resolution to BP and CT. Moreover FFT can also be used by assuming that the GPS and LEO satellites are in a circular orbit, considerably reducing computational cost. The FSI method is therefore widely used.
In the standard retrieval program at CDAAC, the FSI method only retrieves the L1C 25 signal below some sewing altitude. The bending angle of L2 and the upper part of the bending angle profile is inverted using GO. The final profile is constructed by combining (sewing) the RH-based and GO-based profiles at an appropriate altitude ( Sokolovskiy and C. Rocken, personal communication, 2006). The sewing altitude of FSI-based and GO based profiles depends on strong fluctuation and abrupt noise increases in L2 Doppler, or large deviations between L1 and L2 Doppler. Because CDAAC focuses on reducing multipath effects, especially below the tropopause, the sewing points are usually defined between 10 and 20 km. However, atmospheric 5 profiles with better vertical resolution are sometimes required for scientific and practical applications in meteorology, the atmospheric sciences, and so on. In this study, we use a different sewing scheme from CDAAC. Specifically, in our version the sewing altitude is much higher (30 km). Note that features such as optimization of the bending angle using a model profile, Abel inversion, and further derivation of atmospheric parameters 10 are the same as with their counterparts in the GO inversion scheme. Figure 1 shows a flowchart of the retrieval procedure used in this study, which is modified from that developed at CDAAC. L1 and L2 excess phases are the input of the GO inversion block. The L1 and L2 Doppler shifts are computed for a sewing criterion that outputs to the Sewing block, which combines the GO and RH based profile. Another 15 output is L1 and L2 GO-based bending angles, which help calibrate the ionospheric effect in the combined bending angle profile. Note that because the RH method is only used below the sewing point, where the L2 signal is usually considered unavailable, ionospheric calibration for the RH-based profile is computed from extrapolation of the L2 bending angles. Ionospheric calibration in higher altitudes is calculated in the same 20 manner as in the CDAAC version, while below 30 km the RH-based profile is corrected using the GO-based L1 and L2 bending angles.

Validation of GPS RO temperature profiles with radiosondes
When the sewing height is elevated, it may introduce noise with short vertical scales from several error sources, such as ionospheric perturbations or receiver system noise, 25 and such noise can contaminate the atmospheric excess phase. It is noteworthy that the atmospheric excess phase decreases exponentially with altitude, and thus noise 2077 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | effects become more serious at higher altitudes. Consequently, validation of retrieved results is particularly necessary at higher altitudes.
Here, we verify the FSI-based temperature profiles through comparison with radiosonde results with good height resolution. We obtained 54 radiosonde profiles launched at Kuching (1.5 • N, 110.4 • E) by the Malaysian Meteorological Department. 5 We analyzed the original data records of VAISALA radiosondes (type RS-80), and obtained a temperature profile with a height resolution of about 100 m. Figure 2  Malaysia. The temperature profiles are shown in Fig. 3. The solid lines in Fig. 3 are the FSI retrieval for #49 and #50 in our study, and the dashed lines are taken from CDAAC and analyzed with the GO method. Since the effect of humidity on the refractive index becomes negligible above the tropopause, located at about 16 km, we use the dry temperature of GPS-RO. 15 In Fig. 3 the profiles above the tropopause exhibit small-scale temperature perturbations in the FSI and radiosonde profiles, while the GO profile shows a rather smooth height structure. A sharp change in the temperature gradient is clearly seen near the tropopause. In the stratosphere wave-like structures are evident at 18-20 km altitude, showing a downward phase progression from #50 to Kuching and toward #49. This 20 can be interpreted as a manifestation of an atmospheric gravity wave that propagated northward passing over these sites. This example provides proof that the FSI method is sufficient to analyze temperature perturbations with vertical resolution comparable with a radiosonde, even in the stratosphere.
Next the radiosonde launch sites were less than 400 km and with time separation less than 30 min. The number of compared pairs ranges from 210 to about 350, depending on the month and altitude range. Standard deviations and the mean temperature differences (GPS-RO minus radiosonde) are shown in Fig. 4 for the GO and FSI methods in November 2006. The 5 mean absolute temperature difference at 10-30 km for FSI is less than 0.5 K for all three months. Note that the results are not shown for October and December, as the tendency is very similar to those in Fig. 4. There is no obvious temperature bias for either FSI or GO, and the accuracy of the FSI method seems to be similar to that of GO. The average standard deviation for FSI does not vary significantly with altitude. At 10 high altitudes, where the number of comparisons is similar between the FSI and GO profiles, the standard deviation and the mean are very similar. At CDAAC most of the GO profiles are replaced with the FSI results below about 16 km, so as the number of GO profiles decreases the number of comparisons below 16 km also decreases. 15 Here, we analyze vertical wave number spectra of temperature fluctuations, and compare them with a theoretical model and radiosonde observations. Through statistical investigation of the spectra, we aim at identifying the vertical resolution of FSI retrievals, and their dependence on the altitude region in the stratosphere.

Vertical wave number spectrum
First, we briefly introduce the fundamental concept of the vertical wave number spec-20 trum of saturated gravity waves (Smith et al., 1987). Earlier studies reported that temperature profiles in the upper troposphere and the stratosphere fluctuate with vertical scales ranging from tens of meters to several kilometers. These fluctuations are believed to be superpositions of various atmospheric gravity waves with different wave lengths (e.g., Dewan, 1979;VanZandt, 1982). Spectral analysis is adopted to express Introduction regardless of the altitude range and geographical location (VanZandt, 1982). Observational results indicate that the asymptotic logarithmic slope of the vertical wave number spectra is close to −3 for large wave numbers, while the spectral density is limited by the saturation values determined by the background Brunt-Vaisala frequency (Smith et al., 1985;Fritts et al., 1988;Tsuda et al., 1989). A formula for the saturated grav-5 ity wave spectrum for the normalized temperature perturbations can be expressed as follows for a large wave number m, where N is the Brunt-Vaisala frequency, g is the acceleration of gravity and T and T are the temperature perturbation and the background temperature, respectively (Tsuda 10 et al., 1991). This relation was validated by intensive radiosonde observations (Tsuda et al., 1994). It is noteworthy that Eq.
(1) provides the maximum of the spectral density for saturated gravity waves because of wave dissipation due to convective instability (e.g., Fritts and Alexander, 2003). The vertical wave number spectra were analyzed using GPS RO temperature profiles 15 from GPS/MET, which demonstrated similar characteristics of the asymptotic logarithmic slope with a model spectrum of saturated gravity waves (Steiner and Kirchengast, 2000;Tsuda and Hocke, 2002). These results are based on GO retrieval, however, which may not fully resolve small vertical scales. In this study we analyze the vertical number spectra of temperature fluctuations using FSI results. 20 Figure 5 shows the spectra from GPS RO with both GO and FSI methods compared to the theoretical model. A spectrum from the US radiosonde data is also plotted. The Brunt-Vaisala frequency N is computed from the mean temperature and pressure profiles with radiosondes. (Note that the unit for m in Fig. 5 is cycles m −1 .) Spectra from both radiosonde and FSI in Fig. 5 are approximately parallel with the model except for 25 the first few spectral points where gravity waves are not saturated, and therefore the spectrum does not follow Eq. (1). The FSI spectrum closely matches that from radiosondes at wave numbers as large as about 2 × 10 −3 cycles m −1 . Beyond 2 × 10 −3 cycles m −1 , the FSI spectrum tends to show a gradual (flat) slope, and the spectral density becomes larger than the model values. This indicates the appearance of a noise floor for large m. The spectra retrieved by GO shows a much smaller spectral density than the model, particularly for large m. 5 This is caused by smoothing with a filter in the inversion program, which is adjusted to a height resolution of about 1.5 km by considering the Fresnel radius in the stratosphere. For statistical validation, we calculated the spectra by averaging over all December 2006 COSMIC GPS RO profiles at 20-30 km in the latitude ranges 20 • -40 • S and 20 • -40 • N. The results using both FSI and GO are shown in Fig. 6. The spectra from these two latitude bands show a very similar shape. The number of averaged profiles is much larger than in Fig. 5, making the spectra much smoother, but the noise floor does not disappear.
The noise floor is related to the vertical resolution of the FSI method. Because a bend from −3 to 0 (white noise) in the logarithmic spectral slope appears at around 15 3 × 10 −3 cycles m −1 , we could resolve waves with a vertical wave length of about 300 m. The vertical resolution of the FSI retrieval seems to be about 150 m, i.e., half the shortest detectable wave length.
Earlier studies estimated the vertical resolution of RH methods, including FSI, using numerical simulations. Mortensen et al. (1999) showed the resolution to be about 20 100 m, while Gorbunov et al. (2004) indicated a better resolution of about 60 m. Those studies, however, assumed an ideal case without additional errors such as ionospheric noise or horizontal inhomogeneities, thus giving inferred vertical resolutions better than our estimate based on observed results.
It is important to note that vertical resolution can also vary with altitude. Figure 7a  respectively. The noise floor magnitude increases with altitude, and is most evident at 23-30 km in Fig. 7. This suggests that the noise floor in the spectra at 20-30 km (see 2081  Fig. 6) could be mainly attributable to the high altitude region.

Upper limit of the FSI profiles
The atmospheric excess phase decreases exponentially with altitude, as it is basically related to atmospheric density. While system noise may not depend on altitude, moreover, ionospheric effects increase at higher altitudes. The quality of the FSI profile, 5 therefore, tends to degrade at higher altitudes. The main purpose here is to define the maximum available height of FSI profiles for studies of small vertical scale geophysical perturbations such as atmospheric gravity waves. We utilize the behavior of the noise floor in the vertical wave number spectra as an index of deterioration of the FSI retrieval. We analyze the wave number spectra for 10 individual GPS RO profiles using the same procedure as in Figs. 5-7, but shorten the height range for spectral analysis from 7-10 km to 2 km and concentrate on the behavior of the spectra at large (>10 −3 cycles m −1 ) wave numbers. If the spectral density noticeably exceeds the model value given in Eq. (1), then origin of the temperature perturbations are most likely not gravity waves, but they are affected by other noise. 15 We frequently found that the spectra at 26-28 km and 28-30 km altitudes show noisy behavior, with the spectral density obviously exceeding theoretical predictions.
For the Fourier Integral Operator (FIO) method including FSI, Gorbunov et al. (2006) proposed a technique that can estimate signal tracking errors for bending angles. This technique is based on the width of the running spectra of the transformed wave field 20 multiplied by a reference signal. The reference signal is computed by smoothing the phase of the transformed field. In this method, however, the smoothing parameter should match the vertical scale of the noise. In our case, noise may have various wave lengths in different profiles, and thus it is difficult to use a constant smoothing parameter for all measurements. 25 Lohmann (2006)  work on each profile without external data other than a signal-to-noise ratio (SNR) threshold. In the standard inversion program at CDAAC, estimated errors of the bending angle are calculated using this method. But this estimation does not distinguish perturbations with different wave lengths.
Here we propose a simple procedure concentrating on the fluctuation with different 5 wave lengths in the bending angle profiles, and then determine the upper limit for the FSI profile. Considering saturated gravity wave theory, the normalized temperature perturbations T /T 0 should not exceed a certain threshold. Similarly, the ratio of the bending angle perturbation α to α o , the optimized bending angle with GO, can be limited by a threshold. Provided the bending angle perturbation caused by system 10 noise remains the same magnitude throughout an occultation event, the ratio α /α o exponentially increases with altitude. We therefore assume that the height variation of the ratio α /α o can be used as an index to judge contamination by noise. When noise dominates the ratio increases exponentially with altitude, but it remains nearly constant when the signal instead comes from actual atmospheric perturbations.
15 Figure 8 shows an example of α /α o versus height for the profile C001.2006.335.00.50.G23. We consider two vertical scale ranges, 0.1-0.5 km and 1-2 km, and estimate α from the standard deviation of bending angle perturbations at each scale range. The ratio for short vertical scales rapidly increases above about 26 km, but it is less than about 1% for longer scales in the entire height 20 range. For this particular case, we should limit the maximum height to 26 km to study small-scale atmospheric perturbations.
According to Jensen et al. (2003), the accuracy of FSI is 1 K. We take this value to infer a threshold. Considering the standard atmospheric model, the maximum bending angle fluctuation corresponding to the temperature perturbation with 1 K amplitude and 25 vertical wave length of 0.5 km becomes about 2.1%. For the estimate of standard deviation, we take 1.4% as the root mean square value of the sinusoidal temperature perturbations. We statistically analyzed the distribution of the upper boundary in terms of local time at latitudes 20 • S-20 • N, where ionospheric noises seem to be larger. We show the results for November and December 2006 in Fig. 9a and b, respectively, as the amount of GPS RO data fluctuates depending on local time. We only focus on the local hours where numerous profiles existed, specifically ranges 05:00-16:00 LT (local time) 5 in November, and 09:00-14:00 and 21:00-3:00 LT in December. Note that longitudinal asymmetry of the ionosphere is greatest at sunrise and sunset. There is little difference between the two cases, however, suggesting that ionospheric effects are not significant in determining the upper limit of FSI retrieval. To summarize, Fig. 9 indicates that FSI profiles can be used up to about 28 km when studying temperature perturbations with 10 vertical wave lengths as short as about 0.5 km.

Concluding remarks
Because of its advantages in vertical resolution, accuracy, and self-calibration, the GPS RO technique is quite promising for areas such as the atmospheric and ionospheric sciences, global and mesoscale weather forecasting, climate change monitor- 15 ing, and boundary layer research. Among the various retrieval algorithms for GPS-RO, the GO and FSI methods are commonly used. Because of limitations due to the Fresnel radius for GO, its vertical resolution is limited to about 1.5 km in the lower stratosphere, while FSI is not affected by such limitations. We applied FSI up to 30 km on COSMIC data and obtained temperature profiles with 20 a better vertical resolution than those published by CDAAC with GO. We compared the FSI profiles with radiosonde obtained in Malaysia, which showed good consistency with perturbations caused by atmospheric waves. Further comparisons were done with North American radiosonde data from October to December 2006. Statistical tests indicate that the temperature errors for FSI are less than 0.5 K, and the temperature 25 bias is equivalent to GO retrieval. We studied the vertical resolution of small-scale temperature fluctuations using a vertical wave number spectrum. Spectral analysis showed that a noise floor tends to appear for wave numbers larger than about 2-3 × 10 −3 cycles m −1 . The vertical resolution of the FSI method in the lower stratosphere is inferred to be about 100-200 m.
Since the atmospheric excess phase decreases exponentially with altitude, various forms of noise can easily contaminate GPS-RO retrieved data for high altitudes. A case study also showed noisy fluctuations in the bending angle profile above about 26 km altitude. We aimed at defining the upper limit for FSI profiles. By separating perturbations of the bending angle at different wave lengths, we found that small-scale 10 perturbations are more sensitive to noise. We thus used a normalized perturbation of the bending angle as an index to determine the upper boundary. Assuming a temperature precision of 1 K, the upper limit for the available profiles is estimated to be about 28 km. Distribution of the upper boundary does not seem to be affected by ionospheric conditions. 15 Although FSI retrieval provides a height resolution in the lower stratosphere (below about 28 km) superior to other satellite measurements, horizontal resolution is not significantly improved. GPS RO therefore observes only part of the horizontal wave number spectra of temperature perturbations, and so the effects of an observational filter must be considered in interpreting the analyzed characteristics of atmospheric

2097
A Self-archived copy in Kyoto University Research Information Repository https://repository.kulib.kyoto-u.ac.jp