Interactive comment on “ Comparison of aerosol LIDAR retrieval methods for boundary layer height detection using ceilometer backscatter data

In this paper, three algorithms for estimating the boundary layer height (BLH) with a Vaisala CL31 ceilometer are assessed: the manufacturer’s aerosol gradient method, the cluster analysis method of Toledo et al. (2014), and the Haar-wavelet method, described e.g. in Cohn and Angevine (2000) or Davis et al. (2000), as well as in many publications later on. The ceilometer is located near the urban area of Houston, Texas, and the ceilometer BLH retrievals are compared with BLH retrievals from collocated radiosonde measurements. In total, 48 radiosonde launches over more than 4 years are considered, all during daytime (mostly in the early afternoon) and cloud-free periods.

1. Introduction p.2, l.12 ", facilitating the monitoring of the nocturnal stable layer, internal aerosol layers and the nighttime residual layer" What about the ML? p.2, l.14 replace: "in the need for determining the most reliable and accurate method" 2. Data an Instrumentation p.2, l. 25-26 "Over 80 radiosonde launches and CL31 backscatter profiles beginning in January 2011 through March2015 are used." Why not say 85 (see Sect. 2.2) directly?Also, this gives me a launch every 19 days during this time period.Please me more specific (regular launches (i.e.every 19 days)?concentrated on some periods, some seasons?) p.2, l.27 replace: "The effect of cloud signals is analyzed separately for each method" p.2, l.27-28 "this data set includes data from the NASA DISCOVERAQ Texas campaign" What kind of data?Radiosonde measurements?Please specify.

Vaisala CL31
p.3, l.4 Please be more precise: "The backscatter coefficient \beta or the backscattering cross section per unit volume of air is related to the received power with the following formula:" I prefer the formulation "is related to the received power" instead of "can be calculated from", because calculating \beta from eq (1) is not straightforward.The CL31 gives you a (scaled) attenuated backscatter coefficient (i.e.C_Vaisala \beta \tauˆ2) in 10ˆ-9 C4 mˆ-1sr-1 units, with C_Vaisala around 1. I accept that you call this aerosol backscatter or simply backscatter throughout the paper, as long as you remain consistent.Side comment: the liquid cloud calibration method of O'Connor et al. ( 2004) is a method that is used to calculate C_Vaisala.
p.3, l.6 I am missing the background term (i.e.P(x,\lambda) = . . .+ B) in eq. ( 1) and in its description in lines 7-10.Please add.Please be more specific.Add how many launches happened e.g. between 12:30 and 13:30 (or some other period that you define around 13:00 CST), and how many after that midday period.
3. Boundary Layer Height Retrieval Methods p.3, l.28 remove "therefore a constant concentration of backscatter" Indeed, "concentration of backscatter" does not make sense.Also, the growth of aerosols due to swelling does not change their concentration, but changes the C5 backscattered signal (i.e.constant concentration does not mean necessarily constant backscatter).
p.4, l.1-2 "within 10 minutes before or after radiosonde launch." Please be more specific.Do you choose the closest in time in the 20min time-window?the mean aerosol-BLH in the 20min time-window?the closest in altitude in the 20min time-window?Or else?
3.1 Skew-T log-P Diagram for Radiosonde Boundary Layer Heights p.4, l.5 replace: "an unstable boundary layer is identified by having a dry adiabatic lapse rate greater than " p.4, l.8 replace: "where the temperature profile" p.4, l.9 replace: "where relative humidity or dew point temperature profiles sharply decrease as seen in the skew-T log-P diagram in Fig. 1b" p.4, l.11 be more specific: "when comparing ceilometer and radiosonde derived BLHs (both manually) using " Did you retrieve the BLH from the skew-T log-P diagram manually for each day?Please specify.Please comment on this quantitative difference in the paper, in Sect.4.1 for example.
I actually do not expect a correlation coefficient of much more than 0.9 when comparing aerosol gradient BLH with radiosonde BLH, especially if the ML is not fully developed, since aerosol gradients are only a proxy for the thermodynamic BLH, and are not expected to correlate perfectly with it.
3.2 Vaisala Corporation Aerosol Backscatter Gradient p.4, l.18 replace: "The temperature correction of -10 is an algorithm setting that adjusts the shape" p.4, l.24 "Therefore, a manual analysis of the algorithm's three resulting layers (Fig. 2) is required in order to prevent the incorrect identification of other aerosol layers".
Did you do this manual analysis in this study?Please specify, since not entirely clear here.
p.4, l.25 "The algorithm gives three maximum negative gradients every 1-minute, these are averaged to 10 minutes for radiosonde comparison".
How are they averaged?Layer by layer?Or by clustering all points inside the 10 minutes into 3 clusters/layers, and then taking the mean of these clusters?Or else?
3.3 Cluster Analysis p.4, l.30-31 be more specific "The BLH is typically identified as the (temporal) variance local maximum" p.5, l.8 replace: "the greater the EZ depth the greater the overestimation of the BLH" p.5, l.9 punctuation, missing point at the end: "whereas aerosol gradient methods can give multiple results."

C7
Comment: the cluster method gives a unique BLH only if you choose the number of clusters to be k=2.

Data Processing for Cluster Analysis and Application
p.5, l.14 replace "Due to increasing noise in backscatter profiles with height" with "Because the range correction needed to invert Eq. 1 increases noise in backscatter profiles with height, " p.5, l.17 be more specific: "The moving time average" p.5, l.21In Eq. ( 2), the power 2 should be applied on the brackets, not on the averaged profile.
Also, you already used P as the received power in Eq. ( 1).Use an another letter for this variable.
p.5, l.22 try to add space and correct: "where P(z,space ti) is the averaged LIDAR backscattered signal at time ti and height z, and P is the averaged profile from" p.5, l.22-23 "from N number of profiles" How much is N here?The number of profiles corresponding to 10min?Please specify.
p.5, l.24 "K-means is a data-partitioning algorithm that assigns observations" Explain that the observations are 3D-points where the first dimension is the range, the second the backscattered signal and the third the variance, and that these 3D-points are standardized (Toledo et al. ( 2014)).
p.5, l.30 punctuation, missing point at the end: "Step 4. Compute the average of the observations in each cluster to obtain new centroid locations." p.6, l.3 replace "Cluster analysis will typically divide a well-mixed boundary layer into two clusters, one below a peak in variance corresponding the center of the EZ," with C8 "By choosing k=2, cluster analysis will typically divide a well-mixed boundary layer into two clusters, one below a peak in variance corresponding to the center of the EZ," Indeed, the cluster method gives two clusters only if you choose the number of clusters to be k=2.
p.6, l.5 "will cause the cluster analysis to assign clusters using other criteria" I do not understand what you mean with "using other criteria".Replace with "will cause the cluster analysis to assign clusters somewhere else" Indeed, there may be more than on big aerosol gradient (e.g.top of BL and top of lofted aerosol layer).C9 p.6, l. 24-26 change the order of the 2 sentences, i.e. change: "Lower dilation values create numerous CWTC local minimums (Fig. 4b) at heights of smaller aerosol gradients in the measured profiles.Larger values create large local minimums (Fig. 4c and 4d) at the height of the biggest aerosol gradients in the backscatter profile (Fig. 4a)." with "Larger values create large local minimums (Fig. 4c and 4d) at the heights of the biggest aerosol gradients in the backscatter profile (Fig. 4a).Lower dilation values create numerous CWTC local minimums (Fig. 4b) at heights of also smaller aerosol gradients in the measured profiles." The reason is that at the location of a strong negative gradient, the CWTC with a small dilation factor will also have a (strong) local minimum there.
This is why, in some studies, for robustness reasons, the local minima over the averaged CWTC profile (averaged over multiple dilations, say here 30m to 300m) are searched for.Indeed, the gradient you seek at the top of the BL, under the assumptions you stated in the beginning of Section 3, has a peak in the Wavelet transform at multiple dilations.See for e.g.Cohn and Angevine (2000) Fig. 2 b).
As you already mentioned, going to high dilations however does not apply when seeking the top of the NSL, where only small dilations should be used.
p.6, l.29-30 be more specific, i.e. replace "A higher value of 300m (Fig. 4d) for the dilution factor a is applied for daytime BLHs to identify the sharp transition between ML and FT" with "A higher value of 300m (Fig. 4d) for the dilution factor a is applied for daytime BLHs and the strongest CWTC local minimum is used to identify the sharp transition between ML and FT" C10 p.7, l.7 replace "t test" with "t-test" p.7, l.10 "where X is the mean of the aerosol BLH samples, \mu is the radiosonde BLHs mean, S is the standard deviation of samples," p.7, l.13 replace "p value" with "p-value" p.7, l.13 replace "or 95% confidence" with "or 5% significance level" p.7, l.20 replace: "following precipitation or during periods of high wind speeds." p.8, l.1 "(not statistically significant; p < 0.05)." Not statistically significant means p > 0.05, please recheck this and all other statements you make about statistically significance.
p.8, l.5-6 replace: "when the algorithm did not find strong enough gradients in the backscatter profile" p.8, l.8-9 reformulate, i.e. replace "This is due to the assumption in the methodology of using aerosol gradients to detect the BLHs and thermal parameters to detect radiosonde BLHs." with "This is due to the difference of assumptions in the methodologies, using aerosol gradients to detect the BLHs on one side and thermal parameters to detect radiosonde BLHs on the other side."

Cluster Method Results
p.8, l.12 What about the statistical significance for this method?
p.8, l.13-14 "From the 45 comparisons performed, 13.3% showed the algorithm mistakenly finding maximum variance peaks not corresponding not corresponding to the C11 BL but to noise or other aerosol layers" You claim later in this paragraph 16 cases with noise and 5 cases with lofted aerosol layers, which makes (16+5)/45=46.7%and not 13.3%.Please change accordingly.
p.8 l.19, replace "CL31 displays a significant increase in noise with increasing altitude." Note that all lidars in general, not only the CL31, display an increase in noise with height, due to the range-correction.
p.8, l.21-22 "These were not algorithm errors but instead due to the implicit assumptions in using aerosol backscatter for BLH detection (constant concentration of backscatter within the BL and a negative gradient in backscatter corresponds to the top of the BL)".
Change "constant concentration of backscatter" by "constant backscatter" (see my comment on p.3, l.28)I am not satisfied with the formulation: "these were not algorithm errors but instead due to the implicit assumptions".I trust that you did not make a programming error.Still, these errors are due to limitations of the cluster algorithm, and not due to violations of your basic assumptions.
Indeed, assume you have constant backscatter within the BL and a negative gradient in backscatter at the top of the BL.Now suppose there is a lofted aerosol layer on top of the BL and that your measured backscatter profile is affected by noise.This does not violate your basic assumptions of constant backscatter in the BL and gradient at the top, but still perturbs the algorithm.
The cluster method works best to retrieve the top of the residual layer or of the fully developed mixed layer, assuming no clouds, high signal-to-noise ratio and no lofted aerosol layers above the BL.The noise aspect in particular, affecting the temporal variance, seems to have been an important issue here, and you could for e.g. say that this issue should be kept in mind and tackled when using the cluster-or other C12 variance-based algorithms on the CL31 (even if you say that later in the conclusions, p.10, l.22).

Wavelet Method Results
p.9, l.4 Add: "capturing the BLH.Also, it was less affected by noise than the gradient method or the cluster method." Remark: These conclusions are not surprising, since the wavelet method with 300m dilation means taking the difference of two 150m-averaged backscatter signals.With this, you increase the total amount of signal considered and decrease the variability.The drawback is that with a big dilation of 300m, you risk to detect the top of a lofted aerosol layer instead of the top of the BL, which happens 12.5% of the cases, as you mentioned.
4.4 Retrieval with cloud signals p.9, l.6 soften the appreciation: "The identification of the BLH is more difficult in the presence of clouds".
Indeed, you say yourself in p.9, l.13-15 that in case of cloud the gradient and wavelet methods often compare well with the thermodynamic BLH.
p.9, l.6-7What about precipitation or fog events?Are they included into this category of BLH retrieval with cloud signals?Please specify.
p.9, l.26 "The presence of clouds limits the detection of the BLH for all methods due to either the extinction of aerosol backscatter signals above the cloud, or the presence of low clouds mistakenly identified as the BLH." In case of cumulus or stratocumulus, even if you have extinction above the cloud, the error you make with your ceilometer BLH detection compared to the thermodynamic BLH is acceptable and not a major limitation.The presence of low clouds mistakenly identified as the BLH is similar to the cases where you have a lofted aerosol layer above C13 the BL, in the sense that the algorithms get "distracted" by those (strong) gradients and or variances and fail to pick the correct gradient at the top of the BLH.This is an attribution problem (linked to the processing of the measurement), but not a problem of the ceilometer measurements.
p.9, l.27-28 soften the appreciation: "Hence the removal of profiles with cloud signals is preferred for the automatic retrieval of the BLH".
The removal of profiles with cloud signals is per se not needed for the automatic retrieval of the BLH.Clouds introduce simply additional uncertainty.The algorithms you use in this paper do not really tackle the attribution problem: retrieved BLHs are here simply defined as the strongest gradient in the profile, or the first gradient from the ground, or the first separation point between two clusters, etc.There are more robust algorithms that have been developed in the last few years and advances have been made in this attribution problem issue, see my comments below for the "conclusions" section.
p.9, l.29 be more specific: "since the moving time averaging performed before the application" 5 Summary and Conclusions p.10, l.2-3 be more specific: "were compared to over 40 cloud-free daytime radiosonde profiles" p.10, l.9 replace: "corresponding to the top of the BL by calculating the wavelet transform" p.10, l.20 "Profiles were also mistakenly divided into relative variance concentrations rather than a peak in variance," I do not quite understand the term "relative variance concentrations".Replace with: "Profiles were also mistakenly divided due to the increasing noise with height rather C14 than a peak in variance," p.10, l.21-22 "Previous knowledge of the BL will aid in identifying these algorithm errors, but is not necessary to obtain a BLH with this method." Not very clear.Replace with: "This method being automatic, previous knowledge of the BL aids in identifying such algorithm errors, but is otherwise not necessary."p.10, l.24 "without previous knowledge of the BL required to derive at a BLH." Not very clear.Replace with: "without previous knowledge of the BL required, as this method is also automatic." p.11, l.1 replace: "signal and a cluster for cloud-free signal", i.e. signal without ending s, because it is along one backscatter profile.
p.11, l.5 be more specific: "more sensitive due to the moving time averaging applied" p.11, l.3-5 "Limited detection of the BLH in aerosol profiles with cloud signals is seen for all methods due to either the extinction of aerosol backscatter signals above the cloud, or the presence of low clouds mistakenly identified as the BLH" You did not really show this "limited detection of the BLH in aerosol profiles with cloud signals".(85-48)/85 = 43.5% of the measurements were not analyzed because of cloud presence, which is not negligible.I would be very interested to see a scatter plot exactly like the one in Figure 5, but for the 37 days with clouds.You said yourself in Section 4.4 that the methods tested here, especially the gradient and the wavelet methods, correlated quite well with the skew-T log-P derived BLH, so why actually exclude them from the analysis?p.11, l.8-10 "The errors found with this method were due to the basic assumptions used to derive BLH from aerosol backscatter concentrations rather than errors with the C15 algorithm itself." Again, as in p.8, l.21-22,I am not satisfied with this formulation.Replace with: "The errors found with this method were due to lofted aerosol layers, low-level clouds and differences in determining BLHs using aerosols and thermodynamically using radiosondes." p.11, l.10-13 "For this reason, the use of this automated method is recommended for BLH observations with careful determination of the dilation values needed for individual instrument and locations.It is suggested to employ the wavelet method in future studies, in particular for long-term boundary layer evolution studies and spatial analysis of the BL using multiple LIDAR aerosol backscatter measurements." More results are needed in order to accept these two phrases, as I will explain below.Please reformulate: "In order to use this method on other instruments and locations, dilation values should be determined carefully and individually.Out of the three methods tested in this study, it is suggested to employ the wavelet method in future studies, in particular for long-term boundary layer evolution studies and spatial analysis of the BL using multiple LIDAR aerosol backscatter measurements." What do you mean by "evolution" in "boundary layer evolution studies"?Seasonal diurnal cycle?Annual cycle of the e.g.13:00 CST value?Please specify.
Here you did not show any "climatological results" (i. ).Such a more robust method has not been investigated or commented here, and your suggestion to use the Haar-wavelet method (out of all existing methods) is thus not justified enough.
Note that there is a review study of (Haeffelin et al., 2012), where several state-ofthe-art algorithms at that time on three different ceilometers, including the CL31, are compared.
Here you did not do a review of all state-of-the-art algorithms, but compared three rather simple methods (one manual, two automatic) with radiosonde measurements, and showed that the wavelet method with a large dilation is a robust method on the CL31 to derive the height of strong gradients in the aerosol backscatter profiles, and that using the temporal variance should be done with precaution on this ceilometer.This is valuable information for future work.You also showed that there are cases C17 where the aerosol layer based and thermodynamic based PBLHs do not match perfectly, which would be a further example giving input to the discussion on how well both methods compare.To my knowledge, it is also the first time the Cluster method, which is based on a good idea but seems quite sensitive to noise, is evaluated with the CL31 ceilometer, and you showed also that it was a potentially interesting technique to detect cloud layers.
Remove the word "Density" in the z-label "Aerosol Backscatter Density" of the figure, because "backscatter density" does not make sense and it will be more consistent with the "Aerosol Backscatter" that you write in all other plots.
Change the legend: "Aerosol backscatter time series for September 26, 2013." Interactive comment on Atmos.Meas.Tech.Discuss., doi:10.5194/amt-2016-340,2016. C18 p.3, l.15 Please specify which version of the firmware of the CL31 you use and what noise_h2 setting do you use.With the noise_h2 setting turned off, the CL31 backscat- sure (at least I need to see it) that you are able to see the full diurnal cycle of the ML with the wavelet method as implemented here.For example, during the morning hours, by taking the strongest local minimum, you risk to still detect the top of the RL instead of the top of the developing CBL.But for the annual cycle of the CBL at say 13:00 CST (where you could roughly assume that the CBL is fully developed), I could imagine it works reasonably well.There are also other automatic methods than the wavelet method or the cluster method, or in a more general sense gradient or variance methods.I think in particular of the "fitting an idealized backscatter profile to observed profiles" method ofSteyn et al. (1999), which is a publication you cited.In addition, you might be confronted to physical inconsistencies if you only take the largest peaks in gradient or variance (especially during the development of the CBL in the morning, you risk to jump back and forth between the still existing RL and the CBL).More robust algorithms aiming to tackle these issues exist (among others, gradient-variance coupling with temporal height tracking (Martucci et al. (2010), coupling with a BL-model (di Giuseppe et al. (2012), Biavati (dissert., 2014)), coupling with surface measurements (Pal et al. (2013)), coupling with graph theory (de Bruine et al. (AMTD, 2016)) e. do we recover a diurnal and annual cycle with this automatic wavelet method and do they compare qualitatively well with the cycles depicted in theHaman et al. (2012)study?)Thiswouldbe the expected second part of the validation of the algorithm (wavelet method) you preconize in your conclusions, now that you have good comparison results with the radiosonde measurements, and especially if you suggest this algorithm as suitable for PBL studies on multiple LIDARs.C16I am not