With 13.4 GW of new capacity installed in 2021 in the U.S., 1 wind energy is playing a major role in the transition toward the netzero emission goal by 2050. 2 However, the increased wind power harvesting capabilities are challenged by the scattering of available wind resources. To address this challenge, wind turbines are typically clustered into wind farms installed over sites characterized by variable terrain characteristics. As a consequence, compact wind farm layouts can give rise to highly complex flow features, such as wake interactions, [3][4][5] terrain-related effects, 6 internal boundary layers, 7 and farm-to-farm interactions. 8,9 A detailed quantification of the extension and magnitude of the velocity reduction of the incoming wind field while approaching a wind turbine rotor, which is referred to as induction zone, is of primary importance for a correct prediction of the power harvested by a turbine. The latter is, therefore, related to a location where the wind velocity is not affected by the pressure field induced by the turbine rotor. This task can be addressed only by quantifying the velocity deficit occurring by approaching the rotor as done, for instance, through an analytical model. 10 For wind farm operations, the incoming velocity field is also affected by the mutual interactions of multiple rotor-induced flow modifications, which are exacerbated as wind turbines are closely arranged in a wind farm layout. [11][12][13][14] Therefore, besides the singleturbine induction, a cumulative effect on the incoming wind field is generated for the entire wind plant, which is termed blockage. 15 Major experimental and numerical efforts have been performed to understand this wind farm phenomenon. For instance, two-way couplings between the first row of a wind farm and downstream turbines were identified by wind tunnel experiments with subsequent effects on the power produced by the upstream row due to blockage-related velocity reductions. 16 This phenomenon was empirically modeled as a function of the wind farm layout, in particular of the streamwise (s x ) and spanwise (s y ) spacing among turbines, as well as of the number of rows.
The global blockage effect on power production was estimated for the Lillgrund offshore wind farm by comparing wind speed distributions before and after the wind farm construction. 17 Annual energy production (AEP) losses due to wind farm blockage were estimated to be in the range between 2% and 4%. A Reynolds-averaged Navier-Stokes (RANS) model was used to interpret the experimental results collected by meteorological towers in the upstream region of three wind farms, before and after the construction of the power plants. 18 Specifically, a 3.4% velocity reduction was found 2D upstream of the first row (where D is the turbine diameter), while 1.9% velocity reduction was quantified between 7D and 10D upstream, associated with lower energy production from the first row with respect to the isolated-turbine case. Comparing the mean velocity field evolving around a single turbine, a single isolated row, and an eight-row wind farm with infinite width along the spanwise direction, it was found that the upstream row benefits both from close spanwise spacing (due to flow speedups caused by mass conservation) and from large streamwise spacing. 19 However, the authors evidenced that, for an isolated row of wind turbines, the velocity reduction within 2D upstream of the turbine row is larger than what was predicted with the analytical model proposed by Ref. 10. Furthermore, stronger velocity reduction was observed as the spanwise spacing between rotors reduces. By contrast, for a finite spanwise array of wind turbines, the increase in power production is different depending on the resolved domain size in the vertical direction and the presence of an inversion layer aloft the atmospheric boundary layer. 20 Finally, the impact of atmospheric stability on the induction zone and global blockage was investigated through numerical simulations as well, 20 with particular regards to stable conditions 19,21 and the possible presence of gravity waves. 22 Another important feature of a wind farm flow is the presence of speedup regions in the area between neighboring rotors with limited spacing. Due to mass conservation, a local increase in the freestream wind speed arises due to the confinement generated by the rotors and their wakes. The presence of these speedups has been observed both through numerical simulations 19,[23][24][25][26] and experiments, 11,12 and it has been shown to have a direct impact on the overall power production. [26][27][28] Similar results were obtained through large eddy simulations (LES) of wind farms with constant streamwise spacing but different spanwise spacing, and different spanwise offsets between consecutive wind turbine rows. 24 From the analysis of the time-averaged streamwise velocity, it was observed that speedups were enhanced for a decreasing spanwise distance. Furthermore, the overall power production was maximized for a layout where the wind turbines in the second row were installed within the observed speedup regions. Speedups were also identified at the first row of an onshore wind farm for varying wind directions leading to the generation of channel regions between consecutive wind turbines. 26 The occurrence of modifications of the incoming wind field induced by wind turbine rotors, such as single-turbine induction, global blockage, and speedups, calls for further assessments through real-scale flow measurements around utility-scale wind turbines. In this realm, light detection and ranging (LiDAR) anemometry offers compelling features to probe the flow both at the single turbine and at the wind farm level. In particular, the typical spatial resolution (from 10 m up to 100 m) and sampling rate (few Hertz) of modern LiDARs allows the design of several scans able to capture bi-and threedimensional features of wind farm flows. 13,[29][30][31][32][33][34] Focusing on the induction zone of an isolated turbine, three continuous-wave LiDARs were recently deployed to quantify a velocity reduction of 1%À3% at a distance of 1D upstream of an isolated onshore wind turbine, and up to 18% at the turbine position for wind turbine operating under high thrust and power coefficients. 35 Similarly, the induction zone of a single wind turbine was investigated both through profiling LiDAR and wind tunnel measurements. 36 After quantifying a velocity reduction of 8% at 0:64D upstream of the rotor, the authors found significant correlations between time series of power and upstream velocity measured 0:25D above hub height, thus suggesting that remote sensing data collected at upstream locations can efficiently be leveraged for proactive turbine control purpose. Nonetheless, numerical simulations highlighted the impact of terrain elevation changes upstream of a wind farm on the accuracy for the quantification of global blockage. 37 Regarding global blockage, scanning pulsed LiDARs were deployed to perform plane position indicator (PPI) scans at hub height reaching up to 40D upstream of an offshore wind farm under different thermal stratification and turbine operating conditions. 15 The strongest velocity reduction due to blockage (À4% between 25D and 5D upstream the array) was found to occur for stable conditions and operative conditions corresponding to the high thrust coefficient of the wind turbine rotors (C t >0:8), while no evidence of blockage was found for convective conditions and low-C t operations. The impact of atmospheric stability on wind farm blockage was further highlighted by data collected by five profiling LiDARs positioned in proximity of an onshore wind farm. 38 For the site under investigation for the present work, previous measurements performed with wind profiling LiDARs confirmed the primary role of atmospheric stability on the interaction between the incoming atmospheric boundary layer and a wind turbine rotor. 39 Specifically, under convective atmospheric conditions, velocity reductions of up to 3% were observed at an upstream distance of 1.5D while they were negligible at a further upstream distance of 3D. On the other hand, for more unusual flow conditions, such as during the occurrence of low-level jets and negative shear, the rotor-induced modifications consisted of speedups above the rotor top tip and velocity reductions in the lower half of the turbine rotor. 40 Although extensive efforts have been done to characterize global wind farm blockage, single-turbine induction zone and speedups for onshore sites, a detailed characterization of these phenomena for different operative and atmospheric conditions through real-scale field experiments is still lacking. In the present work, the experiment introduced in Ref. 39 is leveraged to quantify induction zone and speedup regions under a broader range of atmospheric conditions. In particular, mean velocity and turbulence intensity (TI) are quantified through scanning and profiling wind LiDAR measurements within a domain encompassing three onshore utility-scale wind turbines. The high spatial resolution adopted by the scanning LiDAR allows to quantify induction zone and speedup regions at the single-turbine level within the wind farm; such a resolution is obtained through the optimal design of LiDAR scans via LiDAR Statistical Barnes Objective Analysis (LiSBOA) technique. 4,41 The same algorithm is also used in post-processing to Journal of Renewable and Sustainable Energy ARTICLE pubs.aip.org/aip/rse retrieve the hub-height mean velocity field from the LiDAR data. The induction and speedup regions are systematically quantified for different wind turbine spacing associated with varying wind directions, as for previous laboratory experiments, 16,19 and for different regimes of ambient turbulence intensity. 15,38,39,42 The remainder of this paper is organized as follows: in Sec. II, the experimental campaign (Sec. II A), scanning LiDAR strategy (Sec. II B), data selection (Sec. II C), and quality control (Sec. II D) are presented. In Sec. III, the results are organized in profiling-vs scanning-LiDAR data assessment (Sec. III A), mean flow variability with atmospheric turbulence intensity (Sec. III B) and turbine spacing (Sec. III C). Finally, conclusions are reported in Sec. IV, while in the Appendix detailed analysis of the LiSBOA design of the volumetric scan is presented.
The experimental data presented in this study were collected before and after the construction of four onshore wind turbines characterized by a diameter D ¼ 127 m, hub height H ¼ 89 m, and rated power close to 3 MW. The aerial view of the experimental site is reported in Fig. 1(a). As appears from the topography, some changes in altitude are present but, overall, the site is classified as flat according to the IEC international standards. 43 In this work, two distinct right-handed reference frames are used: the first one, referred to as the "global reference frame," is aligned along west-east and south-north directions with the vertical coordinate increasing moving away from the ground; the second one, referred to as "hub reference frame" and indicated by ðx; y; zÞ is oriented along streamwise (with x increasing moving downstream), transverse (y), and vertical (z) directions, respectively, and it is obtained by aligning the x-axis with the downwind direction.
During the pre-construction phase (from June to August 2020), four ground-based Windcube v1 profiling LiDARs (PLs hereinafter) were deployed in the proximity of the future turbine locations to investigate the site climatology. [38][39][40] PL data are available at z ¼ 40, 64, 89 m (hub height), 114, 140, 153, 160, 180, and 200 m, where z is the height above the ground. The PLs are utilized to reconstruct the incoming wind at each height through a six-beam Doppler Beam Swinging (DBS) technique. 44 For a given height, the DBS technique adopts six beams with constant elevation angle (28 ) equally spaced between 0 -and 360 -azimuth angle to retrieve the 10-min-averaged radial wind speed at each azimuth. Subsequently, the magnitude of the horizontal wind component (U) and direction (H) are obtained by fitting a sinusoidal function onto the time-averaged radial wind speed across azimuth angles. This operation is repeated for all the heights to obtain vertical profiles of wind speed and direction [UðzÞ; HðzÞ, respectively].
In addition to U(z) and HðzÞ, the incoming wind is characterized by the shear exponent (a). The latter is obtained by fitting a power-law function over each wind speed profile obtained from the DBS scans: 13,30,45
where the subscript "R" refers to reference values assumed at the lowest available height. Each vertical profile of wind direction is further characterized through the wind direction value at hub-height (H hub ) and veer (dH=dz), both of them calibrated through a first-order polynomial fitting of the wind direction profile obtained from the DBS scans:
In the remainder of this study, large veer and shear coefficients will be used to rule out non-canonical boundary layer flows from further analysis (see Sec. II C for more details).
For the pre-construction phase, a prevailing hub-height wind direction from south-west (%240 ) is estimated from the reconstructed During the construction of the wind turbines (August 2020 to September 2020), all the remote sensing instruments deployed for this experiment were co-located for a period of 35 days to perform DBS scans with the aim of characterizing the instruments accuracy. The LiDAR data are quality-controlled and averaged over 10-min periods to then perform a linear-regression of data collected from the different LiDARs. The obtained cross correlation coefficient is always within the range ½0:998; 1:002 for all the measured heights with a deviation from unit slope smaller than 1.5%, and intercept smaller than 0.2 m/s. 39 All these statistical values are deemed suitable to investigate wind farm blockage and speedup zones.
Finally, each turbine was monitored through a Supervisory Control And Data Acquisition (SCADA) system collecting the 10min-averaged hub-height wind speed, power, turbulence intensity (TI), thrust coefficient (C t ), power coefficient (C p ), revolutions-perminute (RPM), and set point of each turbine. For the data postprocessing, SCADA data are interpolated over the PL04 time stamp to monitor inflow and turbine operating conditions in a synchronized fashion.
The SL is deployed to perform optimally designed volumetric scans over a region encompassing turbines T06, T07, and T08 [shaded area in Fig. 1(a)]. Each volumetric scan is made of a sequence of PPI scans probing the flow over a prescribed azimuth range with increasing the elevation angles (/). The angular limits chosen in this work are 10 h 100 for the azimuthal direction (assuming zero azimuth aligned with north) and / 75 for elevation angle. The SL gate length and sampling rate have been set equal to 50 m (0:39D) and 2 Hz, respectively.
The azimuth resolution (Dh) and the number of PPI scans (N PPI ) within the above-mentioned limits are selected by applying the LiDAR Statistical Barnes Objective Analysis (LiSBOA) algorithm. 4,41 The latter allows reconstructing time-averaged flow statistics within a threedimensional Cartesian volume with user-defined resolution. However, due to the finite resolution of the SL, statistics of the measured velocity field can be reconstructed only down to a fundamental halfwavelength vector [indicated by Dn ¼ ðDn x ; Dn y ; Dn z Þ in the hub reference frame], while spatial fluctuations with smaller halfwavelengths will be filtered out. To reconstruct meaningful flow features both in the induction zone and the speedup region at hub height, in this work the fundamental half-wavelengths are chosen as Dn x ¼ 1D; Dn y ¼ 0:5D, and Dn z ¼ 0:8D, while the grid resolution is set equal to 0:25Dn. A sensitivity analysis on the choice of the fundamental half-wavelength, which is reported in the Appendix, reveals a low sensitivity of the LiSBOA-reconstructed mean velocity field to the choice of the fundamental half-wavelength vector around the selected value.
Several combinations of Dh and N PPI have been tested to optimally reconstruct the chosen spatial wavelength within the domain ðx; y; zÞ ¼ ½À2; 4D  ½À4; 4D  ½À0:25; 0:25D centered around T07 at hub height (see the Appendix for further details). Eventually, the chosen volumetric scan parameters are Dh ¼ 2 and N PPI ¼ 5 corresponding to the following elevation angles / ¼ ½3:4 ; 4:4 ; 6:1 ; 10:1 ; 28 , and to a volumetric scan time of s s ¼ 129 s.
After collecting the data through the volumetric scan, the LiSBOA algorithm is leveraged to reconstruct streamwise mean velocity and TI on the same Cartesian grid chosen during the scan design phase. 4 To validate the mean flow reconstruction through the selected scan, the non-dimensional LiSBOA mean flow will be first assessed against PL data (Sec. III A). Subsequently, the non-dimensional LiSBOA mean velocity and TI will be evaluated for different inflow conditions (cf. with Secs. III B and III C) in terms of hub-height cross sections of the three-dimensional reconstructed flow.
In addition to the statistical uncertainty (addressed in Sec. II C), one source of error that must be accounted for when the induction zone is quantified via SL measurements is the instrumental error; 15 this uncertainty is introduced by fluctuations of the LiDAR's scanning head performing volumetric scans (resulting in uncertainty of the probed position at a certain radial distance), Earth curvature over the measured domain (resulting in uncertainty of the vertical position of the scanned points) and, if present, motions of the structure supporting the LiDAR. For the present experiment, the scanning LiDAR was fixed to the ground without any intermediate support. Furthermore, the error due to the Earth curvature over the measured radial distance (2050 m) is quantified as %0:66 m (thus negligible when compared to the turbine diameter and hub-height). Finally, the standard deviation of the azimuth and elevation angles with respect to the nominal positions is of 0:017 and 0:0012 , respectively, which correspond to a maximum deviation from the nominal position of 0.61 and 0.04 m, respectively, at a distance of 2050 m, thus negligible for the present purpose.
After collecting wind data from different instruments (PLs, SCADA, and SL) during the post-construction phase, the data are synchronized on a common time stamp. Specifically, PL and SCADA data are interpolated via the nearest-point algorithm onto the SL time stamp limited to 65 min. After interpolating all the data on the same time stamp, only region-II, non-derated operating conditions are down-selected. Thus, for each turbine only flow realizations characterized by hub-height wind speed (U hub , obtained from SCADA) within ½0:6; 1U rated (where U rated is the rated wind speed) have been considered in order to exclude both turbine de-rated conditions and operating conditions forced by the controller for low incoming wind speed. Finally, only conditions with power less than the rated value (P rated ) and RPM greater than the turbine set-point value are retained for further analyses. The latter ensures that curtailed turbine operations are ruled out from further analysis. The selected operative conditions are visualized as a function of the normalized hub height wind speed (U hub =U rated ) in Fig. 2 for further analysis. The significant reduction of available data is due to the highly occurring wind farm curtailment throughout the duration of the experiment.
The variability of the selected dataset is further explored considering the bi-dimensional normalized probability density functions (pdfs) reported in Fig. 3. In particular, the pdf of incoming wind speed at hub height normalized with the rated wind speed (U 1;hub =U rated , where U 1;hub is evaluated from PL04 data) vs veer [Fig. 3(a)] unveils significant occurrence of high-veer events at low wind speed (U 1 0:6 U rated ). It is noteworthy that such conditions mostly occur during nighttime and Journal of Renewable and Sustainable Energy ARTICLE pubs.aip.org/aip/rse are associated with the presence of low-level jets and negative-shear velocity profiles. 40 However, these events are characterized by low incoming velocity (U 1;hub < 0:6U rated ); thus, based on the abovementioned data filtering, they are not considered for further analysis.
In addition, wind veer and shear seem to be practically uncorrelated, as shown in Fig. 3(b).
The bi-dimensional pdf of hub-height wind speed vs hub-height wind direction [Fig. 3(c)] confirms the presence of two peaks highlighted by high occurrence values roughly located around 45 and 240 . Thus, following the results shown in Fig. 3(c), in the remainder of this study the selected data in region II will be binned based on the incoming wind direction at hub height within the range 200 H hub 270 with a bin size of 10 (as prescribed by Ref. 46). The bin sizes are reported in Fig. 3(c) with vertical red dashed lines. The bi-dimensional plot of TI (evaluated from SCADA data) vs hubheight wind speed [Fig. 3(d)] reveals events occurring all across the considered TI range. Thus, simultaneously with the wind direction binning, the selected data in region II will be subdivided also in terms of TI in the range 0% to 15% with a bin size of 5% [reported in Fig. 3 (d) with red dashed lines].
As a result of the aforementioned binning criteria, 21 nonoverlapping wind direction and TI conditions are identifiable. For each bin, the radial wind speed collected by the SL is injected into the quality control process described in Sec. II D. However, only bins including a minimum number of volumetric scans (N t ) should be considered for further analyses to limit the statistical uncertainty on the mean, which is a crucial parameter when the mean velocity in the induction zone must be estimated. To this aim, the statistical uncertainty on the mean for auto-correlated signals has been used 47 [Eq. (A2) in the Appendix]. Assuming a streamwise velocity standard deviation of r u ¼ 0:1U rated , a minimum hub-height freestream velocity of U 1;hub ¼ 0:6U rated , and an integral timescale of 24 s (see the Appendix), the non-dimensional statistical uncertainty on the mean (e II ) calculated from Eq. (A2) is estimated a priori as e II 0:056U 1;hub if N t ! 9 (cf. with the Appendix). Thus, only data bins featuring a minimum number of nine volumetric scans have been considered for further analysis. Based on this criterion, ten data bins (reported in Table I) are available for the remainder of this study. From Table I, homogeneous conditions across bins are observed in terms of C t , C p , incoming wind speed, and diameter-based Reynolds number (Re D ¼ U 1 D=, where ¼ 1:41 Â 10 À5 m 2 =s is the kinematic viscosity).
For each bin reported in Table I, the SL scans are examined to individuate possible outliers in the data, i.e., erroneous quantification of the radial wind speed (v r ) due to instrumental errors. To this aim, a dynamic filtering algorithm is used to remove instantaneous radial wind speed samples characterized by low occurrences. 48 Specifically, a spatiotemporal interval of ðDx; Dy; Dz; DtÞ ¼ ð1D; 1D; 0:5D; 112 sÞ in the global reference frame is defined around each SL physical point. For each point, all the instantaneous (v r , CNR) pairs (where CNR is the carrier-to-noise ratio) are utilized to estimate a normalized bi-variate probability density function centered about the space-time averaged values of v r and CNR. All the samples characterized by low occurrence ( < 0:1%) are treated as outliers and marked by Not-a-number (Nan) values. An example of the outcome of this step is visualized in Fig. 4(b) where, with respect to Fig. 4(a), velocity samples in the proximity of turbines T06, T07, and T08 are removed for this reason. Finally, the outliers located within the boundaries of the SL domain are replaced through the "inpaint_nans" function available in Matlab. 49 After assessing the quality of the radial wind speed, the nondimensional streamwise velocity (u=U 1 ) is calculated from each v r sample [Fig. 4(c
where ðh; /Þ represent the azimuth and elevation angles utilized by the SL to probe the generic position ðx; y; zÞ, t is the time stamp of the SL, and H ðz; tÞ is the incoming wind direction obtained fitting Eq.
( Bin ID 3), the dimensional (U) and non-dimensional (U=U 1 ) mean velocity fields are calculated through the LiSBOA algorithm as a time average of SL streamwise wind speed realizations. 4,42 Similarly, the standard deviation of the streamwise velocity fluctuations (r u ) and turbulence intensity based on shear-corrected non-dimensional SL data (TI SL ) are estimated as r u ðx; y; zÞ ¼ h u ðx; y; z; tÞ À U ðx; y; zÞ
i.e., obtained by subtracting the local time-averaged flow from the instantaneous streamwise wind speed evaluated over the SL physical points. The LiSBOA mean flow, originally calculated on the Cartesian points, is estimated over the SL points through linear interpolation and subtracted from the instantaneous values. 41 It is noteworthy that, for the data analysis detailed in Sec. III, r u from Eq. (4a) has been leveraged to estimate the statistical uncertainty on the mean (e II ) through Eq. (A2).
Before delving deeper into the analysis of the wind field over the induction zone and speedup regions for different wind conditions, it is important to assess the mean velocity field reconstructed through LiSBOA at hub height against PL01 and PL03 mean velocity data collected within the SL measurement domain. In particular, PL01 and PL03 are placed upstream of T07 and between T07 and T08, respectively, for the considered wind sector (200 À270 ); thus, they may be used as further quantification of the induction and speedup effects investigated in this study.
For each bin of quality-controlled data reported in Table I, the instantaneous non-dimensional velocity sampled by the SL is injected into the LiSBOA algorithm to reconstruct mean velocity and standard deviation through Eq. (4a) within the Cartesian domain limits ðDX; DY; DZÞ=D ¼ ½À2; 4  ½À4; 4  ½À0:25; 0:25; as stated in Sec. II B, a fundamental half-wavelength vector of Dn ¼ ½1 D; 0:5 D; 0:8 D is chosen. The resulting hub-height mean velocity is linearly interpolated over the PL01 and PL03 locations and finally compared against the PL mean velocity records at hub height (here named U PL01 and U PL03 , respectively). To obtain a meaningful comparison, the latter are averaged in time considering only PL flow realizations simultaneous to the SL time samples injected into the LiSBOA algorithm for each bin.
The linear regressions of the PL01 and PL03 data with the SL data are reported in Figs. 5(a) and 5(b), respectively, for all the bins listed in Table I. For both panels, the abscissa reports the hub-height mean velocity estimated by LiSBOA at the considered PL position, while the ordinate reports the reference mean velocity value measured Journal of Renewable and Sustainable Energy ARTICLE pubs.aip.org/aip/rse by the PLs. Overall, good agreement is found between reference and LiSBOA-estimated values of mean wind speed over the PL locations, with normalized mean absolute errors (NMAEs) of 5.8% and 2.0% of the incoming wind speed, thus lower than the estimated statistical uncertainty on the mean for the chosen scanning strategy (7.4%, cf. with the Appendix).
Although consistent mean velocity estimates are obtained from both regression analyses (R-square values of 0.75 and 0.98 for PL01 and PL03 references, respectively), a better agreement between SL and PL03 is observed with respect to SL vs PL01 data (slope of 0.957 vs 0.754, and intercept of 0.018 U rated vs 0.178 U rated ). The larger data scattering observed for the PL01 assessment is thought to be due to the enhanced flow three-dimensionality in the induction region, which is harder to reproduce both through SL and PL velocity measurements.
For flat terrain, the freestream TI at hub height is generally correlated with the atmospheric stability regime. 42 Thus, SL data from data bins 8, 9, and 10, which are characterized by the same wind direction and incoming TI within 0%À5%; 5%SSÀ10%, and 10%À15%, respectively (Table I), are injected into the LiSBOA algorithm to evaluate the hub-height non-dimensional mean velocity and TI SL (Fig. 6). The two-dimensional mean velocity fields [Figs. 6(a From Fig. 6, it is observed that with increasing TI wake recovery occurs over shorter downstream distances due to the enhanced turbulent fluxes [Figs. 6(d)-6(f)], as already documented in the literature. 13,30,42 Furthermore, for the speedup regions located on both sides of T07, the magnitude of the velocity decreases with increasing freestream TI [Figs. 6(a)-6(c)]. In the hub reference frame, the speedup region is characterized for different levels of the incoming TI both through the position of the maximum speedup (x max ; y max ) and through the centroid of the area where U ! U 1 . In particular, the streamwise (x C ) and spanwise (y C ) positions of the centroid are calculated as The mean velocity along the streamwise direction, x, is then extracted at two different spanwise locations situated along the T08 axis and within the T07-T08 speedup region, whose transects are reported in Figs. 6(a)-6(c) with continuous and dashed lines, respectively. The selection of T08 instead of T07 is due to the fact that, for the chosen bins, the SL-T08 direction is better aligned with the wind direction (260 À270 ) rather than the SL-T07 direction, thus SL measurements can provide a better accuracy in probing the streamwise velocity component. The mean velocity along the selected transects is reported in Fig. 7 in the hub reference frame along the x-direction centered at the turbine location (x T ). The shaded areas represent the statistical uncertainty on the mean calculated via Eq. (A2) with r u obtained from Eq. (4a). In particular, the mean wind speed within the streamwise range ½À2 D; 0 from T08 is depicted in Fig. 7(a), while the velocity deficit within 1D is further quantified in Table II.
Comparing the velocity reduction for different levels of the incoming TI, bin 8 with the lowest TI range shows the strongest velocity reduction as ðxÀx T Þ 0:5D (cf. also with Table II), in agreement with previous scanning wind LiDAR experiments. Journal of Renewable and Sustainable Energy ARTICLE pubs.aip.org/aip/rse
Both Table II and Fig. 7(a) reveal a good agreement between the considered data bins and the induction-zone model proposed by Ref. 10, which is expressed as
where a is the axial induction factor related to the thrust coefficient (C t ) by the classic momentum theory: C t ¼ 4að1ÀaÞ, and C t ¼ T= ð0:5qU 2 1 AÞ (T being the rotor thrust force). For each bin of the LiDAR data, a is estimated based on the mean C t value for turbine T08 conditionally averaged across the time realizations included in each bin. No sensible difference for the predictions obtained with the model of Eq. ( 6) is found across bins 8, 9, and 10 since the mean C t values are very similar, thus only the prediction for bin 10 is reported in Fig. 7(a) for the sake of clarity. However, the experimental mean velocity profiles show stronger velocity reduction as ðxÀx T Þ>À1D with respect to Ref. 10, in analogy with the results of Ref. 19 for similar spanwise turbine spacing (s y ¼ 1:26À3:93D vs s y ¼ 1:93D for the present bins). One possible explanation might be that the Ref. 10 model is based on the vortex sheet theory applied to an isolated wind turbine, while the current wind farm layout is characterized by four wind turbines closely spaced along the same row (s y ¼ 1:93D). Thus, the global wind farm blockage effect (not accounted for by the Ref. 10 model) could be responsible for the stronger mean velocity reduction observed from the experimental data. Another discrepancy between the Ref. 10 model and the current experimental scenario is the significant presence of shear in the wind speed profile, 20,21 which is not accounted for by the reference model.
The current findings about the speedup under different ambient TI conditions are slightly different than the results of Ref. 39, who found the strongest speedup (15%) under daytime, convective conditions for westerly wind. However, it should be noticed that in that study, TI data are binned over two larger ranges (TI < 10% and 10% TI 25% for night and daytime, respectively), which might lead to a more heterogeneous flow conditions in terms of incoming TI and, in turn, increase the uncertainty of the reconstructed mean flow. By contrast, in the current study TI is clustered into three bins (0 < TI 5%; 5% < TI 10%; 10% < TI 15%), thereby providing a more accurate data binning according to the background flow conditions.
Focusing on the wake region of T08 in Fig. 7(b), a faster wake recovery is observed for bin 10, as expected since it is characterized by the highest TI level, for ðxÀx T Þ ! 2D. In contrast, no TI-related effects are observed in the near-wake region (xÀx T < 2D). This confirms that the effects of the incoming TI on the downstream evolution of wind turbine wakes are mainly localized in the far-wake as a result of the enhanced turbulent fluxes and mixing. 30,42,50,51 The transect in the speedup region [Fig. 7(c)] reveals a significant dependence of the velocity profile on the ambient TI within xÀx T < 2D. Specifically, as TI increases, the velocity magnitude in the speedup region decreases. The speedup reduction is further quantified by the magnitude of the maximum velocity value across the speedup region [Fig. 8(a)] as a function of the ambient TI. As anticipated from Figs. 6 and7, the speedup effect observed in the T07-T08 speedup region decreases from 5:7% 6 1:4% of U 1 for bin 8 with the lowest TI level to 0:3% 6 0:7% of U 1 for bin 9 and becomes negative, i.e., no speedup, for bin 10 with the highest TI level (À1% 6 0:9% of U 1 ).
The centroid position of the speedup region is reported in Figs. 6(a TABLE II. Non-dimensional velocity extracted at streamwise locations between 0 and 1D upstream of the reference turbine (T08) for different freestream TI. Uncertainty intervals represent the statistical uncertainty on the mean calculated through Eq. (A2). insensitive to the ambient TI regime, as evidenced in Fig. 8(b) where they are referred to the T07 spanwise location and normalized with the T07-T08 spanwise distance [s y ¼ 1:96D, also reported in Fig. 6(f)]. By contrast, the streamwise positions with respect to T07 location [Fig. 8(c)] decrease transitioning toward higher freestream TI (x max ¼ 1:37 s x for the range TI ¼ 0À5%; 0:21 s x for the range TI ¼ 5À10%, and 0:25 s x for the range TI ¼ 10À15%).
The feature that the speedup magnitude is enhanced and detected further downstream from the rotor with reducing incoming-wind TI suggests that the occurrence of a speedup is not only associated with the pressure field induced by the turbine rotor, but also with the associated wake velocity field. In particular, with increasing incoming-wind TI, turbine wakes recover faster. 30,42,51 Therefore, the flow confinement induced by adjacent rotors and wakes is observed further downstream for reduced incoming-wind TI.
Given the wind farm layout reported in Fig. 1(a) and assuming that the turbine rotors are perpendicular to the incoming wind direction, any change in wind direction (H hub ) affects the streamwise (s x ) and spanwise (s y ) spacing among wind turbines in the hub reference frame, which are fundamental parameters driving the flow evolution across a wind farm. 16,19 It is noteworthy that the terms "streamwise spacing" and "spanwise spacing" here refer to the distances between adjacent wind turbines [cf. with Fig. 9(f)]. This definition is physically meaningful only when both s x and s y are strictly positive, and this nomenclature should not be mistaken with the analogous used for multiple-rows and multiple-columns wind farm configurations. 16,28 The following results will be presented both in terms of the incoming wind direction and the related horizontal area between turbines (s x s y ). For the wind farm under investigation, the spanwise spacing decreases as the wind direction shifts from south to west, while the streamwise spacing increases. Overall, the horizontal area between turbines (s x s y ) increases as the wind direction shifts from south to west.
For this analysis, bins 3, 5, 6, and 9 of Table I are selected as they are characterized by the same incoming-wind TI, wind directions within the sectors 230 À240 ; 240 À250 ; 250 À260 , and 260 À270 , respectively, and horizontal areas between turbines of s x s y =D 2 ¼ 0:78; 1:66; 2:46; 2:70, respectively. For each bin, non-dimensional mean velocity and TI SL around the turbines are depicted in Fig. 9. Focusing on the mean velocity [Figs. 9(a)-9(d)], a significant velocity decrease in the speedup region between T07 and T08 is observed increasing the horizontal area between turbines. Furthermore, the LiSBOA domain for bins 3 and 5 [related to Figs. 9(a) and 9(b), respectively] allows for comparing the speedup across different turbines. As a result, the most downstream turbine (T08 among those probed by the SL) experiences the strongest speedup, consistently with the numerical results of Ref. 26. Regarding the TI SL distribution [Figs. 9(e)-9(h)], no clear dependence is observed on the turbine spacing, which confirms that the TI evolution downstream the turbines is mainly governed by the ambient value of TI.
The mean flow within the speedup region between T07 and T08, as well as the wake center of the reference turbine, are further . Notably, for bins 3, 5, and 6 the direction SL-T07 is more aligned with the incoming wind, whereas for bin 9 the direction SL-T08 provides less misalignment; the reference turbine is selected accordingly. The so-obtained streamwise profiles of mean velocity are visualized in Fig. 10. The mean velocity distribution within the induction region is reported in Table III and Fig. 10(a). Good agreement is found between the reconstructed velocity profiles and the Ref. 10 model as ðxÀx T Þ À1D, while the latter seems to underestimate the velocity deficit closer to the rotor. As mentioned in Sec. III B, this might be due to a global blockage effect, which is not taken into account by Ref. 10. However, the statistical uncertainty associated with the reconstructed mean flow (reported in Table III) does not allow for further considerations about the effect of the spanwise spacing on the induction zone. Thus, we believe that further research is needed to assess the impact of turbine spacing on the single induction zone for a real wind farm layout. Observing Fig. 10(b), no dependence is found between the turbine spacing and the wake recovery associated with the reference turbine, while, as anticipated by Figs. 9(a)-9(d), a decreasing speedup between the turbines is associated with increasing values of s x s y =D 2 [Fig. 10(c)]. This effect is further quantified through the evaluation of the maximum speedup magnitude across different spacing and reported in Fig. 11(a). Here it is confirmed that an increased horizontal area between turbines (s x s y ) is associated with a reduced speedup at hub height. By contrast, the streamwise location of the centroid of the speedup region (x C ) reduces as the horizontal area between turbines increases [cf. with Fig. 11(c)], while the streamwise position of the maximum speedup is less sensitive to the area between turbines. Thus, it can be concluded that the speedup region reduces both in magnitude and extension as the area between turbines increases.
In this work, the flow developing at the inlet zone of an onshore wind farm has been investigated through one scanning wind LiDAR (SL) and three profiling LiDARs (indicated as PL01, PL03, and PL04) to quantify the mean streamwise velocity and hub-height turbulence intensity based on SL data, TI SL , within induction and speedup regions at hub height for different incoming wind conditions and ambient TI. Three turbines within the wind farm have been selected to position the profiling LiDARs within the main regions of interest both before and FIG. 10. Mean streamwise velocity at hub height for different incoming wind directions: (a) induction zone; (b) wake region; and (c) speedup region. Darker colors correspond to the wind shifting toward west (increasing the horizontal area between turbines). Uncertainty intervals are evaluated through Eq. (A2) with r u from Eq. (4a). after the construction of the wind turbines. In particular, PL04 was placed %10D upstream of the turbine row along the prevailing wind direction (240 from north), while PL01 was positioned within the induction zone (%1D upstream of the turbines row), and PL03 was located in the space between two consecutive turbines. The inflow data have been characterized from PL04, then classified into 10 nonoverlapping bins based on wind direction and TI of the incoming wind. The mean flow reconstructed through the LiDAR Statistical Barnes Objective Analysis (LiSBOA) algorithm from the scanning LiDAR data has been positively assessed against the profiling LiDAR data in terms of hub-height wind speed.
To investigate the effects of atmospheric stability on the turbine induction zone, the mean velocity for a selected wind sector and different incoming-wind TI has been investigated. First, it has been observed that the induction-zone model proposed by Ref. 10 systematically under-estimates the velocity deficit, which might be due to the mutual interaction of neighboring wind turbines and global blockage, which are disregarded by the considered model, or the effect of wind shear on the induction zone. Furthermore, wind conditions with incoming TI smaller than 5% have been found to provide the largest velocity reduction upstream of the turbine rotors, consistently with previous scanning wind LiDAR experiments. 15 The speedup region is characterized by a reduced intensity as the incoming-wind TI increases (from 5.7% of the asymptotic velocity for TI 5% conditions to be negligible for TI>10%). These effects have been associated with the enhanced turbulent mixing due to the increased TI and faster wake recovery.
The analysis of the hub-height mean velocity field for varying horizontal area between adjacent turbines, s x s y (simulated considering data bins for different wind sectors and the same ambient TI) has shown no remarkable effects neither on wake recovery nor on the induction zone upstream of the turbines. In contrast, turbine spacing has a significant effect on the speedup region extension and the velocity speedup. In particular, the speedup is found to decrease as the horizontal area between turbines increases. It should be noticed, however, that the present results hold for a single row made by four turbines. Results may be qualitatively or quantitatively different changing the number of turbines or adding rows in the downstream region.
The data that support the findings of this study are protected by a non-disclosure agreement.
In this study, the SL volumetric scan parameters are selected based on the LiDAR Statistical Barnes Objective Analysis (LiSBOA) algorithm. 4,41 As mentioned in Sec. II B, the statistical variability in space can be reconstructed only down to a fundamental halfwavelength vector Dn ¼ ðDn x ; Dn y ; Dn z Þ. Thus, assuming a certain half-wavelength vector dictated by physical flow features, it is desirable to reduce the portion of the domain under-resolved by the volumetric scan. According to the LiSBOA formulation, the percentage of the Cartesian domain where the selected half-wavelength vector is resolved is quantified by the Petersen-Middleton 52 where V is the sphere volume with radius 3r (r is the standard deviation of the Gaussian weights) centered on the generic Cartesian point (x, y, z) and n is the number of LiDAR points within the sphere. Notably, Dd ! 0 means perfect sampling, while Dd ! 1 implies under-sampling and, thus, local aliasing. The percentage of the Cartesian grid where Eq. (A1) is not satisfied (here named e I ) is the first cost function to minimize for an efficient volumetric scan design. The second cost function utilized for the optimal design of the volumetric scans is the statistical uncertainty on the mean (e II ). The latter is quantified for auto-correlated signals as 47 e
where s s is the sampling time of a single volumetric scan (129 s), s is the integral timescale of the flow, r u is the guessed standard deviation of the streamwise velocity (here assumed equal to 0:1U rated , where U rated is the turbine rated wind speed) and N t is the number of time realizations. In this work, the integral timescale is assumed as the ratio between a typical flow length scale and an advection velocity. 53,54 Assuming the latter as 5 m/s and an integral length scale of 120 m, 55 an integral timescale of s ¼ 24 s is therefore considered. It is understandable that, for a given spherical sector, a volumetric scan featuring low angular resolution (thus larger values of e I ) allows for a higher number of repetitions (N t ) within a certain time, which represents an efficient way to decrease e II . On the other hand, high-angular resolution scans are inevitably slow in time and, thus, associated with larger statistical uncertainty on the mean. Hence, a well-designed volumetric scan provides a trade-off between spatial and temporal resolutions, IV. For each combination of LiSBOA parameters, the non-dimensional LiSBOA mean velocity at hub-height is evaluated on the LiSBOA grid and assessed on the same grid against an independent Delaunay triangulation, which is a wellestablished statistical reconstruction technique for wind turbine flows. 30,41,[56][57][58][59] The normalized mean absolute error (NMAE) between the two reconstructed non-dimensional mean flows over the same Cartesian domain is assumed as error metric to quantify the effect of different LiSBOA fundamental half-wavelengths on the reconstructed flow. The SL data from bin 2 of Table I have been used for this analysis.
The NMAE is reported in Fig. 13 against the percentage of domain loss due to aliasing (e I ) calculated for each combination of Dn and r via Eq. (A1), assuming the SL azimuth resolution and number of PPI chosen for the experiment (Dh ¼ 2 and N PPI ¼ 5, respectively). The current choice of reconstruction parameters (Dn x ¼ 1D; Dn y ¼ 0:5D; Dn z ¼ 0:8D, and r ¼ jDnj=4, reported with red symbols in Fig. 13) leads to a NMAE of 0.094 with respect to the Delaunay algorithm. Notably, increasing r and Dn x [Fig. 13 The discrepancy between LiSBOA-and Delaunayreconstructed flows is further analyzed in terms of streamwise distributions of the difference between the two mean flows at Journal of Renewable and Sustainable Energy ARTICLE pubs.aip.org/aip/rse hub-height within the induction zone; the error is evaluated changing the streamwise half-wavelength within the interval reported in Table IV and keeping Dn y and Dn x constant (equal to 0.5 and 0.8, respectively), as well as r ¼ 1=4. The result is reported in Fig. 14.
Here, it is observed that, within the induction zone, the error between the two reconstruction algorithms is insensitive to the choice of streamwise half-wavelength and it is bounded between 610% up to x ¼ x T , then it increases past the rotor. This might be due to the larger spacing of the LiDAR physical points within the wake region leading to a coarser Delaunay triangulation, which ultimately causes larger discrepancies between the two reconstruction algorithms.
L. Zhan, S. Letizia, and G. Valerio Iungo, "Lidar measurements for an onshore wind farm: Wake variability for different incoming wind speeds and atmospheric stability regimes," Wind Energy 23, 501-527 (2020).
Journal of Renewable and Sustainable Energy ARTICLE pubs.aip.org/aip/rse