Towards improved turbulence estimation with Doppler wind lidar VAD scans

The retrieval of turbulence parameters with profiling Doppler wind lidars (DWL) is of high interest for boundarylayer meteorology and its applications. The DWL measurements extend beyond the observations with meteorological masts and are comparably flexible in their installation. Velocity-azimuth display (VAD) type scans can be used to retrieve turbulence kinetic energy (TKE) dissipation rate through a fit of measured azimuth structure functions to a theoretical model. At the elevation angle of 35.3◦ it is also possible to derive TKE. We show in this study how modifications to existing methods allow 5 to retrieve TKE and its dissipation rate even with a small number of scans, how a simple correction for advection improves the results at low altitudes and that VAD scans at different elevation angles with the same instrument provide comparable results of TKE dissipation rate after all filters and corrections. For this purpose, data of two experiments are utilized: First, measurements at the Observatory Lindenberg – Richard-Aßmann Observatory (MOL-RAO) are used for validation of the DWL retrieval with sonic anemometers on a meteorological mast. Second, distributed measurements of three DWL during the CoMet campaign 10 are analyzed and compared to in-situ measurements of the DLR Cessna Grand Caravan 208B. The comparison to in-situ instruments shows that the methods to improve turbulence retrievals from VAD scans introduced in this study are effective, especially at low altitudes and for narrow cone angles, but it also shows the limits of turbulence measurement with state-ofthe-art DWL in low turbulence regimes.

1. Abstract: The first 3 or 4 sentences could be probably be reduced to a single sentence in favor of allowing for a more quantitative summary later in the abstract. As it stands, the abstract lacks sufficient substance. The author should incorporate more hard results from the comparisons with sonic anemometers.
We believe that we need to explain the basic idea of the measurements with the introductory sentences even in the abstract but can still add more substance to it in the end. A modified version is given in the revised manuscript. 5 2. page 4 line 3: Not everyone will know where Upper Silesia is (including myself until I looked it up), I suggest ". . . were installed in Upper Silesia, in southern Poland (or where ever), . . . " We will add the country Poland in parantheses, the exact location is indicated on a map of Europe in Figure 2. 3. page 4 line 15: change ". . . and were finally fixed. . . " to ". . . and were finally choosen. . . " Ok. From the lidar manufacturer Halo Photonics, we were only provided with the official value of 1.5 µm. If the reviewer can give us a reference for the value he suggests, we will very much like to correct it in the table. On the other hand, this value is not very important for the study and if serious doubts about this value persist, we would rather not mention it at all. 15 5. Section 3.2: The author should include the equation for the measured azimuth structure function -since this is key for the dissipation rate retrieval methods.
We will include the equation just before Eq. 8 in the revised manuscript: (2) 20 6. Equation 3: The condition that ϕ = 35.3 • should be made more explicit to prevent possible misused.
The sentence introducing Equation 3 explicitly states that it is for the special case of ϕ = 35.3 • .
We add the suggested addition to the equation. 7. Page 6, lines 4-6: The author states that the "TKE dissipation rate is estimated through a fit of the measured secondorder structure function of horizontal velocity to the theoretical . . . " This statement implies that the observations are 25 adjusted to fit the model, when in fact it's the other way around, i.e. the model parameters are adjusted in order to fit the observations. Please rephrase.
We apologize for the mistake in language and correct it in the revised manuscript.
8. Page 6, lines 23-24: Similar to last comment. The author states that "A fit of the azimuth structure to the equation . . . " again implies that the observations are being adjusted to fit the model, when in fact it's the other way around. Please rephrase.
We apologize for the mistake in language and correct it in the revised manuscript. 9. Page7, line 1: The author states that "Scanning with Doppler lidar in a VAD implies a volume averaging of radial velocities in longitudinal and transverse directions." Aside from the grammatical errors, this statement is not generally true because transverse averaging is not an issue for step-stair scans, only for continuous motion scans. The author 5 should briefly mention the two different types of scans in their introduction. Also, the author should define what they mean by longitudinal and transverse (i.e. along the beam, and orthogonal to the beam).
Velocity azimuth display (VAD) is a term from radar technology where continuous motions of the azimuth motor are the standard.
Step-and-stare scans like Doppler-Beam-Swinging techniques are not considered in this manuscript. In any case, even step-and-stare scans with pulsed DWL have a transverse averaging effect due to the pulse-averaging over a 10 certain accumulation time. We add to the manuscript that transverse averaging is regarded for scans with a continuous motion. We also define longitudinal and transverse in the revised manuscript.
10. Page 7 line 5: change ". . . radial wind speed. . . " to " . . . radial velocity. . . ". Wind speed is a (positive) scalar, velocity is a vector. In this sentence your talking about the radial component of the velocity vector. "Radial wind speed" makes no sense. 15 We change the wording in the revised manuscript.
11. Page 7 starting at line 7: The discussion here is a bit disjointed and difficult to follow. Equations 5-7 should be listed after the sentence on line 5 (starting with "It is based on the decomposition. . . "). As it is, these equations are introduced without any corresponding text. One suggestion might be: "In Smalikho and Banakh (2017), this method has been combined with the E89-method to yield TKE, and the momentum 20 fluxes. It is based on the decomposition of radial velocity variance into its subcomponents, i.e.
σ 2 a = σ 2 r − σ 2 t (4) where σ 2 L is the lidar measured variance, σ 2 a is the lidar measured variance without instrumental error, σ 2 e is the instru- 25 mental error variance, and σ 2 t is turbulent broadening of the lidar measurement. In Smalikho and Banakh (2017), all of these variances and corresponding structure functions are calculated for single azimuth angles and then averaged." We agree that the modifications make the text easier to follow and incorporate the changes in the revised manuscript.
12. Page 7, line 17: Recommend changing "Substituting σ 2 e in Eq. 7 with Eq. 8 yields:" to "Combining Eq. 7 with Eq. 8 yields:" 30 We change the sentence in the revised manuscript according to the suggestion of the referee.
13. Page 7, lines 18-23, including equations 9, 10 and 11: There is a dependence on the separation distance on the right side of equation 9 that presumably cancels such that the right side is effectively constant, i.e. independent of separation distance. This is a subtle point that is not made by the author. Also, in equation 10, the author has substituted Ψ l with Ψ 1 without any explanation or justification. Please explain.
We agree that an explanation is lacking here. Since the instrumental error σ 2 e is assumed to be a constant offset of azimuth 5 structure function D a (ψ l ) and the lidar measurement D L (ψ l ), l can actually be chosen arbitrarily in the TKE equation.
It is set to l = 1 because potential random errors like unstationary flow will be least effective for small separation angles.
14. Page 8, line 10-11: The author states ". . . from VAD scans with other elevation angles as well." You should be a bit more specific here, since readers may not know what you mean by "other elevation angles." I assume you're referring to elevation angles different from 35.3 • . 10 We change the sentence to explicitly say "elevation angles different from 35.3 • .
15. Page 8, line 13-15: The author states "The value of l = 9 is chosen following the example of Smalikho and Banakh (2017) and corresponds to l∆θ = 9 • as it was found to be suitable in all conditions in that study." The discussion up to this point had been fairly general. Now, suddenly the author is referring to a very specific VAD scan. The author should be a bit more specific as to which scan (and which experiment) they are referring to.

15
Here we define the upper separation angle that will be used in the further manuscript and in all experiments. The separation angle is not referring to a specific VAD scan and in our opinion can be introduced here. We rephrase in the revised manuscript to state that this separation angle was found by Smalikho and Banakh (2017) to be a reasonable value in the ABL. It is illustrated by an example structure function in Figure 3.

21.
Page 10 line 5: The author states "(here: g=360 for all azimuth angles)". The reference to a specific value of g here is a bit perplexing. Please explain.
In this study we work with VAD scans with 1 • azimuthal resolution, which yields 360 values per scan. Since we introduce a general method here, we will remove this information at this point of the text.

Equation 20:
The summation is over j, but there is no dependence on j in the quantity being summed. Please explain.
This is a mistake. The summation is over the variable m (index of the azimuth angle within one VAD scan).
24. Page 10, line 16-17: The author states "Measured PDFs of the variables . . . are fit to the model PDFs to obtain an estimation of the corresponding standard deviations σ 1 , σ 2 and σ 3 and probability of bad estimates P 1 , P 2 and P 3 ." This statement implies that the observations are fit to the model. In other words, the observations are tweaked to get agreement with the model. That's certainly not what is happening. Please rephrase. 10 We apologize for the confusion in language and rephrase in the revised manuscript. A first guess of the standard deviations is needed to find the ±3.5σ region for the integral over the PDF. All the details of this method are given in . Since this method is not essential for this study, we do not want to expand too much on it in this manuscript. It is mostly relevant if better measurements in low-signal conditions are targeted. 20 26. Section 3.2.3: In this section the author throws down a series of equations without adequate discussion. The authors need to do a better job explaining their line of reasoning.
We thank the referee for their recommendation on improving the section and incorporate it in the revised manuscript.
27. Appendix C, lines 17: The author introduces the quantity Xj (i.e. a 1D vector), and then in equation C1 it is indicated to be a 2D matrix, i.e. Xij. Please explain. 25 j is the subsample index and i is the index for each element of the subsample. We clarify this in the revised manuscript.
28. Appendix D: I find no mention of the "FSWF-retrieval" in the paper (did I miss it?). Please elaborate and provide relevant citations.
Filtered sine-wave fitting is introduced with the corresponding reference in Section 3.2.1, but without giving the abbreviation FSWF. We add this in the revised manuscript. 1. The method to retrieve ε from sonic anemometer (and research aircraft) measurements is based on averaging 2 min estimates of ε over 30 min sampling time. Taking the average means that the distribution of ε 2min is assumed to be Gaussian, which is not necessarily true as the magnitude of ε 2min may vary over several scales of magnitude. When the 5 distribution is not Gaussian, the average will introduce a bias to the ε 30min estimate. Therefore, authors should check the shape of the distribution of ε 2min values during each 30 min period and choose an estimate for ε 30min that is more representative for the distribution of ε 2min values, such as the median of these values.
We thank the referee for this comment which is very valid and should be considered in the evaluation of sonic anemometer data. An arithmetic mean of dissipation rates ε is not the best solution given the exponential character of the variable.

10
Instead of using the median as suggested, we however believe now that for this study it is more reasonable to calculate the structure function over the full half-hour period and estimate ε from it.  calculated 2-minute periods because they wanted to see bursts of turbulence on shorter timescales. In our case we are however comparing to lidar measurements that are averaged over half-hour periods, so that a comparison can be best made with the same period for the calculation of the structure function of sonic anemometer data. 15 We did investigate the difference between median of 2-minute estimates, mean of 2-minute estimates and 30-minute estimates and can confirm that the referee is right that the 2-minute mean is skewed towards larger values compared to the the median approach. A systematic error can however also occur when the median underestimates the dissipation rate within the half-hour (see Fig. 1). We will not present these results in the revised manuscript because we believe that 30-minute stucture function estimates are the right choice for this analysis.

20
In any case, for all possible estimates of ε from sonic anemometers, the differences are small enough to not change the conclusions that are made for the comparison to the lidar retrievals. In Fig. 2 below, we show the differences. The best correlations can be found with the 30-minute structure function estimate. The median estimate shows the highest bias for the S17 method, but the general conclusions remain the same. be considered as robust as sonic anemometer data. In fact, the explanation for using aircraft measurements is provided only in Section 5 on lines 8-9 of page 21. This explanation should be given already in the abstract but also in the Introduction (Section 1) and maybe also in Section 3.3.
We add the term "research aircraft" in the abstract explicitly and also explain that the research aircraft is a unique 2. Section 2, page 2, lines 33-34: "data from two different sites and sets of instruments": This is not clear: you use data from four DWLs (three in Upper Silesia, one at Falkenberg), from two sonic anemometers and from one research aircraft.
Although you introduce the measurements in two Figures (Figures 1 and 2), it does not change the fact that you have 5 several type of instruments and sites, not just two of each.
We change the text to only state that measurements from two different experiments are analyzed. The details of the sites and the instrumentation is given in detail throughout the section: "In this study, data from two different experiments are analyzed. Both of the experiments and the instrumentation is introduced in this section." 3. Section 2.1., page 3, line 19: remove "and infrared gas analysers LI7500 (LiCor Inc.)" because you do not use this data.
We remove the information about the gas analyzer in the text: "Continuous turbulence measurements (20 Hz sampling frequency) using sonic anemometer of type USA-1 (METEK GmbH) are performed at the 50m and 90m levels of the tower and have been used for validation purposes. The instruments are mounted at the tip of the booms pointing towards South." 4. Section 2.2, the first two paragraphs on page 4: the description of the CoMet mission is too detailed and not relevant 15 to the topic of this paper, as here the aim is not to investigate CO2 or methane but to develop DWL data processing methods. Please, provide here only such information that is relevant for the present study.
We believe that the information about the CoMet campaign is relevant, because it puts the DWL measurement into context for the special issue to which this manuscript has been submitted. This manuscript is important for the CoMet research community, so that we want to make a connection to the overall project. Ok.
(b) The acronyms of the Doppler lidars in Upper Silesia region are misleading: for the research aircraft you use the 5 acronym "DLR" (e.g. in the abstract but also on line 8 of page 4) and for Doppler lidars you have introduced the acronym "DWL". Why do you introduce here another acronym for DWL, i.e., "DLR"? Please, use only one acronym for Doppler wind lidars throughout the paper.
We change the naming of the lidars in the CoMet campaign to DWL#1-DWL#3 in the revised manuscript. The tense of the verb are all changed to past tense. .
See above answer to major comment. It is actually Eq. 3, but the way this is written is confusing and we rephrase in the revised manuscript.
10. Page 7, line 20: Why does Ψ l change to Ψ 1 here? Please, explain what it means that l=1?
We agree that an explanation is lacking here. Since the instrumental error σ 2 e is assumed to be a constant offset of azimuth structure function D a (ψ l ) and the lidar measurement D L (ψ l ), l can actually be chosen arbitrarily in the TKE equation. 25 It is set to l = 1 because potential random errors like unstationary flow will be least effective for small separation angles. 11. Page 7, line 25: Typo: "Kolmogorov-Obhukov spectrum" ok.
14. Page 10, Equation 20: There is no index j in the equation (which is included in the summation). Moreover, are there some 5 parenthesis missing?
The index j is a mistake, it should be m. There are no paranthesis missing.
15. Page 12, Section 3.3: Again, it is not clear why research aircraft data is used: is it used a) because it gives more reliable results than DWLs and therefore can be used for validation of DWL data or b) is it used because it would be interesting to know how good the research aircraft data is compared to DWLs (and sonics)? 10 We consider turbulence measurements with flow probes on research aircraft a well-established method which provides more reliable measurements than a DWL since it does not depend on many assumptions except the Taylor's hypothesis.  17. Page 14, line 5: "modified version W19 introduced in this study", maybe you should use acronym W20 for the method introduced in this paper? 20 We changed the acronym throughout the manuscript.
18. Page 15, Figure 6: Does the biases in (b) and (e) include all points or only those after the advection filter?
All points except the grey points below the threshold are used. We will clarify in the caption. We think it is important to use all points to evaluate the effect of the advection correction.
19. Page 17, Figure 8: Why there is an oscillatory pattern in TKE bias as a function of horizontal wind speed? The oscillatory 25 pattern is more significant than the differences between the methods, and therefore it should be discussed. Could you also provide here the amount of data in each bin, maybe by adding another y-axis for that?
There is no physical reason for an oscillatory pattern in TKE vs. horizontal wind speed. The shape of the curve is merely a specialty of the dataset where more values of certain bins fall in nighttime hours with lower TKE and thus also lower absolute error in TKE. We add the number of data points in the bin in Figure 8 of the manuscript as shown here in Fig. 4 30 20. Page 18, line 9: "Here, it shows that the difference between S17 and W19 only occurs at the very lowest level" -this cannot really be seen from Figure 12.  It shows that for DLR1, there is a small difference between the two methods, but it is not very well visible. We change the discussion of this figure in the revised manuscript. 21. Page 19, Figure 11: Another horizontal axis with a km scale would be nice, because in the text you give the length of the flight path in kilometers.
We change the plot to give the x-axis in kilometers instead. We rescaled the figure and changed the color of the greyscale for the S17-method slightly for better readability of the plot.
23. Page 20, lines 16-18: "The advection effects are most relevant at the lowest measurement heights where the spatial 10 separation of lidar beams along the VAD cone ∆y are small." This is true, but what could perhaps be also mentioned is that the advection speed increases with height because of the logarithmic wind profile.
We add this thought to the discussion. In our observation, the effect of increasing wind speed with height is however not as strong as the effect of a counteracting larger ∆y.
24. Page 21, lines 1-3: "dissipation rates of values smaller than 10 −3 m 2 s −3 are underestimated by the lidars, likely because 15 the small scale fluctuations that are carrying much of the energy in these cases, cannot be resolved any more." This can be true, but you should still check the method to retrieve dissipation rates from sonic anemometers as mentioned in the previous comments.
As shown above we have evaluated the sonic anemometer retrieval and agree that there can be differences between different methods to obtain a representative value of ε within the half-hour period. However, the strong underestimation 20 of lidar measurements below 10 −3 m 2 s −3 is much more significant and found in any case. 25. Page 21, lines 8-10: this information should be provided much earlier in the manuscript. This is not just a result but also the motivation to use research aircraft data in the first place.
More information on the motivation to use aircraft data is added to Section 2.2. in the revised manuscript.

Relevant changes to the manuscript
We list here the relevant changes to the manuscript:

Abstract:
-Text modifications in response to referee comments. More focus on results of the study.

5
-Text modifications in response to referee comments.

Section 2:
-Text modifications in response to referee comments.

Section 3:
-Text and equation modifications in response to referee comments.

Section 4:
-Text and equation modifications in response to referee comments.
-Added subplot to Figure 5. - Figure 8&10: Added number of data points per bin in the plot.
- Figure 12: Changed line colors and scaling.

300
-Text modifications in response to referee comments.

Appendix:
-Text modifications in response to referee comments.

Introduction
The observation of turbulence in the atmosphere and in particular the atmospheric boundary layer (ABL) is of great importance for basic research in boundary-layer meteorology as well as in applied fields such as aviation, wind energy (van Kuik et al., 2016;Veers et al., 2019) or pollution dispersion (Holtslag et al., 1986).
A wide range of instruments are used to measure turbulence: sonic anemometers are nowadays the most popular in-situ in-5 strument which can be installed on meteorological masts and provide continuous data of the three-dimensional flow and its turbulent fluctuations (Liu et al., 2001;Beyrich et al., 2006). For in-situ measurements above the height of towers, airborne systems are applied such as manned aircraft (Bange et al., 2002;Mallaun et al., 2015), remotely-piloted aircraft systems (RPAS, van den Kroonenberg et al., 2011;Wildmann et al., 2015) or tethered lifting systems (TLS, Frehlich et al., 2003) which can be equipped with turbulence probes such as multi-hole probes or hot wire anemometers. A different category of instru-10 ments are remote-sensing instruments such as radar, sodar and lidar which can measure wind speeds and allow the retrieval of turbulence based on assumptions of the state of the atmosphere and the structure of turbulence. In this study, we focus on ground-based Doppler wind-lidars (DWL), which have become increasingly popular in boundary-layer research because of their light weight ::: ease :: of :::::::::: installation, invisible and eye-safe lasers, their reliability and high availability which is only restricted by clouds/fog and rain or very low aerosol content in the atmosphere.

15
A variety of methods already exist to retrieve turbulence from DWL measurements. They can be categorized according to the respective scanning strategy applied: the simplest scanning pattern is a constant vertical stare to zenith, which allows to obtain variances of vertical velocity and estimates of turbulence kinetic energy (TKE) dissipation rate (O'Connor et al., 2010;Bodini et al., 2018). More complex are conical scans (velocity azimuth display, VAD) with continuous measurements along the cone (Banakh et al., 1999;Smalikho, 2003;Krishnamurthy et al., 2011;Smalikho and Banakh, 2017). These scans include informa-20 tion on the horizontal wind component as well. A simplification of VAD scans are Doppler-beam-swinging (DBS) methods, that reduce the number of measurements taken along the cone to a minimum of 4-5 beams and thus increase the update rate for single wind profile estimations (Kumer et al., 2016). Both, VAD and DBS are popular scanning strategies that are applied in commercial instruments. Kelberlau and Mann (2019a, b) introduced new methods to obtain better turbulence spectra from conically scanning lidars by corrections for the scanner movement. Significantly different scanning strategies are vertical (or 25 horizontal) scans which can also provide vertical profiles of turbulence (Smalikho et al., 2005), but even allow deriving twodimensional fields of TKE dissipation rate . Multi-Doppler measurements require more than one lidar with intersecting beams, but do not need assumptions on homogeneity to measure turbulence at the points of the intersection directly (Fuertes et al., 2014;Pauscher et al., 2016;. For operational or continuous monitoring of vertical profiles of turbulence in the ABL, VAD or DBS scans are most suitable. At an elevation angle of 35.3 • , a VAD scan 30 allows to retrieve TKE, its dissipation rate, integral length scale and momentum fluxes according to a method that was first developed for radar by Kropfli (1986) and adapted for lidar later by Eberhard et al. (1989) using the variance of radial velocities along the scanning cone. Further improvements of this method have been implemented by Smalikho and Banakh (2017) and , which also account for lidar volume averaging effects. ::: We :::::::: introduce ::::::::::: modifications ::: on ::: the ::::::::: estimation :: of ::::::: structure ::::::: function ::: and :::::::: variances :: in ::::: order :: to ::: be ::: able :: to ::::::: retrieve ::::::::: turbulence ::::::::: parameters :::: from :: a :::::: smaller ::::::: number :: of :::: VAD :::::: scans. For conditions with significant advection, the method is not applicable and causes large ::: can ::::: cause : errors, especially at low altitudes where the cone diameter of the VAD scan is small. With this study, we propose a method to significantly reduce this error. We also apply the turbulence retrieval to VAD scans with 75 • elevation angle, which still allows to retrieve TKE dissipation rate. In this case the advection correction is particularly important. The experiments that were carried out are explained in Sect. 2. The and an outlook are given in Sect. 5.

Experiment description
In this study, data from two different sites and sets of instruments ::::::::: experiments : are analyzed. Both of the sites ::::::::: experiments : and 10 the instrumentation is introduced in this Section :::::: section.

The MOL-RAO Falkenberg field site
The Continuous turbulence measurements (20 Hz sampling frequency) using the eddy-covariance method with sonic anemometer-thermomete :::: sonic ::::::::::: anemometer :: of :::: type USA-1 (METEK GmbH) and infrared gas analysers LI7500 (LiCor Inc.) are performed at the 50m 30 and 90m levels of the tower and have been used for validation purposes. The instruments are mounted at the tip of the booms pointing towards Southwith the LI7500 behind the USA-1. The DWL measurements are particularly helpful to support the CoMet measurements by providing wind information which is essential to derive emission flux estimates from passive remote sensing (Luther et al., 2019) or in-situ measurements of mass concentrations (Fiehn et al., 2020). The DWL wind information can also be used to validate modeled wind of the transport models for greenhouse gases. The lidars were remotely operated during the whole CoMet campaign period from 16 May 2017 to 17 June 2017 and were continuously measuring. The locations of the three lidars were planned to cover the whole region of 15 interest and were finally fixed :::::: chosen based on logistical constraints.
The lidars were operating in VAD modes with two different elevation angles. Since the focus for the CoMet campaign was on continuous wind profiling and a good height coverage was desired, the lidars were programmed to perform VADs with an elevation angle of 75 • (see Tab. 1, VAD75) for a longer period, i.e. 24 scans (≈29 minutes), followed by only six scans (≈7 minutes) at 35 :::: 35.3 • elevation (VAD35) for turbulence retrievals.

VAD turbulence measurements
Methods to retrieve turbulence parameters from VAD scans are well-known and a variety of different methods exist. The method we refine in this study is based on the theory that was originally described by Eberhard et al. (1989) for lidar measurements. The variance of radial velocities σ 2 r depends on the range gate distance R, the azimuth angle θ and the elevation angle 15 ϕ. It is calculated from the measured radial wind speeds :::::::: velocities V r : v r (R, θ, ϕ, t) = V r (R, θ, ϕ, t) − V r (R, θ, ϕ) From a partial Fourier decomposition (see App. A) and for the special case of ϕ = 35.3 • a simple equation for E TKE is derived: :::::: . (3) In this equation, σ 2 r is the mean of the variance of radial velocites over all azimuth angles. In the following, we will refer to this method as E89-retrieval.
The retrieval method for ε using Eq. 16 and TKE using Eq. 11 will be referred to as S17 in the following.

Modifications for small number of scans
The VAD at ϕ = 35.3 • during the CoMet-campaign was not run continuously, but only six individual scans are performed successively before switching back to the VAD at ϕ = 75 • as described in Sect. 2.2. This means that only six data points 20 are available to calculate variance and mean of the radial wind speeds :::::::: velocities : at each azimuth angle, which cannot be considered a solid statistic. We introduce two modifications of data processing to overcome this problem which are based on the assumptions of stationary and homogeneous turbulence.

Practical implementation of the ensemble average
In Eq. 2, V r (R, θ, ϕ) can be calculated as the arithmetic mean of radial wind speeds ::::::: velocities : at specific azimuth angles: where N is the number of scans. Instead of this approach, we suggest to use the reconstructed radial wind speed ::::::: velocity from the retrieved wind field over all individual scans as the expected value in the variance calculation. For the retrieval of the three wind components (û,v,ŵ), filtered sine-wave fitting :::::: (FSWF) : is applied (Smalikho, 2003). The reconstructed radial wind speeds ::::::: velocities :V are then used as the expected value in the variance calculation: V =ŵ(R) sin ϕ +v(R) cos ϕ cos θ +û(R) cos ϕ sin θ (18) With this approach, all measurement points in the VAD with the same elevation angle are used to obtain the expected value 5 V r and thus, a better statistical significance is achieved. This method has also been proposed in Smalikho and Banakh (2017) as a practical implementation of Eq. 2.

Averaging of variances
In (Smalikho and Banakh, 2017) :::::::::::::::::::::::: Smalikho and Banakh (2017), the variances of the lidar measurements are defined as the average of variances at individual azimuth angles: The variances σ 2 r (θ m ) are variances of a subsample of radial wind speeds ::::::: velocities : of the VAD (i.e. those at a specific azimuth angle θ m ). We use a simple relation between the variances of subsamples and the total variance of a dataset (see Appendix C). Applying this to the radial wind speed ::::::: velocity variances yields: where k is the number of samples at each azimuthal angle, g is the number of subsamples (here: g = 360 for all azimuth 5 angles) and n is the number of total samples in the dataset (n = gk). Since the mean of the radial wind speed fluctuations v r = 0 :::::: velocity :::::::::: fluctuations :: is :::: zero by definition, it is :: Eq. ::: 21 ::::::: becomes:

Filtering of bad estimates
Improvements of turbulence estimates in low signal conditions can be achieved with filtering of bad estimates as described in 10 Stephan et al. (2018). This approach is not based on the calculation of the azimuth structure function from measured radial wind speeds :::::::: velocities, but uses probability density functions (PDFs) and their corresponding standard deviations. The model PDF is defined as a Gaussian function with a filter term P : where P is the probability of bad estimates of x, σ is the standard deviation of the PDF and B v is the velocity bandwidth of 15 the lidar. Measured PDFs of the :: For ::: the :::::::: measured : variables v r (R, θ), ∆v r (R, θ + ∆θ) and ∆v r (R, θ + l∆θ)are fit to the model PDFs to obtain an estimation of : , ::: the :::::: best-fit ::::: model ::::: PDFs ::: are ::::: found :: to ::: to ::::: obtain : the corresponding standard deviations σ 1 , σ 2 and σ 3 and probability of bad estimates P 1 , P 2 and P 3 . However, since the PDFs cannot be assumed Gaussian in atmospheric turbulence, the standard deviations are finally calculated as the integral over the measured PDFs in the range ±3.5σ according to .

20
Replacing σ 2 L with σ 2 1 , D L (ψ 1 ) with σ 2 2 and D L (ψ l ) with σ 2 3 in Eqs. 11 and 16 yields: As suggested in , Eqs. 24 and 25 are only used if P > 0. We also introduced a quality control to discard 25 any measurements with P > 0.5 for best results. In practice, this method will thus only be applied in some conditions when the signal is weak and can extend the range of vertical profiles to some degree.
The effects of advection on the turbulence estimation is largest in the lowest levels of the VAD-scans, because ∆y is small 15 compared to U ∆t in this case. The retrieval method including the filtering for bad estimates, and the advection correction is referred to as W19 :::: W20 : in the following.

Quality control filters
In order to fulfill the assumptions that are made with regards to the turbulence model and the turbulence retrieval method, the data is filtered according to the criteria given in Smalikho and Banakh (2017): R ω s | V | For the purpose of evaluating the methods in a broad range, we set mild criteria for Eq. 32 and 34 using Equations 32 and 33 are criteria that require the integral length scale L v to be larger than the sensing volume of the lidar in transversal (∆y) and longitudinal (∆z) direction. Unfortunately, there is no independent measurement of L v at all heights of the VAD scan, so that it is derived from the lidar measurement itself as L v = 0.3796 E 3/2 ε (Smalikho and Banakh, 2017). The filter criteria in Eq. 34 is a filter for conditions with significant advection which distorts the measured structure functions and is only applied if the method described in Sect.3.2.3 is not used.

5
Except for the retrieval method WS19 :::: W20, which uses the filtering of bad estimates, we set fixed CNR filter thresholds adapted to the lidar type. Since the turbulence retrievals are very sensitive to bad estimates, we set the CNR thresholds to conservative values that are given in Tab. 1.
An overview of all retrieval methods and their characteristics and filters that are applied is given in Tab. 2.

Turbulence estimation from airborne data 10
The estimation of turbulence parameters from the wind measurement system on the DLR Cessna Grand Caravan 208B (Mallaun et al., 2015) is done very similarly to the in-situ estimations from the sonic anemometer. TKE is calculated from the sum of variances as described in Sect. 3.3. Dissipation rate is also calculated from the second order structure function, but with different bounds for the time lag. For the flight data we use τ 1 = 0.2 s and τ 2 = 2 s, corresponding to approximately 13-130 m lag at 65 m s −1 mean airspeed.

15
To evaluate the heterogeneity of turbulence due to changing land use along the flight legs of more than 50 km length, we divided the legs into sub-legs of 6.5 km (i.e. 100 s averaging time) and calculate turbulence for each leg individually. The location of the legs and the sublegs is shown in Fig. 2. 14 4 Validation

Comparison to sonic anemometer
Best possible validation of the methods introduced in Sect. 3 can be performed with the lidar in close proximity to the meteorological mast such as at the measurement site at MOL-RAO. The sonic anemometers at 50 m and 90 m on the mast almost coincide with measurement levels of the lidar at 52 m and 93.6 m respectively. Since the VAD-retrieval with elevation angle 5 of 35.3 • yields TKE as well as its dissipation rate, both turbulence parameters can be compared to values obtained from the sonic anemometers. In this section we will evaluate the methods described in Sect. 3.2 and in particular the validity of the assumptions made in Sect. 3.2.1 and the efficiency of the advection correction described in Sect. 3.2.3.

Validation of modified variance
In Sect. 3.2.1 we introduced two modifications on the calculation of the averaged lidar radial wind speed :::::: velocity : variances.

Comparison of lidar retrievals
The MOL-RAO dataset allows us to compare the retrievals without consideration of lidar volume averaging (E89 and S17A) to the S17 retrieval and its modified version W19 :::: W20 : introduced in this study. For this purpose, the individual retrieval results 25 are compared to the sonic anemometer estimates of TKE and its dissipation rate ε. Figure 6 shows the scatter plots for TKE at the two measurement levels. For each method, the coefficient of determination R 2 c of the linear regression between sonic measurement and lidar retrieval is given, as well as a bias which is calculated as b = (y − x) for TKE and b log = (log 10 y − log 10 x) for ε. We find that with the E89-method, TKE is systematically underestimated, as expected. In contrast to that, the S17-method yields slightly overestimated TKE-values, if no advection filter is applied (light red dots), but a good agreement with the sonic 30 anemometer in the absence of advection (red dots). The overestimation of TKE is larger for the lower level at 50 m compared to the 90-m level which we attribute to the smaller averaging volume ∆y. Our refined version of the retrieval, including the advection correction improves the results slightly and especially for the high turbulence cases (corresponding to high wind speeds) at the 50-m level, as the scatterplots show. Figure 7 gives the scatterplot comparison of ε-retrieval from S17A, S17 and W19 :::: W20 : with the sonic anemometer respectively. E89 does not provide an estimate for ε. Even more clearly than for TKE, the underestimation of the method without consideration of lidar volume averaging (S17A) is found. Also, a now positive 5 bias of b = 0.09 m 2 s −2 of lidar estimates with the S17-method at the 50-m level is reduced with the advection correction to b = 0.06 m 2 s −2 . It is evident from the ε-estimates that all lidar retrievals underestimate turbulence significantly compared to the sonic anemometer in the low turbulence regimes, and in particular for values smaller than 10 −3 m 2 s −3 ::::::::::::: 2 · 10 −3 m 2 s −3 , which is why these values have been excluded for the estimation of biases (grey dots).

Evaluation of advection error 10
To evaluate the error that is caused by advection in the S17-retrieval, all data that was collected at MOL-RAO was binned into wind speeds with a bin-width of 1 m s −1 . The mean absolute error between the lidar retrievals and the sonic anemometers at the respective level is calculated and shown in Figure 8 for TKE and ε. Although the averaged errors of TKE are small in general, it shows that the W19-method ::::::::::: W20-method does reduce the error in comparison to the S17 method at the 50-m level. The error for ε at the 50-m level increases with wind speed, but less for the W19-method ::::::::::: W20-method. Hardly any improvement is found 15 for the already small errors of TKE and ε at 90 m.

Comparison of elevation angles
VAD-scans with 35.3 • allow to retrieve TKE as well as momentum fluxes using the methods described in Sect. 3. The disadvantage compared to VAD scans with larger elevation angles is that at the same range of the lidar line-of-sight measurement, lower altitudes are reached. If the limit of range is not given by the ABL height in any case, this can lead to significantly lower data availability at the ABL top. Another advantage of greater elevation angles is that the horizontal area that is covered with 5 the VAD and thus the footprint of the measurement is much smaller than with low elevation angles. From the theory derived in Sect. 3.3, we see that the dissipation rate retrieval does not depend on the elevation angle and can thus be obtained from VAD scans with 75 • elevation angle as well if the assumptions of isotropy and homogeneity hold. However, since the VAD at 75 • has the more narrow cone, the separation distances ∆y at respective measurement heights are smaller and thus the sensitivity to advection errors is expected to be larger. The measurements with two different elevation angles at the CoMet-campaign 10 allow to compare dissipation rate retrievals for both kinds of VAD scans with the restriction that they are not simultaneous, but sequential. Figure 9 shows the comparison of both types of VADs in scatter plots of three measurement heights for the whole campaign period and DLR ::::: DWL#1. In general, a large scatter is found between the two types of VAD which can be attributed to the different measurement times, different footprint and heterogeneous terrain. Applying a filter for significant advection as described in Sect. 3.2.4 removes most of the measurement points at the 100 m level in the 75 • VAD. Without the filter (grey and red 5 points), large, systematic overestimation against the VAD at 35.3 • is found, which can still be seen at 500 m, but is not found at 1000 m any more. With the advection correction of W19 :::: W20, the systematic error is reduced, but the random errors remain.
As for the comparison with the sonic anemometer, we evaluate the error in dependency of wind speed by binning the data in wind speeds between 0 m s −1 and 10 m s −1 (see Fig. 10). A clear trend is found for the 100 m-level which can be significantly 10 reduced with the W19 :::: W20 : method. A very small difference between the two elevation angles at the 500-m level is only reduced very little in the W19 :::: W20 : method.

Comparison to D-FDLR
At CoMet, no meteorological tower with sonic anemometers on levels that could be compared to the Doppler lidars was available. Instead, the aircraft D-FDLR was operating with a turbulence probe and provided in-situ turbulence data. On 5 June 2018, the aircraft was flying a so-called 'wall'-pattern, with long, straight and level legs at five altitudes (800 m, 1000 m, 1100 m, 1300 m and 1600 m). At least the lowest three levels of this flight allow a comparison to measurement levels of the 5 top levels of lidar measurements on this day. Figure 11 shows the measured TKE of D-FDLR in the five flight levels. Only at the lowest two light levels (i.e. 800 m and 1000 m), significant turbulence is measured, with strong variations along the flight path. Figure 12 shows the measurements from D-FDLR and the lidars DLR ::::: DWL#1 and DLR :::: DWL#3 that were taken between 1400 and 1530 UTC as vertical profiles. The solid red line gives the average of all sub-legs along the 50 km flight of D-10 FDLR, and the shaded areas give the range between the minimum and the maximum at each height. It shows nicely how the measurements of turbulence at the DLR :::: DWL#3-site are significantly lower than at DLR :::: DWL#1, which we attribute to the lake fetch. The D-FDLR measurements of TKE almost match with the DLR ::::: DWL#1 site and are all higher than at the DLR :::: DWL#3-site, which fits to the environmental conditions of heterogeneous land-use. Figure 12a also gives the comparison between E89-, S17-and W19-estimates :::::::::::: W20-estimates : of the same dataset. Here, it shows that the difference between S17 15 and W19 :::: W20 : only occurs at the very lowest level, but the underestimation of the E89-method is found up to 750 m. In dissipation rate estimates, the DLR :::: DWL#1 measurements are at the low end of the range that was measured with D-FDLR and the lake-site measurements are even smaller. The estimates from 35.3 • -scans and 75 • -scans agree very well, especially at the higher levels, which shows that the assumption of isotropy and homogeneity seem to hold. The presumed underestimation of ε of lidar retrievals compared to in-situ measurements at absolute values of 10 −3 m 2 s −3 is consistent to what was found for the comparison to sonic anemometer measurements. This single case of airborne measurements compared to lidar retrievals at higher altitudes can however not provide any statistical validation. Figure 13 shows measurements of ε retrieved with the W19-method :::::::::: W20-method : for the VAD scans with 75 • elevation and all three lidars. It shows that the growth of the boundary-layer with its increased turbulence can be nicely captured by the lidars.

5
There are some differences between the three locations, especially lower turbulence close to the ground at the DLR :::: DWL#3location and a higher boundary layer at the DLR :::: DWL#2 location. More studies will be necessary in future, analyzing the data of the whole campaign to improve the understanding of land-atmosphere interaction in this case.

Conclusions
In this study, we used four different methods to retrieve turbulence from the same data obtained through lidar VAD scans.
The MOL-RAO experiment allowed us to show that methods which do not account for the lidar volume averaging effect underestimate turbulence compared to sonic anemometers at 50 m and 90 m systematically. :::: This ::: has :::: been :::::: shown ::: for ::: the :::: first :::: time :::: with :::: such :: a ::: big ::::::: dataset. : The S17-method tackles this problem, but introduces an overestimation in our dataset. Parts 5 of this overestimation can be attributed to advection, which distorts the retrieval of the azimuth structure function and the transverse filter function in the lidar model. The advection effects are most relevant at the lowest measurement heights where the spatial separation of lidar beams along the VAD cone ∆y are small. :: is ::::: small. :::: This ::::: effect :: is ::::::: stronger :::: than ::::::::: increasing ::::: wind ::::: speeds :: at :::::: higher ::::::: altitudes ::: in ::: our ::::::::::: observations. : We propose a correction for this issue and show here that our method reduces the systematic errors compared to the sonic anemometers at the 50-m level. It is also shown :: To ::::::: confirm ::: that ::::::::: advection :: is ::: the 10 ::::: reason ::: for :::: this ::::::::::: improvement ::: we ::::: show that the bias increases with wind speed. With all retrievals, dissipation rates of values smaller than 10 −3 m 2 s −3 are underestimated by the lidars, likely because the small scale fluctuations that are carrying much of the energy in these cases, cannot be resolved any more. A remaining piece of uncertainty are the lidar parameters ∆p and T w which are given by the manufacturer for the lidar type, but could potentially differ for individual lidars. Exact knowledge about these parameters could reduce the uncertainty of the model functions F and A (see App. B) and thus improve the corrections 15 of the volume averaging effects. It is conceivable that the observed overestimation of the S17 (and W19 :::: W20) based TKE can partly also be attributed to these uncertainties.
The aircraft measurements that were carried out during the CoMet campaign were used to show the agreement of the lidar retrievals with in-situ measurements at higher altitudes. : It :: is ::: the :::: first :::: time :::: that :::: these ::::: lidar :::::::::::: measurements :::: have :::: been ::::::::: compared Figure 13. Diurnal cycle of TKE dissipation rate on 5 June 2018 at the three lidar locations calculated with the W19-method ::::::::: W20-method. :: to ::::: in-situ ::::::: aircraft :::: data. : Unfortunately, only measurements of one day allowed a comparison and the spatial separation of the measurements introduces additional uncertainty. It was found that TKE estimates from lidar and aircraft compare rather well, but the small values of dissipation rates at these heights are underestimated by the lidar to a similar order of magnitude as for low turbulence conditions in the sonic anemometer comparison. Dedicated experiments will be necessary in future to provide more comprehensive validation datasets for turbulence retrievals with lidar VAD scans. Given the larger separation 5 distances ∆y of the lidar beams at higher altitudes, the assumption that l∆y L v is more likely to be violated. Airborne in-situ measurements are the best way to validate the assumptions and the lidar retrievals in these cases.
The CoMet dataset was also used to show that with VAD-scans with larger elevation angle (here: 75 • ) can be used to retrieve TKE dissipation rate with the same method as for VAD-scans with 35.3 • and the results are comparable. For this narrow VAD cone scans, we showed that the advection correction is much more important than for lower elevation angles and strong scans with different elevation angle can in future potentially help to analyze horizontal heterogeneity in the boundary layer and its impact on the calculation of area-averaged fluxes.
Data availability. The data are available from the author upon request.

10
A partial decomposition of Eq. A1 yields: u 2 + v 2 + 2 tan 2 ϕ w 2 = 1 π cos 2 ϕ These equations provide the basis for the retrieval of TKE and momentum fluxes from lidar VAD measurements.

24
Appendix B: Lidar filter functions Theoretical models for the spectral broadening of lidar measurements (F (∆y)) and the structure function(A(l∆y)) are derived in Banakh and Smalikho (2013) from the two-dimensional Kolmogorov spectrum for lidar measurements of turbulence: Θ(κ z , κ y ) = C 3 (κ 2 z + κ 2 y ) −4/3 1 + 8 3 · κ 2 y κ 2 z + κ 2 y (B1) where H is the longitudinal and H ⊥ the transverse filter function of lidar measurements in the VAD scan: H (κ 1 ) = exp −(π∆pκ 1 ) 2 sin (π∆Rκ 1 ) π∆Rκ 1 2 (B4) where ∆p is derived from the FWHM pulse width τ p , ∆R from the time window T w and ∆y from the VAD azimuth increment For statistically independent subsamples X j with size k j , the total variance of the dataset can be derived as follows:
where n = k j . Eventually, for equally sized subsamples one obtains: Appendix D: Validation of wind retrieval The FSWF-retrieval is used to obtain the three-dimensional wind vector from the lidar VAD scans. The results of the retrieval 10 is compared to the sonic anemometers at 50 m and 90 m and shown in Fig. D1. To show the distortion of the mast, no data have been removed in the retrieval of wind speed and wind direction for this figure.
26 Figure D1. Scatter plot of horizontal wind speed retrieved from lidar measurements compared to sonic anemometer measurements at 50 m (a) and 90 m (b) and wind direction (c and d). 27