Analyzing the turbulent planetary boundary layer by remote sensing systems: the Doppler wind lidar, aerosol elastic lidar and microwave radiometer

The planetary boundary layer (PBL) is the lowermost region of troposphere and is endowed with turbulent characteristics, which can have mechanical and/or thermodynamic origins. This behavior gives this layer great importance, mainly in studies about pollutant dispersion and weather forecasting. However, the instruments usually applied in studies of turbulence in the PBL have limitations in spatial resolution (anemometer towers) or temporal resolution (instrumentation aboard an aircraft). Ground-based remote sensing, both active and passive, offers an alternative for studying the PBL. In this study we show the capabilities of combining different remote sensing systems (microwave radiometer – MWR, Doppler lidar – DL – and elastic lidar – EL) for retrieving a detailed picture on the PBL turbulent features. The statistical moments of the high frequency distributions of the vertical wind velocity, derived from DL, and of the backscattered coefficient, derived from EL, are corrected by two methodologies, namely first lag correction and $-\mathrm{2}/\mathrm{3}$ law correction. The corrected profiles, obtained from DL data, present small differences when compared with the uncorrected profiles, showing the low influence of noise and the viability of the proposed methodology. Concerning EL, in addition to analyzing the influence of noise, we explore the use of different wavelengths that usually include EL systems operated in extended networks, like the European Aerosol Research Lidar Network (EARLINET), Latin American Lidar Network (LALINET), NASA Micro-Pulse Lidar Network (MPLNET) or Skyradiometer Network (SKYNET). In this way we want to show the feasibility of extending the capability of existing monitoring networks without strong investments or changes in their measurements protocols. Two case studies were analyzed in detail, one corresponding to a well-defined PBL and another corresponding to a situation with presence of a Saharan dust lofted aerosol layer and clouds. In both cases we discuss results provided by the different instruments showing their complementarity and the precautions to be applied in the data interpretation. Our study shows that the use of EL at 532 nm requires a careful correction of the signal using the first lag time correction in order to get reliable turbulence information on the PBL.

1 Introduction

The planetary boundary layer (PBL) is the atmospheric layer directly influenced by the Earth's surface that responds to its changes within timescales of around an hour (Stull, 1988). This 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 moments after sunrise, the ground surface temperature increases due to the positive net radiative flux (R_n). This process intensifies the convection, where there is an ascension of warm air masses, causing the downward displacement of colder air masses and consequently causing the convective boundary layer (CBL) or mixing layer (ML). This 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 R_n 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 CBL leads to the development of two layers, namely a stably stratified boundary layer called the stable boundary layer (SBL) close to the surface and the residual layer (RL) that contains features from the previous day's ML and is just above the SBL.

Knowledge of the turbulent processes in the CBL is important in diverse studies, mainly for atmospheric modeling and pollutant dispersion, since turbulent mixing can be considered the primary process by which aerosol particles and other scalars are transported vertically in the 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 an analysis provides information about the field of turbulent fluctuation as well as a description of the mixing process in the PBL (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 this equipment restricts the analysis to regions close to surface. Aircrafts have also been used in atmospheric turbulence studies (e.g., Lenschow et al., 1980, 1994; Williams and Hacker, 1992; 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 an adequate range are necessary in order to provide more detailed results along the day throughout the whole thickness of the PBL.

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 differential absorption lidar – DIAL; e.g., Wulfmeyer, 1999; Kiemle et al., 2007; Wulfmeyer et al., 2010; Turner et al., 2014; Muppa et al., 2016), temperature (rotational Raman lidar; e.g., Behrendt et al., 2015) and aerosol (elastic lidar; e.g., Pal et al., 2010; McNicholas and Turner, 2014). This allows for 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 PBL region. In addition, when different lidar systems operate synergistically, as in Engelmann et al. (2008), for example, 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 of characterizing the PBL turbulence by DL. Pal et al. (2010) have shown the feasibility of retrieving information on the PBL turbulence from high-order moments of elastic lidar operating at 1064 nm. These approaches are even more attractive when considering facilities of networks, e.g., the European Aerosol Research Lidar Network (EARLINET; Pappalardo et al., 2014), microwave radiometer network (MWRNET; Rose et al., 2005; Caumont et al., 2016) and European Research Infrastructure for the observation of Aerosol, Clouds and Trace gases (ACTRIS) CLOUDNET (Illingworth et al., 2007). For these reasons, and keeping in mind the wide spread of elastic lidar systems operated at other wavelengths like 532 or 355 nm, it would be a worthy test of the feasibility of these other wavelengths in the characterization of the turbulent PBL behavior.

The use of simple techniques, applied to the aforementioned remote systems, provides robust and similar information on the PBL height (PBLH) during the convective period (see, for example, Moreira et al., 2018) or complementary information when the CBL is substituted by the presence of the SBL and the RL (Moreira et al., 2019). Thus, the combination of information obtained from the active remote sensing systems, Doppler lidar (DL) and elastic lidar (EL), acquired with a temporal resolution close to 1 s, and that provided by the microwave radiometer (MWR) can provide a detailed understanding about different features of the PBL, like its structure (CBL versus SBL and RL), the height of the layers, and the rate of growth of the PBLH and turbulence.

In this study we show the feasibility of obtaining clear insight into the PBL behavior using a combination of active and passive remote sensing systems (EL, DL and MWR) acquired during the Sierra Nevada Lidar AerOsol Profiling Experiment (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 EL at 532 nm, considering the widespread use of lidar systems based on laser emission at this wavelength in different coordinated networks, like EARLINET (Pappalardo et al., 2014) and LALINET (Latin American Lidar Network; Guerrero-Rascado et al., 2016). In addition, this study shows the variety of applications that can be done with EARLINET data, applying some simple changes to the data acquisition procedures.

This paper is organized as follows. Descriptions of the experimental site and the equipment setup are presented in Sect. 2. The methodologies applied are introduced in Sect. 3. Section 4 presents the results of the analyses using the different methodologies. Finally, conclusions are summarized in Sect. 5.

2 Experimental site and instrumentation

The SLOPE-I campaign was performed from May to September 2016 in southeastern Spain in the framework of the 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 IISTA-CEAMA and in situ measurements operating at different altitudes in the northern slope of the Sierra Nevada, around 20 km away from IISTA-CEAMA (Bedoya-Velásquez et al., 2018; Román et al., 2018). The IISTA-CEAMA station is part of the EARLINET (Pappalardo et al., 2014) since 2005 and, at present, is an ACTRIS station (https://www.actris.eu/DataServices/ObservationalFacilities/AccesstoObservationalFacilities.aspx, last access: 21 January 2019). The research facilities are located at Granada, a medium-sized city in southeastern Spain (Granada; 37.16^∘ N, 3.61^∘ W; 680 m a.s.l.), surrounded by mountains and with Mediterranean continental 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 a different kind of aerosol particles the originated locally and is affected by medium–long-range particles 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 the EARLINET. This system operates with a pulsed Nd:YAG laser, with its frequency doubled and tripled by potassium dideuterium phosphate crystals, emitting at wavelengths of 355, 532 and 1064 nm and 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 N₂), 408 (from H₂O) and 607 nm (from N₂). 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 in Guerrero-Rascado et al. (2008, 2009).

The Doppler lidar (Halo Photonics, model Stream Line XR) is also operated at IISTA-CEAMA. Since May 2016, this system works in a continuous and automatic mode. 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/, last access: 21 January 2019). Since November 2011, this system operates in an automatic and continuous mode at IISTA-CEAMA. 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: the K band (water vapor; frequencies of 22.24, 23.04, 23.84, 25.44, 26.24, 27.84 and 31.4 GHz) and V band (oxygen; frequencies of 51.26, 52.28, 53.86, 54.94, 56.66, 57.3 and 58.0 GHz). From these bands, it is possible to obtain profiles of water vapor and temperature with the 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 min. Relative humidity is monitored with an accuracy of ±3 %, and air surface temperature is acquired with an accuracy and precision of 0.6 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 W m⁻²), and the calibration factor of stability has been periodically checked against a reference CM-11 pyranometer (Antón et al., 2012).

3 Methodology

3.1 MWR data analysis

The MWR data are analyzed combining two algorithms, the parcel method (PM; Holzworth, 1964) and the temperature gradient method (TGM; Coen et al., 2014), in order to estimate the PBL height (PBLH_MWR) in convective and stable situations, respectively. The different situations are differentiated by comparing the surface potential temperature (θ(z₀)) with the corresponding vertical profile of θ(z) up to 5 km. The cases in which all the points in the vertical profile have values larger than θ(z₀) are labeled as stable, and the TGM is applied. Otherwise the situation is labeled as unstable, and the PM is applied. The vertical profile of θ(z) is obtained from the vertical profile of T(z) using the following equation (Stull, 2011):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\theta }\left(z\right)=T\left(z\right)+\mathrm{0.0098}\cdot z,\end{array}$$

where T(z) is the temperature profile provided by MWR, z is the height above sea level, and 0.0098 K m⁻¹ is the dry adiabatic temperature gradient. A meteorological station co-located with the MWR is used to detect the surface temperature (T(z₀)). In order to reduce the noise, θ(z) profiles were averaged providing a PBLH_MWR value at 30 min intervals. This methodology of PBLH detection was selected as the reference due to the results obtained during a performed intercomparison campaign between MWR and radiosonde data, where 23 radiosondes were launched. High correlations were found between PBLH retrievals provided by both instruments in stable and unstable cases. Further details are given by Moreira et al. (2018).

3.2 Lidar retrieval of the PBLH

The simple processing of DL and EL data allows for the estimation of the CBL height. Moreira et al. (2018) have discussed this issue in depth, while Moreira et al. (2019) have used the complementarity of the data obtained from distinct remote sensing systems in order to distinguish the sub-layers during the period when the SBL and RL substitute for the CBL as well as during complex situations, like the presence of dust layers.

The PBLH obtained from DL data (PBLH_Doppler) is estimated from the variance threshold method. In this method the PBLH_Doppler is attributed to the height where the variance of vertical wind speed (${\mathit{\sigma }}_{\mathrm{w}}^{\mathrm{2}}$) is lower than a determinate threshold, which was adopted as 0.16 m² s⁻² (Moreira et al., 2018). For the PBLH_Doppler, calculations were selected from a time interval of 30 min. Concerning the PBLH obtained from EL (PBLH_Elastic), the variance method is applied. This method assumes the maximum of the variance of range corrected signal (${\mathit{\sigma }}_{\mathrm{RCS}}^{\mathrm{2}}$) as PBLH_Elastic (Moreira et al., 2015). The ${\mathit{\sigma }}_{\mathrm{RCS}}^{\mathrm{2}}$ is obtained from a time interval of 30 min.

3.3 Lidar turbulence analysis

Both lidar systems, DL and EL, gathered data (q(z,t)) with a temporal resolution of 2 s. Then, the data are averaged in 1 h packages from which the mean value is extracted ($\stackrel{\mathrm{‾}}{q}\left(z\right)$). This mean value is subtracted from each q(z,t) profile in order to estimate the vertical profile of the fluctuation for the measured variable ($q^{\prime }\left(z,t\right)$; i.e., vertical velocity for the DL):

$$\begin{array}{}\text{(2)}& {\displaystyle }{q}^{\prime }\left(z,t\right)=q\left(z,t\right)-\stackrel{\mathrm{‾}}{q}\left(z\right).\end{array}$$

Then, from $q^{\prime }\left(z,t\right)$, it is possible to obtain the high-order moments (variance – σ², skewness – S – and kurtosis – K) as well as the integral timescale (τ, 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, M_ij:

$$\begin{array}{}\text{(3)}& {\displaystyle }{M}_{ij}=\underset{\mathrm{0}}{\overset{t_{\mathrm{f}}}{\int }}\left[q^{\prime }\left(z,t\right){]}^i\right[q^{\prime }\left(z,t+t_{\mathrm{f}}\right){]}^j\mathrm{d}t,\end{array}$$

where t_f is the final time, and i and j indicate the order of autocovariance function.

However, it is necessary to consider that the acquired real data contain instrumental noise, ε(z). Therefore, Eq. (3) can be rewritten as

$$\begin{array}{}\text{(4)}& {\displaystyle }{M}_{ij}=\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}[q\left(z,t\right)+\mathit{\epsilon }(zt\left){]}^i\right[q\left(z,t+\mathit{\tau }\right)+\mathit{\epsilon }(zt+\mathit{\tau }){]}^j\mathrm{d}t.\end{array}$$

The autocovariance function of a time series with zero lag results in the sum of the variances of the atmospheric variable and its ε(z). Nevertheless, atmospheric fluctuations are correlated in time, but the ε(z) 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) made the suggestion of obtaining the corrected autocovariance function, M₁₁(→0), from two methods, namely first lag correction or $-\mathrm{2}/\mathrm{3}$ law correction. In the first method, M₁₁(→0) is obtained directly by the subtraction of lag 0, ΔM₁₁(0), from the autocovariance function, M₁₁(0). In the second method, M₁₁(→0) is generated by the extrapolation of M₁₁(0) at the first nonzero lags back to lag zero ($-\mathrm{2}/\mathrm{3}$ law correction). The extrapolation can be performed using the inertial subrange hypothesis, which is described by the following equation (Monin and Yaglom, 1979):

$$\begin{array}{}\text{(5)}& {\displaystyle }{M}_{\mathrm{11}}\left(\to \mathrm{0}\right)=\stackrel{\mathrm{‾}}{{q^{\prime }}^{\mathrm{2}}(z,t)}+C{t}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where C represents a parameter of the turbulent eddy dissipation rate. The high-order moments and τ corrections and errors are shown in Table 1 (columns 2 and 3, respectively).

Table 1 Variables applied to statistical analysis (Lenschow et al., 2000).

Figure 1 Procedure to remove the errors of autocovariance functions. M₁₁(→0) is corrected autocovariance function errors, M₁₁(0) is autocovariance function without correction and ΔM₁₁(0) is error of autocovariance function.

The same procedure of analysis is applied in studies with DL and EL, with the main difference being the tracer used by each system, which is the fluctuation of vertical wind speed (w^′) for DL and aerosol number density (N^′) for EL. DL provides w(z,t) directly, and therefore the procedure described in Fig. 2 can be directly applied. Thus, the two corrections described above are applied separately, and τ and high-order moments with and without corrections can finally be estimated.

On the other hand, the EL does not provide N(z,t) 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, therefore adopting the following relation (Pal et al., 2010):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{par}}\left(z,t\right)\approx N\left(z,t\right)Y\left(z\right)\Rightarrow {\mathit{\beta }}_{\mathrm{par}}^{\prime }\left(z,t\right)=N^{\prime }\left(zt\right),\end{array}$$

where β_par and ${\mathit{\beta }}_{\mathrm{par}}^{\prime }$ represent the particle backscatter coefficient and its fluctuation, respectively, and Y(z) does not depend on time.

Considering the lidar equation,

$$\begin{array}{}\text{(7)}& {\displaystyle }{P}_{\mathit{\lambda }}\left(z\right)=P_{\mathrm{0}}{\displaystyle \frac{c{t}_{\mathrm{d}}}{\mathrm{2}}}AO\left(z\right){\displaystyle \frac{{\mathit{\beta }}_{\mathit{\lambda }}\left(z\right)}{z^{\mathrm{2}}}}{e}^{-\mathrm{2}{\int }_{\mathrm{0}}^z{\mathit{\alpha }}_{\mathit{\lambda }}\left(z^{\prime }\mathrm{d}{z}^{\prime }\right)},\end{array}$$

where P_λ(z) is the signal returned from distance z at time t, z is the distance (m) from the lidar of the volume investigated in the atmosphere, P₀ is the power of the emitted laser pulse, c is the light speed (m s⁻¹), t_d is the duration of laser pulse (ns), A is the area (m²) of telescope cross section, O(z) is the overlap function, α_λ(z) is the total extinction coefficient (due to atmospheric particles and molecules; (km)⁻¹) at distance z, β_λ(z) is the total backscatter coefficient (due to atmospheric particles and molecules; (km sr)⁻¹) at distance z and the subscript λ represents the wavelength. The two-path transmittance term related to α(z) is considered nearly negligible at 1064 nm (Pal et al., 2010). Thus, it is possible to affirm that

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathrm{RCS}}_{\mathrm{1064}}\left(z\right)=P{\left(z\right)}_{\mathrm{1064}}{z}^{\mathrm{2}}\cong G.{\mathit{\beta }}_{\mathrm{1064}}\left(z\right),\end{array}$$

and, consequently, that

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{RCS}}_{\mathrm{1064}}^{\prime }(z,t)\cong {\mathit{\beta }}_{\mathrm{1064}}^{\prime }\left(z,t\right)={{\mathit{\beta }}^{\prime }}_{\mathrm{par}}\left(z,t\right)=N^{\prime }\left(z,t\right),\end{array}$$

where RCS₁₀₆₄ and RCS${}_{\mathrm{1064}}^{\prime }$ are the range corrected signal and its fluctuation, respectively, G is a constant, and the subscripts represent the wavelength.

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

In this way, Pal et al. (2010) have shown the feasibility of using EL operating at 1064 nm for describing the atmospheric turbulence. However, keeping 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 use RCS₅₃₂ fluctuations to determine turbulence, following the procedure described in Fig. 3. This EL methodology is very similar to that described earlier for DL.

Figure 3 Flowchart of data analysis methodology applied to the study of turbulence with elastic lidar.

We perform the validation of the RCS₅₃₂ in analyses about turbulence using EL, following the procedure described in Fig. 3, which is basically the same methodology described earlier for DL.

4 Results

4.1 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, for the first step, in order to ascertain the applied methodology and our data quality, we performed the error treatment of DL data, as described in Fig. 2. For the DL analysis we selected the period 08:00–09:00 UTC on 19 May, the same day that will be presented in Case Study 1. This day is characterized by a well-defined PBL.

Figure 4 illustrates the autocovariance function, generated from w^′, 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 PBLH_MWR (1100 m a.g.l. in our example).

Figure 4 Autocovariance function (ACF) of w^′, obtained from Doppler lidar at three different heights on 19 May 2016 from 08:00–09:00 UTC in Granada.

Figure 5 (a) Vertical profile of integral timescale (${\mathit{\tau }}_{w^{\prime }}$). (b) Vertical profile of variance (${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$). (c) Vertical profile of skewness ($S_{w^{\prime }}$). All profiles were obtained from Doppler lidar data on 19 May 2016 from 08:00–09:00 UTC in Granada.

Figure 5a illustrates the comparison between the integral timescale (${\mathit{\tau }}_{w^{\prime }}$) without correction and the two corrections cited in Sect. 3.2. Except for the first height bins, below the PBLH_MWR, the profiles have little differences as well as small error bars. Above the PBLH_MWR the first lag correction presents higher differences in relation to the other profiles at around 1350 m.

Figure 5b and c show the comparison of variance (${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$) and skewness ($S_{w^{\prime }}$), respectively, with and without corrections. The profiles corrected by $-\mathrm{2}/\mathrm{3}$ law correction 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 PBLH_MWR, mainly the ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ ($S_{w^{\prime }}$ only in the first 50 m). Therefore, considering high signal-to-noise ratio (SNR) conditions, although the presence of ε can slightly change the value of high-order moments, it is not enough to distort the observed phenomena, as shown by the impact of the corrections applied.

For EL we use the same procedure for the correction and error analysis that we apply to the DL data. The same day was chosen (19 May), however the period selected is between 12:00 and 13:00 UTC, due to the incomplete overlap of MULHACÉN.

In this sense, we studied the influence of noise at two wavelengths: 1064 nm, which has been previously analyzed by Pal et al. (2010), as presented in the Sect. 2 and adopted as a reference (considering the rather low impact of molecular signal and the two-way transmittance shown in Eq. 9), and 532 nm was only used 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 show the autocovariance function, obtained from ${{\mathrm{RCS}}^{\prime }}_{\mathrm{1064}}$ and ${{\mathrm{RCS}}^{\prime }}_{\mathrm{532}}$, respectively, at three distinct heights. As expected, ε increases with the range, principally above the PBLH_MWR. However, the wavelength 532 nm is more influenced by the noise, which can be verified by the higher peak at lag 0 in Fig. 7, in comparison with peaks at same lag in Fig. 6.

Figure 6 Autocovariance of RCS^′₁₀₆₄ obtained from MULHACÉN elastic lidar data for three different heights on 19 May 2016 from 12:00–13:00 UTC in Granada.

Although the level of influence of ε in each wavelength depends on the SNR of them (which is associated to technical factors such as laser output power, filters and the type of detectors), considering the proposed methodology, evaluating the composition of each wavelength is also important. The large contribution of ${\mathit{\beta }}_{\mathrm{Molecular}}^{\mathrm{532}}$ to the total β at 532 nm, in comparison with the behavior at 1064 nm, can influence the results obtained from this wavelength, because our methodology is based on the use of ${\mathit{\beta }}_{\mathrm{Aerosol}}^{\prime }$. 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 SNR values at this wavelength. Since we used the elastic lidar technique, we could not calculate aerosol extinction profiles, but an estimation of these transmittances was done on the basis of the Klett method (Klett, 1985). With this method, a constant lidar ratio value was constrained for each profile using the aerosol optical depth (AOD) derived from a collocated AERONET sunphotometer (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 May).

Figure 8a, b, c and d show the vertical profiles of ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ (a.u. – arbitrary units), ${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$, $S_{{\mathrm{RCS}}^{\prime }}$ and kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$), respectively, obtained at 1064 nm, with and without the corrections described in Sect. 3.2. In general, the corrections do not affect the profiles generated from 1064 nm data in a significant way, so the higher influence of corrections is observed in the $K_{{\mathrm{RCS}}^{\prime }}$ profile, which is underestimated in some regions. In the Fig. 9a, b, c and 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 the first lag correction, make these profiles very similar to those generated from the wavelength 1064 nm, so the same phenomena can be observed in both.

Therefore, in spite of the larger attenuation expected at 532 nm wavelength, which reduces the SNR of the profiles in comparison with 1064 nm, the application of the proposed corrections, mainly the first lag, reduces this influence significantly 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 Sect. 4.2. The first lag correction was adopted as default due to it generates more relevant results in comparison with the $-\mathrm{2}/\mathrm{3}$ law correction, providing a more careful analysis.

Figure 7 Autocovariance of RCS^′₅₃₂ obtained from MULHACÉN elastic lidar data for three different heights on 19 May 2016 from 12:00–13:00 UTC in Granada.

Figure 8 (a) Vertical profile of variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$). (b) Vertical profile of integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$). (c) Vertical profile of skewness ($S_{{\mathrm{RCS}}^{\prime }}$). (d) Vertical profile of kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$). All profiles were obtained from MULHACÉN elastic lidar data on 19 May 2016 in Granada from 12:00–13:00 UTC.

Figure 9 (a) Vertical profile of variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$). (b) Vertical profile of integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$). (c) Vertical profile of skewness ($S_{{\mathrm{RCS}}^{\prime }}$). (d) Vertical profile of kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$). All profiles were obtained from MULHACÉN elastic lidar data on 19 May 2016 in Granada from 12:00–13:00 UTC.

4.2 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 PBL. 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.

Table 2 Products and their respective meaning, provided by each system.

4.2.1 Case study I: clear sky situation

In this case study we use measurements gathered with DL, the MWR and the pyranometer over 24 h. The EL was operated under an operator-supervised mode between 08:20 to 18:00 UTC.

Figure 10a shows the integral timescale obtained from DL data (${\mathit{\tau }}_{w^{\prime }}$). The gray area represents the region where it is not possible to analyze the turbulent process from our DL data, either because of the low SNR values, which result in null values of the ${\mathit{\tau }}_{w^{\prime }}$, or due to regions where the ${\mathit{\tau }}_{w^{\prime }}$ is not null but is lower than the acquisition time of the DL. However, the gray area is located almost entirely above the PBLH_MWR (white stars).

Figure 10 (a) Integral timescale obtained from Doppler lidar data (${\mathit{\tau }}_{w^{\prime }}$), (b) variance obtained from Doppler lidar data (${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$), (c) skewness obtained from Doppler lidar data ($S_{w^{\prime }}$), (d) net radiation obtained from pyranometer data (R_n), and  (e) air surface temperature (blue line) and surface relative humidity (RH – orange line) were both obtained from surface sensors. All profiles were acquired on 19 May 2016 in Granada. In (a), (b) and (c), black lines and white stars represent air temperature and PBLH_MWR, respectively.

The ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ has low values during the entire period when the SBL is present (Fig. 10b). Nevertheless, as air temperature begins to increase (around 07:00 UTC), the ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ increases along with the PBLH_MWR. The ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ reaches its maximum values in the middle of the day, when we also observe the maximum values of air temperature and PBLH_MWR. The combination of the ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ and PBLH_MWR provides us with better comprehension about the PBLH growth speed, so in the moments where high values of ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ are observed, it means higher values of turbulent kinetic energy (TKE), which favor the fast ascension of PBLH.

The skewness of w^′ ($S_{w^{\prime }}$) is shown in Fig. 11c. The $S_{w^{\prime }}$ describes the distribution of the turbulent velocities. Thus a positive $S_{w^{\prime }}$ implies strong but narrow updrafts surrounded by weaker but more widespread downdrafts, and vice versa for negative $S_{w^{\prime }}$. 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 $S_{w^{\prime }}$. Nevertheless, as air temperature increases (transition from stable to unstable period), $S_{w^{\prime }}$ values begin to become larger. Air temperature begins to decrease around 18:00 UTC, and there is a reduction of $S_{w^{\prime }}$, so the generation rate of convective turbulence decreases. Therefore, if the turbulence cannot be maintained against dissipation, then the CBL becomes an SBL covered by the RL. Thus, the reduction observed in the PBLH_MWR is due to the detection of SBL height.

Figure 11 Time–height plot of RCS obtained on 19 May 2016 in Granada. Pink stars represent the PBLH_MWR, black stars represent the PBLH_Elastic and blue stars represent the PBLH_Doppler.

Figure 10d shows the values of net surface radiation (R_n) that are estimated from solar global irradiance values using the seasonal model described in Alados et al. (2003). The negative values of R_n are concentrated in the stable region. The R_n begins to increase around 06:00 UTC and reaches its maximum in the middle of the day. Comparing Fig. 8c and d, we can observe similarity among the behavior of $S_{w^{\prime }}$ and R_n, so the joint analysis of these variables reinforce the characterization of this PBL as a surface-heating-driven CBL.

Figure 10e presents the values of surface air temperature and surface relative humidity (RH). Air surface temperature has a pattern of increase and decrease similar to that observed in R_n and $S_{w^{\prime }}$. On the other hand, the RH is inversely correlated with temperature.

Figure 12 Statistical moments obtained from 532 nm wavelength data of elastic lidar (MULHACÉN) in Granada from 13:00 to 14:00 UTC on 19 May 2016. From left to right is the variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$), integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$]), skewness ($S_{{\mathrm{RCS}}^{\prime }}$) and kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$).

Figure 11 shows the RCS₅₃₂ profile obtained from 08:00 to 18:00 UTC. At the beginning of the measurement period (08:20 to 10:00 UTC), observing the presence of a thin residual layer (around 2000 m a.s.l.) is possible, and later from 13:00 to 18:00 UTC, a lofted aerosol layer is evident. In this plot there are the PBLH_MWR values (pink stars), the PBLH_Doppler values (blue stars), obtained from the maximum of ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ (Moreira et al., 2018), and the PBLH_Elastic values (black stars), obtained from the maximum of ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ (Moreira et al., 2015). In the initial part of measurement, all profiles have similar behavior. However due to the distinct PBLH definition and tracer applied by each profile, the differences increase as CBL becomes more complex, e.g., with the presence of the lofted aerosol layer at 14:00 UTC. The joint observation of the results provided by these three methods can provide us information about the sub-layers in the PBL, both in convective and stable situations. Due to low variability of PBLH, the period between 13:00 and 14:00 UTC has been selected from the high-order moments to be analyzed.

Figure 12 presents the statistical moments generated from the RCS^′ of wavelength 532 nm which were obtained between 13:00 and 14:00 UTC. The red line in all graphics represents the PBLH_Elastic (2200 m a.s.l.), and the blue one represents the average value of PBLH_MWR (2250 m a.s.l.), both obtained between 13:00 and 14:00 UTC.

Due to the presence of a decoupled aerosol layer at 13:30 UTC, the average values of PBLH_Elastic and PBLH_MWR have a difference of around 500 m. The ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ has small and practically constant values between 1000 and 1400 m, evidencing the homogeneity of aerosol distribution in this region. From 1400 m the value of ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ begins to increase, reaching a positive peak at PBLH_MWR, which represents the entrainment zone (region characterized by an intense mixing between air parcels coming from the CBL and free troposphere – FT, causing a high variation in aerosol concentration). The PBLH_Elastic values observed at approximately 2900 m demonstrate an inherent difficulty of the variance method for detecting the PBLH in the presence of several aerosol layers (Kovalev and Eichinger, 2004). Above the PBLH_Elastic the values of ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ decrease slowly due to location of the lofted aerosol around 2500 m. However, above this aerosol layer the value of ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ is reduced to zero, indicating a large homogeneity in aerosol distribution at this region, which is expected, because the aerosol concentration at the FT is negligible in this case. The integral timescale obtained from RCS^′ (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$) has values higher than EL time acquisition throughout the CBL, evidencing the feasibility of studying turbulence using this elastic lidar configuration. The skewness values obtained from RCS^′ ($S_{{\mathrm{RCS}}^{\prime }}$) give us information about aerosol motion. The positive values of $S_{{\mathrm{RCS}}^{\prime }}$ observed in the lowest part of profile and above the PBLH_Elastic represent the updrafts of aerosol layers. The negative values of $S_{{\mathrm{RCS}}^{\prime }}$ indicate the region with low aerosol concentration due to clean air coming from the FT. This movement of the ascension of aerosol layers and the descent of clean air with a zero value of $S_{{\mathrm{RCS}}^{\prime }}$ at PBLH (characteristic of the CBL growing) was also detected by Pal et al. (2010) and McNicholas and Turner (2014). The kurtosis of RCS^′ ($K_{{\mathrm{RCS}}^{\prime }}$) determines the level of mixing at different heights. There are values of $K_{{\mathrm{RCS}}^{\prime }}$ larger than 3 in the lowest part of profile and around 2500 m, showing a peaked distribution in this region. On the other hand, values of $K_{{\mathrm{RCS}}^{\prime }}$ lower than 3 are observed close to the PBLH_Elastic, therefore this region has a well-mixed CBL regime. Pal et al. (2010) and McNicholas and Turner (2014) also detected this feature in the region near the PBLH. In Fig. 13 the high-order moments obtained at the same period described above are shown, however they are shown from the 1064 nm data (our reference wavelength). It is possible to observe a similarity between the profiles obtained from each wavelength, so the same phenomena observed in the profiles generated from 532 nm and described above are also detected in the profiles obtained from the reference wavelength.

Figure 13 Statistical moments obtained from 1064 nm wavelength data of elastic lidar (MULHACÉN) in Granada from 13:00–14:00 UTC on 19 May 2016. From left to right is the variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$), integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$), skewness ($S_{{\mathrm{RCS}}^{\prime }}$) and kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$).

The results provided by DL, pyranometer and MWR data agree with the results observed in Figs. 12 and 13. In the same way, the analysis of high-order moments of RCS^′ fully agree with the information in Fig. 10. Thus, the large values of $S_{{\mathrm{RCS}}^{\prime }}$ and $K_{{\mathrm{RCS}}^{\prime }}$ detected around 2500 m a.s.l, where we can see a lofted aerosol layer, suggest the ascent of an aerosol layer and the presence of a peaked distribution, respectively.

4.2.2 Case study: dusty and cloudy scenario

In case study measurements with DL, the MWR and pyranometer expand over 24 h, while EL data are collected from 09:00 to 16:00 UTC.

Figure 14a shows ${\mathit{\tau }}_{w^{\prime }}$. Outside of the period from 13:00 to 17:00 UTC, the greatest part of gray area is situated above the PBLH_MWR (white stars), thus DL time acquisition is enough to perform studies about turbulence in this case.

Figure 14 (a) Integral timescale from Doppler lidar data (${\mathit{\tau }}_{w^{\prime }}$), (b) variance from Doppler lidar data (${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$), (c) skewness from Doppler lidar data ($S_{w^{\prime }}$), (d) net radiation from pyranometer data (R_n), and (e) air surface temperature (blue line) and surface relative humidity (RH – orange line) from surface sensor data. All profiles were obtained in Granada on 8 July 2016. In (a), (b) and (c), black lines and white stars represent air temperature and PBLH_MWR, respectively.

${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ has values close to zero during the entire stable period (Fig. 14b). However, when air temperature begins to increase (around 06:00 UTC), the ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ also increases and reaches its maximum in the middle of the day. The higher values of PBL growth speed are observed in the moments where ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ reaches its maximum values. In the late afternoon, as air temperature decreases, the values of ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ (and consequently the TKE) decrease gradually until they reach the minimum value associated with the SBL. Figure 14c shows the profiles of $S_{w^{\prime }}$. The main features of this case are the low values of $S_{w^{\prime }}$, the slow increase and ascension of positive $S_{w^{\prime }}$ values, and the predominance of negative $S_{w^{\prime }}$ values from 12:00 to 13:00 UTC. The first two features are likely due to the presence of the intense Saharan dust layer (Fig. 15), which reduces the transmission of solar irradiance and consequently the absorption of solar irradiance at the surface, generating a weak convective process. From Fig. 16 we can observe the presence of both middle-altitude clouds and very intense dust layers from 12:00 to 15:00 UTC. This combination contributes to the intense negative values of the ${\mathrm{S}}_{w^{\prime }}$ observed in this period until around 2 km, because, as mentioned previously, the${\mathrm{S}}_{w^{\prime }}$ is directly associated with the direction of turbulent movements, which can be characterized as cloud-top longwave radiative cooling during this period (Ansmann et al., 2010).

Figure 15 Time–height plot of RCS obtained from MULHACÉN elastic lidar data on 8 July 2016 in Granada. Pink stars represent the PBLH_MWR, black stars represent the PBLH_Elastic and blues stars represent the PBLH_Doppler.

The influence of the Saharan dust layer can also be evidenced by the R_n pattern (Fig. 14d), which maintains negative values until 12:00 UTC and reaches a low maximum value (around 200 W m⁻²). The observation of $S_{w^{\prime }}$ and R_n 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 of a case of the cloud-top longwave radiative cooling in the CBL. The air surface temperature and the RH (Fig. 14e) 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 RH are detected in coincidence with the maximum daily value of PBLH_MWR.

As mentioned before, Fig. 15 shows the RCS profile obtained from 09:00 to 16:00 UTC in a complex situation, with presence of the 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 PBLH_MWR, PBLH_Doppler and PBLH_Elastic, respectively. Due to the presence of dusty layers and clouds, the difference between the methods is more evident, mainly in the PBLH_Elastic, which uses aerosol as tracers. This method only produces results close to the others at 15:00 UTC, when the dust layer is mixed with the CBL.

Figure 16 Statistical moments obtained from 532 nm wavelength data of elastic lidar (MULHACÉN) in Granada from 11:00–12:00 UTC on 8 July 2016. From left to right is the variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$), integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$), skewness ($S_{{\mathrm{RCS}}^{\prime }}$) and kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$).

Figure 16 illustrates the statistical moments of RCS^′ of 532 nm wavelength obtained from 11:00 to 12:00 UTC. The ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ profile presents several peaks due to the presence of distinct aerosol sub-layers. The first peak is coincident with the value of PBLH_MWR. The value of PBLH_elastic is coincident with the base of the dust layer. This difficulty of detecting the PBLH in presence of several aerosol layers is inherent in the variance method (Kovalev and Eichinger, 2004). However, the joint observation of PBLH_MWR and PBLH_elastic enables us to characterize and distinguish the several sub-layers. The values of ${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$ are higher than EL acquisition time all along the PBL, evidencing the feasibility of EL time acquisition for studying the turbulence of PBL in this case. The $S_{{\mathrm{RCS}}^{\prime }}$ profile has several positive values, due to the large number of aerosol sub-layers that are present. The characteristic inflection point of $S_{{\mathrm{RCS}}^{\prime }}$ is observed to be coincident with the PBLH_MWR, which confirms the agreement between this point and the PBLH. From the analysis of $S_{{\mathrm{RCS}}^{\prime }}$ and $S_{w^{\prime }}$, it is possible to justify these phenomena from the mixing process demonstrated in the earlier case study. The $K_{{\mathrm{RCS}}^{\prime }}$ predominately has values lower than 3 below 2500 m, thus showing how this region is well mixed, as can be seen in Fig. 16. Values of $K_{{\mathrm{RCS}}^{\prime }}$ larger than 3 are observed in the highest part of profile, where the dust layer is located.

In order to show the feasibility of 532 nm wavelength, in the Fig. 17 are presented the high-order moments obtained between 11:00–12:00 UTC from 1064 nm wavelength data. Although the error of ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{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 (Figs. 16 and 17).

Figure 17 Statistical moments obtained from 1064 nm wavelength data of elastic lidar (MULHACÉN) in Granada from 11:00–12:00 UTC on 8 July 2016. From left to right is the variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$), integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$), skewness ($S_{{\mathrm{RCS}}^{\prime }}$) and kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$).

Figure 18 Statistical moments obtained from 532 nm wavelength data of elastic lidar (MULHACÉN) in Granada from 12:00–13:00 UTC on 8 July 2016. From left to right is the variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$), integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$), skewness ($S_{{\mathrm{RCS}}^{\prime }}$) and kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$).

Figure 19 Statistical moments obtained from 1064 nm wavelength data of elastic lidar (MULHACÉN) in Granada from 12:00–13:00 UTC on 8 July 2016. From left to right is the variance (${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$), integral timescale (${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$), skewness ($S_{{\mathrm{RCS}}^{\prime }}$) and kurtosis ($K_{{\mathrm{RCS}}^{\prime }}$).

Figure 18 shows the RCS^′ 532 nm wavelength high-order moments obtained from 12:00 and 13:00 in presence of cloud cover. The method based on the maximum of ${\mathit{\sigma }}_{{\mathrm{RCS}}^{\prime }}^{\mathrm{2}}$ locates the PBLH_Elastic at the cloud base, due to the high variance of RCS^′ generated by the clouds. The ${\mathit{\tau }}_{{\mathrm{RCS}}^{\prime }}$ presents values larger than EL time acquisition, therefore this configuration enables us to study turbulence by EL analyses. The $S_{{\mathrm{RCS}}^{\prime }}$ has few peaks, due to the mixing between CBL and dust layer that generates a more homogenous layer. The highest values of $S_{{\mathrm{RCS}}^{\prime }}$ are observed in regions where there are clouds, and the negative values (between 3500 and 4000 m) occur due to presence of air from the FT between the two aerosol layers (Fig. 15). The inflection point of the $S_{{\mathrm{RCS}}^{\prime }}$ profile is observed in the PBLH_MWR region. The $K_{{\mathrm{RCS}}^{\prime }}$ profile has low values in most of the PBL, demonstrating the high level of mixing during this period, where the dust layer and PBL are combined. The higher values of $K_{{\mathrm{RCS}}^{\prime }}$ are observed in the region of clouds. In the same way as the previous analysis, the high-order moments of the period mentioned above were calculated for the wavelength of 1064 nm (Fig. 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 the same phenomena can be observed, demonstrating the viability of the 532 nm wavelength in the proposed methodology.

5 Conclusions

In this paper we perform an analysis of the PBL turbulent features from three different types of remote sensing systems (DL, EL and MWR) and surface sensors during the SLOPE-I campaign. We applied two kind of corrections to the lidar data: the first lag correction and $-\mathrm{2}/\mathrm{3}$ law correction. The corrected DL statistical moments showed little variation with respect to the uncorrected profiles, denoting a rather low influence of the noise on high SNR conditions. The EL high-order moments were obtained from two wavelengths: 1064 nm, adopted as reference, and 532 nm, adopted in order to verify the viability to use the last one in turbulence analysis. From this comparison, it 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 this influence, so the same phenomena can be observed in the high-order moments provided from both wavelengths.

The case studies present two kind of situations: a well-defined PBL and a more complex situation with the presence of Saharan dust layer and some clouds. In both cases, it 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 compliment 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 to retrieve a detailed picture on the PBL turbulent features. In addition, the feasibility of using the analyses of high-order moments of the RCS collected at 532 nm at a temporal resolution of 2 s offers the possibility of using the proposed methodology in networks such as EARLINET or LALINET with a reasonable additional effort.