Tilted lidar profiling: Development and testing of a novel scanning strategy for inhomogeneous flows

The ability to reproduce realistic inflow conditions has been identified as the leading cause of errors for numerical wind turbine and farm models.^1,2 With this lesson in mind, a large portion of the suite of instruments deployed during the American WAKE experimeNt (AWAKEN)^3,4 has been devoted to the characterization (i.e., “profiling”) of the kinematic and thermodynamic properties of the atmospheric boundary layer when unaffected by wind turbines, or upwind of turbines. In this regard, the ground-based scanning pulsed Doppler lidars (simply lidars, hereafter) are the most versatile and powerful instruments among the AWAKEN fleet to estimate the inflow in terms of wind velocity field as a function of time and space. The reconstruction of the three-dimensional wind vector from line-of-sight (LOS) measurements detected by a single instrument is inherently ill-posed from a mathematical standpoint and thus requires strong assumptions to be carried out. The most commonly adopted lidar profiling techniques assume horizontal homogeneity of the velocity statistics of interest within the scanned volume. This assumption generally holds true over flat terrain or offshore in the absence of strong baroclinicity within the lidar scan volume. This assumption is also expected to apply at the AWAKEN site, which is characterized by relatively flat terrain and low vegetation, provided that the lidar scans are performed sufficiently far away from sources of heterogeneity, such as tall wind turbines. However, one of the inflow lidars at AWAKEN had to be deployed within a radius of 2.5 rotor diameters from an active turbine, where horizontal homogeneity of the wind flow may not apply (e.g., in the turbine wake and in the induction zone⁵) and the local wind field is not representative of the undisturbed inflow. This challenge prompted the investigation of innovative scanning techniques that could avoid the flow distortion caused by nearby turbines and still allow for a well-posed estimation of the undisturbed inflow for a wide range of wind directions, which is the focus of this work. In particular, we seek an unsupervised profiling technique capable of accurately reconstructing the mean flow and the Reynolds stresses (hereafter, velocity statistics).

The origin of current velocity profiling techniques dates back to radar research in the 1960s. The main idea is to retrieve the velocity components in the Cartesian reference frame from LOS wind velocity probed at different locations by knowing the scan geometry. Probert-Jones⁶ reported a method to estimate horizontal winds from a two-beam scan with low elevation and assuming horizontally homogeneous flow. Lhermitte⁷ proposed the vertical conical scan, named vertical-azimuth-display (VAD), that is used to estimate the horizontal and vertical wind component at each height from the magnitude of the LOS velocity as a function of the azimuth. Browning and Wexler⁸ expanded this technique by including linear wind gradients and perfecting the retrieval algorithm to use all the measurements. Fixed-beam techniques were still employed during that period (e.g., Balsey and Gauge,⁹ Strauch et al.¹⁰) but the error analysis by Koscienly et al.¹¹ proved that VAD is preferable since it exhibits the lowest theoretical error. The VAD became a standard in later years and was transferred with minimal adaptation to lidars,¹² extended to nonlinearly varying wind fields,¹³ and now represents the most commonly used wind reconstruction techniques for commercial profiling lidars.¹⁴ The uncertainty of VAD in horizontally homogeneous flow is comparable to that of traditional anemometers in measuring mean velocity,¹⁵ and profiling lidars based on VAD principles gained recognition among the wind industry as reliable instruments to perform power performance tests of wind turbines.⁵ 

Due to the importance of second-order point statistics (i.e., Reynolds stresses and other turbulence metrics) for meteorology and wind energy, many efforts have been made to estimate such quantities from VAD-like scans. Already, Lhermitte¹⁶ and Wilson¹⁷ proposed a method to decompose the variance of the radial velocity into single Reynolds stress components. Examples of the application of this theory were produced by Kropfli¹⁸ and Eberhard et al.¹⁹ using two different kernels to extract Reynolds stresses. Frisch²⁰ enhanced this technique by allowing Reynolds stresses to vary linearly in space. More recently, Mann et al.²¹ applied this approach to estimate vertical momentum flux, while Sathe et al.²² optimized a six-beam scan that yields Reynolds stresses through a simple matrix inversion rather than a decomposition. Frehlich et al.²³ and later on Smalikho and Banakh²⁴ also added an analytical correction for the effect of probe averaging and instrumental noise.

Another technique that has been successfully used to estimate second-order statistics includes either a paired^25–27 or a single, wind-aligned range-height-indicator (RHI) scan.^28–31 The method is also called z-bin since radial velocity data are analyzed in thin vertical layers where horizontal homogeneity is assumed. The difference is that the multiple-RHI method uses Reynolds stress decomposition (similar to Eberhard et al.¹⁹) and works regardless of the wind direction, while the single-RHI method requires good alignment with the wind direction (within $10°$– $20°$⁠) and is aimed at calculating the streamwise velocity mean and variance. Often, the vertical velocity is also neglected for low elevations. Indeed, when active alignment of the scan with the wind is possible and horizontal homogeneity holds true, even a single lidar stare can yield important turbulence information.^32–34 It is worthy of mention that more advanced and novel applications of RHIs that do not require the assumption of flow homogeneity have been lately proposed for the reconstruction of dissipation rate³⁵ and boundary layer depth³⁶ in complex terrain.

Finally, an attractive and more simple alternative to the former methods is the estimation or second-order statistics from a straightforward application of the eddy-covariance method to high-frequency (which means calculated once per scan cycle) fluctuating wind component.^37–40 

Several sources of error have been identified for lidars used to retrieve vertical profiles of wind speed in horizontally homogeneous flow. Most of the recipes proposed for error quantification of wind statistics boil down to the same principle: sampling error and instrumental noise are propagated on the wind retrieval as a function of the scan geometry and retrieval method. Furthermore, the LOS velocity departs from the local instantaneous value due to probe averaging and time accumulation, producing a bias of mean quantities in regions of strong vertical gradients,⁴¹ and an attenuation of second-order statistics.^39,42 For comparative error analysis using a reference sensor, the representativeness error due to the inherently different acquisition processes and sampling volume of the different instruments also adds up. Given the linearity between LOS velocities and the reconstructed wind vector, the error propagation part is straightforward, and several authors agree that low elevation angles are favorable for the estimation of the horizontal wind speed components, while higher elevations are preferable to measure the vertical velocity.^11,43,44 There is also widespread agreement that the error is fairly independent of the wind direction.^39,41,44 The treatment of the sampling error is instead more complex and deeply influenced by the atmospheric state. In the simplest models, the sampling error is proportional to the variance of the wind components and the reciprocal of the number of samples,^11,25,45 but more advanced approaches also take into account the integral spatiotemporal scales of turbulence.^19,41,42 Generally, larger absolute errors on the mean components occur in convective/high-turbulence conditions as a consequence of the larger sampling error.^41,45–48 Instrumental noise is instead commonly characterized as a function of signal-to-noise ratio (SNR) and quite small (⁠ $∼0.1$ m s⁻¹) if a robust quality control (QC) is implemented.⁴⁵ 

A special type of error affects the velocity variances obtained through direct calculation from profiling lidars. Some authors report a negative,^21,22,37 positive,⁴⁹ or case-dependent bias.^39,50 Newman et al.⁴⁹ elucidates this dilemma by providing evidence that the direct calculation of the horizontal velocity variances from high-frequency VAD-derived velocities leads to an overestimation due to the decorrelation between LOS velocity measured across the scanning area. In simpler terms, the estimation of Reynolds stress from direct calculation requires not only horizontal statistical homogeneity, but an even stronger instantaneous homogeneity, implying that all the beams probe simultaneously the same “frozen” turbulent eddy. The severity of this effect, named “cross-contamination,” is directly proportional to the distance between lidar beams at a given height and inversely proportional to the integral spatiotemporal scales of the flow. The inconsistent results on this topic can then be related to the strength of cross-contamination and instrumental noise (which can induce more variance) counterbalanced by the variance attenuation due to the probe/time averaging. Significant cancelation of the bias can, therefore, be expected,²⁴ which calls for caution in the interpretation of the results.

In addition to the mentioned sources of error, flow inhomogeneity within the scan volume results in additional uncertainty, being the wind reconstruction methods based on the assumption of statistically (or even instantaneously) homogeneous flow.⁵¹ First and foremost, an inhomegeneous wind field makes the definition of the “true” velocity profile quite elusive⁴¹ and certainly not representative of undisturbed conditions. Moreover, persistent spatial wind gradients can induce biases compared to the local wind field in the retrieval method itself.^11,52 Some authors suggested corrections based on numerical models in complex terrain.^52–55 Others successfully tested more pragmatic modified scanning strategies to mitigate the issue. The reduction of the coning angle for complex terrain scenarios has produced conflicting results.^50,52 Partial VAD or “arc scan” strategies can be used to restrict the sampling area to narrower (and thus more likely homogeneous) flow regions,^23,46,48,56 which, however, leads to a worse-conditioned retrieval than the full scan.⁴⁷ 

In this study, we address the problem of characterizing the wind profile in proximity of regions of horizontally inhomogeneous flow by adapting traditional conical scanning techniques. Our attention is focused on conical scans, due their low sensitivity to wind direction and widespread use in commercial lidars. In particular, the scan pattern is “tilted” to avoid the region where flow inhomogeneity is expected to occur. The new scans are tested by comparing their retrieval to the wind statistics obtained from a sonic anemometer. Furthermore, their application for the AWAKEN campaign is assessed through Weather Research and Forecasting Model (WRF) and large-eddy simulation (LES) numerical simulations incorporating a virtual lidar. A semiempirical model of the measured systematic and random error is also developed to shed light on the error mechanisms affecting the observations.

The paper is structured as follows: Sec. II outlines the mathematical framework adopted to design the tilted scans; Sec. III describes the field experiment; Sec. IV discusses the error metrics, data quality control and retrieval methods; Sec. V includes the results of the error analysis along with the semiempirical error model; Sec. VI quantifies the relevance of tilted profiling for AWAKEN through numerical simulations; and Sec. VII discusses the potential and limitations of the proposed scans and error model.

The modification of a conical profiling scan pattern to avoid a horizontally inhomogeneous flow region (IFR) can be formalized as a constrained optimization problem. For simplicity, we consider the reference system defined in Fig. 1 where a lidar is placed at the origin of the Cartesian reference system x = 0, y = 0, z = 0, where x and y in this case are generic horizontal axes, and z is the vertical axis. Without loss of generality, the IFR can be thought of as a volume bounded by $x<xmin, z<zmax$ and unbounded on the y axis. We can identify two possible scenarios: (1) if $xmin≤0$⁠, the lidar is not placed within the IFR, and a well-posed profiling is possible down to the lowest range gate [Fig. 1(a)]; (2) if $xmin>0$⁠, the lidar is within the IFR, and thus profiling is possible only above a certain $zmin$⁠, which is also a function of the scan geometry and that can be considered as an additional constraint for this more extreme case. The green and red areas in Fig. 1 indicate the regions where profiling is possible and ill-posed, respectively. An additional constraint on the minimum elevation $βmin$ (where β is the elevation, measured from horizontal) can also be enforced to limit the scan volume and stay clear of low obstacles.²² 

Two scanning strategies are redesigned within the framework defined: the six-beam scan²² and the Doppler-beam-swinging (DBS), a simplified version of the VAD used in commercial lidars (see Fig. 2). For each of these three scanning strategies, we create four different scans, based on the four sets of constraints in Table I. This process produces four scans, a regular one with the reference “untilted” geometry and used as a baseline, plus three with increasingly pronounced tilt, the latter being defined as the offset of the scan axis from the vertical direction which depends on the location of the IFR. The design of the tilted versions of these two families of traditional scans requires different approaches that are described in the following.

The main purpose of the six-beam scan is the reconstruction of the Reynolds stress tensors from the variance of the LOS velocities through a direct matrix inversion. The optimal scan pattern proposed by Sathe et al.²² is the result of the minimization of an objective function, F, representing the mean of the error variance of the reconstructed Reynolds stresses divided by the error variance of the individual LOS variance (assumed to be constant).

Tilted variations of the six-beam scan are obtained by adding the aforementioned constraints to the original optimization process. The minimization of F is carried out through a sequential least squares quadratic programming method with the constraints of Table I to identify optimal sets of azimuths (measured counterclockwise starting from the positive x direction) and elevations. The optimization is initialized 20 times with random azimuth and elevation angles to minimize the probability of the optimizer being attracted by local minima. A converged solution is always obtained within a few seconds. The optimal configurations are reported in Table II and also visualized in Fig. 2 (top row). Two elements substantiate the goodness of the optimization: first, the unconstrained scan corresponds to the regular six-beam scan, and second, the scans are symmetrical in the y direction, which is expected given the symmetry of the constraints along the same axis. The value of the objective function is $F=10.2,18.4,52,2299$ for the regular, slightly, moderately, and severely inclined scans, respectively. This variability suggests how severe the error amplification is expected to be for each scan.

The coordinates of the four tilted DBSs are given in Table III. All the scans are also illustrated in Fig. 2. It is noteworthy that a similar approach can be used to design any inclined VAD-like scan.

In order to assess the performance of the new scans, a Halo XR+ Streamline lidar has been deployed $∼70$ m southwest (typically upwind) of the base of the 135-m M5 meteorological (met) tower at the National Renewable Energy's (NREL's) Flatirons Campus near Boulder, Colorado, USA, from May to August 2022. Figure 3 shows a topographic map of the site, which is located approximately 5 km east of the Front Range of the Rocky Mountains. The experimental site is characterized by changes in altitude of less than 15 m. Local changes in the mean flow speed due to terrain should be less than 0.2 m s^–1 according to Wang et al.⁴⁸ Local orography can trigger unsteady flow features such as rolls,⁵⁷ although their presence is expected to similarly impact the velocity statistics of the lidar and the met tower. In light of these considerations, the flow probed by lidar and met tower is assumed to be horizontally statistically homogeneous, which allows for a well-posed comparison between the two instruments regardless of the scanning pattern.

Assuming ergodicity, turbulent statistics are calculated as time statistics using a 10-min averaging window. Also, following a widely accepted nomenclature, the instantaneous velocity components in the streamwise (x), lateral (y), and vertical (z) direction are indicated by u, v, and w, respectively. Capitalized U, V, and W indicate the instantaneous velocity components in the fixed frame of references west-east, south-north, and vertical (i.e., $W≡w$⁠) and will be useful in the derivation of the error model (Sec. V). Fluctuations on top of the mean are indicated with a prime symbol (e.g., $u′$⁠), and the time-averaged quantities are indicated with an overline (e.g., $u¯$⁠). Four quantities are selected for the comparison between lidar and met tower: the horizontal mean wind speed (in a vector average sense), $u¯$⁠; mean wind direction $θ¯w$ (which is $0°$ for northerly wind and is measured clockwise); the turbulent kinetic energy, $TKE=0.5(u′2¯+v′2¯+w′2¯)$⁠; and the turbulent shear stress, $u′w′¯$⁠. The selected variables provide a concise representation of the atmospheric state in terms of mean flow (⁠ $u¯$ and $θ¯w$⁠) and diagonal (TKE) and off diagonal (⁠ $u′w′¯$⁠) Reynolds stresses.

The lidar executed a sequence encompassing all the eight scan patterns shown in Fig. 2, tilted upstream toward the prevalent wind direction, i.e., 275° from north. Each scan is run for a duration of 10 min for a number of repetitions ranging from 21 for the slowest scan (the moderately inclined six-beam) to 64 for the fastest (the severely inclined DBS). The lidar operates in step-stare mode where each single beam represents the accumulation of 5000 individual shots, and it is collected once at each angular location. The gate length is 18 m with gate overlapping, which results in a different Doppler shift estimation every 1.5 m. The minimum range of the lidar is around 100 m; thus, colocated records from lidar and sonic anemometer are available only at 119 m above the ground. The anemometer is an ATI “K” type sonic anemometer and with a nameplate accuracy of 0.01 m s⁻¹ and an acquisition frequency of 20 Hz,⁵⁸ enabling it to characterize first- and second-order turbulence moments in the atmospheric boundary layer. To minimize the differences in the results due to imperfect synchronization between the different instruments, the target flow quantities are obtained by averaging high-frequency sonic data over the actual duration of each lidar scan. Overall, 6'079 10-min periods (about 42 days) are available for this study.

The details of the error metrics and the methods adopted to process raw met tower and lidar data are reviewed in this section, with special emphasis on the quality control of the sonic anemometer data, the retrieval of velocity statistics from the lidar, and the correction for tower distortion.

For the error analysis that follows, the four target turbulence quantities (⁠ $u¯, θ¯w$⁠, TKE, and $u′w′¯$⁠) retrieved from the lidar are compared to the sonic anemometer at 119 m in order to quantify the accuracy of the different scans and identify the different sources of error. In this framework, the sonic anemometer is taken as representative of the “true” wind characteristics. To this aim, the following error metrics are used:

• the mean error or bias, defined as:

  $Δf=⟨flidar−fsonic⟩,$

  (2)

  where f is the specific turbulence quantity from either lidar or sonic anemometer and the bracket indicate the average over a selected subset of observations (e.g., clusters of wind speed, direction, etc.). The bias is an informative metric to identify systematic discrepancies between the two instruments.

• The variance of the error (or alternatively, the standard deviation), defined as

  $Δ2f=⟨(flidar−fsonic−Δ)2⟩,$

  (3)

  which quantifies the time-varying or random errors. The form provided in Eq. (3) is particularly advantageous since it is consistent with the uncertainty that is estimated through linear error propagation;

• The linear regression coefficients (intercept and slope) between $flidar$ and $fsonic$⁠;

• The Pearson correlation coefficient, (4) x

  $ρ,$

  (4)

  between $flidar$ and $fsonic$⁠.

Sonic anemometer data are preprocessed through the data pipeline described in Clifton.⁵⁸ An inspection of the time series of the raw signals at 20 Hz revealed the presence of ubiquitous spikes (i.e., isolated points with nonphysical values) and ramps (sudden jumps in the signals), which required additional quality control (QC) to be removed. A detailed account of the QC applied to the present dataset is provided in  Appendix B. Overall, 96% of sonic velocity data and 90% of sonic temperature data at 119 m passed the QC.

Subsequently, tilt and roll in the anemometer are corrected through the planar fit method,⁵⁹ which identifies constant values of tilt and roll for the whole dataset. This rotation method is more suitable for a comparison with the lidar, as it allows for a nonzero mean vertical velocity, which is indeed possible in some instances, especially at this site. The algorithm yields a significant tilt of $3.8°$⁠, which has been attributed to an imperfect leveling of the instrument, in that the median slope of the terrain in the zoomed area shown in Fig. 3 is too small (⁠ $∼1.2°$⁠) to explain such flow inclination.

Finally, the 10 min temporal averaging and the eddy-covariance method are used to calculate mean velocity vector and Reynolds stresses from the quality-controlled data. Only 10-min periods with data availability higher than 80% are considered to ensure appropriate data coverage. The use of a 10-min period for Reynolds decomposition is standard in the wind energy industry,⁶⁰ although longer periods are typically used in boundary-layer meteorology.⁶¹ Throughout the work, except for Sec. IV D, the flow statistics sonic anemometer are taken as the benchmark for the quantification of the error of the lidar. Although this is a common practice, we acknowledge that, even in the case of perfect horizontal homogeneity, a zero error is hardly attainable due to the different sampling error and the larger vertical sampling volume of the lidar compared to the sonic anemometers.

The velocity retrieval applied to the lidar data transforms LOS velocity into corresponding mean velocity vector and Reynolds stresses. Before proceeding with the flow reconstruction, outliers in the raw lidar data are flagged through the dynamic lidar filter.⁶² The random noise affecting the LOS measurements is then characterized as a function of SNR by applying the autocorrelation function method⁶³ to a set of 33 h of vertical stares. The noise in the quality-controlled observations is estimated by interpolating the obtained SNR vs noise curve on the SNR detected during the experiment at 119 m. The noise has a median and 95th percentile of 0.12 and 0.27 m s⁻¹, respectively. In comparison, the characteristic standard deviation of the wind components from the sonic,⁴⁸  $23TKE$⁠, is greater than the noise 98.7% of the time. This indicates that the contribution of noise to the second-order moments is negligible, thanks to the close range of the measurements used in this work. The results in Sec. V will confirm that the error is dominated by the sampling error with no detectable contribution from instrumental noise. Therefore, noise will not be considered for the present study, although we acknowledge is may play a more important role for observations collected at further ranges.

The components of the mean velocity vectors are obtained for all the scans by simply averaging over the scan duration the “instantaneous” (i.e., collected over a single scan sweep) velocity estimated field, U. An important property ensues from the linearity of Eq. (6): the mean velocity retrieved from the lidar, regardless of the scan, only requires the assumption of statistical horizontal homogeneity. This requirement is significantly more lax than the assumption of instantaneous homogeneity that one makes when interpreting U as the estimator for the instantaneous velocity field. The knowledge of the mean flow is sufficient to estimate the wind direction and rotate the mean velocity vector, $U¯$⁠, into the final mean velocity vector aligned with the wind, $u¯$⁠. It is worth noting that Eq. (6) assigns equal weights to all the measurements that in this case were affected by generally negligible noise. To reconstruct winds at further ranges, where the SNR is expected to be lower and more heterogeneous across the scanning volume, it may be better to use an inversion method that weights in the uncertainty of the LOS data due to noise.⁴⁵ This is especially true for tilted profiling scans, where the range (and thus the SNR) corresponding to a certain height above the ground may differ significantly from one beam to the other, especially in the far range.

Unlike the calculation of the mean flow, the estimation of Reynolds stresses is carried out differently for each scan type depending on the number of available independent beams. For the DBS scans, which have less than six beams, the eddy-covariance method is applied directly to the instantaneous estimation of the velocity field, U, thus assuming instantaneous homogeneity.

After the Reynolds stress tensor in the fixed frame of reference is known, TKE is readily obtained as half of its trace (which is invariant to rotation), whereas the streamwise shear stress $u′w′¯$ requires a projection and rotation of the stress tensor into the wind-aligned frame of reference.²² 

A careful data filtering and correction is carried out to mitigate the influence of several factors that can jeopardize the comparison between the lidar and sonic anemometer, especially the flow distortion due to the met tower. Preliminarily, data collected during precipitation events as recorded by the sensors installed at the met tower base are discarded (8.1% of cases). Furthermore, instances when either the met tower or the lidar are within the wake of one the three research wind turbines active at the site are excluded using the method formalized in the IEC standards,^5,  Appendix A (2.7% of cases). Finally, a more elaborated procedure is implemented to minimize the effect of the distortion induced by the met tower and the supporting structure on the flow experienced by the anemometer. Significant flow distortion of a twin met tower at the site was already identified by Clifton et al.⁶⁴ and tailored CFD simulations (unpublished work) but not comprehensively documented. For this work, we adopt the method reported in McCaffrey et al.⁶⁵ who used conventional profiling lidar data from a lidar placed $∼40$ tower widths away from the tower base as representative of the undisturbed inflow. In our study, more than 98% of the center points of the lidar beams at 119 m across all the scans are located more than $∼40$ tower widths away from the met tower, and we can therefore consider the overall collection of lidar measurements negligibly affected by flow distortion due to the tower. Before proceeding to the analysis, it is necessary to ensure that the innovative scanning geometry did not introduce any significant error. In this regard, we inspect the results of a preliminary error analysis between the lidar and sonic anemometer prior to any corrections and for all wind directions (Fig. 4, red dots). The high correlation (⁠ $ρ∼0.99$⁠) and good error metrics, in general, confirm that the lidar is a trustworthy source to estimate the mean flow, regardless of the scanning geometry. The whole lidar dataset can then be used as the reference to quantify the flow distortion due to the met tower. Just for clarity, this section will be the only instance where the lidar, instead of the met tower, is used as ground truth.

The flow distortion is quantified by evaluating the percentage difference of the mean wind speed and direction between the sonic anemometer and lidar for non-overlapping $15°$-wide bins in the wind direction. Normalizing the difference in velocity by the incoming wind speed is customary for analyzing flow distortion around met masts^5,65 and it is based on the assumption that, after such normalization, a high Reynolds number flow should be relatively insensitive to the undisturbed wind speed (Reynolds similarity⁶⁶). Bin-averaged percentage flow distortions as a function of wind speed and direction (not reported here) showed that Reynolds similarity holds true for the present dataset, except for values of wind speed below 2 m s⁻¹, that have therefore been excluded in this section. The statistical convergence is checked by evaluating 95% confidence interval using the bootstrap method.^67, Figure 5 shows the outcome of such a procedure, both in the form of a linear plot [Fig. 5(a)] and a vectorial representation [Fig. 5(b)]. The most evident feature is a 15% velocity deficit in correspondence to the expected tower wake (⁠ $θ¯w=98°$⁠). Furthermore, a more subtle but statistically significant 5% slowdown occurs for wind approaching the tower from the anemometer's side, a probable signature of induction. The vectorial plot provides evidence that a significant flow deflection also occurs for winds blowing perpendicular to the mounting boom. For those wind directions, the anemometer likely senses the streamline expansion caused by the tower. All the foregoing effects appear to have a clear physical explanation and are also reported by McCaffrey et al.,⁶⁵ albeit for a different tower. In the remainder of the paper, the sonic anemometer data from the 60°-wide sector around the tower expected wake location (shaded area in Fig. 5 as suggested by Clifton and Lundquist⁶⁸) are excluded (9.4% of occurrences). The rest of the measurements are corrected for tower blockage and deflection by subtracting the differences shown in Fig. 5. The effect of the correction on the global error can be seen in Fig. 4 (blue data and trend line). No tower correction is applied to second-order statistics outside of the wake, since the effect of tower induction on turbulent statistics was not detected.

After describing the overall methodology and QC of sonic and lidar data, we present the results of the comparison between these two instruments, preceded by an overview of the atmospheric conditions encountered during the field test.

In general, the atmospheric conditions during the experiment were compatible with historical observations conducted in the same season.⁶⁹ In general, the climatology at the experimental site is heavily affected by the presence of the Rocky Mountains to the west that, in conjunction with global and mesoscale weather patterns, creates a markedly directional behavior of the wind. Clifton and Lundquist⁶⁸ identified four different wind regimes, characterized respectively by northerly, southerly, westerly weak, and westerly strong winds. The mean flow during the three months of the lidar deployment characterized based on the data collected by the sonic anemometer at 119 m (Fig. 6). The wind rose [Fig. 6(a)] reveals equally probable northerly, westerly, and southerly winds, with more sustained wind events from the west and west-northwest. Diurnal cycles of TKE and turbulent shear stress and their 95% bootstrap confidence intervals are shown in Figs. 6(b) and 6(c) separately for three wind sectors, similar to those defined by Clifton and Lundquist⁶⁸ with the westerly weak and westerly strong clusters merged. These flow quantities exhibit the expected daily patterns, with low TKE and (negative) turbulent shear stress indicative of stable conditions during the nighttime (∼01:00–12:00 UTC) and enhanced TKE and shear stress (unstable conditions) during daytime. Wind direction plays a role as well. Specifically, the strongest daytime winds (from the west) are associated with higher TKE and shear stress, all suggesting the occurrence of significant mechanical generation of turbulence. There is an unexpected reversed sign of the shear stress for northerly and southerly winds, which is not fully understood. Although this may indicate non-canonical flow conditions, their relatively large statistical uncertainty over the limited dataset prevents further analysis.

In this first part of the actual comparison between lidar and sonic anemometer, the overall agreement between the wind statistics estimated by the two instruments is first quantified by evaluating the error metrics for the eight different scans separately and for the whole dataset. Figures 7 and 8 show the linear regression and error metrics for mean wind speed and direction, respectively. The level of agreement is excellent and in line with the literature reviewed in Sec. I for regular scans¹⁵ and slightly better than what reported for arc scans operating in a similar way (see Figs. 6 and 13 in Wang et al.⁴⁸). Most importantly, the error is weakly dependent on the scan type and tilt, which points to the important conclusion that the mean flow can be accurately estimated regardless of the scan geometry. A systematic negative bias in wind speed, $Δu¯<0$⁠, seems to occur for all cases, but its magnitude is so small that it could not be investigated further, as it overlaps with the sonic anemometer's instrumental accuracy. Also, some sparse outliers appear in the wind direction estimate, although further analysis (Sec. V D) showed that they occurred during calm wind events when the wind direction itself is mathematically ill-conditioned and not of interest.

While the mean quantities are robust to changes in scan geometry, second-order statistics are not (Figs. 9 and 10). The main takeaways are

• TKE and $u′w′¯$ are fairly predicted with regular scans, but the random error increases rapidly when the scans are tilted; in some cases, TKE can even become negative for the six-beam, a pitfall of these techniques that was previously reported.⁴⁵ 

• The error metrics indicate better agreement for TKE than $u′w′¯$⁠. We note that Mann et al.²¹ shows a better correlation for their comparison of lidar-based $u′w′¯$ with sonic anemometer; however, their statistics were likely less affected by statistical uncertainties since they calculated 30-min statistics while we use 10-min windows,⁶³ which is the industrial standard for wind energy.

• Inclined DBS scans systematically overpredict the TKE (and magnitude of shear stress, although not clear from this global analysis) as indicated by the large slope and bias, a behavior that is investigated in Sec. V C.

To conclude, this global error analysis indicates that the six-beam scans up to moderate tilt provide an excellent reconstruction of the mean flow and are also devoid of significant biases on second-order statistics. We, therefore, conclude that the six-beam scan with slight or moderate tilt is the most suitable for tilted profiling. This early evidence will be further substantiated in Secs. V C and V D, where we investigate the causes of the trends just observed.

In Sec. IV C, it was shown how the Reynolds stresses from the DBS scans are obtained by direct application of the eddy-covariance method to the least-square estimate of velocity component by assuming horizontal homogeneity of the instantaneous flow. We also noted in Sec. V B how second-order statistics are systematically overestimated from the DBS. Similar biases were already observed in previous studies^39,49 and attributed to the breakdown of the assumption of instantaneous homogeneity, a phenomenon which is commonly referred to as cross-contamination.

To validate the proposed framework, the Reynolds stresses for both the DBS hypothesis and wide scan approximation are estimated by replacing the $U′nU′p¯$ term in Eq. (16) with the Reynolds stresses obtained from the sonic measurements every 10 min. Since bias is defined as the mean error, it is crucial to choose the appropriate averaging dimensions to properly characterize its behavior. In fact, calculating the bias for each scan but for all wind conditions, as done in Figs. 9 and 10, may mask the dependence of bias from specific quantities. Ideally, we want to characterize the bias as a function of the lowest number of variables that significantly affect its magnitude. To single out the environmental variables that are descriptive of the bias, we perform a permutation feature importance ranking based on a random forest model. The algorithm takes as inputs six variables (⁠ $u¯, θ¯w$⁠, TKE, $u′w′¯$⁠, Ri, and scan tilt) as candidate descriptors of the outputs. We consider two separate models, the first having as target the signed error on TKE, and the second the signed error on shear stress, between lidar and met tower for the DBS scans only. We tune five hyperparameters (the number of trees in the forest, the maximum depth of the trees in the forest, the maximum number of features considered when looking for the best split, the minimum number of samples for a split, and the minimum number of samples in a leaf) with a randomized search under a fivefold cross-validation, where the ranges of the hyperparameters are chosen such that the models do not overfit the training data. Once the hyperparameters have been tuned, the two final random forests (one for each target variable) are fitted on all the available data, and the permutation importance is assessed on them. The results are provided in Fig. 12 and clearly indicate that the signed error (and thus the bias) is driven by TKE and tilt for the TKE itself, and by wind direction and tilt for the shear stress. In light of these findings, the biases are calculated separately for each scan tilt and in 10 equally populated bins of TKE for TKE and 20°-wide wind sectors for shear stress.

Figure 13 shows the polar plot of the obtained bias of TKE and shear stress for the DBS scans from experiment (continuous lines) and the proposed model (dashed lines). The agreement is fair, thus suggesting that the observed bias in the measurements is mainly due to the failure of the instantaneous homogeneity hypothesis and that the wide scan approximation applies at the selected height of 119 m. Both observation and model capture the positive bias of TKE that is exacerbated by the tilt and high TKE. The magnitude of the bias of shear stress also increases with tilt, but is strongly directional with a positive sign for westerly flow (where the scan is tilted) and negative for easterly flow. Interestingly, $u′w′¯$ is unbiased only for winds perpendicular to the scan tilt direction, where the beams along the flow directions are symmetric. A deeper inspection of Eq. (16) revealed that the estimation of Reynolds stresses from DBS benefits from symmetry properties of the scan geometry, which result in a more diagonal-dominant system that minimizes cross-contamination.

In this section, we showed how Eq. (16), despite its simplicity, can provide an accurate estimate of the bias of second-order statistics due to cross-contamination in a profiling scan. A weakness of this framework is that it requires the knowledge of the true Reynolds stresses, which is the reason why it was used diagnostically in this section. A prognostic use of Eq. (16) would require some sort of prediction of $Un′Up′¯$⁠. Alternatively, the bias could be calculated parametrically for different turbulent states (e.g., turbulence intensity, shear, etc.).

Unlike the bias, the random error, quantified through the error standard deviation, is ubiquitous and affects all the scans. The feature selection based on the random forest is used again to identify the main drivers of the random error. The process is analogous to the one described in Sec. V C, with the only difference that the output variable is the absolute value of the error after the bias of the DBS quantified in Sec. V C is removed. This preconditioning of the error ensures that (1) error cancelation does not occur and (2) the effect of bias is excluded. The resulting importance scores are shown in Fig. 14 for all the wind statistics. For the wind speed absolute error [Fig. 14(a)], the TKE and scan tilt are slightly more important than the other variables and are then selected as descriptors for the random error on wind speed in the rest of this section. The error on wind direction [Fig. 14(b)] delivers a clearer picture, where wind speed emerges as the main error driver, followed by TKE and scan tilt. On the other hand, second-order statistics [Figs. 14(c) and 14(d)] exhibit a strong correlation with the TKE and scan tilt. It is noteworthy that the scan type did not emerge as significant feature for the description of the error standard deviation, as it could be already surmised from Figs. 7 and 8. This indicates that after the bias of the DBS is removed, the leftover random error of this scan is similar to that of the six-beam.

The feature selection method will guide the analysis of the random error on wind statistics that follows. More specifically, the error standard deviation, which is the primary metric for the random error, is calculated over subsets of data clustered based on the most important features. The discussion of the experimental results is also paired with the predictions from a semiempirical random error model along the line of the one presented for the bias in Sec. V C. This model has limited prognostic value (i.e., it has unknown validity outside of the present test site and climatology), but offers important insight into the underlying mechanism behind the observed error.

To assess the validity of Eq. (18) for the present dataset, we investigate the dependence of the error standard deviation of the horizontal velocity components on the TKE and scan tilt (Fig. 15). The dimensional error [Figs. 15(a) and 15(b)] approximates well a polynomial function and increases with the tilt. Even more interestingly, when dividing the error by the geometrical factor, $∑jNLOS∑k3Ai,j+2Aj,k2$⁠, to isolate its sampling component, the obtained “geometry-normalized error” curves collapse nicely on a single curve for both $U¯$ and $V¯$ [Figs. 15(c) and 15(d)]. This indicates that Eq. (18) is indeed sufficiently accurate to describe the error observed in the present dataset and that the sampling error can be approximated as $f(TKE)∼c1TKEc2$⁠. The dashed lines in Figs. 15(a) and 15(b) simply represent the error predicted by the present model after multiplying back by the geometrical factor. It is noteworthy that the error does not plateau to a constant value for low TKE, as would happen if the TKE-independent instrumental noise was significant. This seems to point to the fact that, as a consequence of the relatively close range at which measurements were collected (z = 119 m) and with the strict lidar data QC in place, instrumental error is dwarfed by the sampling counterpart, as already concluded in Sec. IV C.

After the error on the individual velocity components is known, Eq. (C5) can be used to estimate the resulting error on wind speed (Fig. 16, as function of TKE and tilt) and direction (Fig. 17 as a function of wind speed, TKE, and tilt). The agreement of the combined Eqs. (17) and (C5) is slightly worse than seen for $U¯$ and $V¯$ in Fig. 15 due to the intrinsic nonlinearity of the wind speed and direction functional form, but overall satisfactory. This analysis also shows that the well-known high error of wind direction for low wind speed is simply a direct consequence of the error propagation from $U¯$ and $V¯$ to $θ¯w$ and it is a mathematical feature that has little to do with the actual flow physics.

Finally, the omnidirectional errors on TKE and shear stresses can be obtained from the errors on Reynolds stresses in the fixed reference frame [see Eq. (C9)]. The results are shown in Fig. 19. The accuracy of the error model is comparable to that on the individual Reynolds stresses, also because of the linear relationship between these two flow variables and the $UjUk¯$ terms that were originally modeled.

The performance of the semiempirical error model in predicting all the error standard deviations calculated over the selected bins of important features is summarized in Fig. 20. The correlation is above 0.94 for all variables and was obtained through the use of the six tuning parameters shown previously (two for the mean flow, two for the Reynolds stresses of six-beam, and two for the Reynolds stresses of DBS). How generalizable these tuning parameters are in describing the dependence of the sampling error on TKE remains an open question and subject of future work.

Regardless of its universality, the present framework for the modeling of the error on wind statistics provides valuable insight into its driving factors, thanks to the complete separation of the effect of scan geometry and sampling error. While the latter is approximated through an empirical polynomial function of TKE (which happened to be quite general for different flow statistics and scans), the geometrical factors are calculated in closed form from matrices A, A⁺, and M, which are fully determined by the scan geometry and readily evaluated numerically. It is also noteworthy that random error is negatively affected by tilt due to the larger geometrical factor as it happened for the bias of the DBS (Sec. V C). One could argue that this is also related to the loss of symmetry, orthogonality and extreme elevation angles that occur when moving from regular to tilted scans. However, the prohibitive complexity of a closed-form solution of Eqs. (18) and (19) hindered a further breakdown of the error sources.

The development of the tilted profiling techniques is mainly motivated by the need to install a ground-based scanning lidar to measure wind speed vertical profiles close to operating turbines at AWAKEN,^3,4 where blockage and wakes are expected to be significant. Specifically, the candidate location of the main southerly inflow site “A1” was within 2.5 rotor diameters, D, of the H05 turbine of the King Plains wind farm.⁴ Two options were contemplated: option 1 was a deployment at the turbine base, which would have had the advantage of easing construction efforts and power access; option 2 was a deployment 1D south of H05. To identify the best scanning strategy for site A1, the field tests of tilted profiling scan were complemented by a numerical analysis of the AWAKEN site. The numerical part of the study leverages the virtual lidar tool presented by Robey and Lundquist⁴¹ deployed within the wind field simulated through state-of-the-art high-fidelity simulations by Sanchez Gomez et al.⁷² The simulations employ a nested mesoscale–microscale approach using the WRF model,⁷³ with NREL 5-MW wind turbines represented using a generalized actuator disk approach^74–76 in an LES environment. For the simulated case, the wind direction is from the south, which is the prevalent wind direction at the site, and the wind speed at 90 m is in region II of the NREL 5-MW turbine, corresponding to stable/high-thrust conditions, to maximize wake and blockage effect. The virtual lidar includes the time-resolution of the beams and averaging over the probe volume using a range gate of 120 ns or 18 m and pulse full-width, half-maximum of 320 ns or 48 m. The simulations include a high-resolution topography (1/3 arc s), which significantly impacted the flow upstream of the farm. Rather than assuming flow homogeneity upstream of the farm to assess the effect of the wind farm thrust, two sets of simulations are run using the same precursor and same setup: a case without the turbines (“no farm”), and a case with the turbines (“farm”) operating at nominal thrust. The difference between farm and no farm is a well-posed estimate of the wake and blockage effect on the undisturbed flow.

Virtual lidar experiments are carried out for 10 min in the farm velocity field using all 8 scan configurations under investigation. The ground-truth representative of the undisturbed flow are equivalent virtual measurements obtained in the no farm flow and with an ideal lidar that samples the scan region without any probe averaging or LOS projection. The difference between “lidar/farm” and “ideal/no farm” will encompass the effects of wakes, blockage, and lidar acquisition and retrieval error.

An overview of the farm-no farm flow at hub height is shown in Fig. 21 and highlights three types of flow distortions caused by the wind farm within the region of interest: (1) wakes behind the turbines, (2) a sharp gradient caused by individual turbine induction $∼2D$ upstream of the farm, and (3) a more gentle gradient attributable to wind farm induction or blockage that extends beyond 15D upstream.

The cropped area around turbine H05 is the region where the virtual lidar operates, and it is further expanded in Fig. 22. In these plots, we superpose the location of the lidar beams at hub height for all the scans and the deployment option 1 (top) and option 2 (bottom). We notice that even for the most tilted scan, the lidar probes a region deeply within the wind farm blockage area, thus ruling out the possibility that tilted profiling could be used to mitigate this large-scale flow deceleration. However, tilting the scan appears to be an effective strategy to avoid the near wake and the local induction zone where strong gradients occur.

The virtual lidar data analysis provides a thorough quantification of the potential benefit of tilted profiling. This is done by visualizing the difference of mean wind speed and direction between lidar estimates with farm and ideal sampling without farm at different heights within the rotor layer. Figure 23 highlights significant benefits on wind speed moving from regular profiling (first column) to a moderate and severe tilt for deployment option 1, where the lidar sits at the tower base, with a threefold reduction of the negative bias due to near wake and induction zone. The improvement is limited in the case of a lidar deployed 1D upstream, as most of the effects from turbine H05 have vanished. It is noteworthy that a zero error is never observed, which is expected since the difference between farm and no farm does not vanish up to 15D (see Fig. 21). This points to the fact that measuring a truly undisturbed inflow (i.e., the inflow that would occur without the turbines) is nearly impossible in practice. Even by placing a sensor at a large distance from the wind farm, the sampled flow would likely be a poor representation of the conditions experienced by the turbines due to spatial patterns induced by terrain and mesoscale circulation.

Figure 24 regarding wind direction error conveys a more complex scenario, where significant wind direction biases occur for the lidar at x = 0 for regular DBS and for the most tilted scans. The first effect is likely due to wake or induction effects occurring at the specific sampling locations of the DBS that are wrongfully interpreted as wind direction shift by the retrieval. The source of the bias for high tilt is consistent among all scans at both deployment locations and is due to an actual flow turning induced by the turbines 1–2D upstream.

The virtual lidar analysis points to the fact that tilted profiling is mostly beneficial in deployment option 1. It also indicates a slightly superior performance of the six-beam scan in mitigating wind direction artifacts for this specific flow. These results, paired with those from the field test, suggest that a moderately tilted six-beam is a suitable option for a deployment close to a turbine. In fact, this scan delivers an excellent wind speed (⁠ $Δ2=0.31$ m s⁻¹, $ρ=0.993$⁠) and wind direction (⁠ $Δ2=12.8°, ρ=0.993$⁠) prediction. TKE is also well estimated (⁠ $Δ2=0.61$ m² s⁻², $ρ=0.78$⁠) although it has some challenges in retrieving shear stress (⁠ $Δ2=0.31$ m² s⁻², $ρ=0.44$⁠) which requires further investigation. The moderate tilt is sufficient to at least halve the effect of turbine-induced flow distortion for the deployment option 1. The tilted DBS scans perform well for the estimation of the mean flow, but are affected by unacceptable bias of second order statistics (Sec. V C) and similar random error to the six-beam.

The present work proposes a modified version of traditional lidar profiling techniques that allow for a tilted scanning of the flow. The scan tilt is an additional degree of freedom offered to the lidar user to intentionally avoid regions that are expected to be contaminated by flow disturbances, such as those induced by buildings, topographic elements, or wind turbines. The design of tilted versions of the six-beam scan uses a constrained optimization approach, whereas tilt is introduced in the DBS by means of geometric modifications.

A three-month field comparison carried out at the NREL Flatirons Campus of wind statistics vs a collocated sonic anemometer showed that mean horizontal winds are reliably reconstructed by all the proposed scan geometries and at all tilts, with negligible biases and correlations exceeding 0.98 for both wind speeds and direction. Conversely, TKE and shear stresses are very sensitive to the scan tilt, and their accuracy deteriorates dramatically past moderate tilts.

A deeper analysis of the mean error revealed that the tilted DBS scans suffer from strong biases in TKE and shear stress. The proposed analytical model explains such bias as a consequence of the breakdown of the assumption of instantaneous flow homogeneity or cross-contamination, which generally causes an overestimation of the magnitude of the Reynolds stresses. The bias is exacerbated by tilt due to the loss of symmetry, orthogonality, and uniform and moderate elevation angles that characterize the geometry of the regular DBS.

On the other hand, the analysis of the random component of the error for all the scans was also coupled with a semiempirical model of error standard deviation that facilitated the interpretation of the observed trends and achieved a global correlation around 0.94 with observations. The error standard deviation on all wind statistics can be decomposed into a geometrical amplification factor, which is a function of the scan geometry only and completely known, and a sampling error due to the lidar finite sampling frequency and scan duration. The latter showed a monotonic trend as a function of ambient TKE, with higher error in more turbulent flows. Moreover, the sampling error, which is significant especially for Reynolds stresses over a relatively short 10-min window, is quite uniform across different scans and turbulent quantities, although it is magnified differently for different tilts and flow statistics based on the geometrical factors. Model and observations agree that larger random error amplification occurs at high tilts.

High-fidelity WRF-LES simulations of the AWAKEN site enhanced with a virtual lidar simulator shed light on the benefits of tilted lidar profiling for that experimental campaign. Tilted profiling can be useful for reducing the impact of the turbine wake and induction zone for a logistically convenient deployment at tower base, although it does not avoid the broader wind farm induction zone.

This work builds on previous studies of narrow VAD and arc scans and highlights the effectiveness and limitations of unconventional profiling techniques that represent a pragmatic and actionable alternative to more sophisticated computational fluid dynamics-based flow corrections that have been adopted mostly for sites with complex terrain. The proposed model for the DBS bias and the generalized random error model also proved to be effective at describing the observations. The main shortcoming of these models is the unknown applicability for other sites and environments, as part of their formulation is data-driven and was tuned and validated on the present dataset. Moreover, important lidar acquisition and processing issues such as the probe and time averaging and the instrumental noise have been neglected. Future work would require a generalization of the model for different flows, an estimation of the contribution of lidar-specific errors, and an analytical estimation of the sampling error, which would allow for a closed-form quantification of the measurement uncertainty. Furthermore, a more careful uncertainty-aware inversion algorithm could be introduced to compensate the inhomogenoeus SNR across the scanning volume for tilted scans. Also, we focused on readapted traditional conical scan types which exhibited smaller errors than arc scans.⁴⁸ However, future research should explore completely novel scanning patterns specifically designed for tilted profiling that could outperform the scans proposed in this study.

This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Additional work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 and Pacific Northwest National Laboratory. PNNL is operated by Battelle Memorial Institute for the Department of Energy under Contract No. DE-AC05-76RL01830. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award No. DE-SC0021110. This research has been supported by the U.S. National Science Foundation (Grant No. AGS-1565498). This research was performed using computational resources sponsored by the U.S. Department of Energy's Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory.

The authors have no conflicts to disclose.

Stefano Letizia: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Rachel Robey: Data curation (supporting); Software (equal); Supervision (equal); Validation (supporting); Writing – review & editing (equal). Nicola Bodini: Software (supporting); Supervision (equal); Validation (supporting); Writing – review & editing (equal). Miguel Sanchez Gomez: Software (supporting); Supervision (equal); Writing – review & editing (equal). Julie Lundquist: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal). Raghavendra Krishnamurthy: Writing – review & editing (equal). Patrick Moriarty: Funding acquisition (equal); Writing – review & editing (equal).

The lidar data collected during the pre-campaign test are available at www.a2e.energy.gov at channel awaken/nwtc.lidar.z02.b0, while wind statistics from the AWAKEN campaign collected through the moderately inclined six-beam scan at site A1 are available at channel awaken/sa1.lidar.z03.c1. M5 met tower data are also downloadable from https://wind.nrel.gov/MetData/135mData/M5Twr/. Python scripts to calculate the bias due to cross-contamination [Eq. (16)] and the geometrical error standard deviation factors [Eqs. (18) and (19)] are available at https://github.com/StefanoWind/Tilted_lidar_error.

A tailored QC is applied to correct and/or reject segments of sonic data affected by spikes and ramps and described in this appendix. The overall QC process is illustrated in Fig. 26, and it is applied to 10 min segments of data collected at 20 Hz. The process includes the following steps:

1. Raw segments for each velocity component and sonic temperature, T_s, are normalized by dividing them by the typical standard deviation calculated as the median of the standard deviation across all heights and for the whole period; calculating the median for the whole dataset mitigates the influence of data artifacts (especially spikes) that can exaggerate the standard deviation value at this stage.

2. An iterative despiking algorithm based on an adaptive median filter⁷⁷ is applied to each normalized signal iteratively for a maximum of 10 times; at each iteration, the detected spikes are replaced through linear interpolation; signals showing residual spikes after 10 iterations are discarded as a whole.

3. 10-min segments containing more than 10% of spikes are discarded to limit possible artifacts due to the interpolation.

4. For the same reason, signals presenting more than five consecutive spikes are discarded.

5. Finally, signals presenting linear ramps larger than 10 times the 95th percentile of the typical gradient in each signals are discarded; this is done since the despiking algorithm is vulnerable to ramps.

An example of the data artifacts that are filtered through the outlined QC pipeline is given in Fig. 27, where an uncorrupted sonic anemometer signal is contaminated by isolated spikes, consecutive spikes, and a 4-point ramp.

In this appendix we provide the mathematical details of the random error of the lidar retrievals based on linear error propagation.