Mesospheric gravity wave momentum flux estimation using hybrid Doppler interferometry

Mesospheric gravity wave (GW) momentum flux estimates using data from multibeam Buckland Park MF radar (34.6 S, 138.5 E) experiments (conducted from July 1997 to June 1998) are presented. On transmission, five Doppler beams were symmetrically steered about the zenith (one zenith beam and four off-zenith beams in the cardinal directions). The received beams were analysed with hybrid Doppler interferometry (HDI) (Holdsworth and Reid, 1998), principally to determine the radial velocities of the effective scattering centres illuminated by the radar. The methodology of Thorsen et al. (1997), later re-introduced by Hocking (2005) and since extensively applied to meteor radar returns, was used to estimate components of Reynolds stress due to propagating GWs and/or turbulence in the radar resolution volume. Physically reasonable momentum flux estimates are derived from the Reynolds stress components, which are also verified using a simple radar model incorporating GW-induced wind perturbations. On the basis of these results, we recommend the intercomparison of momentum flux estimates between co-located meteor radars and verticalbeam interferometric MF radars. It is envisaged that such intercomparisons will assist with the clarification of recent concerns (e.g. Vincent et al., 2010) of the accuracy of the meteor radar technique.


Introduction
There has recently been particular interest in the use of specular returns from all-sky interferometric meteor radar to measure the gravity wave (GW)-driven vertical fluxes of horizontal momentum (herein momentum fluxes) in the mesosphere-lower thermosphere/ionosphere (MLT/I; ∼ 80-100 km altitude) (e.g.Antonita et al., 2008;Clemesha and Batista, 2008;Beldon andMitchell, 2009, 2010;Clemesha et al., 2009;Fritts et al., 2010aFritts et al., , b, 2012a, b;, b;Vincent et al., 2010;Placke et al., 2011aPlacke et al., , b, 2014Placke et al., , 2015;;Andrioli et al., 2013aAndrioli et al., , b, 2015;;Liu et al., 2013;de Wit et al., 2014bde Wit et al., , a, 2016;;Matsumoto et al., 2016;Riggin et al., 2016).This has largely arisen from a need to obtain improved spatial coverage in the parameterization of GWs and their associated momentum transport in climate models of the whole atmosphere (e.g.Kim et al., 2003;Ern et al., 2011), and the suggestion (Hocking, 2005) that such measurements can be made accurately in the MLT/I (with integration times of the order of 2 months) using relatively low-cost commercial meteor radar systems.Nevertheless, there are concerns over the accuracy and precision of the momentum flux estimates from this technique (e.g.Vincent et al., 2010).In this paper, a previously established technique which makes use of partial reflections from the mesosphere at medium frequency (MF) is re-visited and contrasted with the meteor technique.Our aim in doing this has been to determine if interferometric MF radars, in particular those which have a small antenna aperture and only transmit a vertical beam, are viable candidates for verifying momentum flux estimates from meteor radars.
Direct measurements of momentum fluxes in the MLT/I were pioneered by Vincent and Reid (1983), in a study that Published by Copernicus Publications on behalf of the European Geosciences Union.
utilized partial reflection returns from the Buckland Park MF radar (34.6 • S, 138.5 • E).Their experiment consisted of transmitting a broad, vertically directed beam, and applying Doppler beam steering (DBS) (Woodman and Guillen, 1974) to narrower, fixed receive beams, offset from the zenith in the cardinal (initially east and west) directions (herein referred to as a "complementary" beam arrangement).The momentum fluxes were then estimated from the difference in the beams' mean square radial velocities.Similar approaches have since been applied to the same system (e.g.Fritts and Vincent, 1987;Reid and Vincent, 1987;Murphy andVincent, 1993, 1998) and to other high-frequency (HF) and very high-frequency (VHF) radars at various sites (e.g.Fukao et al., 1988;Reid et al., 1988;Fritts and Yuan, 1989;Fritts et al., 1990Fritts et al., , 1992;;Sato, 1990Sato, , 1993Sato, , 1994;;Tsuda et al., 1990;Wang andFritts, 1990, 1991;Hitchman et al., 1992;Nakamura et al., 1993;Murayama et al., 1994;Placke et al., 2014Placke et al., , 2015;;Riggin et al., 2016).It was obvious in some of these studies, especially those involving radars with relatively broad ( 3 • ) transmit and receive beams, that the aspect sensitivity (or Bragg anisotropy; Muschinski et al., 2005) of the partially reflecting scatterers illuminated by the transmit beam needed to be measured and/or accounted for, or else the apparent receive beam zenith angles would overestimate the true values (see, e.g., Reid and Vincent (1987); Murphy and Vincent (1993) for two such approaches).Thorsen et al. (1997) introduced an extension of the Vincent and Reid (1983) approach for radars with an interferometric capability, which accounted for the "brightness distribution" (the normalized angular and Doppler-frequency power spectral density; see, e.g., Woodman, 1997) in the radar receive beam(s).The authors applied this approach to a broad, vertically transmitted MF radar beam and relied on geophysical variability in the brightness distribution to obtain a sufficient number of radial velocity-pointing direction pairs to solve for the mean winds and wind covariances.An important assumption they made was that the "true" measured fluctuations in the weighted angle of arrival exceeded those due to statistical estimation errors.Despite the apparent validity of this assumption for the radar utilized, and the retrieval of momentum fluxes that appeared to be physically reasonable, no such studies incorporating interferometric radar techniques at MF have since been published.Hocking (2005) later demonstrated the application of the Thorsen et al. (1997) technique to returns from specular meteor echoes, for estimation of momentum fluxes.There have been a number of concerns raised over the accuracy and precision of the estimates that have since been reported.Those concerns of relevance to this paper, which are related only to the wind field and scatterer location characteristics (and are valid for both the meteor and partial reflection approaches), may be summarized as follows: 1.The assumption of statistical stationarity of the wind and wave field over the volume spanned by the scatter-ers.As discussed by Reid (1987), this assumption is independent of the beam configuration (or brightness distribution) used to sample the wind field.Additionally, if the beam configuration is anything other than complementary, the momentum flux estimates will only converge to the true value for horizontal wave scales much larger than the beam separation.This is clearly more problematic for meteor observations, where angle-of-arrival (i.e. the "effective receive beams" of the radar) distributions peak at large offzenith angles (typically 40-50 • ) and are often asymmetric in azimuth for integration times less than a few hours.An example of an approach to alleviate the former effect can be found in the design of the SAAMER and DrAAMER systems (Fritts et al., 2010b(Fritts et al., , 2012b)), in which a larger transmitting aperture was used to increase the signal-to-noise (SNR) of returns (and hence number of detections) at small zenith angles.While reducing the total number of detections, such an approach also has the advantage of increasing the relative contribution of the perturbation component of the vertical wind to the radial velocity, as well as reducing errors in meteor height estimates.
2. The required integration time for statistical significance of the estimates.Kudeki and Franke (1998) showed that, for a perfect dual-beam "anemometer", at least 16 days of integration is required to reliably estimate a momentum flux in the stratosphere, if one assumes that the flux represents a fraction of around 1 % of the mean horizontal and vertical variance.Accounting for the effects of measurement noise and finite spatial correlation of the wind fluctuations (as would be relevant in the case of using a radar sampling multiple, separated volumes), Thorsen et al. (2000) came to a similar conclusion.Also, in their modelling assessment of the meteor technique, Vincent et al. (2010) argued that integration times of over 1 month are necessary for accurate momentum flux estimation, in approximate agreement with the findings of Kudeki and Franke (1998) and Thorsen et al. (2000) for the dual-beam technique.
A contrary argument was put forward by Fritts et al. (2012a).They independently assessed the meteor technique, and argued that the large-amplitude monochromatic gravity waves often observed in the MLT/I region would result in stronger correlations between component velocity fluctuations (i.e.larger ratios of momentum flux to mean variance) than Kudeki and Franke (1998) assumed, and that the required integration time would be reduced from 16 days by around a factor of 20.
Additionally, Riggin et al. (2016) suggested that the use of longer integration times may reduce the measured correlation between vertical and horizontal mo-tions, and hence lead to systematic underestimation of momentum fluxes.
3. The contamination of momentum flux estimates by temporal wind shear (leading to overestimation of the true fluxes).Some authors (e.g.Fritts et al., 2010bFritts et al., , 2012a, b) , b) have opted to remove the temporal shear imposed by tides and planetary waves by fitting these components to the Cartesian winds, and then subtracting their radial projection from the radial velocities prior to momentum flux estimation.Andrioli et al. (2013aAndrioli et al. ( , b, 2015) ) also attempted to remove the same shears by fitting them, directly evaluating the contribution of the fitted components to the momentum fluxes, and then subtracting these from those estimated previously.
This paper presents the application of hybrid Doppler interferometry (HDI) (Holdsworth and Reid, 1998) and the Thorsen et al. (1997) estimator to momentum flux determination from the large-aperture Buckland Park MF radar.In HDI, Doppler beams are created and steered in hardware, and multiple-receiver interferometry is used on reception to form the actual receive beams.Therefore, the paper represents the first application of the Thorsen et al. (1997) estimator to an MF radar with beam steering applied on transmission, as well as on reception.The use of off-zenith transmission in HDI allows for greater beam directivity and hence SNR in the received off-zenith beams (for a given receiver array geometry).Additionally, the HDI analysis yields a radial velocity estimate at an effective beam position (EBP) corresponding to the peak of a fitted Gaussian brightness distribution.This means that the aspect sensitivity in the partially reflecting scattering surfaces is directly accounted for in the velocity estimation.
Section 2 describes the experiment configuration utilized on the Buckland Park MF in this study, and the array characteristics at the time (transmit/receive polar diagrams and data acceptance rate).Section 3 evaluates the ability of the radar to determine momentum fluxes by simulating the radar's sampling of a model GW-perturbed wind field with a predetermined spatiotemporal distribution of EBPs based on real data.Section 4 presents momentum flux estimates from experiments conducted between July 1997 and June 1998 using the Buckland Park MF radar.Discussion and conclusions follow.

HDI implementation on the Buckland Park MF radar
The Buckland Park MF radar, located about 36 km NNW of Adelaide, South Australia, operates at a frequency of 1.98 MHz.It consists of 89 crossed half-wave dipole antennas, each aligned ∼ 4 • west of north, arranged on a rectangular grid with a circular outline of diameter ∼ 1 km.The basic antenna spacing is 3λ 5 (∼ 91.4 m) (where λ is the radar wavelength).A detailed description of the hardware, configurable experiments and analyses can be found in Reid et al. (1995) and Holdsworth and Reid (2004).
The HDI experiments presented here were conducted across four multi-day campaigns during July, September, and October 1997 and June 1998 (see Fig. 4 for an indication of data availability during these periods).The radar utilized 20 channels for transmission and 16 channels for reception.For both transmit and receive modes, each channel was connected to three antennas on the array.The channels were phased to sequentially form transmit and receive beams in five directions: one vertical, and four at small off-zenith angles in the cardinal directions.As shown in Fig. 1 (the antenna configuration for June 1998), 60 antennas on the northsouth-oriented array were used for transmission, and a total of 48 antennas on either orientation for reception.The antenna configuration was varied slightly during and between experiments, as necessitated by antenna outages.
Modelled far-field transmit and receive polar diagrams for the Fig. 1 configuration are shown in Fig. 2. In both cases, the signal phase has been progressed by about 48 • per column of the array (increasing to the right), so as to steer the beam approximately 13 • west of the zenith.As in the real experiment, the three antennas constituting each group have been driven with an equivalent phase, which has been evaluated for the centroid of each group.The main lobe of the transmit beam is approximately circular with a half-power half-width (HPHW) of about 5.7 • .There appears to be some influence of side/grating lobes in the opposite sector of the sky, though given their power relative to the main lobe and their large zenith angles, it is not expected that they will result in ambiguous radial velocities.The receive polar diagram (which is used for beam synthesis in HDI) is slightly wider, with a half-width of around 10 • .There are side/grating lobe influences in the opposite sector of the sky, though again these are not expected to be of major concern.
The sampled range gates encompassed the height range 50-102 km (daytime) and 74-102 km (overnight), in 2 km bins.The transmit pulse's half-power half-width (HPHW) was 2 km.The pulse repetition frequency was fixed at 100 Hz for daytime observations, and 20 Hz for night time, with 20 and 4 software coherent integrations applied, respectively.The time series recorded for each beam contained 560 points (i.e. a record length of 112 s), and the beam direction was changed every 2 min.The beam sequence was [vertical, north, east, south, west], with an off-zenith angle of 12 • used for 1997 experiments, and 13 • for 1998.
HDI was applied to analyse each 112 s raw data record.Briefly, this involved the synthesis of beams in software (using the post-statistics steering method; Kudeki and Woodman, 1990) across an 11 × 11 grid of positions centred on the nominal transmitted beam position, bounded by the e −1 width of the receive beam in the cardinal directions.The EBP was set to the position of the peak of a 2-D Gaussian fit- For ease of viewing, antenna elements have not been drawn to scale.Each bold triangle indicates a group of antennas that were connected to a given transmit or receive channel: "TR" denotes both transmit and receive, and "T" transmit only.Insufficient good power estimates with synthesized beam angle to estimate EBP ted to the distribution of signal power across this grid.The beam was then resteered to the EBP, and power, SNR, radial velocity, and spectral width estimates were subsequently estimated using standard Doppler analysis (Woodman and Guillen, 1974).An example of the EBPs determined with HDI is shown in Fig. 3. Substantial fluctuations in the EBP relative to the transmitted beam directions are clear, as is an effect of atmospheric aspect sensitivity in the lower range gates.Further discussion of these points is taken up in Sect. 5.
The scenarios under which the analysis failed on any given data record are summarized in Table 1.Acceptance rates (and hence data availability) based on these error codes for the four conducted experiments are shown in Fig. 4. It is clear that acceptance rates of greater than 75 % were obtained during the day across all range gates between 76 and about 90 km.At night more than 75 % of data were rejected below about 80 km; the accepted data in this region have still been included in the analysis presented here, and so it is acknowledged that the reported momentum flux estimates may  be biased towards the daytime values (as an aside, note that this is a result of diurnal changes in ionospheric reflectivity (and hence received signal SNR), rather than an artefact of the HDI analysis).The October 1997 and June 1998 results are also still analysed as a single experiment, despite the clear radar outages.

Simulation of momentum flux estimation
A simple computer model has been created to obtain a qualitative assessment of the accuracy and precision of the momentum flux estimates from the five-beam Doppler, vertical beam HDI, and meteor radar techniques.The model propagates monochromatic gravity and tidal waves over the field of view of a "radar", which samples radial wind velocities at positions and times corresponding to real records of the re-spective radar techniques.The approach used has parallels to those used in the following previous works:  4. Nicolls et al. (2012), who applied an approach similar to that used in Fritts et al. (2010a) to narrow-beam fixedlook phased array radars.
5. Vincent et al. (2010), who used a Monte Carlo-based simulation of monochromatic gravity waves with random phases, propagation directions and amplitudes to assess the ability of a meteor radar to measure mean winds and momentum fluxes.
6. Murphy (1992), who used a time-varying wind field to simulate the effects of pointing angle variations on momentum flux extraction with an MF Doppler radar.

Simulation description
The model is based around the following workflow (where necessary, the individual steps are described in more detail in subsequent subsections): 1. Specify a wind field analytically.
2. Acquire an ensemble of scattering positions and times (and add Gaussian-distributed uncertainties to a copy of the positions).
3. Evaluate radial projections of the wind velocity at the "correct" positions/times (measurement noise is effectively added to the radial velocities due to uncertainty in the scattering position).
4. Evaluate the momentum flux components by inverting the fluctuating radial velocity components at the noiseinfluenced positions/times.
5. Evaluate the momentum flux components by computing the covariances of the wind field at a fixed position directly above the simulated radar.
6. Loop back to 1, and repeat for a different realization of the wind field/different uncertainties in the scattering positions.
7. Investigate the mean and standard deviation of the differences between the results of 4 and 5 (herein referred to in the text and in Figs.6-9 as "biases").

Wind field specification
The wind field in the model is comprised of a fixed "mean flow" background velocity (with speed v 0 , bearing ϕ 0 , and no height variation or vertical component), and a superposition of linear waves (which can resemble gravity, Rossby, or tidal waves).It is parameterized in space and time as the velocity vector: where v 0 = [v 0 sin ϕ 0 , v 0 cos ϕ 0 , 0] is the fixed background velocity, n is the number of included waves, and ω i are the vectors of component wave amplitudes, wave vectors, and angular frequencies respectively for the ith wave, r = [x, y, z] is the Cartesian position vector, and t is the time since some arbitrary zero.
Prior to a simulation, the background wind vector, the number of waves to include, and the horizontal perturbation amplitude v h , propagation direction ϕ, ground-based phase speed c p , and ground-based period T for each wave are specified.The remaining parameters in (Eq. 1) are computed using known dispersion and polarization relations (see, e.g., Fritts and Alexander, 2003).The radial component of the wind field is then evaluated at a set of EBPs and times corresponding to those recorded in real samples of Buckland Park MF Doppler/VHF meteor radar measurements (see Fig. 5).For simplicity, the wave parameters are evaluated at a fixed height.A "flat-Earth" coordinate system is also assumed, so that r j = z tan θ j sin φ j , z tan θ j cos φ j , z , where j represents the position index, and θ and φ correspond to the zenith and azimuth angles of the scattering locations, respectively.

Spatiotemporal sampling configurations
The basis of all the scattering positions used in the model is shown in Fig. 5.The scattering positions for the fivebeam Doppler (herein 5BD) technique were obtained from the June 1998 experiment discussed in Sect.2, from the range gate centred on 88 km.The temporal order of the points was preserved, in hope to best account for the effects on the beam position of structures in electron density propagating over the radar's field of view.
Only data with zero error code (see Table 1) were selected, and so some of the points in the resulting time series were missing.To avoid further complicating the results of this simulation with the effects of missing data, an attempt has been made to fill in these gaps and hence make the off-vertical beam dataset "continuous" (i.e. to have the four beams present in each 10 min steering cycle).To do this, the nominal azimuth for each of the off-vertical beams was firstly subtracted from the azimuth of the EBP recorded for each beam.These "wrapped positions" were then assigned to the four off-vertical beams in temporal order, and were "unwrapped" by re-adding the nominal azimuth of the beam the position had been assigned to.The vertical beam positions were simply assigned to subsequent 2 min records in temporal order, again removing the effect of temporal gaps.
A configuration utilizing solely vertical beam data from the 5BD experiment (herein V5BD) was also considered in the model.This resembles an analysis that could be applied on systems with no (practical or otherwise) capability for beam steering on transmission.Vertical transmission in that experimental case was only applied for 2 min per 10 min steering cycle, and so a 10 min analysis interval has been used to represent it here.
A slight variation on the V5BD was also included, which was based on data from experiments consisting of a solely vertical transmitted beam (herein VBD).The HDI-derived EBPs in this case were obtained from experiments run between 20 June 1997 and 15 July 1997, again from the range gate centred on 88 km and from a beam with a half-power half-width of 5.7 • .It should be noted that this technique only differs from the V5BD in that in uses a sampling interval of 2 (rather than 10) min, and also obviously employs a different sample of EBPs from a vertical beam.
A fourth simulated technique based on MF radar Doppler returns was intended to emulate that used in older (and the only yet reported in the literature) MF Doppler experiments for momentum flux estimation, which did not incorporate direct estimates of the EBP (instead they were calculated based on aspect sensitivity estimates; see Sect. 1).In the model, this technique uses radial velocities evaluated at the same positions and times as those in the 5BD technique, but the velocities are assumed to be from scatterers located at fixed zenith angles in the appropriate Cardinal directions.The fixed zenith angle used here was 9 • (recall that the transmitted beam was steered to a zenith of 12 • ).Herein, this technique is referred to as "conventional five-beam Doppler" (C5BD).
Finally, a technique based on 55 MHz all-sky meteor radar data was included.For this case, data from the Buckland Park meteor radar recorded during May 2014 at heights between 88 and 90 km were used.The radar obtained a peak count rate of around 60 h −1 over this height interval and period.Off-zenith beams from 5BD, but without EBP information (inferred from aspect sensitivity)

EBP error
Errors of the form -where xj and yj are the errors along the xand y-direction cosines respectively, 1 and 2 are two numbers drawn from a Gaussian distribution with unit variance, and σ is the desired standard deviation of the distribution -were added to the direction cosines of each scattering position (where θ j and φ j are the zenith and azimuth angles of the scattering positions, respectively).A σ value of 1 • was used to represent the meteor technique, and 1.5 • for techniques based on the Buckland Park MF radar (justification of these values is provided in Spargo (2016), p. 62-63 -though it is highly likely that the latter is an overestimate of the true value).The wind field was evaluated at the points without the error, and in the inversions described in the next subsection, they were interpreted to be at locations corresponding to the positions with the added error.

Estimation of the Reynolds stress components
The spatial distribution of radial velocities at this point will contain contributions from the gravity and tidal waves specified in Eq. (1).Contributions from the latter are not desired in the Reynolds stress component estimates here.To remove them, a method similar to that devised in Andrioli et al. (2013a) is applied; it firstly involves estimating the "background mean" wind field by partitioning the radial velocity (and corresponding EBP) data into non-oversampled windows of width one hour.Wind velocities are estimated using a standard least-squares formulation (e.g.Vandepeer and Reid, 1995).The tidal wave contributions are removed from each radial velocity by subtracting from them a timedependent radial projection of a least-squares fit y to wind time series, of the form where T is an n-element array of periods (in this work, [1/3, 1/2, 1, 2] days) and is an n-element array of phases (in this case giving the time at which the ith component maximizes).
It should be acknowledged that Fritts et al. (2010aFritts et al. ( , b, 2012a, b) , b) used an "S-transform" Gaussian wavelet fit to estimate diurnal and semidiurnal tidal components from meteor radar time series.While this approach does allow for amplitude and phase modulation of the tidal components in the time series, the Andrioli et al. (2013a) approach is applied here on the basis of the ease of its implementation.
For the 5BD, V5BD, VBD, and meteor techniques, the wind field variance and covariance components are then calculated from the residuals by the inversion technique described in Thorsen et al. (1997) (their Eq. 15).Note that the selection of the length of data over which this inversion is to be applied is clearly a compromise of increasing time resolution (shorter windows) and correctly sampling longer-period and/or larger-scale fluctuations (favouring longer windows).The effect of the data length on the bias in the returned momentum flux component estimates is discussed in the next sections.
To emulate the C5BD technique, momentum flux components are calculated using the Vincent and Reid (1983) estimator.

Test cases
Two different wind field "configurations" have been used in the model; they are summarized in Table 2.The first case considers a single gravity wave of horizontal amplitude 20 ms −1 propagating to the north-east.In this case, the horizontal and vertical perturbation components are out of phase (as m ∝ − √ k 2 + l 2 , and k, l > 0) -i.e.u w , v w < 0. The second case considers an ensemble of 37 smalleramplitude waves propagating in uniformly distributed directions in the eastern sector; here only u 2 , v 2 , w 2 , and u w will take on non-zero values, with u 2 v 2 ≈ 1 and w 2 ≈ 15.The waves will have a net propagation direction due east, so u w < 0. The two configurations are intended to emulate the limiting cases observed in real mesospheric wind fields: the first obviously a case in which a single welldefined monochromatic wave dominates the spectrum, and the second in which a spectrum of equal-amplitude waves from an isotropic source propagate with component directions opposite to that of the background wind.
With exception to the final case discussed (Fig. 9), the simulations were performed over a sequence of gravity wave periods so as to test the sensitivity of the techniques to waves of differing scales.A total of 200 realizations were performed at each period in order to obtain a distribution of the measured momentum flux component biases.At the start of each realization, the initial phases of all waves considered were assigned a random value in the interval [0, 2π ).In the case of the second configuration, the periods of each of the 37 waves in a given realization were varied by obtaining periods from the equation T i = R 2 × T i , with R 2 being selected from a uniform distribution with bounds ( 3 4 , 5 4 ).This was essentially done to reduce the correlation distance of the wave field.The phase speeds of the waves in the second configuration were determined as indicated in Table 2.
In each realization, diurnal and semidiurnal tides with fixed amplitudes of 15 and 20 ms −1 respectively were also superposed onto the wind field.Their horizontal wavelengths and phase speeds were adjusted to resemble those of real atmospheric tides (i.e. with horizontal wavelengths equal to the ratio of the Earth's circumference and the tidal mode num-ber, and phase speeds such that a full cycle of a given tidal component would be completed in the ratio of 24 h and the tidal mode number).The spatial variability of the tides was included so that the tide-induced bias of the momentum flux estimates (particularly those derived from the meteor technique) could be inferred (it is assumed that the subtraction of the fit in Eq. 3) will remove most of the temporal variability).

Single gravity wave
The biases for the 5BD and meteor technique momentum fluxes for a wind field with a single gravity wave (Case 1), plotted against the horizontal wavelength of the wave (for a fixed phase speed of 50 ms −1 ), are shown in Fig. 6.Each panel shows a mean value of the "true" components across all the periods examined.Each result is the average of all 2 h blocks in a 48 h time series, over 200 wave field realizations and scattering position errors.The error bars shown depict the standard deviation in the bias of the estimates across all samples (herein, the "accuracy" of the results will be taken to refer to the size of the mean bias, and the "precision" to the size of the standard deviation of the bias).
On the whole, it is clear for such a wave field that both techniques will statistically measure the sign of momentum flux components correctly, and to a similar level of accuracy.The 5BD generally obtains better precision, although it has a tendency to underestimate these components, especially at low and high horizontal wavelengths.This occurs as a result of a failure of the technique to sample full wave cycles.It indicates the requirement for continuous sampling windows (or "integration times") much longer than the maximum gravity wave period under investigation, if unbiased estimates of the covariances at those periods are sought.The obvious downside to this approach is the required assumption for stationarity of the wave field for the duration of the sampling window (this is more likely to be satisfied for a shorter window)., 2, 3, 4, 6, 8, 12, 24, and 48 h, with colours as shown in the key), for a spectrum of waves propagating in the eastern sector (Case 2).Average biases are shown in the upper panels, and the standard deviations of the biases in the lower panels.
The 5BD technique also appears to obtain very low accuracy and precision when the horizontal wavelength considered is such that wave period does not sufficiently exceed the time taken (10 min) for a full five-beam cycle (note that T = λ/c p (where λ is the horizontal wavelength) and so the wave period matches the experiment sampling time for horizontal wavelengths of around 30 km).The fact that the true flux is substantially underestimated at horizontal wavelengths shorter than this indicates an aliased sampling of the wave field at these wavelengths.

Effects of different integration times
Figures 7 and 8 explore the effect of varying the integration time, for Case 2 (which corresponds to an ensemble of waves propagating in the eastern sector).In both of these figures, the upper panels show the mean biases in u w and v w as a function of average horizontal wavelength of the waves in the ensemble (200 realizations of horizontal wavelength and initial wave phases per average horizontal wavelength), and the lower panels show the standard deviation of the bias in the samples used to compute those means.
Results for the 5BD technique are shown in Fig. 7.It is clear that an underestimation of the true value of the non-zero u w occurs at large horizontal wavelengths for shorter integration times.Increasing the integration time generally reduces the bias, though the level of improvement diminishes rapidly once windows much longer than the maximum grav- ity wave period considered are used (note that the average phase speed used in this case is 50 ms −1 , and so the maximum average wave period shown in Fig. 7 is about 200 min).
It is also clear that there is little effect of integration time on the precision of the results for the 5BD technique.
In contrast, results for the meteor technique in Fig. 8 have accuracies exhibiting less dependence on integration time (with biases only readily apparent for the 1 h window), but whose standard deviations are highly dependent on the integration time.In fact, the best precision is obtained for the shortest window used.It is also worth noting that the meteor technique's precision is worse than the 5BD's, and also worsens with average horizontal wavelength of the waves (whereas the 5BD's improves).At present, we do not have an explanation for either of these features.
Clearly, the selected integration time in the analysis is a compromise between the desired accuracy and precision of the techniques.On the basis of our desire to obtain the most accurate possible results for the 5BD and other similar techniques, 48 h windows are adopted for the analysis presented herein.

Spectrum of gravity waves
Figure 9 compares the momentum flux estimate biases of the meteor technique and the four Doppler techniques discussed in Sect.3.1.2.In contrast to the previous section, a "wide" spectrum of gravity waves is used in this case, with the periods of the waves selected from a uniform distribution with bounds (6180) min (and with the propagation direction of the waves spanning the entire eastern sector as in Case 2).The results from a total of 10 000 realizations (each with differ-ent scattering position errors, wave periods, and initial wave phases) are shown.The first two entries in the upper-right corner of each panel indicate the mean and standard deviation of the corresponding distribution.A third statistic, given by 1.4826 × MAD, where MAD is the median absolute deviation of the distribution, is also given; it corresponds approximately to the standard deviation of the distribution subject to the outlier robustness of the MAD, and if the distribution itself were Gaussian (Rousseeuw and Croux, 1993).This has been shown along with the "true" standard deviation, since a few large outliers were present in some of the bias samples -especially those of the Doppler techniques.Standard deviations greatly exceed MAD values for these cases.
The results imply that it is possible to measure momentum flux terms of the correct sign with the V5BD and VBD techniques (for a wave field with a realistic non-zero momentum flux), albeit with less accuracy than both the 5BD and meteor techniques and less precision than the 5BD technique.Like the 5BD technique, both of these techniques have also shown a tendency to underestimate the non-zero u w term in the model.
The C5BD results show poorer precision than the VBD (but greater than the V5BD) and are also substantially biased.The technique overestimates the non-zero u w term and also clearly estimates a non-zero v w term, for which the corresponding true value in the model is very close to zero.
These modelling results have shown, at least qualitatively, the pleasing result that the momentum flux estimation and tide removal procedures employed work to a satisfactory level on all variations of the Doppler technique tested.We particularly stress the finding that EBPs from the V5BD experiment exhibit a spatial variability sufficient to estimate these terms reliably.We also note again that this result incorporates a modelled EBP uncertainty that is likely much larger than that in reality but that is ultimately extremely difficult to quantify (Klövekorn, 1992).

Analysis procedure
The analysis performed here on HDI radial velocities follows the methodology applied to Doppler data in the previous section, making use of the 5BD, V5BD and C5BD techniques to estimate momentum fluxes.
The C5BD technique required that pattern scale data (derived from the full correlation analysis (FCA) Briggs, 1985) exist in the routine analysis run concurrently with the Doppler experiments.Each of the analyses also required the pre-determination of tidal components in the measured velocities, which were projected onto the radial velocities (if they were Cartesian components, as in the first and third techniques) and subtracted from them prior to evaluating any stress terms.An outline of the procedure used to determine the mean winds and momentum flux components is given below.
1. Partition the 2 min resolution data from a given campaign, as well as the concurrent FCA, into nonoversampled 2 h blocks (steps 2 to 7 pertain to each separate 2 h block).
2. Using the pattern scale and axial ratio information derived from the FCA, calculate circularly averaged values of the aspect sensitivity parameter, θ s (see, e.g., Lesicar and Hocking, 1992).
3. Calculate the altitude of each Doppler measurement (from the recorded line-of-sight ranges) using two different techniques: a. Simply assume that the zenith angle of the return is equal to the HDI-derived EBP (θ e ).
b. i. Interpolate a θ s value from the averaged SA-FCA data at an altitude corresponding to the range of the Doppler measurement (assuming θ a = θ e , where θ a is the transmitted beam direction).The interpolant should be determined from a block of data centred on the time of interest, with a block width of the order of a few days.In this study, a block length of 14 days has been used, given the high variance observed in θ s on short timescales.ii.Use the interpolant to calculate an FCA-derived θ e .iii.Use the acquired θ e to calculate the "true" altitude of the Doppler measurement from its original range.
4. Using the two sets of altitudes calculated in Step 3, partition each block into 2 km width bins, with the lowest bin starting at 70 km and the highest at 96 km.For brevity, call the set of bins pertaining to those measurements with altitudes derived from the HDI-based θ e "A", and those from the second method "B".Vandepeer and Reid (1995) on A, estimate the mean horizontal and vertical winds.
than 10 measurements exist across all available beams, consider the wind estimate for this block as "missing".
If calculating winds based on a single beam (e.g. the vertical beam), only compute this if at least three measurements exist.Refer to these velocity estimates as v A .
6. Evaluate mean radial velocities for each nominal beam direction in B. Remove points more than 3 standard deviations from the means, and recalculate the means.Refer to them as V rad Bi (for the ith beam direction).If there are fewer than two points in a given beam, consider this measurement as "missing" and do not perform any further analysis on it.
7. Estimate mean horizontal and vertical winds from B, using data from pairs of off-zenith beams of opposite azimuths.Adjust θ a only for the local value of θ s (i.e.do not apply any correction to the apparent azimuthal angle).Refer to the wind estimates as v B .
8. Re-partition the 2 min resolution Doppler and FCA data into 48 h blocks, with the centre of each block displaced by 6 h from the adjacent one (the remaining steps pertain to each separate 48 h block).
9. Perform a least-squares fit for (Cartesian) tidal/planetary wave components in v A .The fit should be performed over a window encompassing data in the vicinity of the block currently being analysed.In this study, a window width of 4 days (centred on the current block) was used.
10. Subtract a radial projection of the fitted components from the individual radial velocity records in A.
11. Estimate the variance and covariance components from the residuals using the Thorsen et al. (1997) inversion.Again, scale the system of equations with the radial velocity variances.
12. Repeat steps 9 and 10 for the mean radial velocity time series V rad Bi .
13. Simultaneously solve for u w and v w using the Vincent and Reid (1983) estimator.

Momentum fluxes
Unweighted average profiles of the momentum flux components for the four campaigns and the three Doppler techniques are shown in Figs. 10 and 11.The results shown have been evaluated at 2 km resolution, from 76 to 94 km.The uncertainties shown correspond to the standard error in the mean at each height evaluated over each campaign.
In general, the 5BD and C5BD results show the best level of agreement, with especially good agreement at heights where acceptance rates are high.A noteworthy result is the very similar vertical structure from the three techniques' measurements of u w during the June 1998 campaign around 80-90 km.While the V5BD results do show large departures from the those of the other two techniques, it is encouraging to see some level of qualitative agreement.We again stress the point that no transmission beam steering has been used to acquire the V5BD results.
As an aside, the better agreement between the 5BD and C5BD techniques (relative to those from the V5BD) also lends support to the simulation results presented in Fig. 9.However, there does not appear to be substantial qualitative evidence that the C5BD technique systematically overestimates the covariance components, as the results in Fig. 9 predict.

Body forces and Coriolis torques
In an attempt to verify the validity of these experimental results, we (following the approaches of, e.g., Placke et al., 2015 andReid andVincent, 1987) have computed the body forces arising from the vertical divergence of the densityweighted wind field covariance, and have compared them to the Coriolis torque due to the perpendicular mean wind.These quantities should be equal when zonally averaged.Mathematically, the relation is expressed as where ρ is the atmospheric density and f is the Coriolis parameter.Preliminary calculations of these quantities for the four campaigns are shown in Figs. 12 and 13.The densities used were derived from the NRLMSISE-00 model for Buckland Park's location, evaluated at a time resolution of 1 day for an entire year and height resolution of 2 km (with a corresponding density extracted for each height and time spanned by the radar data).Following Reid and Vincent (1987), a function has been fitted to the density-weighted covariance profiles so as to reduce the effects of measurement noise on their derivatives.A quartic polynomial was used here, and was found to adequately replicate most of the major features in the profiles.It also led to smaller least-squares residuals than lower-order polynomials, and lacked the spurious "edge effects" associated with higher orders.The fit was weighted by the inverse square standard error of the individual covariance estimates.
The body force uncertainty σ bf (expressed in the error bars on the profiles of Figs. 12 and 13) was evaluated via the equation where R is the covariance matrix of the least-squares fit parameters (e.g.Tellinghuisen, 2001;Markwardt, 2009), , where F is the analytical form of the fitted function and α i are the fit parameters, and T is the transpose operator.It was assumed that the NRLMSISE-00-derived densities had zero uncertainty.
The vertical structures of the mean inferred body forces and Coriolis torques show few similarities, in both the zonal and meridional planes.Some of these discrepancies may be explained by noting that the relation between the two quantities is only valid for a zonal average.Additionally, the results presented here are centred on the winter months; during this time, planetary waves can propagate into the MLT/I, and may both contribute to the body force and change the "local" mean wind in such a way as to filter gravity waves from the wave spectrum (Andrews et al., 1987).Only the body force contributions from gravity waves have been considered here.
Nevertheless, the inferred body forces are generally large enough to balance the Coriolis torque due to the orthogonal wind.Their senses are also consistent with what is expected: for example, in the July 1997 and June 1998 cases, the body forces are predominantly westward in the heights of the highest acceptance rates.This is consistent with a deceleration of eastward MLT/I winds during winter.
Other studies of the same intercomparisons have had mixed conclusions, which is not surprising given that good local agreement between the quantities is not necessarily expected.Reid and Vincent (1987) compare Coriolis torques and inferred body forces using measurements derived from the Buckland Park MF radar, essentially encompassing all seasons and, as in this study, having measurements taken as part of dedicated campaigns lasting several days.They noted that the zonal body force was usually of the correct sense to balance the Coriolis torque due to the meridional wind, though the agreement varied from excellent (e.g.May 1982, their Fig.12h) to poor (e.g.July 1982, their Fig.12j).Hall et al. (1992) used 17 days of measurements from the Saskatoon MF radar (summer, 1989) to perform a similar comparison, and found good agreement between the zonal body force and meridional torque, but poor agreement in the orthogonal plane.Frame et al. (2000) also considered this intercomparison, using 1 month of data from the Buckland Park and Christchurch (New Zealand) MF radars (May 1992).They obtained good agreement between the gravity-wave-driven body forces and Coriolis torques in the zonal and meridional directions at Christchurch but inconsistent results at Buckland Park.In a more recent study employing DBS on the Saura MF radar, Placke et al. (2015) note good agreement between body forces and Coriolis torques during summer in the MLT/I, but not during winter (as done here, attributing the winter result to planetary wave contributions to the momentum flux).

Discussion and conclusions
This study has suggested, through the use of both synthesized and real observations, that vertical beam MF radars can measure momentum fluxes in the MLT/I to an acceptable degree of accuracy and precision.This implies that the EBP distribution brought about by refractive index irregularities propagating through a fixed vertical beam volume should contain a sufficient spread of radial velocity-pointing direction pairs to solve for the wind field covariances.In particular, we believe this sheds new light on the Thorsen et al. (1997) study, which attempted to make the same assessment through the use of the small-aperture interferometric Urbana MF radar (which transmitted a solely vertical beam).We speculate that the main reason this study has not been followed up is because of the lack of confidence in the community around making momentum flux estimations without steering narrow beams about the zenith, as done in the original studies beginning with Vincent and Reid (1983) (the C5BD in the present paper).We have compared these two approaches in this paper, and have largely found their results to be consistent.It is also clear that a substantial EBP spread occurs in the 5BD results for the same reason as in the VBD; however, as the 5BD-C5BD intercomparison suggests, this is of little consequence, even if the received beam positions are not measured interferometrically.
We thus conclude that vertical beam interferometric MF radars are viable candidates to use for testing momentum flux estimates from interferometric meteor radars, over which there are well-known concerns regarding accuracy and precision, as discussed in Sect. 1.However, we stress that there is more work to be done concerning the prediction of the wave field conditions under which 5BD and V5BD-like techniques will show large discrepancies.
This study did consider basic testing of the meteor technique as well in a simulated setting, using a synthetic GW field consisting of a superposition of monochromatic waves, in much the same way as in previous studies such as Vincent et al. (2010), Fritts et al. (2012a), andAndrioli et al. (2013a).The drawback of the present and these older studies is the lack of realism of the simulated wave field (i.e. its finite spatial correlation, and transient features) and subsequent simulation of the response to this of the typically large EBP distributions of the meteor technique.Placke et al. (2011a) went a step further in this context by using a wind field output from a mechanistic model, but ultimately did not consider a realistic spatiotemporal distribution of meteors in sampling that wind field.In performing the simulations shown in Fig. 8 in the present study, we noticed that the bias standard deviation in the meteor technique's covariance estimates was highly sensitive to the spectral width of the wave field used (i.e. the frequency spanned by the superposed waves, for a given average horizontal wavelength).In the case of using a wave field in which all the superposed waves had equal wavelengths, the bias standard deviation in both the covariance estimates increased by a factor of around 2, relative to that shown in Fig. 8. Clearly, the technique's precision is highly dependent on the correlation length of the wind field (as the assumption in item 1 of Sect. 1 implies), and so a simulation which incorporates this to a realistic extent, along with the spatiotemporal sampling characteristics of a meteor radar, is needed to more fully understand this technique's limitations.
A potential problem with all Doppler beam-steering measurements of momentum fluxes in the presence of nonuniform volume scatter (and hence EBP distributions that may deviate away from an idealized complementary beam arrangement) concerns the extent to which accurate estimates of the vertical wind velocity perturbation can be obtained.A good contemporary review of the well-known biases inherent in volume-scatter-derived wind measurements in the presence of correlations between refractive index fluctuations and the underlying dynamics is provided by Fritts et al. (2012c).This study used a numerical algorithm to compute off-zenith backscatter from simulated Kelvin-Helmholtz instabilities in the mesospheric region (Franke et al., 2011), and revealed the biases in the obtained Doppler spectra as a function of the stage of turbulence development.Simulation-based ap-proaches like this clearly underpin future investigations of the validity of the partial reflection Doppler techniques employed in the present paper.

Figure 1 .
Figure1.Antenna configuration used in the June 1998 HDI experiments.Each thin vertical line denotes an approximately north-southoriented antenna.For ease of viewing, antenna elements have not been drawn to scale.Each bold triangle indicates a group of antennas that were connected to a given transmit or receive channel: "TR" denotes both transmit and receive, and "T" transmit only.

Figure 2 .
Figure 2. Modelled transmit (left) and receive (right) polar diagrams for the antenna configurations in Fig. 1, and used for the HDI experiments.Phasing has been applied to steer the beam away from the zenith; see text for details.The models were produced in EZNEC v. 5.0.63.

Figure 3 .
Figure 3. Normalized 2-D histograms of the EBPs estimated by HDI at selected range gates for the June 1998 experiment.The thin black circles approximately denote the half-power half-width contours of the transmitted beams (5.7 • in all cases shown).
1.Fritts et al. (2010a), who evaluated the abilities of a new meteor radar on Tierra del Fuego (53.8 • S, 67.8 • W) to measure gravity wave momentum fluxes.The wind fields simulated included mean winds, diurnal and semidiurnal tides, and propagating gravity waves with variable phase angles and time-varying amplitudes, and the scattering locations used were based on observed meteor distributions.2.Fritts et al. (2012a), who employed the same tests asFritts et al. (2010a) on a new meteor radar on King George Island (62.1 • S, 58.7 • W), and three more conventional meteor radars.

Figure 4 .
Figure 4. Acceptance rate profiles for the four HDI experiments.

Figure 5 .
Figure 5.A 48 h sample of scattering position data used in the model (black: 5BD; red: VBD; blue: meteor.)See text for details.

Figure 6 .
Figure 6.A comparison of the biases in the model covariance terms extracted from the 5BD (black) and VHF meteor radar (red) techniques, averaged over a 48 h period, with individual bin widths (or integration times) of length 2 h, for the single wave case.The error bars show the standard deviation in the bias determined over 200 realizations of the initial gravity wave phase/scattering position errors.

Figure 7 .
Figure 7.A comparison of the biases in the model momentum flux terms extracted from the 5BD technique for different integration times (1, 2, 3, 4, 6, 8, 12, 24, and 48 h, with colours as shown in the key), for a spectrum of waves propagating in the eastern sector (Case 2).Average biases are shown in the upper panels, and the standard deviations of the biases in the lower panels.

Figure 8 .
Figure 8.As per Fig. 7, but for the all-sky meteor technique.

Figure 10 .
Figure 10.Vertical profiles of the momentum flux components obtained from the July 1997 (upper panels) and September 1997 (lower panels) campaigns.The error bars shown correspond to the standard error in the mean of the samples at each height.

Figure 12 .
Figure 12.Vertical profiles of the zonal and meridional body forces (determined by three independent techniques) and corresponding accelerations due to Coriolis torques, obtained from the July 1997 (upper panels) and September 1997 (lower panels) campaigns.

Table 1 .
Error codes for the HDI analysis.

Table 3 .
Summary of different Doppler techniques referred to in this paper.See Sect.3.1.2for more details on each technique.