Errors in radial velocity variance from Doppler wind lidar

A high-fidelity lidar turbulence measurement technique relies on accurate estimates of radial velocity variance that are subject to both systematic and random errors determined by the autocorrelation function of radial velocity, the sampling rate, and the sampling duration. Using both statistically simulated and observed data, this paper quantifies the effect of the volumetric averaging in lidar radial velocity measurements on the autocorrelation function and the dependence of the systematic and random errors on the sampling duration. For current-generation scanning lidars and sampling durations of about 30 min and longer, during which the stationarity assumption is valid for atmospheric flows, the systematic error is negligible but the random error exceeds about 10 %. 1 Motivation and approach Coherent Doppler lidars (hereafter called lidars) are increasingly being deployed to measure flow in the atmospheric boundary layer (ABL) particularly for applications to wind engineering (Banta et al., 2013). Accordingly, uncertainties in lidar-derived mean wind velocity estimates have been well characterized (Wang et al., 2016; Lindelöw-Marsden, 2009) and methods and procedures have been developed for error reduction and uncertainty control (Clifton et al., 2013; Gottschall et al., 2012). However, use of lidar for turbulence measurements, while possible (Newman et al., 2016; Branlard et al., 2013; Mann et al., 2010), is less established (Sathe et al., 2015; Sathe and Mann, 2013). Two methods are commonly used to derive the second-order moments (i.e., velocity variances and momentum fluxes) of turbulent flow from lidar data (Sathe and Mann, 2013). The first method initially estimates the three orthogonal wind components from radial velocities acquired through a scan geometry (e.g., the conical scan) and then uses the eddy covariance method to derive the turbulence statistics. The second method involves three steps: (1) obtaining radial velocities to form a time series at different locations through a scan geometry, (2) estimating radial velocity variance from the time series obtained in (1), and (3) deriving the second-order moments based on their relationships with the radial velocity variance. Estimates from both methods have errors relative to the true or expected values defined as the ensemble means over all possible realizations, and these errors are commonly quantified in terms of the systematic and random components. Systematic errors are consistent deviations from the true values in all estimates, and they are also referred to as biases. Random errors are varying deviations from the true values arising for unknown reasons and are commonly modeled as random variables following Gaussian distributions. The first method can suffer from cross-contamination errors (biases) due to the correlation between the radial velocities used for wind velocity estimates (Sathe et al., 2011). Therefore, the second method is considered to be more suitable for lidar turbulence measurements (Newman et al., 2016; Sathe et al., 2015; Sathe and Mann, 2013; Mann et al., 2010) Using the second method mentioned above, errors in the estimated variances and momentum fluxes are accumulations of the following three types of errors in the estimated radial velocity variance: – measurement errors caused by radial velocity estimator uncertainty and atmospheric turbulence (Frehlich, 1997); – attenuation errors due to the volumetric averaging effect of lidar measurement (Sathe and Mann, 2013; Mann et al., 2010); Published by Copernicus Publications on behalf of the European Geosciences Union. 4124 H. Wang et al.: Errors in radial velocity variance from Doppler wind lidar – sampling errors as a result of estimating the ensemble mean by the time average (Mann et al., 2010; Lenschow et al., 1994). Measurement errors can be estimated and potentially removed using either the autocovariance or spectrum of the measured radial velocities (Frehlich, 2001; Lenschow et al., 2000). Alternatively, the signal-to-noise ratio (SNR) can also be used to approximate the measurement error (Pearson and Collier, 1999). Correction for the attenuation error requires knowledge of three dimensional spatial statistics (Mann et al., 2010; Frehlich, 1997), and the accuracy of this correction depends on the suitability of the selected turbulence spectrum model and its parametrization (Sathe and Mann, 2013; Mann et al., 2010). The sampling error is well understood for turbulence statistics estimated from sonic anemometer measurements (Lenschow et al., 1994), but has not yet been studied for the radial velocity variance estimated from lidar measurements. Mann et al. (2010) call for a thorough sampling error analysis in order to understand the difference between radial velocity variance and momentum fluxes derived from lidars and sonic anemometers. The size of sampling error is a function of the sampling interval and duration (Lenschow et al., 1994), both of which are determined by lidar scan geometries which can, in turn, be optimized to minimize the uncertainty in the estimated turbulence statistics (Sathe et al., 2015). Improved understanding of sampling errors in lidar-derived radial velocity variance estimates is a necessary prerequisite to the development of robust techniques that will enable the widespread use of lidar for high-fidelity turbulence measurements. Accordingly, the objectives of this work are to improve the characterization of radial velocity variance error properties and to develop tools to quantify and reduce these errors. The approach taken and the format of this paper are as follows. The theoretical framework used herein to quantify errors in radial velocity variance from lidar measurements leverages the theory developed to characterize uncertainties in statistical moments estimated from a time series of sonic anemometer measurements in Lenschow et al. (1994), and is modified to incorporate the effect of volumetric averaging and the slow sampling rate of lidars (Sect. 3). The theoretical findings are then validated using empirical estimates obtained from measurements of a Galion lidar and three co-deployed sonic anemometers (Sect. 4). The theoretical framework is then used to investigate the sampling duration required to obtain a predefined error magnitude (Sect. 5).


Motivation and approach
Coherent Doppler lidars (hereafter called lidars) are increasingly being deployed to measure flow in the atmospheric boundary layer (ABL) particularly for applications to wind engineering (Banta et al., 2013).Accordingly, uncertainties in lidar-derived mean wind velocity estimates have been well characterized (Wang et al., 2016;Lindelöw-Marsden, 2009) and methods and procedures have been developed for error reduction and uncertainty control (Clifton et al., 2013;Gottschall et al., 2012).However, use of lidar for turbulence measurements, while possible (Newman et al., 2016;Branlard et al., 2013;Mann et al., 2010), is less established (Sathe et al., 2015;Sathe and Mann, 2013).Two methods are commonly used to derive the second-order moments (i.e., velocity variances and momentum fluxes) of turbulent flow from lidar data (Sathe and Mann, 2013).The first method initially estimates the three orthogonal wind components from radial velocities acquired through a scan geometry (e.g., the conical scan) and then uses the eddy covariance method to derive the turbulence statistics.The second method involves three steps: (1) obtaining radial velocities to form a time series at different locations through a scan geometry, (2) estimating radial velocity variance from the time series obtained in (1), and (3) deriving the second-order moments based on their relationships with the radial velocity variance.Estimates from both methods have errors relative to the true or expected values defined as the ensemble means over all possible realizations, and these errors are commonly quantified in terms of the systematic and random components.Systematic errors are consistent deviations from the true values in all estimates, and they are also referred to as biases.Random errors are varying deviations from the true values arising for unknown reasons and are commonly modeled as random variables following Gaussian distributions.The first method can suffer from cross-contamination errors (biases) due to the correlation between the radial velocities used for wind velocity estimates (Sathe et al., 2011).Therefore, the second method is considered to be more suitable for lidar turbulence measurements (Newman et al., 2016;Sathe et al., 2015;Sathe and Mann, 2013;Mann et al., 2010) Using the second method mentioned above, errors in the estimated variances and momentum fluxes are accumulations of the following three types of errors in the estimated radial velocity variance: measurement errors caused by radial velocity estimator uncertainty and atmospheric turbulence (Frehlich, 1997); attenuation errors due to the volumetric averaging effect of lidar measurement (Sathe and Mann, 2013;Mann et al., 2010); Published by Copernicus Publications on behalf of the European Geosciences Union.
H. Wang et al.: Errors in radial velocity variance from Doppler wind lidar sampling errors as a result of estimating the ensemble mean by the time average (Mann et al., 2010;Lenschow et al., 1994).
Measurement errors can be estimated and potentially removed using either the autocovariance or spectrum of the measured radial velocities (Frehlich, 2001;Lenschow et al., 2000).Alternatively, the signal-to-noise ratio (SNR) can also be used to approximate the measurement error (Pearson and Collier, 1999).Correction for the attenuation error requires knowledge of three dimensional spatial statistics (Mann et al., 2010;Frehlich, 1997), and the accuracy of this correction depends on the suitability of the selected turbulence spectrum model and its parametrization (Sathe and Mann, 2013;Mann et al., 2010).The sampling error is well understood for turbulence statistics estimated from sonic anemometer measurements (Lenschow et al., 1994), but has not yet been studied for the radial velocity variance estimated from lidar measurements.Mann et al. (2010) call for a thorough sampling error analysis in order to understand the difference between radial velocity variance and momentum fluxes derived from lidars and sonic anemometers.The size of sampling error is a function of the sampling interval and duration (Lenschow et al., 1994), both of which are determined by lidar scan geometries which can, in turn, be optimized to minimize the uncertainty in the estimated turbulence statistics (Sathe et al., 2015).Improved understanding of sampling errors in lidar-derived radial velocity variance estimates is a necessary prerequisite to the development of robust techniques that will enable the widespread use of lidar for high-fidelity turbulence measurements.Accordingly, the objectives of this work are to improve the characterization of radial velocity variance error properties and to develop tools to quantify and reduce these errors.
The approach taken and the format of this paper are as follows.The theoretical framework used herein to quantify errors in radial velocity variance from lidar measurements leverages the theory developed to characterize uncertainties in statistical moments estimated from a time series of sonic anemometer measurements in Lenschow et al. (1994), and is modified to incorporate the effect of volumetric averaging and the slow sampling rate of lidars (Sect.3).The theoretical findings are then validated using empirical estimates obtained from measurements of a Galion lidar and three co-deployed sonic anemometers (Sect.4).The theoretical framework is then used to investigate the sampling duration required to obtain a predefined error magnitude (Sect.5).

Preliminaries
A brief description of lidar measurements is given below.For more details see Mann et al. (2008) and Sathe and Mann (2012).Symbols introduced hereafter are explained in the nomenclature in Appendix B.
x 2 Based on the streamline coordinate system in Fig. 1, a wind velocity vector is defined as where u 1 , u 2 , and u 3 are streamwise, transverse, and vertical wind components respectively at position x = (x 1 , x 2 , x 3 ).Without loss of generality, we can assume that the mean streamwise velocity U 1 has been removed, and u only consists of the fluctuating components of turbulence with zero means.A lidar measures the radial velocity (v r ) from the Doppler frequency shift induced by the motion of scatterers along the line of sight (LOS), where the orientation of LOS is defined by the unit directional vector: where θ is the azimuth angle that increases clockwise from being zero in a positive x 2 direction and φ is the elevation angle relative to the x 1 -x 2 plane (Fig. 1).The radial velocity is the projection of wind velocity on the LOS and is defined as The Galion lidar used here is a pulsed lidar that measures radial velocity with the step-stare technique in this application.Each radial velocity is measured over a dwell time of approximately 1.0 s, during which spectra of a large number of returned signals are averaged to improve the measurement accuracy.When operated with a scan geometry, the lidar steers its transceiver to probe at different sets of θ and φ.Hence, the sampling interval of two consecutive measurements at one location depends on the dwell time, the scan geometry, and the mechanical design of the lidar.The shortest sampling interval can be achieved with the staring scan for which the lidar measures with fixed θ and φ; that is, the sampling interval is close to the dwell time.The radial velocity (v R ) is estimated from an averaged spectrum acquired over a range gate; hence, it is a weighted average of radial velocities along the LOS: where s denotes the range gate location and s is the range distance on the LOS.Here we use the subscript R to denote an average quantity and r to denote a point quantity.Note that v R defined in Eq. ( 4) is the expected value of a radial velocity measurement.The actual measured value differs from v R due to measurement errors as discussed in Frehlich (1997) and Sect. 1.The weighting function Q in Eq. ( 4) can be approximated by a Gaussian function (Kristensen et al., 2011): where the standard deviation σ Q is a measure of the volumetric averaging size, which is 15.4 m for the Galion lidar used herein (Wang et al., 2016).The covariance of v r (R r ) and v R (R R ) along the x 1 direction is defined as follows (Mann et al., 2008): where s and s denote the range distance, R ij is the velocity tensor for the ith and j th velocity components, n i and n j are the components of n, and e 1 and r 1 are the unit vector and the separation distance on the x 1 axis, respectively.Note that i, j = 1, 2, 3 and summation is assumed over repeated indices.We can map r 1 to the temporal lag τ through the frozen turbulence hypothesis (Taylor, 1938); r 1 = U 1 τ .The spatial autocorrelation function of v r (ρ r ) and v R (ρ R ) separated by (r 1 , 0, 0) is then defined respectively as where µ 2,r = R r (0) is the ensemble variance of v r , and µ 2,R = R R (0) is the ensemble variance of v R .Due to the averaging given in Eq. ( 4), R R (0) < R r (0) and the ratio R R (0)/R r (0) decreases when the size of volumetric averaging increases (i.e., σ Q /L 1 increases, where L 1 is the integral length scale of u 1 ) (Mann et al., 2008).However, R R (r 1 ) = R r (r 1 ) when r 1 is sufficiently large (e.g., r 1 L 1 ) because values of R r (r 1 ) with large r 1 are determined by eddies of large sizes that are not affected by the averaging.As 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Examples of turbulence statistics of point radial velocity (v r ) and averaged radial velocity (v R ) derived from Eqs. ( 6) and (7) using the isotropic turbulence model (Pope, 2000) and the von Kármán spectra (Burton et al., 2011) with streamwise velocity U 1 = 8 m s −1 , turbulence intensity = 12.5 %, and lidar elevation angle = 10 • .The covariance and autocorrelation functions for v r (R r and ρ r ) and v R (R R and ρ R ) for azimuth angle θ = 30 • are shown in (a) and (b), respectively, as functions of r 1 /L 1 , where r 1 is the spatial lag in streamwise direction and L 1 is the integral length scale of streamwise velocity.R R and ρ R are presented in terms of σ Q /L 1 , where σ Q represents the size of volumetric averaging (see Eq. 5).
The effect of volumetric averaging on the integral length scale of radial velocity is presented in (c) in terms of the relationship between L R /L r and σ Q /L 1 for different θ values, where L r and L R are the integral length scales for v r and v R , respectively.illustrated in Fig. 2, with the assumption of isotropic turbulence, R R (r 1 ) starts at a value lower than R r (r 1 ) at r 1 = 0 and gradually converges to R r (r 1 ) as r 1 increases.Because R r (0) and R R (0) are the denominators in Eqs. ( 8) and (9), respectively, ρ R ≥ ρ r and (ρ R − ρ r ) increases with increasing As a result, if we denote L r and L R respectively as the integral length scales for v r and v R , L R > L r and L R /L r increases with increasing σ Q /L 1 (Fig. 2).Although quantitative evaluation of R R /R r or (ρ R − ρ r ) as a function of r 1 requires the knowledge of velocity tensors R ij , the qualitative H. Wang et al.: Errors in radial velocity variance from Doppler wind lidar statement above is also true for non-isotropic turbulence and has implications for the errors of radial velocity variance estimates as demonstrated in the next section.

Errors in radial velocity variance
In the following, the analysis and notation used are from Lenschow et al. (1994), and we use the word "error" to refer to the sampling error of radial velocity variance estimated from time series unless otherwise stated.
Radial velocity variance is estimated from time series of radial velocity that is related to u(t) which is a stationary and ergodic time series generated from a Gaussian process characterized by the mean wind speed U 1 and covariance tensors.Assuming that the mean has been removed from u(t), it can be shown that both v r and v R are also from a stationary and ergodic time series of a Gaussian process that has the following properties: zero ensemble mean for both v r and v R ; ensemble autocorrelation ρ r (τ ) and ρ R (τ ).
The brackets • denote the ensemble average.The measured radial velocity ( vR ) from a lidar will deviate from the expected value (v R ) as a result of the measurement error (Frehlich, 2001).The bias introduced by the measurement error is negligible for a lidar measurement with a high pulse accumulation number (e.g., 20 000 for the Galion lidar) (Frehlich, 2001(Frehlich, , 1997)).The random measurement error in vR is modeled by an independent Gaussian random variable with zero mean and variance σ 2 R , which is inversely proportional to the SNR (Pearson and Collier, 1999;Frehlich and Yadlowsky, 1994).Turbulence can cause σ 2 R to increase by an amount that is proportional to the range gate size and the turbulent kinetic energy (TKE) dissipation rate (Frehlich, 1997).Frehlich et al. (1994) show that the amount of increase due to turbulence is about 20 % for a 48 m range gate size.For the Galion lidar with 30 m range gate size, the contribution from turbulence to σ 2 R is no more than 10 % with a relative weak SNR and a TKE dissipation rate on the order of 0.01 m 2 s −3 (based on the method in Frehlich, 1997).Therefore, it is a good approximation to estimate σ 2 R from the SNR data using Eq. ( 5) in Pearson and Collier (1999) for the Galion lidar (Wang et al., 2016).
Assuming zero bias in the measurement error, it can be shown that the radial velocity variance ( μ2,R (T )) estimated from lidar measurements over a sampling period of T can be represented as where µ 2,R (T ) is the estimate of µ 2,R from a time series over a period T without measurement errors.According to Lenschow et al. (1994), where r,R (T ) is random sampling error that is characterized by a Gaussian distribution with zero mean and variance denoted as E 2 r,R (T ).
Both systematic and random errors are results of estimating ensemble means by time averages over a finite time period; therefore, E s,R (T ) and E 2 r,R (T ) must be functions of T .Note that σ 2 R in Eq. ( 10) is independent of T and can be estimated and removed using the SNR data as discussed above.Hence, we assume no measurement error in the following derivation of E s,R (T ) and E 2 r,R (T ) from a disjunct time series of v R .The derived results for v R can also be used for v r .As discussed in more detail in the subsequent text, the difference between v r and v R in terms of the relative sampling errors derives solely from the difference between ρ r and ρ R .
Radial velocities from a pulsed lidar naturally form a discrete time series (see Sect. 2).For example, when a pulsed lidar is configured to estimate turbulence statistics with the six-beam method of Sathe et al. (2015), the lidar samples six locations sequentially, and therefore the sampling interval (δt) at one location is 6 times the sampling interval per one lidar measurement (i.e., δt > 6 s and is at least 60 times slower than a sonic anemometer sampling at 10 Hz).Hence, the radial velocity variance from a time series of v R acquired over a period T with a sampling interval δt is estimated from where t i is the time stamp of measurement, N = 1 + T /δt is the sample number, and µ 1,R (T ) is the ensemble mean estimate: The radial velocity variance estimated from Eq. ( 12) has a systematic error E s,R given by Note that µ 2 1,R (T ) is the sample mean variance; hence the systematic error is always negative, and its relative magnitude is given by Box et al. (2015): where The variance of the random error r,R (T ) of radial velocity variance estimated from Eq. ( 12) is defined as Applying the Isserlis relation of a Gaussian process (Lenschow et al., 1994), it can be shown that the relative sampling random error variance is where The error expressions derived in Eqs. ( 15) and ( 18) are valid for v r variance with the replacement of the autocorrelation and the variance by ρ r and µ 2,r , respectively.It is clear from Eqs. ( 15) and ( 18) that the errors are functions of the sampling interval, sampling duration, and autocorrelation function.The averaging defined in Eq. ( 4) affects the shapes of the radial velocity covariance and autocorrelation function as illustrated using the isotropic turbulence model (Fig. 2).To investigate the impact of anisotropic turbulence under realistic atmospheric conditions, wind vectors are statistically simulated from the Risø Smooth-Terrain (SMOOTH) spectrum model using TurbSim (https://nwtc.nrel.gov/TurbSim).The simulation domain is 725×60×12 m (x 1 × x 2 × x 3 ) centered at 80 m above the ground, and has a horizontal and vertical resolution of 1.0 and 0.5 m, respectively.In these simulations, the mean wind speed is 8 m s −1 and the time interval is 0.125 s.Point radial velocities (v r ) are calculated using Eq.(3) for varying LOS orientations in the horizontal plane relative to the x 1 direction (denoted by β in Fig. 1) with a fixed elevation angle of 10 • .Lidar-equivalent radial velocities (v R ) are derived by averaging v r on the LOS within ±30 m from the domain center using the weighting function in Eq. ( 5).As expected, variance reduction occurs for all LOS orientations.The difference between R R and R r in Fig. 3 is similar to that in Fig. 2, and can be explained using the structure function (D r ) defined as for a given time lag (τ ) and used to represent the energy of eddies of sizes that are smaller than the scale U 1 τ (Pope, 2000).Volumetric averaging only attenuates the energy of eddies of small sizes (Mann et al., 2008).When τ is small relative to the integral timescale τ 0 , volumetric averaging causes D r (τ ) to decrease; hence, R R decreases more slowly than R r with respect to τ per Eq. ( 21).When τ/τ 0 is large (e.g., τ/τ 0 > 0.25), R R = R r because volumetric averaging has little effect on eddies of scales of large τ .The difference between R R and R r results in ρ R > ρ r for all time lags and LOS orientations (Fig. 3).Simulations conducted with the Kaimal (IECKAI) spectrum model TurbSim (https: //nwtc.nrel.gov/TurbSim)produce similar results.Errors associated with the autocorrelation functions ρ r /ρ R derived from the simulated v r /v R from both models are shown in Fig. 4.Although the difference between the errors of variance of v r and v R varies with LOS orientation and turbulence model (i.e. the turbulence structure), errors associated with v R are consistently higher than those related to v r , indicating that volumetric averaging increases errors associated with radial velocity variance estimates.Thus both the statistically simulated wind data and physical reasoning provide evidence that volumetric averaging increases the autocorrelation of radial velocity and inflates the errors in radial velocity variance estimates.Further confirmation will be provided in the next section, where data from a www.atmos-meas-tech.net/9/4123/2016/Atmos.Meas.Tech., 9, 4123-4139, 2016  (Barthelmie et al., 2016).During the experiment, a Galion lidar was configured with 20 kHz pulse repetition frequency and 1.0 s dwell time to scan at four elevation angles (4.8, 10.0, 15.2, and 20.6 • ) with a fixed azimuth angle of 349 • such that the seventh range gate of the lidar sampled at 20, 40, 60 m (and 80 m) above the ground where three Gill Windmaster Pro sonic anemometers were installed on a slender meteorological mast and sampled at 10 Hz.The sampling interval of the lidar at each elevation angle is about 7.5 s, which is similar to the sampling interval of the six-beam technique used for turbulence measurement in Sathe et al. (2015).The measurements from the sonic anemometers are used to describe the atmospheric turbulence conditions and evaluate the accuracy of the lidar measurements.
The lidar conducted automatic cleaning at the beginning of each hour, resulting in a 60 s gap in the measurements; thus analysis presented here uses a sample period of one hour.Measurement errors were corrected in lidar radial velocity variance estimates using the SNR data according to Eq. ( 5) in Pearson and Collier (1999).Radial velocity variance from the sonic data is estimated by extrapolating the autocovariance of the first 10 time lags to lag zero, assuming that autocovariance is linearly proportional to τ 2/3 for the first 10 lags (Lenschow et al., 2000).Comparison of the hourly mean and variance of radial velocities from the lidar and sonic anemometers indicates good agreement (correlation coefficient = 0.998) with the exception of periods around 18 May 2015 23:00 UTC when the measurements were in the wake of a wind turbine located 60 m southwest of the Galion (Fig. 5).
The variance is consistently higher (on average by 19 %) for the sonic radial velocity than the lidar radial velocity because of the expected attenuation in variance caused by volumetric averaging.
Stationarity is the fundamental assumption required to obtain theoretical estimates of the errors (as in Sect.3).Thus, the hourly time series of radial velocity from both the lidar and sonic anemometers are evaluated for stationarity using the approach of Foken and Wichura (1996).Each hourly time series is evenly divided into 12 subsets.If the mean of the variance of the subsets deviates by less than 30 % from the variance of the hourly time series, the time series is considered to be stationary.Among all the hourly time series obtained at the three heights, 34 concurrent pairs of lidar and sonic data pass the stationarity test (Fig. 5); therefore, they are used to derive the empirical estimates of errors.

Error estimation method
The systematic error e s (T n , T ) and random error variance e 2 r (T n , T ) associated with sampling duration T n derived from a time series of length T are estimated using the following stationary bootstrap method (Politis and Romano, 1994), where the sample numbers associated with T and T n are denoted as N and N n , respectively.
-A new time series of size of N is constructed by resampling blocks of the original time series.To keep the new time series stationary, the sizes of the blocks are randomly drawn from a geometric distribution, and the locations (the start of each block) are randomly drawn from a discrete uniform distribution on (1, 2, . .., N ).
The only parameter to be specified is the optimal mean block size for the geometric distribution, which is found by minimizing the mean squared error of the estimate of the sample mean variance (Politis and White, 2004).
-A subset of size of N n is randomly selected from the new time series and its mean and variance, denoted as µ 1 (T n,i ) and µ 2 (T n,i ), respectively, are recorded.The subscript i denotes that it is the ith resampling.
-Sequences of µ 1 (T n,i ) and µ 2 (T n,i ) are acquired after repeating the two steps above for N b times.The variance of the sample mean µ 2 1 (T n , T ) is approximated as -Per the definition in Eq. ( 14), the systematic error e s (T n , T ) can be calculated as -To calculate the random error variance, the expected value µ 2 (T n , T ) is first estimated by adding the systematic error to µ 2 , which is approximated by µ 2 (T ), i.e., Then the random error variance is derived with the following equation:

Observed errors
Two methods are used to estimate the errors associated with different sampling durations after the means are removed from hourly time series of the point radial velocity from sonic anemometers (v r,sonic ) and lidar-measured averaged radial velocity (v R,lidar ).The first method, denoted as the M ρ method, is based on Eqs. ( 15) and ( 18) and the autocorrelation function derived from measurements (Fig. 6a).The observed autocorrelation functions show the postulated effect of volumetric averaging on the autocorrelation function.For all the hourly time series studied here, the autocorrelation of v R,lidar at time lag one (τ = 7.5 s) is significantly higher than that of v r,sonic (Fig. 7).At time lag two (τ = 15 s), the autocorrelation of v R,lidar is still larger than that of v r,sonic , but the difference is not always statistically significant.Beyond time lag two (τ > 15 s), the difference between the autocorrelation of v R,lidar and v r,sonic vanishes.Because integral timescales of streamwise velocity calculated from the sonic data are all below 30 s, we assume that the autocorrelation function values of both v r,sonic and v R,lidar are zero for time lags larger than 60 s.The autocorrelation-based systematic error and random error variance will be denoted as e s,ρ and e 2 r,ρ , respectively.The second method, denoted as the M b method, uses the stationary bootstrap method described in Sect.4.2, and the resultant systematic error and random error variance will be denoted as e s,b and e 2 r,b , respectively (see the example given in Fig. 6).
Consistent with expectations, error estimates from both M ρ and M b methods are higher for v R,lidar than v r,sonic (Fig. 8) due to the difference in autocorrelation functions.The random errors estimated from both methods in general decrease with height except for a few cases related to the M ρ method (Fig. 9).The systematic errors show a similar trend with height (not shown here).This is possibly the result of increasing turbulence integral length scale with height.For a pulsed lidar, increasing integral length scale can cause the amount of attenuation of radial velocity variance to decrease (Mann et al., 2010) and consequently the errors to decrease (as demonstrated in Fig. 2).Both methods give very similar estimates of systematic error, although there are some hours for which the M ρ method produces higher systematic errors than the M b method (Fig. 10).The median of e s from v R,lidar is 1.5/0.9% for T = 30/55 min.Estimates of random error from the M ρ method are always higher than the M b method, due in part to the negative bias in the M b method (Politis and White, 2004) (Fig. 10).For v R,lidar , the median of e 2 r from v R,lidar is 1.6/0.9% and the median of e r (standard deviation) is 12.7/9.5 % for T = 30/55 min, respectively.Despite the expected difference between errors from the two methods, both methods consistently confirm the trend of decreasing error with increasing sampling duration (Fig. 11), and the relatively close agreement of results from the two approaches offers empirical support for the relationships between the er- rors and the autocorrelation function of radial velocity as described by Eqs. ( 15) and ( 18) for ergodic and stationary time series.

Discussions
On the basis of empirical evidence presented in Sect.4, in the following we use the theoretical framework presented in Sect. 3 to describe how the error in estimating radial velocity variance from lidar measurement that results from (i) autocorrelation function, (ii) sampling duration, and (iii) sampling interval can be minimized.The first factor (autocorrelation function) is determined by the underlying wind field, the lidar LOS orientation, and the size of volumetric averaging (i.e., σ Q in Eq. 5), and needs to be specified to estimate uncertainty in the variance estimates.It typically approximates an exponential but the precise functional form and the rate of decay varies depending on the flow field.Therefore, in practice, error reduction can only be achieved by adjusting the other two factors: sampling duration (T ) and sampling inter- .The relative variance of random error e 2 r for the estimated radial velocity variance from the observational data at three heights ordered by their occurrence time (x axis).Errors in (a) and (c) are estimated using the autocorrelation method (M ρ ) for 30 and 55 min of sampling duration, respectively.Errors in (b) and (d) are estimated using the bootstrap method (M ρ ) for 30 and 55 min of sampling duration, respectively.r and the sampling duration for radial velocity variance estimated from the lidar data using the autocorrelation method (gray lines) and the stationary bootstrap method (box plots).All errors are normalized by the errors estimated from the autocorrelation method with a sampling duration of 55 min. .Contours of the relative variance e 2 r of random errors of radial velocity variance from lidar measurements estimated with the isotropic turbulence model (Pope, 2000) and the von Kármán spectra (Burton et al., 2011) as a function of the sampling duration (T ) and the sampling interval (δt) normalized by the integral timescale (τ 0 ) for four different β values, where β is the angle between the LOS and the wind direction.The weighting function (Eq.5) representing the volumetric averaging has a standard deviation σ Q = 0.2 L 1 , where L 1 is the streamwise integral length scale.The other input parameters include the elevation angle φ = 10 • and the mean wind speed U 1 = 8 m s −1 .val (δt), with assumption or knowledge of the autocorrelation function ρ(τ ).
Atmospheric turbulence is rarely isotropic, and for all the hours presented in Sect. 4 the three wind components had non-equal variance.However, here we use the isotropic turbulence model to model the autocorrelation function, noting that it is always possible to find an integral length scale to reproduce the observed autocorrelation function of radial velocity from lidar measurements with the isotropic turbulence model (Pope, 2000) and the von Kármán spectra (Burton et al., 2011) (e.g., Fig. 6).Therefore, we argue that with a proper length scale, the isotropic turbulence model can generate an autocorrelation function that can be used to approximate non-isotropic turbulence conditions, justifying the use of the isotropic turbulence model here.
Based on the isotropic turbulence model and the von Kármán spectra, the relative systematic error (e s ) is negligible in comparison to the random error.Hence, only the analysis of random error variance e 2 r or standard deviation (e r ) is given below.The magnitude of e 2 r is not sensitive to the sam-pling interval (δt).Because the integral timescale (τ 0 ) has the order of magnitude of 10 s and δt = 1-10 s, it is likely that δt/τ 0 < 1, implying that in practice the sampling interval is a minor factor on error reduction (Fig. 12).Increasing the size of volumetric averaging in terms of σ Q can cause e 2 r to increase with a rate that decreases to nearly zero when the sampling duration increases and the LOS moves from being parallel to perpendicular to the wind direction (i.e., β from 0 to 90 • ) (Fig. 13).Therefore, it is also a minor factor on error reduction when the sampling duration is long.The LOS orientation relative to the wind direction (β) naturally affects the properties of random errors because the timescale of radial velocity varies with β.For example, for the isotropic turbulence model used here, streamwise velocity has the largest timescale and as a result, errors are large when the LOS is aligned with the wind direction (i.e., β = 0 • ) (Figs. 12a and  13a).In general e 2 r is not sensitive to β, but it must be acknowledged that the effect of β on e 2 r will change if a different turbulence model is applied.The factor that has the most significant effect on error reduction is the sampling du-  (Pope, 2000) and the von Kármán spectra (Burton et al., 2011) as a function of the sampling duration (T ) normalized by the integral timescale (τ 0 ) and the volumetric averaging size (σ Q ) normalized by the streamwise integral length scale (L 1 ) for four different β values, where β is the angle between the LOS and the wind direction.The other input parameters include the sampling interval δt = 0.5τ 0 , the elevation angle φ = 10 • , and the mean wind speed U 1 = 8 m s −1 .ration (Figs. 12 and 13).Thus optimization of lidar operation for retrieval of radial velocity variance can be considered through the lens of "how long is long enough?" (Lenschow et al., 1994).Provided that the optimum six-beam configuration proposed in Sathe et al. (2015) (i.e., φ = 45 • ) is applied to the Galion lidar for which σ Q = 15.4 m under neutral atmospheric conditions over flat terrain, both the systematic and random errors are almost independent of LOS orientation and surface roughness length, which is used to predict the integral length scale and turbulence intensity (Fig. 14).The systematic error is lower than 1.0 % when T > 30 min (Fig. 14).The standard deviation (e r ) of random errors can be reduced from 12.0 to 9.0 % by increasing the sampling duration from 30 to 60 min (Fig. 14), which is consistent with the observed standard deviation 12.7-9.5 % in Sect.4.3.The standard deviation remains higher than 6.0 % when T increases to 120 min (Fig. 14).Note that the volumetric averaging increases e r by less than 1.0 % because of its relatively small size (σ Q /L 1 < 0.1).The implication of the analysis above is that, given that 0.5-1.0h is usually the length over which the stationarity assumption is valid in the ABL (Larsén et al., 2016), the random error in radial velocity variance estimates will likely be around 10.0 % and it will be difficult to estimate radial velocity variance with random errors lower than 5.0 %.
The relatively high sampling error discussed above, in combination with measurement and attenuation errors described in Sect. 1, will propagate through and result in a bias and random error in the estimated turbulence parameters.The bias is mainly due to the attenuation error from volumetric averaging.The contribution from the measurement error to the bias can be ignored given that the SNR and the number of spectra used for radial velocity estimation are high (i.e., the dwell time is sufficiently long).The random error in the estimated turbulence parameters is mainly from the random sampling error.Thus, we can consider only attenuation and sampling errors in error quantification and scan geometry optimization for lidar turbulence measurements.The error quantification method is given in Appendix A for lidar-derived second-order moments (i.e., velocity variances and momentum fluxes), and is applied to evaluate the uncertainty in the second-order moments esti- r with respect to the sampling duration (T ) in (a) and (b), respectively, predicted from the isotropic turbulence model (Pope, 2000) and the von Kármán spectra (Burton et al., 2011) at 80 m height under neutral conditions for surface roughness lengths (z 0 ) from 10 −3 to 10 0 m and LOS orientations β from 0 to 90 • , where β is the angle between the LOS and the wind direction.The blue lines with squares and the dark lines are the mean error values from the entire range of z 0 and β for the averaged and point radial velocity variance, respectively.The red shaded areas denote the range of errors of the averaged radial variance.The equation used to estimate the integral length scale can be found in Wang et al. (2016) with a mean wind speed of 7 m s −1 and Coriolis parameter of 10 −4 s −1 .The elevation angle φ = 45 • .mated with the optimal scanning geometry of the six-beam method from Sathe et al. (2015) using the turbulence conditions, attenuation errors, and sampling errors observed during the PEIWEE experiment.The derived uncertainties of both the streamwise velocity variance and the vertical momentum flux, in terms of the relative standard deviation, decrease with increasing sampling duration (Fig. 15).The uncertainty values are similar for the streamwise velocity variance derived from the six-beam method and the sonic data (Fig. 15); they are around 10 % for T = 30 min and 8 % for T = 55 min.However, the vertical momentum flux derived from the six-beam method is much higher than that from the sonic data (Fig. 15).For a 30 (55) min sampling duration, the median uncertainty is 38 (28) % for the vertical momentum flux derived from the six-beam method and 11 (8) % from the sonic data.Such large uncertainties in the derived turbulence statistics can introduce uncertainty when verifying lidar turbulence measurements with sonic data.Using the median uncertainty values in Fig. 15 and the linear model for verification, statistical simulations show that the coeffi- Relative standard deviations for the sixbeam method are calculated using the method in Appendix A and the observed attenuation and sampling errors in Sect.4, and those for sonic data are based on the derived turbulence statistics from the sonic data described in Sect.4.1 and error quantification methods given in Lenschow et al. (1994).The sampling duration is denoted by the colors of box plots (see the legend).
cients of determination (i.e., R 2 ) for the vertical momentum flux are 0.552 ± 0.073 for a 30 min sampling duration and 0.694 ± 0.065 for a 55 min sampling duration.The uncertainty is much lower for the streamwise velocity variance.
The R 2 values are 0.905 ± 0.020 for a 30 min sampling duration and 0.946 ± 0.012 for a 55 min sampling duration.

Concluding remarks
Use of lidar for estimation of turbulence fields if realized could revolutionize atmospheric boundary layer characterization studies and has applications to many fields.Accurate radial velocity variance estimates are necessary (but not sufficient) to obtain robust turbulence statistics from lidar.The accuracy of radial velocity variance estimates and their relationship to pseudo-point measurements from sonic anemometers are determined by (i) the applicability of the stationarity assumption, (ii) the effect of volumetric averaging on the radial velocity autocorrelation function, (iii) the sampling interval, and (iv) the sampling duration.Of these factors, (i), the stationarity assumption, is determined only by atmospheric conditions, but it is most likely to be achieved within the period of 1 h in environments where the surface conditions are homogeneous.The second factor, (ii), the volumetric averaging, is dictated by the probe length that is determined by the lidar properties; it causes the radial velocity autocorrelation function to increase, and thus increases errors in radial velocity variance estimates.Large probe length can result in high errors.The third factor, (iii), the sampling interval, is determined partly by the scan geometry which is needed to sample radial velocities with different LOS orientations to reconstruct the wind field, and partly by the lidar configurations of e.g., the dwell time of each measurement and the scanning speed.Errors are not sensitive to the sampling interval because the sampling interval for lidar turbulence measurement is commonly smaller than the turbulence integral timescale.The last factor, (iv), the sampling duration, which together with the sampling interval determines the number of samples available for radial velocity estimates, can only be extended to the limit implied by the stationarity assumption, but in principle, as sampling duration increases, the errors associated with the radial velocity variance decrease.
Given these constraints on radial velocity variance estimates, this paper uses theories and empirical observations to show that for sample periods for which stationarity can reasonably be asserted (approximately 1 h), the systematic error can be reduced to a level lower than 1 %, and the standard deviation of random errors will be around 10 %.These errors will propagate through to estimation of turbulence statistics from lidar measurements, and thus provide a fundamental limit on the likely accuracy of those estimates.

Figure 1 .
Figure 1.The schematic of lidar line of sight (LOS) orientation in the streamline coordinate system

Figure 3 .
Figure 3. Covariance (R r and R R on the left) and autocorrelation functions (ρ r and ρ R on the right) of the point radial velocity (v r ) and averaged radial velocity (v R ) derived from statistically simulated time series using TurbSim with the SMOOTH model and the default parameter values.Covariance and autocorrelation functions are presented for four different LOS orientations here as indicated by the angle of LOS relative to the mean wind direction (β).The time lag (τ ) is normalized by τ 0 , which is the first time when ρ r = 1/e at β = 0.A fixed elevation angle of 10 • is used.

Figure 4 .
Figure 4.The relative systematic error (e s ) and random error variance e 2 r derived from the autocorrelation functions of the simulated point radial velocity (ρ r ) and averaged radial velocity (ρ R ) from the SMOOTH model in (a) and (c) and from the IECKAI model in (b) and (d) (TurbSim) as functions of the angle between the LOS and mean wind direction (β).The sampling duration is 100τ 0 , where τ 0 is the first time when the point radial velocity autocorrelation crosses 1/e at β = 0 • .

Figure 5 .
Figure5.Time series of (a) hourly mean and (b) variance of radial velocity (v R ) from the lidar (markers) and the sonic anemometers (lines) and time series of (c) hourly mean wind speed and (d) direction from the sonic anemometers at three different heights.The mean and variance of radial velocity from the hourly time series classified as stationary using the method ofFoken and Wichura (1996) are shown by the filled markers.The dashed line in (d) gives the wind direction under which the sonic anemometers are in the center of the wind turbine wake.

Figure 6 .
Figure 6.An example of the autocorrelation and errors in radial velocity variance estimate using data from the hour starting at 19 May 2015 04:00 at 40 m height.The autocorrelation functions derived from the lidar and sonic measurements are presented in (a).Systematic errors (e s ) and random errors e 2 r are presented in (b) and (c), respectively.Errors associated with both lidar and sonic anemometer data are estimated using the autocorrelation function from measurements (denoted by M ρ ) and the stationary bootstrap method (denote by M b ) described in Sect.4.2.

Figure 7 .
Figure 7.Comparison of the values of radial velocity autocorrelation function from sonic data (ρ sonic ) and lidar data (ρ lidar ) at the first time lag δt = 7.5 s.The 95 % confidence interval of ρ lidar is indicated by the error bar.

Figure 8 .
Figure8.Comparison of systematic errors estimated from sonic data (e s,sonic ) and from lidar data (e s,lidar ) in (a) for a sampling duration of 30 min and (c) a sampling duration of 55 min, and comparison of the random error estimated from sonic data e 2 r,sonic and from lidar data e 2 r,lidar in (b) for a sampling duration of 30 min and (d) a sampling duration of 55 min.The method used to estimate the errors is indicated by the subscript in the legend, where ρ denotes the M ρ method using the autocorrelation function and b denotes the M b method using the stationary bootstrap method.The solid lines are the 1 : 1 lines.
Figure9.The relative variance of random error e 2 r for the estimated radial velocity variance from the observational data at three heights ordered by their occurrence time (x axis).Errors in (a) and (c) are estimated using the autocorrelation method (M ρ ) for 30 and 55 min of sampling duration, respectively.Errors in (b) and (d) are estimated using the bootstrap method (M ρ ) for 30 and 55 min of sampling duration, respectively.

Figure 10 .
Figure 10.Comparison of lidar radial velocity systematic errors e s,ρ and e s,b estimated using the autocorrelation (M ρ ) method and the stationary bootstrap (M b ) methods, respectively, for (a) a sampling duration of 30 min and (c) a sampling duration of 55 min, and comparison of the lidar radial velocity random errors e 2 r,ρ and e 2 r,b estimated using the M ρ method and the M b method, respectively, for (b) a sampling duration of 30 min and (d) a sampling duration of 55 min.The measurement heights are shown in the legend.

Figure 11 .
Figure 11.Relationships between (a) the systematic error (e s ) and the sampling duration and (b) the random error e 2r and the sampling duration for radial velocity variance estimated from the lidar data using the autocorrelation method (gray lines) and the stationary bootstrap method (box plots).All errors are normalized by the errors estimated from the autocorrelation method with a sampling duration of 55 min.
Figure12.Contours of the relative variance e 2 r of random errors of radial velocity variance from lidar measurements estimated with the isotropic turbulence model(Pope, 2000) and the von Kármán spectra(Burton et al., 2011) as a function of the sampling duration (T ) and the sampling interval (δt) normalized by the integral timescale (τ 0 ) for four different β values, where β is the angle between the LOS and the wind direction.The weighting function (Eq.5) representing the volumetric averaging has a standard deviation σ Q = 0.2 L 1 , where L 1 is the streamwise integral length scale.The other input parameters include the elevation angle φ = 10 • and the mean wind speed U 1 = 8 m s −1 .

Figure 13 .
Figure13.Contours of the relative variance e 2 r of random errors of radial velocity variance from lidar measurements estimated with the isotropic turbulence model(Pope, 2000) and the von Kármán spectra(Burton et al., 2011) as a function of the sampling duration (T ) normalized by the integral timescale (τ 0 ) and the volumetric averaging size (σ Q ) normalized by the streamwise integral length scale (L 1 ) for four different β values, where β is the angle between the LOS and the wind direction.The other input parameters include the sampling interval δt = 0.5τ 0 , the elevation angle φ = 10 • , and the mean wind speed U 1 = 8 m s −1 .

Figure 14 .
Figure14.Variation of the systematic error (e s ) and variance of random error e 2 r with respect to the sampling duration (T ) in (a) and (b), respectively, predicted from the isotropic turbulence model(Pope, 2000) and the von Kármán spectra(Burton et al., 2011) at 80 m height under neutral conditions for surface roughness lengths (z 0 ) from 10 −3 to 10 0 m and LOS orientations β from 0 to 90 • , where β is the angle between the LOS and the wind direction.The blue lines with squares and the dark lines are the mean error values from the entire range of z 0 and β for the averaged and point radial velocity variance, respectively.The red shaded areas denote the range of errors of the averaged radial variance.The equation used to estimate the integral length scale can be found inWang et al. (2016) with a mean wind speed of 7 m s −1 and Coriolis parameter of 10 −4 s −1 .The elevation angle φ = 45 • .

Figure 15 .
Figure15.Box plots of the relative standard deviation of the streamwise velocity variance u 2 and vertical momentum flux u w , assuming that they are derived from the six-beam method fromSathe et al. (2015) (denoted by the subscript lidar) and the sonic data (denoted by the subscript sonic) based on the observed hourly time series in Sect.4.1.Only those hourly data with u w < −0.1 m 2 s −2 are included.Relative standard deviations for the sixbeam method are calculated using the method in Appendix A and the observed attenuation and sampling errors in Sect.4, and those for sonic data are based on the derived turbulence statistics from the sonic data described in Sect.4.1 and error quantification methods given inLenschow et al. (1994).The sampling duration is denoted by the colors of box plots (see the legend).
μ2,R (T ) Estimated variance of radial velocity from lidar measurements over a sampling duration T (m 2 s −2 ) μ2,RVector of estimated radial velocity variance from lidar measurements (m 2 s −2 ) µ 2,r Ensemble variance of point radial velocity (m 2 s −2 ) ρ R Autocorrelation of averaged radial velocity ρ r Autocorrelation of point radial velocity Random sampling error covariance matrix of μ2,R (m 4 s −4 ) m Error covariance matrix of m (m 4 s −4 ) σ Q Standard deviation of radial velocity weight function (m) σ R Standard deviation of random measurement error of v R (m s −1 ) τ Autocorrelation or autocovariance time lag (s) τ 0 Integral timescale (s) φ Lidar elevation angle ( • ) A Diagonal matrix formed by a vector of a a Attenuation factor of radial velocity variance due to averaging D r Structure function (m 2 s −2 ) E 2 r,R Variance of random sampling error in variance estimates using averaged radial velocities (m 4 s −4 ) E s,R Systematic sampling error in variance estimates using averaged radial velocities (m 2 s −2 ) random sampling error in variance estimates using averaged radial velocities e s Relative systematic sampling error e s,ρ Estimated e s using autocorrelation function e s,b Estimated e s using the bootstrap method e s,R Relative systematic sampling error in variance estimates using averaged radial velocities e s,r Relative systematic sampling error in variance estimates using point radial velocities g Vector of the coefficients that relates m to µ 2,r G Matrix of the coefficients that relates m to µ 2,r K Matrix defined as (G T G) −1 G T I Identify matrix L 1 Integral length scale of the streamwise velocity (m) L R Integral length scale of the averaged radial velocity (m) L r Integral length scale of the point radial velocity (m) m Vector of the components of the velocity covariance tensor (m 2 s −2 ) m Estimate of m (m 2 s −2 ) www.atmos-meas-tech.net/9/4123/2016/Atmos.Meas.Tech., 9, 4123-4139, 2016