Parameters of internal gravity waves in the mesosphere-lower thermosphere region derived from meteor radar wind measurements

A procedure of revealing parameters of internal gravity waves from meteor radar wind measurements is presented. The method is based on dividing the measuring volume into different parts and, using wavelet analysis, calculating the phase progression of frequency peaks in the vertical and horizontal direction. Thus, the distribution of vertical and horizontal wavelengths and directions of IGW energy propagation, using meteor radar data, has been obtained. The method was applied to a 4-month data set obtained in July and August, 1998 and 1999. As expected, the majority of waves have been found to propagate upwards, although a considerable number seem to propagate downwards as well. High-frequency (intrinsic periods T* of less than 2 h) waves are dominating. The distribution of waves over the course of an average day is only weakly structured, with weak maxima in the morning and evening.


Introduction
Studies of the wind regime in the mesosphere/lower thermosphere (MLT) regions have shown the essential influence of internal gravity waves (IGW) on dynamic processes and energy transport and dissipation in these height layers (e.g.McLandress, 1998).Propagating from below, but partially also originating in the middle and upper atmosphere, they can produce turbulence and substantially deposit momentum and energy, and they can influence the general circulation, thermal regime, and composition of the middle and upper atmosphere (Gavrilov et al., 2000(Gavrilov et al., , 2002)).
During the last decades, radar measurements have successfully been used to provide information about IGW.These measurements have been carried out, e.g., using the Japanese Correspondence to: Ch.Jacobi (jacobi@uni-leipzig.de) middle and upper atmosphere radar at Shigaraki (Yamamoto et al., 1987;Gavrilov et al., 1996Gavrilov et al., , 1997Gavrilov et al., , 2002)), medium frequency (MF) radar (Vincent and Reid, 1983;Manson and Meek, 1988;Manson et al., 1999Manson et al., , 2003;;Gavrilov et al., 1995), and even the low frequency (LF) ionospheric E-region drift measuring method (Gavrilov and Jacobi, 2004).In addition, meteor radar measurements have been used to derive IGW information (Kalchenko et al., 1985;Xiong et al., 2003;Manson et al., 2004).However, thus far, detailed parameters of IGW, such as horizontal wavelength and the energy propagation direction, have not been available from typical meteor radars, only from more sophisticated systems, such as Mesosphere-Stratosphere-Troposphere (MST) radars with a controllable beam direction, and large MF radars have provided those parameters (Vincent and Reid, 1983;Manson and Meek, 1988).
In this paper results of a study carried out using an automatic goniometer on the meteor radar at the Kharkiv National University of Radioengineering (49 • 30 N, 36 • 51 E) are presented, with the purpose to obtain the individual time-spatial parameters of IGW and to determine features of their statistical distribution.The automatic goniometer of the Kharkiv meteor radar station (AG MRS) is intended to study dynamic parameters of the Earth's atmosphere at altitudes 80-105 km with the opportunity of revealing the high-altitude structure of zonal wind motions.AG MRS, operating on 31.1 MHz, consists of a five-element antenna array that is directed to the east, and makes Doppler measurements of the radial drift velocity, angular data (azimuth and elevation), and range, i.e. slant distance up to the reflecting area of meteor trails (Kashcheyev and Oleynikov, 1994).

Technique of data processing
The algorithm to process the meteor radar data to obtain information about IGW parameters is based on the one described by Gavrilov (1981) and Kalov and Gavrilov (1985).The volume that is surveyed by the meteor radar is divided   into several regions of 6-km height to ensure a sufficient number of meteors in each sub-volume.These sub-volumes are constructed for the entire height layer under consideration (77-103 km, i.e. centres of bins from 80-100 km), with a height shift of 2 km, i.e. with a 4-km overlap.Similarly, in the horizontal direction the zone between the 80-and 220km distance is divided into sub-volumes of d i =60 km (in some cases with a small meteor number, d i =80 km) width (in the viewing direction of the radar), with a horizontal step of 10 km or 20 km.We discard waves with a horizontal (vertical) wavelength of more than λ h =1500 (λ z =150) km.The construction of the sub-volumes is outlined in Fig. 1.The typical total meteor count per day is 5500-6500.The separate altitude/distance bins contain between 100-1000 meteor echoes per day, with an approximate average of 500-700 echoes per bin.These numbers are dependent on the position, i.e. near the centre of the measuring volume (height near 90 km, distance near 150 km) more meteor echoes are detected in comparison with the edge bins.We use horizontal zonal wind velocities, which were calculated from projecting the horizontal line-of-sight winds on the east/west direction.This procedure gives real estimates of the zonal winds, if the meteor distribution is symmetric with respect to the east-west direction.Therefore, we use only the zonal component both for wind estimation, and for gravity wave estimation.This means that the horizontal wavelengths are not uniquely determined, depending on the north-south component.This has to be taken into account when interpreting the data, assuming an isotropic wave spectrum this would mean that the most probable horizontal wavelengths and phase speeds are overestimated by a factor ∼ √ 2 (Vincent and Reid, 1983).The winds have been bandpass filtered before further processing.Useful results can be obtained with a very few (>20) meteor contribution to an average wind value.
To detect IWG, an algorithm that includes a wavelet analysis to obtain the amplitude and phase spectra is used.A sliding window low-pass filter is applied to suppress random errors in the high frequency component of the meteor radar   wind time series in each sub-volume.The size of the sliding time window is 10 min, so that the frequency components with periods less than 5 min will be suppressed.From the time series 5-min averages are calculated; this allows detecting IGW with periods of more than about 30 min.A continuous wavelet transform with complex 8-order Gaussian or complex Morlet wavelet is carried out.These complex mother wavelets have been found to be the most suitable for IGW modelling.An example of fitting the wavelet to the time series is shown in Fig. 2. Using the horizontal wind velocity wavelet spectra for separate range zones d i , an average wavelet spectrum is computed for each altitude layer h i. , and times and frequencies of wavelet spectrum maxima are determined.Only data are accepted, where the horizontal phase change of the wavelets satisfied a linear dependence.This procedure was repeated for each height bin.As a criterion for accepting the data, the horizontal wavelengths should not deviate more than their standard deviation from the mean wavelength.
The next step in processing the data is a robust leastsquares fitting for detecting phase dependences.It is usually assumed that errors can be described by a normal distribution, and that extreme values are rare.However, a simple least-squares fitting is sensitive to outliers, because squaring the residuals magnifies the effects of these extreme data points.To minimize the influence of outliers, we use iteratively a reweighed least-squares robust fitting with bisquare weights.Robust fitting allows for a reduction in the miscalculation of the horizontal and vertical lengths of the gravity waves with outliers in a data set (Holland and Welsch, 1977;Huber, 1981;Street et al., 1988;DuMouchel and O'Brien, 1989).
The range-phase dependence of the IGW is constructed from the spectra in each range and height zone using a robust least-squares fitting, with a check of the hypothesis of a linear regression for the range-phase dependence data (Brandt, 1975).The horizontal wavelength and phase speed of the   IGW is then calculated.The altitude area, in which the calculated horizontal wavelengths do not deviate more than the r.m.s.deviation from their mean value, is accepted as the spatial area where the wave is present.A vertical IGW structure check is applied similarly.For each altitude layer h i the mean values of phase and altitude are calculated, from which the vertical phase dependence is constructed, again using a robust least-squares fitting and the hypothesis of linear regression.Finally, the dispersion relation is checked, i.e. it is tested whether the IGW phase velocity satisfies the following expression (Gavrilov and Medvedev, 1997): with c is the horizontal phase speed, N is the Brunt-Väisälä frequency, and H is the scale height.Taking H =6000 m and N 2 =5.5•10 −4 s −2 gives an upper limit of 280 ms −1 for the phase speed.An initial number of 1716 gravity waves was measured, 14 of these have been rejected through the dispersion relation test, 50 waves were rejected, because their horizontal wavelength was longer than 1500 km, while 21 waves were rejected because their vertical wavelength was longer than 150 km; these maximum wavelength values were chosen as upper wavelength limits.Additionally, we rejected 107 waves or 224 waves, because their wave periods or intrinsic wave periods, respectively, were longer than 5 or 6 h.In total, 416 waves were rejected, and 1300 waves were accepted.

Results
The above-described algorithm was applied to the meteor radar wind data that was obtained during 2 summer campaigns in July-August 1998 and 1999.Background winds and long-period waves during the 1998 time interval have been presented in Jacobi et al. (2001).Height-time crosssections of the background zonal prevailing winds are shown in Fig. 3.Note that the profiles have been calculated using a sliding 4-day window, so that most of the quasi two-day wave signal (Jacobi et al., 2001) has been averaged out.Using data of these 124 days, 1300 waves and their time-spatial parameters are determined.
Histograms of vertical and horizontal IGW wavelength are shown in Fig. 4. For convenience, positive (negative) wavelength values correspond to upward (downward) or eastward (westward) propagation of IGW energy, respectively.The figure has been calculated using all data between 80 and 100 km, while division of the height range into subranges did not alter the structure of the histograms.On average, the prevailing wind changes with height, from westward winds in the mesosphere to eastward winds in the lower thermosphere (see Fig. 3).From this one would expect a change in the preferred propagation direction towards the west.The fact that we cannot see this may indicate that either our method is not sensitive enough, or the disturbance of the mean profile through waves (e.g. the strong quasi-two-day wave, see Jacobi et al., 2001) and tides is so strong that the simple division into two layers does not show the effect clearly.The most probable values of the vertical wavelength λ z are 10-30 km, for any direction of IGW propagation.Upward propagating waves are more frequent (65% of all waves) than downward propagating (35%), which is in accordance with theory (Andrews et al., 1987).The dominant horizontal wavelengths lie between λ h =100-300 km.This is in accordance with the overall results by Manson and Meek (1988), however, their summer data tend to show somewhat shorter wavelengths than what we measured.This may be explained by an overestimation of our wavelengths, since we only measured to the zonal component, which defines the upper limit of the wavelengths, while the real values are shorter in the case of a north-south component.Vincent and Reid (1983) also presented shorter wavelengths using one direction only, however, this difference is partly due to the limits of our system that is not able to detect wavelengths shorter than 60 km.
In Fig. 5, histograms of IGW intrinsic periods and the time of day of their occurrence are shown.The histogram of the IGW daily cycle (left panel) has two weak maxima at early morning and early evening.The histogram of IGW periods (right panel) displays an exponential dependence.Highfrequency IGW with periods of 0.8-2 h are prevailing.Compared to Manson and Meek (1988), the spectrum is shifted towards longer wavelengths, which may be due to our system limits.The probability of occurrence is seen that IGW decreases with increasing period, as is also shown in the literature.
Histograms of the observed vertical and intrinsic horizontal phase velocity are shown in Fig. 6.The vertical phase speeds are generally small and the probability is approximately decreasing with phase velocity.The horizontal phase speeds are most probably found near 50 m/s, which is well in accordance with the results from literature (Vincent and Reid, 1983;Manson and Meek, 1988).There is no pronounced difference between the frequency of occurrence of eastward and westward propagating waves (54% vs. 46%).
As with the horizontal wavelengths, division of the vertical measuring volume does not show a clear effect.Therefore, in Fig. 7 intrinsic phase speeds are shown vs. the background wind.The upper and lower lines denote the 280 ms −1 limit of the phase speed, as described in Sect.2, while the two slant lines near the middle of the panel show the region of discarded waves, with periods of more than 5 h.The dashed line shows zero intrinsic phase speed.As expected, the intrinsic phase speeds decrease with background wind, so that preferably negative (westward) waves are found when the prevailing wind is positive (eastward), and vice versa.Considering the mean wind changes, however, we cannot distinguish between prevailing wind, tidal, and planetary wave changes.
In Fig. 8 dependences of horizontal intrinsic phase velocities and wavelengths on the IGW intrinsic period are shown.IGW phase velocities decrease with period.This, at the lower edge of the data, is due to the system limits |c-u| min = λ h,min /T*, with T* as the intrinsic period, while the decrease in the maximum phase speeds with period is qualitatively similar to that reported by Manson and Meek (1988).Our absolute maximum values are larger, which may be owing to the influence of waves with a north-south component.The horizontal wavelengths tend to increase with period (right panel of Fig. 8), which is in agreement with literature results.
The vertical wave characteristics are shown in Fig. 9.The vertical phase speed is also decreasing with period.Again, the lower limit is given by the system, however, it is not uniquely determined by T*, but is also depending on the background wind.The vertical wavelength vs. T* is presented in the right panel of Fig. 9.No clear connection between these parameters is visible.The vertical wavelengths tend to decrease with strong background winds (Fig. 10) due to the dispersion relation.
The obtained distributions of time-spatial IGW parameters agree with the results of other authors (Vincent and Reid, 1983;Manson and Meek, 1988;Gavrilov et al., 1996Gavrilov et al., , 1997;;Kalov and Gavrilov, 1985;Gavrilov and Medvedev, 1997), who have presented results about the horizontal IGW structure.Their distribution of horizontal IGW phase velocity, distribution of horizontal wave number, and also their dependence on IGW period, have similar character and are quantitatively of the same order of magnitude as has been shown here.However, Kalchenko et al. (1985) and Gavrilov et al. (1996Gavrilov et al. ( , 1997) ) showed a more uniform IGW frequency distribution, while the received distribution of IGW periods shown in Fig. 5, displays a clear domination of highfrequency IGW, i.e. of short-period waves.

Conclusions
In this paper a procedure of revealing the parameters of internal gravity waves from meteor radar wind data using wavelet analysis is described.The use of this algorithm has allowed us to identify waves and determine time-spatial IGW parameters.The distribution of horizontal and vertical wavelengths and directions of IGW energy propagation using meteor radar data has been obtained.However, our limitation of the meteor radar interferometric system to the zonal component means that the horizontal wave structure is not uniquely determined, but phase speeds and wavelengths must be considered as upper limits of the true values.Nevertheless, considering the fact that meteor radars with interferometric meteor position detection (also in both horizontal directions) are frequently and increasingly used in upper mesosphere/lower thermosphere studies provides the potential of analysing IGWs and their impact on the background circulation on a global scale.The majority of waves have been found to propagate upwards, although the portion of downward propagating waves is nonnegligible.High-frequency waves are dominating.The distribution of waves in the course of an average day is only weakly structured, with weak maxima in the morning and evening.
The data set used is limited to the 2 campaigns in the summers 1998 and 1999.Therefore, in this study, no information on the seasonal cycle of gravity wave activity and the potential variation of IGW parameters with season could be obtained.Further studies will be performed using a larger set of AG MRS data.

Figure 1 :
Figure 1: Schematic view of the division of the measuring volume into subvolumes of height hi and length di.Note that for the calculation we constructed overlapping subvolumes not shown in the figure.

Fig. 1 .
Fig. 1.Schematic view of the division of the measuring volume into sub-volumes of height h i and length d i .Note that for the calculation we constructed overlapping sub-volumes not shown in the figure.

Figure 2 :
Figure 2: Example of 5-min time series and fitted wavelet function.

Figure 3 :
Figure 3: Height-time cross-sections of background zonal winds during July and August1998, 1999.Data are calculated using a 96 h window each.

Fig. 3 .
Fig. 3. Height-time cross sections of background zonal winds during July and August 1998, 1999.Data are calculated using a 96-h window each.

Figure 5 :Fig. 5 .Figure 6 :Fig. 6 .
Figure 5: Histograms of IGW occurrence during the day (left panel) and period distribution (right panel) of IGWs.The upper and lower limits in the left panel are due to the resolution limitations of measurements and data analysis procedure.