LiSBOA: LiDAR Statistical Barnes Objective Analysis for optimal design of LiDAR scans and retrieval of wind statistics. Part II: Applications to synthetic and real LiDAR data of wind turbine wakes

The LiDAR Statistical Barnes Objective Analysis (LiSBOA), presented in Letizia et al., is a procedure for the optimal design of LiDAR scans and calculation over a Cartesian grid of the statistical moments of the velocity field. The LiSBOA is applied to LiDAR data collected in the wake of wind turbines to reconstruct mean and turbulence intensity of the wind velocity field. The proposed procedure is firstly tested for a numerical dataset obtained by means of the virtual LiDAR technique applied to the data obtained from a large eddy simulation (LES). The optimal sampling parameters for a scanning Doppler pulsed wind LiDAR are retrieved from the LiSBOA, then the estimated statistics are calculated showing a maximum error of about 4% for both the normalized mean velocity and the turbulence intensity. Subsequently, LiDAR data collected during a field campaign conducted at a wind farm in complex terrain are analyzed through the LiSBOA for two different configurations. In the first case, the wake velocity fields of four utility-scale turbines are reconstructed on a 3D grid, showing the capability of the LiSBOA to capture complex flow features, such as high-speed jet around the nacelle and the wake turbulent shear layers. For the second case, the statistics of the wakes generated by four interacting turbines are calculated over a 2D Cartesian grid and compared to the measurements provided by the nacelle-mounted anemometers. Maximum discrepancies as low as 3% for the normalized mean velocity and turbulence intensity endorse the application of the LiSBOA for LiDAR-based wind resource assessment and diagnostic surveys for wind farms.


Introduction
The use of Doppler light detection and ranging (LiDAR) technology for wind energy applications has largely increased over the last decade (Clifton et al. 2018;Veers et al. 2019). Thanks to the achieved measurement accuracy, simpler and cost-effective deployments compared to traditional met-tower instrumentation, this remote sensing technique is now included in the international standards as a reliable tool for performance diagnostic of wind turbines and wind resource assessment (International Electrotechnical Commission 61400-12-1 2017). Nonetheless, due to the limited spatio-temporal resolution and the distribution of the sample points in a spherical reference frame, the reconstruction of wind statistics from LiDAR samples still presents several challenges (Sathe et al. 2011;Newman et al. 2016).
In the companion paper (Letizia et al.), we presented a revisited Barnes objective analysis (Barnes 1964) for the calculation of wind statistics from scattered LiDAR data, which is referred to as LiDAR Statistical Barnes Objective Analysis (LiSBOA). This procedure enables the estimation over a Cartesian grid of the mean, variance and even higher-order central statistical moments of the radial velocity field probed by a scanning Doppler pulsed wind LiDAR. The LiSBOA performs also adequate filtering of small-scale variability in the mean field and mitigation of the dispersive stresses on the higher-order statistics provided that the algorithm is tuned based on the characteristics of the flow under investigation and the data collection strategy is optimally designed through the LiSBOA.
The LiSBOA capability to estimate statistics of an ergodic turbulent velocity field makes it a suitable tool for the analysis of wind turbine wakes and the resource assessment of sites characterized by heterogeneous wind conditions, such as in presence of flow distortions induced by a complex terrain. Over the last decade, wind LiDARs have been used to investigate wind turbine wakes; for instance, Käsler et al. (2010) and Clive et al. (2011) measured the velocity deficit past utilityscale wind turbines, while Bingöl et al. (2010) used a nacelle-mounted LiDAR to detect wake displacements and validate the dynamic wake meandering model (Larsen et al. 2008). Fitting of the wake velocity deficit was successfully exploited to extract quantitative information about wake evolution from LiDAR measurements (Aitken and Lundquist 2014; Wang and Barthelmie 2015;Kumer et al. 2015;Trujillo et al. 2016;Bodini et al. 2017).
A deeper understanding on the physics of turbine wakes was achieved by calculating temporal (Trujillo et al. 2011;Iungo et al. 2013b;Iungo and Porté-Agel 2014;Kumer et al. 2015;Machefaux et al. 2016;Van Dooren et al. 2016) or conditional (Aubrun et al. 2016;Machefaux et al. 2016;Garcia et al. 2017;Bromm et al. 2018;Iungo et al. 2018;Zhan et al. 2019Zhan et al. , 2020 statistics of the velocity collected through LiDAR scans performed at different times. Using this approach, Iungo and Porté-Agel (2014) detected a significant dependence of the wake recovery rate on atmospheric stability based on time-averaged volumetric LiDAR scans. The same concept was expanded by other authors using ensemble statistics (Machefaux et al. 2016;Carbajo Fuertes et al. 2018;Zhan et al. 2019Zhan et al. , 2020. Kumer et al. (2015) carried out a comparison between instantaneous, 10 minutes and daily-averaged velocity and turbulence intensity fields around utility-scale wind turbines, highlighting the presence of persistent turbulent wakes. Trujillo et al. (2011) used a nacelle-mounted LiDAR to quantify meandering-induced wake diffusion and added turbulence from statistics calculated over 10-minute periods.
Second-order statistics are of great interest in wind energy. Iungo et al. (2013b) used velocity time-series extracted from LiDAR fixed scans performed downstream of a 2-MW wind turbine to detect enhanced turbulence intensity in the proximity of the wake shear layers. More recently, temporal statistics over 30-minute periods allowed for the identification of turbulent wake shear layers from both numerical (Fuertes Carbajo and Porté-Agel 2018) and experimental (Carbajo Fuertes et al. 2018) velocity fields. Aubrun et al. (2016) attempted to characterize the turbulence intensity using bin statistics, even though achieving values higher than expected, i.e. larger than 50%. Zhan et al. (2019) used clustered data of wake velocity fields to retrieve a proxy for the standard deviation of wind speed in the wake of utility-scale turbines. These authors reported significant variability in the wake turbulent statistics depending on the atmospheric stability regime and operative conditions of the wind turbines.
In the light of great relevance for the wind energy applications of the statistical analysis of wind LiDAR data, for this work, the LiSBOA procedure is applied to virtual and real LiDAR measurements of wind turbine wakes. The scope of this study is dual: first, assessing the capabilities provided by the LiSBOA for the optimal design of the LiDAR scanning strategy by maximizing the statistical accuracy of the measurements and coverage of the sampling domain with the prescribed spatial resolution; second, showing the potential of the LiSBOA to reconstruct mean velocity and turbulence intensity fields from LiDAR data to unveil important flow features of wind turbine wakes.
With these aims, the LiSBOA is initially applied to virtual LiDAR data generated by a LiDAR simulator scanning the wake of a turbine modeled through the actuator disk approach in LES environment. This numerical test case allows for an extensive error analysis enabling the quantification of the LiSBOA accuracy. Then, real LiDAR data collected in the wakes generated by four 1.5-MW wind turbines are analyzed through the LiSBOA. Specific wake features, such as the high-speed jet around the nacelle and the turbulent shear layers, as well as perturbations induced by the complex topography, are detected. Finally, to provide a quantitative comparison with the data retrieved by means of traditional anemometers, the LiSBOA is employed to calculate mean velocity and turbulence intensity fields of the wakes generated by four 1-MW turbines interacting each other.
The remainder of the manuscript is organized as follows: section 2 reports the results of the virtual LiDAR data analysis, along with the description of the LiDAR simulator. Section 3.a provides a description of the site and the experimental setup of the field campaign. In section 3.b, the scan design and the reconstruction of the statistics of the non-interacting wakes are discussed, while section 3.c presents the results of the comparison between nacelle anemometer statistics and LiSBOA for the multiple interacting wakes. Finally, conclusions are drawn in section 4. The paper uses symbols introduced in the companion paper Letizia et al., which the reader is encouraged to review for a better understanding of the present manuscript.

LiSBOA validation against virtual LiDAR data
The LiSBOA algorithm is applied to a synthetic dataset generated through the virtual LiDAR technique to assess accuracy in the calculation of statistics for a turbulent flow. For this purpose, a simulator of a scanning Doppler pulsed wind LiDAR is implemented to extract the line-of-sight velocity from a numerical velocity field produced through high-fidelity large-eddy simulations (LES). Due to their simplicity and low computational costs, LiDAR simulators have been widely used for the assessment of post-processing algorithms of LiDAR data and scan design procedures Stawiarski et al. 2015;Lundquist et al. 2015;Mirocha et al. 2015).
A main limitation of LiDARs is represented by the spatio-temporal averaging of the velocity field, which is connected with the acquisition process. Three different types of smoothing mechanisms can occur during the LiDAR sampling: the first is the averaging along the laser beam direction within each range gate, which has been commonly modeled through the convolution of the actual velocity field with a weighting function within the measurement volume (Smalikho 1995;Frehlich 1997;Sathe et al. 2011). The second process is the time-averaging associated with the sampling period required to achieve a back-scattered signal with adequate intensity (O'Connor et al. 2010;Sathe et al. 2011), while the last one is the transverse averaging (azimuth-wise or elevation-wise) occurring in case of a scanning LiDAR operating in continuous mode (Stawiarski et al. 2013).
These filtering processes lead to significant underestimation of the turbulence intensity (Sathe et al. 2011), overestimation of integral length-scales (Stawiarski et al. 2015), and damping of energy spectra for increasing wavenumbers (Risan et al. 2018).
Three versions of a LiDAR simulator are implemented for this work: the simplest one is referred to as ideal LiDAR, which samples the LES velocity field at the experimental points through a nearest-neighbor interpolation. This method minimizes the turbulence damping while retaining the geometry of the scan and the projection of the velocity onto the laser beam direction. The second version of the LiDAR simulator reproduces a step-stare LiDAR, namely the LiDAR scans for the entire duration of the accumulation time at a fixed direction of the LiDAR laser beam. Two filtering processes take place for this configuration: beam-wise convolution and time averaging.
To model the beam-wise average, the retrieval process of the Doppler LiDAR is reproduced by means of a spatial convolution : where n is the LiDAR laser-beam direction. A triangular weighting function φ(s) was proposed by Mann et al. (2010): where ∆r is the gate length. The former expression is valid assuming matching time windowing, i.e. gate length equal to the pulse width, and the velocity value is retrieved based on the first momentum of the back-scattering spectrum. Despite its simplicity, Eq.
(2) has shown to estimate realistic turbulence attenuation due to the beam-wise averaging process of a pulsed Doppler wind LiDAR (Mann et al. 2009). Furthermore, a time-averaging occurs due to the accumulation time necessary for the LiDAR to acquire a velocity signal with sufficient intensity and, thus, signalto-noise-ratio. This process is modeled through a window average within the acquisition interval of each beam. For the sampling of the LES velocity field in space and time, a nearest-neighbor interpolation method is used.
The third version of the LiDAR simulator mimics a pulsed LiDAR operating in continuous mode and performing PPI scans, where, in addition to the beam-wise convolution and time-averaging, azimuth-wise averaging occurs due to the variation of the LiDAR azimuth angle of the scanning head during the scan. The latter is taken into account by adding to the time average an azimuthal averaging among all data points included within the following angular sector: where r is the radial distance from the emitter, while θ and β are the associated azimuth and elevation angles, respectively. The subscript p refers to the p-th LiDAR data point. Following the suggestions by Stawiarski et al. (2013), the out-of-plane thickness, ∆z, is considered equal to the length of the diagonal of a cell of the computational grid.
As a case study, we use the LES dataset of the flow past a single turbine with the same characteristics of the 5-MW NREL reference wind turbine (Jonkman et al. 2009). The rotor is three-bladed and has a diameter D = 126 m. The tip-speed-ratio of the turbine is set to its optimal value of 7.5.
A uniform incoming wind with freestream velocity of U ∞ = 10 m s −1 and turbulence intensity of 3.6% is considered. The rotor is simulated through an actuator disk with rotation, while the drag of the nacelle is taken into account using an immersed boundary method (Ciri et al. 2017). More details on the LES solver can be found in Santoni et al. (2015). The computational domain has dimensions (L x × L y × L z = 9D × 6D × 10D) in the streamwise, spanwise, and vertical directions, respectively, and it is discretized with 512 × 256 × 384 grid points, respectively. A radiative condition is imposed at the outlet, freeslip is enforced on the top and bottom while periodicity is applied in the spanwise direction. Ergodic velocity vector fields are available for a total time of T = 750 s. Fig. 1a shows a snapshot of the streamwise velocity field over the horizontal plane at hub height obtained from the LES. The respective data of the radial velocity obtained from the three versions of the LiDAR simulator by considering a scanning pulsed wind LiDAR deployed at the turbine location and at hub height highlight the increased spatial smoothing of the radial velocity field by adding the various averaging processes connected with the LiDAR measuring process, namely beam-wise, temporal and azimuthal averaging (Fig. 1).
The optimal design of the virtual LiDAR scan is based on the procedure outlined in Letizia et al. (see specifically the flowchart provided in Fig. 7). For the estimation of the flow characteristics, the azimuthally-averaged mean and standard deviation of streamwise velocity, as well as the integral time-scale are considered (Fig. 2). The use of cylindrical coordinates is justified by the axisymmetry of the statistics of the wake velocity field generated by a turbine operating in a uniform velocity field (Iungo et al. 2013a;Viola et al. 2014;Ashton et al. 2016).
The streamwise LES velocity field shows the presence of a higher-velocity jet surrounding the To reconstruct the mentioned flow features, the fundamental half-wavelengths in the spanwise and vertical directions selected for this application of the LiSBOA are ∆n 0,y = ∆n 0,z = 0.5D.
Furthermore, considering the streamwise elongation of the isocontours of the flow statistics shown in Fig. 2a, a conservative value of the fundamental half-wavelength in the x-direction ∆n 0,x = 2.5D is selected.
The relevance of the selected ∆n 0 is tested by evaluating the 3D energy spectrum of u/U ∞ and u 2 /U 2 ∞ in the physical and in the scaled reference frame (see Letizia et al.,Eq. (6)). The spectra are azimuthally averaged by exploiting the axisymmetry of the wake. The spectra in the physical features lying outside will be damped. Numerical integration of the 3D energy spectrum shows that 94% of the total spatial variance of the mean is contained within that sphere, which ensures that the energy-containing modes in the mean flow are adequately reconstructed with the selected parameters.
In For the optimization of the LiDAR scan, the LiDAR angular resolution, ∆θ, is evenly varied, for a total number of 7 cases, from 0.75 • , corresponding to a total number of scans L = 2, to 4 • , which corresponds to L = 41. The four values of σ recommended in Table 2 of Letizia et al. to achieve a response of the mean D m (∆ñ 0 ) = 0.95 are considered here. In Fig. 4, markers indicate the different σ and, thus, the response of high-order statistical moments, D 0 (∆ñ 0 ). The Pareto front shows that increasing ∆θ from 0.75 • up to 2.5 • drastically reduces the uncertainty on the mean ( I I ) by roughly 70 %, but does not affect significantly data loss consequent to the enforcement of the Petersen-Middleton constraint ( I ). For larger angular resolutions, the statistical significance improves just marginally, but at the cost of a relevant data loss. For ∆θ > 2.5 • , in particular, I becomes extremely sensitive to σ, with the most severe data loss occurring for small σ (i.e. small R max ). The Pareto front also shows that to achieve a higher response for the higher-order statistics, D 0 (∆ñ 0 ), generally entails an increased data loss and/or statistical uncertainty of the mean. This analysis suggests that the optimal LiDAR scan for the reconstruction of the mean velocity field should be performed with ∆θ = 2.5 • and σ = 1/4 or 1/6.  From a more technical standpoint, the error on the mean velocity field, u/U ∞ , appears to be relatively insensitive to the type of LiDAR scan, with the spatial and temporal filtering operated by the step-stare and continuous LiDAR being even beneficial in some cases. In contrast, the error on the turbulence intensity exhibits a more consistent and opposite trend, with the continuous LiDAR showing the most severe turbulence damping. This feature has been extensively documented in previous studies, see e.g. Sathe et al. (2011).
This error analysis confirms that the optimal configurations selected through the Pareto front (i.e. and data loss ( I = 33% and 37%, respectively).
The 3D fields of mean velocity and turbulence intensity for first optimal configuration, (i.e. ∆θ = 2.5 • , σ = 1/4, m = 5), are rendered in Figs. 6 and 7, respectively. Furthermore, in Fig. 8, azimuthally-averaged profiles at three downstream locations are also provided for a more insightful comparison. The mean velocity field is reconstructed fairly well regardless of the type of LiDAR scan, due to the careful choice of the fundamental half-wavelength, ∆n 0 , for this specific flow. On the other hand, the reconstructed turbulence intensity is highly affected by the LiDAR processing, which leads to visible damping of the velocity variance for the step-stare and even more for the continuous mode. The ideal LiDAR scan, whose acquisition is inherently devoid of any space-time averaging, allows retrieving the correct level of turbulence intensity for locations for x ≥ 4D, while in the near wake it struggles to recover the thin turbulent ring observed in the wake shear layer.
Indeed, such short-wavelength feature has a small response for the chosen settings of the LiSBOA, in particular ∆n 0 and σ (see Fig. 3b of the companion paper Letizia et al.). On the other hand, any attempt to increase the response of the higher-order moments, for instance by reducing the fundamental half-wavelengths or decreasing the smoothing and the number of iterations, would result in higher data loss and less experimental points per grid node.
Finally, Figs. 9 and 10 show u/U ∞ and u 2 /u over several cross-flow planes and for all the combinations of σ − m tested for the ideal LiDAR and the optimal angular resolution. For the mean velocity, the most noticeable effect is the increasingly severe data loss consequent to the reduction of σ, which indicates σ = 1/4 -m = 5 as the most effective setting. The turbulence intensity exhibits, in addition to the data loss, a moderate increase in the maximum value for smaller σ, which is due to the higher response of the higher-order statistics (see Table 2 in Letizia et al.).
However, this effect is negligible in the far wake were the radial diffusion of the initially sharp turbulent shear layer results in a shift of the energy content towards scales with larger ∆n, that are fairly recovered even for σ = 1/4.

a. Site description and experimental setup
LiDAR data collected during an experimental campaign carried out at an onshore wind farm are used to assess the potential of the LiSBOA algorithm for wind energy applications. The measurements were collected during a long-term experimental campaign conducted at a large wind farm located in North-East Colorado (Fig. 11). This wind park encompasses 221 Mitsubishi 1-MW and 53 General Electric 1.5-MW wind turbines. More technical specifications of the wind turbines are provided in Table 1 Table 2.
The atmospheric stability is characterized through the Obukhov length ( where ρ ref = 1.225 Kg m −3 is the reference density at the sea level, U SCADA is the 10-minute average of the wind speed measured by the nacelle-mounted anemometers, while the local air density ρ met is calculated from the meteorological data according to the international standard (International Electrotechnical Commission 61400-12-1 2017). Another important parameter derived form the SCADA data is the turbulence intensity at the rotor, defined as: where U SD,SCADA is the standard deviation of wind speed over 10-minute periods.
The two LiDARs performed a great variety of scans during the campaign, based on the specific phenomena under investigation. For the present analysis, we focus on the 3D reconstruction of non-interacting wakes using the high-resolution data collected with the Halo Streamline XR LiDAR and the 2D reconstruction of multiple overlapping wakes detected by the Windcube 200S.

b. Application of the LiSBOA to volumetric LiDAR data
The present section aims to explore the potential of the LiSBOA for the optimal design of a LiDAR experiment, data post-processing, and reconstruction of 3D flow statistics. The dataset used in the this section was collected on October 11, 2018 over the farm region shown in Fig. 11b through a StreamLine XR LiDAR. The goal of the experiment is to investigate the evolution of multiple turbine wakes advected over complex terrain. Fig. 14 shows the site of the deployment and the relative distances between the LiDAR and the turbine hubs. The deployment location was chosen to scan the wakes generated by the wind turbines B16-B19 for south-south-east wind directions.
The LiDAR was deployed off a county road that connects the plateau with the surrounding plains, with a consequent difference in altitude between the instrument and the base of the turbines of about 40 m. To probe the wake region of turbines B16-B19 (Fig. 12b)  For the selection of the optimal azimuthal angular resolution of the LiDAR scan, the LiSBOA is applied to produce a Pareto front for six possible angular resolutions, ∆θ, between 0.25 • and 4 • , and four values of the smoothing parameter, σ = [1/4, 1/6, 1/8, 1/17]. As shown in Fig. 15, the optimal LiDAR scan is that with angular resolution ∆θ = 1 • and σ = 1/4 or σ = 1/6. Generally, an increasing ∆θ entails a reduction of the standard deviation of the mean, I I , yet values higher than ∆θ = 1 • do not lead to significant reductions of I I while worsening the data loss, I , indicating a larger number of grid points not satisfying the Petersen-Middleton constraint.
In Fig. 15, the values of the cost function I and I I calculated from the LiDAR data after the quality-control process (Beck and Kühn 2017) are also reported for the optimal angular spacing of the LiDAR ∆θ = 1 • . It is noteworthy that there is negligible difference between the values calculated before and after the quality control of the LiDAR data, indicating that the data loss due The wind and power data, which are recorded by the SCADA (Fig. 18), confirms that the turbines experienced fairly homogeneous inflow conditions, with differences in power capture 5% smaller than the rated value. The values of normalized velocity together with the performance curves ( Fig. 13) indicate that the turbines were operating in region II of the power curve for the whole interval of interest.
Since statistical stationarity is an important assumption for the LiSBOA applications, adequate post-processing of the LiDAR data is needed to avoid effects on the reconstructed flow statistics due to the wind variability. Specifically, the wind speed variability is corrected by making the line-ofsight velocity non-dimensional with the wind speed measured from the met tower #1. Furthermore, scans performed when the wind direction was outside of the range θ w ± ∆θ w , with ∆θ w = 10 The effect of the combination σ − m on higher-order statistics is investigated by extracting the turbulence intensity at different cross-stream planes. The optimal pairs σ − m identified by the Pareto front analysis (Fig. 15), viz. σ = 1/4 -m = 5 and σ = 1/6 -m = 2, are tested here. One may expect that due to the difference in the response of the high-order moments of the fundamental mode between the two pairs, D 0 (∆ñ 0 ), the first case would exhibit a significantly lower u 2 /u with respect to second one. However, as shown in Fig. 21, the peaks of turbulence intensity are quite similar between the two cases. The main difference between the two reconstruction processes is a smoother distribution of u 2 /u for σ = 1/4 − m = 5. The similarity between the two cases is due essentially to two reasons: first, the smallest energy-containing length scales of the turbulence intensity field (i.e. shear layer thickness) are larger than the selected fundamental mode ∆n 0,y = ∆n 0,z = 0.5D; second, the larger number of points per grid node averaged for the σ = 1/4 case, leads to a higher variance due to the reduction of the bias of the estimator of the variance, which partially compensates the lower theoretical response. Summarizing, this sensitivity analysis suggests that the choice of the σ − m pair cannot be based purely on the theoretical response, since it does not take into account non-ideal effects deriving for the discrete and non-uniform data distribution. Instead, an a posteriori analysis of the statistics retrieved is recommended to select the best σ − m values.
Turbine-wake statistics are extremely sensitive to the width of the selected wind sector can be observed in the shape of the wakes among different wind turbines, with the wake of turbine B19, in particular, showing the velocity deficit and turbulence peak that are displaced above the hub height. Turbine B19 is also the only one facing a slightly inclined terrain (see Fig. 22), which may have caused a skewed inflow.

c. Application of the LiSBOA to interacting wind turbine wakes
An assessment of the accuracy of the LiSBOA in the calculation of mean wind speed and turbulence intensity is now provided for LiDAR measurements performed during the occurrence of wake interactions. To this aim, point-wise measurements provided by the nacelle-mounted anemometers and saved in the SCADA data of four closely spaced Mitsubishi wind turbines, roughly aligned with the wind direction, are compared with the statistics obtained from the postprocessing of the LiDAR data with the LiSBOA. The incoming wind is characterized by averaging measurements collected from all the anemometers and wind vanes installed on both met towers, which are located 12 km and 10.4 km away from the leading turbine F01 (Fig. 25). The Obukhov length is calculated from both sonic anemometers indicating a stable stratification regime. The SCADA data exhibits the typical signature of multiple wake interactions with reduced wind speed and power for downstream turbines, while turbulence intensity is enhanced, in particular for the F02 and F04 wind turbines.
The optimal design of the LiDAR scan is performed considering six values of ∆θ and four values of σ. The obtained Pareto front is shown in Fig. 26, which indicates ∆θ = 0.5 • and σ = 1/3, 1/4 or 1/6 as the optimal scanning parameters. The equivalent velocity retrieved by the LiDAR is made non-dimensional with the freestream velocity provided by the met-towers. The wind direction range is set to ∆θ w = 10 • , resulting in a total measuring period of 150 minutes. Data points lying above the top-tip or below the bottom-tip heights are excluded for this data analysis. The dynamic filter technique is used to reject corrupted LiDAR data, producing a total of 544,000 quality-controlled LiDAR samples over 1,327,000 collected LiDAR data within the selected wind-direction range.
The LiSBOA is carried out on a grid with resolution dx = 0.25∆n 0 , using the combination smoothing parameters -number of iterations σ = 1/6 -m = 1, which is, among the allowable combinations, the one providing the largest response of the higher-order moments. The obtained velocity and turbulence intensity fields over the horizontal plane at hub height are displayed in Fig. 27. The velocity deficit of F02 appears slightly larger than that detected behind the unwaked turbine F01, which is most probably due to the wake superimposition. An even deeper velocity deficit can be observed behind F03, which operates in a partially waked condition for this specific wind direction. Downstream of the third turbine, the wake deficit build-up saturates, confirming results from previous studies on close wake interactions (Barthelmie et al. 2010;Chamorro and Porté-Agel 2011). Finally, the relatively fast recovery of the wake of the trailing turbine, F04, can be ascribed to the enhanced mixing due to the wake-generated turbulence. Indeed, Fig. 27b shows significant wake-generated turbulence increasing past the leading turbine that reaching its maximum at a distance of 1D downstream of the rotor of F03. Interestingly, wake-generated turbulence is concentrated on the sides of the wake of F01, which experiences undisturbed flow, while it spreads among the whole wake region for the downstream turbines. This feature might be related to the presence of coherent wake vorticity structures in the near wake of turbine F01 (Iungo et al. 2013a;Viola et al. 2014;Ashton et al. 2016), while further downstream, the perturbed inflow promotes the breakdown of such coherent structures leading to more homogeneous turbulence. Finally, the large velocity deficit/high turbulence detected in the wake of F03 may be a consequence of the mentioned partial wake interaction, which exposes the rotor to a non-homogeneous flow resulting during off-design operations.
From a more quantitative standpoint, the incoming wind conditions experienced by each turbine are characterized to perform a direct comparison with the nacelle-anemometer data. To this aim, the mean velocity and turbulence intensity profiles are extracted from the LiDAR statistics at a distance of 1D upstream of the rotors over a segment spanning the whole rotor diameter. The sampling location is chosen based on previous studies (Politis et al. 2012;Hirth et al. 2015), since 1D is generally considered the minimum distance upstream of the rotor where the influence of the induction zone can be neglected for normal operative conditions. The averaged values of u/U ∞ and u 2 /u of each upstream profile are then used for the comparison with the respective values recorded through the SCADA.
A well-posed comparison of the wind statistics obtained from the LiSBOA, the SCADA and metdata requires two important elements: firstly, the statistical moments compared have to be equivalent; secondly, both the LiSBOA and the SCADA data must be representative of the freestream conditions experienced by each turbine.
Regarding the first issue, the mean field obtained through the LiSBOA, u, can be expressed as: where . T is the average calculated over the whole sampling period of 150 minutes, while . T is the 10-minute average performed by the SCADA and the met-tower acquisition system. U SCADA and U met are the 10-minute averaged velocities recorded from the SCADA and met-tower, respectively, while the symbol ∼ indicates statistical equivalence.
Similarly, for the comparison between the velocity variance calculated through the LiSBOA and the respective values recorded through the SCADA, we have the following relationship: where u andû are the velocity fluctuations with zero mean calculated over the time periods T and T, respectively. The parameter U 2 SD, SCADA is the velocity variance recorded by the SCADA over the periodT of 10 minutes.
To ensure that the SCADA mean and standard deviation of velocity are representative of the undisturbed wind conditions at each rotor, these velocity statistics are corrected for the flow distortion induced by the turbine through appropriate nacelle transfer functions (NTF), which converts the velocity statistics measured at the nacelle of a wind turbine to the corresponding freestream values measured from a met-tower located nearby. The IEC standard 61400-12-2 (International Electrotechnical Commission, 61400-12-2 2013) prescribes to calculate the NTF from the bin average with bin size 0.5 m s −1 of the velocity measured by a reference anemometer as a function of the nacelle wind speed. In the present work, besides correcting the mean wind speed as indicated by the IEC standards, a linear correction of the wind speed standard deviation is also applied, as suggested by Argyle et al. (2018). We adopted as reference anemometer that installed at 69 m above the ground on met tower #2. The SCADA data of Mitsubishi turbines H05 and H06, both falling in the range of distances from the met-tower recommended by the IEC 61400-12-1 (International Electrotechnical Commission 61400-12-1 2017), are used. Only the unwaked wind sectors calculated based on the same standard are considered. The described layout is shown in Fig.   28, while Fig. 29 shows the result of this analysis. There is a high correlation between the velocity measured by the met-tower and the nacelle-mounted anemometer (ρ = 0.976). Nevertheless, the NTF of the velocity reveals consistently lower values occurring at the nacelle compared to the met-tower, with a peak at 20 m s −1 . Concerning the standard deviation of velocity, the agreement between SCADA and met-tower data is significantly lower (ρ = 0.828), yet a linear correction can be still calculated with acceptable significance (error on slope and intercept are 0.0038 and 0.0034 with 95% confidence).
The results of the comparison between LiSBOA and SCADA are provided in Fig. 30. The mean velocity is accurately captured and confirms that F02 and F04 are the turbines mainly affected by the upstream wakes. The slightly higher momentum impinging F03 is mostly due to the imperfect alignment of that rotor with the upstream turbine wakes, which creates a condition of partial-wake interaction. A slightly larger discrepancy between LiSBOA and SCADA data is observed for the turbulence intensity, with a maximum difference of ∼ 3% for F03. Nonetheless, the main trend is well reproduced and the overall agreement is satisfactory. The observed difference in turbulence intensity can be related to several factors, such as turbulence damping due to the LiDAR measuring process and LiSBOA calculations, the accuracy of the NTF, estimate of the streamwise velocity from the LiDAR radial velocity or vertical dispersive stresses.
The effect of the sampling location upstream of the turbines in the LiSBOA field is investigated by quantifying the discrepancy of the LiSBOA statistics with respect to the reference SCADA values for all the turbines through the 95-th percentile of the absolute error, AE 95 . Fig. 31 shows AE 95 as a function of the distance upstream where the incoming flow is extracted from the LiSBOA statistics. For the mean velocity, it is confirmed that the value suggested by the literature (x = −1D) is sufficiently far from the rotor to limit the effects of the induction zone on the definition of the reference freestream velocity. Furthermore, the rotor thrust does not seem to have noticeable effects on the incoming turbulence, being the induction zone essentially devoid of significant turbulent fluctuations due to the loads of the turbine blades. The discrepancy between the turbulence intensity retrieved through LiSBOA and SCADA steeply increases for sampling locations further than 2D from the rotor.
Summarizing, the satisfactory agreement between LiSBOA and SCADA data achieved in the present study indicates the proposed procedure as a promising candidate for wind resource assessment, especially for complex terrains, and investigations of the intra-wind-farm flow.

Conclusions
The LiDAR Statistical Barnes Objective Analysis (LiSBOA) has been applied to three different cases of wind turbine wakes to estimate the optimal LiDAR scanning strategy and retrieve mean velocity and turbulence intensity fields.
The first dataset objectively analyzed has been generated through the virtual LiDAR technique, namely by numerically sampling the turbulent velocity field behind the rotor of a 5 MW turbine obtained from a large eddy simulation (LES). The 3D mean normalized streamwise velocity and turbulence intensity fields have shown a maximum error with respect to the LES statistics of about 4%. Wake features, such as the high-velocity stream in the hub region and the turbulent shear layer at the wake boundary, have been accurately reconstructed. This analysis has also confirmed that the optimal scanning strategy identified by the LiSBOA has been that producing the most accurate flow statistics.
Subsequently, the LiSBOA has been used to process real LiDAR data collected for a utility-scale wind farm. For the first test case, the statistics of the wakes of four non-interacting 1.5-MW turbines placed at the brink of an escarpment have been reconstructed. The optimal LiDAR scanning strategy has been selected through the LiSBOA, while the mean velocity and turbulence intensity fields retrieved through the LiSBOA have offered a detailed insight of the wake morphology.
Furthermore, a sensitivity analysis of the wind direction range has confirmed the robustness of the data selection and quality control methods.
Finally, the complex velocity field arising from the interaction of four 1-MW turbines has been analyzed by calculating first and second-order moments on the horizontal plane. The mean velocity and turbulence intensity extracted 1D upstream of the rotors have agreed well with the values provided by the nacelle anemometers, with maximum discrepancies as low as 3%.
The applications of the LiSBOA discussed in this work aims to showcase the potential of the Syst., XVIII (4), 529-542, doi:10.1017/cbo9780511921247.018.  F . 15. Pareto front for the design of the optimal LiDAR scan for the reconstruction of the wakes generated by the wind turbines B16-B19. The markers highlighted in red correspond to the respective parameters obtained from the actual LiDAR data after the quality control process.

LIST OF TABLES
F . 17. 3D LiDAR scans of five wind turbines: a) 10-minute average wind speed measured from the anemometers installed at 50-m and 80-m height on the met-tower #1; the error bar represents the standard deviation over 10 minutes; the shaded area represents the interval selected for the LiSBOA application; b) 10minute average wind direction in geophysical reference system measured from the vanes installed at 50 m and 80 m on met-tower #1; c), d) and e) equivalent velocity fields measured with PPI scans; the green arrow is oriented as the mean wind direction measured by the met-tower #1, while the black cross indicates the LiDAR location.
F . 19. 3D rendering of the normalized mean equivalent velocity field reconstructed with ∆θ w = 10 • . The three isosurfaces represent u/U ∞ = 0.45, 0.6 and 0.75, while the color maps represent cross-sections of the mean velocity field over the respective planes reported in the rendering. The dashed circles correspond to the rotor swept area of turbines B16-B19 (from left to right) projected onto the specific cross-plane.
F . 20. 3D rendering of the turbulence intensity field reconstructed with ∆θ w = 10 • . The two isosurfaces represent u 2 /u = levels of 20% and 30%, while the color maps represent cross-sections of the turbulence intensity field over the respective planes reported in the rendering. The dashed circles correspond to the rotor swept area of turbines B16-B19 (from left to right) projected onto the specific cross-plane.
F . 22. Rotor-averaged streamwise mean velocity and turbulence intensity as a function of the downstream distance from the turbine and associated altitude profile.
F . 25. Probability density functions of the met and SCADA data recorded from 21:00 to 1:00 MDT on the night between September 5 and 6 2018: a) wind speed from met-towers; b) wind direction from met-towers; c) inverse Obukhov length from our sonic anemometers; d) normalized wind speed from SCADA; e) turbulence intensity from SCADA; f) normalized power from SCADA. F . 28. Met-tower and turbines selected for the nacelle transfer function estimation. The directions highlighted in grey represent the valid wind sectors unaffected by turbine wake interactions. The dashed circles bound the allowed range of distances from the tower in compliance with IEC standard 61400-12-1 (International Electrotechnical Commission 61400-12-1 2017).