Studies of gravity wave propagation in the mesosphere observed by MU radar

Mesospheric data were analyzed by a composite method combining phase and group velocity tracing technique and the spectra method of Stokes parameter analysis to obtain the propagation parameters of atmospheric gravity waves (AGW) in the height ranges between 63.6 and 99.3 km, observed using the MU radar at Shigaraki in Japan in the months of November and July in the years 1986, 1988 and 1989. The data of waves with downward phase velocity and the data of waves with upward phase velocity were independently treated. First, the vertical phase velocity and vertical group velocity as well as the characteristic wave period for each wave packet were obtained by phase and group velocity tracing technique. Then its horizontal wavelength, intrinsic wave period and horizontal group velocity were obtained by the dispersion relation. The intrinsic frequency and azimuth of wave vector of each wave packet were checked by Stokes parameters analysis. The results showed that the waves with intrinsic periods in the range 30 min–4.5 h had horizontal wavelength ranging from 25 to 240 km, vertical wavelength from 2.5 to 12 km, and horizontal group velocities from 15 to 60 m s −1. Both upward moving wave packets and downward moving wave packets had horizontal group velocities mostly directed in the sector between directions NNE (north-north-east) and SEE in the month of November, and mostly in the sector between directions NW and SWS in the month of July. Comparing with mean wind directions, the gravity waves appeared to be more likely to propagate along with mean wind than against it. This apparent prevalence for downstream wave packets was found to be caused by a systematic filtering effect existing in the process f phase and group velocity tracing analysis: A significant portion of upstream wave packets might have been Doppler shifted out of the vertical range in phase and group velocity tracing analysis.


Introduction
It was proved in a recent simulation study of various analysis methods (Lue and Kuo, 2012) that the errors in intrinsic period and horizontal propagation direction obtained by traditional methods (e.g., hodograph analysis, Stokes parameters analysis) were unacceptably large when data consisted of both upward propagating waves and downward propagating waves.During the past 30 years a number of groups have tried to measure gravity waves' propagation parameters in the mesosphere with ground-based radar at a single location.The horizontal phase velocities in the range of 20-90 m s −1 with wave periods in the range of 10 min to 10 h and horizontal wavelength of 40-1000 km were reported by several researchers (Vincent and Reid, 1983;Meek et al., 1985;Manson and Meek, 1988;Nakamura et al., 1993).The horizontal propagation direction has also been analyzed by traditional methods (Ebel et al., 1987;Manson and Meek, 1988;Tsuda et al., 1990;Nakamura et al., 1993;Gavrilov et al., 1997).But because data of upward propagating waves and downward propagating waves were not separately treated, no consistent pattern of horizontal propagation direction could be summarized from their results.In this paper we shall present the results of analyzing mesospheric data using a composite method combining the phase and group velocity tracing Published by Copernicus Publications on behalf of the European Geosciences Union.
technique (Kuo et al., 1998(Kuo et al., , 2007(Kuo et al., , 2008(Kuo et al., , 2009) ) and spectral method of Stokes parameters analysis (Vincent and Fritts, 1987;Eckermann and Vincent, 1989).Remarkably, quite often upward waves and downward waves coexisted at the same heights in the data of this study.So we must emphasize that the separation of upward waves and downward waves is essential in studying atmospheric gravity wave (AGW) propagation.However, the window to separate upward waves and downward waves was found to cause an unbalanced filtering on upstream wave and downstream wave.Its consequence will be discussed in this study.
2 Data and analysis procedure

Data
Wind velocity data were taken in 2-4 sampled days of each month of 1986, 1988 and 1989 by the MU radar (35 • N, 136 • E) at Shigaraki, Japan.Since the radar signal returned from mesosphere is only available during daytime, we chose data from 09:00-15:00 LT for analysis, during which the signal was the strongest and the amount of missing data was relatively small.In every inter-pulse period the radar antenna beam was steered sequentially toward the vertical and oblique directions at a zenith angle of 10 • (vertical → north → east → south → west).The beam width was 3 degrees, and the aspect sensitivity was negligible at the 10 degree zenith angle.The other observation parameters were observation range z = 63.6-99.3km; vertical resolution z = 300 m; and temporal resolution t was 147 s for 1986 and 1988 data, 222 s for July 1989 data, and 210 s for November 1989 data.There were many missing data due to insufficient signal power or time breaks during the experiment operation.Before doing data analysis, each missing data was filled by interpolation from its nearest neighbor good data, and the interpolation process was repeated until all the missing data within the data set were filled.A round of interpolation process was as follows: for missing data (denoted by 999) at time step I and height step K, if the nearest neighbor data at (I − 1, K) and (I + 1, K) were good, then the missing data at (I , K) would be replaced by the average of the data at (I − 1, K) and (I + 1, K).Otherwise, we would try to interpolate by the other nearest neighbor data at (I , K − 1) and (I , K + 1).If that interpolation also failed, we would try to interpolate by data at (I + 2, K) and (I − 2, K), then by (I , K + 2) and (I , K − 2), . . ., if the last try of interpolation by data at (I , K − 4) and (I , K + 4) also failed, we would leave the data at (I , K) as temporarily missing and would go on to repair other missing data until all data were exhausted.If missing data still existed after one round of interpolation process, we would go on to repeat another round of interpolation.We required all missing data in each data set under study to be successfully repaired.All data sets in this study had a small number of missing data that remained missing af- ter one round of interpolation, and were successfully repaired at the second round of interpolation.
The dual beam method (described in detail in Sect.2.1 of Kuo et al., 2008) is based on the assumption that east beam and west beam (respectively north and south beam) detect one wave at the same altitude and the same phase.The eastward velocity u and the vertical velocity w then can be determined from the Doppler velocities measured by the east beam V E and the west beam V W by Eq. ( 1): The velocities v and w can be similarly converted from the measurements of the north beam and south beam.The error made due to the assumption of constant phase is given by Eq. (2) (taken from Eq. (A5) in Appendix A of Kuo et al., 2008): where η = 2π z • tan 10 • λ x , z is the height and λ x is the projection wavelength of the wave along east-west (or northsouth) direction.Equation ( 2) is very useful to estimate the measurement error of a wave with known projection wavelength λ x along the dual beam line, and at height z.For example, at the altitude of 66 km (86 km), the minimum projection wavelength is 101.3 km (131.9 km) to guarantee the dual beam error to be smaller than 25 %.The dual beam errors for a gravity wave measured by east-west dual beam and north-south dual beam were generally different because the projection wavelengths were different.It was found in this study that the error at mesospheric height was so big that the vertical velocities obtained by east-west beams and north-south beams were often inconsistent with each other to such an extent that they had opposite sign at many heights.So instead of applying vertical flux of horizontal momentum analysis (which needs vertical velocities), we applied Stokes parameters analysis (which requires no vertical velocity) to help identify the azimuth of the most probable gravity wave packet propagation.In this study, 364 wave packets were measured by east-west dual beam and 319 wave packets were measured by north-south dual beam, totaling 683 wave packets analyzed.

Propagation parameters of wave packets calculated by dispersion relation
Each data set was separated by double-Fourier transformation over height and time (see Sect.A1 of Appendix A in Kuo et al., 2009), using windows defined in Table 1, into a data set of waves with downward phase velocity and a data set of waves with upward phase velocity for wave packet analysis.The windows in Table 1 set the upper-and lower-limits for vertical wavelength and observed wave period to be analyzed; its consequence on the outcomes of analysis will be discussed later.
Figure 1 shows wave structures after decomposition into upward and downward propagating waves.Several distinct wave structures are observed.The vertical phase velocity v pz and ground-based period τ can directly be estimated from the plot and are determined by the method of phase velocity tracing described in Appendix A2 of Kuo et al. (2009).The resulting values for v pz are indicated by phase lines (along the center package of waves).The waves also show distinct amplitude maxima.These are interpreted as upward (Fig. 1a) or downward (Fig. 1b) propagating wave packets.Applying group velocity tracing (Kuo et al., 2009), the corresponding group velocities v gz are determined (marked by energy lines across the packages of waves).It should be noted that this is based on two assumptions: first that the chief cause of the amplitude structure is the shape of the wave packet (i.e., that 765 Fig. 1b A partial height-time plot of ( ) 2 V δ of upward-phase propagating gravity waves converted from zonal velocity data set observed by MU radar at Shigaraki Japan on July 8 1986 (same data set of Fig. 1a). the wave packet is not distorted by, e.g., interference of different waves), and second that the vertical structure of the wave packet is horizontally slowly varying and thus that no additional altitude shifts are generated while the wave propagates obliquely through the radar beam.The wave packets are narrow (sometimes less than one vertical wavelength), which makes these assumptions the chief uncertainty cause of the determined group velocities.The observed wave frequency σ and vertical wave number m were consequently obtained by σ = 2π τ and m = σ v pz , respectively.
For every wave packet the best estimate of the wave vector is determined based on the dispersion relation and the Doppler shift between ground-based and intrinsic frequency and ground-based and intrinsic phase velocity.For an assumed horizontal phase speed and propagation azimuth ( ṽph , φph ), the horizontal wave vector k, ˜ and, using the full dispersion relation, the corresponding values ω, m, ṽpz and ṽgz , are calculated.The calculation is based on background atmospheric parameters of N = 2.09 × 10 −2 s −1 (5 min BVperiod; Gossard and Hooke, 1975), f = 8.31 × 10 −5 s −1 (21 h inertial period) at the latitude of the MU Radar site and H = 6 km (scale height).The background wind velocity ( ū, v) is obtained from the radar measurements by averaging over all heights and over all consecutive sampled days' (2-4 days') data for each respective month.Comparing the calculated values ṽpz and ṽgz with the actually measured values for v pz and v gz , the following cost functions are defined: (3) These cost functions are evaluated on a grid of phase speeds 2, 3, . . . , 198, 199, 200 (m s −1 ) and azimuths φ of 0.5, 1.5, 2.5, . . . , 357.5, 358.5, 359.5 (degree).P and G form tilted surfaces with basically opposite slope and thus S forms a clear minimum defining the best estimate values for ( ṽ * ph , φ * ) and (k * , * ), respectively.It should be noted in this context that, while the measurement error of v gz is larger than the error of v pz , P and G are equally weighted.A different weight would shift the location of the minimal value for S.
For ideal values of v pz and v gz , the cost function S should become zero for the best estimate (k * , * ).In reality the measured v pz and v gz have measurement errors and the minimum value for S is non-zero.In this study all cases were tested with notable amplitudes, where in principle the phase and group velocity tracing techniques could be applied, i.e., in order not to miss any prospective gravity wave packet, any three consecutive packages of waves which could be identified as wave packets were analyzed.This also included cases where the shape of the wave packet was less clearly pronounced than in the examples shown in Fig. 1, and also short wave events for which the group velocities could be identified only with much larger error ranges.This problem was enhanced by the application of the double Fourier-transform applied for separating upward and downward wave packets.The cost functions P and G were therefore also used to discern physically plausible cases from un-plausible ones by retaining only events where both P and G remained smaller than 0.15.From a total 683 cases, 471 were rejected and 212 were retained on this basis, including 141 observed by eastwest dual beams and 71 by north-south dual beams.Among these 212 potential wave packets, 113 were upward wave packets and 99 were downward wave packets.There are two solution-wave vectors (if existing) symmetric with respect to the mean wind velocity vector u h for each wave packet (see Sect. 2.2 of Kuo et al., 2009).Solving the mid-frequency approximation of the dispersion relation (see Appendix A) in order to determine (k * , * ) from v pz and v gz also shows that in principle two solutions, with horizontal wave vectors symmetric with respect to the mean wind velocity, exist if the measurement errors are sufficiently small.Eventually we must decide which one of the two solution-wave vectors is more likely to be the true solution.A possible method to help identify this true solution is Stokes parameters analysis.

Spectra method for Stokes parameters analysis
Stokes parameters analysis proposed by Vincent and Fritts (1987) provided intrinsic period and azimuth of a gravity wave.But the calculation of one of the parameters was not straightforward because it involved a 90 • rotation from the zonal perturbation velocity.Therefore, Eckermann and Vincent (1989) developed a spectral method for Stokes parameters analysis which can be summarized by the following equations: The meaning of each term and the fact that this method alone cannot separate upward and downward waves were explained in Lue and Kuo (2012).The range of summation over the vertical wave number m is to be properly selected to estimate the characteristic intrinsic period τ , azimuth ϕ as well as the degree of polarization dof the wave packet by the following equations: It is clear that Eq. (6a) cannot distinguish between the major axis orientation ϕ and ϕ ± 180 • , implying that Stokes parameters analysis cannot distinguish between eastward (northward) wave and westward (southward) wave.
It was concluded in the study of Lue and Kuo (2012) that the resulting azimuth from both methods of phase and group velocity tracing and Stokes parameters analysis were consistent with the statistical average of the azimuths of the component waves of the wave packet, and their resulting periods from Stokes parameters analysis tended to correspond to   1.
We propose that the wave packet-related averaging Stokes parameters P , Q, D, Ī in Eqs.(5a, b, c, d), and consequently ϕ and τ in Eqs. ( 6a) and (6d), respectively, can be obtained by averaging over a time interval t 0 − T < t < t 0 + T and summing over vertical wave numbers in the range m 0 − M < m • Z 2π < m 0 + M, with the time center t 0 coinciding with the time center of the corresponding wave packet, and the center vertical wave mode number m 0 given by m 0 = Z λ z , where λ z is the characteristic vertical wavelength of the wave packet and Z is the vertical range (35.7 km) under analysis.M = 1, 2, or 3, and T equals to one or half characteristic period of the wave packet.If the scale (T , M) was changed, these averaging values will also change correspondingly.Since these averages are statistical quantities, they can be used only as a reference to decide which one between the two symmetric solutions of the dispersion relation is more likely to be the true solution.In this study, we calculated the averages over six scales (T , M) = ( τ, 1), ( τ, 2), ( τ, 3), τ 2, 1 , τ 2, 2 and τ 2, 3 ; here τ is the characteristic wave period of the corresponding wave packet.For each calculation, azimuth ϕ and intrinsic period τ obtained from Stokes parameters analysis were compared with the corresponding azimuth φ and intrinsic period τ obtained from phase and group velocity tracing analysis.Then we picked the solution wave vector with its φ closest to ϕ as the wave vector of the virtual gravity wave packet if it also satisfied following conditions: and dual beam error < 25 %.
If such conditions ( 7) and ( 8) are not satisfied for each calculation, then the corresponding wave packet would be removed.A total 81 of the 212 possible gravity wave packets (38.2 %) were identified to be virtual gravity wave packets satisfying polarization relation with condition (7) and dual beam error condition (8).Among these 81 virtual wave packets, 39 were propagating upward and 41 were downward.

Examples of wave packet analysis by composite method
Let us have a look at some examples of wave packet analyses by the technique of phase and group velocity tracing, as shown in Fig. 1a, which is a partial height-time plot of (δV ) 2 (representing local power of the fluctuation velocity, details see Sect.A2 of Appendix A in Kuo et al., 2009) of upward propagating waves converted from zonal fluctuation velocity obtained on 8 July 1986.The corresponding downward propagating wave packets are shown in Fig. 1b.The locations (time and heights) of several upward wave packets (from Fig. 1a) and downward wave packets (from Fig. 1b) were indicated by their phase lines (along the packages of waves) and energy lines (across the packages of waves) in Fig. 2, which revealed that some upward wave packets (marked by red) and downward wave packets (marked by blue) coexisted at the same time and same height.The mean wind velocity had a magnitude of 14.46 m s −1 with azimuth angle of −93.02The observed characteristic wave period, vertical phase velocity, and vertical group velocity of the 3 wave packets from left to right (denoted by WP1, WP2, WP3, respectively) in Fig. 1a measured by phase and group velocity tracing technique were summarized in the 2nd-4th rows and the 2nd-4th columns in Table 3.The vertical wavelengths of 12.03, 9.83 and 12.72 km were obtained by definition, λ z = τ ob ×v pz , for WP1, WP2 and WP3, respectively.Then the observed vertical wave number, m = 2π λ z , was obtained for each wave packet.To search for the most probable horizontal wave vector k, ˜ to fit the measured values of v pz and v gz , any combination pair of ṽph , φ from the following values was used to calculate the characteristic gravity wave parameters for each measured wave packet: 2, 3, . . ., 198, 199, 200  For each pair k, ˜ = sin φ • 2π τ ṽph , cos φ • 2π τ ṽph , the corresponding intrinsic frequency ω, vertical wave number m, and vertical phase velocity ṽpz were obtained from Doppler equation ( ω = 2π τ − k • ū − ˜ • v), dispersion relation and definition for phase velocity, ṽpz = 2π m τ , respectively.The corresponding vertical group velocity ṽgz was obtained by a group velocity formula (see Eq. (3a) in Kuo et al., 2009).With the restriction of Eq. (3), we obtained the characteristic intrinsic period, horizontal wavelength and azimuth angle of the three wave packets as summarized in 5th-7th columns and 2nd-4th rows of Table 3.
The uncertainty of azimuth can be solved with the help of Stokes parameters analysis.First of all, we noticed that phase and group velocity technique tends to yield wave parameters corresponding to the high frequency part of the wave packet's component waves, while the Stokes parameters analysis tends to yield a result corresponding to the low frequency part of the wave packet's component waves (Lue and Kuo, 2012).So we used a time window of 30 min to 1 h (6th to 12th frequency mode) and wavelength window of 5.95 to 35.7 km (1st to 6th wave number mode, same as the window for velocity tracing analysis) to separate upward phase velocity waves from downward phase velocity waves.Then we applied the spectral method of Stokes parameters analysis to calculate the Stokes parameters from Eqs. (5a)-( 5d) and ( 6a)-(6d), and obtained intrinsic wave periods and azimuth angles for the three wave packets as summarized in 5th-7th rows, 5th and 7th columns of Table 3. Comparing with the results of phase and group velocity tracing technique, we decided that the most probable azimuths of the three wave packets were WP1: ϕ az = −85 • ; WP2: ϕ az = −100 • ; and WP3: ϕ az = −106 • .So their projection horizontal wavelengths along east-west line were 121, 98 and 131 km for WP1, WP2 and WP3, respectively.Substituting the east-west projection wavelength and the height of each wave packet into Eq.(2), we obtained the dual beam error of WP1, WP2 and WP3 to be 22.43, 32.06 and 21.03 %, respectively.All these three wave packets satisfy dispersion relation with Eq. (3) and polarization equation with Eq. ( 7), but WP2 did not meet the upper limit of 25 % for dual beam error (Eq.8) set to select wave packets for statistical analysis, so WP2 was rejected for further analysis.Among the three downward wave packets in Fig. 1b, only the leftmost packet was identified as a virtual gravity wave packet satisfying all the conditions, and the other two wave packets were rejected by dispersion relation with Eq. (3).

Results on AGW propagation parameters
Figure 3 shows the three-dimensional plot of horizontal wavelength λ h and vertical wavelength λ z versus intrinsic period τ in of the virtual gravity wave packets with upward group velocity (marked by circles) and with downward group velocity (marked by crosses) in the months of November.The corresponding plot for the month of July is presented in Fig. 4. The year of observation is distinguished by different colors (blue for 1986, red for 1988 and green for 1989).The ranges of intrinsic periods of both upward moving wave packets and downward moving wave packets in the month of November (Fig. 3) were 0.5-4.5 h, and their horizontal wavelengths and the vertical wavelengths ranges were 50-240 km and 2.5-12 km, respectively.τ λ λ with their projections on the ( , )   in h τ λ plane were drawn for easy viewing.Fig. 3. Three-dimensional plots of relation among horizontal wave length λ h , vertical wavelength λ z and intrinsic wave period τ in of upward and downward moving virtual gravity wave packets in the month of November."o" stands for upward-, while "+" for downward propagating virtual gravity wave packets; blue denotes 1986; red denotes 1988; and green denotes 1989.The vertical color lines connecting the points of τ in, λ h , λ z with their projections on the (τ in , λ h ) plane were drawn for easy viewing.
results in the month of July (Fig. 4) were quite different: their ranges of intrinsic periods and horizontal wavelengths were reduced to 0.5-2 h and 25-150 km, respectively, while the vertical wavelength range 2.5-14 km was similar to that in the month of November.The seasonal difference in τ in and λ h were caused by the difference in the mean wind's strength.November wind was significantly stronger than July wind, so the intrinsic wave period spreading in the month of November was much wider than that in the month of July.The relation between horizontal wavelength (vertical wavelength) and intrinsic wave period is that the wavelength grows with intrinsic wave period, consistent with the mesosphere study by Nakamura et al. (1993) and the study of lower stratosphere and troposphere by Kuo et al. (2009).
The horizontal wave vectors of the virtual gravity wave packets in the month of November (July) are plotted in Fig. 5a (Fig. 5b), the color lines indicating only the direction (no amplitude) of the corresponding mean wind vectors.The black circles in Fig. 5a and b from inside to outside represent horizontal wavelengths of 200, 100 and 50 km, respectively.The wave packets located outside the middle circle have horizontal wavelengths smaller than 100 km; however, the projection wavelengths of these wave packets, measured by east-west beams (or north-south beams), along the eastwest (or north-south) direction, are larger than 100 km.The results showed that the wave vectors were more likely to lie in the forward sector of the mean wind velocities than in the backward sectors.Such trend is even more convincing in the corresponding plots of horizontal group velocity plots (Fig. 6a and b) where color lines represent mean wind velocities (including direction and amplitude).Our lower stratosphere and troposphere study (Kuo et al., 2009) showed the same trend with much better statistics.It can be seen from Fig. 6a and b that the azimuths of horizontal group velocities were mostly distributed in the sector between NNE and SEE in the month of November (Fig. 6a), while in the sector between NW and SWS in the month of July (Fig. 6b).Clearly, wave packets apparently tend to propagate more likely along mean wind than opposite to it.

Doppler shift effect on wave propagation and momentum flux
Conventionally, horizontal propagation directions of gravity waves were referred by their wave vectors rather than by their group velocities.So we shall investigate the results on horizontal wave vectors in this study.In a windless situation, dispersion relation does not depend on azimuth angle of the horizontal wave vector, and the horizontal propagation direction is uniquely determined by polarization equation alone.When background wind is not negligible, dispersion relation and polarization relation combined determine the propagation direction.Table 4a.The dependence of the vertical phase velocity (2nd column), vertical group velocity (3rd column), vertical wavelength (4th column), the intrinsic period (5th column), vertical flux of zonal momentum to energy ratio (6th column), and vertical flux of meridional momentum to energy ratio (7th column) on the azimuth (1st column) of a gravity wave with period τ = 45 min and horizontal phase speed v ph = 35 m s −1 (horizontal wavelength λ h = 94.5 km) assumes that the wave is propagating in an eastward wind with wind speed u h = 30.9m s −1 ( azimuth angle ϕ az = 90 • ) in the mesosphere over the MU radar site.The results were obtained by gravity wave dispersion relation.which is the scattering plot of the horizontal momentum flux for both upward waves and downward waves for the month of November (July).These time-and height-averaged momentum fluxes were obtained for each sampled day.Then the downward flux of horizontal momentum was multiplied by −1 to yield the correct direction of the corresponding horizontal momentum.Figure 8b shows that the horizontal momentum of both upward waves and downward waves in the month of July were dominantly pointed against mean wind, completely contradicting Fig. 7b.The statistics of propagation direction revealed in Fig. 8a is also inconsistent with that of Fig. 7a, but to a lesser degree.
To uncover the effect of mean wind on the gravity wave's propagation parameters, we calculated dependences of the wave's intrinsic period, vertical phase and group velocities on the azimuth angle with the existence of mean wind.Some results are presented in 2nd-5th columns of Tables 4a and  b.The corresponding horizontal momentum flux to energy ratio was also calculated (by the formulas presented in Appendix B), and the results are listed in 6th and 7th column of Tables 4a and b.Existence of an eastward wind (90 • of azimuth) with wind speed of 30.9 m s −1 was assumed in these calculations.This mean wind speed corresponded to the mean wind of November 1988 in this study.In the first case (Table 4a), we consider a gravity wave with wave period of 45 min and horizontal phase speed of 35 m s −1 ; while in the second case (Table 4b), the same wave period of 45 min but higher horizontal phase speed of 60 m s −1 were considered.We calculated 7 azimuth angles from downstream (2nd row in Tables 4a and b) to upstream direction (last row in Tables 4a and b) with a step size of 30 • .The results showed that, as the azimuth angle's difference between mean wind and the gravity wave increased (from 2nd to bottom row), the vertical wavelength (4th column, ignore negative sign, which indicated downward propagation) increased while the intrinsic period (5th column) decreased.If the vertical wave length was larger than the observation range (35 km in this mesosphere studies, see Table 1), or the intrinsic period was smaller than B-V period (5 min), the wave packet under study would be systematically eliminated by dispersion relation.Even worse, the phase and group velocity tracing technique tends to yield the high frequency part of the composition waves of the wave packet, so the largest vertical wavelength obtained in this study was about 14 km (see Figs. 3 and  4).Therefore, many upstream wave packets such as those in the bottom 3 rows of Table 4a and bottom 4 rows of Table 4b Fig. 5a Plot of meridional wave number .vs. zonal wave number of virtual gravity wave packets in the month of November.'o' stands for upward, while '+' for downward propagating virtual gravity wave packets; Blue color denotes the year of 1986; red color denotes the year of 1988; and green color denoted the year of 1989.Colored straight lines represent the respective directions (not amplitude) of background mean wind.The inner, middle and outer black circles correspond to horizontal wavelengths of 200, 100, and 50 km respectively.
Upward and downward wave packets / November Fig. 5a.Plot of meridional wave number vs. zonal wave number of virtual gravity wave packets in the month of November."o" stands for upward, while "+" for downward propagating virtual gravity wave packets; blue denotes 1986; red denotes 1988; and green denotes 1989.Colored lines represent the respective directions (not amplitude) of background mean wind.The inner, middle and outer black circles correspond to horizontal wavelengths of 200, 100, and 50 km, respectively.
would be systematically filtered by this mean wind.There is no doubt that the effectiveness of upstream waves filtering depends on the strength of background wind.This filtering effect may explain why the horizontal propagation direction of gravity wave packets were dominated by downstream wave packets in this study of mesosphere and our previous study of lower stratosphere and troposphere (Kuo et al., 2009).The vertical flux (of zonal momentum) to energy ratio shown in the 6th column of Tables 4a and b suggest that the Doppler effect seems to upgrade (downgrade) the momentum flux of upstream waves (downstream waves).This effect might have contributed to the dominance of upstream momentum in Fig. 8b, and might also explain why most conventional momentum flux studies suggest that upstream waves dominate.

Summary and discussion
Horizontal wind velocities measured by the Shigaraki MU radar are analyzed by a space-time Fourier transform and separated into upward and downward propagating waves.
After back-transformation into the spatiotemporal domain, wave packets are identified and the vertical phase velocity v pz and the vertical group velocity v gz are determined by phase and group velocity tracing (Kuo et al., 2009).The latter relies on the assumption that the propagation of the wave packet is the dominant factor in determining the apparent shape of the wave packages in the vertical-time cross sections.From the full gravity wave dispersion relation and Doppler shift relations, a complete set of wave parameters is estimated by optimization.The depth of the minimum in the cost function is used to distinguish whether candidate wave packets are physically plausible.In this step, out of a total 683 candidate wave packets, 212 events were retained as potential wave packets.To make sure these parameters also satisfied polarization relation, the vicinity of each possible gravity wave packet was further analyzed by spectral method of Stokes parameters analysis to obtain intrinsic wave period and azimuths.Among these 212 possible gravity wave packets, only 81 were accepted as virtual gravity wave packets because their periods and azimuths obtained by the spectral method of Stokes parameters analysis were consistent with the corresponding values obtained by phase and group velocity tracing technique, satisfying Eqs. ( 7) and (8).Our composite method not only takes care of the dispersion relation and the polarization relation but also takes care of dual beam error.The result emphasizes the dominant role of mean wind in the dispersion relation.It was found in this study that gravity waves apparently propagated more likely along mean wind than opposite to it.Very probably many upstream wave packets would have been Doppler shifted out of the observation window of phase and group velocity tracing analysis if they ever existed.
It was concluded in a simulation study of various analysis methods (Lue and Kuo, 2012) that upward waves and downward waves must be separated and analyzed independently by any method of analysis.It was quite often in this study that upward waves and downward waves coexisted at the same height, as demonstrated in Fig. 2, which justifies the necessity of up-and down waves separation.This could also explain one of the reasons why our results were inconsistent with others' results.The simulation study also revealed that phase and group velocity tracing technique tends to display the characteristics of the high frequency part of component waves of the wave packet, while hodograph analysis (hence Stokes parameters analysis also) tends to display the characteristics of low frequency component waves.Furthermore, the unpublished result in that simulation study showed that momentum flux analysis also tends to display the characteristics of low frequency component waves.Therefore, phase and group velocity tracing technique will experience more serious upstream waves filtering than Stokes parameters analysis and momentum flux analysis.This may also account for the difference in intrinsic periods between this analysis and the analyses by Tsuda et al. (1990) and Nakamura et al. (1993).velocity along the wind direction equals wind velocity, dynamic instability will arise, and its wave amplitude will grow and break.Therefore, some downstream waves will be filtered by background wind while upstream waves will be free to propagate to higher atmosphere.This theory is supported by the results of vertical flux of horizontal momentum by conventional method, which did not separate upward waves from downward waves before momentum flux analysis.But in a study of vertical flux of horizontal momentum in the lower stratosphere and troposphere made by Kuo et al. (2008), in which upward waves and downward waves were separately treated, the zonal momentum of both upward waves and downward waves in the quiet day with relatively small momentum flux were found to be dominantly eastward (i.e., propagating along the mean wind).However, in the active days with relatively large momentum flux, the upward flux was dominated by westward propagating waves (i.e., propagating against mean wind) while the downward flux was dominated by eastward propagating waves (i.e., propagating along the mean wind).In other words, in the active days both upstream waves and downstream waves were detected to be dominating in different environments.Furthermore, in the active days, downward flux was about 50 % larger than upward flux, implying that downstream waves prevail over upstream waves.This result (Kuo et al., 2008) does not agree with the conventional analysis because of the process of separating upward waves from downward waves before momentum flux analysis.The upward-downwardwaves separation must be done by double Fourier transformation over a height and time window, leading to unbalanced filtering effect.The result of this study does not rule out the possibility of a prevailing upstream gravity waves propagating against mean wind, as conventional theory insists.Now we have a dilemma in gravity wave propagation analysis: On one hand, upward waves and downward waves must be separately analyzed to guarantee accuracy; on the other hand, the height-time window in double Fourier transformation to do up-and down-waves separation yields an unbalanced filtering effect between upstream waves and downstream waves, leading to a misleading result.This dilemma needs to be solved.
Another interesting result in this study is that there were many wave packets in the mesosphere with downward group velocities, whose number was almost comparable to that of wave packets with upward group velocities (see Fig. 5a and  b).The comparison of upward waves and downward waves in terms of wave packet number and mean wave energy are summarized in Table 5, which shows that the upward waves clearly dominate over downward waves in terms of wave energy, but are comparable in terms of wave packet's number.If the source of gravity wave is in the lower atmosphere, how were these downward gravity wave packets generated?We speculate that elastic scattering (ES) and parametric subharmonic instability (PSI) might be responsible.Because atmospheric density decreases with height, the amplitude of an upward moving gravity wave will grow with height.When  the amplitude grows to a level that PSI occurs, the wave may break into two oppositely propagating waves, one upward and one downward.Another possibility is that a small wave with faster upward group velocity catches up to a much larger wave but with slower upward group velocity; the smaller waves may suffer an ES to become a downward wave.An indirect proof of the existence of PSI and ES processes seems to have been revealed in our previous momentum flux analysis of the MU data in the lower stratosphere and troposphere (Kuo et al., 2008): In a quiet day, both upward flux and downward flux were dominated by eastward momentum, implying ES might be responsible; in the active days, upward flux was dominated by westward momentum while downward flux was dominated by eastward momentum, implying PSI might be responsible.Direct proof would require a detailed analysis of upward-and downward wave packet pairs.This is our current interest of research.

Fig. 1a A
Fig.1a A partial height-time plot of ( ) 2 V δ Fig. 1b.A partial height-time plot of (δV ) 2 of upward-phase propagating gravity waves converted from zonal velocity data set observed by MU radar at Shigaraki, Japan, on 8 July 1986 (same data set as Fig.1a).

Fig. 4 Fig. 4 .
Fig.4 Same as Fig.3 except for the month of July.
Figure 7a (7b) displays the scattering plots of horizontal wave vectors of the possible gravity wave packets in the month of November (July) satisfying Doppler relation and dispersion relation with Eq. (3).There were always two solutions symmetric with respect to the mean wind direction, as demonstrated in Fig.7a and b.Clearly, the distributions were dominated by downstream wave packets in both figures, inconsistent with most conventional results of momentum flux studies.For the purpose of direct comparison with

Fig. 7a (
Fig. 7a (7b), the corresponding time and height averages of the vertical flux of zonal momentum ( u w ) and meridional momentum ( v w ) were analyzed as shown in Fig.8a(8b), which is the scattering plot of the horizontal momentum flux for both upward waves and downward waves for the month of November (July).These time-and height-averaged momentum fluxes were obtained for each sampled day.Then the downward flux of horizontal momentum was multiplied by −1 to yield the correct direction of the corresponding horizontal momentum.Figure8bshows that the horizontal momentum of both upward waves and downward waves in the month of July were dominantly pointed against mean wind, completely contradicting Fig.7b.The statistics of propagation direction revealed in Fig.8ais also inconsistent with that of Fig.7a, but to a lesser degree.To uncover the effect of mean wind on the gravity wave's propagation parameters, we calculated dependences of the wave's intrinsic period, vertical phase and group velocities on the azimuth angle with the existence of mean wind.Some results are presented in 2nd-5th columns of Tables4a and b.The corresponding horizontal momentum flux to energy ratio was also calculated (by the formulas presented in Appendix B), and the results are listed in 6th and 7th column of Tables4a and b.Existence of an eastward wind (90 • of

875Fig.
Fig.5b Same as Fig.5a except for the month of July.
Fig.6b Same as Fig.6a except for the month of July.

Fig. 7b
Fig.7b Same as Fig.7a except for the month of July.
Fig.8a Plot of zonal momentum flux .vs. meridional momentum flux of upward waves (marked by open circle) and downward waves (marked by cross) for the month of November data.All downward fluxes of horizontal momentum were multiplied by − 1 in this plot.

Fig. 8a .
Fig. 8a.Plot of zonal momentum flux vs. meridional momentum flux of upward waves (marked by open circle) and downward waves (marked by cross) for November data.All downward fluxes of horizontal momentum were multiplied by −1 in this plot.

Table 1 .
Wave periods bands and wave length bands corresponding to the frequency-and wave number-windows used to sample the wave packets for phase and group velocity tracing analysis.

Table 2 .
Same as Table1except that the window was used for Stokes parameters analysis.lowfrequencypart among the component waves, while the result from phase and group velocity tracing technique tended to correspond to the high frequency part among the component waves of the wave packet.Therefore, we used windows defined in Table2to generate data for the Stokes parameters analysis on the corresponding wave packets obtained from the windows defined in Table the

www.ann-geophys.net/31/845/2013/ Ann. Geophys., 31, 845-858, 2013 Fig.6a Plot
of meridional group velocity .vs. zonal group velocity of virtual gravity wave packets in the month of November.'o' stand for upward, while '+' for downward propagating virtual gravity wave packets; Blue color denoted the year of 1986; red color denoted the year of 1988; and green color denoted the year of 1989.Colored straights lines represent the respective mean wind velocities.Plot of meridional group velocity vs. zonal group velocity of virtual gravity wave packets in the month of November."o" stands for upward, while "+" for downward propagating virtual gravity wave packets; blue denotes 1986; red denotes 1988; and green denotes 1989.Colored lines represent the respective mean wind velocities.

Table 5 .
Mean energy u 2 + v 2 , the number of potential wave packets N, and their up/down ratio for the different years and months.•• • represents average over all heights, and over-bar represents time average over all sampling days of the same month.The up/down ratio represents the ratio of the quantities of upward waves to the downward waves.