**Sub-Mode Aerosol Volume Size Distribution and Complex Refractive Index from the Three-Year Ground-Based Measurements in Chengdu China**

**Chi Zhang 1,2, Ying Zhang 1, Zhengqiang Li 1,\*, Yongqian Wang 3, Hua Xu 1, Kaitao Li 1, Donghui Li 1, Yisong Xie <sup>1</sup> and Yang Zhang <sup>3</sup>**


Received: 9 January 2019; Accepted: 22 January 2019; Published: 26 January 2019

**Abstract:** Chengdu is a typical basin city of Southwest China with rare observations of remote sensing measurements. To assess the climate change and establish a region aerosol model, a deeper understanding of the separated volume size distribution (VSD) and complex refractive index (CRI) is required. In this study, we employed the sub-mode VSD and CRI in Chengdu based on the three years observation data to investigate the sub-mode characteristics and climate effects. The annual average fraction of the fine-mode aerosol optical depth (AODf) is 92%, which has the same monthly tendency as the total AOD. But the coarse-mode aerosol optical depth (AODc) has little variation in different months. There are four distinguishing modes of VSD in Chengdu; the median radii are 0.17 μm ± 0.05, 0.31 μm ± 0.12, 1.62 μm ± 0.45, 3.25 μm ± 0.99, respectively. The multi-year average and seasonal variations of fine- and coarse-mode VSD and CRI are also analyzed to characterize aerosols over this region. The fine-mode single scattering albedos (SSAs) are higher than the coarse-mode ones, which suggests that the coarse-mode aerosols have a stronger absorbing effect on solar light than the small-size aerosol particles in Chengdu.

**Keywords:** Chengdu; aerosol; sub-mode volume size distribution; sub-mode complex refractive index; remote sensing

#### **1. Introduction**

As one component of the terrestrial atmosphere, aerosol is an important factor in global climate change, with direct effects and indirect effects. Direct effects include changing the radiation balance of the Earth-atmosphere system by absorbing and scattering sunlight [1,2]. On the other hand, aerosol serves as cloud condensation nuclei (CCN) involved in cloud microphysics processes, referred to as indirect effects [3]. Moreover, as one major constituent of haze, aerosol particles endanger public health [4].

Chengdu, located in the central region of the Sichuan basin (Chengdu Plain) with a population of ~16 million (the resident population of Chengdu was counted in 2016), is the economic center and transportation hub in southwest China. Due to the special topography of the basin, the average wind speed is low [5,6]. Coupled with the development economy in Chengdu, a large amount of pollutants associated with regional industrial emissions are easy to accumulate. Abundant water vapor

is conducive to the formation of aerosol particulates, which has made Chengdu become one of the regions in China with serious haze pollution [7–10]. The special topography and climate in Sichuan basin leads to a kind of wet, rainy and cloudy weather condition [11]. Therefore, it is difficult to obtain a large number of measurements through optical remote sensing observations in Chengdu. Sampling method of atmospheric particulates is one kind of Chengdu air pollution study [7,9,10,12], which fails to acquire the aerosol's larger-scale spatial characteristics. The spatial characteristics of aerosols can be obtained by remote sensing, one kind of non-contact method [13–15]. For instance, some studies analyzed the spatial-temporal distribution of aerosol optical depth (AOD) in Sichuan basin using satellite remote sensing [16,17]. Besides, aerosol optical parameters obtained by the ground-based remote sensing observation were statistically analyzed in this area [18–21].

Although the above-mentioned research studies focus on the total-column aerosol properties by remote sensing methods, it still remains a big challenge to obtain the characteristics of sub-mode aerosol (fine- and coarse-mode) in Chengdu. The fine- and coarse-mode particles are related to the different sources and compositions. Fine-mode aerosols are mainly dominated by anthropogenic emissions and are related to the fog or low-altitude cloud dissipated events [22,23]. Coarse-mode particles are determined by the sea salt, dust and other natural sources. Therefore, the knowledge of aerosol sub-mode properties plays a role in the research on regional climate and the improvement of the aerosol model used in the satellite retrieval algorithms.

In this study, the sub-mode volume size distribution (VSD) and complex refractive index (CRI) in Chengdu area were investigated based on three years of measurements of the Sun-sky radiometer Observation NETwork (SONET). Section 2 presents the observation data and the method to separate VSD and CRI into sub mode (fine and coarse mode). Section 3 shows the seasonal characteristics of sub-mode VSD and CRI. The sub-mode optical properties (AOD, single scattering albedo (SSA)) and the climate effect of sub-mode aerosol in Chengdu are discussed in Section 4. The results are summarized in Section 5.

#### **2. Data and Method**

#### *2.1. Observation Site and Data*

SONET Chengdu site (Figure 1 red dot) is located in Chengdu, Sichuan province, a developing urban area. The CE318 (Cimel Electronique, Pairs France) is an automatic instrument for long-term continuous observation of direct solar radiation and diffused sky radiation in the Chengdu site with nine individual spectral channels: 340 nm, 380 nm, 440 nm, 500 nm, 675 nm, 870 nm, 1020 nm, 1640 nm and 940 nm. The parameters include AOD, Angstrom Exponent (AE), VSD, CRI and SSA, and so on [24–26]. To ensure data quality, SONET instruments are calibrated by a routine process of laboratory and field calibration experiments. The direct sun measurements are calibrated every year in a field, located on Ling Mountain (~1600 m, MSL), and compared with a master instrument that is regularly calibrated by Langley plot method with high precision [27,28]. For the calibration of sky radiance measurements, SONET adopts the method of vicarious/transfer calibration, which is different from the AErosol RObotic NETwork (AERONET) calibration method [29].

The accuracy assessment of SONET products is based on the Distributed Regional Aerosol Gridded Observational Network (DRAGON)–Korea–United States Air Quality Study (KORUS-AQ) 2016 campaign [30]. The aerosol optical and microphysical parameters were inverted by SONET and AERONET algorithm at the same time. The differences of two network products can reflect the accuracy and data acquisition of SONET products. It turned out that the Level 1.5 data amount of two kinds of products are all above 85%, and the average AOD difference between SONET and AERONET is 0.002 ± 0.0001, less than AERONET AOD uncertainty [28]. The VSD difference between SONET and AERONET is 1.5% ± 26% (radius range from 0.1 to 7 μm) and 18% ± 85% (radius more than 7 μm) [28]. The real part of CRI has a difference of 0.007 ± 0.04 from AERONET, and the imaginary part has a difference value of 18% ± 46%. The average difference on SSA is slightly higher, but other parameters are close to or less than AERONET normal uncertainties [28].

The products of SONET are graded into three level (Level 1.0, Level 1.5, Level 2.0), defined following the AERONET data level protocols of version 2.0. In detail, Level 1.0 is raw data calculated from measurement and calibration coefficient. The Level 1.5 is based on Level 1.0 with automatic cloud screening procedures [31]. Level 2.0 has the additional application of pre- and post-calibration coefficient and expert checking. This paper is based on Level 2.0 data of SONET Chengdu site from June 2013 to December 2016.

ˊ **Figure 1.** Topographic map of Sichuan province and location of the Sun-sky radiometer Observation NETwork (SONET) Chengdu site (marked with red dot). Data source is ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model).

#### *2.2. Method*

In this paper, the sub-mode VSDs and CRIs were retrieved by the aerosol products [24,26], following by Zhang et al. (2016) and Zhang et al. (2017). The method developed by Cuesta et al. (2008) is employed to separate VSD into a single Log-Normal Modes (LNM). Each LNM parameters have three parameters: the modal concentration, geometric standard deviation and modal radius [32,33]. The sub-mode VSD can be modeled by the following function:

$$\frac{\text{dV}(\mathbf{r})}{\text{dlnr}} = \sum\_{\mathbf{i}=1}^{\text{m}} \frac{\mathbf{C}\_{\mathbf{i}}}{\sqrt{2\pi}|\ln \sigma\_{\mathbf{i}}|} \exp\left[-\frac{1}{2} \left(\frac{\ln \mathbf{r} - \ln \mathbf{r}\_{\mathbf{i}}}{\ln \sigma\_{\mathbf{i}}}\right)^{2}\right] \tag{1}$$

where Ci μm3/μm2 is the volume modal concentration, ri(μm) is the median radius, σ<sup>i</sup> is standard deviation, m is the total number of modes, dV/dlnr μm3/μm2 is aerosol volume particle size distribution.

The modes of VSD with radius less than 1 μm can be considered as the fine mode and others belong to the coarse mode.

CRI can describe the scattering and absorption properties of atmospheric particulates. Most research analyzed aerosol total-columnar CRI, but fine- and coarse-mode particles are associated with different composition and source of pollution. So, in this study we recalculated the complex refractive indices for both of fine and coarse mode, following Zhang et al. (2017). For the calculation of the sub-mode CRI, we also choose the same radius limit (1 μm) as the sub-mode VSD. The separated results of CRI have their own fine and coarse modes in different wavelength as follow:

$$\mathfrak{n}\_{\mathfrak{f}/\mathbb{c}}(\lambda) = \mathfrak{n}\_{\mathfrak{f}/\mathbb{c}}\,\lambda = 440\,\text{,}\,675\,\text{,}\,870\,\text{,}\,1020\,\text{nm}\tag{2}$$

$$\mathbf{k}\_{\mathbf{f},\mathbf{c}}(\lambda) = \begin{cases} \mathbf{k}\_{\mathbf{f},\mathbf{c}} \, \, 440 & \lambda = 440 \text{ nm} \\\ \mathbf{k}\_{\mathbf{f},\mathbf{c}} & \lambda = 675,870,1020 \text{ nm} \end{cases} \tag{3}$$

where n is real part of CRI, k is imaginary part of CRI, λ denotes the standard wavelengths of AERONET products, the subscripts f and c represent the fine and coarse modes, respectively.

The input parameters are the VSD, spectral AOD, and absorbing AOD. The initial guesses of sub-mode CRIs are from the inversion CRIs of measurements [34]. The effective CRIs are corresponding to each VSD bin, following the volume average rule [35]:

$$\mathbf{n}(\mathbf{r}) = \frac{n\_f V\_f(r) + n\_c V\_c(r)}{V\_f(r) + V\_c(r)} \tag{4}$$

$$\mathbf{k}(\lambda, \mathbf{r}) = \frac{k\_f V\_f(r) + k\_c V\_c(r)}{V\_f(r) + V\_c(r)} \tag{5}$$

Then the fine- and coarse-mode CRIs are found by iterative fitting of the input AODs and the calculated AODs by the CRIs (Equations (4) and (5)) and VSDs.

With regard to the test of error estimation, the error of real part of CRI is less than 0.046 and that of imaginary part is less than 0.003 in three typical modal (WS: water soluble, BB: biomass burning, DU: dust). As this algorithm applied to AERONET measurements, the total uncertainties are Δnf/c = 0.11, Δkf/c = 78% by considering all possible input of AERONET parameter errors together.

#### **3. Results**

#### *3.1. Fine- and Coarse-Mode AOD*

O'Neil et al. (2003) developed a spectral deconvolution algorithm (SDA) that utilizes spectral total extinction AOD data with the assumption of bimodal aerosol size distribution to infer the fine and coarse mode contributions to atmospheric AOD. SONET employs the algorithm's ability to separated coarse and fine mode AOD that used in Figure 2 [25,36]. As illustrated in Figure 2, we find that the fine-mode AOD (AODf) has the same monthly tendency with the total AOD. The annual AODf percent is 92%, which varies from 86% to 96%. Nevertheless, the coarse-mode AOD (AODc) has little variation in different month. It is demonstrable that the fine-mode aerosols are the principal pollutant in Chengdu area, which lead to the change of AOD. In summer, rainfall affected by the typical subtropical monsoon basin climate is significantly higher than that in winter [37]. The increasing precipitation can play a role in the removal of atmospheric pollutants. As a result, AOD decreases in summer, which reaches the lowest value of 0.64 in June. However, the average AOD in summer is 0.89, which still acts as a pollutant. There are two main reasons: firstly, the average surface wind speed all over the basin is low under the control of subtropical anticyclone in summer [38], which makes aerosol transport to other regions difficult. In addition, high temperature can cause the increasing formation of secondary organic aerosol particles [39]. In winter, the large variation of temperature from day to night leads to the rapid condensation of water vapor, which is conducive to the concentration of water and particulate matter. And the cold air is not conductive to the diffusion of aerosol because of the enclosed basin [16]. Therefore, AOD is the highest in winter (Winter: AOD = 1.12).

ˊ **Figure 2.** Monthly mean of aerosol optical depth (AOD), AODf, AODc at 440 nm and Angstrom Exponent (440–870 nm) in SONET Chengdu site from 2013 to 2016.

AE (440–870 nm) changes little with the value more than 1.0, which is similar to the previous research [12,20]. Convective precipitation occurs frequently in summer that the majority of coarse particles have been eliminated by wet settlement, but the fine particles remained in atmosphere with long suspension times [20,40–42]. Also, in Figure 2, it also can be illustrated by that the difference between AODf and AOD is lowest in summer. The average AE is less than 1.2 in each month of spring. The main reason for the high AOD and low AE in spring is the long-distance transport of dust pollution from North China [43,44]. In addition, another reason may be the wind speed increasing, by which it can be speculated that large particles are emitted to the atmosphere or transported from other regions by strong wind [20,45].

The AOD can be separated into fine- and coarse-mode and the different mode has the obviously distinguishing extinction for solar light in Chengdu. Therefore, it is important to explore the more sub-mode properties to research on the climate effect of different mode aerosols.

#### *3.2. Sub-Mode VSD*

VSD is one of the important aerosol microphysical properties. Although the particle size changes constantly, atmospheric aerosol particles exist in three-mode VSD stably. According to particle radius, the modes can be divided into nuclear mode (less than 0.1 μm), accumulation mode (0.1–1 μm) and coarse mode (larger than 1 μm) [46]. The process of nuclear mode is concentrated in nanometer scale, and the accumulation mode and coarse mode can be shown in VSDs [47,48]. According to Eck et al. (2012), large fine mode-dominated aerosols (submicron radius) were observed after the fog or low-altitude cloud dissipated events. As cloud condensation nuclei or ice nuclei, the smaller coarse mode-dominated aerosols (supermicron radius) are involved in the fog/cloud formation and dissipation [49]. Li, et al. (2014) found that an unusual increase of submicron fine modes is an important mechanism for haze growth in the polluted region. As identified in Figure 3, there are four modes of the aerosol VSD in Chengdu area and the corresponding parameters are listed in Table 1. The different color lines are the average of each mode and the values of N in the legend are the amounts of data involved in averaging. The modes can be clearly distinguished. As the median radius is less than 1μm, there are two peaks: the fine mode and submicron fine (SMF) mode. Furthermore, there are two peaks as the radius more than 1 μm: supermicron coarse (SMC) mode and coarse mode. In the figure, the fine mode has an almost equal amount of data as the coarse mode, but the fine-mode with one addition, which is one record with only one-peak VSD. The SMC mode has the least data amount of all (*N* = 133). The median radius of each mode is obviously different. As to the fine mode, the average median radius is 0.17 μm ± 0.05 and the volume modal concentration is 0.10 ± 0.07, that is, the largest

value of all modes. The median radii of the SMF modes are mainly in the range of 0.2–0.5 μm, which can indicate fog dissipation and haze growth, as previously mentioned. The SMC mode has a median radius of 1.62 μm ± 0.45 and the volume modal concentration of 0.07 ± 0.05. The coarse mode has an average median radius of 3.25 μm ± 0.99.

**Figure 3.** The average sub-mode volume size distribution in Chengdu.


**Table 1.** The parameters of sub-mode volume size distribution.

#### *3.3. Fine- and Coarse-Mode VSD and CRI*

The primary focus of Section 3.2 shows that the VSD can be separated into four distinguishing modes. However, in general, the natural aerosols, which are predominately coarse mode particles (r > 1 μm), and combustion-produced and anthropogenic emissions particles, which are predominately fine mode particles (r < 1 μm), of various mixed relative fractions are the mixtures in the aerosols. Furthermore, the fine mode and coarse mode particles are from different components. Also, at the same time, the CRI can only be separated into fine- or coarse-mode due to technical and precision limitations. Therefore, we focus on the fine- and coarse-mode VSD and CRI in Chengdu in this section.

In Figure 4, we present the multi-year average separated CRIs and the breakdown results of VSDs. It can be seen the fine- and coarse-mode VSD are well separated. The pictures (a) and (b) in the first row are the average values retrieved from the ground-based Sun-sky radiometer over multiple years. The second and third rows are the average of fine- and coarse-mode VSD and CRI. The sub-mode real part of CRI has non variation of wavelengths referring to Equation (2). The total CRI and sub-mode CRI are from different algorithms [21,32,34], but the sub-mode CRI is related to the volume modal concentration in the iterative algorithm of estimation of CRI for fine and coarse mode [34].

The fine-mode volume modal concentration is clearly higher than the coarse-mode one. The fineand coarse-mode real parts of CRI exhibit little difference (nf = 1.43, nc = 1.46). However, the coarse mode has a lower imaginary part of CRI at 440 nm than that of fine mode (kf440 = 0.0106, kc440 = 0.0072). At longer wavelengths, the imaginary part of CRI has little variation between fine- and coarse-mode (kf = 0.0121, kc = 0.0112).

**Figure 4.** The separated volume size distribution and sub-mode complex refractive index of multi-year average. The data source of (**a**) and (**b**) for total modes are from the inversion algorithm of SONET, and others (**c**–**f**) are from the sub-mode algorithms. (**a**) Total volume size distribution (**b**) Total complex refractive index in four wavelengths; (**c**) Fine-mode volume size distribution; (**d**) Fine-mode complex refractive index; (**e**) Coarse-mode volume size distribution; (**f**) Coarse-mode complex refractive index; The solid lines in figure (**b**), (**d**), (**f**) are the real parts of complex refractive index, and the dash lines are the imaginary parts. The corresponding parameters are listed in Table 2.

Figure 5 shows the fine-mode VSDs and CRIs in different seasons. The corresponding parameters are listed in Table 2. The typical bimodal or multimodal VSDs in all seasons imply a fine-coarse mixed–size distribution in the Chengdu area, similar to the urban-industrial aerosol type [50]. The fine-mode volume concentration is higher in summer, followed by winter, spring and autumn (Table 2). Furthermore, the fine-mode median radii are higher in summer (0.21) and winter (0.21). This can be explained by the increasing precipitation in summer and the low wind speed in winter that would lead to the high relative humidity and hygroscopic growth [22,23].

ˊ **Figure 5.** Seasonal variation of fine-mode volume size distribution (VSD) and complex refractive index (CRI) in SONET Chengdu site. The grey line refers to the total parameters. The corresponding parameters are listed in Table 2.

ˊ



The CRI can reflect the aerosol chemical composition. The real part indicates the aerosol refractivity and scattering characteristics. Specifically, the real part of CRI of water is generally considered to be 1.33 [51], much lower than other dry matters. Therefore, the real part can also reflect the water content in aerosol. The total real part of CRI in summer is the lowest, which is associated with the high humidity and water content in summer. Moreover, the fine-mode real part of CRI is also the lowest (nf = 1.38). As shown in Table 2, the fine-mode real parts of CRIs are all lower than the coarse-mode ones, which suggests that the water content in fine-mode particles play a leading role in Chengdu aerosols. ˊ

The imaginary part of CRI is related to the absorption characteristics of aerosol. The fine-mode imaginary part of CRI generally indicates that the fine-mode absorption component, that is, Black carbon (BC) or Brown carbon (BrC) [52].

Figure 6 presents the coarse-mode VSDs and CRIs in different seasons. The coarse-mode volume concentration is obviously less than the fine-mode. In particular, the volume concentration in summer gets the lowest (Vc = 0.069), which is associated with the wet removal of coarse particles. In summer, the coarse-mode real part is the lowest value over the four seasons, but also higher than the fine-mode one in Figure 5 that suggests the coarse particles are weakly hygroscopic (nf = 1.38, nc = 1.41). For all seasons, the coarse-mode imaginary part of CRI is quite constant (Table 2). In contrast, the fine-mode imaginary part has great seasonal variation. That demonstrates the coarse particle components are relatively stable. The spectral difference of the imaginary part (between 440 nm and longer wavelength) has a discrepancy between fine- and coarse-mode. It should be mentioned that the imaginary part reflects the aerosol absorbing property, and its spectral pattern can reveal the relative fractions of absorbing aerosols, that is, BC and DU [52–54].

**Figure 6.** ˊSeasonal variation of coarse-mode VSD and CRI in SONET Chengdu site. The grey line refers to the total parameters. The corresponding parameters are listed in Table 2.

#### **4. Discussion**

In order to evaluate the overall analysis scheme, a numerical experiment was used to assess the performance of sub-mode results. In Figure 7, we illustrated the recovery of AOD in four wavelengths by the sub-mode VSD and CRI. It can be seen that the correlation coefficients are all larger than 0.98. The absolute deviations in four wavelengths are 0.04, 0.02, 0.03 and 0.03, corresponding relative standard deviations are 0.03, 0.03, 0.06, 0.09, respectively. These biases are basically close to the claimed uncertainties of SONET products (AOD), which demonstrate that the sub-mode VSD and CRI results are acceptable in understanding of optical closure.

**Figure 7.** Recovery of AOD at different wavelength based on the separated volume size distribution and sub-mode complex refractive index. *N* = 380 and curves (at 675 nm, 870 nm, 1020 nm) are shifted for a better viewing.

For further details of the climate effect, the sub-mode SSAs were taken into account. SSAs could reflect the aerosol absorption and scattering of solar lights, which is an important parameter in climate modeling [55,56]. Hansen et al. (1997) noted that a change in SSA from 0.9 to 0.8 can change the radiative forcing from negative to positive depending on the reflectance of the underlying surface and the altitude of the aerosols. Moreover, strongly absorbing aerosols may have a large impact on the regional climate and heating the atmosphere [57,58].

In Figure 8, the sub-mode SSA is calculated by the separated VSD and sub-mode CRI under ignoring the influence of nonsphericity on dust aerosols. The total SSAs are all larger than 0.9 in the four seasons and the fine-mode SSAs are closed to the total SSAs that indicates the scattering properties are mainly dominated by fine-mode aerosol. Higher SSA at 675 nm could indicate the main absorption component BrC [52,53], but the high value of SSA indicates the absorbing effect is obvious less than scattering effect of solar lights. The SSA spectral trends from 675 nm to 870 nm show different patterns: the increasing pattern is usually corresponding to the absorbing coarse particles (e.g., Dust); the decreasing pattern is usually corresponding to the absorbing fine particles (e.g., BC, BrC) [52]. The fine-mode SSAs decrease and the coarse-mode SSAs increase with the increase of wavelength, which indicates the different absorbing component in different modes. Furthermore, the coarse-mode SSAs are lower than that of the fine-mode that suggests the absorption of coarse-mode particles is likely stronger than that of the small size aerosol particles. In this regard, the absorbing component in large-size aerosols (such as mineral dust, biomass burning, etc.) is possibly more than that of small-size aerosols.

**Figure 8.** The sub-mode single scattering albedo in different season based on the separated volume size distribution and sub-mode complex refractive index (ignore the influence of nonsphericity on dust aerosols). The grey, red, blue lines are the whole-, fine- and coarse-mode single scattering albedos (SSAs), respectively.

#### **5. Conclusions**

In this study, we investigated the sub-mode VSD and CRI for fine- or coarse-mode in the Chengdu area. The annual average of AODf percentage is 92%, which has the same monthly tendency with the total AOD, but AODc has little variation in different months. The typical bimodal or multimodal VSDs employ a fine-coarse mixed–size distribution in Chengdu area. There are four distinguishing modes of VSD in Chengdu that the median radii are 0.17 μm ± 0.05, 0.31 μm ± 0.12, 1.62 μm ± 0.45, 3.25 μm ± 0.99, respectively.

The fine-mode annual average volume modal concentration is clearly higher than the coarse-mode one. The fine-mode volume concentration and median radius are higher in summer and winter. In particular, the coarse-mode volume concentration gets lowest in summer (Vc = 0.069), which is associated with the wet removal of large-size particles. For multi-year average results of CRI, the fine-mode and coarse-mode real parts show little difference. However, the coarse mode has a lower imaginary part at 440 nm than the fine-mode (kf440 = 0.0106, kc440 = 0.0072). At longer wavelengths, the imaginary part of CRI has little variation between fine and coarse mode (kf = 0.0121, kc = 0.0112). In summer, both fine- and coarse-mode real parts get the lowest value respectively because of the high humidity (nf = 1.38, nc = 1.41). For all seasons, the coarse-mode imaginary part of CRI is quite constant, but the fine-mode imaginary part has great seasonal variations. It indicates the coarse particle components are relatively stable.

In order to assess the performance of the sub-mode results, we illustrated the recovery of AOD by the sub-mode VSD and CRI. It can be seen that all the correlation coefficients are larger than 0.98. The sub-mode SSAs are calculated by sub-mode VSD and CRI under the condition of neglecting the non-sphericity. The total SSAs are all larger than 0.9 in the four seasons and the fine-mode SSAs are closed to the total SSAs that indicates the scattering properties are mainly dominated by fine-mode aerosols. The coarse-mode SSAs are lower than fine-mode, which suggests the absorbing effect of coarse-mode particles is likely stronger than that of the small size aerosol particles.

**Author Contributions:** Z.L. conceived and designed this study, and participated in drafting and revising the article. C.Z. substantially contributed to the analysis and interpretation of data, and drafted the articles. H.X. reviewed and edited the manuscript. Y.Z. and Y.W. undertook a part of the instrument maintenance job. D.L., K.L. and Y.X. contributed to the calibration of the sun-sky radiometer and analyses of the observation data. Y.Z. do formal analysis and investigation.

**Acknowledgments:** This work was supported by the National Key Research and Development Program of China under Grant 2016YFE0201400, the National Natural Science Foundation of China under Grant 41671364,41671367, 41771396, and Science and Technology Service network initiative (STS) Project of Chinese Academy of Sciences (KFJ-STS-QYZD-022).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Summertime Urban Mixing Layer Height over Sofia, Bulgaria**

#### **Ventsislav Danchovski**

Department of Meteorology and Geophysics, Faculty of Physics, Sofia University, Sofia 1164, Bulgaria; danchovski@phys.uni-sofia.bg; Tel.: +359-2-8161-413

Received: 7 November 2018; Accepted: 9 January 2019; Published: 17 January 2019

**Abstract:** Mixing layer height (MLH) is a crucial parameter for air quality modelling that is still not routinely measured. Common methods for MLH determination use atmospheric profiles recorded by radiosonde but this process suffers from coarse temporal resolution since the balloon is usually launched only twice a day. Recently, cheap ceilometers are gaining popularity in the retrieval of MLH diurnal evolution based on aerosol profiles. This study presents a comparison between proprietary (Jenoptik) and freely available (STRAT) algorithms to retrieve MLH diurnal cycle over an urban area. The comparison was conducted in the summer season when MLH is above the full overlapping height of the ceilometer in order to minimize negative impact of the biaxial LiDAR's drawback. Moreover, fogs or very low clouds which can deteriorate the ceilometer retrieval accuracy are very unlikely to be present in summer. The MLHs determined from the ceilometer were verified against those measured from the radiosonde, which were estimated using the parcel, lapse rate, and Richardson methods (the Richardson method was used as a reference in this study). We found that the STRAT and Jenoptik methods gave lower MLH values than radiosonde with an underestimation of about 150 m and 650 m, respectively. Additionally, STRAT showed some potential in tracking the MLH diurnal evolution, especially during the day. A daily MLH maximum of about 2000 m was found in the late afternoon (18–19 LT). The Jenoptik algorithm showed comparable results to the STRAT algorithm during the night (although both methods sometimes misleadingly reported residual or advected layers as the mixing layer (ML)). During the morning transition the Jenoptik algorithm outperformed STRAT, which suffers from abrupt changes in MLH due to integrated layer attribution. However, daytime performance of Jenoptik was worse, especially in the afternoon when the algorithm often cannot estimate any MLH (in the period 13–16 LT the method reports MLHs in only 15–30% of all cases). This makes day-to-day tracing of MLH diurnal evolution virtually impracticable. This problem is possibly due to its early version (JO-CloVis 8.80, 2009) and issues with real-time processing of a single profile combined with the low signal-to-noise ratio of the ceilometer. Both LiDAR-based algorithms have trouble in the evening transition since they rely on aerosol signature which is more affected by the mixing processes in the past hours than the current turbulent mixing.

**Keywords:** mixing layer; urban area; ceilometer; radiosonde

**PACS:** 01.30.-y; 01.30.Ww; 01.30.Xx

#### **1. Introduction**

The effect of air quality on human health is a serious problem, especially in densely populated areas. Hence, a lot of effort is being made to better understand the processes controlling pollution levels, particularly in numerical modelling. Key input parameters of these models are meteorological variables, which are needed to be identified in order to calculate the production, diffusion, transport and scavenging of atmospheric pollutants. These harmful substances are dispersed vertically within the mixing layer (ML) due to its inherent turbulence. According to Seibert et al. [1], ML is ". . . the layer adjacent to the ground over which pollutants or any constituents emitted within this layer or entrained into it become vertically dispersed by convection or mechanical turbulence within a time scale of about an hour". However, one should bear in mind that there are situations when time-scales of the dominant processes (such as diabatic processes like radiative cooling in the evening transition, or unsteadiness of pressure gradients, or intermittent turbulence due to breaking gravity waves, just to name a few) are much longer so that the ML is unsteady [2]. Obviously, near-ground pollution levels will depend on the mixing layer height (MLH) since it constrains the dispersion volume. Thus, the MLH is vitally important to be identified especially in urban areas where pollution sources and inhabitants are much greater [3–10]. Moreover, urban MLH can be characterized by enormous temporal and spatial variability due to inhomogeneity in surface roughness and heating in cities [11]. Therefore, MLH is worthwhile to be continuously monitored and also compared with parametrizations in numerical weather and/or pollution prediction models [12–14].

Despite its importance, MLH is not a part of routine measurements. Furthermore, because it is associated with the spatial distribution of turbulence, we need turbulence profiles to determine MLH. Consequently, TKE (turbulent kinetic energy)-based criteria (MLH is marked by the level where TKE drops below a predefined threshold) are often used in numerical models with turbulence closure of order 1.5 or higher for MLH determination [1,15]. Moreover, profiles of the TKE and its dissipation rate can be measured by remote sensing instruments [16,17]. Therefore, Doppler LiDARs [18] and sodars [19,20] can serve as "turbulence profilers" but the former are quite expensive and the latter have limited vertical range. Fortunately, vertical profiles of non-reactive scalar meteorological variables should be nearly constant with height within a well-mixed boundary layer [21], so we can detect the MLH by looking for abrupt changes in the uniformly distributed profiles of these tracers [22,23].

Regardless of the wide variety of remote sensing methods, the most used instrument for MLH detection is the radiosonde, which is still used as a reference. The derivation of the MLH from the radiosonde profiles of the atmospheric temperature, humidity, and wind dates back to 1960s [24]. Moreover, these radiosonde-based methods are still used independently [25] or as a reference for validating MLH measurements from remote sensing instruments [26]. However, the radiosonde also has some drawbacks as it measures atmospheric properties along its flight, which is slant instead of vertical, due to horizontal wind. Therefore, the radiosonde profiles do not coincide with rising thermal or vertical profiles derived from the remote sensing instruments. Additionally, the radiosondes' main limitation is their coarse temporal resolution since they are usually launched no more than twice a day.

The necessity of continuous MLH monitoring can be met by operating ground-based remote sensing instruments. A comprehensive review of existing techniques for MLH determination through ground-based remote sensing instruments, along with their advantages and limitations, can be found in Wiegner et al. [27] and Emeis et al. [28]. It is worthwhile to note that individual disadvantages of each instrument in the MLH diurnal cycle determination can be overcome if apparatuses are used together [29].

One should note that relying on ground-based remote sensing instruments in MLH estimation cannot provide good spatial representativeness, especially over areas with non-homogeneous land-use and/or complex topography. Fortunately, this lack of information can be filled if space-based remote sensing instruments are used. Among them are the Cloud-Aerosol LiDAR with Orthogonal Polarization (CALIOP) [30] and the Moderate Resolution Imaging Spectroradiometer (MODIS) [31], which are the most used for determination of the atmospheric boundary layer height over continents and oceans [32–35]. The radio occultation method based on global position system signals can provide vertical profiles of the refractivity index that can be used in MLH retrieval [36,37].

In recent years, laser-based remote sensing instruments, especially automatic LiDARs and ceilometers (ALC), have become more affordable and widely used in the field of atmospheric research, particularly for MLH determination [38–45]. We should also note the considerable efforts made in the COST Action ES 1303 TOPROF [46] which provides standards for calibrated profiles of the aerosols, winds, temperature, and humidity to fill the observational gap in the lower troposphere. These quality controlled observations are delivered in near real-time through the EUMETNET Composite Observing System (EUCOS) network E-PROFILE [47] to the national weather services in order to improve numerical weather prediction.

Different LiDAR-based methods for MLH retrieval from the range corrected signal are summarized in Haeffelin et al. [48]. These methods are the basis of many proprietary [49,50] and in-house [51–54] algorithms explicitly designed for MLH retrieval from ceilometers' data. Possible performance improvement of these LiDAR-based techniques can be achieved by monitoring diurnal variations in Radon-222 [55] or it can be used alone to evaluate MLH [56].

The objectives of the present study are to evaluate the performance of a proprietary algorithm as well as a popular, freely available algorithm in detecting of MLH from ceilometer data over an urban area. Both methods are evaluated against MLHs retrieved using radiosonde profiles as a reference. The structure and evolution of the mixing layer over Sofia in summertime is also discussed. To highlight the advantages and disadvantages of both algorithms, an analysis was performed in the summer when MLH is high enough to minimize the negative effects due to incomplete overlapping in the near-field range of the ceilometer. The paper is organized as follow: measurement sites and specifications of the instruments used for observations, as well as details about the collected data and methods applied to determine the MLH, are in Section 2. Inter-comparison of the three radiosonde-based methods is in Section 3.1. Verification of the MLH retrieved from the ceilometer compared to the MLH measured from the radiosonde is in Section 3.2. Diurnal evolution of the MLH over Sofia derived from the ceilometer data by both the proprietary and the freely available algorithm, as well as a discussion of their main benefits and drawbacks and suggestions for performance improvement, are in Section 3.3. Statistical analysis of the MLH diurnal cycle is discussed in Section 3.4. The article ends with summary of findings.

#### **2. Data and Methodology**

Sofia is the largest and the most densely populated city in Bulgaria with roughly 1,400,000 inhabitants. The city is located in a valley that is almost fully encircled by mountains; therefore, the micro- and meso-scale processes, as well as the ML dynamics, are heavily influenced by both complex orography and urban territory. To perform our analysis of urban MLH we used 3 months of intensive measurements, from 1 June until 31 August 2015. The data used in this study were obtained from a continuously operating ceilometer, Jenoptik CHM 15k (in 2014, the company G. Lufft Mess- und Regeltechnik GmbH acquired the product segment of ceilometers from Jenoptik and now the ceilometer is known as Lufft CHM 15k), and a balloon sounding launched on a daily basis at 12:00 UTC (14:00 LT).

The CHM15k (firmware version 0.63) is operated by the Department of Meteorology and Geophysics, Sofia University. The ceilometer is situated in the city centre on the territory of the University Astronomical Observatory in the park "Borisova gradina" (Figure 1). The CHM15k is an eye-safe biaxial LiDAR system equipped with an Nd: YAG solid-state near-infra-red laser operating at 1064 nm. It emits pulses with an energy of 8 μJ and a repetition frequency of 5–7 kHz. The ceilometer provides data with a vertical resolution of 15 m, the maximum height of the signal is 15,000 m and the temporal resolution is set up to 60 s. As CHM15k is a biaxial LiDAR it suffers from incomplete overlap in the near range since only a small portion of the laser beam gets into the receiver field of view. According to the manufacturer, the overlapping is ∼1% at 15th bin (225 m), ∼10% at 24th bin (360 m), ∼50% at 40th bin (600 m), ∼90% at 57th bin (855 m), ∼99% at 78th bin (1170 m) and full overlap is achieved at 120th bin (1800 m). Further details about the instrument can be found in [57].

The MLH was retrieved from the ceilometer profiles by supposing that aerosol concentration is rapidly adapted to the thermal stratification of the ML and that aerosol loading above the city is not dominated by advection. The manufacturer's software includes a proprietary algorithm for automatically deriving the MLH every minute (software JO-CloVis version 8.80) [58]. Because the Jenoptik algorithm is proprietary not much is known about it. Haeffelin et al. [48] reported that the

algorithm uses vertical derivatives and wavelet transforms on the range-corrected signal to identify local minima which are used as MLH, however, it is not specified which version is referred to and there is no information on how signal-to-noise ratio is enhanced in the pre-processing. A freely available Structure of the Atmosphere (STRAT) algorithm [52], which is designed for the retrieval of aerosol vertical profiles in the atmospheric boundary layer and free troposphere, is used for comparison. In contrast to the Jenoptik algorithm that is based on vertical gradients of backscatter signal in a single profile, STRAT uses both temporal and vertical gradients (*Gt* and *Gv*) by using Sobel 2-D derivation operators. The global gradient is calculated as *G* = *Gt* <sup>2</sup> + *Gv* 2. The edges in backscatter are kept if *G* is greater than predefined thresholds. Additionally, edges in low signal-to-noise ratio zones are rejected. Finally, minimum (450 m due to overlapping) and maximum (here, 3000 m during the day and 1500 m at night are used) allowed heights are applied and three global gradients—the strongest, the second strongest, and the lowest-height—are reported as MLH candidates [48]. The MLH is then determined as the lowest-height candidate at night, during day quality control based on relative change in backscatter around each candidate is performed and the first existing one in the line strongest, second strongest and lowest is selected as MLH [18].

**Figure 1.** The locations of the ceilometer and the radiosonde indicated by a blue and a purple triangle in the Sofia valley. (Source of the map is Google LLC).

The atmospheric sounding system is the Vaisala's MW41 which is located in the Central Aerological Observatory on the territory of the National Institute of Meteorology and Hydrology, which is about 4.4 km south-east from the ceilometer (Figure 1). In this study, low resolution radiosonde data are used, which are freely available at Integrated Global Radiosonde Archive (IGRA) [59]. Archived radiosonde data consist of atmospheric parameters recorded at mandatory and significant levels which were used to restore the atmospheric profiles by linear interpolation.

Following de Haij et al. [60], three different MLH detecting algorithms were applied to the radiosonde data. The Bulk Richardson (Ri) method is based on the Richardson number which is the ratio of thermally and mechanically driven turbulence. According to this method, MLH is the level where the bulk Richardson number exceeds predefined threshold values [61–63]. In this study, the commonly recommended value of 0.21 was used. It is worth noting that the Ri method is suitable for both convective and stable conditions. In the parcel method [24,64], the MLH is determined by extending dry-adiabatically surface temperature to its intersection point with temperature profile. However, this method provides reliable results only for unstable convective boundary layer as it neglects wind shear effects on vertical mixing. The last method for MLH determination from radiosonde data is the lapse rate method [21,65]. It is based on threshold values of vertical gradients of potential temperature (*θ*) and relative humidity (RH). Adhering to de Haij et al. [60], negative gradients of RH and a gradient of *θ* >2 K/km were used as the basis for this study. As the selected critical value of potential temperature gradient is more or less subjectively chosen, the performance of lapse rate values of 0.5, 1, 1.5, 2.5, 3, 3.5, and 4 K/km was also tested.

As main synoptic-scale systems are associated with the suppression or stimulation of parcel ascending, it is interesting to examine their role on mixing layer height [66]. Therefore, the difference (Δ*p*) of surface layer atmospheric pressure (*p*) and its smoothed value (*psmooth*, which is obtained by low pass filter with cut-off 6 days) is calculated by Equation (1). Then Δ*p* is standardized by Equation (2), i.e., the Δ*p* is rescaled to have a mean of zero (subtraction of the mean value Δ*p*) and a standard deviation of one (division by the standard deviation *σ*Δ*p*).

$$
\Delta p = p - p\_{smooth} \tag{1}
$$

$$
\Delta p\_{std} = \frac{\Delta p - \overline{\Delta p}}{\sigma\_{\Delta p}} \tag{2}
$$

Finally, the atmospheric pressure (atm.press) is classified as "Low" if Δ*pstd* is smaller than −0.5 while it is marked as "High" if Δ*pstd* is higher than 0.5. If the Δ*pstd* values are greater than −0.5 but less than 0.5, atmospheric pressure is marked as "Normal".

#### **3. Results and Discussions**

#### *3.1. Inter-Comparison of Radiosonde-Based MLH Retrieval Methods*

The three aforementioned radiosonde-based algorithms—Richardson, parcel and lapse-rate were applied on the dataset for a total of 92 days (for 28 days the atmospheric pressure was "High", for 43 "Normal", and for 21 "Low"). MLH values were successfully estimated at 92, 92, and 81 days, respectively. The estimated MLHs were then compared against one another on Figure 2. The perfect correlation between the Richardson and parcel method indicates that in summer at 14:00 LT (12:00 UTC) the urban mixing layer over Sofia is dominated by thermally driven turbulence. It is a fairly expected result since the study period took place in summer and radiosonde launching occurred in the early afternoon. The box-plot shows that slightly higher MLH values are related to prevailing low atmospheric pressure and that when atmospheric pressure is marked as normal or high, MLHs are slightly decreased; however, the observed difference is not statistically significant (Wilcoxon–Mann–Whitney test with a significance level of 5% was performed).

The lapse rate method shows worse alignment with the Richardson and parcel methods, therefore, we tested how a threshold value of vertical gradient *θ* influences concurrence with the other two approaches. We found that lapse rate values of 1, 1.5 and 2 K/km perform similarly and Pearson correlation coefficients with respect to the Ri method are about 0.89. However, the correlation diminishes if smaller or higher threshold values are used. It is also worth mentioning that a negative vertical gradient of the relative humidity is not changed because it agrees with the mixing layer conception (the Earth's surface is the water vapour source and free atmosphere is low in humidity, so the humidity gradient should be negative at MLH). Keeping in mind that the Ri method incorporates both mechanical and buoyancy production of turbulent mixing we choose it as a reference in the following analysis. Performance of the parcel and lapse rate methods using the set of critical values mentioned above were evaluated against the Richardson method and is summarised in Table 1.

**Figure 2.** Inter-comparison of the three radiosonde-based MLH methods. The correlation matrix (**a**) shows correlation coefficients in the upper-right triangle, the diagonal shows a histogram of each method, and the lower-left triangle shows scatter-plots and linear regression lines with corresponding 95% confidence intervals. The box and whisker plot (in the style of Tukey) is on plot (**b**). The box lines correspond to the 25, 50 and 75 percentiles. The lower and upper whiskers represent the lowest values still within 1.5 IQR (inter-quantile range) of the lower quartile, and the highest values still within 1.5 IQR of the upper quartile. The data beyond the end of the whiskers signify outliers and are plotted as black dots. White dots indicate mean values. In both figures, atmospheric pressure is color-coded as "High" (blue), "Low" (yellow), "Normal" (green).


**Table 1.** Skill scores (MD—mean deviation; RMSD—root-mean-square deviation, r—Pearson correlation coefficient, slope—linear regression slope, intercept—linear regression intercept) of parcel and lapse rate methods compared against the Richardson method as a reference in the MLH determination.

#### *3.2. Inter-Comparison of MLHs Derived from Ceilometer and Radiosonde Data*

MLHs calculated from radiosonde data are often used for reference since they are based on the thermodynamic structure of the lowest atmosphere that directly reflects changes in the surface forcing. However, since routine balloon launching usually occurs only twice a day so it does not allow for MLH diurnal evolution to be tracked. Low-cost ceilometers that provide backscatter power profiles are a tempting alternative because they operate continuously.

To evaluate the overall performance of the ceilometer-based methods in the MLH determination, the calculated values are compared against the Richardson method estimates from the radiosonde data. Since the radiosonde in Sofia is launched once a day at 12:00 UTC (14:00 LT), the ceilometer-retrieved MLHs from within a 20 min timespan are averaged and used in the comparison. After this procedure the size of the STRAT's datasets at "High", "Normal" and "Low" atmospheric pressure is reduced to 18 (64%), 24 (56%) and 10 (48%) days, respectively. The Jenoptik algorithm successfully estimates

MLHs in 13 (46%), 17 (40%) and 11 (52%) days at "High", "Normal" and "Low" atmospheric pressure, respectively. In other words, both ceilometer-based algorithms cannot estimate MLHs in about half of the days with "Low" atmospheric pressure. The percentage of the Jenoptik-retrieved MLHs becomes even lower at "Normal" and "High" pressure, while the performance of STRAT is slightly increased. The left and right panels of Figure 3 show a correlation matrix and box and whiskers plots of the MLH determined by STRAT, Jenoptik and Ri methods at different atmospheric pressures. It is evident that both LiDAR-based algorithms tend to underestimate MLH compared to radiosonde (Richardson). We should bear in mind that MLH estimation from the ceilometer and the radiosonde data rely on different tracers, which may contribute to the observed discrepancy. When optically thick clouds or rain are presented the backscatter signal can be strong enough to saturate the ceilometers receiver so the cloud base or somewhere under the cloud within the rain column is reported as the MLH. To prove the hypothesis, the data was spited to rainy (if nonzero ceilometer's precipitation index is registered from 11:20 to 11:40 UTC) and dry cases. The analysis showed that the difference between Ri and STRAT, and between Ri and Jenoptik are statistically non-significant (t-test with a significant level of 0.05 is performed) in rainy days. In the rest of the days the Ri method reports about 180 m (750 m) higher MLH than STRAT (Jenoptik) and the results are statistically significant. Additionally the role of low clouds was tested. In days with low clouds (if cloud base height is < 1500 m) the difference between Ri and STRAT and between Ri and Jenoptik are evaluated as a statistically non-significant. In the rest of the days radiosonde estimates are 212 m and 493 m higher than STRAT and Jenoptik respectively and t-test showed that both results are statistically significant. The observed underestimation of the MLH by the ceilometer could be attributed to the difference in land surface type. The ceilometer is situated in the park (while the radiosonde is in a built-up area) so one can expect that some of the solar energy is consumed during evapotranspiration; therefore the rest of the energy that would produce the thermally driven turbulence, and thus MLH raising, is reduced. To prove this, the number of consecutive days with no precipitation (ceilometer's precipitation index is used for the classification) are used to split the data into categories. It was found that the difference between Ri and STRAT methods (the Jenoptik is not included since it shows significant deviation from the Richardson method, see Figure 3) increases with the number of consecutive droughty days, which is supposed to be a result of the lack of available water for evaporation in a built-up zone. Additionally, if the drought period is long, MLH becomes higher and the ceilometer-based method experiences difficulties that are supposed to be a result of diminished backscatter due to the increased volume for aerosol dispersion (Figure 4).

**Figure 3.** A correlation matrix (**a**) and Tukey's box and whiskers plot (**b**) of radiosonde-(Richardson) and LiDAR-based (STRAT and Jenoptik) algorithms for MLH detection. Conventions are the same as in Figure 2.

**Figure 4.** The dependence of drought duration (in number of dry days) on the mean MLH determined by Richardson and STRAT methods.

Skill scores of LiDAR-based algorithms against the Ri method are listed in Table 2. It can be seen that the average underestimation of MLH by STRAT and Jenoptik is ∼160 m and ∼660 m respectively. The STRAT-estimated values of MLH are reasonably comparable with those retrieved from radiosonde profiles, but Jenoptik's performance is quite unpromising and needs further clarification.

**Table 2.** Skill scores (MD—mean deviation; RMSD—root-mean-square deviation, r—Pearson correlation coefficient, slope—linear regression slope, intercept—linear regression intercept ) of aerosol-based algorithms (Jenoptik and STRAT) compared against the Richardson method as a reference in the MLH determination.


#### *3.3. Diurnal Evolution of the MLH Determined by the Ceilometer—A Case Study*

To elucidate the above-mentioned ceilometer's capacity to track the MLH diurnal cycle, a case study is first considered. In Figure 5 diurnal evolution of the range-corrected ceilometer signal (PR2) on July 24 is presented along with MLHs determined according to STRAT and Jenoptik algorithms. Radiosonde-derived MLH by the Richardson method is also plotted for comparison.

**Figure 5.** Time-height cross section of the ceilometer's range-corrected backscatter power (PR<sup>2</sup> in arbitrary units) on 24 July 2015. The MLH retrieved from ceilometer's data by Jenoptik and STRAT algorithms are marked by magenta triangles and red circles, respectively (for clarity, the Jenoptik MLHs are plotted with the same temporal resolution as STRAT—10 min). Radiosonde-based MLH according to the Ri method is presented by black "x" marks.

As observed, the range corrected ceilometer's signal reveals some characteristic features in the MLH diurnal evolution. The backscatter power within the first 500–700 m is high in the early hours of the night which can be associated with mechanically mixed aerosols within the nocturnal boundary layer. As seen, the layer was shrinking and at ∼8:00 LT (less than 2 h after sunrise) it had disappeared. One may expect a new convective layer to be identifiable at that moment but we should keep in mind that the ceilometer has virtually zero overlapping in the first ∼200 m (overlapping is <1%) so that the first signs of the rising thermals are visible at ∼9:00 LT. Above the nocturnal layer, there is a zone with decreased signal that is capped by a high backscatter layer, which most likely outlines aerosol burden air in the residual layer, or it is a result of advection at that elevation. The ceilometer's signal also depicts the daytime evolution of the MLH with its typical growth due to the solar heating of the surface. After sunrise, thermals start forming and rising due to positive buoyancy. These updrafts produce turbulent mixing so that the diminished vertical backscatter within ML in the afternoon results from an increased volume for aerosol dispersion. An enhanced signal close to the ML top in the afternoon can be attributed to hygroscopic growth of aerosols due to increased relative humidity. As can be seen, MLH reached its maximum (∼2250 m) at ∼16:00 LT and remains almost constant until ∼19:00 LT. In the evening, the thermals cease to form (in the absence of cold air advection), allowing turbulence to decay in the formerly well mixed layer. A new nocturnal layer starts forming and overhead air associated with the new residual layer becomes decoupled from the mechanical source of turbulence on the ground. However, the evening transition period is non-stationary as heat fluxes decrease over a few hours after sunset so that the aerosol vertical distribution does not respond to surface forcing within an hour [2]. At that part of the day the ceilometer profiles are mostly a result of the turbulence dynamics in the recent periods, therefore, they do not reveal the present ML but its history. That problem is inherent to all remote sensing instruments which use aerosol backscatter to trace the ML but can be overcome if a "true turbulence profiler" is used. As observed, the MLH determined by the Ri method is ∼1940 m which corresponds very well to the aerosol distribution depicted by the ceilometer backscatter signal at the moment of balloon launching (at 13:30 LT that day).

It is also noticeable that the STRAT algorithm plausibly represents the diurnal evolution of the MLH. However, in the time interval from sunrise (6:09 LT) to approximately 10:30 LT, which corresponds to the morning transition period, STRAT misleadingly reports an overhead backscatter gradient (associated with the residual layer) as MLH instead of the one closest to the ground. Similar behaviour is found across all days and seems to be due to the layer attribution technique implemented in STRAT. According to Haeffelin et al. [48] the algorithm reports the strongest, second strongest, and the lowest gradients in backscatter and then, depending on the local time, it constructs a diurnal evolution of its "best estimate" (used here as MLH) which is the lowest gradient during the night and the strongest gradient during the day. Thus, the STRAT method reports abrupt changes in MLH around sunrise and sunset instead of smooth transitions from the nocturnal to convective boundary layer, and vice versa. Possible improvement of layer attributions and representations of the MLH diurnal evolution can be achieved through the use of statistical analysis [67] or graph theory [68].

It can be seen that the overall consistency of the MLHs reported by the Jenoptik algorithm with the observed aerosol distribution and evolution is relatively poorer than the consistency of the STRAT's MLHs. However, Jenoptik outperforms STRAT in the morning transition, although neither method can track the MLH from 8 to 9 LT when the MLH is in the zone of incomplete overlapping. The performance of the Jenoptik method during daytime is much worse and it cannot represent the MLH evolution. The method cannot report MLH from 13 to 17 LT and it significantly underestimates ML depth around noon and in the late afternoon. It should be noted that STRAT also locates these aerosol gradients at intermediate levels (Figure 6) but reports them as the lowest and the second strongest candidates, which are then successfully filtered out by the attribution procedure in the algorithm. This worsened performance of the Jenoptik method is likely to be a result of the immaturity of the outdated version of the algorithm used. Additionally, the Jenoptik method operates in real-time so it is likely to use only the current backscatter profile without taking into account previous measurements. Therefore,

the signal-to-noise ratio (SNR) will be lower, which may result in poor performance compared to STRAT. Consequently, the poorer daytime performance of the Jenoptik can be attributed to reduced backscatter signal within the increased depth of the MLH (and enlarged volume for aerosol dispersion) and augmented background signal due to the higher sun elevation angle. Data shows the STRAT method also has similar troubles with backscatter gradient detection from ∼13:00 LT to ∼15:00 LT when only a few MLHs are reported. However, the process of smoothing incorporated within the algorithm enhances the SNR, enabling the MLH evolution to be tracked against the Jenoptik algorithm.

#### *3.4. Diurnal Evolution of the MLH Determined by Ceilometer—A Statistical Analysis*

To compare the performance of both LiDAR-based algorithms we first make the datasets comparable. Since the Jenoptik algorithm has a 1-min resolution but STRAT's temporal resolution is 10 min, Jenoptik-derived MLHs are averaged in 10-minute intervals. The availability of STRATand Jenoptik-derived MLHs after applying the described procedure is presented in Figure 7. As seen, MLH data availability of both methods show similar patterns related to the diurnal cycle. The STRAT algorithm manages to estimate MLHs in about 70–95% of the cases but in the afternoon its availability drops to 50–70% with minimum of ∼45% at 14 LT. In contrast, the Jenoptik method provides MLHs in about 60–85% of the cases but in the afternoon it hardly reaches even 35–40% with a minimum of ∼15% at 15 LT. The observed diurnal pattern in MLH availability in both LiDAR-based methods is closely related to decreased SNR due to reduced aerosol concentration (due to increased volume for aerosol dispersion) and increased background signal (due to higher solar radiation) in the afternoon. Neither of the two applied algorithms show a clear atmospheric pressure dependency.

**Figure 7.** Diurnal evolution of the availability of MLH determined by STRAT (**a**) and Jenoptik (**b**) algorithms at "Normal", "Low" and "High" atmospheric pressure in summer of 2015.

The aforementioned features of both LiDAR-based techniques are also visible if all data for the MLH daytime progress are summarised and presented as box-plots (Figure 8). From midnight to

7 LT STRAT and Jenoptik algorithms provide comparable MLHs and most of the estimated values are in the range of 500–1000 m. However, the Jenoptik also shows several quite large values marked as outliers (most of them in the range 1500–3500 m) which are a result of improper selection of high aerosol layers that cannot be related to the near-ground turbulence. When atmospheric pressure is "Normal" the Jenoptik algorithm also reports a few quite low MLH in the ceilometer's incomplete overlapping zone which should be treated as incorrect values (most likely they are result of multiple scattering). The morning transition is marked by STRAT as an abrupt jump that is a result of its layer attribution criterion, while the Jenoptik represents the transition less steeply. Daytime performance of both algorithms is, thus, easily distinguishable. The MLHs retrieved by Jenoptik are often in the first 1 km and rarely reach 2 km. As was noted, the algorithm tends to report mid-level gradients that are also marked by STRAT lowest-height and/or second strongest gradients. However, in STRAT these mid-level gradients in the ML are successfully filtered out by the successive layer attribution. The daily maximum of MLH (sometimes more than 2000 m) is registered in the late afternoon (∼18:00 LT), a few hours before sunset and, more importantly, during peak car traffic, which can help against excessive concentrations of air pollution. The evening transition is hard to be correctly traced by the LiDAR's backscatter profile as the aerosol signature is more related to turbulent mixing in the past than the current state. Therefore, although showing different behaviours, it is difficult to designate one of the two methods as more reliable. It is worth noting that there are a large number of outliers in the retrieved MLHs by the Jenoptik algorithm; most of them are related to high aerosol layers due to advection or residual layers at night. As seen, both methods report lower daytime MLHs in "Low" atmospheric pressure, especially Jenoptik algorithm, whose estimates do not reach 1 km in 50% of cases.

**Figure 8.** Diurnal cycle of the MLH over Sofia determined by STRAT (red) and Jenoptik (magenta) algorithms as a box and whiskers plot (in Tukey's style) at "High", "Low" and "Normal" atmospheric pressure in summer of 2015.

#### **4. Conclusions**

In this paper, MLHs derived by different algorithms over three summer months from radiosonde and ceilometer data were analysed and compared. It was shown that the Richardson and parcel methods produce identical MLHs which indicates that ML is primarily thermally driven while the lapse-rate method underestimates the MLH. It was also found that threshold values for potential temperature higher than 2 K/km or smaller than 1K/km deteriorate the agreement of lapse rate with the Richardson and parcel methods. Based on performed comparison against the Richardson method, it was shown that the ceilometer tends to underestimate the MLH. The observed discrepancy was mainly attributed to different land surface types where the instruments are situated and the distance between the sites. Additionally, the proprietary algorithm has difficulty with the low SNR of the ceilometer and frequently cannot report any MLH. In contrast, STRAT handles this better in that it incorporates SNR enhancement. It was shown that the ceilometer-derived aerosol profiles provide consistent with expected MLH information, which can be used to trace the urban MLH dynamics during day. However, the Jenoptik algorithm has difficulty (low availability of reported MLH due to reduced SNR when ML is high) primarily due to the early version of the software used in this work. The hampered tracking of the MLH by the proprietary algorithm may also be a result of the real-time operation on a single profile without making use of the previously collected data. The primary issues of both LiDAR-based techniques were identified as layer attribution, particularly at night and during transition periods when high aerosol layers were mistakenly used by the algorithms. It was underlined that incomplete overlapping of the ceilometer impacts the detection of low MLH at night. Based on the performed statistical analysis it was shown that the STRAT algorithm reconstructs expected MLH dynamics during the day, with maximums in the late afternoon. On the other hand, the Jenoptik method rarely reports MLH values in daytime, which embarrasses the tracking of the MLH diurnal evolution.

**Funding:** This research and publication costs were funded by the Bulgarian National Science Fund grant number DM 04/1 2016.

**Acknowledgments:** This study would not be possible without TOPROF—European COST action ES1303 and the advices and recommendations of all TOPROF members. The author is also grateful to NOAA's National Centers for Environmental Information for providing the IGRA. Acknowledgements are due to all contributors to the R project. The author would also like to thank anonymous readers whose valuable comments and corrections significantly improved paper quality.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Atmosphere* Editorial Office E-mail: atmosphere@mdpi.com www.mdpi.com/journal/atmosphere

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-2961-5