Spectral correction of turbulent energy damping on wind LiDAR measurements due to spatial averaging

. Continuous advancements in pulsed wind LiDAR technology have enabled compelling wind turbulence measurements within the atmospheric boundary layer with probe lengths shorter than 20 m and sampling frequency of the order of 10 Hz. However, estimates of the radial velocity from the back-scattered LiDAR signal are inevitably affected by an averaging process within each probe volume, generally modeled as a convolution between the true velocity projected along the LiDAR line-of-sight and an unknown weighting function representing the energy distribution of the laser pulse along the probe length. 5 As a result, the spectral energy of the turbulent velocity ﬂuctuations is damped within the inertial sub-range, thus not allowing to take advantage of the achieved spatio-temporal resolution of the LiDAR technology. We propose to correct the turbulent energy damping on the LiDAR measurements by reversing the effect of a low-pass ﬁlter, which can be estimated directly from the power spectral density of the along-beam velocity component. LiDAR data acquired from three different ﬁeld campaigns are analyzed to describe the proposed technique, investigate the variability of the ﬁlter parameters and, for one dataset, assess 10 the corrected velocity variance against sonic anemometer data. It is found that the order of the low-pass ﬁlter used for modeling the energy damping on the LiDAR velocity measurements has negligible effects on the correction of the second-order statistics of the wind velocity. In contrast, the cutoff wavenumber plays a signiﬁcant role in spectral correction encompassing the smoothing effects connected with the LiDAR probe length. Furthermore, the variability of the spatial averaging on wind LiDAR measurements is investigated for different wind speed, turbulence intensity, and sampling height. The results conﬁrm 15 that the effects of spatial averaging are enhanced with decreasing wind speed, smaller integral length scale and, thus, for smaller sampling height. The method proposed for the correction of the second-order turbulent statistics of wind-velocity

, where " is the mean velocity calculated over a 5-minute sub-period, is the number of nonoverlapping sub-periods, and &'& is the mean of the entire velocity signal. For each dataset, the maximum relative absolute error for the mean is reported in the following table and in 11.67 At L 214, it is now reported: "The statistical steadiness of the LiDAR signals is estimated for both first-and second-order statistics. For the mean velocity, the absolute percentage error is calculated as follows: where " is the mean velocity as a function of height, z, calculated for the j-th sub-period with a duration of 5 minutes, sub-period, &'& is the mean wind velocity as a function of height for the entire velocity signals, while is the total number of sub-periods generated without overlapping". At L 224, it is reported: "For quality control purpose, signals with ( ≥ 15% or ! ≥ 40% are usually rejected (Foken and Wichura, 1996;Foken et al., 2004). The parameters ( and ! are calculated for each range gate and their maximum values for the selected datasets are reported in Table 2".

Point 4:
In your reply to point 27, you wrote that the periodogram method was replaced with Welch's algorithm and that similar results were obtained. In Matlab, using Welch's algorithm without overlapping and a single window is similar to applying the periodogram method but with a Hamming window instead of a rectangular window. It seems that in your data analysis, you have applied a single window (I might be wrong here). Therefore, this might explain the lack of substantial improvement in the power spectral density (PSD) estimates. Nevertheless, the periodogram method is known to be a poor power spectral density estimator. I recommend trying two to three segments with Welch's algorithm and 50% overlapping. This could significantly change the outcome of your fitting algorithm. R: We are thankful to the Reviewer for this further comment. As reported in the manuscript, before the calculation of the power spectral density, each signal is high-pass filtered with a cut-off wavenumber )' = 10 *+ m *% (L 317). Subsequently, the power spectral density for the Celina and SLTEST data are calculated with the Welch's algorithm without window overlapping and with a sub-window length inversely proportional to )' . In particular, the number of sampling points within a single window is given by: whereis the sampling rate and )' is the detrending frequency calculated from )' through Taylor's frozen-turbulence hypothesis (Taylor, 1938). Therefore, the number of non-overlapping windows has been calculated as: , = floor ; ./. , =, were ./. is the total number of samples of the velocity signal. The number of windows for each dataset is reported in the following For the XPIA dataset, only one window has been used due to the short duration of the velocity acquisition (14.5 minutes) and the need to characterize the lowest frequency available. Finally, a smoothing procedure in the wavenumber domain is applied following the Savitztky-Golay filter (L270). As pointed out by the Reviewer, the specific algorithm used to evaluate the spectra may affect the result of the correction procedure. In the manuscript, it is stated that the main parameter determining the amount of variance recovery is the cutoff wavenumber, .0 . Thus, we conducted a sensitivity analysis of .0 on the algorithm used for the calculation of the spectra. For this analysis, we leveraged the SLTEST dataset. In particular, for each LiDAR gate, the spectrum of the velocity signal is calculated through the following algorithms: 1. Welch algorithm with 3 windows and 50% overlapping; 2. Welch algorithm with 3 windows, 50% overlapping and smoothing; 3. Welch algorithm with , windows and no overlapping; 4. Welch algorithm with , windows, no overlapping and smoothing (as done in the manuscript). An example of the resulting spectra is reported in the figure below for the LiDAR signal collected at z= 10 m.
Subsequently, the correction procedure is carried out for each LiDAR velocity signal and spectral procedure. It is evident that the resulting cutoff wavenumber is practically insensitive to the different spectral procedures used. In your reply to point 29, the authors wrote that measurement noise in the lidar velocity records is unlikely to produce the big overestimation observed in Fig 5. I understand that my initial suggestion was maybe too vague. Therefore, I have reproduced such an overestimation in fig. 1 by filtering an idealized velocity spectrum using functions which emulates the presence or absence of noise at high frequencies. I have actually applied two low-pass filters with a different order for the sake of simplicity. In the left panel of fig. 1, one could assume that the red curve has a larger amplitude than the blue one because of measurement noise. The correction of the velocity spectrum is done by using the ratio H(f r ) between the black curve and the blue curve. Therefore, the blue curve is properly corrected. However, the noise in the red curve is massively amplified when applying the same correction, because H(f r ) erroneously assumed no measurement noise. The right panel of fig. 1 shows similar behavior as in the right panel of your Fig. 5. Even if the random noise is small, the application of an idealized spectral correction will likely produce a velocity spectrum with an unacceptable level of noise in the high-frequency range. The misalignment of the scanning beam with the instantaneous wind direction might also contribute to the overestimation of the corrected power spectral density in the high-frequency range. As stated by the authors, the different probe volume length could also affect the high frequencies velocity fluctuations, but I am not sure to what extent. I believe that different probe volume lengths are likely to affect more the low-frequency velocity fluctuations than the high-frequency ones. In this regards, your correction algorithm outperforms those based on an idealized spatial filtering function because it accounts for the random error in the spatially filtered velocity spectrum. R: We agree with the Reviewer about the impact of the instrument noise on the high-frequency portion of the spectrum. In the manuscript at L 304, it is now reported: "A possible explanation for the poor performance of these deconvolution models could be the presence of residual noise in the data, which is not accounted for in the models of Eq. 5 and 6. Other factors, like the different probe length used for the XPIA campaign ( = 50m) in contrast to = 30 used in the original study of this deconvolution model , are thought to have marginal effects on the estimation of the energy spectrum at higher frequencies."

Introduction
Over the last decades, wind Doppler Light Detection and Ranging (LiDAR) technology has provided compelling features to perform wind turbulence measurements within the atmospheric boundary layer (ABL) for different scientific and industrial pursuits, such as air quality, meteorology (Spuler and Mayor , 2005;Emeis et al. , 2007;Bodini et al. , 2017), aeronautic Doppler wind LiDARs were assessed against other measurement techniques, such as sonic anemometers and scanning Doppler wind radars, during the eXperimental Planetary boundary layer Instrumentation Assessment (XPIA) campaign Debnath et al. , 2017a, b;Choukulkar et al. , 2017;Debnath , 2018).
Different scanning strategies can be designed to characterize different properties of the ABL velocity field through LiDAR measurements (Sathe and Mann , 2013), while the highest spectral resolution is achievable by maximizing the sampling fre-30 quency and measuring over a fixed line-of-sight (LOS). Provided the use of a probe length, l, sufficiently small to perform wind measurements within the inertial sub-layer, e.g. at a height from the ground, z, with l < 2πz (Banerjee et al. , 2015), turbulence statistics of the wind velocity field can be retrieved through fixed scans while providing a spectral characterization of the inertial sub-layer (Iungo et al. , 2013). 3D fixed-point measurements can be performed by retrieving the radial velocity measured simultaneously by three or more LiDARs intersecting at a fixed position (Mikkelsen et al. , 2008;Carbajo et al. , 35 2014). In Mann et al. (2009), auto-and cross-spectral densities for the three velocity components were estimated through multiple scanning-LiDAR measurements.
Besides the easier deployment compared to the installation of classical meteorological towers, wind LiDARs tailored for investigations on atmospheric-turbulence currently provide probe volumes smaller than 20 m along the direction of the laser beam and sampling frequency higher than 1 Hz, which are welcomed features for studies on ABL turbulence.

40
A Doppler wind LiDAR allows probing the atmospheric wind field utilizing a laser beam, whose light is back-scattered in the atmosphere due to the presence of particulates suspended in the ABL. The velocity component along the laser-beam direction, denoted as radial or LOS velocity, is evaluated from the Doppler shift of the back-scattered signal. A pulsed Doppler wind LiDAR, like those used for the present work, emits laser pulses to perform quasi-simultaneous wind measurements at multiple distances from the LiDAR as the pulses travel in the atmosphere. The wind measurements performed over each probe volume 45 can be considered as the convolution of the actual wind velocity field projected along the laser-beam direction with a weighting function representing the radial distribution of the energy associated with each laser pulse. Therefore, LiDAR measurements can be considered as the result of low-pass filtering of the actual velocity field, where the characteristics of the low-pass filter are functions of the energy distribution of the laser pulse over the probe volume, probe length, and accumulation time (Frehlich et al. , 1998;Sjöholm et al. , 2009;Held and Mann , 2018).

50
A reduced variance of the wind velocity is generally measured with a Doppler wind LiDAR compared with that measured through a sonic anemometer due to the laser-pulse averaging and different size of the measurement volume. For single-point measurements performed with a Windcube 200S LiDAR and azimuth angle of the laser beam set equal to the mean wind direction, a variance reduction of 8% was predicted for a gate length of 25 m, while it was increased up to 20% for a gate length of 100 m (Cheynet et al. , 2017).

55
Attenuation of the measured turbulent kinetic energy due to the averaging over each probe volume can be corrected through a spectral transfer function introduced in Mann et al. (2009). For fixed scans, by leveraging the Taylor's frozen-turbulence hypothesis (Taylor , 1938;Panofsky and Dutton , 1984), the velocity energy spectrum is recovered through the deconvolution of the radial velocity with the weighting function representing the energy of the laser pulse. The critical part of this correction 2005; Lindelöw , 2008;Mann et al. , 2009). As it will be shown in this paper, corrections performed through this deconvolution procedure often do not provide a satisfactory accuracy for wind turbulence measurements.
Another method for spatial-averaging correction of wind LiDAR measurements was proposed in Brugger et al. (2016).
By assuming a linear averaging over each range gate and a Gaussian shape of the energy along the laser pulse, this method estimates the corrected velocity variance and the outer scale of turbulence by leveraging the von Kármán model of the secondorder structure function for the streamwise velocity (Von Kármán , 1948). The spatial averaging is included directly in the calculation of the structure function, following the work by Frehlich et al. (1998). In Brugger et al. (2016), compelling results were achieved comparing corrected LiDAR data with simultaneous and co-located data collected with an ultrasonic anemometer. However, these authors noticed residual systematic errors in the LiDAR corrected data, which might be related to the assumptions of the laser-pulse shape, the linear averaging process, or the von Kármán model of the structure function, 70 which was originally formulated for isotropic neutrally-stratified turbulence (Von Kármán , 1948).
In this work, a semi-empirical procedure is proposed to correct the damping of turbulent kinetic energy associated with wavelengths comparable to the LiDAR probe length for turbulent velocity measurements collected within the atmospheric surface layer (ASL), which is defined as the lower portion of the ABL where momentum and thermal fluxes are assumed to be constant (Stull , 1988). The ASL height can be quantified through the analysis of the turbulent fluxes or the variance of 75 the streamwise velocity as a function of height (Gryning et al. , 2016). In contrast to the above-mentioned methods for the correction of the streamwise velocity variance for LiDAR spatial averaging Sjöholm et al. , 2009;Brugger et al. , 2016), the correction method proposed in this paper does not require any a-priori information about the technical specifications of the used LiDAR systems, such as probe length or shape of the laser pulse. The proposed method allows to correct the second-order statistics of the streamwise velocity from spatial averaging by inverting the effects of a low-pass 80 filter, whose characteristics are directly determined from the power spectral density (PSD) of the LiDAR measurements. It is noteworthy that this method leverages surface layer similarity (Stull , 1988), thus it can only be applied for wind LiDAR measurements collected within the ASL.
The remainder of this paper is organized as follows: the theoretical aspects of the correction procedure are discussed in §2, while in §3 the experimental campaigns performed to collect the various LiDAR datasets are described. In §4, an assessment 85 of the proposed correction procedure is performed against sonic anemometry, while in §5 the correction procedure is tested for various LiDAR datasets. In §6, the spatial-averaging effects are investigated by varying the friction velocity, aerodynamic roughness length, and sampling height, thus for different mean wind speed and standard deviation. Finally, concluding remarks are reported in §7.
2 Correction procedure for the LiDAR velocity spectra Surface-layer scaling is typically used for spectral models of the wind speed assuming that the velocity integral length-scale is proportional to the height from the ground, z, and the Reynolds stresses can be normalized with the square of the friction 3 velocity, u τ . A classical approach to model the power spectral density of the streamwise velocity, S u , is the following: where f is frequency, n = f z/U is the reduced frequency, U is the mean advection velocity, and A and B are parameters 95 estimated through the best-fitting of the pre-multiplied energy spectra of the LiDAR velocity signals with Eq. 1. The term φ (≥ 1) represents a dimensionless dissipation for non-neutral atmospheric stability regimes, with φ equal to 1 for neutrallystratified surface-layer flows (Kaimal et al. , 1972). For this work, we only consider near-neutral atmospheric conditions and slight variations connected with atmospheric stability are embedded in the coefficient A.
The spectral model of Eq. 1 is typically referred to as blunt model (Olesen et al. , 1984) or Kaimal model (Kaimal et al. , 100 1972;IEC , 2007;Worsnop et al. , 2017;Risan et al. , 2018), and the parameter A is typically assumed equal to 105 (Kaimal et al. , 1972), later revised to 102 (Kaimal and Finnigan , 1994), and B equal to 33. It is noteworthy that within the inertial sub-layer, the pre-multiplied spectra scale as n −2/3 , while the maximum value occurs for a reduced frequency equal to 1.5/B, which corresponds to the wavenumber k p = 3π/(B z).
Considering the Cartesian reference frame, (x 1 , x 2 , x 3 ), where the coordinates are aligned with streamwise, transverse, and 105 vertical directions, respectively, the wind speed measured by a pulsed Doppler wind LiDAR can be modeled as the convolution between the projection of the wind velocity u = (u 1 , u 2 , u 3 ) along the laser-beam direction, n = (n 1 , n 2 , n 3 ), with a weighting function, φ, representing the energy distribution of the laser pulse within a probe volume (Sjöholm et al. , 2009;Mann et al. , 2009;Cheynet et al. , 2017):

110
where v r is the radial or LOS velocity measured along the laser-beam direction, n, at a radial distance x from the LiDAR. The probe length is l, while s is the radial position within the considered probe volume. The weighting function, φ, is normalized to unit integral. If the Doppler frequency is determined as the first moment of the signal PSD with the background subtracted appropriately, then the weighting function can be expressed as (Banakh and Werner , 2005;Mann et al. , 2009;Cheynet et al. , 2017): For the LiDAR Windcube 200S, the following weighting function can also be used (Lindelöw , 2008;Mann et al. , 2009): In the spectral domain, the Fourier transform of Eq. (3) is:

120
4 while for Eq. 4 is: where k = 2πf /U is the wavenumber evaluated through the Taylor's frozen-turbulence hypothesis (Taylor , 1938). As shown in Mann et al. (2009), the measured velocity spectrum, S L , can be modeled as: 125 where: k = (k 1 , k 2 , k 3 ) is the wavenumber vector and summation over repeated indices is assumed. In Eq. 7 , Φ ij (k) is the spectral tensor obtained as Fourier transform of the Reynolds stress tensor and ϕ(k · n) is the Fourier transform of the convolution function. When the laser beam stares along the mean wind direction with a relatively low elevation angle, namely with n ≈ (1, 0, 0), the PSD of the radial velocity, S L , is equal to the product between the spectrum of the actual radial velocity, S L (k 1 ), and the square of the Fourier transform of the weighting function, ϕ(k 1 ): Equation 8 shows that the spectrum of the measured radial velocity, S L , is equal to the true velocity spectrum,Ŝ L , low-pass filtered with a certain transfer function. In this work, the latter is modeled as: where α and k T h represent the order and cutoff wavenumber, respectively, of a low-pass filter (Ogata , 2010). The symbol 135· is used to differentiate the analytical model of the low-pass filter from its empirical estimate through the ratio between the fitted Kaimal spectrum and the PSD of the LiDAR velocity, ϕ 2 * . These features of the low-pass filter and, thus, of the LiDAR measuring process, are functions of the LiDAR probe length, the elevation angle of the laser beam, the relative angle between wind direction and azimuth angle of the laser beam, accumulation time, and characteristics of the laser pulse. Therefore, it is highly challenging to predict these parameters a priori, while it is advisable to estimate α and k T h directly from the specific 140 LiDAR data under analysis. To this aim, we propose the following iterative procedure to correct the effects of the spatial averaging on wind LiDAR measurements, which is summarized in the flow chart of Fig. 1.
First, the pre-multiplied spectrum of the radial velocity projected in the horizontal mean wind direction is fitted with the spectral model of Eq. 1 only for wavenumbers smaller than k T h,0 = 2π/l. Indeed, we expect to observe significant spatialaveraging effects for turbulent length scales smaller than the probe length, l. For wavenumbers higher than the selected cutoff 145 value, the ratio between the fitted Kaimal spectrum and the PSD of the LiDAR velocity, ϕ 2 * , is calculated to quantify the effect of the energy damping due to the LiDAR measuring process. Subsequently, the LiDAR-to-Kaimal ratio, ϕ 2 * , is fitted with Eq. 9 through a least-square algorithm to estimate the filter order, α, and provide an updated value for the cutoff wavenumber, k T h . This process is iterated until convergence on the parameter k T h is achieved (for this work, the convergence condition imposed is a variation of k T h smaller than 1% of the previous value). If during the iterative process, k T h achieves a value equal or 150 smaller than that corresponding to the spectral peak, k p , then the procedure is arrested and a warning is dispatched indicating that the correction procedure was not successful. This warning condition never occurred for all the data analyzed in this work.
Furthermore, it should be considered that when k T h achieves values close to k p , the part of the velocity spectrum, S u , used for the fitting procedure with Eq. 1 can be so limited to jeopardize the accuracy of the fitting procedure. Once convergence in k T h is achieved, the corrected velocity spectrum,S L (k), is calculated as: It is noteworthy that in contrast to existing models using pre-defined functions to correct the energy damping of the velocity fluctuations, see e.g. Eqs. 5 and 6 (Sjöholm et al. , 2009;Brugger et al. , 2016;Cheynet et al. , 2017), which require information about the LiDAR probe length and the energy distribution over a pulse, the proposed procedure calculates the characteristics of the damping on the LiDAR velocity signals directly from the experimental data, which leads, as it will be shown in the 160 following, to enhanced accuracy in the correction of the LiDAR velocity spectra. On the other hand, the proposed procedure leverages the surface-layer similarity for the Kaimal spectral model for the streamwise velocity and, thus, it can only be applied for wind LiDAR measurements collected within the ASL. 6 The present study is based on wind LiDAR measurements collected from three different experimental campaigns. The first 165 dataset was acquired during the period June 9-24, 2018 at the Surface Layer Turbulence and Environmental Science Test (SLTEST), which is part of the U.S. Dugway Proving Ground facility in Utah (GPS location: 40 • 08 07"N, 113 • 27 04"W , UTC offset −6 h). Characterized by an elevation variability of 1 m every 13 km (Kunkel and Marusic , 2006;Metzger and Klewicki , 2001), this facility is located in the South West of the Great Salt Lake and extends for 240 km and 48 km along North-South and East-West directions, respectively. An aerial view of the SLTEST facility is reported in Fig. 2a. During the 170 experiment, the prevailing wind direction was from North-North-East.
The second field campaign was carried out at a test site in Celina, TX (GPS location: 33 • 17 35.3"N, 96 • 49 17.5"W, UTC offset −5 h), which is a relatively flat terrain with a certain variability in land cover (Fig. 2b). For these two field campaigns, wind velocity measurements were performed with a Streamline XR scanning Doppler pulsed wind LiDAR manufactured by Halo Photonics, whose technical details are reported in Table 1 Source of Google Earth. Black crosses represent the instrument locations. In (c), each LiDAR is labeled with its respective name. and sonic anemometers are analyzed to assess the proposed spectral correction procedure of LiDAR measurements.
For the Celina and SLTEST field campaigns, the regime of the atmospheric stability was monitored through sonic anemometers mounted at a height of 3 m in the proximity of the LiDAR location. The sampling frequency of the sonic anemometer data was 20 Hz, while atmospheric stability was characterized through the Obukhov length calculated as follows (Monin and Obukhov , 1954): where κ = 0.41 is the von Kármán constant, g is the gravitational acceleration, w θ v is the sensible heat flux, θ v is the average virtual potential temperature (in Kelvin) and u τ,S is the friction velocity calculated from sonic anemometer data as (Stull , 1988): To avoid effects of thermal stratification and buoyancy on our analysis, only datasets acquired under near-neutral conditions are considered, which are selected by imposing the threshold: |z/L| ≤ 0.05 (Kunkel and Marusic , 2006;Liu et al. , 2017).
The LiDAR velocity signals undergo a quality control process to ensure statistical significance and accuracy of the measurements. Only datasets with a variability of the 10-minute averaged wind direction within the range ±20 • have been considered to avoid significant offset between the LiDAR azimuth angle and the instantaneous wind direction (Hutchins et al. , 2012).

210
The quality of the LiDAR signals is then checked based on the intensity of the back-scattered signal. For the Windcube 200S LiDAR, the samples with a carrier-to-noise ratio (CNR) higher than -25 dB are selected, while for the Streamline XR LiDAR data are analyzed only if the intensity of the back-scattered signal is higher than 1.01.
The statistical steadiness of the LiDAR signals is estimated for both first-and second-order statistics. For the mean velocity, the absolute percentage error is calculated as follows: where U j is the mean wind velocity as a function of height, z, calculated for the j-th sub-period of the velocity signal with a duration of 5 minutes, U tot is the mean wind velocity as a function of height for the entire velocity signal, while N is the total number of sub-periods generated without overlapping. A similar parameter ε σ is calculated for the second-order statistics (Foken and Wichura , 1996): where CV j is the variance of the j-th sub-period with a duration of 5 minutes, while CV tot is the variance of the signal over the entire period. For Celina and SLTEST campaigns, the parameters ε M and ε σ were calculated for 1-hour periods, while for the XPIA deployment the whole 14.5-minute record was analyzed. For quality control purpose, signals with ε M ≥ 15% or ε σ ≥ 40% are usually rejected (Foken and Wichura , 1996;Foken et al. , 2004). The parameters ε M and ε σ are calculated for 225 each range gate and their maximum values for the selected datasets are reported in Table 2.
Subsequently, a gradient-based procedure is used to remove outliers from the LiDAR radial velocity signals. Specifically, the partial derivative in time of the radial velocity is calculated through a second-order central finite-difference scheme, and velocity samples with absolute partial derivative larger than 15 times the respective median value calculated over the entire signal are marked as outliers and replaced through the inpaint-nans function available in Matlab (D'Errico , 2004). The used threshold 230 value is selected based on a sensitivity analysis. Based on the above-mentioned quality-control procedure, five datasets were selected, whose details are reported in Table 2.
The radial velocity, V r , measured by a Doppler wind LiDAR, as explained in §2, is expressed as: . Description of the selected datasets: Φ is the LiDAR elevation angle, fS is the sampling frequency, l is the probe length, IST is the maximum value of the non-stationary index, uτ is the friction velocity, and z0 is the aerodynamic roughness length. where θ and Φ are the LiDAR azimuth and elevation angles, respectively, θ w is the wind direction, V h and W are the horizontal 235 and vertical wind velocities, respectively. As previously mentioned, for the SLTEST and Celina field campaigns, LiDAR measurements were carried out with the azimuth angle equal to the mean wind direction and very low elevation angles (Table   2). Therefore, we can calculate an approximation of the horizontal wind speed as: which is referred to as horizontal equivalent velocity. In the following, U eq is considered to calculate the streamwise velocity 240 spectrum. Furthermore, the variance of the radial velocity is the first-order approximation of the streamwise-velocity variance given the above-mentioned setup constraints (Eberhard et al. , 1989;Sathe and Mann , 2013).

Assessment of the LiDAR spectral correction against sonic anemometry
In this section, the procedure proposed in §2 to correct the energy damping in the LiDAR velocity measurements due to the energy pulse distribution over the probe volume is assessed against sonic anemometry by leveraging the XPIA dataset, whose 245 characteristics are summarized in Table 2. LiDAR fixed scans were performed with an elevation angle of 5 • to have one range gate in the proximity of a sonic anemometer installed on the BAO tower at a height of 100 m. The LiDAR probe length used for that experiment was equal to 50 m, while the sampling rate was set equal to 2 Hz.
Based on the instantaneous wind direction measured by the sonic anemometer and neglecting the vertical velocity due to the very low elevation angle of the LiDAR laser beam, the horizontal equivalent velocity, U eq , is calculated from the LiDAR radial 250 velocity through Eq. 15 and it is reported in Fig. 3 with a blue line. The LiDAR equivalent velocity, U eq , is then high-pass filtered to remove low-frequency non-turbulent velocity fluctuations, using the following spectral transfer function:  where k co is the cutoff wavenumber, which should be smaller than k p to avoid effects on the spectral peak. The parameter β is set equal to 100 to generate a sufficiently sharp filter across the cutoff wavenumber, k co (Hu et al. (2020)). The PSD of the 255 velocity signal high-pass filtered with a cutoff wavenumber k co = 1.26 × 10 −3 m −1 is reported in Fig. 4(a) with a grey line.
In case significant noise in the velocity spectra is observed in the proximity of the Nyquist wavenumber, see e.g. Debnath (2018), as for this velocity signal, a denoising procedure is then applied to remove possible noise effects on the velocity signals.
Following the wavelet-transform-based procedure proposed by , the velocity signal is decomposed in a 10-level orthogonal wavelet basis. For each level, a soft-threshold selection is applied to the wavelet coefficients to remove those related 260 to noise. The estimated noise-free wavelet coefficients, d jk , are calculated as: where j = 1, ..., 10 is the number of levels in the wavelet basis; k = 1, ..., 2 j and w jk are the coefficients of the discrete wavelet transform of the original signal. T j represents a noise-based threshold for the j th level, that for this work is set to (To et al. ,

2009):
where med(·) stands for the median value. Finally, the denoised signal is reconstructed in time through the modified wavelet coefficients d jk . The spectrum of the denoised LiDAR velocity signal is reported in Fig. 4(a) with a light-blue line.
For modeling purpose, the velocity spectra are then smoothed in the wavenumber domain following the Savitztky-Golay filter (Savitzky and Golay , 1964), by using a second-order polynomial function and windows with the width equal to int[10(160k) 0.5 ], 270 where k is in m −1 and int is rounding to the closest integer number (Balasubramaniam , 2005). The result is an increased level of smoothness moving towards the Nyquist wavenumber. The LiDAR velocity spectrum resulting from the de-noising and smoothing procedures is reported in Fig. 4(a) with a blue line. The PSD of the LiDAR velocity is then fitted through the spectral model (Eq. 1) producing the following fitting parameters: A = 23.3, B = 24.7. The resulting Kaimal spectrum is plotted in Fig. 4(a) with a black line.

275
A deviation of the LiDAR velocity spectrum from the −5/3 scaling of the inertial sub-range is observed due to the LiDAR measuring process over the probe volume. The ratio between the fitted Kaimal spectrum and the LiDAR velocity spectrum, |ϕ * | 2 in Fig. 4(b), is then fitted with Eq. 9 to estimate the low-pass filter of order α and cutoff wavenumber, k T h . For this LiDAR velocity signal, α is equal to 0.774 and k th l/(2π) ≈ 0.8, with an R 2 value of 0.702, which confirms the proposed model is a good approximation for the damping of the velocity fluctuations over the LiDAR probe volume. In Fig. 4(b), the 280 weighting functions of Eqs. 5 and 6 are also reported for a probe length l =50 m. spectrum (black) and −5/3 slope (black dashed); (b) |ϕ * | 2 (blue), |ϕ * | 2 fitted with Eq. 9 (black); predictions from Eq. 5 (bright green) and Eq. 6 (dark green) are reported as well.
To assess the accuracy of the estimated low-pass filter in representing the LiDAR averaging process over a probe volume, first we apply the estimated low-pass filter to the simultaneous and co-located sonic anemometer velocity signal. The horizontal velocity retrieved from the sonic anemometer is first down-sampled with the sampling frequency of the LiDAR measurements, namely 2 Hz, using the Matlab function "decimate" with a finite-impulse response (FIR) low-pass filter with order equal to 10 285 (Weinstein , 1979). The resulting down-sampled velocity signal is reported with a red line in Fig. 3, and the respective PSD in Fig. 5(a). Subsequently, the down-sampled sonic-anemometer signal is low-pass filtered with the filter of Eq. 9 modeled only using the LiDAR data (yellow line in Fig. 3 and in Fig. 5(a)). The comparison in Fig. 3 of the sonic anemometer signal down-sampled and low-pass filtered with the LiDAR raw signal already highlights a very good agreement, which suggests that the energy damping carried out by the LiDAR pulse over the probe volume is well represented through the proposed low-pass 290 filter. This feature is further corroborated by the respective spectra reported in Fig. 5(a). Specifically, the spectrum of the LOS velocity has the same slope in the inertial sub-range of the sonic-anemometer signal down-sampled and low-pass filtered, while some differences are observed for lower frequencies, which are most probably due to the different size of the measurement volume of the two instruments, namely 50 m for the LiDAR and 0.3 m for the sonic anemometer.
The comparison between LiDAR and sonic anemometer data is now presented through a linear regression analysis, which 295 is reported in Fig. 6. The LiDAR horizontal equivalent velocity, U eq , is analyzed against the horizontal wind speed measured by the sonic anemometer before (Fig. 6(a)) and after ( Fig. 6(b)) the low-pass filtering. All the linear regression parameters improve for the low-pass filtered sonic anemometer data: the slope increases from 0.878 to 0.962, the R-square value increases from 0.88 to 0.904, and the correlation coefficient, ρ, increases from 0.88 to 0.904.
We now aim to correct the LiDAR velocity signal from the energy damping due to the laser pulse distribution over the probe 300 volume. First, the LiDAR velocity spectrum is corrected by using the existing models of Eqs. 5 and 6 with l = 50 m, i.e. the used LiDAR probe length. As shown in Figs. 4(b) and 5(b), these correction methods largely over-estimate the turbulent energy for wavenumbers larger than k T h or, in other words, the characteristic length scale should be smaller than the LiDAR probe length to provide reasonable spectral corrections. A possible explanation for the poor performance of these deconvolution models could be the presence of residual noise in the data, which is not accounted for in the models of Eqs. 5 and 6. Other 305 factors, like the different probe length used for the XPIA campaign (l = 50 m) in contrast to l = 30 m used in the original study of this deconvolution model , are thought to have marginal effects on the estimation of the energy spectrum at higher frequencies.
According to the correction technique proposed in this paper, the LiDAR velocity spectrum can now be corrected for the averaging process by reversing the effect of the estimated low-pass filter through Eq. 10. The corrected velocity spectrum is 310 reported in Fig. 5(b) with a black line. The velocity spectrum of the corrected LiDAR velocity signal clearly shows that the expected slope of −5/3 in the inertial sub-range is recovered, while the spectral energy for lower wavenumbers are practically unchanged.
14 For the SLTEST and Celina field campaigns, LiDAR velocity measurements were collected for time periods between 2 and 3 315 hours (see Table 2). The procedure used to obtain the streamwise velocity spectra is the following: for each LiDAR velocity signal, the high-pass filter of Eq. 17 is applied with a cutoff wavenumber k co = 0.001 m −1 to remove low-frequency velocity fluctuations connected with atmospheric mesoscales. The PSD of each velocity signal is then calculated with the pwelch function implemented in Matlab (Welch , 1967) without window overlapping and window width corresponding to k co . Subsequently, smoothing of the velocity spectra is carried out through the Savitzky-Golay filter, as detailed in the previous section 320 (Savitzky and Golay , 1964;Balasubramaniam , 2005).
The mean values and variance of the LiDAR equivalent velocity are plotted in Figs. 7(a) and (b), respectively. For the mean velocity field in Fig. 7(a), while a logarithmic region is generally observed at the lower heights, a noticeable difference in terms of terrain roughness between the SLTEST and Celina sites reflects in different vertical intercepts and respective aerodynamic roughness length. The latter is estimated to be equal to 0.021 mm for SLTEST and, as average, 37 mm for the Celina site.

325
The vertical profiles of streamwise velocity variance are reported in Fig. 7(b) as a function of height. For the datasets collected at the Celina site, a general increase of the velocity variance is observed with increasing height. Specifically, for the dataset Celina1, after achieving a maximum value at height z ≥ 40 m, a quasi-logarithmic reduction of the velocity variance is observed with increasing height, which is in agreement with previous laboratory and numerical studies of canonical boundary layer flows (Kunkel and Marusic , 2006;Meneveau and Marusic , 2013). A logarithmic reduction of the velocity variance with increasing height is also observed for the SLTEST dataset throughout the entire height-range.
For the SLTEST dataset, the PSD of the LiDAR velocity signals acquired at the different gates from 10 m up to 60 m with a vertical spacing of 1 m are plotted in Fig. 8. A departure from the expected −5/3 slope in the inertial sub-range is observed for wavenumbers larger than 0.07 m −1 . The Kaimal spectra, which are obtained by fitting the measured LiDAR velocity spectra only for wavenumbers lower than the respective k T h for each height, are reported in Fig. 9(a) for the lowest and highest range 335 gates. The ratio between the LiDAR and Kaimal spectra, |ϕ * | 2 , is then calculated and fitted through Eq. 9 to estimate the order and cutoff wavenumber of the respective low-pass filter ( Fig. 9(b)). For the lowest gate at height of 10 m, the fitting procedure estimated α = 5.08 and k T h = 0.065 m −1 , while for the highest range gate at height of 60 m α = 4.16 and k T h = 0.066 m −1 .
Since the correction procedure is based on two consecutive best-fit operations, the robustness of the model is assessed for each LiDAR gate through the R-square value of the respective fitting procedure. All the datasets generally show a very 340 good agreement between experimental data and the spectral model, i.e. 0.82 ≤ R 2 ≤ 0.98, with the SLTEST dataset showing the highest level of agreement (R 2 ≥ 0.96). Accuracy in modeling the actual energy damping due to the LiDAR measuring process through the low-pass filter of Eq. 9 is quantified through the R-square value of the fitting procedure of the experimental energy damping, ϕ 2 * , with the analytical model of Eq. 9. The R-square values result to be always larger than 88%, corroborating the good approximation for the proposed model.  The proposed spectral correction of the LiDAR measurements is now applied to all the datasets collected at the Celina and SLTEST sites (see Table 2). As explained in §2, the first step of the proposed procedure consists in fitting each velocity spectrum with the spectral model of Eq. 1. In the right column of Fig. 10, the results of this operation are reported for all the datasets of the SLTEST and Celina field campaigns. The spectral model (depicted with a red line) has been fitted on the uncorrected LiDAR spectrum (blue lines) using a cutoff wavenumber, k T h , estimated for each velocity signal using the iterative 350 procedure illustrated in Fig. 1. In the figures, line colors become darker with increasing height. For the sake of clarity, the fitted Kaimal spectrum is only shown for the highest LiDAR range gate.
The second step of the correction procedure consists of approximating the LiDAR-to-Kaimal spectral ratio with the low-pass filter of Eq. 9. Firstly, the energy ratio is quantified for each LiDAR range gate, as reported in the left column of Fig. 10. We can observe that the ratio always settles about the unit at the lowest range of the spectral domain, while it monotonically reduces 355 from the cutoff wavenumber towards the Nyquist wavenumber, which is an effect of the LiDAR measuring process. By plotting |ϕ * | 2 as a function of (k/k T h ) α , all the estimated transfer functions practically collapse on the same curve for measurements collected at different heights. The latter is then compared with the analogous of Eqs. 5 and 6.
As the last step, to retrieve the corrected LiDAR velocity spectra, the original spectrum is divided by the modeled correction function, |φ| 2 (Eq. 10). These corrected LiDAR velocity spectra are reported on the right column of Fig. 10, where we can 360 observe that the −5/3 slope of the inertial sub-range is always recovered. The LiDAR spectra corrected with the models of Eq. 5 and 6 are also reported to highlight the improved accuracy achieved through the proposed method. In particular, it is observed that the existing models of Eqs. 5 and 6 always under-estimate the spectral energy attenuation; thus, for the actual choices of gate length and sampling rates, a data-driven approach is preferred to correct the LiDAR smoothing effect.
The corrected variance of the LiDAR velocity signals is compared with the respective quantity calculated for the raw LiDAR 365 data in Fig. 11(a). As expected, the wall-normal profile of the variance calculated as integral of the de-convoluted spectrum is considerably larger than the respective value obtained from the convoluted spectrum, indicating that the underestimation related to the spatial averaging is significant. To quantify the effects of the spectral correction on the LiDAR data, the relative percentage increment of variance is calculated from the smallest frequency up to the noise-free high-frequency content as in where σ 2 C and σ 2 U are the corrected and uncorrected, respectively, streamwise velocity variance. The parameter % is reported as a function of height in Fig. 11(b). The underestimation in the velocity variance through the LiDAR measurements seems to change with the wall-normal location for the SLTEST dataset and the highest portion of Celina1 (z > 30m); for the remaining datasets, the percentage error does not change with height. To clarify this aspect, in the next section we will investigate the 375 variability of the effects of the LiDAR spatial averaging for different mean wind speed, standard deviation and sampling height.
For the SLTEST dataset (see Table 2), the correction of the velocity variance obtained with the proposed method (red marker in Fig. 12a) is compared with those obtained from Eq. 5 (dark green symbols), Eq. 6 (light green marker) and the method of  (Risan et al. , 2018). Consistently with the spectra 380 of Fig. 10, the correction methods of Eqs. 6 and 5 under-estimate the effects of spatial averaging of the streamwise velocity variance and they do not allow the complete recovery of the -5/3 slope of the inertial sub-range. The correction based on the structure function proposed by Brugger et al. (2016) leads to a variance distribution larger than what is estimated by the new method based on the Kaimal spectral model, yet with percentage correction in the same order of magnitude (Fig. 12b). (b) percentage correction of the velocity variance, % . The marker colors are black for the raw LiDAR data, red for the proposed method, light green for Eq. 6, dark green for Eq. 5, and blue for the method proposed by Brugger et al. (2016).
We now focus on the variability of the parameters of the low-pass filter of Eq. 9 among the various datasets. Median and 385 interquartile (IQ) range calculated over the height for the various datasets are reported in Table 3 for both order, α, and cutoff wavenumber, k T h , of the low-pass filter of Eq. 9. First, the order of the low-pass filter, α, is found to be roughly constant over the height for the datasets Celina1 and Celina3, while it decreases with height for the SLTEST and Celina2 datasets. Among the various datasets analyzed, the filter order α has values from 2.2 up to 4.9, a variation entailing limited effects in the correction of the spatial averaging on the LiDAR measurements (not shown here for the sake of brevity). The mean IQ range of α among 390 the various datasets is 0.399. The cutoff wavenumber of the low-pass filter, k T h , is practically constant with height (mean IQ range of 0.022) and has a mean value among the various datasets of: k T h l/(2π) =0.163.
6 Variability in spatial filtering with mean wind speed, turbulence intensity and sampling height To investigate effects of the spatial averaging on wind LiDAR measurements for different mean wind speed, turbulence intensity, and sampling height of the velocity signals, synthetic turbulent velocity spectra are generated using the spectral model of 395 Eq. 1, while the energy damping connected with the LiDAR measuring process is estimated through Eq. 9 by using a filter order α = 3 and (k T h l) = 0.95, in analogy to the respective experimental values reported in Table 3. Within the inertial sublayer, namely for heights smaller than about 30% of the surface layer height , the mean streamwise velocity for near-neutral stability conditions can be modeled through the following logarithmic law (Monin and Obukhov , 1954;Stull , 1988): where κ = 0.41 is the Von Kármán constant. Furthermore, we should expect a logarithmic decrease in the velocity variance with increasing wall-normal distance (Townsend , 1976), as follows: where δ is the outer scale of turbulence, namely the boundary layer height, G 1 = 0.98 is the Townsend-Perry constant (Baars 405 and Marusic , 2020), while H 1 might be dependent on the characteristics of the specific boundary layer flow under investigation.
Previous field campaigns performed at the SLTEST site have quantified H 1 = 2.14 .
According to the spectral model of Eq. 1, the velocity variance can be obtained by integrating S u in the spectral domain, which leads to: 23. The probe length, l, is varied from 10 m up to 100 m with a step of 10 m.
The parameter % , which is defined in Eq. 20, is used to quantify the effects of the LiDAR spatial averaging on the variance of the wind velocity. In Fig. 13(a), it is observed that, as expected, % increases with increasing probe length, which indicates that the probe length is the main root cause of spatial averaging. Furthermore, it is noteworthy that for a given value of l, 420 % decreases with increasing u τ and, thus, mean wind speed, U . Indeed, by fixing the probe length, the cutoff wavenumber, k T h , representing the spatial averaging is, in turn, fixed while the reduction of % with increasing U can be explained in the perspective of the Taylor frozen-turbulence hypothesis (Taylor , 1938). In other terms, with increasing U , a turbulent spectrum will shift towards lower wavenumbers (denominator in Eq. 1) and, thus, reducing the percentage of the spectral energy that is filtered for wavenumbers larger than k T h .

425
The second test case, whose results are reported in Fig. 13(b), is performed by varying the aerodynamic roughness length within the range: z 0 = [10 −5 , 10] m, while keeping fixed u τ = 0.5 m s −1 , and sampling height z/δ = 0.3 (B=33). In this case, variations of z 0 affect directly the mean velocity (Eq. 21), while the velocity standard deviation is unchanged (Eq. 22).
Therefore, this study can be considered as a test to investigate the effects on the LiDAR spatial averaging due to variation of the wind turbulence intensity, which might be connected to different site terrain roughness, variations of wind direction, and 430 atmospheric stability regime. Specifically, for a fixed probe length, an increase of aerodynamic roughness length leads to a reduction of the mean velocity for a given height and friction velocity, and a shift of the turbulent spectrum towards higher wavenumbers and, in turn, enhanced effects of the spatial averaging on the LiDAR measurements.
The last case study is performed for a given wind condition, namely with u τ = 0.5 m s −1 and z 0 = 10 −4 m, and by varying the measurement height within the range z/δ = [0.1, 0.6]. The vertical variability of the parameter B = [65, 228] has been 435 selected equal to that measured for the SLTEST dataset, and the respective values of the wavenumber of the spectral peak, k p , are reported in the side panel of Fig. 13(c). For each height, the velocity variance is calculated via Eq. 22 and the corresponding value of A is obtained through Eq. 23. The resulting percentage reduction of variance is plotted in Fig. 13(c). For a fixed probe length, it is observed that % generally decreases with height in a way similar to what has been observed experimentally in Fig. 11(b), and the variations of % are strongly dependent on the variations of the parameter k p . In other words, with increasing height, the general reduction of k p leads to a smaller percentage of the spectral energy of the velocity signal present for wavenumbers larger than k T h , which is fixed once the LiDAR probe length is selected. Therefore, smaller effects of the LiDAR spatial averaging occur with increasing height for a given wind condition and probe length.

Conclusions
Pulsed Doppler wind LiDAR technology is gradually achieving compelling technical specifications, such as probe lengths 445 smaller than 20 m and sampling frequencies higher than 1 Hz, which are instrumental to investigate atmospheric turbulence with length scales typical of the inertial sub-range. However, the emission of a laser pulse over the probe volume to measure the radial velocity entails a spatial smoothing process leading to damping on the measured variance of the velocity fluctuations. Existing methods propose to correct the effects of spatial averaging on LiDAR measurements using as input technical specifications of the used LiDAR systems, such as probe length and pulse energy distribution, which might not be available 450 and, thus often approximated with analytical functions. According to previous works, and also confirmed through this study, existing methods have limited accuracy in correcting the LiDAR velocity fluctuations.
In this work, we have proposed to correct the measured LiDAR velocity signals by inverting the effects of a low-pass filter representing the energy damping on the velocity fluctuations due to the LiDAR measuring process. The filter characteristics, namely order and cutoff wavenumber, are directly estimated from the spectrum of the LOS velocity under investigation. Specif-455 ically, the spectrum of the LiDAR velocity signal is fitted through the Kaimal spectral model for streamwise turbulence only for wavenumbers lower than a cutoff value for which the slope of the LiDAR velocity spectrum is observed to deviate from the expected −5/3 slope typical of the inertial sub-range. The ratio between the LiDAR and the Kaimal spectra is then fitted with the analytical expression of a low-pass filter to estimate the order and cutoff wavenumber. An iterative procedure is proposed to estimate order and cutoff wavenumber of the low-pass filter. The modeled low-pass filter is then reverted on the LiDAR data 460 to correct the LiDAR measurements and produce more accurate second-order statistics and spectra of the streamwise wind velocity. It is noteworthy that the Kaimal spectral model leverages surface layer similarity and, thus, the proposed method can only be used for LiDAR measurements collected within the atmospheric surface layer (ASL). Specifically for this work, we performed fixed scans with low elevation angle (less than 10 • ) and azimuth angle equal to the mean wind direction to achieve good accuracy in the measurements of the streamwise velocity component, high vertical resolution (≈ 1m), measurements up 465 to the ASL height, and high sampling frequency (between 1 and 3.3 Hz).
For this study, the proposed method for correction of the LiDAR data has been applied to datasets collected during three different field campaigns and for one dataset the procedure has been assessed against simultaneous and co-located sonic anemometer data. For this case, it has been shown that the proposed procedure allows us to correct the second-order statistics of the LiDAR data to estimate a velocity variance comparable to that measured by a sonic anemometer. The compelling 470 results obtained for the correction of the second-order statistics of LiDAR data corroborate the advantage of applying the proposed method, which does not require as input any information of the LiDAR system used, such as probe length and energy distribution over the laser pulse. In contrast to existing methods for correction of LiDAR spatial averaging, all the method parameters are directly estimated from the collected LiDAR data. However, the proposed method can only be applied for LiDAR data collected within the ASL.

475
To better understand the role of the cutoff wavenumber and order of the low-pass filter representing the LiDAR energy damping, further analysis has been conducted on synthetic turbulent velocity spectra. This analysis has been performed by varying mean wind speed, turbulence intensity, and sampling height. This analysis has shown that the main parameter for efficiently correct the LiDAR energy damping is the cutoff wavenumber of the low-pass filter, which is mainly affected by the probe length, while the velocity statistics are weakly affected by the filter order. Furthermore, the results have confirmed that 480 for a given probe length, effects of spatial averaging are enhanced with decreasing wind speed, smaller integral length scale and, thus, for a lower sampling height.
Code availability. The Matlab code for the spectral correction of LiDAR velocity signals is available for free download at https://www. utdallas.edu/windflux/ .
Author contributions. MP contributed to the field campaigns at SLTEST and Celina. GVI is the principal investigator of the LiDAR field 485 campaigns, provided guidance for the research strategy and mentored his Ph.D. student. MP and GVI performed the data analysis, examined the results and prepared the paper.
Competing interests. The authors declare that they have no conflict of interest.