A new approach to momentum flux determinations using SKiYMET meteor radars

Abstract. The current primary radar method for determination of atmospheric momentum fluxes relies on multiple beam studies, usually using oppositely directed coplanar beams. Generally VHF and MF radars are used, and meteor radars have never been successfully employed. In this paper we introduce a new procedure that can be used for determination of gravity wave fluxes down to time scales of 2-3h, using the SKiYMET meteor radars. The method avoids the need for beam forming, and allows simultaneous determination of the three components of the wind averaged over the radar volume, as well as the variance and flux components , where refers to the fluctuating eastward wind, refers to the fluctuating northward wind, and refers to the fluctuating vertical wind. Data from radars in New Mexico and Resolute Bay are used to illustrate the data quality, and demonstrate theoretically expected seasonal forcing. Keywords. Meteorology and atmospheric dynamics (Middle atmosphere dynamics; Waves and tides; Climatology)


Introduction
Determinations of gravity wave momentum fluxes are a crucial requirment for understanding middle atmosphere dynamics and energetics.A few measurements of this parameter have been made (Reid and Vincent, 1987;Fritts and Vincent, 1987;Fritts et al., 1992;Murphy and Vincent, 1998;Nakamura et al.,1993Nakamura et al., , 1996;;Tsuda, et al., 1990;Thorsen et al., 1997), but in the context of a full understanding of global morphology, studies have really only just begun.Previous measurements have been primarily made by radar (VHF and MF), but VHF radars have limited seasonal and height coverage, and MF radars with narrow beams are relatively rare.The need exists for a larger suite of instruments, better distributed on a global scale, to supply geographic, seasonal and annual variability of momentum flux parameters.Meteor radars are well suited to fill this niche, since a large number Correspondence to: W. K. Hocking (whocking@uwo.ca) have been developed in recent years.Currently there are almost 30 such radars distributed world-wide, with locations as high as 78 N, to the equator, and on to the south pole, and more are under development.The advantage of this suite of radars is that most are similar in design (typically being SKiYMET radars; Hocking et al., 2001a), making intercomparisons simple and reliable.However, these radars are not narrow-beam systems, but use all-sky coverage.Individual meteors are located to an accuracy of ±1.5 • in zenith and azimuth, and in principle it should be possible to choose a subset of these meteors to emulate a dual-beam system.However, meteor count rates are generally too low to make this possible.Hence to date meteor radars have not been able to contribute to momentum flux studies.
The principle of the primary current momentum flux method relies on two beams orientated in the same vertical plane, but tilted at an angle θ in opposite directions (Reid and Vincent, 1987).For simplicity, assume that the two beams are in the +x and −x directions.In one beam the instantaneous radial velocity measured is (U +u ) sinθ +(W +w ) cosθ , where U is the mean eastward wind, u is the fluctuating eastward component, W is the mean vertical wind and w the fluctuating vertical wind.In the other beam the radial velocity is −(U +u ) sinθ +(W +w ) cosθ .The difference in variance in the two beams is then 4u w cosθ sinθ , or 2u w sin(2θ ), where u w is the flux of horizontal fluctuating momentum in the vertical direction (and also of the vertical fluctuating momentum in the horizontal direction).In other words where v rad1 2 and v rad1 2 are mean square fluctuating radial velocities in the two beams.Similarly we can determine v w using antenna beams pointing to the north and south.An alternative method for momentum flux calculations, using MF broad beam interferometric radars, was introduced by Thorsen et al. (1997).
In this paper, we wish to generalize the dual-beam formulation to deal with cases where narrow radar-beams do not exist, but in which the location of individual scatterers is known to good accuracy.We will show that the formulation that we develop is quite general, and the traditional dual beam method is a special case of this formulation.The new method has parallels with the medium-frequency interferometric method developed by Thorsen et al. (1997), but is different in the details of application, and is new because for the first time meteor echoes are used to measure momentum fluxes.

Method outline
When mean winds are determined by meteor methods, an all-sky fit of all radial velocities is performed, minimizing the quantity , where v radm is the mean radial velocity expected if the winds were uniform in a horizontal plane with values in cartesian coordinates of (U, V , W ). The summation is over all detected positions.We have not added subscripts or a summation index, but these are implied.For a meteor at position (θ, φ) in spherical coordinates, where θ is the angle from zenith and φ is the azimuthal angle anticlockwise form due east, v radm =U sinθ cosφ+V sinθ sinφ+W cosθ .However, in reality the wind is not uniform, and the measured radial velocity v rad usually differs from the value of v radm .
Recognizing that these deviations between v rad and v radm represent true wind variability (primarily due to gravity waves), we now propose to mimimize the quantity We assume first that we have performed a fit of the mean wind, as described above, and at any meteor position we have removed the radial velocity due to the mean wind.Here v rad represents the difference in radial wind between the measured value and the value expected from the knowledge of the mean wind (v rad =v rad −v radm , where we will refer to v radm as the "model" values).Application of Eq. (2) amounts to optimizing the similarity between the measured and modeled variances of radial velocity as a function of time and position.We write that the model radial velocity at position (θ, φ), for assumed fluctuating velocities u , v and w , is where the fluctuating velocity components are assumed to be due to wave and turbulent motions.Squaring this term and substituting into Eq.( 2) means that we must mimimize where the summation is over all detected meteor positions (θ, φ) within a user prescribed height and time interval.Again, we have not specifically added indices.
To minimize , we partially differentiate this function with respect to u 2 , v 2 , w 2 , u v , u w and v w , and set each derivative to zero.For example, if we differentiate with respect to u 2 , we obtain 2 [(v rad ) 2 −(u 2 sin 2 θ cos 2 φ+v 2 sin 2 θ sin 2 φ+w 2 cos 2 θ +2u v sin 2 θ cosφ sinφ+2u w sinθ cosθ cosφ +2v w sinθ cosθ sinφ)] sin 2 θ cos 2 φ=0. (5) Similarly we may differentiate with respect to all 6 parameters.The final result is a matrix equation of the following form, assuming that the parameters u 2 , v 2 etc. are all uniform across the field of view: This equation can readily be inverted to produce an estimate for the 6 parameters u 2 , v 2 , w 2 , u v , u w and v w .This is the principle of the method.
As a check on the validity of this equation, it is relatively easy to show that the traditional dual beam method is a special case of this equation.To see this, imagine that all targets occur in the φ=0 or φ=180 • plane, so sinφ is zero for all targets.In addition, assume all targets occur at a zenith angle of θ , where this time θ is a singular value.Targets may then only occur at (θ, φ) values of (θ, 0) and (θ, 180 • ).The summations in the above equation then become sums over two possible angles.Then multiplying the first row of the first matrix by the first column matrix, and recognizing that cos 2 φ=1 for all cases, gives Dividing through by sin 2 θ gives Assuming N scatters in the beam at φ=0 and M in the beam at φ=180 • , separating out the terms in to cases of φ=0 and φ=180 • , denoting the mean square radial velocity fluctuation in each beam as v rad0 2 and v radπ 2 , and recognizing that the sum is just the product of the number of points and the mean, gives (N+M)u 2 sin 2 θ+(N+M)w 2 cos 2 θ+2Nu w sinθ cosθ cos(0) where overbars indicate averages.Then we may write The second line of the matrix equation contains sinφ in all terms, which is zero for our dual-beam configuration, so is of no use.The third line of the matrix multiplied by the first column matrix gives identical information to Eq. ( 7).The fourth line of of the matrix contains sinφ in all terms and gives no information.Likewise for the last line.Thus the remaining useful equation is found as the fifth line multiplied by the column matrix, or Dividing through by sinθ cosθ gives Again considering N scatterers in the first beam and M in the second, and substituting for the sums with the products of the appropriate means and numbers of points, gives 2Nu 2 sin 2 θ cos(0)+2Mu 2 sin 2 θ cos(180) +2Nw 2 cos 2 θ cos(0)+2Mw 2 cos 2 θ cos(180) +4Nu w sinθ cosθ+4Mu w sinθ cosθ=2Nv rad0 2 cos(0) or We can now substitute for (u 2 sin 2 θ +w 2 cos 2 θ ) from (9) to give which is just Eq. ( 1).
A similar logic may be followed for the case of dual beams at φ=90 and −90 • which will give the standard dual-beam result for beams pointing to the north and south.Clearly our matrix equation is a general equation which encompasses existing techniques.Henceforth we will use the matrix formulation for all our analyses.

Application
In order to test this new theory, we use data from two sitesa mid-latitude site and a polar site.The polar site is at Resolute Bay, Nunavut, Canada, at coordinates of 75N and 95W.The details about this radar have been reported in Hocking et al. (2001b).The radar operates at 51.5 MHz and shares its resources between various recording modes, and runs in meteor mode for typically 60% of the time.Four receiver antennas are used to determine meteor locations by interferometry.In 2001 it was upgraded to employ four separate receivers, rather than multiplexing the signals as had been done previously, and this increased the data rate.The mid-latitude site is at Socorro, New Mexico, in the USA.It is located at latitude 34 N and longitude 107 W, and is a standard uncrossed SKiYMET radar (Hocking et al., 2001a), operating at 35.24 MHz.It uses five receivers for reception.In each case, meteors are located to an accuracy of about ±1.5 • in angle, and an accuracy of about 3 km in height.Position, amplitude, radial velocity and decay time are determined for each meteor, as well as other parameters related to reliability and ambiguity.Only unambiguously located meteors have been used in the analyses presented in this paper.Typically the Resolute Bay radar detects about 6000 meteors per day in summer, and about 1500-2000 in winter.At Socorro, typical count rates are of the order of 3500 meteors per day throughout the year, with a slight decrease in January and February.Occurrence of noise can diminish these rates on any one day.
Good quality data exist for the Resolute Bay site since July 2001, when the receivers were upgraded, but with gaps due to downtime and system failures.The Socorro radar was established at that site in April 2002, having been previously moved from Starfire Optical Range in Albuquerque, NM.Data quality in 2002 was variable, due to an old computer, but since late 2002 a new computer has been installed and downtime has been minimal.
In order to test our theory, we have applied our equations to data from both these sites.For the Socorro radar, peak meteor count rates occur at zenith angles of 50 to 60 • , and for the Resolute Bay radar, wich uses a transmitter beam which concentrates more on meteors closer to overhead, they occur at about 50 • .However, because of the importance of the vertical velocity in these determinations, we have excluded data beyond 45 • from zenith in our analyses.We have also excluded meteors detected at angles closer to zenith than 15 • , as is normal in meteor studies.Standard calculations of momentum fluxes using dual-beam radars (e.g.Reid and Vincent, 1987) employ off-vertical tilts of 10 to 15 • , where the ratio of relative contributions from vertical and horizontal velocities to the radial valocity are approximately in the ratio 4:1.At 30 • the ratio is 1.7:1, and at 45 • it is 1:1.Thus meteor radars are less sensitive to vertical velocity contributions, but they can compensate for this by achieving higher count rates than MF and VHF radars.Ideally, we should concentrate the meteor selection to zenith angles even closer to zero, but then count rates fall to unacceptable values.A new radar is under design which will allow optimal selection closer to vertical, Socorro, NM. 15 March 2004-19 March 2004 (incl.)but for now we wish to evaluate the capabilities of standard SKIYMET meteor radars.
Because the parameters we are deriving are second order terms, considerable care is needed in order to ascertain that the results are reliable.Whereas typically 7 to 9 meteors are required in any height-time bin for reliable estimation of a mean wind, many more will be needed to make reliable estimates of the momentum flux parameters.Our first task has been to determine acceptable limiting parameters.The upper two graphs in Fig. 1 shows determinations of mean winds and the northward velocity variance for a 5-day period at Socorro, using 1-h data bins at a range of 88 km.The bins actually cover 90 min, covering 45 min before and 45 min after the nominal time, and the nominal time is shifted in steps of 1 h.The mean winds show a clear tidal signature.In the second graph, the shaded region represents values less than zero.It is clear that on occasions the values of v 2 fall below zero, which is unphysical, but is not precluded by the equations.The points labelled A, B, C,...H corespond to cases of extremely large variances or negative variances, and all corespond to low count rates.This indicates breakdown of the method in these cases.Investigations have shown that such negative values, and unusually large values, correspond to cases where the number of points used in the analysis was less than 30.The "box-car" type line shown within the shaded region highlights cases where the number of counts is greater than and less than 30.When this function is high, it indicates over 30 points, and when it is low, it indicates less than 30 points.The lower two figures show momentum fluxes, but cases of less than 30 points, and negative values of u 2 and v 2 (which usually coincide with cases of less than 30 points) have been removed.Some values of u w are clearly revealed, but the acceptance rate is low.Not surprisingly, we also found that less erratic values of momentum flux can be determined when the positions of meteors in any time-height bin are uniformly distributed around the sky.
Figure 2 also shows data from Socorro, for the same time of year as in Fig. 1 (March, 2004), but in this case we have used three-hourly data bins, shifted in steps of 2 h.In this case the number of meteors per accepted bin always exceeded 30, and u 2 and v 2 were always positive.We, therefore, consider that the data are more reliable, and henceforth will use such averaging.The reliability is still of course unproven, but will be addressed in the coming paragaphs.
In order to ascertain the importance of knowledge about the vertical velocities, we repeated the above procedures for the cases that we evaluated the mean vertical wind, and the case that we assumed that the mean vertical wind was zero.Our observations showed that the vertical wind variance was affected by this assumption, but the estimates of u 2 , v 2 , u v , u w , and v w were fairly consistent for the two cases.We have adopted the practice of setting w to zero for our subsequent calculations.
The variances and fluxes shown in Fig. 2 show strong variability, with frequent occurrences of strong episodic increases.Typical values of u w are of the order or 10 to 50 m 2 s −2 , with occasional values as high as 200.These values appear to be higher than previously reported values, but Fritts and Vincent (1987), Fig. 14, shows values as high as 30 m 2 s −2 , and these are for a 3-day average.In addition, MF radars are known to underestimate the winds above 90 km altitude, so their values could have been in fact larger.Fritts and Vincent (1987) discussed the fact that intermittent episodic values of high momentum flux should be an expected feature of mesospheric forcing.Fritts and Alexander (2003), Sect.8.1.2,note that values as high as 30-60 m 2 s −2 can easily be expected.The values presented in our own Fig. 2 are therefore large, but not unreasonably so.
In order to further check the validity of our measurements, we have formed monthly averages for all the available recent data from Resolute Bay and Socorro.Results are shown in Figs. 3 and 4. For Socorro (Fig. 3), data used covers the period April 2002 to January 2005.Data were averaged in two-monthly bins in order to make the graphs more compact.The Resolute Bay data cover the period from July 2001 to September 2004.Data were only accepted if count rates are suitable for more than 20 days per month.In some cases this limits us to 2 months of data, and in other cases data were available every year.Error bars are also shown, and are typically ±3 m 2 s −2 at 88 and 91 km.The data at 82 and 85 km had slightly larger errors, in part because of slightly lower meteor count rates, but more importantly because the mo- mentum fluxes showed significant variability on time scales of 10 to 20 days, presumably due to planetary-wave filtering at lower heights.In both Figs. 3 and 4, a least-squares fit straight line is plotted on each graph, for the data points 82 to 91 km.The 94 km data are included, but have not been used in the straight-line fit.The significance of these lines is discussed in the next section.

Discussion
Current theories of gravity wave breaking in the mesosphere ascribe the reversal of mean winds above 80 km to deposition of momentum by gravity waves (Lindzen, 1981;Holton, 1983;Fritts and Alexander, 2003).It is proposed that if, as gravity waves propagate upwards, the vertical flux of horizontal zonal momentum diminishes, then this lost momentum flux must be transferred to the mean flow.Thus a decreasing vertical flux of zonal momentum should cause an acceleration of the mean flow, according to ∂U ∂t u w =− ∂u w ∂z , where the subscript u w on the left hand side indicates that this refers only to the component of the acceleration specifically due to momentum flux deposition.Because of filtering by the mean flow at lower heights, an asymmetry in gravity wave propagation directions ensues, with the consequence that during summer the vertical flux of zonal momentum should decrease with increasing height, causing a mean flow reversal above 80 km altitude.As a result the expected radiatively balanced westward flow becomes a dynamically/thermally balanced eastward flow.In winter the reverse should occur.Both Figs. 3 and 4 are consistent with these scenarios, with increasing momentum fluxes as a function of height in winter (December to February), and decreasing momentum fluxes as a function of height in summer.At Socorro the momentum fluxes show a tendency to decrease as a function of height as early as April (Fig. 3), and this continues into September.In October the fluxes are constant on average with increasing height, and the slope reversal begins in late November and early December.At Resolute Bay the momentum fluxes do not appear to show a decrease with increasing height until well into May, but by July/August the rate of decrease with increasing height is very strong indeed.September and October are variable, and in late November and December the momentum fluxes u w increase again with increasing height.Typical magnitudes of the monthly mean values of the mean flow acceleration are of the order of 1 m 2 s −2 km −1 , with values at Resolute Bay in July and August reaching as high as 2.3 m 2 s −2 km −1 .These may be converted to accelerations in m s −1 day −1 by multiplying by 86.4 (number of seconds in one day divided by 1000).Hence typical accelerations are of the order of ±80 to 100 m s −1 day −1 , and the acceleration over Resolute Bay in July and August can reach monthly average values of almost 200 m s −1 day −1 .These values are not inconsistent with values summarized by Fritts and Alexander (2003), Sect. 8.1.2,where values of 10 to 70 m s −1 day −1 are proposed for midlatitude sites.Measurements have not yet been made at polar sites, so the results discussed here for Resolute Bay are especially timely and important.
The above results are therefore indicative that meteor radars can indeed make useful contributions to measurements of momentum fluxes.The fact that a world-wide network of SKiYMET radars already exists makes this ability especially appealing.However, we also recognize that this is a very preliminary study, and no doubt the abilities of meteor radars in this regard can be refined significantly.We have already discussed the idea of designing a special radar which is more sensitive to meteors from overhead.However, we still advocate that all meteor radars should have interferometric capability, because the huge dynamic range of meteor echo cross-sections means that they can easily be detected in the sidelobes of monostatic and non-interferometric radars.
One further issue that needs to be considered is the relative contributions of spatial and temporal variability.If a meteor radar detects meteors out to 45 • from zenith, this means that the scattering volume of the radar covers a width of typically 180 km.Thus gravity waves with horizontal wavelengths of typically 180 km can contribute to the variance.If such a wave has a vertical wavelength of typically 5 km, then the period of such a wave is approximately 3 h.Hence it does not make sense to calculate the variance over time scales of less than 3 h, since even if the averaging was restricted to periods of say 90 min, the spatial variability would mean that waves with periods of 3 h would still contribute.Thus our choice of a 3 h averaging period and a limit of 45 • zenith angle are a relatively good combination.If averaging intervals of say 90 min are needed, then zenith angles should be restricted to values less than 27 • , so that the horizontal wavelengths and averaging times are commensurate.
It can be argued that traditional momentum flux determinations might be superior in this regard, in that the scattering volumes are separated by only a few tens of kilometers.However, it needs to be also noted that the atmosphere will move during an averaging interval, so that if an integrating period of 3 h is used, and the mean wind is 20 m s −1 , then an atmospheric volume of width 216 km will drift through the beams in that time.Hence the traditional methods are not immune to the effects of interplay between temporal and spatial resolution either, and for integrating periods greater than 3 h, the effects are in fact similar for both meteor and traditional methods.In cases where averaging times as high as 8 h are used (e.g.Fritts and Vincent, 1987;Murphy and Vincent, 1998) the spatial variations can be the main contributor to variability.We would argue that meteor-derived momentum fluxes offer a great deal in the area of momentum flux studies, and the relative large numbers, broad geographic spread, and ease of construction of SKiYMET radars make them excellent candidates for determinations of momentum fluxes.

Conclusions
The theory for determination of gravity wave momentum fluxes in the middle atmosphere using meteor techniques has been outlined.While meteor radars cannot actually image the gravity waves, they can easily measure the radial velocity variances associated with them, and as long as it is accepted that gravity waves are the main cause of this variance, then meteor radars may be used to determine fluxes.Values measured are consistent with earlier measurements by dual-beam methods, although on occasion can be somewhat higher.We believe that these larger values are physically real.Monthly mean forcings measured at Socorro (NM) and Resolute Bay (Canada) show values consistent in form with current theories of gravity wave forcing of the mean circulation, and typical mean flow accelerations are of the order of 50 to 100 m s −1 day −1 during summer and winter (although opposite in sign).Mean summertime forcings at polar sites can be even stronger, reaching as high as 200 m s −1 day −1 .Plans for a new class of SKiYMET meteor radar are underway which will optimize the capability of these radars for momentum flux measurements.

Fig. 1 .
Fig. 1.Examples of mean winds (upper left), northward gravity wave variance (upper right), northward flux of zonal momentum (and also the eastward flux of meridional momentum) (lower left) and vertical flux of meridional momentum (lower right) for Socorro NM, using 90 min data sets.See text for details.

Fig. 2 .Fig. 3 .
Fig. 2. Socorro data using 3-h bins and 2-h steps.Slightly different parameters are shown compared to Fig. 1, just to illustrate that all flux terms are available.Quality has improved significantly compared to Fig. 1.

Fig. 4 .
Fig. 4. Height profiles of mean values of u w for Resolute Bay, for the period July 2001 to September 2004, plotted as two-monthly averages.