An Inter-Comparison Study of Multi- and DBS Lidar Measurements in Complex Terrain

: Wind measurements using classical proﬁling lidars suffer from systematic measurement errors in complex terrain. Moreover, their ability to measure turbulence quantities is unsatisfactory for wind-energy applications. This paper presents results from a measurement campaign during which multiple WindScanners were focused on one point next to a reference mast in complex terrain. This multi-lidar (ML) technique is also compared to a proﬁling lidar using the Doppler beam swinging (DBS) method. First- and second-order statistics of the radial wind velocities from the individual instruments and the horizontal wind components of several ML combinations are analysed in comparison to sonic anemometry and DBS measurements. The results for the wind speed show signiﬁcantly reduced scatter and directional error for the ML method in comparison to the DBS lidar. The analysis of the second-order statistics also reveals a signiﬁcantly better correlation for the ML technique than for the DBS lidar, when compared to the sonic. However, the probe volume averaging of the lidars leads to an attenuation of the turbulence at high wave numbers. Also the conﬁguration (i.e., angles) of the WindScanners in the ML method seems to be more important for turbulence measurements. In summary, the results clearly show the advantages of the ML technique in complex terrain and indicate that it has the potential to achieve signiﬁcantly higher accuracy in measuring turbulence quantities for wind-energy applications than classical proﬁling lidars.


Introduction
Over recent years, lidar (light detection and ranging) technology has quickly penetrated wind-energy applications. Especially in resource assessment, Doppler lidars are now widely used for wind measurements to predict the annual energy production of windfarms. In flat and homogeneous terrain, classical profiling lidars using the Doppler beam swinging (DBS) or velocity azimuth display (VAD) techniques achieve high accuracy [1]. Therefore, in recent years, their use has been adopted into national and international standards and guidelines [2,3].
The application of Doppler lidars in complex terrain, however, remains difficult and is often associated with systematic errors in mean wind-speed estimations [4,5]. If the flow is complex, the homogeneity assumption underlying the DBS and VAD techniques is often violated. First comparison studies have shown promising results in the capabilities of flow models in correcting this complex terrain error at individual sites [4,6]. However, they also indicated a high sensitivity to model parameterisations [6] and remain limited to moderately complex terrain. Moreover, the application of corrections leads to an increased uncertainty.
Another disadvantage of profiling lidars compared to mast-based measurements is that the estimation of turbulence characteristics, which are, e.g., relevant for wind-turbine loads, is not possible with the desired accuracy [7][8][9]. Most prominently, the estimation of the variance of the wind-vector components and the respective spectra suffer from systematic errors due to the effects of cross-contamination and volume averaging [9,10]. Newman and co-workers showed that some of the cross-contamination effect can be mitigated by applying a correction based on the measurement of the vertical variance [11]. However, their experimental results indicated that even after the application of the correction, a large deviation between lidar and mast-based measurements can remain. In general, a turbulence correction is difficult since it depends not just on lidar technology and measurement configuration but also the three-dimensional structure of the turbulence itself and, as a consequence, on atmospheric conditions [8,9]. New scanning strategies thus need to be developed. A promising example is the six-beam method [12,13], where turbulence parameters are directly derived from the radial velocity variances. Although the first encouraging experimental results have been reported in [12], other experiments indicated mixed results [11]. In complex terrain, the six-beam method might also be affected by an inhomogeneous turbulence field.
Probably the most intuitive solution to reduce the potential errors in lidar measurements is to set up multiple lidars such that the various lidar beams intersect at a single point from different directions (muti-lidar technique or ML). This strongly reduces the need for flow homogeneity in lidar measurements. Moreover, intersecting lidar beams probably offer the best approach to derive highly resolved time series of the three-dimensional wind-vector and turbulence statistics [14].
In flat and homogeneous terrain, first experimental results intersecting three lidar beams next to a sonic anemometer showed encouraging results for mean wind speeds [15][16][17] and turbulence statistics [14,16]. Mann and co-authors showed good agreement between second-order turbulence statistics of the radial wind velocities from multiple lidar systems and a sonic anemometer measurement [14] which was corrected using the spectral tensor model of [18]. Fuertes and co-authors achieved good agreement in second-order turbulence statistics between three intersecting lidars and a sonic anemometer to which they applied a simpler, pseudo-spatial filtering approach [16]. However, both studies were performed under idealised conditions in flat and homogeneous terrain and only report on very limited sample sizes. Newman and co-workers have also used three lidar beams intersected in a single point to derive turbulence statistics [19], but no adequate reference was available to assess the accuracy of their measurements. This paper presents measurements performed during the Kassel 2014 Experiment. The aim of the Kassel 2014 Experiment was to experimentally explore and demonstrate the advantages of ML measurements under conditions which are realistic for wind-energy applications in complex terrain. The experiment also served as a test bed for the complex terrain experiments in the New European Wind Atlas which will heavily rely on ML measurements [20]. Here, we present measurements and comparisons of four long-range WindScanners, which intersected at one point, a DBS lidar, and a sonic anemometer at a reference mast at a complex forested site. An analysis of the radial velocity statistics is used to examine the measurement accuracy of the individual instruments as well as the accuracy of the setup process. It is also used to provide insight into the path-averaging effect on turbulence measurements with the WindScanners. An analysis of the first-and second-order statistics of the horizontal wind speed contrasts error sources for the different lidar techniques. The ML-DBS-mast inter-comparison is especially valuable to analyse and quantify the benefits of the ML technique over the standard DBS lidar and experimentally compare their measurement behaviour in complex terrain. The WindScanners were located at distances ranging from directly next to the mast to 3740 m from the reference mast; they reflect a measurement scenario which is "realistic" in windfarm development scenarios in complex terrain.

Wind Vector Reconstruction and Scanning Strategies
Doppler lidars provide an estimate of reflectivity-weighted, radial velocities of aerosol particles. This velocity was applied as a good approximation for the radial wind velocity. In the context of wind profilers, the difference between Doppler velocity and real-wind velocity has been addressed by [21]. For Doppler lidars, the assumption seems justified as most of the backscattering for optical wavelengths is from particles which are small enough to be advected with the wind (e.g., [22]). Moreover, a homogeneous aerosol concentration within the probe volume needs to be assumed at this point. The measured radial velocity (v r i ) is thus a projection of the three-dimensional wind vector (v) onto the lidar beam where n i is the unit vector in the direction of the lidar beam and i denotes the individual beam and its location. In meteorological applications, it is convenient to write this equation as: where u is directed into the mean wind direction during the averaging interval, v is perpendicular to u in the horizontal and w is the vertical component of the wind vector. The projection onto the radial is defined by the elevation angle θ in addition to the angle between the mean wind direction and the beam φ. The azimuth angle (ψ) or beam-pointing angle is often used instead of φ in the literature to define (2). However, in this paper, the use of φ makes the notation in the later chapters more consistent. Assuming homogeneous flow (i.e., v i = v) between the measurement locations or-as is the case for ML configurations-the intersection of the beams in a single location, multiple measurements can be combined to derive v using matrix inversion where M = {n 1 , n 2 , . . .}. For a reconstruction of all components of v, at least three beams need to be combined. The horizontal wind speed (V h ) is then defined as: Classical profiling wind lidars usually rely on a scanning strategy, which employs a conical scan above the instrument and measures v r i at a fixed height and at different azimuth angles. The wind profiler lidar (manufacturer Leosphere, type Windcube v2), which is used as a reference instrument in this study, scans in DBS mode and uses four beams (with θ = 62 • ), which are arranged in two perpendicular planes ( Figure 1). An additional vertical beam is used to measure w directly. The resulting separation between two opposing beams in this setup is roughly equal to the measurement height.
Especially in complex terrain, the assumption of homogeneous flow across the different beams of a profiling lidar is often violated; significant, direction-dependent errors in the estimation of the mean can be introduced [4,5]. The main reason for the observed errors is thought to be the difference in w (the overbar denotes the mean) across the different beam locations [4,5]. Differences in the other components will of course also introduce an error. For the site of the current experiment, a deviation to mast-bound cup anemometer measurements of approx. −4% to +2% in V h has previously been reported [6]. Also the scatter in the inter-comparison was higher than in comparable experiments in flat and homogeneous terrain [1]. reported [6]. Also the scatter in the inter-comparison was higher than in comparable experiments in flat and homogeneous terrain [1]. In contrast to profiling lidars, the ML-technique uses a combination of multiple lidars, which are located around the measurement location and focus their laser beams to intersect in a single location ( Figure 1). Since multiple independent measurements of the radial velocity are made at the point of interest, the assumption of horizontal homogeneity of the mean flow among the different beam locations is not needed. The effect of a potential inhomogeneity within the probe volume of the ML measurement is expected to be much smaller and the beam separation problem inherent in a profiling lidar measurement is potentially overcome. Figure 1. Sketch of the multi-lidar (ML) and Doppler beam swinging (DBS) strategies next to a reference mast with a sonic as used in this study. Note that the angles and positions of the lidars differ from the instrument setup for reasons of illustrative clarity.

Second-Order Statistics of
One of the characteristics of wind-speed measurements from Doppler lidars is that, due to the measurement principle, is retrieved from a measurement within a probe volume. The width of the laser beam is usually negligible compared to the length of the probe volume along the path. The effect of the measurement volume on can thus be simplified to line averaging along the direction of the lidar beam. This can be expressed as where ( ) is the weighting function along the probe volume, is the centre of the probe volume and is the three-dimensional velocity field. Depending on the signal processing and shape of the emitted laser pulse, ( ) will vary and different forms have been suggested (e.g., [23,24]). Here we follow the concept that the weighting function can be estimated from a convolution of the form of the laser pulse ( ) and the gate-weighting function ( ), which results from the size and shape of the sampling gate-i.e., the internal signal processing [25]. Since the effects of line averaging depend on the scale of the turbulent eddies, it is convenient to express the problem in Fourier space. The Fourier transform of (5) can be written as [14,26]: where is the wave vector ( , , ) and ^ denotes the Fourier transform. For the WLS 200s (a scanning lidar manufactured by Leosphere) used in this study, both weighting functions can be approximated by a Gaussian function [27]. ( ) can thus be approximated by where is the standard deviation of ( ) and is the standard deviation of ( ). In contrast to profiling lidars, the ML-technique uses a combination of multiple lidars, which are located around the measurement location and focus their laser beams to intersect in a single location ( Figure 1). Since multiple independent measurements of the radial velocity are made at the point of interest, the assumption of horizontal homogeneity of the mean flow among the different beam locations is not needed. The effect of a potential inhomogeneity within the probe volume of the ML measurement is expected to be much smaller and the beam separation problem inherent in a profiling lidar measurement is potentially overcome.

Second-Order Statistics of v r
One of the characteristics of wind-speed measurements from Doppler lidars is that, due to the measurement principle, v r is retrieved from a measurement within a probe volume. The width of the laser beam is usually negligible compared to the length of the probe volume along the path. The effect of the measurement volume on v r lidar can thus be simplified to line averaging along the direction of the lidar beam. This can be expressed as where ϕ (s) is the weighting function along the probe volume, r is the centre of the probe volume and V is the three-dimensional velocity field. Depending on the signal processing and shape of the emitted laser pulse, ϕ (s) will vary and different forms have been suggested (e.g., [23,24]). Here we follow the concept that the weighting function can be estimated from a convolution of the form of the laser pulse ϕ p (s) and the gate-weighting function ϕ g (s), which results from the size and shape of the sampling gate-i.e., the internal signal processing [25]. Since the effects of line averaging depend on the scale of the turbulent eddies, it is convenient to express the problem in Fourier space. The Fourier transform of (5) can be written as [14,26]:v r lidar (k) =φ (k·n) n·V (k), where k is the wave vector (k 1 , k 2 , k 3 ) andˆdenotes the Fourier transform. For the WLS 200s (a scanning lidar manufactured by Leosphere) used in this study, both weighting functions can be approximated by a Gaussian function [27].φ (k) can thus be approximated bŷ where σ p is the standard deviation of ϕ p (s) and σ g is the standard deviation of ϕ g (s).
In practice, pulsed lidars use an accumulation time to derive the radial wind speed. We thus need to account for the temporal averaging. This can be done by adding the additional term ζ (k 1 ) = sinc (k 1 ·t·u/2), where t is the averaging time [14]. Assuming a homogeneous turbulence field along the probe volume and Taylor's hypothesis of frozen turbulence, we can express the line averaging in (6) on the measured turbulence spectrum, F v r (k 1 ), using the spectral density tensor Φ ij (k) =V (k)V * (k) as [14,26] where * denotes the complex conjugate. Summation over repeated indices is assumed. In analogy to what Kaimal and co-workers [26] suggested for sonic anemometers, we can now define a spectral transfer function between the spectrum measured by the lidar and a point measurement as The attenuation in v r lidar can then be calculated by integrating the spectra over k 1 . One of the difficulties with (9) is that a model for Φ ij (k) [18,26] is needed to estimate T v r (k 1 ). If n is in the direction of x 1 ,φ(k·n) becomes independent of k 2 and k 3 and they can be removed from the integral [26]. In this case, the effect of line averaging can be expressed as This special case has the advantage of removing the need for a model of Φ ij (k). The line averaging effect can now be modelled by applying the theoretical formulation in (10) to a point measurement. If the centre of the lidar range gate and the point measurement are co-located, the only assumption that we have to make is Taylor's frozen turbulence hypothesis along the probe volume of the lidar and the form ofφ (k). In field measurements, the lidar beam will rarely be exactly parallel to the streamlines. In fact, during the Kassel Experiment 2014, the beams of all WindScanner devices were pointed upwards into the atmosphere to some degree. Therefore, (9) would be the exact expression. However, in the later analysis, we assume that (10) is still a good approximation if the angle between the mean streamlines is small.

Second-Order Statistics from Multiple Lidar (ML) Beams
In wind-energy applications V h 2 or u 2 and, in more complex terrain, also v 2 and w 2 are of key interest in a site suitability analysis [28]. The respective spectra are also important as they directly feed into models commonly used for load simulation in the design process of wind turbines [29,30].
In case the lidar beam is not pointing exactly in the direction of the wind component of interest, we need to use a combination of multiple beams (Section 2.1). In this situation, in addition to the line averaging, the location of the different lidar beams to each other becomes important. For a mathematical formulation of the effect of beam separation, it is convenient to define the distance of the lidar measurement from the point for which the measurement is assumed to be representative of where x is the location of the centre of the averaging path. Adding the measurement location to (6) and including the temporal averaging,v r lidar can now be rewritten as [26] v r lidar (k, d) =ζ (k 1 ) e ik·dφ (k·n) n·V (k).
The term e ik·d basically decorrelates the wind field which is sensed at the different measurement locations for high wavenumbers by introducing a phase shift. By combining (3) and (12), a matrix N can now be defined, which accounts for the effects of line averaging and path separation in the vector reconstruction. The entries of N are and its inverse J = N −1 . The vector reconstruction in the Fourier domain then becomeŝ The spectrum of e.g., the u component derived from three lidar beams can be written as where again the summation for repeating indices i and j is assumed. Sathe and co-workers call J "weighting factors," because they define how Φ ij (k) is weighted in the vector reconstruction [9]. The effect of the beam separation on the time series-derived wind vector is twofold. The decorrelation leads to a damping of the contribution of the spectral tensor, Φ ij (k), in the direction of the desired component (e.g., i = j = 1 for u) for high wavenumbers. However, the decorrelation effect is also present for the other components of Φ ij (k). As a result, some of the spectral energy of these components is folded onto F u (k 1 ). This effect is often referred to as cross-contamination in the literature [9,26], because the desired spectrum is contaminated by the other terms of Φ ij (k). As explained in Section 2.2. the effect of the probe volume averaging is mainly attenuation at small wavenumbers.
Using the spectral tensor model of [18], this effect for two different conically scanning lidars was investigated [9,10]. It was shown that the combined effects of beam separation and line averaging strongly depend on the beam configuration and the probe volume, which define J, but also the atmospheric conditions i.e., Φ ij (k). Experimental evidence confirms these results and an underestimation as well as an overestimation of u 2 and v 2 has been reported [8,9,11,12]. It should also be noted that, in complex terrain, the assumption of a statistically homogeneous turbulence field implicitly inherent in (15) is likely to be violated and additional errors might be introduced into the turbulence measurement.
If the range gates which are used for the wind-vector reconstruction intersect at the same location, e ik·d = 1 and the effect of beam separation vanishes. Also, the distance for which a homogeneous turbulence needs to be assumed is greatly reduced. It is interesting to note that the line averaging still introduces some cross-contamination. The effect is small compared to the spectral probe-volume averaging, but may be increased if the probe volume differs between the beams.
In general, the same principles concerning the line averaging and path separation apply for a sonic anemometer as well. However, the sonic anemometer used in this study uses intersecting measurement paths which are 2-3 orders of magnitude smaller than the path averaging of the lidars. Thus these effects are neglected and the sonic anemometer is treated as an approximation of a point measurement in the analysis in this study.

Experimental Setup: The Kassel 2014 Experiment
The data presented in this paper was collected during the Kassel 2014 Experiment at Rödeser Berg-a hill in central Germany close to Kassel-during summer 2014 (3 July-17 August 2014). Five long-range WindScanners (scanning lidars based on the Windcube WLS200s) were set up around the hill (EE, SE, SW, WW, NN in Figure 2). A sixth scanner (MA) was placed on the ridge of the hill directly next to a 200 m high meteorological mast. An additional Windcube v2 (DBS lidar) was operated at a few meters' distance from the mast. The measurements from a sonic anemometer mounted on the mast serve as a reference in this study. around each WindScanner, with the exception of the MA WindScanner, at distances of approx. 100 m, three 2 cm thick surveying stakes were installed, separated by 120° in azimuth. The absolute positions of the WindScanners, stakes, and reference mast top were surveyed using a differential GPS and a theodolite. These positions were acquired with cm accuracy. The stakes were then mapped using the CNR-Mapper [15]. After the calibration, the reference mast top was mapped by each WindScanner to validate the setup procedure, and the static laser beam pointing accuracy and the sensing range accuracy were estimated to be 0.05° and 1 m, respectively. Due to the terrain configuration, reference mast structure and surrounding vegetation, the deployment steps, which are described above, could not be applied completely to the MA WindScanner. Hence, a single surveying stake and the branch of a tree nearby were used for the CNR mapping. Both hard targets were located about 50 m from the WindScanner and separated by 180° in azimuth. However, since the MA system was located directly next to the reference mast, the static pointing accuracy and sensing range accuracy could not be validated by mapping the mast top. The WindScanners were coordinated and synchronized by a remote master computer [15], which communicated with the WindScanners using a 3G network [32].
This paper analyses data from measurements when all WindScanners were staring at a single point close to a Gill HS50 sonic anemometer mounted on the reference mast at a height of 188 m (agl).
The reference anemometer is mounted on a foldable boom with a length of 5.4 m to minimise the effects of the mast structure. The wind sector 100°-160° is excluded from the analysis to remove data which is affected by the mast shadow. According to [2], the influence of the mast structure is estimated to be less than 0.5%. In a wind tunnel test the directional deviation of the Gill HS50 is well below ±1%. The inclinometer accuracy is given as ±0.15°. The direction-dependent error of the measurement of the horizontal wind speed is thus estimated to be below ±1%. The sampling frequency of the Gill HS50 was set to 50 Hz.
Although there is data from six available WindScanners, we restrict the analysis to four locations (EE, SW, SE and MA). These locations have the advantage that they have the longest time series in the analysed scanning mode. Details about the setup of the four scanners can be found in Table 1. The reference mast is erected on a small clearance within a forest which, depending on direction, stretches from 0.4 to 2.5 km around the mast. The wider surroundings are characterised by a patchy landscape of mainly agricultural land use, villages and forested hills ( Figure 2). The height of the hill above the surrounding terrain is approx. 100-200 m depending on direction. More details about the measurement site can be found in [6].
During the deployment of the WindScanners the procedures for calibration of the home position, levelling, and sensing range as outlined in [15] were applied. For the purpose of these procedures around each WindScanner, with the exception of the MA WindScanner, at distances of approx. 100 m, three 2 cm thick surveying stakes were installed, separated by 120 • in azimuth. The absolute positions of the WindScanners, stakes, and reference mast top were surveyed using a differential GPS and a theodolite. These positions were acquired with cm accuracy. The stakes were then mapped using the CNR-Mapper [15]. After the calibration, the reference mast top was mapped by each WindScanner to validate the setup procedure, and the static laser beam pointing accuracy and the sensing range accuracy were estimated to be 0.05 • and 1 m, respectively.
Due to the terrain configuration, reference mast structure and surrounding vegetation, the deployment steps, which are described above, could not be applied completely to the MA WindScanner. Hence, a single surveying stake and the branch of a tree nearby were used for the CNR mapping. Both hard targets were located about 50 m from the WindScanner and separated by 180 • in azimuth. However, since the MA system was located directly next to the reference mast, the static pointing accuracy and sensing range accuracy could not be validated by mapping the mast top. The WindScanners were coordinated and synchronized by a remote master computer [15], which communicated with the WindScanners using a 3G network [32].
This paper analyses data from measurements when all WindScanners were staring at a single point close to a Gill HS50 sonic anemometer mounted on the reference mast at a height of 188 m (agl).
The reference anemometer is mounted on a foldable boom with a length of 5.4 m to minimise the effects of the mast structure. The wind sector 100 • -160 • is excluded from the analysis to remove data which is affected by the mast shadow. According to [2], the influence of the mast structure is estimated to be less than 0.5%. In a wind tunnel test the directional deviation of the Gill HS50 is well below ±1%. The inclinometer accuracy is given as ±0.15 • . The direction-dependent error of the measurement of the horizontal wind speed is thus estimated to be below ±1%. The sampling frequency of the Gill HS50 was set to 50 Hz.
Although there is data from six available WindScanners, we restrict the analysis to four locations (EE, SW, SE and MA). These locations have the advantage that they have the longest time series in the analysed scanning mode. Details about the setup of the four scanners can be found in Table 1.

Data Treatment, Quality Control and Coordinate Systems
Quality control of lidar data is usually done by using a lower carrier-to-noise ratio (CNR) threshold to identify periods when not enough backscatter is received to derive a reliable wind-speed estimate. If a structure is present close to the measurement volume, returns from this hard target can also contaminate the signal. Often this leads to an increased CNR and/or to an increased scatter in the observed CNR values. Careful analysis of the WindScanner data, which was collected during the presented experiment, however, showed that also periods during which the CNR fell into the expected CNR range (here −27.5 and −5 dB were used) and had a "normal" variability sometimes exhibited spurious radial wind-speed statistics. Therefore, two additional filter criteria were imposed on the WindScanner data.
The first sorts the radial wind velocities of each averaging interval according to their magnitude. Outliers are then removed if the gap between them and the rest of the data exceeds 1 m·s −1 . This method is useful in situations with a relatively high radial velocity which is contaminated with hard target returns, as these are often characterised by a radial wind speed close to zero. The filter, however, becomes ineffective if the mean radial wind speed is close to zero-i.e., winds perpendicular to the lidar beam or low wind speeds. The second filter criterion therefore checks the time series during the averaging interval for its temporal consistency. If the time series after the application of the first filter still exhibits a difference between two adjacent measurements which exceeds 1 m·s −1 , the respective period is removed completely from further analysis. This threshold seems justified because of the relatively high sampling frequency (0.5 Hz) and the damping effect of the lidar probe volume on small-scale turbulence. The exact reason for the problems with the hard targets could not fully be clarified. At the beginning of the experiment, the WindScanners were focused within approx. 1 m above the reference sonic. The following points may then have led to interference with the mast structure: instability of the WindScanner systems (sinking into soil after installation, vibration of the systems due to e.g., wind), switching of the staring mode to other patterns inducing pointing error due to backlash of the internal gearing, and/or movement of the mast structure. For future inter-comparison experiments, it is thus recommended to increase the distance to the reference instrument/mast. The problem with the hard targets was not present in the data of the Windcube v2 and only a CNR filter was applied. The values (−22 and −5 dB) slightly differ from the WindScanners as the system configuration is different.
The analysis in this paper is based on 10 min averaging intervals, which is the standard in wind-energy applications. Before the analysis of the wind-vector components, a rotation into the mean wind direction (v = 0 m·s −1 ) is performed. A second rotation forcing w = 0 m·s −1 as done in many other micrometeorological studies (cf. [33]) is not carried out as the aim of the study is mainly a comparison of the different measurement techniques. For this same reason, corrections which are often applied to eddy-covariance measurements (cf. [34]), such as correcting for high-and low-frequency loss, are not applied to the sonic anemometer measurements. In a previous study, Mauder and co-authors reported deviations between a Gill HS and other sonic anemometers, which was attributed to the Gill HS [35]. However, in their study, the internal correction/calibration was turned off, which removes the sensor head correction and the instrument-dependent calibration. In the present study, the internal corrections were turned on and the problems reported in [35] are assumed to not to be present in the collected data.

Radial Velocity Components
In a first step, the radial components of the individual WindScanners (v r lidar ) are compared to the projection of the sonic measurement on the respective WindScanner beam (v r sonic ). This allows an evaluation of the performance of the individual instruments and the accuracy of the setup process. Moreover, some of the complexities in the inter-comparison of the second-order moments, which arise from the combination of multiple lidar beams in three-dimensional space, can be simplified when the radial components of the sonic and lidar wind speed are compared directly. In Section 3.2., we then turn to the wind statistics which are derived from the combination of multiple WindScanner devices.

First-Order Statistics
The inter-comparison of the mean radial components reveals an excellent agreement between the sonic and all WindScanner devices. A two-parametric linear regression between the WindScanners located around the mast and the reference sonic produces slopes (m) and intersects (b) which are very close to 1 and 0 m·s −1 ( Table 2). The goodness of fit for the linear regression for these instruments is R 2 ≥ 0.998. The linear-regression statistics for MA are slightly worse. However, this instrument is pointing vertically into the atmosphere and w is usually small in the atmospheric boundary layer. In fact, v r MA only varied between approx. ±1 m·s −1 and the root-mean-square deviation (RMSD) is comparable to the other WindScanners.
The availability of all WindScanners is affected by hard targeting (mast structure, boom or guying). As a result, the number of valid measurement periods, especially for the MA location, is relatively low. A conclusion on how the measurement range affects the data availability is therefore unfortunately not possible. Table 2. Summary of the inter-comparison of the radial velocities of the WindScanners and the sonic anemometer at 188 m; RMSD indicates the root-mean-square deviation; R 2 is the goodness of the fit, m is the slope and b the intersect of the of the linear regression.

Second-Order Statistics
For the locations SE, SW and EE, the linear-regression statistics between v 2 r lidar (measured by the WindScanners) and the respective variance derived from the projection of the sonic measurement onto the direction of the respective WindScanner beam (v 2 r sonic ) are very similar ( Table 2). The slope indicates a significant underestimation of v 2 r sonic by the WindScanner measurements. This is likely to be due to the relatively large probe volume. In contrast, the slope of the linear regression for the MA system is significantly closer to 1. Here, the smaller range gate size seems to more than compensate for the smaller turbulent length scales which are usually found for w. R 2 is similar for all systems. The values are slightly worse than for the first-order statistics, except for the MA location, where R 2 is slightly higher. The reduced values in R 2 are probably partly caused by the reduced range of the observed variances when compared to the mean values. Another reason for the increased scatter is due to the fact that line averaging does not necessarily scale linearly with the observed variance. It is rather direction dependent on and selective of small wavenumbers. This makes the attenuation dependent on the form of Φ ij (k) (9), which varies with atmospheric conditions-i.e., under stable conditions, a larger portion of the variance will be contained in small eddies as opposed to under unstable conditions. Also, since a finite time interval is used (10 min), some scatter will be introduced by the stochastic distribution of the eddy sizes in the individual period. To further investigate the probe-volume averaging effect on v 2 r lidar periods with wind directions close to parallel (±15 • ) to the azimuth angle of the individual WindScanner beams, SE, SW and EE were selected. The spectral filter function defined by (7) and (10) was applied to v r sonic (Figure 3). The application of the filter function to the sonic data brings the slope of the linear regression between v 2 r sonic and v 2 r lidar closer to unity and slightly improves R 2 for all systems. Interestingly, this improvement is also observed for the EE system, which has quite a large inclination angle into the atmosphere (θ = 22.3 • ).
A comparison of the average spectra of v r of the sonic and WindScanner measurements reveals a clear spectral attenuation for SE, SW and EE spectra from approx. k 1 > 10 −2 m −1 (Figure 4a-c). This corresponds nicely to what we would expect for the spectral transfer function (T v r ) from theoretical considerations using the parameters in Table 1 and Equations (7) and (10) (Figure 4e). In contrast to the other systems, the spectrum of MA does not exhibit a clear spectral attenuation, which again is likely to be due to the smaller probe volume. Instead, the spectral energy at the high wavenumber tail for MA is actually higher than for the reference sonic. This behaviour can also be observed for the SW system (Figure 4c). Also T v r starts to increase and deviate from its theoretical form at appox. k 1 > 10 −1 m −1 for all systems (Figure 4d). This effect becomes worse if measurements with small variances are included (data not shown) and is likely to be attributed to instrumental noise.
introduced by the stochastic distribution of the eddy sizes in the individual period. To further investigate the probe-volume averaging effect on periods with wind directions close to parallel (±15°) to the azimuth angle of the individual WindScanner beams, SE, SW and EE were selected. The spectral filter function defined by (7) and (10) was applied to (Figure 3). The application of the filter function to the sonic data brings the slope of the linear regression between and closer to unity and slightly improves R 2 for all systems. Interestingly, this improvement is also observed for the EE system, which has quite a large inclination angle into the atmosphere (θ = 22.3°). after the application of (10). The sonic data was filtered with ( ) to separate out the effects of the temporal averaging.
A comparison of the average spectra of of the sonic and WindScanner measurements reveals a clear spectral attenuation for SE, SW and EE spectra from approx. > 10 (Figure 4a-c). This corresponds nicely to what we would expect for the spectral transfer function ( ) from theoretical considerations using the parameters in Table 1 and Equations (7) and (10) (Figure 4e). In contrast to the other systems, the spectrum of MA does not exhibit a clear spectral attenuation, which again is likely to be due to the smaller probe volume. Instead, the spectral energy at the high wavenumber tail for MA is actually higher than for the reference sonic. This behaviour can also be observed for the SW system ( Figure  4c). Also starts to increase and deviate from its theoretical form at appox. > 10 m for all systems (Figure 4d). This effect becomes worse if measurements with small variances are included (data not shown) and is likely to be attributed to instrumental noise.   Up until now, only a few experimental studies investigating v 2 r lidar and its spectrum against reference measurements have been undertaken. Most interesting in regards to pulsed lidars, Mann and co-authors reported a similar agreement between the observed and theoretical forms of T v r . However, in contrast to this study, they only used five selected periods and the lidars were set up in the direct vicinity of the reference mast. As this leads to large elevation angles, they had to use a spectral tensor model to derive T v r . Also, they used a simpler, triangular-shaped weighting function. The influence of measurement noise at the high-wave-number end of the spectra has also been reported for continuous wave lidars [36,37] and deserves further investigation for the WindScanner systems. Figure 4e indicates that the concept of line averaging (9 and 10) can be used to explain the majority of the differences in the observed second-order statistics between the sonic and the WindScanners.

First-Order Statistics of the Horizontal Wind Speed
For resource estimation in the wind industry, it is standard to use the scalar horizontal velocity (V h ). This approach is also followed here. Although in many remote sensing applications it is common practice to work with the vector average of the wind velocity, publications concerning the accuracy of lidars in complex terrain often use the scalar mean [4,6]. While vector averaging produces slightly different numbers, the general patterns that can be observed remain the same and the differences are very small.
The comparison of V h measurements of the ML combination using three lidars and the sonic reveals an excellent agreement between the two techniques as would be expected from the comparison of v r (Figure 5a). Also, the ML technique using two lidars and neglecting the vertical component in the vector reconstruction only produces marginally different results for the scanner combination SE/SW (Figure 5b). For the SW/EE combination, the scatter slightly increases and we can see a few points which significantly deviate from the linear relationship (Figure 5d). The observed differences between the different dual-lidar configurations are likely to be at least partly caused by the setup of the instruments.
One characteristic which sets the different dual-lidar measurements apart is the difference in azimuth angle between the WindScanner beams. For SW/EE, the angle is smallest (43 • ) which causes the strongest error propagation from the radial wind speed onto the wind vector when solving Equation (3) (cf. [38]). The angle for SE/EE is slightly larger (49 • ) and the SE/SW combination has an angle difference of 88 • , which is close to the ideal 90 • .
The DBS measurement with the Windcube v2 exhibits more scatter and the linear regression significantly deviates from 1:1 when compared to the sonic (Figure 5e). Especially for low wind speeds, an overestimation of the V hsonic by the DBS lidar can be observed.
In general, the linear-regression statistics for the DBS lidar are similar to what has been reported at Rödeser Berg in earlier measurement campaigns [6] and are similar to what has been reported from other lidar-mast inter-comparisons in complex terrain [39] albeit significantly worse than in flat terrain (e.g., [1]). In contrast, the statistics for the ML combinations are similar to the results which we regularly observe during DBS lidar calibrations in flat and homogeneous terrain. The "bump" in the DBS data at low wind speeds might be associated with turbulence effects which are more important for low wind speeds. Moreover, thermal effects on the flow might create increased vertical velocities (relative to the horizontal wind speed) which are likely to spatially vary in the complex surroundings at Rödeser Berg. In combination with the large scanning volume of the DBS lidar at 188 m this would introduce an increased complex terrain error.
To our knowledge no statistics from a comparably long-term data set for assessing the performance of ML measurements to derive the mean wind speeds have been published so far. Previously published results of mast-lidar inter-comparisons either focused on the evaluation of individual measurement periods [16] or employed scanning strategies that focused more on assessing the spatial heterogeneity of the wind field than on evaluating ML measurements against a reference [17,40]. The ML-mast inter-comparison in the latter two references is slightly worse, but due to the different focus of the scanning strategy, is likely to not be directly comparable. A recently published study by Newman only compares co-located DBS and ML measurements [19]. The sonic anemometer used in this study only measures at significantly lower altitudes than the lidars and is only used to give an indication of the wind field below the lidar measurements. This precludes a quantification of the advantages of either of the techniques over the other in their study or a comparison to the results presented.
One of the typical observations for profiling lidars in complex terrain is that the deviation in V h between the lidar and the mast measurements varies with the shape of the terrain and thus wind direction [4,6]. This behaviour can also be observed for the DBS data collected during the Kassel Experiment 2014 (Figure 6). The pattern clearly reflects the shape of the terrain and is very similar to what has been observed in earlier inter-comparisons between cup anemometers and a DBS lidar at Rödeser Berg [6]. If the flow is directed across the ridge of the hill (aprox. 50 • -100 • and 200 • -270 • ), the DBS lidar underestimates V h . This is a typical observation for convex flow as observed over e.g., ridges or hilltops [4,5]. In contrast, the sectors showing an overestimation of V h by the DBS lidar (approx. 290 • -20 • and 160 • -180 • ) correspond to along the ridge flow. Here, the mast is located on the slope in front of/behind the highest point on the ridge and the terrain is slightly concave.
This might also cause a concave flow pattern which results in an overestimation of the horizontal wind speed by the DBS lidar. Some evidence for this can be found in the change of sign of the flow angle measured by the sonic between 290 • -20 • and 160 • -180 • . Whereas for the first sector (mast in front of the hill top) the median flow angle is 1.7 • (directed upwards), it is −2.2 • (directed downwards) for the second sector (mast behind the hill top). It should be noted, however, that the real flow might be more complex than this simplified explanation. A more detailed discussion and modelling results of directional errors of DBS lidar measurements at Rödeser Berg can be found in [6]. To our knowledge no statistics from a comparably long-term data set for assessing the performance of ML measurements to derive the mean wind speeds have been published so far. Previously published results of mast-lidar inter-comparisons either focused on the evaluation of individual measurement periods [16] or employed scanning strategies that focused more on assessing the spatial heterogeneity of the wind field than on evaluating ML measurements against a reference [17,40]. The ML-mast inter-comparison in the latter two references is slightly worse, but due to the different focus of the scanning strategy, is likely to not be directly comparable. A recently published study by Newman only compares co-located DBS and ML measurements [19]. The sonic anemometer used in this study only measures at significantly lower altitudes than the lidars and is only used to give an indication of the wind field below the lidar measurements. This precludes a quantification of the advantages of either of the techniques over the other in their study or a comparison to the results presented.   increase the comparability to other complex terrain measurements reported in the literature. Also, bins with n < 5 are not displayed; error bars denote ± one standard deviation.

Second-Order Statistics of the Horizontal Wind Vector Components
This section presents the second-order statistics (variances and spectra) of the horizontal wind components and . The second-order statics for the component can be found in Section 3.1.2. (MA WindScanner). Figure 7 displays the scatter plots and linear-regression statistics for the different ML combinations and the DBS lidar when compared to the sonic measurements. The differences between the DBS and ML, but also among the different ML combinations, are quite large. Almost identical linear-regression statistics can be observed for the SW/SE/EE and SW/SE combinations (Figure 7a,b). The R 2 values of the linear regression for ′ and ′ are similar to the statistics of (Table 2 and Figure 3). Also the slope indicates a similar underestimation of ′ (and ′ ) of about Figure 6. (a) mean of the directional deviation in V h between the sonic and the different lidar configurations; (b) same as (a) but for the SE/EE combination with and without the correction for w; see text for details; periods during which V h sonic < 4 m·s −1 were excluded from the comparison to increase the comparability to other complex terrain measurements reported in the literature. Also, bins with n < 5 are not displayed; error bars denote ± one standard deviation.
On average, the DBS-mast deviations are slightly shifted (1%-2%) towards an increased V hlidar (relative to V hsonic ) when compared to the earlier inter-comparison to cup anemometry. This might be caused by the different reference anemometers or the different DBS instruments used in the present study, although the reason could not be fully clarified.
The deviation between V hlidar and V hmast for the ML measurements is smaller than for the DBS system for most sectors and WindScanner combinations ( Figure 6). Especially for the WindScanner combinations SW/SE/EE and SE/SW, the differences are small for all directions and only vary between approx. −2% and 2%. For the SE/EE combination, a relatively large positive deviation of the ML system in comparison to the sonic wind speed can be observed between approx. 350 and 25 • , which even exceeds the DBS-sonic deviations. This is likely due to the fact that within this range, the angle between the mean wind direction and the beam directions of SE and EE are large and the difference in beam directions is small (cf. [38]). This configuration leads to an increased error propagation in the estimation of u that largely dominates V hlidar -i.e., the azimuth angles and flow direction are unfavourable for the retrieval of V hlidar . For flow angles which are directed in a more parallel manner to the orientation of the lidar beams, the error propagation for v is expected to be large. The effect on V hlidar , however, is then much smaller.
A second effect which can affect the quality of the measurements of SE/EE combination is the relatively large elevation angle of the EE system compared to the other systems ( Table 1).
The assumption, w = 0 m·s −1 , thus results in a contamination of V hlidar by w. Due to the low availability, an inclusion of the MA system which could correct for the influence of w in the directional analysis was not possible. Therefore, a direction-dependent correction factor was calculated using w sonic . Equation (3) is first solved using v r SE , v r EE and w sonic . In a second step, (3) is solved only using v r SE and v r EE . The ratio of the two results is then applied as a correction factor to V hlidar . The correction reduces the deviations in the 350 • -25 • sector. The effect is, however, relatively small. Therefore, the assumption w = 0 m·s −1 does not seem to be the main cause of the observed deviations. Some of the deviations between the sonic and the lidar measurements might also be caused by the flow distortion of the mast. However, we estimate this effect to be small (c.f. Section 2.4.).

Second-Order Statistics of the Horizontal Wind Vector Components
This section presents the second-order statistics (variances and spectra) of the horizontal wind components u and v. The second-order statics for the w component can be found in Section 3.1.2. (MA WindScanner). Figure 7 displays the scatter plots and linear-regression statistics for the different ML combinations and the DBS lidar when compared to the sonic measurements. The differences between the DBS and ML, but also among the different ML combinations, are quite large. Almost identical linear-regression statistics can be observed for the SW/SE/EE and SW/SE combinations (Figure 7a,b). The R 2 values of the linear regression for u 2 and v 2 are similar to the statistics of v 2 r (Table 2 and Figure 3). Also the slope indicates a similar underestimation of u 2 (and v 2 ) of about 20%. Looking at the spectra of SW/SE/EE and SW/SE, one can clearly see the effects of line averaging which causes the reduction in u 2 and v 2 (Figure 8a,b). As for v r , an attenuation in the spectral density from approx. k 1 = 10 −2 m −1 for u and v can be observed. The line averaging effect thus also causes the reduced u 2 and v 2 in SW/SE/EE and SW/SE.
The scatter in the inter-comparison between the SW/EE and SE/EE is larger than for SW/SE/EE and SW/SE and the slope of the linear regression shifts close to 1. One of the reasons for the increased scatter in the ML combinations using two WindScanners including the EE system might be that the EE system has an inclination angle which is not close to 0 • . In the wind-vector reconstruction for two scanners, the w component is assumed to be 0 m·s −1 . Some of w 2 will thus contaminate the measurement of u 2 and v 2 . Moreover, as discussed earlier, the direction of the flow in relation to the beam directions will also influence the error propagation for the individual components, especially if the angle between the beams is small [38]. In the experimental setup it is, however, difficult to separate these effects as they are likely to be correlated.
This observation suggests that for the ML method, an adequate positioning of the lidar devices is even more important for measuring turbulence quantities than for the mean wind speed.
In general, random measurement errors in the time series will always lead to an increased variance. Therefore, both the contamination by the w component and the unfavourable angles in the lidar combination lead to increased u 2 and v 2 in the lidar measurement. The slope values close to the 1:1 line are thus likely a result of the combination of these measurement errors and the line averaging, which act in opposing directions. 20%. Looking at the spectra of SW/SE/EE and SW/SE, one can clearly see the effects of line averaging which causes the reduction in ′ and ′ (Figure 8a,b). As for , an attenuation in the spectral density from approx. = 10 m for and can be observed. The line averaging effect thus also causes the reduced ′ and ′ in SW/SE/EE and SW/SE.   The scatter in the inter-comparison between the SW/EE and SE/EE is larger than for SW/SE/EE and SW/SE and the slope of the linear regression shifts close to 1. One of the reasons for the increased scatter in the ML combinations using two WindScanners including the EE system might be that the EE system has an inclination angle which is not close to 0°. In the wind-vector reconstruction for two scanners, the component is assumed to be 0 m·s −1 . Some of ′ will thus contaminate the measurement of ′ and ′ . Moreover, as discussed earlier, the direction of the flow in relation to the beam directions will also influence the error propagation for the individual components, especially if the angle between the beams is small [38]. In the experimental setup it is, however, difficult to separate these effects as they are likely to be correlated.
This observation suggests that for the ML method, an adequate positioning of the lidar devices is even more important for measuring turbulence quantities than for the mean wind speed.
In general, random measurement errors in the time series will always lead to an increased variance. Therefore, both the contamination by the component and the unfavourable angles in the lidar combination lead to increased ′ and ′ in the lidar measurement. The slope values close to the 1:1 line are thus likely a result of the combination of these measurement errors and the line averaging, which act in opposing directions. to the mast; solid lines are the sonic spectra; dashed lines are the lidar spectra; sonic time series have been aggregated to 0.5 Hz (a and b) and 0.89 Hz before calculation of the spectra; only periods with u sonic > 4 m·s −1 and u 2 sonic > 0.2 m 2 ·s −2 were used in the spectral averaging; u * is the friction velocity computed from the sonic anemometer measurements; the black line indicates the theoretical −2/3 slope in the inertial subrange; n indicates the number of periods used for averaging.
The comparison between the DBS system and the sonic indicates a clear overestimation of u 2 and v 2 by the DBS system (Figure 7e and Table 3). Also, the R 2 values are lower than for all ML combinations. The spectra of the DBS system show a significant contamination of the variance measurements across almost the whole wavenumber range. Only for very large k 1 are the spectra for the u component of the DBS lidar and the sonic close together. For the u component, a local maximum can be observed around approx. k 1 = 10 −1.8 m −1 . The v spectrum of the DBS lidar is higher than the sonic spectrum for the whole wave number range and a (smaller) second maximum can be observed around k 1 = 10 −1.4 m −1 . For higher wave numbers, the spectra then quickly drop off.
The exact shape of the measured spectra of a Windcube does not only depend on the atmospheric conditions and measurement heights, but will also vary with the flow direction towards the scanning geometry of the Windcube [10]. It is thus difficult to explain all features of the measured spectra.
The strong positive deviations from the sonic spectra in this range are likely to be caused by the cross-contamination effects between u 2 , v 2 and especially by the contamination of the horizontal variances by w 2 . The short region where the spectra are almost horizontal at the high wavenumber end is likely to be caused by the algorithm which derives the wind vector from the measurements of v r . In analogy with the internal software of the Windcube v2, in this paper the wind vector was calculated whenever a new measurement of v r was available-i.e., at a frequency of approx. 0.89 Hz. One full rotation, however, takes five measurements of v r and thus acts similarly to a moving average filter [41]. Table 3. Summary of the linear-regression statistics for u 2 and v 2 for different lidar configurations; values in brackets are statistics for which periods with u sonic < 4 m·s −1 have been excluded; R 2 is the goodness of the fit, m is the slope and b the intersect of the of the linear regression. The relatively strong overestimation in the measured variances by the Windcube compared to some other studies using similar instruments (Windcube v1 and v2) [9,19] is likely to be related to the great measurement height (188 m) and thus the large distance between the different beams. In parts, it might also be caused by a dominance of unstable conditions which are more frequent during the summer season and the inhomogeneous flow caused by the terrain complexity.

SE SW EE
In Table 3, we also provide the statistics for periods when u sonic > 4 m·s −1 , as this filter criterion is often used in applications and studies in the wind-energy sector, thereby allowing a more direct comparison to other experiments. For the ML combinations SW/SE/EE, SW/SE and SE/EE, there is only a small change in the linear-regression statistics. For the Windcube v2 and the SW/EE, which have lower R 2 values, R 2 significantly improves after the application of the u sonic > 4 m·s −1 criterion.

Conclusions
Two of the biggest challenges for ground-based remote sensing in wind energy are the systematic errors in the mean wind speed which can be introduced by complex flow in complex terrain and the inability of traditional profiling devices to accurately measure turbulence. This paper experimentally explores the improvements in these two areas which can be made using the multi-lidar (ML) technique. For this purpose, results from the Kassel 2014 Experiment are presented, during which multiple WindScanners and a Doppler beam swinging (DBS) lidar were operated surrounding a reference mast in complex terrain.
The comparison of the radial wind velocities between the WindScanner measurements and a reference sonic shows that, given an accurate setup, an excellent agreement with the reference sonic could be achieved for individual systems, even if they are positioned several kilometres from the mast. Combination of multiple lidars to an ML measurement demonstrate that this also holds true for the derivation of the horizontal wind speed. ML wind speeds were considerably closer to the mast reference than the DBS lidar at this complex terrain site. The significantly reduced scatter and the smaller directional deviation in comparison with the reference sonic clearly demonstrate the advantages of focusing multiple lidar beams at one point over the classical DBS-profiling method. Correction methods (based on e.g., flow models) can be avoided, which can significantly reduce the uncertainty in applications like wind resource assessments or site calibrations. Exploiting the flexibility of the WindScanner technology, the measurements presented here could be extended to a scan of multiple measurement locations within a planned windfarm site. Especially in complex terrain this could significantly reduce the uncertainty associated with the horizontal extrapolation of the wind resource. For research applications like flow model validation, the ML technique removes the constraints of a mast measurement for an accurate measurement of the horizontal wind speed.
The ML measurements of v 2 r , u 2 and v 2 (for SW/SE/EE and SW/SE) showed a high correlation to the reference sonic. However, while larger scale turbulence was accurately captured, due to the relatively large measurement volume and the associated path averaging, attenuation for higher wave numbers could be observed. This results in an overall underestimation of the variances. An obvious way forward would be to reduce the measurement volume. With currently available lidar technology, however, this comes at the cost of a reduced measurement range and, therefore, reduces the flexibility of ML measurements. An alternative is to develop correction methods for the variance which is lost due to the line averaging. The good agreement between observed and theoretical spectral transfer functions is encouraging for this approach. First-order corrections could be based on model spectra as often done for eddy-covariance measurements (e.g., [42]). Alternatively, the correction could be based on an extrapolation of the spectra in the inertial sub-range [43,44]. This, however, requires that a part of the inertial sub-range is not affected by line averaging. More sophisticated approaches could exploit e.g., the spectral broadening of the received backscatter [25,45,46]. The results for the variance measurements of the DBS lidar clearly illustrate the errors introduced by cross-contamination in this technique and reiterate the unreliability of turbulence measurements from profiling lidars.
From a practical perspective, it is interesting that the good results of the ML method could also be maintained if only two lidars were used. This reduces the additional costs, which are one of the disadvantages of the ML method, when compared to the DBS lidar (or other profiling lidar). Differences between the different dual-lidar ML systems are likely to be mainly caused by the setup or, more precisely, the different angles between the WindScanners. While the mean wind speed is only slightly affected, the variances of the horizontal wind-speed components seem to be more sensitive to an unfavourable setup of the WindScanners. A detailed analysis of the effect of the setup on the accuracy of ML measurements is part of ongoing investigations.