Analyzing the Atmospheric Boundary Layer by high-order moments obtained from multiwavelength lidar data: impact of wavelength choice

The lowest region of the troposphere is a turbulent layer denominated Atmospheric Boundary Layer (ABL) characterized by high daily variability due to the influence of surface forcings. This is the reason why detecting systems with high spatial and temporal resolution, like lidars, have been widely applied for researching this region. In this paper, we present a comparative analysis on the use of lidar backscattered signals at three wavelengths (355, 532 and 1064 nm) to study the ABL investigating the high-order moments, which give us information about the ABL height (derived by the variance method), 5 aerosol layers movement (skewness) and mixing conditions (kurtosis) at several heights. Previous studies have shown that 1064-nm wavelength, due to the predominance of particle signature in the total backscattered atmospheric signal and practically null presence of molecular signal (which can represent noise in high-order moments), provides an appropriate description of the turbulence field and thus, in this study, it was considered as a reference. We analyze two case studies, which show us that the backscattered signal at 355 nm, even after applying some corrections, has a limited applicability for turbulence studies 10 using the proposed methodology due to the strong contribution of the molecular signature to the total backscatter signal. This increases the noise associated to the high-order profiles and, consequently, generates misinformation. On the other hand, the information on the turbulence field derived from the backscattered signal at 532 nm is similar to that obtained at 1064 nm due to the appropriate attenuation of the noise, generated by molecular component of backscattered signal, by the application of the corrections proposed. 15 Copyright statement. TEXT


Introduction
The Atmospheric Boundary Layer is the part of the troposphere that is directly or indirectly influenced by the Earth's surface (land and sea), and responds to gases and aerosol particles emitted at the Earth's surface and to surface forcing at time scales less than a day. Forcing mechanisms include heat transfer, fluxes of momentum, frictional drag and terrain-induced flow modification. The height of this layer (ABLH) varies from hundreds of meters until some kilometers due to the intensification 5 or reduction of convective or mechanical processes with additional contribution from orographic effects. The ABL presents a daily pattern controlled by the energy balance at the Earth's surface. Thus, after sunrise the positive net radiative flux (R n ) induces the raise of surface air temperature that initiates the convective process, which is responsible for the growth of the so-called Mixing Layer (ML) or Convective Boundary Layer (CBL). This layer grows along the day extending the region affected by the convective process until around midday, when it reaches their maximum development. Slightly before sunset, 10 the decrease of the incoming solar irradiance at the surface results in a radiative cooling of the Earth's surface. This cooling affects the closest air layer, diminishing the convective process. In this way, the CBL disappears and two new layers characterize the ABL, a stable and stratified layer denominated Stable Boundary Layer (SBL) at the bottom and the Residual Layer (RL) over the last one with characteristics of the previous day's ML (Stull, 1988).
The turbulent features of the ABL are relevant in air quality and weather forecasting and thus are worthy of study. As a 15 rule, the turbulent processes are treated as nondeterministic and, therefore, the turbulence is characterized by its statistical properties. Thus, high order statistical moments are used to generate information about the turbulent fluctuation field, besides a description about mixing processes in the ABL . ABL turbulence has been commonly studied by means of anemometer towers (e.g., Kaimal and Gaynor, 1983) and aircrafts (e.g., Lenschow et al., 1980;Williams and Hacker, 1992;Lenschow et al., 1994;Stull et al., 1997;Andrews et al., 2004;20 Vogelmann et al., 2012). Nevertheless, the first ones have a use restricted to regions near the surface, due to their limited vertical range. Aircrafts offer an alternative approach that allows extending the analyses to higher atmospheric layers, but conversely, they have a reduced time window, thus limiting the period of analysis. Due to the large variability of the ABL characteristics along the day, the use of systems endowed with high spatial and temporal resolution allow studies with a higher degree of details. Consequently, remote sensing systems (mainly lidars) become an important tool in ABLH detection 25 (Martucci et al., 2007;Pal et al., 2010;Wang et al., 2012), as well as, in turbulence studies (Lagouarde et al., 2013(Lagouarde et al., , 2015. In addition, the different lidar techniques offer the possibility of analyses with several variables, such as vertical wind velocity by Doppler lidar (Lenschow et al., 2000;Lothon et al., 2006;O'Connor et al., 2010), water vapor mixing profile by Raman lidar or Differential Absorption lidar (DIAL) (Wulfmeyer, 1999;Kiemle et al., 2007;Wulfmeyer et al., 2010;Turner et al., 2014;Muppa et al., 2016), temperature by rotational Raman lidar (Hammann et al., 2015) and aerosol number density by elastic lidar 30 or High Spectral Resolution lidar (HSRL) McNicholas and Turner, 2014). Therefore, a wider range of results can be obtained, especially when different types of systems are synergistically used, as shown by Engelmann et al. (2008) who combine elastic and Doppler lidar data for deriving the vertical aerosol flux. Pal et al. (2010) have shown that it is feasible the use of elastic lidar measuring at a high acquisition rate for characterizing the atmospheric turbulence. In particular, they have shown that the fluctuation of the Range Corrected Signal (RCS) at 1064 nm is a proxy for the fluctuation of the particle concentration, due to predominance of particle signature (β par ) in the total backscattered signal at this wavelength, and, thus, it can be used for observing the turbulent aerosol movements in the CBL.
However, if other wavelengths are used in this kind of analysis, the effects of molecular backscatter coefficient (β mol ) and 5 atmospheric extinction (α) must be considered. In this work, we perform a comparative analysis regarding the use of three different wavelengths, namely 355, 532 and 1064 nm (the last one adopted as reference), to obtain the high-order moments, i.e. variance (σ 2 ), skewness (S) and kurtosis (K), and also the integral time-scale (τ ). Moreover, it was analyzed the interference of noise ε and β mol over the high-order moments and τ obtained from each one of the considered wavelengths, in order to quantify how such factors can influence the correct interpretation of the statistical variables. The goal of this study is to 10 show the viability of the proposed methodology for studying the turbulence by computing the high-order moments of the backscattered signal at different wavelengths. We pay special attention to the advantages and limitations of each wavelength analyzed considering the importance of the proposed correction schemes. This paper is organized as follows. The measurement site and the experimental set up are introduced in Section 2. The methodology is described in Section 3. The comparisons and case studies are analyzed in Section 4. Conclusions are given in Section 5. America, with a population of approximately 12 million of inhabitants, and endowed a subtropical climate where winter is mild (15 • C) and dry, while summer is wet and has moderately high temperatures (23 • C) (IBGE, 2017). The São Paulo Lidar station (SPU) has a coaxial ground-based multiwavelenght Raman lidar system operated at LEAL. The system operates with a pulsed Nd:YAG laser, emitting radiation at 355, 532 and 1064 nm, a laser repetition rate of 10 Hz and a laser beam pointing to zenith direction. The pulse energy (and stability) of each wavelength are 225 mJ (2mJ) at 355 nm, 400 mJ (4 mJ) at 532 nm, 25 and 850 mJ (6 mJ) at 1064 nm. The MSPI lidar detects three elastic channels at 355, 532 and 1064 nm and three Raman-shifted channels at 387 nm, 408 nm (corresponding to the shifting from 355 nm by N 2 and H 2 O) and 530 nm (corresponding to the Rotational Raman shifting from 532 nm by N 2 (Veselovskii et al., 2015)). This system is equipped with photomultipliers Hamamatsu R7400. The SPU lidar reaches full overlap at around 300 m a.g.l. (Lopes et al., 2018). This system operates with a temporal and spatial resolutions of 2 s and 7.5 m, respectively.

Methodology
The turbulence study is based on the observation of the fluctuation q (t) of a determined variable (q) in the time t. The values are obtained as follows: firstly q(t) are averaged in packages that cover a certain time interval, from which the mean value (q) is extracted. Then, such value is subtracted from each q(t) value, providing the fluctuation q (t) as demonstrated in the equation below by Reynold's decomposition (de Arruda Moreira et al., 2019): In the analysis performed with elastic lidar systems, the variable of interest is the aerosol number density (N ), from which we obtain its fluctuation (N ) by the equation 1. However, elastic lidar systems do not provide directly the value of N . Therefore, considering the validity of Mie-theory (where the aerosol backscatter coefficient is linked to the backscatter efficiency, particle radius (r) and the number of particles with radius r we can write the 2, under several assumptions. The premises adopted 10 here are (i) the variation of aerosol size with the relative humidity can be neglected, (ii) the atmospheric volume probed is composed by similar types of aerosol particles and (iii) the fluctuations of the aerosol microphysical properties are smaller than the fluctuations of the total number density in the volume probed by the lidar. More details about these assumptions can be found in Pal et al. (2010). Feingold and Morley (2003) and Titos et al. (2016) demonstrated the relation between relative humidity and hygroscopic growth, so that, such effects can start at 80% RH. The two cases presented in this work were gathered 15 in winter, the driest season of São Paulo. In particular, RH was below 80% in both days (see section 4). Such value is lower than the RH threshold to hygroscopic effects indicated by the two papers above mentioned. Consequently, ignoring the hygroscopic growth and assuming similar types of aerosol throughout the atmospheric column, the following equation can be used: 20 where β aer and β aer represent the particle backscatter coefficient and its fluctuation, respectively. The variable z is the height above the ground, t is the time and Y is a variable that does not depend on time.
The lidar equation is defined as follows Weitkamp (2005): , and A is the effective area of the telescope receptor [m 2 ], η is a variable related to the efficiency of the lidar system and O(λ, z) is the laser-beam receiver-field-of-view overlap function. The most important quantities are β(λ, z), which is the total backscatter coefficient, due to atmospheric molecules, β mol (λ, z), and aerosol β aer (λ, z), in other words, β(λ, z) = β mol (λ, z) + β aer (λ, z) [(m.sr) −1 ] at distance z, and α(λ, z) is the total extinction coefficient, due to atmospheric molecules, α mol (λ, z), and aerosols α aer (λ, z), in other words, α(λ, z) = α mol (λ, z)+α aer (λ, z) [(m) −1 ] at distance z. If the wavelength 5 1064 nm is used, we can neglect the influence of the extinction coefficient α(λ, z) provided by aerosol, the Rayleigh scattering generated by atmospheric molecules) and the β mol (λ, z) . Therefore, the equation 4, for the wavelength of 1064 nm, can be rewritten as follows: where RCS 1064 is the Range Corrected Signal, G is a constant and the subscribed indexes represent the wavelength and the 10 particles. Then, applying Reynold's decomposition (Eq. 1) over Eq. 5, the following equation is derived: Our purpose is to evaluate the use of other wavelengths when the effects of molecular backscatter coefficient (β mol ). The interest is based on the best performance of the technology for detecting wavelengths in the VIS and UV and on the extended

High-order moments
The high-order moments used in this study are obtained from RCS (z, t), generated by equation 1, where RCS(z) represents the 1-hour average package of RCS(z, t) data. From this, the high order moments, variance (σ 2 ), skewness (S) and kurtosis 20 (K) are obtained as demonstrated in the first column of Table A1 , as well as, their corrections and errors in the second and third columns of the same table, respectively. In table A2 are presented the physical meaning of each high-order moment in the context of the proposed analysis The integral time scale (τ ) is an important prerequisite in turbulence studies. It guarantees that the most part of the horizontal variability of the turbulent eddies is detected with good resolution, enabling the solution of inertial subrange and dissipation 25 range in the spectrum and autocorrelation function, respectively . τ must be larger than the temporal resolution of the analyzed time series (SPU Lidar station time acquisition is 2s). In the same way of high-order moments, such variable is obtained from RCS (z, t) as shown in the first column of Table A1.

Error analysis
The high-order moments and τ generated from RCS (z, t) can also be obtained from the following autocovariance function M ij , which has its order represented by the sum of the subscript i and j , according to the following equation: where tf means final time. However, it is important to consider the influence of instrument noise ε(z, t) in the RCS (z, t) 5 profile. Therefore, M ij can be rewritten as follows: Although atmospheric fluctuations are correlated in time, ε(z, t) is random and uncorrelated with the atmospheric signal, therefore ε(z, t) is only associated with lag 0. Consequently, it is possible to obtain the corrected autocovariance function, M 11 (→ 0), removing the error ∆M 11 (0) of the uncorrected autocovariance function M 11 (0), as demonstrated in the equation 10 below: Based on this concept, Lenschow et al. (2000) proposed two methods to correct for the noise influence: -First lag correction: the lag 0 (∆M 11 (0)) is directly subtracted from the uncorrected autocovariance function M 11 (0), generating M 11 (→ 0).

15
--2/3 law correction: A new lag 0 value is obtained by the extrapolation of M 11 (0) to the firsts non-zero lags back to lag zero, using the inertial subrange hypothesis (Monin and Yaglom, 2013): where C represents a parameter of turbulent eddy dissipation rate. In this study, we also used the first five points after lag 0 to perform this correction. In Table A1 the second and third columns present the corrections and errors, respectively, of high-order 20 moments and τ .

Results
In this section we present two case studies, applying the methodology described in section 3, in order to perform a comparative analysis about the influence of β mol , and ε in the high-order moments and τ obtained from different wavelengths (355, 532 and 1064 nm). In this case study we gathered measurements from 13:00 to 19:00 UTC. Figure A2 shows the time-height plot of RCS 532 during this period. This case is composed by two distinct periods, in the first two hours there is a RL with an underlying shallow CBL. Nevertheless, in the last part of the second hour the CBL quickly grows and it mixes with RL forming a fullydeveloped ABL, with its top situated between 1500 and 1600 m from 15:00 to 19:00 UTC. The black dotted box, between 15 17:00 and 18:00 UTC represents the period selected to perform the statistical analysis.
In order to check the hypothesis proposed by Pal et al. (2010), which assumes that there is not particle hygroscopic growth and that the same type of aerosol is present in the entire atmospheric column in the ABL region, were analyzed the relative humidity and mixing ratio profile retrieved from radio-sounding measurements (http://weather.uwyo.edu/upperair/sounding.html), launched at the Campo de Marte Airport (São Paulo, Brazil), which is about 10 km away from the SPU lidar system. Figure A3- 20 A and A3-B shows the relative humidity and mixing ratio profiles, respectively, measured on 26 th July 2017 at 12 UTC. Both, relative humidity and mixing ratio can be considered constant below 1500 m, with mean values of 67 ± 8% and 7.6 ± 0.9 g/kg, respectively. Since there are no large variation of water vapor mixing ratio and relative humidity values in this region, we assume that this case is not affected by particle hygroscopic growth. In addition, the AERONET Sunphotometer (Holben et al., 1998a) data from the São Paulo station were retrieved in order to check the aerosol type, as can be seen in the figure A3-C. 25 According to Eck et al. (1999), the Ångström Exponent (AE) can be a useful tool to distinguish different types of atmospheric aerosols. Figure A3-C shows the aerosol AE time series for the case study of 26 th July 2017. The AE was calculated at the spectral range 340-440 nm and 440-675 nm using AERONET (Holben et al., 1998b) products from Level 1.5 version 3 data.
For this measurement period the percentage variation of AE was no more than 3% in both cases. Therefore, there are no considerable changes during the whole measurement period, which is a strong indication that there is no aerosol type change 30 throughout the day.
In figure A4 is presented the Signal-to-Noise-Ratio (SNR) profile of the raw lidar signal, as calculated by Heese (2010), of the three wavelengths (1064 nm (red line), 532 nm (green line), and 355 nm (violet line)) during the analyzed period.
All wavelengths have values of SNR higher than 1 (the threshold for good quality) below the ABLH (dotted blue line) with predominance of values lower than 1 in the Free Troposphere (FT), what was expected due to the strong reduction of aerosol concentration in such region. Although the three wavelengths have similar SNR profiles, close to ABLH the difference among 5 then become more evident, principally the fast decreasing of the 355 nm and the high values of 532 nm. Figure A5 shows the autocovariance function (ACF), obtained between 17:00 and 18:00 UTC for the wavelengths 355 (ACF 355 ), 532 (ACF 532 ) and 1064 nm (ACF 1064 ) at 1000 m agl and 1700 m agl. Thus, from the comparison of the figures A2 and A5 it is possible to observe that the altitude chosen at 1000 m (red line) is situated below the top of CBL, while the altitude chosen at 1700 m (light green line) is in the FT. As expected, the ε, which is represented by the peak on the lag 0 of 10 the autocovariance function (A5), increases with height for all the wavelengths due to reduction of aerosol load with height.
ACF 355 has the lowest intensity (around 90% smaller those of ACF 532 and ACF 1064 ) and it is clearly much more affected by the magnitude of ε that represents approximately 25% of ACF 355 , while for ACF 532 and ACF 1064 the noise represents around 10% of the respective autocovariance. Figure A6 presents all statistic variables, their respective corrections and errors (shadows), generated from the methodology 15 described in section 3, for data acquired between 17:00 and 18:00 UTC.
The variance profiles, σ 2 RCS (z), with and without corrections for all wavelengths are represented in Figure A6, from 1 to 9. The low and almost constant values of uncorrected σ 2 RCS1064 (z) from the bottom until around 1000 m of altitude demonstrates an almost constant distribution of aerosol particles in this region, as can be seen in Figure A6.1. Above 1000 m of altitude, the value of uncorrected σ 2 RCS1064 (z) increases, reaching its maximum peak at around 1600 m. This peak represents the 20 Entrainment Zone, the region where a mixing occurs between air parcels coming from the CBL and FT. According to Menut et al. (1999), there is an intense variation of aerosol concentration during this process, generating a maximum in the uncorrected σ 2 RCS1064 (z), which represents the ABLH. Above the ABHL, the aerosol concentration is considerably lower than in CBL and, thus, the uncorrected σ 2 RCS1064 (z) is reduced to practically zero. This methodology to estimate the ABLH is named Variance Method or Centroid Method and it was described by Hooper and Eloranta (1986) and Menut et al. (1999), respectively. The 25 main limitations of this method are its applicability only for CBL, and the ambiguous results in complex cases, like as the presence of several aerosol layers (Emeis, 2011). In such situations more sophisticated methods like as Wavelet , PathfinderTURB (Poltera et al., 2017) and POLARIS (Bravo-Aranda et al., 2017) are recommended.
The uncorrected σ 2 RCS532 (z), presented in Figure A6.4 is rather similar to uncorrected σ 2 RCS1064 (z), including the position of maximum peak. Nevertheless, although uncorrected σ 2 RCS355 (z), presented in Figure A6.7, also has the maximum peak 30 situated at around 1600 m of altitude, the profile is nosier than the profiles obtained from the other wavelengths and, therefore, it is not possible to identify the regions with uniform aerosol distribution as evidenced in uncorrected σ 2 RCS1064 (z). Although the σ 2 RCS355 (z) is nosier than another ones, there is a low difference among the ABLH estimated from the three different wavelengths (lower than 10%).
The correction 2/3, shown in Figures A6.2, A6.5 and A6.8, does not cause significant changes in the uncorrected profiles. On the other hand, the first lag correction changes significantly the profiles, thus σ 2 RCS532 (z) becomes very similar to σ 2 RCS1064 (z), while σ 2 RCS355 (z) continues with some differences, mainly in the region below the ABLH, as can be seen in Figures A6.3, A6.6 and A6.9.
The integral time scale profiles τ RCS (z), with and without corrections τ corr RCS (z) and τ unc RCS (z), respectively, calculated 5 for the three wavelengths are presented in the Figure A6, from 10 to 18. The τ unc RCS (z) presents values larger than SPU Lidar station time acquisition showed as black dotted line, in the region below ABLH at all wavelengths, as can be seen in Figures   A6.10, A6.13 and A6.16. The largest values of τ unc RCS (z) correspond to 1064 nm, while the lowest values are computed for 355, which is practically half of those obtained with the reference wavelength, 1064 nm. The low value for the τ unc RCS (z) at 355 nm can be associated to the influence of the noise in the signal retrieved at this wavelength. The application of the correction 10 2/3 does not cause significant changes in the profiles, while the first lag correction changes significantly the profiles mainly in the region below the ABLH, as can be checked in Figures A6.11, A6.14 and A6.17, and in Figures A6.12 The kurtosis profile K RCS is the most complex high-order moment presented in this study and, consequently, in such profiles the differences among the three wavelengths are more evident. In the context of our analysis, the values of K RCS are 35 indicators of the mixing degree at each altitude, as well as, of the intermittence of turbulence caused by large eddies. In reason of some technical limitations of our lidar system, it is possible to resolve eddies only until a predetermined size. Therefore, in regions where turbulence is performed in too small scales, our system cannot solve these eddies. The kurtosis equation presented in the table A1 represents the kurtosis of a Normal Distribution, which is equal 3 (Bulmer, 1965), consequently such value is applied as threshold in the analyses performed in this paper. Values lower than 3 represents a well-mixed region, 5 indicating a flatter distribution in comparison with a normal distribution, thus the turbulence caused by large eddies can be characterized as frequent. In contrast, values higher than 3 indicates a peaked distribution in comparison with a Gaussian distribution, in other words, there is an unusual variation in the RCS (z, t), which represents a low degree of mixing, and the presence of an infrequent large eddies turbulence .
The K unc RCS at 532 and 1064 nm have some differences in the region below 1300 m of altitude, where the profile at 1064 nm 10 only shows values higher than 3, representing a region with low degree of mixing, while the K unc RCS obtained from 532 nm is composed by values higher and lower than 3. From 1300 m to 3500 m of altitude, the profiles of these two wavelengths are very similar, with values lower than 3 in the region below the ABLH, characterizing a well-mixed region, a peak of values higher

Case Study II: 19 th July 2018
In this case study measurements were gathered with the SPU Lidar station from 12:00 to 21:00 UTC. Figure A8 shows the time-height plot of RCS 532 during this period. In the beginning of measurement it is possible to observe the presence of an ascending CBL covered by a RL, which has the top situated at around 1300 m of altitude. At approximately 15:30 UTC the CBL breaks up the RL and becomes fully-developed, thus, its growth speed is reduced and the value of top height maintains practically constant (1600 m) from 17:00 UTC until 21:00 UTC. The black dotted box in Figure A8 represents the chosen period to perform the statistical analysis (18:00 -19:00 UTC).
In the same way of Case Study I, the hypothesis proposed by Pal et al. (2010) is validated from the profiles presented in 5 Figure A9. The profiles of relative humidity and mixing ratio, presented in the Figure A9-A and A9-b, respectively, do not have large variations in the CBL below 1200 m of altitude. In addition, the aerosol optical depth related Ångström Exponent time series did not show considerable changes during the whole measurement period, as can be seen in Figure A9-C. For this measurement period the percentage variation of AE was no more than 4% and 3% in the spectral range 340-440 nm and 440-675 nm, respectively. Therefore, there are no considerable changes during the whole measurement period, which is a 10 strong indication that there are no aerosol type change throughout the day and the atmospheric conditions are not propitious for particle hygroscopic growth events. FT region all profiles have values lower than 1, as expected. Figure A11 shows a comparison among the ACF obtained from the three wavelengths 1064 nm (left), 532 nm (center) and 355 nm (right), between 18:00 and 19:00 UTC, at two heights 1000 m (red line) and 1700 (green line). In the same way of case Study I, the region above ABLH (green line) is more influenced by noise than the region situated below this height (red line).
The intensity of ACF 532 and ACF 1064 are very similar, although the presence of noise in the first one, which is 40% and 46%, 20 below and above ABLH, respectively, is higher than in the last one, 27% and 30%, below and above ABLH, respectively. The ACF 355 presents a lower intensity value in comparison with the other two wavelengths, and a strong presence of noise below and above the ABLH, 50% and 67%, respectively.
The three high order moments and τ RCS , both corrected by the first lag correction and obtained between 18:00 and 19:00 UTC, are presented in figure A12. The τ corr RCS for all wavelengths has values higher than 2s from the bottom of profile until the 25 first meters above the ABLH elastic with maximum of σ 2 RCS (z). Although the values obtained from 1064 nm and 532 nm are almost twice as large as the values generated from 355 nm. In the same way of Case I, although there are some differences among the maximum of the [σ 2 RCS (z)], they do not influence significantly the ABLH estimation, so that, the difference among the ABLH obtained from each wavelength is lower than 10%. The positive values of S corr RCS (z) of 1064 nm indicate the presence of aerosol updrafts from the bottom of the profile until around 750 m of altitude. From this height until the ABLH, the S corr RCS (z) 30 is characterized by negative values, which represents a region with entrainment of clean FT air into the CBL. In the same way of case study I, there is an inflection point at ABLH, which reproduces the transition from negative to positive values, the last ones indicating the presence of aerosol updraft layers in the first 200 m above the ABLH. Such behavior in the region of ABLH also was observed by Pal et al. (2010) and McNicholas and Turner (2014) and it can be considered characteristic of convective regime. The S corr RCS (z) obtained from the wavelengths 1064 and 532 nm presents an identical pattern of behavior, demonstrating 35 the occurrence of the same phenomenon. The S corr RCS (z) obtained from the wavelength 355 nm, in the same way of the previous case study, does not exhibit the behavior observed in the reference wavelength, presenting only positive values in the whole profile. Therefore, it is not possible to identify variations in the aerosol dynamic using 355 nm.
The K corr RCS (z) obtained from the wavelength 1064 nm presents values higher than 3 from the bottom until around 1300 m of altitude, characterizing a region with low degree of mixing. From 1300 m until the ABLH the K corr RCS (z) has values lower 5 than 3, that characterize this region as showing a large degree of mixing and more evidently the presence of turbulence. Such behavior occurs mainly due to of entrainment of cleaner air. A few meters above the ABLH, the K corr RCS (z) has a great peak, which occurs due to rare aerosol plumes penetrating at this region. Such behavior also was observed in case study I, as well as, by Pal et al. (2010) and McNicholas and Turner (2014). Above the ABLH the profile has values only higher than 3, however, as τ corr RCS (z) decreases to values close to zero and low values of SNR of the RCS'are characteristic of this region, it is not possible to extract conclusive information from K corr RCS (z). In the same way of the comparison performed with other variables, the K corr RCS (z) obtained from the wavelength 532 nm presents similar behavior to profile obtained from 1064 nm, thus, the same phenomenon can be observed. On the other hand, the K corr RCS (z) obtained from the wavelength 355 nm does not allow observing the behavior detected in the profile obtained from the reference wavelength, because along the whole profile the K corr RCS (z) at 355 nm presents values higher than 3.
15 Figure A13 shows the composition signal of β aer and β mol , retrieved during the analyzed period of this case study (18:00 -19:00 UTC) using the Klett-Fernald-Sasano inversion (Klett, 1983, 1985Fernald, 1984;Sasano and Nakane, 1984), at each one of the three wavelengths, as well as the β ratio calculated using the backscatter profile of aerosol and molecular component (Bucholtz, 1995). From figure A13-1 it is possible to observe that the backscattered signal at 1064 nm has a predominance of β aer , with almost null values of β mol . The composition of the backscattered signal at 532 nm is shown in figure A13-3. 20 Although, the component β mol has values higher than that ones observed in wavelength 1064 nm, the component β aer is predominant in the backscattered signal composition. The backscattered signal at 355 nm, presented in figure A13-5, unlike the other wavelengths, is predominantly composed by β mol and has a low percentage of β aer .
From the results obtained in both case studies, it is possible to observe the influence of the wavelength in the proposed methodology. The wavelength 1064 nm, considered as our signal reference, has a negligible influence of component molecular, 25 therefore the backscatter signal retrieved at 1064 nm can be considered approximately equal to the backscatter signal retrieved only by the aerosol contribution, β 1064 ≈ β aer . Before, taking into account the approximation demonstrated in equation 5 (RCS 1064 ≈ β 1064 ), we can conclude that the range corrected signal retrieved from a lidar at 1064 nm can be considered, in an good precision, approximately equal to the backscatter signal retrieved at the same wavelength for aerosol components , RCS 1064 ≈ β aer . Such relation enables the observation of behavior of aerosol plumes from high order moments. In the case of 30 wavelength 532 nm, β 532 is composed by β aer and β mol (β 532 = β aer532 + β mol532 ), however, as shown in the Figures A8 and   A13, there is a predominance of β aer . Although the high-order moments profiles obtained from the wavelength 532 nm are noisier than that one generated from the reference wavelength data, the phenomena observed from the 1064 nm data also can be observed in 532 nm data, mainly after the application of first lag correction. Consequently the wavelength at 532 nm can be used in the proposed methodology providing satisfactory results. On the other hand, the backscatter at 355 nm is predominantly composed by β mol and has a small percentage of β aer , as presented in figures A8 and A13.
This fact justifies the low quality observed in the results retrieved using the wavelength of 355 nm. As established in equation 3, the turbulent variable is directly associated with β aer , but due to low contribution of this component in the backscatter signal at 355 nm, the supposition established in equation 6 cannot be applied. Consequently, the high-order moments obtained from 5 the proposed methodology are noisier and the value of τ RCS (z) is almost half of the value obtained from the reference wavelength, both due to influence of β mol that presents the stronger contribution to the total backscatter coefficient at this wavelength. Therefore the behavior observed in the high-order moments profiles generated from the 1064 nm wavelength data can be detected partially, or even totally suppressed as the complexity of high-order moments increase. In the both case studies were possible to observe that from the third order moment (skewness) the results obtained from the wavelength 355 nm provide 10 misinformation.

Conclusions
In this paper we performed a comparative analysis about the use of different wavelengths (355, 532 and 1064 nm) in studies about turbulence. The data were acquired with an elastic lidar, from the SPU Lidar station of LALINET, by measurements gathered with high frequency (0.5 Hz) along July 2017 to July of 2018. The RCS provided by this system was used to calculate 15 high-order moments (variance, skewness and kurtosis) and the integral time scale, which were applied to characterization of aerosol dynamics. Based on previous studies, the wavelength 1064 nm was adopted as reference due to predominance of β aer .
Two case studies (26 th July 2017 and 19 th July 2018) were performed in order to verify the proposed methodology, as well as, the applicability of each wavelength. In both cases, the results obtained from 1064 nm wavelength demonstrate as the high-order moments can support a detailed analysis of the ABL region. In addition, it is remarkable the values of τ RCS in 20 the region below the ABLH, demonstrating the viability of the proposed methodology. The high-order moments obtained from the wavelength 532 nm are slightly more influenced by the noise than the results obtained from the reference wavelength (the value of noise can be observed by the ACF 532 . However, the same phenomena observed in the high-order moments profiles generated from the 1064 nm wavelength can be observed in that one generated from the wavelength 532 nm, mainly with the application of first lag correction. On the other hand, the high-order moments obtained from 355 nm have a strong presence of 25 noise and, thus, from the third order moment (skewness) the phenomenon presented in the high-order moments obtained from 1064 nm wavelength cannot be observed in 355 nm high-order moments profiles.
The analysis of the backscatter signal at each wavelength shows that for both case studies β aer is a predominant contribution at 532 nm, while β mol is predominant at 355 nm. In this way, the high-order statistics become noisier at 355 nm, and it cannot be applied in the proposed methodology. In contrast, the predominance of β aer at 532 nm implicates that this wavelength 30 provides results similar to that obtained at 1064 nm, especially after the application of first lag correction. Consequently, the 532 nm wavelength can be used to apply the proposed methodology, providing results similar to that obtained from 1064 nm wavelength.
The results obtained in this paper show the viability of the proposed methodology and its applicability to the 532 nm wavelength, due to the similarity with results derived at 1064 nm and the evidence of a low ε influence. On the other hand, the wavelength 355 nm does not provide satisfactory results in such methodology due to predominance of molecular signal in its composition. However, a better assessment of the molecular backscatter at 355 can reduce the influence of the noise caused by molecular signal and improve the results obtained from the data generated from this channel. In addition, the high-5 order moments obtained from the SPU Lidar station using an elastic lidar data provided us detailed information about some phenomenon in the ABL, allowing us a better comprehension about the aerosol dynamics.   A1. Variables applied to statistical analysis of turbulence in APBL region (Lenschow et al., 2000). The sum of subindex of autocovariance function Mij represents the order of it.            Figure A13. Total (aerosol and molecular) backscatter profile and backscatter ratio retrieved using Klett-Fernald-Sasano inversion technique for 1064, 532 and 355 nm, respectively, for data retrieved on 19 th July 2018 from 18:00 to 19:00 UTC.