Studies of vertical fluxes of horizontal momentum in the lower atmosphere using the MU-radar

We study the momentum flux of the atmospheric motions in the height ranges between 6 and 22 km observed using the MU radar at Shigaraki in Japan during a 3 day period in January 1988. The data were divided by double Fourier transformation into data set of waves with downward- phase- velocity and data set of waves with upward-phase-velocity for independent momentum flux calculation. The result showed that both the 72 h averaged upward flux and downward flux of zonal momentum were negative at nearly each height, meaning that the upward flux was dominated by westward propagating waves while the downward flux was dominated by eastward propagating waves. The magnitude of the downward flux was approximately a factor of 1.5 larger than the upward flux for waves in the 2~7 h and 7~24 h period bands, and about equal to the upward flux in the 10–30 min and 30 min–2 h period bands. It is also observed that the vertical flux of zonal momentum tended to be small in each frequency band at the altitudes below the jet maximum (10~12 km), and the flux increased toward more negative values to reach a negative maximum at some altitude well above the jet maximum. Daily averaged flux showed tremendous variation: The 1st 24 h (quiet day) was relatively quiet, and the fluxes of the 2nd and 3rd 24 h (active days) increased sharply. Moreover, the upward fluxes of zonal momentum below 17 km in the quiet day for each period band (10~30 min, 30 min~2 h, 2~7 h, and 7~24 h) were dominantly positive, while the corresponding downward fluxes were dominantly negative, meaning that the zonal momentum below 17 km in each period band under study were dominantly eastward (propagating along the mean wind). In the active days, both the upward fluxes and downward fluxes in each frequency band were dominantly negative throughout the whole altitude range 6.1–18.95 km.


Introduction
Gravity wave momentum fluxes observed by radars have been used to estimate the body force due to gravity waves in the lower stratosphere and in the mesosphere and lower thermosphere (Fritts and Alexander, 2003). Technically, the symmetric-beam radar method of Vincent and Reid (1983) has been widely applied to measure momentum flux in both the lower atmosphere (e.g. Fukao et al., 1988;McAfee et al., 1989;Sato, 1990Sato, , 1993Fritts et al., 1990;Thomas et al., 1992;Prichard and Thomas, 1993;Worthington and Thomas, 1996) and the mesosphere (e.g. Tsuda et al., 1990;Reid and Vincent, 1987). Mean momentum fluxes in the lower stratosphere over the MU radar site (at 35 • N latitude) were ∼0.1-0.3 m 2 s −2 , with peak value of ∼1 m 2 s −2 and implied mean flow forcing of ∼1 ms −1 d −1 (Fritts et al., 1990). Momentum fluxes in the mesosphere and lower thermosphere have been measured more widely. At extratropical sites, mean values of momentum fluxes of ∼1-10 m 2 s −2 have been measured with mean forcing in the range ∼10-70 ms −1 d −1 that generally oppose the mean winds (Reid and Vincent, 1987;Fritts and Yuan, 1989;Tsuda et al., 1990;Nakamura et al., 1993). At the tropical Jicamarca radar site (12 • S) (Hitchman et al., 1992), momentum fluxes of 3766 F. S. Kuo et al.: Studies of vertical fluxes of horizontal momentum proposed by Vincent (1984) was widely applied to distinguish what fractions of the gravity wave spectrum are propagating upward and downward. Usually, wave amplitudes and modulus of momentum fluxes are larger for atmospheric gravity waves having upward energy fluxes and propagate opposite to the mean wind in the upper middle atmosphere (Nakamura et al., 1993;Gavrilov et al., 1995Gavrilov et al., , 1996Gavrilov et al., , 1997Hall et al., 1995). We notice that all these fluxes were estimated by conventional statistical method of analysis, upward flux data and downward flux data were not separated before momentum flux calculation, and none of these measurements had directly measured the group velocity of the wave propagation which would determine the flux direction.
It is well known that the energy and momentum of any wave propagation are transported at its group velocity, and that the vertical group velocity of an atmospheric gravity wave may be in either the opposite or the same direction of its vertical phase velocity depending on the relation of its horizontal phase propagation speed v ph and the projection mean wind velocity u h (Kuo et al., 2003;Kuo and Röttger, 2005;Gavrilov et al., 1997). We call a wave packet a type 1 wave packet if its group velocity and phase velocity have opposite sense of vertical propagation, and a type 2 wave packet if they have same sense of vertical propagation. From the linear theory of the gravity wave dispersion relation, a wave packet will appear as a type 1 wave packet if v ph >u h . Otherwise (u h >v ph ), it will appear as a type 2 wave packet. Here u h is the projection wind velocity along the horizontal direction of phase propagation (so, v ph ≥0 always, and u h can be either positive or negative). Therefore direct measurement of the vertical group velocity is required to gain the information on the vertical energy transport by gravity waves in the atmosphere. A method called phase-and group-velocity tracing was developed (Kuo et al., 1998) to measure the vertical phase-and group-velocities of any kind of wave propagation. And this method was further proved by numerical simulation to be reasonably accurate in measuring the phaseand group-velocities of the projected wave packet motion along any oblique radar beam (Kuo et al., 2007).
The technique of velocity tracing was applied to analyze the data observed by SOUSY-Svalbard Radar in October 2000 in the altitude range between 2.4 km and 17.4 km (Kuo et al., 2003), in which hundreds of wave packets with characteristic wave periods in the range 15∼45 min were identified from Range-Time plot of (δV ) 2 : About 85% of them were type 1 and about 15% were type 2 wave packets, revealing the typical characteristics of the propagation of atmospheric gravity waves in the linear regime. Then, for another set of data taken during the period between 30 March and 24 April 2000 in the altitude range between 2.85 km and 14.85 km also using the SOUSY Svalbard Radar, with the simultaneous information on the background wind also available, we not only calculated the horizontal wave-number component of a wave packet along the direction of background wind velocity by the dispersion equation, but also deduced the angle between the horizontal wave vector and its background wind velocity by the polarization equation of gravity wave (Kuo and Röttger, 2005). Again, most of the wave packets (76%) with characteristic wave period of ∼70 min in the Range-Time plot of (δV ) 2 were type 1 wave packets, and the rest (24%) were type 2 wave packets. In this study we will combine this wave packet analysis technique with the conventional method of momentum flux measurement to estimate the body forcing due to gravity waves in the lower atmosphere. The upward flux and downward flux will be treated separately, and the vertical group velocities of as many as possible sampled wave packets will be measured.

Data and analysis procedure
The data of the wind velocity were taken from 08:18 on 5 January to 15:09 on 8 January 1988 by the MU radar (35 • N, 136 • E) at Shigaraki, Japan. 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=5∼24 km, range resolution z=150 m, time resolution t=150 s. There were some missing data points 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. Figure 1 shows the 72 h averaged (upper panel) and 24 h averaged wind velocity profiles (lower panel). A strong zonal wind jet can be seen in the height range 10-12 km with large vertical shear both above and below this jet. The meridional wind was evidently much weaker (than zonal wind) but its day-to-day variation (in terms of relative magnitudes) was larger than that of the zonal wind. The vertical wave number and frequency spectra obtained from this data set had been studied and presented by Kuo et al. (1992).
To discuss the analysis method, let u, v and w denote the wind velocities in the zonal-, meridional-and verticaldirection respectively as usual. The eastward oblique beam with zenith angle θ=10 • would detect a Doppler velocity V E , while the westward oblique beam would detect a Doppler velocity V W given by Conversely, the wind velocities u and w can be obtained from the measurement values of V E and V W by following equations, The velocities v and w can be similarly converted from the measurements of the north beam and south beam. Equations (2a) and (2b) represent respectively the average (or, linear interpolation) of the eastward velocities u and vertical velocities w obtained by the eastward and westward oblique beams at the same height (but different locations). These values of u and w can be regarded as the measurements along the vertical axis. Since the eastward-and westward-oblique beams pointed to different phases of the wave, the error of measurement by the dual beam method should be seriously considered. The error in dual beam method and its effect on the momentum flux analysis were discussed in Appendix A, and it was found that if the horizontal wavelength measured along the azimuthal direction of the relevant beam pair is a factor of 8 larger than the horizontal spacing between the radar volumes of the 2 opposite beams, then the errors are less than 8%. The dual beam method was found to be valid in this study.
The wind velocity u consists of a mean windū and a perturbation velocity component u . The perturbation velocity u is resulted from the combined contributions from both the wave packets with upward group velocity and the wave packets with downward group veloity: where u U (u D ) is contributed by the wave packets with upward (downward) group velocity. Similarly, By conventional definition of the vertical flux of zonal momentum, a positive (negative) value of the time average u w may represent either an upward flux of eastward (westward) momentum or a downward flux of westward (eastward) momentum. Such uncertainty of expression always exists unless either the direction of the zonal momentum or the vertical group velocity of the wave packet is known. Physically, the vertical flux of zonal momentum is a component of a tensor defined as EZ =P E v gz , where P E represents the zonal momentum density and v gz represents the vertical group velocity. Both P E and v gz are vector components: P E is positive (negative) for eastward (westward) momentum, and v gz is positive (negative) when the wave packet is moving upward (downward). Technically, we are able to separate the waves with upward phase velocity from the waves with downward phase velocity by double Fourier analysis (see Sect. B1 of Appendix B). Since our previous studies revealed that the majority of wave packets were type 1 wave packets (Kuo et al., 2003;Kuo and Röttger, 2005), we may approximately assume that the wave packets with downward phase velocity have upward group velocity, and vice versa. This assumption is necessary because the momentum is transported at its  "open circle" for the 1st day zonal wind, "line" for the 2nd day zonal wind, "star" for the 3rd day zonal wind; "dot" for the 1st day meridional wind, "cross" for the 2nd day meridional wind, and "square" for the 3rd day meridional wind. The starting time of observation was 12:28 on 5 January 1988.
group velocity (not phase velocity), but our method of wave separation can not separate the group velocity-upward waves from the group velocity-downward waves. Therefore we emphasize that the validity of such assumption should be statistically confirmed by wave packet analysis case by case. So the separation of upward flux from downward flux becomes possible by our analysis procedure, and the vertical flux of the zonal momentum can be written as, The first term in the right hand side of Eq. (6) is an upward flux of eastward-(westward-) momentum if it is positive (negative), and the second term is a downward flux of The third and fourth terms are physically meaningless, because u U and w D (u D and w U ) are definitely not the velocity components of the same wave packet, and they will contribute some random error to the momentum flux calculation. It was found in this study that these two error terms are negligibly small when the wave activity was strong. We will analyze the first and second term separately, and all these discussions apply equally well to the vertical flux of meridionalmomentum, which is given by, The continuous data of 72 h starting from 12:28 on 5 January 1988 in the altitude range 5.75 km∼23.3 km was selected for the momentum flux analysis. First of all the phase-upwardpropagating waves and the phase-downward-propagating waves were separated for the data sets of u , v and w , and the procedure of the wave separation and the velocity tracing technique were described briefly in the previous studies (Kuo et al., 2003;Kuo and Röttger, 2005), and it is put in Appendix B in this paper to help the readers. So the zonal wind data set {u (z, t)} was divided into upward-group-velocity (downward-phase-velocity) subset {u U (z, t)} and downward-group-velocity (upward-phasevelocity) subset {u D (z, t)}, using a vertical wave number window 3.55×10 −4 m −1 <k z <8.52×10 −3 m −1 (corresponding to 17.7 km>λ>0.7375 km wave length band), coupled with several frequency windows σ 1 >σ >σ 2 (corresponding to wave periods of τ 1 <τ <τ 2 ). So was {v (z, t)} divided into {v U (z, t)} and {v D (z, t)}; and {w (z, t)} into {w U (z, t)} and {w D (z, t)}. The frequency windows corresponding to wave period bands of 10∼30 min, 30 min∼2 h, 2∼7 h, 7∼24 h, and the whole period band 10 min∼24 h were used in this study.

Mean momentum flux profiles
The 72 h averaged profiles of upward-flux (downward-flux) of zonal (marked by square)-and meridional (marked by cross)-momentum for each frequency band, obtained by the  first term (second term) of Eqs. (6) and (7) Table 1. Evidently the meridional momentum fluxes were negligibly small comparing with the corresponding zonal momentum fluxes in the three higher frequency bands (10-30 min, 30 min-2 h and 2-7 h period bands) as can be seen from Figs. 2 and 3. The meridional momentum flux was comparable to (but still significantly smaller than) the zonal momentum flux only in the lowest frequency band (7-24 h period band). So we shall discuss only the   Table 1). It is also observed that the vertical flux of zonal momentum tended to be small in each frequency band at the altitudes below the jet maximum (10∼12 km), and the flux increased toward more negative values to reach a negative maximum at some altitude well above the jet maximum. It is noticed from Table 1 that the peak value of the downward flux in each of the two lower frequency bands (2∼7 h and 7∼24 h period bands) was larger than that of the upward flux approximately by a factor of 1.5∼2, while the downward flux and upward flux of zonal momentum in the two higher frequency bands (10∼30 min, 30 min∼2 h period bands) were approximately equal.

Daily momentum flux profiles
Consecutive 24 h averaged profiles of vertical fluxes of zonal and meridional momentum for 10-30 min period band, with 1.5 km running mean, were shown in Fig  "Open circle" for the 1st day zonal momentum; "line" for 2nd day zonal momentum; "star" for 3rd day zonal momentum; "dot" for 1st day meridional momentum; "cross" for 2nd day meridional momentum; "open square" for 3rd day meridional momentum.
tum fluxes in the 1st 24 h ("open circle" for zonal momentum and "dot" for meridional momentum) were negligibly small comparing with the successive daily profiles. For the convenience of discussion we shall refer the 1st 24 h as a quiet period, and the next 48 h as an active period. In the active period, both the upward flux and downward flux of zonal momentum ("line" for 2nd day and "star" for 3rd day) in each frequency band were dominantly negative throughout the whole altitude range 6.1-18.95 km. The peak of the upward flux of the 30 min-2 h period band (top panel of Fig. 5) evidently moved up from the altitude of 13.1 km ("line") to 15.5 km ("star") in 24 h with a speed of upward movement  Fig. 7) was dominantly positive, while the corresponding downward flux ("star") was negative. In other words, the zonal momentum below 17 km in each frequency band was dominantly eastward. This characteristic is tremendously different from that of the active period. To compare the flux magnitudes between the quiet and the active periods, we listed the 24 h mean fluxes of the 1st 24 h (quiet day) and the 3rd 24 h (most active day) in Table 2 and Table 3, respectively. It can be seen that the magnitudes of the fluxes of zonal momentum in the quiet day  ( Table 2) were smaller than the corresponding fluxes in the most active day (Table 3) by a factor 4∼7 in the two lower frequency bands (2∼7 h, 7∼24 h period bands), and by a factor 17∼70 in the two higher frequency bands (10∼30 min, 30 min∼2 h period bands). In other words, the magnitudes of the fluxes of low frequency waves (in 2∼7 h, 7∼24 h period bands) increased by a factor of 4∼7, while the magnitudes of the fluxes of high frequency waves (in 10∼30 min, 30∼2 h period bands) increased by a surprisingly large factor of 17∼70 on the way from quiet phase to active phase. Such extraordinarily large scale of amplification in the higher frequency waves might be resulted from the wave breaking of larger scale waves.

Vertical profiles of induced drag
Mean profiles of induced zonal and meridional accelerations were computed from the vertical momentum flux profiles for each frequency band using the relation (8a) and (8b): whereū andv are the mean zonal and meridional velocities respectively. The atmospheric density at the altitude z is taken as ρ 0 (z) =1.29 exp −z H kg m 3 , with scale height H =8 km. Since the numerical derivative computation can be irrationally amplified by random error, the momentum flux profiles were repeatedly smoothed with a 2.25 km running mean until all the small scale random fluctuations in the flux profiles were completely removed. The results of 72 h profiles were shown in Fig. 8a, and the 24 h average profiles were shown in Fig. 8b and c.
The 72 h average profiles of density weighted momentum flux ρ 0 u U w U +u D w D (marked by square) and ρ 0 v U w U +v D w D (marked by cross) were shown for each frequency band (1st ∼4th profiles from left) and for the sum  of all wave period bands (the right-most profile) in the upper panel of Fig. 8a. The corresponding inferred acceleration profiles were shown in the lower panel of Fig. 8a. The upper panel (of Fig. 8a) revealed a tendency for the zonal momentum flux to be negative and divergent at lower levels and negative and convergent at greater heights, contributing to eastward inferred acceleration at lower altitudes and westward accelerations at higher altitudes as shown in the lower panel. Such features were also revealed in March 1986 data obtained by MU radar (Fritts et al., 1990). But the inferred total zonal acceleration in this study revealed much more violent variation over height: it was estimated to be ∼+9 m/s/day at 12-13 km and ∼−20 m/s/day at 18 km with smooth variation  at intermediate heights. Between 8 and 13 km, the acceleration decreased with decreasing height from +9 m/s/day to +3 m/s. The magnitude of this zonal acceleration is about a factor of 10, and the peak amplitude of the mean (density weighted) zonal momentum flux (in whole 10 min ∼24 h period band) were about a factor of 5 larger than the corresponding value of March 1986 data (Fritts et al., 1990).
The daily profiles of density-weighted zonal momentum flux (upper panel) and meridional momentum flux (lower panel) in the 10 min-24 h period range were shown in Fig. 8b, and the corresponding inferred acceleration profiles were shown in Fig. 8c. The 1st day data was represented by "open circle", the 2nd day by "line" and the 3rd day by "star". It is noticed that both the zonal momentum flux (open circle, upper panel of Fig. 8b) and the inferred zonal acceleration in the quiet day (open circle, upper panel of Fig. 8c) were slightly smaller than the corresponding values of March 1986 data.

On the relation between momentum flux and wave propagation
The polarization relation between the amplitudes of zonal component U and the vertical component W is known to be (Gossard and Hooke, 1975), where θ 1 = tan −1 ωk 2 z , θ 2 = tan −1 m , and ∼ =3.2×10 −5 m −1 ; N ∼ =1.05×10 −2 s −1 (V-B frequency), 2 z ∼ =8.31×10 −5 s −1 (inertial frequency); ω, k, and m are the intrinsic frequency, zonal-, meridional-, and verticalwave number component respectively, and their relations with the observed wave frequency σ are governed by the gravity wave dispersion equation. The zonal and vertical components of fluctuation velocity u and w of a wave mode with frequency σ , wave numbers k, and m are given by, The vertical flux of zonal momentum is defined as the time average of u w , Similarly, the vertical flux of meridional momentum is obtained as, v w = ω ω 2 −4 2 z N 2 −ω 2 m 2 + 2 ω 2 2 +4 2 z k 2 V 2 × 1 2 cos (θ 3 +θ 2 ) (12) where, θ 3 = tan −1 ω −2 z k . Evidently the sign of u w and v w are decided by cos (θ 1 +θ 2 ) and cos (θ 3 +θ 1 ), respectively. We have made numerical calculation of Eqs. (11) and (12) along with the dispersion equation under following conditions: The mean wind velocity is assumed to be u=45 m/s and v=10 m/s; The characteristic wave periods of different wave packets under study were 15 min, 40 min, 2 h and 8 h; The characteristic horizontal phase speed under consideration were 30, 45 and 60 m/s, and the propagation directions expressed in terms of azimuth angle varied from 7 • to 357 • by a step size of 10 • . The result of calculation shows that the sign of u w and v w are always the same as the sign of  v gz × v gx and v gz ×v gy respectively for type 1 wave packets, where v gz ,v gx and v gy represent the vertical-, zonal-and meridional-group velocity of the wave packet under study. This result confirms our earlier statement that the direction of the vertical flux is defined by the direction of the vertical group velocity and a positive (negative) value of the time average u w may represent either an upward flux of eastward (westward) momentum if the v gz >0 or a downward flux of westward (eastward) momentum if v gz <0. Therefore the momentum flux profiles observed in the quiet day as described in Sect. 2.3 indicated that the wave packets tended to propagate in the same direction of the mean wind. More experiments are needed to check whether this is always favored by nature. It should also be noted that the above discussion is based on idealized assumptions using one wave in contrast to the nature with a superposition of different waves.

Wave packet analysis of vertical group velocity
Atmospheric gravity wave is a highly localized phenomenon with its energy and momentum transported by wave packet. The energy flux (a vector) transported by a wave packet is the product of its energy density (a scalar) multiplied by its group velocity (a vector). Similarly, a component of the momentum flux (a tensor) transported by a wave packet is the product of the corresponding component of its momentum density (a vector) multiplied by the proper component of its group velocity (a vector). For instance, the vertical flux of zonal momentum is the product of the zonal momentum density multiplied by the vertical group velocity. The mean group velocity of the gravity waves may be estimated by averaging the group velocities of all the sampled wave packets, and the group velocity of a wave packet can be measured by the technique called phase-and group-velocity tracing (Kuo et al., 1998, see also Sect. B2 of Appendix B). Previous analyses of SOUSY-Svalbard data (Kuo et al., 2003;Kuo and Röttger, 2005) repeatedly showed that the characteristic wave period (wave length) of a wave packet was solely determined by the high frequency (large wave number) end of the frequency (wave number) window. Changing the low frequency (small wave number) end of the window has very small effect on the measurement of the characteristic wave period (wave length) of the wave packet. Based on this understanding, we have sampled and measured the characteristic wave period τ , vertical phase and group velocities (v pz and v gz ) of the wave packets of zonal wind, using frequency windows corresponding to the period bands of 10 min-1 h, 30 min-2 h, 2-7 h, and 7-24 h, and wave number windows corresponding to vertical wave length bands of 0.7375-17.7 km, 1.475-17.7 km, 2.95-17.7 km, 5.9-17.7 km and 8.85-17.7 km. For the convenience of discussion, we listed the corresponding wave period band and vertical wave length band associated with each window in Table 4. It should be understood that due to the process of phase and group velocity tracing, the characteristic wave periods of the sampled wave packets could not cover short periods less than 20 min because of the time resolution of this data set (2.5 min), and the vertical wave length was limited to the range 1.5 km∼14 km (see Table 5) because of the height resolution (∼150 m) and the observation range (∼18 km).

Examples of velocity tracing analysis of specific wave packet
A partial Range-Time plot of (δV ) 2 of zonal velocity converted from the subset {u U (z, t)} and {u D (z, t)}, which were obtained using window T2Z1 (defined in Table 4), by the method of phase-and group-velocity tracing is shown in the upper-and lower-panel of Fig. 9, respectively. Both panels revealed that the strongest wave packet activity appeared in the height range around 15 km altitude, consistent with the heights of the strongest upward momentum flux (profile b,  upper panel of Fig. 2) and downward momentum flux (profile b, lower panel of Fig. 2). Also shown in Fig. 9 is the velocity measurement of the two most energetic wave packets (denoted by 1 and 2) in each panel. According to the method of phase-and group-velocity-tracing, the slope of the straight line along the middle patch (called phase line) equals the vertical phase velocity v pz ; while the slope of the straight line across the three neighboring patches (called energy line) equals the effective vertical group velocity v gz , and the horizontal distance between two consecutive phase lines equals one half of the characteristic wave period τ . The results of measurements were: v pz =−0.613 m/s, v gz =+0.148 m/s and τ =37.9 min for wave packet (1) (2) in the lower panel. Not only these four energetic wave packets, but also the absolute majority of the wave packets shown in Fig. 9 were type 1 wave packets as expected. The measured characteristic wave periods (≈39 min) were much closer to the period (30 min) corresponding to the high frequency end than the corresponding period (2 h) of the low frequency end of the window (T2Z1) which was used to do the wave separation. The same conclusion is also true in the wave length measurement. Such windowing effect had been explained in a previous paper (Kuo et al., 2003). The most interesting phenomenon is that the upward wave packet (1) (in upper panel) and the downward wave packet (1) (in lower panel) occurred simultaneously at the same height, indicating that these two wave packets might be created by one strong wave event. The statement is also true for wave packets (2) in both panels.
The most interesting example of the wave packet measurement is presented in Fig. 10. The wave packets in the upper panel were obtained from {u U (z, t)} using window T2Z4 (defined in Table 4), while the wave packets in the lower panel were obtained from {u D (z, t)} also using the window T2Z4. The results of the measurements are: v pz =−4.397 m/s, v gz =+4.474 m/s and τ =48.1 min for the sampled packet in the upper panel, and v pz =+5.32 m/s, v gz =−4.763 m/s and τ =44.2 min for the sampled packet in the lower panel. Again, the characteristic wave periods (wave lengths) were close to the period (wave length) corresponding to the high frequency (large wave number) end of the window. We also notice that these two wave packets occurred simultaneously at the same heights. The theoretical analysis using dispersion equation yielded that: the  upward moving wave packet (upper panel of Fig. 10) propagated in the north-west-north direction with azimuth angle of 346 degrees, and the downward moving wave packet (lower panel of Fig. 10) propagated in the opposite direction (southeast-south) with azimuth angle of 166 degrees. Their respective horizontal wave lengths and intrinsic periods were: 54.9 km and 44.9 min for the upward moving wave packet and 63.6 km and 46.8 min for the downward moving wave packet. So, these two oppositely propagating wave packets had nearly equal wave lengths and intrinsic periods. They were seemingly resulted from a parametric sub-harmonic instability, and its energy source might be related to the strong zonal wind jet.

Statistics of vertical group velocity of upward moving wave packets
We are interested in the vertical group velocities of only the upward moving wave packets in this study because only the upward flux of wave momentum and energy may possibly be transported to the middle and upper atmosphere. Figure 11 shows the scattering plot of characteristic vertical wave length vs. wave period for the wave packets obtained by windows T2Z1, T3Z1, T4Z1, T2Z4, T3Z4, T4Z4 and T1Z4. And the scattering plots of vertical group velocity vs. vertical phase velocity of the wave packets obtained by windows T2Z1, T2Z3, T2Z4 and T2Z5 were shown in Fig. 12. Among all the 1476 wave packets (listed in Table 5) we have analyzed, only 86 were type 2 wave packets (6%), and 94% of them were type 1 packets. This statistics confirms our presumption that the majority of wave packets are type 1. The  mean as well as the standard deviation for the measured values are summarized in Table 5, which clearly revealed that shorter period waves moved much faster than the longer period waves vertically on the average. The characteristic wave period, vertical wave length and group velocity of other wave packets can be estimated by interpolation or extrapolation from Table 5. The vertical propagation of gravity wave is known to be modulated by the background wind (Doppler shift effect) in such a way that it moves faster vertically when it propagates opposite to the mean wind. Even if the vertical group velocity of a wave packet is not constant on its way up to the upper atmosphere, we still found some useful information from Fig. 12: a wave packet originated from the lower atmosphere can reach ionosphere E-region (100 km altitude) in less than five hours if its mean vertical group velocity is larger than 5 m/s. Figure 12 clearly shows that there were many sampled wave packets with characteristic wave period ∼40 min moving faster than 5 m/s vertically. By extrapolation, we believe that there were many wave packets with periods ∼15 min moving much faster than 5 m/s vertically.

Summaries and discussions
We have studied the momentum flux of gravity wave by a technique slightly different from the conventional method. Conventionally, the vertical flux of zonal momentum is calculated by the left hand side of Eq. (6). Since u w =w u , u w can be interpreted either as a vertical flux of zonal momentum or a zonal flux of vertical momentum (Gossard and Hooke, 1975). If u w is positive, it can be interpreted either as an upward (downward) flux of eastward (westward) momentum or a eastward (westward) flux of upward (downward) momentum. Likewise, if u w is negative, it can be interpreted either as an upward (downward) flux of westward (eastward) momentum or an eastward (westward) flux of downward (upward) momentum. In contrast to the conventional method, a double Fourier transformation must be applied, before calculating the momentum flux by the method in this study, to separate the phase-upward propagating waves and the phase-downward propagating waves. We assumed that the former contribute to the downward flux of zonal momentum (2nd term of rhs of Eq. 6) and the later contribute to the upward flux of zonal momentum (1st term of the rhs of Eq. 6) because absolute majority of wave packets were type 1. Such assumption is necessary because we are unable to separate the group velocity upward waves from the group velocity downward waves directly. Then the direction of the horizontal momentum can be decided depending on the sign of the flux value. This process of analysis had led us to the following findings: During the quiet day, both upward flux and downward flux were dominated by eastward momentum (propagating along with the mean wind). It seems that the gravity waves prefer to propagate in the same direction of the mean wind under quiet condition, but we need more experiments to check whether this is a characteristic of nature. During the active days, the upward flux was dominated by westward momentum (propagating opposite to the mean wind) and the downward flux was dominated by eastward momentum (propagating along with the mean wind). And, two oppositely propagating wave packets were often observed to be simultaneously created at the height of maximum momentum flux (see Figs. 9,10). The fact that the amplification factor (17∼70) of higher frequency waves growing from quiet phase to active phase was much larger than the amplification factor (4∼7) of the lower frequency waves might be resulted from the wave breaking of larger scale waves. To support this conjecture, a pair of oppositely propagating wave packets (see Fig. 10) was found to be seemingly resulted from a parametric sub-harmonic instability, and its energy source might be related to the wind jet. A detailed analysis of the propagation directions of all the wave packets sampled in this study and the possible non-linear processes in the wind jet vicinity is our current subject of research. Group velocity measurement gave us a starting point to estimate the time required for a wave packet to propagate from lower stratosphere to upper atmosphere. Based on the measured characteristic wave period, vertical phase-and groupvelocities of a wave packet, its horizontal wavelength λ h and the azimuth angle φ of the propagation direction can be estimated by the dispersion equation and polarization equation. Then, the traveling time T of this wave packet from the starting point to the upper atmosphere (if this wave packet was not stopped by filtering, wave breaking or wave-wave interaction) can be calculated by numerical integration: T = dt, where dt=dz v gz (z) and the vertical group velocity v gz (z) at height z depends on the vertical profiles of the horizontal mean wind and the BV frequency. The scattering plot in Fig. 12 revealed that there were many wave packets with characteristic period ∼40 min having vertical group velocity larger than 5 m/s. They are good candidates to be further investigated. Kuo et al. (2007) had proved by the method of wave packet analysis that QP echoes were most likely resulted from the modulation by upward propagating gravity waves with periods less than 10 min. Then simultaneous observation of ionosphere E-region and lower atmosphere to track down the source of QP echoes directly may be meaningful.

Fig. A1.
There is some limitation on the wave packet analysis. The time resolution (2.5 min) limits the characteristic wave period of the sampled wave packets to be no smaller than 25 min, and the total observation range (∼18 km) limits the characteristic vertical wave length to be no larger than 14 km. Since the wave packet with smaller wave period and larger vertical wave length moves faster vertically, and they would be sooner than others to reach upper atmosphere (if they can reach before breaking or disappearing by energy dissipation). In our analysis, the very fast moving wave packets (e.g. wave packets with period ∼15 min) could not be sampled because of the limitation, but they can be estimated indirectly through extrapolation from the sampled packets.

Error estimation of velocity measurement by dual beam method
Equation (2a) and (2b) in the text represent respectively the average (or, linear interpolation) of the eastward velocity u and vertical velocity w obtained by eastward and westward oblique beams at the same height (but different locations). We may approximately assign such value as the velocity at the symmetric center C of the two locations (E and W) detected by the two beams at the same height z (see Fig. A1). This symmetric center is located along the vertical axis on the earth. Let's denote the coordinate of point C by (0, z), point E by (x, z) and point W by (−x, z). Then consider a gravity wave at height z with zonal wavelength λ x . If the Sketch A1 Sketch A2 Then, the zonal velocities contributed by this wave at point E and W can be represented by Eqs. (A2) and (A3), respectively, where η=2π x λ x is the phase difference between points C and E. By Taylor's series expansion of cos (η+ η) and cos (η− η), we obtain the average of u (E) and u (W ) to be, The percentage difference between the true velocity u (C) at point C and the velocity u obtained by the dual beam method is, Some numerical values yielded by Eq. (A5) are listed in Table A1, which reveals that if the zonal wavelength λ x is a factor of 8 larger than the horizontal distance 2x between points E and W, the percentage error resulted from the dual beam method will be less than 8%, and the error decreases with increasing value of the ratio λ x x. The geometrical meaning of this numerical calculation can be further understood by a sketch of a sinusoidal wave (Fig. A2). As Fig. A2 shows, when a sinusoidal curve (of one wavelength) is divided into 8 equidistance sections, each section is approximately a straight line. Therefore, the linear interpolation is a reasonable approximation if λ x 2x>8 for sinusoidal wave. The above statements are equally true for the dual beam measurement along any azimuthal direction. Namely, if the horizontal wavelength measured along the azimuthal direction of the relevant beam pair is a factor of 8 larger than the horizontal spacing between the radar volumes of the 2 opposite beams, then the errors are less than 8%.
In the case of this study, x=z· tan 10 • , so the phase difference η increases with increasing height. At the altitude z=15 km, the distance between the two observation locations E and W is 2x=5.25 km. Imagine a wave packet with characteristic horizontal wave length λ h =30 km to propagate at this height in the direction with azimuth angle of ϕ=45 • . The characteristic zonal wavelength of this imagined wave packet is λ x =λ h / sin ϕ=42.4 km, satisfying the condition of λ x 2x>8. At altitudes of 2, 5, 10, 15, 20 and 25 km, the useful minimum wavelength λ min (required by dual beam method) along the azimuthal direction of the relevant beam pair can be estimated by Eq. (A6), and the results are listed in Table A2.
To check the validity of the dual beam method for the calculation of momentum flux in this study, we had analyzed the wave packets from the 30 min-2 h, 2-7 h and 7-24 h period bands of this study by the dispersion equation to obtain their characteristic intrinsic frequencies and horizontal wave lengths λ h , under the assumption that the mean zonal wind velocity was 37 m/s and the mean meridional wind velocity was 8 m/s. The mean characteristic horizontal wave lengths (regardless of their propagation direction) were found to be 62 km, 264 km, and 1028 km for the sampled wave packets from the 30 min-2 h period band, 2-7 h period band and 7-24 h period band, respectively. They were large enough to confirm statistically that the absolute majority of the sampled wave packets satisfied the condition of λ x 2x>8. So, the dual beam method (Eqs. 2a and 2b in the text) is valid to determine the gravity wave values for u and w in the present case of study.

B1 Separation of upward waves and downward waves
Consider an oscillation profile of the meridional (zonal) wind velocity of atmospheric motion ψ (z, t) as a function of height z and time t. First we make finite Fourier transformation of ψ (z, t) over the height coordinate z at each fixed time t to obtain the coefficients A n (t) and B n (t) of the function cos (nz) and the function sin (nz), respectively: ψ (z, t) = n {A n (t) · cos (nz) + B n (t) · sin (nz) } Then, finite Fourier transformations are made on each A n (t) and B n (t) over time t: A n (t) = σ A (1) nσ · cos (σ t) + B (1) nσ · sin (σ t) (B2a) B n (t) = σ A (2) nσ · cos (σ t) + B (2) nσ · sin (σ t) , where n and σ represent the vertical wave-number and the observed frequency of the component wave respectively. Substituting Eqs. (B2a) and (B2b) nσ · cos (σ t + nz) + B (1) nσ + A (2) nσ · sin (σ t + nz) .
The first term in Eq. (B3) represents the contributions from all the waves with upward phase velocities, and the second term comes from the contributions of all the waves with downward phase velocities. Thus, the original oscillation profile ψ (z, t) can be divided into phase-upward profile U (z, t) and phase-downward profile D (z, t): nσ · cos (σ t − nz) nσ · cos (σ t + nz) + B (1) nσ + A (2) nσ · sin (σ t + nz) .
Notice that the summations are made over selected frequency window and wave number window case by case. We shall analyze U (z, t) and D (z, t) instead of the original data set ψ (z, t). A Range-Time plot of the local power of fluctuation of (z, t) (representing either U (z, t) or D (z, t)) can be obtained by following procedure (Kuo et al., 1998): at each time step t i =i t and height z k =k z, a time series of N data points (Nis chosen to be 15 in this paper) is defined by, k = (z k , t ) , =i− N 2 , i− N 2 +1, ..., i+ N 2 .
First, to remove contribution from very long period wave, a linear trend of this time series (Eq. B5), y=at+b, is obtained by the method of least square fitting, then a new time series (Eq. B6) is obtained from Eq. (A5) by subtracting its linear trend: Then, the time series (Eq. B6) is Fourier analyzed to obtain the power P ki of the first harmonic. Considering P ki as an index of the strength of the wave activity at height z k and time t i , we plot the matrix [P ki ] as an Range-Time plot for investigation (We use contour plot to represent the Range-Time plot in this paper). Numerical simulation study (Kuo et al., 1998) showed that when a wave packet is propagating, packets with high regularity would exist in the Range-Time plot. It was shown that the same phase points of the characteristic wave component of the wave packet (patch) form a straight line called phase trajectory (phase line), while the centers of the neighboring packets (patches) form another straight line called energy trajectory (energy line). It was also noticed that the former determines the direction of the phase progression and the latter constitutes the direction of energy propagation. Therefore the slope of the phase trajectory is just the phase velocity and the slope of the energy trajectory is the group velocity of the wave packet. Also, the horizontal distance between two consecutive packets (corresponding to the peaks and the dips in the oscillation profile) equals one half of the period of the characteristic wave component of the packet.