Momentum flux determination using the multi-beam Poker Flat Incoherent Scatter Radar

In this paper, we develop an estimator for the vertical flux of horizontal momentum with arbitrary beam pointing, applicable to the case of arbitrary but fixed beam pointing with systems such as the Poker Flat Incoherent Scatter Radar (PFISR). This method uses information from all available beams to resolve the variances of the wind field in addition to the vertical flux of both meridional and zonal momentum, targeted for high-frequency wave motions. The estimator utilises the full covariance of the distributed measurements, which provides a significant reduction in errors over the direct extension of previously developed techniques and allows for the calculation of an error covariance matrix of the estimated quantities. We find that for the PFISR experiment, we can construct an unbiased and robust estimator of the momentum flux if sufficient and proper beam orientations are chosen, which can in the future be optimized for the expected frequency distribution of momentumcontaining scales. However, there is a potential trade-off between biases and standard errors introduced with the new approach, which must be taken into account when assessing the momentum fluxes. We apply the estimator to PFISR measurements on 23 April 2008 and 21 December 2007, from 60–85 km altitude, and show expected results as compared to mean winds and in relation to the measured vertical velocity variances.


Introduction
Geophysical vertical flux of horizontal momentum measurements were pioneered by Lhermitte (1983) and Vincent and Reid (1983), the former for studying the spectral variance in underwater tidal channels and the latter for assessing the vertical structure of wave motions in the atmosphere. In particular, Vincent and Reid (1983) developed a technique for ground-based radars to estimate radial velocity variances and corresponding momentum fluxes in the mesosphere and lower thermosphere (MLT) using dual co-planar narrow radar beams. This technique has since been widely used by HF and VHF radars (e.g., Reid and Vincent, 1987;Fritts and Vincent, 1987;Fukao et al., 1988;Reid et al., 1988;Fritts and Yuan, 1989;Fritts et al., 1990;Sato, 1990Sato, , 1993Sato, , 1994Tsuda et al., 1990;Fritts, 1990, 1991;Hitchman et al., 1992;Nakamura et al., 1993;Murayama et al., 1994;Vincent, 1993, 1998). Similar approaches have been extended to MF broad beam interferometric radars (Thorsen et al., 1997), incoherent scatter radars (ISRs) operating at VHF and UHF Fritts et al., 2006;Janches et al., 2006), and Doppler lidar systems (Acott, 2009). More recently, Hocking (2005) has generalized the approach for application to meteor radars, a technique that has been employed in several studies of gravity waves (GWs) (e.g., Antonita et al., 2008;Clemesha et al., 2009;Placke et al., 2011b,a), validated by Fritts et al. (2010) in applications of the SAAMER meteor radar in studies of mean winds, tides and GWs, and used by Vincent et al. (2010) in an assessment of GW momentum flux measurements by meteor radars. In this paper, we outline a technique for multi-beam radars to estimate the vertical flux of  Hocking (2005), this technique is more applicable to the case of fixed-look directions determined a priori, rather than to the case of line-of-sight wind measurements from random meteor trail locations. This technique is particularly applicable to phased array radars such as the Advanced Modular Incoherent Scatter Radar (AMISR) class of radars.

Measurements
ISRs are able to measure the line-of-sight (LOS) or radial ion/neutral velocity within a small volume in the MLT. A component velocity measurement j , made by an arbitrary beam, can be written as: where v = (v e , v n , v z ) is the vector velocity in geographic coordinates and a j = (a j,e , a j,n , a j,z ) = (cos θ j sin φ j , cos θ j cos φ j , sin θ j ) is a simplified straightline geometry vector with φ j and θ j being the azimuth (east of north) and elevation angles, respectively. (Note that throughout this paper, we will use the notation a j,e , a j,n , a j,z to refer to the eastward, northward, and vertical components for a given beam j .) The measurement has an associated error, e j , with expected value e j = 0.
Single instrument ISR experiments typically consist of measurements v j distributed in altitude, space and/or time depending on the system, its steering dexterity, etc. We will consider cases applicable to the pulse-to-pulse beam steering of the Poker Flat Incoherent Scatter Radar (PFISR).

Vector velocities
As described by Heinselman and Nicolls (2008) and , an arbitrary number of component velocity measurements can be used to estimate the velocity vector. This approach can be expressed in matrix form as a 1,e a 1,n a 2,z a 2,e a 2,n a 2,z . . . . . . . . . a j,e a j,n a j,z A weighted linear least squares estimator for v is then, where ||·|| W represents a weighted 2 -norm with weight matrix given by W = C −1 , where C is the covariance matrix of the measurements. In this case, with uncorrelated errors, C is simply a diagonal matrix with diagonal elements given by e 2 j , the estimated variances of the measurements, where the notation · denotes an expected value, i.e., Solution to Eq. (4) is given by the normal equations (e.g., Aster et al., 2005;Tarantola, 2005), This solution is equivalent to the maximum a posteriori (MAP) estimator. If A T C −1 A is ill-conditioned, or if the application of a physical constraint is appropriate, a regularization/constraint term can be added. For example, this approach has been applied to this problem by Heinselman and Nicolls (2008) and Butler et al. (2010). The error covariance matrix of the solution without a constraint term is given by

Velocity variances and momentum fluxes
An estimator for the radial velocity perturbation iŝ where we assume that after removal of the mean value, v j = 0. The variance of this quantity may be written as where measurement errors have been assumed to be uncorrelated.
The quantity v i v j may be related to component fluxes: v i v j = v 2 n a i,n a j,n + v 2 e a i,e a j,e + v 2 z a i,z a j,z + v n v e a i,e a j,n + a j,e a i,n + v n v z a i,n a j,z + a j,n a i,z + v e v z a i,e a j,z + a j,e a i,z .
For the case of i = j (variance measurement), 2 v n v e a j e a j,n + 2 v n v z a j,n a j,z + 2 v e v z a j,e a j,z .
The vertical fluxes of horizontal momentum correspond to the two last terms, v n v z and v e v z , which we would like to estimate. Several techniques exist in estimating these quantities, including choosing look directions such that unwanted terms in Eq. (13) cancel (such as the co-planar beam method of Vincent and Reid (1983)) or by formulating a matrix of scaled variances, applicable to arbitrary pointing/meteor trail drift measurements (Hocking, 2005). Note that, throughout this section, it is implicitly assumed that the time-averaged quantities are constant over the probing volume -an assumption of horizontal homogeneity. This is assumed to be true of the background wind fields (driven by tidal and other lowfrequency motions) as well as the amplitudes of the wave fields. However, covariance measurements involve the estimation of the expected value of non-homogenous quantities (i.e., v i , the velocity perturbation of a single measurement); thus, one must be cautious when including these terms using distributed measurements. This will be discussed in more detail in the simulation of realistic wave fields (Sect. 3.3). More generally, if the full covariance matrix of measurements is able to be formed, then this can be related to the covariance matrix of the wind field by a standard covariance matrix transformation, Here, for a number of measurements N meas (corresponding, for example, to the number of look directions or beams), is the N meas × N meas covariance matrix of the measurements, the observation matrix A has been defined by Eq.
(2), and the wind covariance matrix is given by A solution to Eq. (14) can be found directly, While useful, this approach is more readily solvable by transforming the two-dimensional and D matrices into onedimensional matrices of observational and estimated quantities. For simplicity and convenience, in our case we simply solve a linear least squares problem for 6 unknowns, v 2 defining D v to be the vector of variances and co-variances. For N meas measurements/beams/look directions, the observation matrix v will be size N meas (N meas + 1)/2. The geometry matrix B is size N meas (N meas + 1)/2 × 6, with elements given by the factors in Eq.
The inclusion of the covariance terms is a unique aspect of this approach and distinguishes it from previous techniques. If those terms were unable to be estimated or were ignored, than the approach would reduce to that of Hocking (2005), except that Hocking (2005) formed vectors of weighted variances and a corresponding geometry matrix, such that the system is exactly determined. That approach seems appropriate for meteor scattering for arbitrary locations and times over a large field-of-view. In the way we have formulated it, the problem may be over-or underdetermined and is an appropriate formulation for fixed experimental beam pointing, wherein the B matrix will remain fixed for any given experiment and where the covariance matrix of the measurements can be computed. We will show via simulation that including the covariance terms reduces the standard errors significantly in the presence of non-zero measurement uncertainties. The estimator using the same approach as in the previous section iŝ A great advantage of this approach is that the a posteriori error covariance matrix (errors on the wind variances and momentum fluxes) can be readily estimated. This matrix is given by Here, vv is the covariance matrix of the measurements v , which can be readily computed. The general covariance of any entry in the covariance matrix is given by, = v iv m v jv n + v iv n v jv m + v iv m e 2 j,n δ j,n + v iv n e 2 j,m δ j,m + v jv m e 2 i,n δ i,n + v jv n e 2 i,m δ i,m + e 2 i,m e 2 j,n δ i,m δ j,n + e 2 i,n e 2 j,m δ i,n δ j,m where we have used the fourth moment theorem for multivariate normally distributed variables to expand the higher order moments, and δ i,j is the Kronecker delta (1 for i = j , 0 otherwise). For example, we see that if i = j = m = n, then and for i = m, j = n, From these expressions, the full covariance matrix of the measurements can be computed, which can be propagated to compute the variance of the estimator. Misestimation of the covariance matrix and violation of the linear least squares assumptions of normality can lead to biased estimates of the momentum fluxes. As a result, in practice, it may be better to use a diagonal prior covariance matrix and only use the estimated covariance matrix for estimates of the errors on the results. As a final note, the bias in Eq. (10) can cause biases in the estimated quantities, particularly if the covariance terms are included as described here. While the biases on the variance estimates can be removed if they are known, we have found that it may be better to estimate them as well. The approach described above can be augmented to account for the biases. In this approach, the forward model is more correctly  Table 1. described to include the bias terms and Eq. (14) is modified such that where C is the diagonal measurement error covariance matrix given by Eq. (5). Equation (18) can then also be modified to add N meas additional unknowns ( e 2 1 , e 2 2 ...) with a corresponding modification in the B matrix.

Simulation description
In this section, we use the principles developed in Sect. 2 to perform a simulation of momentum flux extraction. For these simulations, we use two experiment geometries, as indicated in Table 1 and plotted at 80 km altitude in Fig. 1. Beam map #1 is the experiment geometry used by the Dregion PFISR measurements of . This experiment consists of seven look directions, which are listed and numbered in Table 1, and include a vertical look direction. Beam map #2 is a theoretical beam map with the same number of positions (7), but with the maximum off-zenith angle reduced to 5 • , for reasons that will become apparent in the simulation of realistic wave fields.
Given these beam directions, we can directly compute the a posteriori covariance matrix of the estimator using Eq. (21), assuming a diagonal a priori covariance matrix (in general, this is not a good assumption, as the a priori covariance matrix can be highly non-diagonal given both the properties of the wind field and the fact that measurements are reused to compute the full covariance matrix; this exercise is merely to show the influence of the observation matrix). In Fig. 2, we utilise this calculation to investigate the importance of each beam from Table 1 (beam map #1) in the estimation procedure. In the future, a similar procedure could be used to optimize experiment pointing geometry. The output error matrix of the estimated wind covariances is estimated after removing each beam sequentially from the estimation procedure. The fractional increase in the variance of each estimated quantity is indicated in Fig. 2 (black dots). The red dots show the same calculation, but for the procedure neglecting the measurement covariance terms that were included in the estimator in Sect. 2.
In general, we see that v 2 e and v 2 n are most sensitive to the removal of beams (5,6) and (2,7), respectively. This is not surprising, as these beams have the lowest elevation angles and are the most non-redundant look directions. The quantity v 2 z is very well determined and not that sensitive to any beams. This is somewhat surprising, given the presence of a vertical beam, but the fact that all beams have fairly high elevation angles means that they are all sensitive to vertical motions (and, note that this calculation does not take into account anything about the magnitude of the component variances, expected to be much smaller in the vertical direction). The horizontal wind covariance term, v e v n , is the most poorly determined parameter and roughly equally sensitive to all beams other than the vertical beam. Note that the chosen experiment geometry pointed beams in roughly cardinal directions, which could be improved/optimized in future experiments with consideration about sensitivity to the dominant Using beam map #1, estimate of the fractional increase in the variance of each estimated wind variance/covariance term after removing individual beams (indicated by the x-axis) from the estimation procedure, for (black) procedure including covariance terms and (red) procedure neglecting covariance terms. This process assumes a diagonal measurement error covariance matrix. The absolute value of the variance (for a measurement error covariance matrix given by the identity matrix) using all beams is indicated in the upper left of each panel for the procedure including the covariance terms. The ratio of the absolute values using all beams is indicated by the horizontal red dashed line. This turns out to be an overestimate, due mostly to the covariance of the measurements. scales and frequencies. The momentum flux terms, v e v z and v n v z seem fairly well-determined and largely sensitive to individual beams pointed in the westward and southward directions, respectively, again due to the non-redundancy of these beams.
Comparing the results, including the covariance terms outlined in Sect. 2 (black) and neglecting those terms (red), we see that in general including the covariance terms leads to less sensitivity in removing individual beams. The absolute magnitude of the variances is much smaller for the estimates including the covariances -a factor of ∼10 in root-meansquare error for the term v e v n , and a factor of ∼1.5-3 for all other terms. However, these estimates neglect the role of nondiagonal measurement error covariance matrix, which will be included next when examining the results of simulations.
To investigate the statistical properties and any biases associated with momentum flux extraction using the technique developed, we performed two types of simulations: the first (Sect. 3.2), to investigate the statistical aspects of momentum flux extraction, and the second (Sect. 3.3) to simulate realistic wave fields.

Simulation 1: statistical properties of momentum flux extraction
For the purposes of an initial simulation designed to investigate the statistical properties of momentum flux extraction, we use the beam configuration described in the previous section with an assumed, non-essential background wind environment v. We then perform the following steps for N 1 iterations of the simulation: -Generate a 3 × 3 matrix, L, with elements drawn from a uniform distribution from −5 to 5 m s −1 .
-Use this matrix to generate a symmetric positivesemidefinite covariance matrix L T L which now represents the wind variances ( v 2 e , v 2 n and v 2 z , diagonal elements) and covariances ( v e v n , v e v z and v n v z , off-diagonal elements). This approach will generate a wind covariance matrix with variance terms with a median value of ∼25 m 2 s −2 and ranging up to ∼50-70 m 2 s −2 . Covariance terms will have a median value close to 0 m 2 s −2 , ranging from ∼ −50 to 50 m 2 s −2 .
-For N 2 iterations, generate a time series of wind vectors with length N 3 points (corresponding to the number of samples used to estimate the variances, covariances and wind vectors), drawn from the multi-variate normal distribution with covariance matrix L T L and mean v. Note that the resulting uncertainties on the wind variances and covariances will be highly dependent on N 3 , as shown and discussed later.
-For each look direction, generate observed line-of-sight wind components using Eq. (9) with errors drawn from the zero-mean normal distribution with variance e 2 .
-Utilise Eq. (20) to estimate the momentum flux and wind component variances, as well as estimated measurement errors ê 2 j . -Estimate the means and standard deviations of the variance and covariance terms over the N 2 iterations, which can be used to compute the standard error on the estimate as well as any biases on the results.
We have run simulations as described above using N 1 = N 2 = 1000 and N 3 = 60, with σ e = e 2 equal to 0, 2 and 5 m s −1 . Results are shown in Fig. 3, where each panel shows a histogram of the difference between the true value and the estimated value (averaged over the N 2 iterations). This represents a histogram of the bias of the estimator. Here, blue lines correspond to the σ e = 2 m s −1 case and red lines correspond to the σ e = 5 m s −1 case. Solid lines correspond to the standard momentum flux approach, neglecting the covariance terms, whereas dashed lines correspond to the technique described in this paper. The variances, v 2 e , v 2 n and v 2 z , are reasonably well determined for the small error case, especially in the vertical direction. As already mentioned, this is due to the high elevation angles of the beams and the presence of a vertical beam. However, the variances are clearly biased (in all three directions), which is most clearly seen on examination of the v 2 z histogram where the red curves are offset by 25 m 2 s −2 (and the blue curves by 4 m 2 s −2 ). As discussed in Sect. 2, there is an expected bias in the line-of-sight variance by e 2   that propagates to the wind variances. This bias is removed in the approach described (dashed lines), because an estimate of the measurement error variance was removed from the measurement variances. The bias, however, does not affect the covariance terms. For these terms, we see unbiased results in all cases, relatively poorly determined for the v e v n term and reasonably well determined for the vertical flux results.
The approach using the covariance terms (dashed lines) significantly outperforms the traditional approach, as seen by the standard deviation of the biases -the distribution of the biases is much narrower. This is especially true for the v e v n term, which is due to the additional information provided by the covariances.
The errors in the approach are quantified in more detail in Fig. 4, where we plot histograms of the standard deviations of the estimates in a similar format to Fig. 3. Standard deviations were computed using the N 2 iterations and represent the error in the estimate. Note that we earlier described an approach to estimate errors from the covariance matrix of the measurements, and we have confirmed that the standard deviations plotted in Fig. 4 agree with these esti-mates. Errors increase significantly as σ e increases, but are much lower, for all parameters, for the technique using the covariance terms. For the momentum flux terms, the errors decrease by about 25 % for the σ e = 2 m s −1 case and about 40 % for the σ e = 5 m s −1 .
The conclusion that we draw from these simulations is that with (1) appropriate beam geometry and (2) sufficiently precise measurements, we can make unbiased estimates with reasonable errors of the vertical flux of horizontal momentum. These errors can quickly become untenable if these criteria are relaxed. Errors are significantly reduced if the full covariance matrix of the measurements can be computed, as is the case for PFISR measurements.
Given that the momentum flux is a measure of the correlation between horizontal and vertical wind variations, its determination is statistical in nature: a single measurement, even with no instrument errors, has at least 100 % uncertainties, determined by the geometric mean of the horizontal and vertical wind variances (e.g., v 2 n v 2 z ) and the fractional variance embodied in the momentum flux term (e.g., x v 2 z (where results for v e v z and v n v z have been merged). Blue and green lines with errorbars show mean and 1σ spread for the black and red points, respectively, and gray line corresponds to v 2 x v 2 z /N 3 . In all panels, black dots are from the simulation neglecting the covariance terms, where red dots are from the approach described in the text. Kudeki and Franke, 1998;Thorsen et al., 2000;Vincent et al., 2010). For the case of this simulation, we have assumed that N 3 = 60 statistically independent points exist with which to compute the line-of-sight variances and corresponding wind variances and co-variances. In Fig. 5, this issue is demonstrated. Figure 5 (left) shows the fractional error in estimating the horizontal momentum flux as a function of the correlation between the horizontal and vertical wind components, for the case of no measurement errors. In this case, the approach described here including the covariance terms agrees exactly with the traditional approach (or, that of Hocking, 2005, for arbitrary pointing). When the correlation is near 1, the fractional error approaches √ 2/ √ N 3 = 0.18. As the correlation decreases, the fractional error increases dramatically, demonstrating the importance of this quantity (e.g., see Kudeki and Table 1) (top row) and beam map #2 (bottom row). Black curves show the solutions with only the variance terms included in the estimator, and red terms show results including the full covariance matrix. Simulations are with σ e = 2 m s −1 .
Franke, 1998). The right-hand panel shows the absolute error for this simulation, where it is quite clear for any given value of wind variances, v 2 x v 2 z /N 3 (gray line) is a lower bound on the error (note that v x here is either component of the wind field, v e or v n ). The true error can be larger than this, depending on the magnitude of the correlation between the horizontal and vertical fluctuations. The importance of independent samples should be clear, as the only way to reduce the uncertainties on the momentum flux estimates is to incoherently average in time or space. Adding measurement errors (lower two rows of Fig. 5) clearly makes things worse, as one would expect. However, including the covariance terms significantly reduces the fractional and absolute errors -∼30-40 % in the σ e = 2 m s −1 and ∼60-70 % in the σ e = 5 m s −1 .

Simulation 2: realistic wave field
In this section, we extend the simulations to include a realistic wave field, mainly to investigate the effects of realistic horizontal variations in the perturbed radial velocity fields, similar to the approach taken by others (e.g., Fritts et al., 2010Fritts et al., , 2012Vincent et al., 2010). We simulate a case which is similar to Case 1 of Fritts et al. (2010) and Fritts et al. (2012), in which we have fixed mean, diurnal tide and semi-diurnal tidal amplitudes for the background winds in addition to imposed gravity wave variations.
In our modified Case 1, we have two gravity waves, with fixed amplitudes as given in Table 2, propagating orthogonally to one another with a horizontal wavelength and period drawn from a uniform random distribution from −500 to 500 km, and 6 to 45 min, respectively. These waves will lead to fixed wind variances and momentum fluxes, given in Table 2. For the simulation, we follow the following steps for N 1 iterations of the simulation to generate a realistic time series of measurements: -Draw a wavelength and period for each of the two waves from a uniform random distribution, as described.
-For N 2 iterations, generate observed line-of-sight wind components using the background winds, horizontally varying gravity wave fields and measurement errors drawn from the zero-mean normal distribution with variance σ 2 e .
-Estimate the momentum flux and wind component variances for each of those iterations.
-Estimate the means and standard deviations of the variance and covariance terms over the N 2 iterations, which  can be used to compute the standard error on the estimate as well as any biases on the results.
The time series of the measurements corresponds to a twohour interval with samples every 3 min, similar to the measurements that will be presented later. The simulation is performed at a single altitude of 80 km. It is crucially important to accurately estimate the covariance terms, for which a direct covariance calculation is not suitable because of scale-dependent propagation delays across the field-of-view. To do so, we employ a method wherein we take the Fourier transform of the crosscovariance function, or the cross-power spectral density, between measurement i and measurement j . We high-pass filter the cross-spectrum in the frequency domain to eliminate low frequency motions with periods larger than an hour, which we are not able to address using the types of measurements targeted with this technique. We then integrate the cross-spectral density to obtain the magnitude of the covariance. The sign of the covariance is taken from the zero-lag of the cross-covariance function. Short-wavelength, spatially aliased waves may thus contribute to biased results, hence our desire to simulate a breadth of wave scales.
The results for Case 1 are shown in Figs. 6 and 7. Figure 6 shows the bias in the momentum fluxes and vertical velocity variance (average estimate minus true value) for the solutions including the variance terms in the estimator (black points) and the solutions including the covariance terms as outlined in this paper. Data are binned in terms of horizontal wavelength of the waves, since this was found to be a strong controlling factor on the nature of the bias. The top row of Fig. 6 utilises beam map #1 as indicated in Table 1. At large wavelengths, there is a clear bias in both methodologies, approaching ∼5-10 m 2 s −2 in the momentum fluxes. At large horizontal wavelengths, the vertical velocity variance is unbiased as a result of accounting for the biased lineof-sight velocity variances. At small-horizontal wavelengths (<∼100 km) the method including the covariance terms becomes significantly oppositely biased. This is clearly an effect of the chosen geometry and scale-dependent propagation delays across the field-of-view. The bottom row of Fig. 6 shows the same results for beam map #2, which has small off-zenith angles. Both the magnitude of the bias as well as the wavelength at which biases begins to deviate are significantly reduced. The importance of chosen beam geometry on the accuracy of the technique seems very clear. It should also be emphasized that the nature of any biases likely depends strongly on the wave field present at any given time. Figure 7 shows a similar plot for the standard errors on the estimates and demonstrates a tradeoff between accuracy and precision. While eliminating the covariance terms from the approach results in a less-biased solution, standard errors are significantly reduced by retaining those terms. Similarly, choosing a beam pattern with smaller off-zenith angles increases the accuracy of the estimates, but also results in larger standard errors. In the cases where biases are smaller than expected errors (such as in the case of larger measurement errors), it would make sense to trade increased precision for accuracy. In any case, for a fixed experimental beam geometry, both solutions could be computed and compared to obtain the most quantitative picture possible.

Observations on 23 April 2008
We have applied the technique described in this paper to PFISR data from 23 April 2008, which has been described and analysed in detail by . A large-scale, large-amplitude inertia-gravity wave was observed over Poker Flat during this period. For our purposes here, we have averaged incoherent scatter radar spectra in ∼3 min integrations, corresponding to a sufficiently short average to capture the expected frequency range of velocity variability, and long enough to ensure that samples are approximately statistically independent. Spectra were fit as described in  for spectral parameters including the line-of-sight drift velocity, taken to be equal to the component of the neutral wind along the radar line-of-sight. Lineof-sight velocities were averaged over 4 km in altitude (and analysed every 1 km). The covariance matrix of the measurements were computed over a two-hour period, corresponding to 38 independent variance estimates (which can be equated to N 3 in the simulations in Sect. 3.2). As shown in the simulations and by Kudeki and Franke (1998), we expect the minimum error on the zonal and meridional momentum flux estimates over the two-hour periods to correspond to v 2 e v 2 z /38 and v 2 n v 2 z /38, respectively. In Fig. 8, we show the two-hour averaged winds and momentum fluxes in the zonal (top) and meridional (middle) directions, plotted every hour. The short vertical wavelength wave activity is clearly evident in the zonal and meridional winds. The momentum flux is well-behaved and is roughly anti-correlated with the wave perturbations, increasing in amplitude between 75 and 85 km. In the lower panel, we have plotted the averaged vertical velocities, along with the vertical velocity variances. The two-hour averaged vertical velocities fluctuate up to about 0.25 m s −1 (but more typically quite a bit less than that), and the variances for this dataset are quite small.
In Fig. 9, we show averaged results for the 11-h period from 13:00-24:00 UT on 23 April. The top panels of Fig. 9 show the averaged winds and momentum fluxes for this 11-h period, determined by averaging each 2-h estimate of the wind vector and momentum fluxes. Errorbars on the winds (black lines) correspond to standard deviations over the 11-h period and are meant to indicate the variability of the background winds embodied in Fig. 8. The winds are predominantly westward and northward over this period, increasing with altitude and peaking at an altitude of about 80 km. The momentum fluxes are indicated by red lines, with errorbars corresponding to standard errors on the mean using the 2-h estimates over the 11-h period. The average momentum estimates of the variance in each 2-h average, the dashed red lines show the expected (minimum) root-mean-square error on the momentum flux, again averaged over the 11-h period. The solid lines show the averaged errors determined from the mean-square error of the least-squares fitting process using the covariance matrix of the measurements. The dotted lines show the standard deviation of the measurements over the 13:00-24:00 UT period. The curves indicate expected errors on the 2-hr momentum flux estimates of 0-2 m 2 s −2 . The curves have very similar altitude trends, with the estimated errors from the covariance matrix of the measurements somewhat underestimating the variability of the data and the expected minimum errors.

Observations on 21 December 2007
In Figs. 10 and 11, we show results in the same format for 21 December 2007 for the period from 18:00-24:00 UT when signal-to-noise ratios were sufficient to determine winds and momentum fluxes. In this case, data were averaged for 2 min before fitting for the radar spectra, and the analysis was done with 6-km altitude averaging (and computed every 3 km) to ensure robust estimates. In this case, vertical velocity excursions were significantly larger than in the previous case, with vertical velocity variances approaching ∼10 m 2 s −2 for individual two-hour windows, but more typically quite a bit less than that. The momentum fluxes are correspondingly considerably larger (peaking at ∼25 m 2 s −2 ), with fluctuations that mimic the dynamic background wind variations. However, in this case, it is important to realise that the horizontal and vertical wind variations are sufficiently large that it may not be possible to robustly determine the momentum fluxes for a given two-hour window. As shown in the lower panels of Fig. 11, the expected averaged minimum root-mean-square error on the momentum flux is near 5 m 2 s −2 . The averaged results in Fig. 11 show that the average momentum flux for this 6-h period is between −10-0 m 2 s −2 in the zonal direction and ∼ ±10 m 2 s 2 in the meridional direction (varying with altitude). These average results should be significant given the statistical uncertainties and demonstrate larger momentum fluxes on this day associated with the enhanced vertical velocities and variances.
The predominant flux appears to be southwestward, opposite the direction of the background winds. Comparisons to the variance-only measurements show generally good agreement in the meridional direction, but indeterminate results in the zonal direction.

Conclusions
We have demonstrated a technique to estimate the covariance matrix of the wind field, including wind variances and momentum fluxes, using arbitrary, distributed measurements of line-of-sight winds. This approach is applicable especially to systems with flexible, but fixed pointing, such as phased array radars like the Poker Flat ISR, and is useful especially for estimates of higher-frequency and shorter lived wave inputs rather than long-term averages.
The approach is simply a generalization of previous techniques to include covariance terms. We demonstrated via simulation that the approach produces robust estimates of the wind variances and covariances that outperform the direct extension of the co-planar beam technique in the presence of finite measurement uncertainties because of the inclusion of covariance terms in the forward model. Accurate estimation of the covariance terms is crucial to the approach, and scaledependent propagation delays over the field-of-view of the system must be accounted for when estimating them. Misestimation of these terms and issues associated with experiment geometry, may lead to biased estimates of the momentum fluxes and wind variances. However, while complicating the approach, including the covariance terms provides a significant reduction in standard errors that may outweigh any potential biasing by misestimation. The covariance terms can easily be excluded to provide a less biased, but likely more error-prone, estimate of the momentum fluxes and wind variances.
Future work can optimize the beam positions, spatial averaging and time integration to take advantage of the expected range of momentum-containing scales. In addition, other sources of potential biases, such as beam-pointing accuracy, should be investigated.