A six-beam method to measure turbulence statistics using ground-based wind lidars

. A so-called six-beam method is proposed to measure atmospheric turbulence using a ground-based wind lidar. This method requires measurement of the radial velocity variances at ﬁve equally spaced azimuth angles on the base of a scanning cone and one measurement at the centre of the scanning circle, i.e.using a vertical beam at the same height. The scanning conﬁguration is optimized to minimize the sum of the random errors in the measurement of the second-order moments of the components ( u,v,w ) of the wind ﬁeld. We present this method as an alternative to the so-called velocity azimuth display (VAD) method that is routinely used in commercial wind lidars, and which usually results in significant averaging effects of measured turbulence. In the VAD method, the high frequency radial velocity measurements are used instead of their variances. The measurements are performed using a pulsed lidar (WindScanner), and the derived turbulence statistics (using both methods) such as the u and v variances are compared with those obtained from a reference cup anemometer and a wind vane at 89 m height under different atmospheric stabilities. The measurements show that in comparison to the reference cup anemometer, depending on the atmospheric stability and the wind ﬁeld component, the six-beam method measures between 85 and 101 % of the reference


Introduction
Wind lidars are being used significantly for wind energy applications. They measure mean wind speeds with great accuracy, and are very useful tools in the measurement of wind profiles (Smith et al., 2006;Kindler et al., 2007;Peña et al., 2009;Wagner et al., 2011). New recommended practices are being defined for wind resource assessments (Clifton et al., 2013). However their use in measuring atmospheric turbulence has not yet been established, particularly with the commercial lidars (Sathe et al., 2011b). The main reason is that for a commercial lidar, the measured lidar data is processed using the so-called velocity azimuth display (VAD) method, where the measurements of the radial velocity (also called the line-of-sight velocity) at different azimuth angles are combined to deduce the wind field components. For the mean wind speed estimation, the VAD method produces negligible errors. For turbulence statistics the VAD method produces significant systematic errors (Sathe et al., 2011b;Sathe and Mann, 2012) mainly due to two reasons; one is the filtering of the smaller scales due to the large size of the probe volume within which the radial velocity is measured, and second is the contamination by the two-point correlation between the components of the wind field.
In this article we present a so-called six-beam method that significantly improves the measurement of turbulence in comparison to the VAD method. This method uses the variances of the radial velocities from six different lidar beams, five of which are at equally spaced azimuth angles on the base of a scanning cone and one beam is vertical. These variances are then combined in order to deduce the second-order moments of the wind field. A framework for this method was originally proposed by Lhermitte (1969), which was used by Wilson (1970) and Kropfli (1986) for radar studies, and subsequently by Eberhard et al. (1989) and Mann et al. (2010) for lidar studies of turbulence measurements. In their studies only the covariances were estimated, either by combining several measurements of the radial velocity variances from several lidar beams (Eberhard et al., 1989) at equally spaced azimuth angles and one elevation angle, or by using only two lidar beams (Mann et al., 2010). In the present work, six beams are used, five at an elevation angle of 45 • and one vertical that enable us to also deduce the variances.
The ideas to measure turbulence using remote sensing instruments have evolved, albeit slowly, since the pioneering works on radar meteorology (Lhermitte, 1962;Browning and Wexler, 1968). Based on the VAD scanning, Lhermitte (1969) was the first (to our knowledge) to suggest a technique of deducing turbulence using the measurements of the variance of the radial velocity. Subsequently Wilson (1970) was the first to conduct an experiment using a pulsed Doppler radar and deducing turbulence in the convective boundary layer (0.1-1.3 km). Only turbulence scales larger than the pulse volume but smaller than the scanning circle could be measured since all the data from a single scan was used. Also, no comparison with any reference instrument was carried out, and hence, the reliability of the radar measurements could not be verified. Kropfli (1986) extended the study of Wilson (1970) to also include the turbulence scales larger than the scanning circle by using the data from multiple scans. Although the method was developed for Doppler radar studies, it could also be used for Doppler lidar studies. Eberhard et al. (1989) was the first to perform turbulence studies using a lidar following the methods of Wilson (1970); Kropfli (1986). Gal-Chen et al. (1992) also used the variances of the radial velocities to deduce turbulence, but with a different scanning configuration. In all of the aforementioned studies with a Doppler lidar (or radar), the probe length was quite significant (of the order of 100 m), which perhaps was the reason to restrict these studies to the convective boundary layer. However if the turbulence measurements were desired close to the ground then they would be subjected to a significant amount of probe volume averaging. It was perhaps this reason that the focus on turbulence research with lidars shifted to understanding the probe volume averaging effect and providing potential solutions to compensate for it (Frehlich, 1994(Frehlich, , 1997Frehlich et al., 1994Frehlich et al., , 1998Frehlich et al., , 2006Frehlich and Cornman, 2002;Frehlich and Kelley, 2008;Banakh et al., 1995aBanakh et al., , b, 1996Banakh et al., , 1999Banakh et al., , 2010Banakh and Smalikho, 1997a, b;Banakh and Werner, 2005;Smalikho, 1995;Smalikho et al., 2005;Mann et al., 2010;Branlard et al., 2013). Even with the development of the modern lidar systems, where the probe lengths have shrunk to about 30 m for a pulsed lidar, significant amount of averaging still remains in the turbulence measurements within the surface layer, where the wind turbines operate (Mann et al., , 2010Sjöholm et al., 2009;Sathe et al., 2011b;Sathe and Mann, 2012). A detailed review of the state of the art with respect to turbulence measurements using ground-based wind lidars can be found in .
Unfortunately within the wind energy sector, turbulence measurements are being deduced using the VAD scanning method, which results in a significant amount of filtering of turbulence, and contamination by the two-point correla- tion between the components of the wind field (Sathe et al., 2011b). In this work we attempt to significantly improve the turbulence measurements compared to those obtained by the VAD method, by extending the previously developed ideas of using the radial velocity variances (Lhermitte, 1969;Wilson, 1970;Kropfli, 1986;Eberhard et al., 1989;Mann et al., 2010), but restricting them to using only six beams.
The structure of this article is divided into the following sections. Section 2 gives a detailed explanation of the six-beam technique. The optimum six-beam configuration, which is one of the main contributions of this article is also described in detail. In order to verify our method, turbulence measurements using the pulsed lidar WindScanner were performed and compared with a reference cup anemometer at a height of 89 m. The site description for the measurements is given in Sect. 3, whereas the results are described in Sect. 4. Discussions and conclusions are made in Sects. 5 and 6, respectively.

Theory of six-beam configuration
The instantaneous velocity field is characterized as a vector v = (u, v, w), and turbulence is characterized as the components of the Reynolds stress tensor, where the diagonal terms are the variances of the respective wind field components and the off-diagonal terms are the covariances, denote ensemble average, and denotes fluctuations around the average.  Fig. 1, at a given instant of time if we assume that a lidar measures at a point, and that the lidar beam is inclined at a certain zenith angle φ (in some literature the complement of φ is used, which is called as the elevation angle α = 90 • −φ) from the vertical axis, and makes an azimuth angle θ with respect to the axes in the horizontal plane, then the radial velocity (also called as the line-of-sight velocity) can be mathematically written as, where v r is the radial velocity measured at a point, n = (cos θ sin φ, sin θ sin φ, cos φ) is the unit directional vector for a given φ and θ , and d f is the distance at which the measurement is obtained. In Eq.
(2), we have implicitly assumed that v r is positive for the wind going away from the lidar axis, the coordinate system is right-handed, and u is aligned with the x 1 axis in a horizontal plane, i.e. from west to east. In reality, a lidar never receives backscatter from exactly a point, but from all over the physical space. Fortunately the transverse dimensions of a lidar beam are much smaller than the longitudinal dimensions, and for all practical purposes we can consider that the backscatter is received only along the lidar beam axis. Mathematically the radial velocity can be represented as the convolved signal, where v r is the weighted average radial velocity, ϕ(s) is any weighting function integrating to one that depends on the type of lidar, i.e. a continuous wave (c-w) lidar or a pulsed lidar, and s is the distance along the beam from the measurement point of interest. From simple geometrical considerations the radial velocity variance can be written as a function of the components of R (Lhermitte, 1969;Eberhard et al., 1989), v r 2 = u 2 sin 2 φcos 2 θ + v 2 sin 2 φsin 2 θ + w 2 cos 2 φ (4) + 2 u v sin 2 φ sin θ cos θ + 2 u w sin φ cos φ cos θ where v r 2 is the radial velocity variance. From Eq. (4) we can see that for a given θ and φ, if we have six measurements of v r 2 then there are six unknowns to be determined, which in a matrix form can be written as, where is a vector of the components of R (because R is symmetric, we only need six components), M is a 6 × 6 matrix of the coefficients of that consist of different combinations of θ and φ (see Eq. 4), and S is a vector of measurements of v r 2 at different θ and φ (where the suffices denote measurements from beam 1 to 6). In principle we can then estimate using the relation = M −1 S, where −1 denotes matrix inverse. It is interesting to know beforehand whether the measurements from the six beams on only one zenith angle are adequate, i.e. whether we can have six θs and only one φ.
From fundamental algebra we understand that Eq. (5) will have a finite solution if and only if det M = 0, where det denotes the determinant of a matrix. In other words M should not be a degenerate matrix. From the properties of determinants we know that if any two rows (or columns) of a matrix are identical then its determinant is zero. Also, if the elements of any row (or column) are increased (or decreased) by equal multiples of the corresponding elements of any other row (or column), the value of determinant is unchanged. If we use only one φ at different θ , and add the first two columns of M, we get the first and the third columns of M to be multiples of each other, which according to the property of determinants implies det M = 0. Thus M becomes degenerate if we use only one φ, and thus need v r 2 measurements from more than one φ.
We are then confronted with the challenge of obtaining an optimum combination of θ and φ. Measured S is stochastic, and the random error of will depend on the particular choice of the θ s and φs. We thus choose the objective function such that the sum of the random errors of the components of are minimized. For simplicity, we neglect the probe volume filtering effect in the derivation of the optimum combination, but including that will not change the optimum configuration.
It is to be noted that in Eq. (5), we have neglected the random error in the lidar estimate of the radial velocity. For the particular lidar used in this work the uncorrelated noise in the velocity estimation is exceedingly small and will not have any effect for the analysis presented in this paper. The very low noise level can be seen from plots of spectra in Figs. 2, 5, and 7 of Sathe and Mann (2012).

Formulation of the objective function
Equation (4) is valid for the mean wind direction aligned in the x 1 direction. Following standard meteorological conventions, let us consider the mean wind direction to be at an angle with respect to the north, i.e. x 2 axis as shown in Fig. 2. At first we derive an objective function for the wind aligned with the x 1 axis, and then extend the derivation to the coordinate system aligned with the mean wind direction.

Mean wind aligned with the x 1 axis
If we consider that δ is the random error on , and δS is the random error on S, then Eq. (5) can be written as, We can thus write, If we consider the sum of the error variances δ T δ , where T denotes matrix transpose, then the objective is to minimize the sum of the error variances of the components of . Taking the transpose and multiplying by Eq. (7) we get, The task now is to simplify Eq. (8) such that it can be represented as a function of θ and φ only. If we assume that the random errors in the variances of the radial velocities are independent of each other, and that the error variance for each radial velocity variance is 2 s , we get, where Tr is the trace of a matrix. The objective function is to minimize Eq. (9).

Coordinate system aligned with the mean wind direction
In order to align the coordinate system with the mean wind direction, we need to apply coordinate transformations on any tensors that are defined in the original coordinate system. The vector v rotated in the mean wind direction has to be multiplied by a transformation matrix T given as, In the coordinate system aligned with the mean wind direction, we then get in matrix form, where R r is the Reynolds stress tensor in a coordinate system aligned with the mean wind direction. If we denote r as the vector of the components of R r , then we can write, Using Eqs. (6) and (7), we can write, Following the same procedure as in Sect. 2.1.1, we get Equation (14) states that the error variance is dependent on the mean wind direction. In order to make it independent of the mean wind direction, we assume a uniform distribution of the mean wind direction, and estimate the averaged ratio of the error variance. Thus the directionally averaged ratio is, We can also switch the order between integration and matrix trace, i.e. either we can estimate the trace first and then the integration or vice-versa. Thus, Solving the integral we get, The objective is to minimize Eq. (17), subject to the constraints that θ varies between 0 and 360 • and φ varies between 0 and 45 • . The limit of 45 • for φ is arbitrary, and is based on the considerations of statistical homogeneity in a horizontal plane. Depending on the type of the terrain the range of φ could thus be increased or decreased, i.e. if a terrain is horizontally homogeneous over a very large extent, then φ could be greater than 45 • and vice-versa.

Optimizing the objective function
Equation (17) represents a non-linear optimization problem with 12 unknown variables, i.e. six unknowns are zenith angles φ, and the remaining six are the azimuth angles θ . Owing to the complexity of the optimization problem, an analytical solution of Eq. (17) is not possible. We thus use numerical methods, where either gradient or direct search methods could be used (Rao, 2009). For gradient methods, it is essential that the objective function is differentiable. However we assume that Eq. (17) is a discontinuous function, and hence, we do not use gradient methods. Thus we optimize Eq. (17) using direct search methods only. The main advantage of using direct search methods is that they can be used for discontinuous and non-differentiable functions. The main limitation of such methods is that the found optimum may only be a local optimum.
Different algorithms such as simplex (Nelder and Mead, 1965), simulated annealing (Ingber, 1993) and random search (Rao, 2009) are tested, which result in the optimum angles as given in Table 1, which shows that the optimum configuration consists of five beams equally spaced on the base of a scanning cone and one vertical beam.   of the random errors in turbulence measurements (Lenschow et al., 1994). The site is about 2 km from the west coast of Denmark. The eastern sector is characterized by flat homogeneous terrain, and to the south is a lagoon. The Wind-Scanner is placed at the UTM zone 32 V 447 188 m E and 6 256 189 m N (WGS84 datum) which is about 41 m away from the met mast in the west direction. Since the wind turbines are to the east of the WindScanner and the met mast, only the measurements from the western sector (225-315 • ) are used. Owing to the sudden change in the surface roughness from sea to land in the western sector, we expect the turbulence structure to be influenced by the development of the internal boundary layer, particularly under different atmospheric stabilities. However we do not expect a significant influence on the flow homogeneity in the horizontal direction around the scanning circle, which is one of the key assumptions of the six-beam method.
The duration of the full cycle of the six-beam measurements from the WindScanner was about 15 s. The period of measurement was between 1 and 28 July 2013, where 764 30 min periods were measured. After filtering for data availability within each 30 min period, where only those periods were chosen with 95 % data, the number of 30 min periods reduced to 625. Finally filtering for wind directions to avoid wakes from the wind turbines and the met mast rendered 401 30 min periods. The available 30 min ensembles are further classified into different atmospheric stabilities, characterized by Monin-Obukhov length L MO based on the intervals given in Table 2 (Sathe et al., 2011a).  L MO is estimated using the eddy covariance method (Kaimal and Finnigan, 1994) from the high frequency (20 Hz) measurements of a sonic anemometer at 80 m, that is mounted on a 116 m tall met mast (UTM zone 32 V 447 647 m E and 6 255 435 m N WGS84 datum) in the southeast direction (see Fig. 3). Mathematically, L MO is given as, where u * is the friction velocity, κ = 0.4 is the von Kármán constant, g is the acceleration due to gravity, θ v is the virtual potential temperature, . . . denotes time average, and w θ v (covariance of w and θ v ) is the virtual kinematic heat flux. u * is estimated as, where u w (covariance of u and w) and v w (covariance of v and w) are the vertical fluxes of the horizontal momentum.
As an initial validation of the accuracy and precision of the WindScanner, the 30 min mean wind speeds were compared with those obtained from the cup anemometer. Figure 4 shows that the WindScanner is very accurate (within 0.1 %) Atmos. Meas. Tech., 8, 729-740 to fit the cup anemometer measurements. The scatter using both methods is comparable to each other, but there is a slightly more scatter using the six-beam method for v 2 .
205 Fig. 6 shows the same as Fig. 5 but under neutral conditions. As for the unstable conditions, the six-beam method measures more turbulence, about 18 % for u 2 and 10 % for v 2 than the VAD method. The scatter using both methods is comparable to each other, with the six-beam method giving a slightly reduced scatter than the VAD method. Fig. 7 shows the same as Fig. 5 but under stable conditions. As for the unstable conditions, the 210 six-beam method measures more turbulence, about 19 % for u 2 and 4 % for v 2 than the VAD method. The scatter using both methods is comparable to each other, but there is a slightly more 13 Figure 5. Comparison of the turbulence statistics under unstable conditions between the WindScanner and the cup anemometer using two methods.
and precise (coefficient of determination, r 2 ≈ 0.9993) in measuring the mean wind speeds. For one 30 min period, the mean radial velocities measured by each of the six beams on the base of the scanning cone were fitted to Eq. (2) in a least squares sense to obtain the 30 min mean wind speed. This procedure was repeated for all 30 min periods. It is to be noted that the mean wind speed obtained using both methods (six-beam and VAD) is identical, since averaging the radial velocity for each beam, and then making a linear fit to obtain the u, v, an w components commute. Such an exercise provided enough confidence to proceed with deducing the turbulence measurements from the WindScanner using both methods.

Turbulence measurements
Two methods are used to deduce the turbulence statistics from the WindScanner measurements: 1. Six-beam method -for each 30 min period the measured S vector is used in combination with Eq. (5) to deduce the vector. Finally the vector is rotated in the mean wind direction for the respective 30 min period.
2. VAD method -within each 30 min period the v r measurements from every single cycle of the six beams are fitted in a least squares sense to Eq. (2) to deduce a 30 min time series of the wind field components. The vector is subsequently computed and rotated in the mean wind direction for the respective 30 min period. Figure 5 shows the comparison of the turbulence statistics derived from the WindScanner measurements using the six-beam and the VAD methods and those obtained from the cup anemometer under unstable conditions. It is clear that the six-beam method measures more turbulence, about 19 % for u 2 and 3 % for v 2 than the VAD method, where the orthogonal least-squares regression is used to fit the cup anemometer measurements. The scatter using both methods is comparable to each other, but there is a slightly more scatter using the six-beam method for v 2 . Thus under all stabilities the six-beam method is closer to the turbulence measurements carried out using the reference cup anemometer. There is however some probe volume averaging using both 215 methods, but is significantly larger for the VAD method.

Discussion
From Figs. 5-7 it is clear that using both methods the WindScanner measures more turbulence under stable conditions than under unstable and neutral conditions. This may be contrary to our intuitive 14 Figure 6. Comparison of the turbulence statistics under neutral conditions between the WindScanner and the cup anemometer using two methods. Figure 6 shows the same as Fig. 5 but under neutral conditions. As for the unstable conditions, the six-beam method measures more turbulence, about 18 % for u 2 and 10 % for v 2 than the VAD method. The scatter using both methods is comparable to each other, with the six-beam method giving a slightly reduced scatter than the VAD method. Figure 7 shows the same as Fig. 5 but under stable conditions. As for the unstable conditions, the six-beam method measures more turbulence, about 19 % for u 2 and 4 % for v 2 than the VAD method. The scatter using both methods is comparable to each other, but there is slightly more scatter using the six-beam method for u 2 .
Thus under all stabilities the six-beam method is closer to the turbulence measurements carried out using the reference cup anemometer. There is however some probe volume averaging using both methods, but this is significantly larger for the VAD method. The probe volume averaging can be observed clearly by comparing the radial velocity spectra, which can be observed in Fig. 6 in Mann et al. (2009), andFig. 4 in Sjöholm et al. (2009).

Discussion
From Figs. 5-7 it is clear that using both methods the Wind-Scanner measures more turbulence under stable conditions than under unstable and neutral conditions. This may be contrary to our intuitive understanding, because usually the turbulence scales are much larger under unstable conditions than under stable conditions . These results are also contrary to what has been observed by Sathe et al. (2011b) at the same site. However, it is to be noted that Sathe et al. (2011b) used the lidar measurements when the wind was blowing from the eastern direction, whereas in this work we use the measurements when the wind is blowing from the western direction. As described in Sect. 4, in the western sector there is a sudden change of roughness due to the transition from sea to land. As a consequence there is a development of the internal boundary layer (IBL). Also the growth of the IBL depends on atmospheric stability, where under unstable conditions the growth will be faster than under stable conditions. Panofsky and Dutton (1984)  under stable conditions . These results are also contrary to what has been observed by Sathe et al. (2011b) at the same site. However, it is to be noted that Sathe et al. (2011b) used the lidar measurements when the wind was blowing from the eastern direction, whereas in this work we use the measurements when the wind is blowing from the western direction. As described in section 4, in the western sector there is a sudden change of roughness due to the transition from 225 sea to land. As a consequence there is a development of the internal boundary layer (IBL). Also the growth of the IBL depends on atmospheric stability, where under unstable conditions the growth will be faster than under stable conditions. Panofsky and Dutton (1984) state that the growth of the 15 Figure 7. Comparison of the turbulence statistics under stable conditions between the WindScanner and the cup anemometer using two methods.
height of the boundary layer is proportional to the drag coefficient u * / u . And it is well known that the drag coefficient is larger for unstable stratification. Consequently the turbulence scales within 230 the IBL will be smaller as compared to those outside of it. It is then interesting to check whether the WindScanner measures more within the IBL under unstable conditions as compared to the stable conditions. corresponding to the maximum spectral energy, it is then clear that the peak of the v-spectra is shifted to the right for unstable conditions as compared to the stable conditions. It is not that clear for the u-spectra, however the shift of scales to larger wavenumbers under unstable conditions can still be observed. Thus L appears smaller under unstable conditions than under stable conditions for the measurements from the western sector used in this work. There is thus more probe volume 240 averaging under unstable conditions than under stable conditions. Hence the WindScanner attenuates to the drag coefficient u * / u . And it is well known that the drag coefficient is larger for unstable stratification. Consequently the turbulence scales within the IBL will be smaller as compared to those outside of it. It is then interesting to check whether the WindScanner measures more within the IBL under unstable conditions as compared to the stable conditions. Figure 8 shows the u and v spectra derived from highfrequency cup anemometer measurements under different stability conditions. If we define the characteristic length scale L as the length scale corresponding to the maximum spectral energy, it is then clear that the peak of the v spectra is shifted to the right for unstable conditions as compared to the stable conditions. It is not that clear for the u spectra, however the shift of scales to larger wavenumbers un-der unstable conditions can still be observed. Thus L appears smaller under unstable conditions than under stable conditions for the measurements from the western sector used in this work. There is thus more probe volume averaging under unstable conditions than under stable conditions. Hence the WindScanner attenuates the turbulence measurements more under unstable conditions than under stable conditions.
Another interesting observation is that using the VAD method the WindScanner does not measure more turbulence than the reference cup anemometer under any stability condition. This, too, does not agree with that observed by Sathe et al. (2011b), even though the same basic pulsed commercial lidar technology was also used in that work. It is likely due to the fact that in Sathe et al. (2011b) only four beams were used as opposed to six beams, and α was 60 • compared to 45 • used in this work. Therefore the turbulence statistics are not directly comparable with those obtained in Sathe et al. (2011b) even though the same basic commercial lidar was used. Due to the application of the least squares technique on the v r measurements in this work, there is significant volume averaging around the scanning circle, which is also observed in Sathe et al. (2011b) for a continuous-wave lidar.

Conclusions
An alternative so-called six-beam method is proposed in place of the standard VAD method to measure atmospheric turbulence using a ground-based wind lidar. The major difference between the two methods is that the six-beam method uses the measurement of the radial velocity variances, whereas the VAD method uses the high frequency measurement of the radial velocity transformed into Cartesian coordinates to deduce turbulence statistics. The scanning configuration of the six-beam method is optimized to minimize the sum of the random errors in the measurement of the components of the R matrix. In comparison to the reference cup anemometer the six-beam method measures between 85 and 101 % of the reference turbulence, whereas the VAD method measures between 66 and 87 % depending on atmospheric stability. The six-beam method thus overcomes partly the problem of significant probe volume averaging that is otherwise observed by the VAD method.
Furthermore two interesting observations have been made in this study. One is that, using both methods the WindScanner measures more turbulence under stable conditions than under unstable conditions, mainly due to the influence of the internal boundary layer (see Sect. 5). The other is that despite using the same underlying pulsed lidar technology as in Sathe et al. (2011b), the VAD method never measures more turbulence than the reference instrument as was observed in Sathe et al. (2011b) (see Sect. 5 for some explanation). It emphasizes the point that the VAD method is highly sensi-tive to the turbulence structure in the atmosphere, and one must avoid using it to measure atmospheric turbulence.
Future studies must certainly focus on tackling the probe volume averaging effect, which will further strengthen the arguments of using the six-beam method. Smalikho et al. (2005) have provided us with such a framework for pulsed lidars, whereas Mann et al. (2010) have demonstrated it for a continuous-wave lidar.