Analyzing the turbulent Planetary Boundary Layer by remote sensing systems : Doppler wind lidar , aerosol elastic lidar and microwave radiometer by

14 The Planetary Boundary Layer (PBL) is the lowermost region of troposphere and endowed with turbulent 15 characteristics, which can have mechanical and/or thermodynamic origins. Such behavior gives to this layer 16 great importance, mainly in studies about pollutant dispersion and weather forecasting. However, the 17 instruments usually applied in studies about turbulence in the PBL have limitations in spatial resolution 18 (anemometer towers) or temporal resolution (instrumentation onboard aircraft). Ground-based remote 19 sensing, both active and passive, offers an alternative for studying the PBL. In this study we show the 20 capabilities of combining different remote sensing systems (microwave radiometer [MWR], Doppler lidar 21 [DL] and elastic lidar [EL]) for retrieving a detailed picture on the PBL turbulent features. The statistical 22 moments of the high frequency distributions of the vertical wind velocity, derived from DL and of the 23 backscattered coefficient derived from EL, are corrected by two methodologies, namely first lag and -2/3 24 correction. The corrected profiles, obtained from DL data, present small differences when compared against 25 the uncorrected profiles, showing the low influence of noise and the viability of the proposed methodology. 26 Concerning EL, in addition to analyze the influence of noise, we explore the use of different wavelengths 27 that usually include EL systems operated in extended networks, like EARLINET, LALINET, MPLNET or 28 SKYNET. In this way we want to show the feasibility of extending the capability of existing monitoring 29 networks without strong investments or changes in their measurements protocols. Two case studies were 30 analyzed in detail, one corresponding to a well-defined PBL and another one corresponding to a situation 31 with presence of a Saharan dust lofted aerosol layer and clouds. In both cases we discuss results provided 32 by the different instruments showing their complementarity and the cautions to be applied in the data 33 interpretation. Our study shows that the use of EL at 532nm requires a careful correction of the signal using 34 the first lag time correction in order to get reliable turbulence information on the PBL. 35

(Lines 96 -101) "One of the goals is to show the feasibility of using  at 532 nm, considering the widespread use of lidar systems based on laser emission at this wavelength in different coordinated networks, like as EARLINET (Pappalardo et al., 2014) and LALINET -Latin American LIdar Network (Guerrero-Rascado et al., 2016).In addition, this study shows the variety of application that can be done with EARLINET data applying some simple changes in the data acquisition procedures."

General comments
This manuscript presents results from the SLOPE campaign in Granada, Spain, in which the objective was to obtain closure between remote sensing and in-situ measurements.For this manuscript, the focus is on characterizing the planetary boundary layer using a Doppler lidar, multi-wavelength lidar (MULHACEN), and a profiling microwave radiometer, all operating at high temporal resolution (2 seconds).The authors investigate the use of fluctuations in aerosol number density from the elastic system (EL), vertical velocity fluctuations obtained from the Doppler lidar (DL), and potential temperature profiles retrieved from the microwave radiometer (MWR), to identify the boundary layer height (PBLH).As stated in the first and second review, the methodology is relevant, and the influence of random error introducing extra noise in higher-order moments has merit and is explored using suitable techniques.
The manuscript has now been improved significantly, with a clearer focus, and explores the impact of applying the elastic lidar methodology at different wavelengths.New Figures 9 and 10 now show that, although backscattering coefficients are wavelengthdependent (molecular and aerosol), the methodology and correction procedure can account for these differences.It is now clear that the higher moments exhibit more correction at 532 nm but the methodology used to derive PBLH from elastic lidar is not unduly sensitive to the wavelength used, at least.
The new supplementary figures are much better, and clearly display where reliable data may be obtained.However, there are still a few issues for the authors to address before the manuscript is suitable for publication.This question was asked previously: "The EL and DL parameters are calculated over 1hour periods.Is this 1-hour timescale suitable during rapidly varying conditions such as during the morning growth of the boundary layer?"This question is asking whether a 1 hour timescale is suitable when, during the morning growth, a particular region may have been calm for 30 minutes, and then strongly turbulent for 30 minutes.What happens to the integral timescale for both EL and DL properties when you include atmospheric regions with very different turbulent attributes (e.g calm and convective) within the same averaging period?With one hour averaging, you will always miss the rapid growth of the CBL.I understand you want to capture the integral time scale, but how can you do this when the turbulence characteristics themselves are changing?
We thank the Reviewer for this comment.Answer Figure 1-A (case study I of the main document) shows a situation with a fast growth of the PBL during 30 min followed by 30 min interval with a PBLH almost constant.In figures 1-B and 1-C we show the analyses of this case split in two 30 min intervals corresponding to the two different situations, fast growth and almost constant PBLH.In both cases the integral time scale is lower than that computed for the whole 1 hour interval and the profiles of Skewness and Kurtosis are rather noisy, thus complicating the observation of determined phenomena, although the profiles are very similar.So it is evident that the analyses of intervals bellow 1 hour are not good enough, the degradation of the analyses increases with the reduction of the considered interval as can be seen in the figures (15 min [Fig.1-D], 10 min [Fig.1-E] and 5 min [Fig.1-F]).In this sense, it is evident that the features of our equipment does not allow the detection, with appropriate quality, of turbulent events with a temporal resolution lower than 1 hour.So when the PBL present faster changes we can only observe the average behavior of the turbulence with this time window.

Minor comments
Lines 42-47.Please check and reformulate these sentences.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 42 -48) "In an ideal situation, some instants after sunrise, the ground surface temperature increases due to the positive net radiative flux (  ).This process intensifies the convection, where there is an ascension of warm air masses, causing the downward displacement of colder air masses and consequently originating the Convective Boundary Layer (CBL) or Mixing Layer (ML).Such layer has this name due to the mixing process generated by the ascending air parcels.Slightly before sunset, the gradual reduction of incoming solar irradiance at the Earth's surface causes the decrease of the positive   and, consequently, its sign change.In this situation, there is a reduction of the convective processes and a weakening of the turbulence." Line 78: Replace 'turbulenc' with 'turbulence'.

Done
Line 86: The convective PBL is the CBL.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Line 87) "…  height () during the convective period …" Line 96: Replace 'realibility' with 'reliability'.Explain why measurements at 532 nm should be more reliable, or suggest removing this sentence.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 97 -99) "… considering the widespread use of lidar systems based on laser emission at this wavelength in different coordinated networks, like as EARLINET (Pappalardo et al., 2014) and LALINET -Latin American LIdar Network (Guerrero-Rascado et al., 2016)." Lines 129-130: Suggest replacing 'This system record the backscattered signal with 300 gates, being the range gate length 30 m, with the first gate at 60 m' with 'This system records the backscattered signal with a range resolution of 30 m in 300 range gates with the first range gate starting at 60 m from the instrument.' Done Line 131: Suggest stating 'The instrument was operated in vertical stare mode with a temporal resolution of 2 s'.Using the phrase 'with respect to the ground surface' could mean that, on a sloping surface, you imply you are pointing normal (90 degrees) to the surface.Pointing vertically is unambiguous and doesn't require the qualifier.
We thank the Reviewer for this comment.In order to clarify this point the suggested change has been applied.Line 168: Replace 'performed campaign of comparison' with 'intercomparison campaign'.

Done
Line 173: Replace 'allow to estimate the CBL height' with 'allows the estimation of the CBL height'.
We thanks the Reviewer for this comment.In order to clarify this point the text has been changed as follow:

Done
Lines 242-244: Suggest rewriting this phrase as it is unclear.You could use ' in this study we evaluate using RCS532 fluctuations to determine turbulence following the procedure described in Figure 3.This EL methodology is very similar to that described earlier for DL.'

Done
Line 257 and elsewhere: Suggest using 'below the PBLH' rather than 'under the PBLH'.

Done in line 264 and line 270
Line 259: Replace 'relation the other' with 'relation to the other'.

Done
Lines 263-265.This statement is only true in high SNR conditions -figure 5 only shows data within 1200 m of the instrument.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 270-271) "…Therefore, considering high Signal-to-Noise Ratio () conditions, although the presence of ..." Lines 272-273: Suggest replacing 'its spread use in observation network with higher reliability than 1064 nm' with 'its widespread use in observation networks'.
We thank the Reviewer for this comment.In order to clarify this point the suggested change has been applied.
Lines 274-275: Suggest replacing 'As expected, in both cases the increase of height produces the increase of ε,' with 'As expected, ε increases with range'.

Done
Lines 278-283: The difference in noise levels between the two wavelengths depends on the SNR at each wavelength, which is more likely to be determined by the laser output power, filters, and detectors used, at the two wavelengths.Higher molecular extinction at 532 nm can then reduce SNR relative to the 1064 nm wavelength, as does separation of the molecular and aerosol backscattering.Figure 8 is not necessary for the manuscript.
We thank the Reviewer for this comment.In order to clarify this point the figure 8 has been removed and the text has been changed as follow: (Lines 286-292) "Although the level of influence of  in each wavelength depends on the  of them (which is associated to technical factors such as laser output power, filters, type of detectors), considering the proposed methodology, to evaluate the composition of each wavelength is also important.The large contribution of   532 to the total  at 532 nm in comparison with the behavior at 1064 nm, can influence the results obtained from such wavelength, because our methodology is based on the use of  ′  .In addition, the larger extinction (due to both aerosol particles and molecules) at 532 nm produces a lower two-way transmittance, resulting in the reduction of the  values at this wavelength." Line 299: SNR is reduced, the noise doesn't increase.
We thank the Reviewer for this comment.In order to clarify this point the figure 8 has been removed and the text has been changed as follow: (Lines 308-310) "…which reduces the  of the profiles in comparison with 1064 nm, the application of the proposed corrections, mainly the first lag, reduces significantly such influence and…" Lines 303-304: Suggest rewriting this phrase as it is unclear.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 312-313) "The first lag correction was adopted as default because it provides better results than the -2/3 law correction." Lines 315-318: It is not clear that the integral time scale can be retrieved throughout the whole PBL, especially if PBLH_MWR is taken as a reference.It is not necessarily true that the grey areas are where T_w' is lower than the DL acquisition time, just that the DL sensitivity is not high enough to measure T_w'.This can be seen in the supplementary material, where SNR is low above 1500 m and the upper portion of the CBL may not be captured during daytime.I would expect T_w' to be just as large here.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 322-325) "Figure 10 (A) shows the integral time scale obtained from  data ( ′ ).The gray area represents the region where it is not possible to analyze the turbulent process using our  data, either because of the low  values, which results in null values of the  ′ , or because the no null  ′ is smaller than the acquisition time of the .However, the gray area is located almost entirely above the   (white stars)." Lines 319-320: This sentence can be removed.

Done
Lines 327-328: I suggest removing this sentence.

Done
Lines 329-320: Not correct, skewness describes the distribution of the turbulent velocities -positive skewness implies strong but narrow updrafts surrounded by weaker but more widespread downdrafts, and vice versa for negative skewness.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: Lines (333-335) "The skewness of ′ ( ′ ) is shown in Figure 11-C.The  ′ describes the distribution of the turbulent velocities.Thus positive  ′ implies strong but narrow updrafts surrounded by weaker but more widespread downdrafts, and vice versa for negative  ′ ." Line 336: Do you mean that that the MWR method is now selecting for SBLH?
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 341-342) "… Thus, the reduction observed in the   is due to the detection of  height." Lines 344-353: Suggest removing these paragraphs.Some of these statements could replace phrases in lines 329-337.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 349-351) "Figure 10-E presents the values of surface air temperature and surface relative humidity ().Air surface temperature has a daily pattern similar to that of   and  ′ .On the other hand,  is inversely correlated with the temperature." Lines 365-367: The PBLH values in Figure 13 don't agree with Figure 12, where there is a large difference between PBLH_MWR and PBLH_Elastic.
We thank the Reviewer for this comment.In order to clarify this point the figures 12 and 13 have been remade as shwon below: The text has been changed as follow: (Lines 366-374) "Due to presence of a decoupled aerosol layer at 13:30, the average values of   and   have a difference of around 500 m.The  ′ 2 has small and practically constant values between 1000 and 1400m, evidencing the homogeneity of the aerosol distribution in this region.Starting at 1400 m the value of  ′ 2 begins to increase, reaching a positive peak at   , which represents the Entrainment Zone (region characterized by an intense mixing between air parcels coming from  and Free Troposphere (), causing a high variation in aerosol concentration).The   observed at approximately 2900 m demonstrate an inherent difficulty of the variance method to detect the PBLH in the presence of several aerosol layers (Kovalev and Eichinger, 2004).Above   the values of  ′  We thank the Reviewer for this comment.In order to clarify this mistake the text has been changed as follow: (Lines 392 -393) "The results provided by , pyranometer and  data agree with the results observed in figures 12 and 13." Line 399.Suggest removing the second phrase of this sentence.

Done
Line 400: Not true according to the figure.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 399-400) "…the greatest part of grey area is situated above the   ..." Lines 411-414: No clouds are observed in the DL data until 1400 UTC, and it is difficult to prove that the negative skewness extends to 3 km in altitude.Are you sure that all of the white regions are cloud before 1400?One at 1230 and one at 1330 maybe.
We thank the Reviewer for this comment.During this period there is the presence of both middle altitude clouds and very intense dust layers.In order to clarify this point the text has been changed as follow: (Lines 411 -415) "From Figure 16 we can observe the presence of both middle altitude clouds and very intense dust layers from 12:00 to 15:00 UTC.Such combination contributes to the intense negative values of  ′ observed in this period until around 2 km, because, as mentioned previously,  ′ is directly associated with the direction of the turbulent movements.The present situation can be considered representative of cloud-top long-wave radiative cooling in the CBL (Ansmann et al., 2010)." Lines 415-420.As above, it is clear that the Saharan dust layer is having an impact, but Rn alone is probably not sufficient to attribute negative skewness directly to clouds.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 417-420) "The observation of  ′ and   between 12:00 and 14:00, as well as, the presence of clouds and geometrically thick dust layers during this same period, reinforces the hypothesis that we have a situation of cloud-top long-wave radiative cooling in the ." Line 423: See above comment.Not all of the high RCS values (white regions) before 1300 can be attributed to clouds, and there is very little attenuation seen in the profile before 1230.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Line 425 -426)

"…the presence of both middle altitude clouds and very intense dust layers …"
Line 463: This is only true at high SNR.
We thank the Reviewer for this comment.In order to clarify this point the text has been changed as follow: (Lines 467-468) "…low influence of the noise in high  conditions." Line 467: Replace 'lower two ways' with 'reduced two-way'.

Done 1 Introduction
The Planetary Boundary Layer () is the atmospheric layer directly influenced by the Earth's surface that responds to its changes within time scales around an hour (Stull, 1988).Such layer is located at the lowermost region of troposphere, and is mainly characterized by turbulent processes and a daily evolution cycle.In an ideal situation, some instants after sunrise, the ground surface temperature increases due to the positive net radiative flux (  ).This process intensifies the convection, where there is an ascension of warm air masses, causing the downward displacement of colder air masses and consequently originating the Convective Boundary Layer (CBL) or Mixing Layer (ML).Such layer has this name due to the mixing process generated by the ascending air parcels.Slightly before sunset, the gradual reduction of incoming solar irradiance at the Earth's surface causes the decrease of the positive   and, consequently, its sign change.In this situation, there is a reduction of the convective processes and a weakening of the turbulence.
In this process the  leads to the development of two layers, namely a stably stratified boundary layer called Stable Boundary Layer () close to the surface, and the Residual Layer () that contains features from the previous day's  and is just above the .
Knowledge of the turbulent processes in the  is important in diverse studies, mainly for atmospheric modeling and pollutant dispersion, since turbulent mixing can be considered as the primary process by which aerosol particles and other scalars are transported vertically in atmosphere.Because turbulent processes are treated as nondeterministic, they are characterized and described by their statistical properties (high order statistical moments).When applied to atmospheric studies such analysis provide information about the field of turbulent fluctuation, as well as, a description of the mixing process in the  (Pal et al., 2010).
Anemometer towers have been widely applied in studies about turbulence (e.g., Kaimal and Gaynor, 1983;van Ulden and Wieringa, 1996), however the limited vertical range of these equipment restrict the analysis to regions close to surface.Aircraft have also been used in atmospheric turbulence studies (e.g., Lenschow et al., 1980;Williams and Hacker, 1992;Lenschow et al., 1994;Albrecht et al., 1995;Stull et al., 1997;Andrews et al., 2004;Vogelmann et al., 2012), nevertheless their short time window limits the analysis.In this scenario, systems with high spatial and temporal resolution and enough range are necessary in order to provide more detailed results along the day throughout the whole thickness of the .
In the last decades, lidar systems have been increasingly applied in this kind of study due to their large vertical range, high data acquisition rate and capability to detect several observed quantities such as vertical wind velocity [Doppler lidar] (e.g.Lenschow et al., 2000;Lothon et al., 2006;O'Connor et al., 2010), water vapor [Raman lidar and DIAL] (e.g.Wulfmeyer, 1999;Kiemle et al., 2007;Wulfmeyer et al., 2010;Turner et al., 2014;Muppa et al., 2015), temperature [rotational Raman lidar] (e.g.Behrendt et al., 2015) and aerosol [elastic lidar] (e.g.Pal et al., 2010;McNicholas et al., 2015).This allows the observation of a wide range of atmospheric processes.For example, Pal et al. (2010) demonstrated how the statistical analyses obtained from high-order moments of elastic lidar can provide information about aerosol plume dynamics in the  region.In addition, when different lidar systems operate synergistically, as for example in Engelmann et al. (2008), who combined elastic and Doppler lidar data, it is possible to identify very complex variables such as vertical particle flux.
Different works (Ansmann et al., 2010;O'Connor et al., 2010) have evidenced the feasibility for characterizing the  turbulence by .Pal et al. (2010) have shown the feasibility for retrieving information on the  turbulence from high high-order moments of elastic lidar operating at 1064.Such approaches are even more attractive when considering facilities of networks, e. g.European Aerosol Research Lidar NETwork (EARLINET) (Pappalardo et al., 2014), Microwave Radiometer Network (MWRNET) (Rose et al., 2005;Caumont et al., 2016) and ACTRIS CLOUDNET (Illingworth et al., 2007).
For these reasons, and having in mind the wide spread of elastic lidar systems operated at other wavelengths, like 532 nm or 355 nm, it would be worthy test the feasibility of these other wavelengths in the characterization of the  turbulent behavior.
The use of simple techniques, applied to the aforementioned remote systems provide robust and similar information on the  height () during the convective period (see for example Moreira et al, 2018), or a complementary information when the  is substituted by the presence of the  and the  (Moreira et al., in preparation).Thus, the combination of information obtained from the active remote sensing systems,  and , acquired with a temporal resolution close to 1 s, and that provided by  can provide a detailed understanding about different features of the , like structure ( versus  and ), height of the layers, rate of growth of the  and turbulence.
In this study we show the feasibility of obtaining a clear insight on the  behavior using a combination of active and passive remote sensing systems (Elastic Lidar [], Doppler Lidar [] and Microwave Radiometer []) acquired during the SLOPE-I campaign, held at IISTA-CEAMA (Andalusian Institute for Earth System Research, Granada, Spain) from May to August 2016.One of the goals is to show the feasibility of using  at 532 nm, considering the widespread use of lidar systems based on laser emission at this wavelength in different coordinated networks, like as EARLINET (Pappalardo et al., 2014) and LALINET -Latin American LIdar Network (Guerrero-Rascado et al., 2016).In addition, this study shows the variety of application that can be done with EARLINET data applying some simple changes in the data acquisition procedures.This paper is organized as follows.Description of the experimental site and the equipment setup are presented in Section 2. The methodologies applied are introduced in Section 3. Section 4 presents the results of the analyses using the different methodologies.Finally, conclusions are summarized in Section 5.

Experimental site and instrumentation
The SLOPE-I (Sierra nevada Lidar aerOsol Profiling Experiment) campaign was performed from May to September 2016 in South-Eastern Spain in the framework of the European Research Infrastructure for the observation of Aerosol, Clouds, and Trace gases (ACTRIS).The main objective of this campaign was to perform a closure study by comparing remote sensing system retrievals of atmospheric aerosol properties, using remote systems operating at the Andalusian Institute of Earth System Research (IISTA-CEAMA) and in-situ measurements operating at different altitudes in the Northern slope of Sierra Nevada, around 20 km away from IISTA-CEAMA (Bedoya-Velásquezet al., 2018;Román et al., 2018).The IISTA-CEAMA station is part of EARLINET (Pappalardo et al, 2014) since 2005 and at present is an ACTRIS station (http://actris2.nilu.no/).The research facilities are located at Granada, a medium size city in Southeastern Spain (Granada, 37.16°N, 3.61°W, 680 m a.s.l.), surrounded by mountains and with Mediterraneancontinental climate conditions that are responsible for cool winters and hot summers.Rain is scarce, especially from late spring to early autumn.Granada is affected by different kind of aerosol particles locally originated and medium-long range transported from Europe, Africa and North America (Lyamani et al., 2006;Guerrero-Rascado et al., 2008, 2009;Titos et al., 2012;Navas-Guzmán et al., 2013;Valenzuela et al., 2014, Ortiz-Amezcua et al, 2014, 2017).
MULHACÉN is a biaxial ground-based Raman lidar system operated at IISTA-CEAMA in the frame of EARLINET research network.This system operates with a pulsed Nd:YAG laser, frequency doubled and tripled by Potassium Dideuterium Phosphate crystals, emitting at wavelengths of 355, 532 and 1064 nm with output energies per pulse of 60, 65 and 110 mJ, respectively.MULHACÉN operates with three elastic channels: 355, 532 (parallel and perpendicular polarization) and 1064 nm and three Raman-shifted channels: 387 (from N2), 408 (from H2O) and 607 nm (from N2).MULHACÉN's overlap is complete at 90% between 520 and 820 m a.g.l. for all the wavelengths, reaching full overlap around 1220 m a.g.l.
(Navas- Guzmán et al ., 2011;Guerrero-Rascado et al. 2010).Calibration of the depolarization capabilities is done following Bravo-Aranda et al. (2013).This system was operated with a temporal and spatial resolution of 2 s and 7.5 m, respectively.More details can be found at Guerrero-Rascado et al. (2008, 2009).
The Doppler lidar (Halo Photonics, model Stream Line XR) is also operated at IISTA-CEAMA.This system works in continuous and automatic mode from May 2016.It operates at 1.5 µm with pulse energy and repetition rate of 100 µJ and 15 KHz, respectively.This system records the backscattered signal with a range resolution of 30 m in 300 range gates with the first range gate starting at 60 m from the instrument.The telescope focus is set to approximately 800 m.The instrument was operated in vertical stare mode with a temporal resolution of 2 s.Furthermore, we operated the ground-based passive microwave radiometer (RPG-HATPRO G2, Radiometer Physics GmbH), which is member of the MWRnet [http://cetemps.aquila.infn.it/mwrnet/].This system operates in automatic and continuous mode at IISTA-CEAMA since November 2011.The microwave radiometer (MWR) measures the sky brightness temperature with a radiometric resolution between 0.3 and 0.4 K root mean square error at 1 s integration time, using direct detection receivers within two bands: K-band (water vaporfrequencies: 22.24 GHz,23.04 GHz,23.84 GHz,25.44 GHz,26.24 GHz,27.84 GHz,31.4GHz) and V-band (oxygenfrequencies: 51.26 GHz,52.28 GHz,53.86 GHz,54.94 GHz,56.66 GHz,57.3 GHz,58.0GHz).From these bands is possible to obtain profiles of water vapor and temperature, respectively, by inversion algorithms described in Rose et al. (2005).The range resolution of these profiles vary between 10 and 200 m in the first 2 km and between 200 and 1000 m in the layer between 2 and 10 km (Navas-Guzmán et al., 2014).
The meteorological sensor (HMP60, Vaisala) is used to register the air surface temperature and surface relative humidity, with a temporal resolution of 1 minute.Relative humidity is monitored with an accuracy of ± 3%, and air surface temperature is acquired with an accuracy and precision of 0.6º C and 0.01º C, respectively.
A CM-11 pyranometer manufactured by Kipp&Zonen (Delft, The Netherlands) is also installed in the ground-based station.This equipment measures the shortwave (SW) solar global horizontal irradiance data (305-2800 nm).The CM-11 pyranometer complies with the specifications for the first-class WMO (World Meteorological Organization) classification of this instrument (resolution better than ±5 Wm −2 ), and the calibration factor stability has been periodically checked against a reference CM-11 pyranometer (Antón et. al, 2012).

MWR data analysis
The MWR data are analyzed combining two algorithms, Parcel Method [] (Holzworth, 1964) and Temperature Gradient Method [] (Coen, 2014), in order to estimate the  Height (  ) in convective and stable situations, respectively.The different situations are discriminated by comparing the surface potential temperature (( 0 )) with the corresponding vertical profile of () up to 5 km.Those cases where all the points in the vertical profile have values larger than ( 0 ) are labeled as stable, and  is applied.Otherwise the situation is labeled as unstable and the  is applied.The vertical profile of () is obtained from the vertical profile of (z) using the following equation (Stull, 2011): where () is the temperature profile provided by ,  is the height above the sea level, and 0.0098 K/m is the dry adiabatic temperature gradient.A meteorological station co-located with the  is used to detect the surface temperature [( 0 )].In order to reduce the noise, () profiles were averaged providing a   value at 30 minutes intervals.This methodology of detection was selected as the reference due to the results obtained during a performed intercomparison campaign between  and radiosonde data, where twenty-three radiosondes were launched.High correlations were found between  retrievals provided by both instruments in stable and unstable cases.Further details are given by Moreira et al. (2018a).

Lidar retrieval of the PBLH.
The simple processing of  and  data allows the estimation of the  height.Moreira et al. (2018), have discussed this issue in depth, while Moreira et al. (in preparation) have exploited the complementarity of the data obtained from distinct remote sensing systems in order to distinguish the sublayers during the period when the  and  substitute the , as well as, in complex situations, like as, presence of dust layers.
The  obtained from  data (  ) is estimated from variance threshold method.In this method the   is attributed to height where the variance of vertical wind speed (  2 ) is lower than a determinate threshold, which was adopted as 0.16 m²/s² (Moreira et al., 2018).For the   calculations was selected a time interval of 30 minutes.In concerning the  obtained from  (  ), the variance method is applied.Such method assumes the maximum of the variance of Range Corrected Signal (  2 ) as   (Moreira et al., 2015).The   2 is obtained from a time interval of 30 minutes.

Lidar turbulence analysis
Both lidar systems,  and , gathered data [(, )] with a temporal resolution of 2 seconds.Then, the data are averaged in 1-hour packages, from which the mean value is extracted [ ̅()].Such mean value is subtracted from each (, ) profile in order to estimate the vertical profile of the fluctuation for the measured variable [ ′ (, )] (i.e.vertical velocity for the ): Then, from  ′ (, ) is possible to obtain the high-order moments (variance (²), skewness () and kurtosis ()), as well as, the integral time scale ( -which is the time over which the turbulent process are highly correlated to itself) as shown in Table 1.These variables can also be obtained from the following autocovariance function,   : where   is the final time,  and  indicate the order of autocovariance function.
However, it is necessary to considerer that the acquired real data contain instrumental noise, ().
Therefore, the equation 3 can be rewritten as: The autocovariance function of a time series with zero lag results in the sum of the variances of the atmospheric variable and its ().Nevertheless, atmospheric fluctuations are correlated in time, but the () is random and uncorrelated with the atmospheric signal.Consequently, the noise is only associated with lag 0 (Fig. 1).Based on this concept Lenschow et al. (2000) suggested to obtain the corrected autocovariance function,  11 (→ 0), from two methods, namely first lag correction or -2/3 law correction.
In the first method,  11 (→ 0) is obtained directly by the subtraction of lag 0, ∆ 11 (0), from the autocovariance function,  11 (0).In the second method  11 (→ 0) is generated by the extrapolation of  11 (0) at firsts nonzero lags back to lag zero (-2/3 law correction).The extrapolation can be performed using the inertial subrange hypothesis, which is described by the following equation (Monin and Yaglom, 1979): where C represents a parameter of turbulent eddy dissipation rate.The high-order moments and  corrections and errors are shown in Table 1 (columns 2 and 3, respectively).
The same procedure of analysis is applied in studies with  and , being the main difference the tracer used by each system, which are the fluctuation of vertical wind speed (′) for  and aerosol number density (′) for .provides(, ) directly, and therefore the procedure described in Figure 2 can be directly applied.Thus, the two corrections described above are applied separately and finally  and highorder moments with and without corrections can be estimated.
On the other hand, the  does not provide (, ) directly.Under some restrictions, it is possible to ignore the particle hygroscopic growth and to assume that the vertical distribution of aerosol type does not changes with time, and to adopt the following relation (Pal et al., 2010): where   and  ′  represent the particle backscatter coefficient and its fluctuation, respectively, and () does not depend on time.
Considering the lidar equation: where   () is the signal returned from distance  at time ,  is the distance [m] from the lidar of the volume investigated in the atmosphere,  0 is the power of the emitted laser pulse,  is the light speed [m/s], is the duration of laser pulse [ns],  is the area [m²] of telescope cross section, () is the overlap function,   () is the total extinction coefficient (due to atmospheric particles and molecules) [(km) -1 ] at distance ,   () is the total backscatter coefficient (due to atmospheric particles and molecules) [(km•sr) - 1 ] at distance  and the subscript  represents the wavelength.The two path transmittance term related to () is considered as nearly negligible at 1064 nm (Pal et al., 2010).Thus, it is possible to affirm that: and consequently: where  1064 and  ′ 1064 are the range corrected signal and its fluctuation, respectively,  is a constant and the subscripts represent the wavelength.
In this way, Pal et al. (2010) have shown the feasibility of using  operating at 1064 nm for describing the atmospheric turbulence.However, having in mind the more extended use of lidar systems based on laser emission at 532 nm in different coordinated networks, e.g., in EARLINET and LALINET around 76% and 45% of the systems include the wavelength of 1064 nm, while 95% of the EARLINET systems and 73% of the LALINET systems operate systems that include the wavelength 532 nm (Guerrero-Rascado et al., 2016), in this study we evaluate using  532 fluctuations to determine turbulence following the procedure described in Figure 3.This  methodology is very similar to that described earlier for .
we perform the validation of the  532 in analyses about turbulence using , following the procedure described in Figure 3, which is basically the same methodology described earlier for .

Error Analysis
The influence of random error in noisy observations rapidly grows for higher-order moments (i.e., the influence of random noise is much larger for the fourth-order moment than for the third-order moment).
Therefore, the first step, in order to ascertain the applied methodology and our data quality, we performed the error treatment of  data as described in Figure 2.For the  analysis we selected the period 08-09 UTC of 19 th May, the same day that will be presented in Case Study 1.This day is characterized by a welldefined PBL.
Figure 4 illustrates the autocovariance function, generated from ′, at three different heights.As mentioned before, the lag 0 is contaminated by noise (), and thus the impact of the  increases together with height, mainly above   (1100 m a.g.l. in our example).
Figure 5-A illustrates the comparison between integral time scale ( ′ ) without correction and the two corrections cited in section 3.2.Except for the first height-bins, below the   the profiles have little differences, as well as small errors bars.Above the   the first lag correction presents higher differences in relation to the other profiles at around 1350 m.For  we use the same procedure for the correction and error analysis that we apply to the  data.The same day was chosen (19 th May), however the period selected is between 12 and 13 UTC, due to the incomplete overlap of MULHACÉN.
In this sense, we studied the influence of noise at two wavelengths: 1064 nm, that has been previously analyzed by Pal et al. (2010) as presented in the section 2 and adopted as reference (considering the rather low impact of molecular signal and the two ways transmittance shown in 9) and 532 nm, just in order to check the feasibility of this wavelength for turbulence studies considering its widespread use in observation networks (Pappalardo et al., 2014;Guerrero-Rascado et al., 2016).Figures 6 and 7 shows the autocovariance function, obtained from ′ 1064 and ′ 532 , respectively, at three distinct heights.As expected,  increases with range, principally above the   .However, the wavelength 532 nm is more influenced by the noise, what can be verified by the higher peak at lag 0 in figure 7, in comparison with peaks at same lag in figure 6.
Although the level of influence of  in each wavelength depends on the  of them (which is associated to technical factors such as laser output power, filters, type of detectors), considering the proposed methodology, to evaluate the composition of each wavelength is also important.The large contribution of to the total  at 532 nm in comparison with the behavior at 1064 nm, can influence the results obtained from such wavelength, because our methodology is based on the use of  ′  .In addition, the larger extinction (due to both aerosol particles and molecules) at 532 nm produces a lower two-way transmittance, resulting in the reduction of the  values at this wavelength.As we used Elastic lidar technique, we could not calculate aerosol extinction profiles, but an estimation of these transmittances was done on the basis of Klett method (Klett, 1985).With this method, a constant lidar ratio value was constrained for each profile using the AOD derived from a collocated AERONET Sun-photometer (Guerrero-Rascado et al., 2008).Using these constrained lidar ratios, the transmittances were calculated together with aerosol backscatter profiles, integrated up to 2.5 km.The estimated two-way transmittance was 0.85 for the case analyzed in this subsection (19 th May).respectively, obtained at 1064 nm, with and without the corrections described in section 3.2.In general, the corrections do not affect the profiles generated from 1064 nm data in a significant way, so that, the higher influence of corrections is observed in the  ′ profile, which is underestimated in some regions.In the figures 9-A, 9-B, 9-C and 9-D we show same high order moments calculated from 532 nm data.As the complexity of moments increases, it is possible to observe the larger influence of the corrections, due to propagation of noise.Nonetheless, the application of the corrections, mainly first lag correction, make these profiles very similar to those generated from the wavelength 1064 nm, so that the same phenomena can be observed in both.
Therefore, in spite of the larger attenuation expected at 532 nm wavelength, which reduces the  of the profiles in comparison with 1064 nm, the application of the proposed corrections, mainly the first lag, reduces significantly such influence and enable the observation of the same phenomena detected in the high-order moments obtained from 1064 nm.Consequently, the wavelength 532 nm will be applied in the analysis presented in section 4.2.The first lag correction was adopted as default because it provides better results than the -2/3 law correction.

Case studies
In this section we present two study cases, in order to show how the products indicated in table 2 can provide a detailed description about the turbulence in the .The first case represents a typical day with a clear sky situation.The second case corresponds to a more complex situation, where there is presence of clouds and Saharan mineral dust layers.

Case study I: clear sky situation
In this case study we use measurements gathered with ,  and pyranometer during 24 hours.The  was operated under operator-supervised mode between 08:20 to 18:00 UTC.
Figure 10 (A) shows the integral time scale obtained from  data ( ′ ).The gray area represents the region where it is not possible to analyze the turbulent process using our  data, either because of the low  values, which results in null values of the  ′ , or because the no null  ′ is smaller than the acquisition time of the .However, the gray area is located almost entirely above the   (white stars).
The  ′ 2 has low values during the entire period when the is present (Figure 10-B).Nevertheless, as air temperature begins to increase (around 07:00 UTC), the  ′ 2 increases together, as well as, the   .
The  ′ 2 reaches its maximum values in the middle of the day, when we also observe the maximum values of air temperature and   .. The combination of  ′ 2 and   provides us a better comprehension about the  growth speed, so that, in the moments where high values of  ′ 2 are observed, it means higher values of Turbulent Kinetic Energy (), which favor the fast ascension of .
The skewness of ′ ( ′ ) is shown in Figure 11-C.The  ′ describes the distribution of the turbulent velocities.Thus positive  ′ implies strong but narrow updrafts surrounded by weaker but more widespread downdrafts, and vice versa for negative  ′ .Consequently, positive values (red regions) correspond with a surface-heating-driven boundary layer, while negative (blue regions) ones are associated to cloud-top longwave radiative cooling.During the stable period, there is predominance of low absolute values of  ′ .
Nevertheless, as air temperature increases (transition from stable to unstable period),  ′ values begin to become larger.Air temperature begins to decrease around 18:00 UTC, and there is a reduction of  ′ , so that, the generation rate of convective turbulence decreases.Therefore, the turbulence cannot be maintained against dissipation, then the  becomes a  covered by the .Thus, the reduction observed in the   is due to the detection of  height.
Figure 10-D shows the values of net surface radiation (  ) that are estimated from solar global irradiance values using the seasonal model described in Alados et al. (2003).The negative values of   are concentrated in the stable region.The   begins to increase around 06:00 UTC and reaches its maximum in the middle of the day.Comparing figures 8-C and 8-D, we can observe similarity among the behavior of  ′ and   , so that, the joint analysis of these variables reinforce the characterization of this  as surfaceheating-driven .Figure 11 shows the  532 profile obtained from 08:00 to 18:00 UTC.At the beginning of the measurement period (08:20 to 10:00 UTC) it is possible to observe the presence of a thin residual layer (around 2000 m a.s.l.), and later from 13:00 to 18:00 UTC it is evident a lofted aerosol layer.In this picture there are the   (pink stars), the   (blue stars), obtained from the maximum of  ′ 2 (Moreira et al., 2018a), and the   (black stars), obtained from the maximum of  ′ 2 (Moreira et al., 2015).In the initial part of measurement, all profiles have similar behavior.However due to distinct  definition and tracer applied by each one, the differences increase as  becomes more complex, e.g. the presence of lofted aerosol layer at 14 UTC.The joint observation of the results provided by these three methods can provide us information about the sublayers in the , both in convective and stable situations.Due to low variability of , the period between 13:00 and 14:00 UTC has been selected to be analyzed from the high order moments.
Figure 12 presents the statistical moments generated from ′ of wavelength 532 nm, which were obtained from 13:00 and 14:00 UTC.The red line in all graphics represent the   (2200 m a.s.l.) and the blue one the average value of   (2250 m a.s.l.), both obtained between 13 and 14 UTC.
Due to presence of a decoupled aerosol layer at 13:30, the average values of   and   shown the high-order moments obtained at the same period described above, however from the 1064 nm data (our reference wavelength).It is possible to observe a similarity between the profiles obtained from each wavelength, so that, the same phenomena observed in the profiles generated from 532 nm and described above, also are detected in the profiles obtained from the reference wavelength.
The results provided by , pyranometer and  data agree with the results observed in figures 12 and 13.In the same way, the analysis of high order moments of ′ fully agree with the information in Figure 10.Thus, the large values of  ′ and ′ detected around 2500 m a.s.l,where we can see a lofted aerosol layer, suggest the ascent of an aerosol layer and presence of a peaked distribution, respectively.

Case study: dusty and cloudy scenario
In this case study measurements with ,  and pyranometer expand during 24 hours, while  data are collected from 09:00 to 16:00 UTC.
Figure 14-A shows  ′ .Outside the period 13:00 to 17:00 UTC, the greatest part of grey area is situated above the   (white stars), thus  time acquisition is enough to perform studies about turbulence in this case.
′ 2 has values close to zero during all the stable period (Figure 14-B).However, when air temperature begins to increase (around 06:00 UTC), the  ′ 2 also increases and reaches its maximum in the middle of the day.The higher values of  growth speed are observed in the moments where  ′ 2 reaches its maximum values.In the late afternoon, as air temperature decrease, the values of  ′ 2 (and consequently the ) decrease gradually, until reach the minimum value associated to the .As mentioned before, Figure 15 shows the  profile obtained from 09:00 to 16:00 UTC in a complex situation, with presence of decoupled dust layer (around 3800 m a.s.l.) from 09:00 to 12:00 UTC and the presence of both middle altitude clouds and very intense dust layers (around 3500 m a.s.l.) from 11:30 to 16:00 UTC.The pink, black and blue stars represent the   ,   and   respectively.Due to the presence of dusty layers and clouds, the difference between the methods is more evident, mainly of the   , which uses the aerosol as tracers.This method only produces results close to the others at 15 UTC, when dust layer is mixed with the .
Figure 16 illustrates the statistical moments of ′ of 532 nm wavelength obtained from 11:00 to 12:00 UTC.The  ′ 2 profile presents several peaks due to the presence of distinct aerosol sublayers.The first peak is coincident with the value of   .The value of   , is coincident with the base of the dust layer.This difficulty to detect the  in presence of several aerosol layers is inherent to the variance method (Kovalev and Eichinger, 2004).However, the joint observation of   and In order to show the feasibility of 532 nm wavelength, in the figure 17 are presented the high-order moments obtained between 11-12 UTC from 1064 nm wavelength data.Although the error of  ′ 2 obtained from 532 nm (pink shadow) is considerably higher than the error of same variable obtained from 1064 nm, all profiles are very similar, so that, the same phenomena can be observed in both graphics (figure 16 and 17).
Figure 18 shows the ′ 532 nm wavelength high-order moments obtained from 12:00 and 13:00 in presence of cloud cover.The method based on maximum of  ′ 2 locates the   at the cloud base, due to the high variance of ′ generated by the clouds. ′ presents values larger than  time acquisition, therefore this configuration enable us to study turbulence by  analyses. ′ has few peaks, due to the mixing between  and dust layer, generating a more homogenous layer.The highest values of  ′ are observed in regions where there are clouds, and the negative ones (between 3500 and 4000 m) occur due to presence of air from  between the two aerosol layers (Figure 15).The inflection point of  ′ profile is observed in   region. ′ profile has low values in most of the , demonstrating the high level of mixing during this period, where dust layer and  are combined.The higher values of  ′ are observed in the region of clouds.In the same way of the previous analysis, the high-order moments of the period mentioned above were calculated for the wavelength of 1064 nm (figure 19).Although there are some differences in the absolute values of some profiles, the high-order moments generated using 1064 and 532 nm have similar profiles, so that, the same phenomena can be observed, demonstrating the viability of 532 nm wavelength in the proposed methodology.

Conclusions
In this paper we perform an analysis about the  turbulent features from three different types of remote sensing systems (, and ) and surface sensors during SLOPE-I campaign.We applied two kind of corrections to the lidar data: first lag and -2/3 corrections.The corrected  statistical moments showed little variation with respect to the uncorrected profiles, denoting a rather low influence of the noise in high  conditions.The  high-order moments were obtained from two wavelengths: 1064 nm, adopted as reference, and 532 nm, in order to verify the viability to use the last one in turbulence analysis.From this comparison, was possible to observe that the wavelength 532 nm is more affected by noise, in comparison with 1064 nm, due to the large contribution of the molecular component and the reduced two-way transmittance at that wavelength.However, the application of proposed corrections, mainly the first lag, can reduce such influence, so that, the same phenomena can be observed in the high-order moments provided from both wavelengths The case studies present two kind of situations: well-defined PBL and a more complex situation with the presence of Saharan dust layer and some clouds.In both cases was possible to identify the events describe in table 2. The combined use of remote sensing systems shows how the results provided by the different instruments can complement one each other, providing a detailed observation of some phenomena, mainly in complex situations.
Therefore, this study shows the feasibility of the described methodology based on the combination of remote sensing systems for retrieving a detailed picture on the  turbulent features.In addition, the feasibility of using the analyses of high order moments of the  collected at 532 nm at a temporal resolution of 2 s offers the possibility for using the proposed methodology in networks such as EARLINET or LALINET with a reasonable additional effort.(Lenschow et al., 2000) Table 2 -Products and their respective meaning, provided by each system

2
decrease slowly due to location of the lofted aerosol around 2500 m." Line 372: Define FT here.Done Line 392: Do you refer to the correct figure here?

Figures 5 -
Figures 5-B and 5-C show the comparison of variance ( ′ 2 ) and skewness ( ′ ), respectively, with and without corrections.The profiles corrected by -2/3 law do not present significant differences in comparison to uncorrected profiles.On the other hand, the profiles corrected by the first lag correction have slight differences below the   , mainly the  ′ 2 ( ′ only in the first 50 m).Therefore, considering high Signal-to-Noise Ratio () conditions, although the presence of  can change slightly the value of high order moments, it is not enough to distort the observed phenomena as shown by the impact of the corrections applied.

Figure 10 -
Figure10-E presents the values of surface air temperature and surface relative humidity ().Air surface temperature has a daily pattern similar to that of   and  ′ .On the other hand,  is inversely correlated with the temperature.
have a difference of around 500 m.The  ′ 2 has small and practically constant values between 1000 and 1400m, evidencing the homogeneity of the aerosol distribution in this region.Starting at 1400 m the value of  ′ 2 begins to increase, reaching a positive peak at   , which represents the Entrainment Zone (region characterized by an intense mixing between air parcels coming from  and Free Troposphere (), causing a high variation in aerosol concentration).The   observed at approximately 2900 m demonstrate an inherent difficulty of the variance method to detect the  in the presence of several aerosol layers(Kovalev and Eichinger, 2004).Above   the values of  ′ 2 decrease slowly due to location of the lofted aerosol around 2500 m.However, above this aerosol layer the value of σ RCS′ 2 is reduced to zero, indicating a large homogeneity in aerosol distribution at this region, what is expected, because the aerosol concentration at the  is negligible in this case.The integral time scale obtained from ′ ( ′ ) has values higher than  time acquisition throughout the , evidencing the feasibility for studying turbulence using this elastic lidar configuration.The skewness values obtained from ′ ( ′ ) give us information about aerosol motion.The positive values of  ′ observed in the lowest part of profile and above the   represents the updrafts aerosol layers.The negative values of  ′ indicates the region with low aerosol concentration due to clean air coming from .This movement of ascension of aerosol layers and descent of clean air with zero value of  ′ at  (characteristic of the  growing)was also detected byPal et al. (2010) andMcNicholas et al. (2014).The kurtosis of ′ ( ′ ) determines the level of mixing at different heights.There are values of  ′ larger than 3 in the lowest part of profile and around 2500 m, showing a peaked distribution in this region.On other hand, values of  ′ lower than 3 are observed close to the   , therefore this region has a well-mixed  regime.Pal et al. (2010) andMcNicholas et al. (2014) also detected this feature in the region nearby the .In figure13 are Figure14-Cshows the profiles of  ′ .The main features of this case are: the low values of  ′ , the slow increase and ascension of positive  ′ values and the predominance of negative  ′ values from 12:00 to 13:00 UTC.The first two features are likely due to the presence of the intense Saharan dust layer (Figure15), which reduces the transmission of solar irradiance, and consequently the absorption of solar irradiance at the surface, generating weak convective process.From figure16we can observe the presence of both middle altitude clouds and very intense dust layers from 12:00 to 15:00 UTC.Such combination contributes to the intense negative values of  ′ observed in this period until around 2 km, because, as mentioned previously,  ′ is directly associated with the direction of turbulent movements.The present situation can be considered representative of cloud-top long-wave radiative cooling in the (Ansmann et al., 2010).The influence of Saharan dust layer can also be evidenced on the   pattern (Figure14-D), which maintains negative values until 12:00 UTC and reaches a low maximum value (around 200 W/m²).The observation of  ′ and   between 12:00 and 14:00, as well as, the presence of clouds and geometrically thick dust layers during this same period, reinforces the hypothesis that we have a situation of the cloud-top longwave radiative cooling in the .Air surface temperature and  (Figure14-E) present the same correlation and anti-correlation (respectively) observed in the earlier case study, where the maximum of air surface temperature and the minimum of  are detected in coincidence with the maximum daily value of   .
, enable us to characterize and distinguish the several sublayers.The values of  ′ are higher than  acquisition time all along the , evidencing the feasibility of  time acquisition for studying the turbulence of  in this case.The  ′ profile has several positive values, due to the large number of aerosol sublayers that are present.The characteristic inflection point of  ′ is observed in coincidence with the   , that confirming the agreement between this point and the .From the analysis of  ′ and  ′ is possible to justify this phenomena from the mixing process demonstrated in the earlier case study.The  ′ has predominantly values lower than 3 below 2500 m, thus shown how this region is well mixed as can see in Figure16.Values of  ′ larger than 3 are observed in the highest part of profile, where the dust layer is located.

Figure 2 -Figure 3 -
Figure 2 -Flowchart of data analysis methodology applied to the study of turbulence with Doppler lidar

Figure 6 -
Figure 6 -Autocovariance of ′  obtained from MULHACÉN elastic lidar data to three different heights on 19th May 2016 at 12-13 UTC in Granada.

Figure 7 -
Figure 7 -Autocovariance of ′  obtained from MULHACÉN elastic lidar data to three different heights on 19th May 2016 at 12-13 UTC in Granada.

Figure 10
Figure 10 -Aintegral time scale obtained from Doppler lidar data [ ′ ], Bvariance obtained from Doppler lidar data [ ′ 2 ], Cskewness obtained from Doppler lidar data [ ′ ], Dnet radiation obtained from pyranometer data [  ], E -Air surface temperature [blue line] and surface relative humidity [ -orange line] both were obtained from surface sensors.All profiles were acquired on 19th May 2016 in Granada.In A, B and C black lines and white stars represent air temperature and   , respectively.

Figure 11 -
Figure 11 -Time-Height plot of RCS obtained on 19 May 2016 in Granada.Pink stars represent the   , black stars represent the   and blues stars represent the   .

Figure 14
Figure 14 -Aintegral time scale from Doppler lidar data [ ′ ], Bvariance from Doppler lidar data [ ′ 2 ], Cskewness from Doppler lidar data [ ′ ], Dnet radiation from pyranometer data [  ], E -Air surface temperature [blue line] and surface relative humidity [orange line] from surface sensor data.All profiles were obtained in Granada on 08 July 2016.In A, B and C black lines and white stars represent air temperature and   , respectively.

Figure 15 -
Figure 15 -Time-Height plot of RCS obtained from MULHACÉN elastic lidar data on 08 July 2016 in Granada.Pink stars represent the   , black stars represent the   and blues stars represent the   .
Yes, exactly.In order to clarify this point the text has been changed as follow: "…gathered data [(, )] with a temporal resolution …"Line 179: Do you mean the PBLH_Doppler is attributed to the height where σw2 drops below a pre-determined threshold?

Table 1 -
Variables applied to statistical analysis