Diurnal Evolution of the Wintertime Boundary Layer in Urban Beijing, China: Insights from Doppler Lidar and a 325-m Meteorological Tower

: The diurnal evolution of the atmospheric boundary layer—the lowermost part of the atmosphere where the majority of human activity and meteorological phenomena take place—is described by its depth. Additionally, the boundary layer height (BLH) and the turbulence intensity strongly impact the pollutant di ﬀ usion, especially during transition periods. Based on integrated observations from a 325-m meteorological tower and a Doppler Wind lidar in the center of Beijing, the entire diurnal cycle of urban BLH in December 2016 was characterized. Results highlight that the Doppler lidar exhibited it is well suited for monitoring convective BLH while it trudges in monitoring stable BLH, while a 325-m meteorological tower provided an important supplement for Doppler lidar under nocturnal boundary layer and heavily polluted conditions. For the diurnal cycle, under light wind condition, the pattern of urban BLH was largely modulated by thermal forcing of solar radiation and may partly be a ﬀ ected by wind speed. While under strong wind condition, the pattern of urban BLH was largely modulated both by thermal forcing and dynamical forcing. The present work also presented evidence for several new features in the morning and afternoon transitions of the urban boundary layer, showing the duration of the morning transition varied between 1 and 5 h, with the largest value occurring under weak wind with high PM 2.5 concentration; while the afternoon transition ranged from 3 to 6 h, which was positively (negatively) correlated to wind speed (PM 2.5 concentration). Our work highlights that weak wind speed (weak dynamic motion) and heavy aerosol pollution (weak thermal forcing due to the e ﬀ ect of cooling) can dramatically a ﬀ ect the evolution of the boundary layer.


Introduction
The atmospheric boundary layer (ABL) height is an important factor in diluting pollutants emitted at the surface, in which layered pollutants vertically disperse through the action of turbulence. Accurate estimation of ABL diurnal variability is critical to predict pollutant concentrations by model computations [1][2][3][4], especially for the ABL over urban environment, so-called the urban boundary layer [5][6][7]. The most common way to determine the boundary layer height (BLH) is by using vertical profiles of meteorological variables, such as the potential temperature, wind speed (WS), and bulk Richardson number (Ri b ), which are generally obtained from radiosonde launches or aircraft observations [8]. However, radiosonde observations are usually infrequent and rarely made in central urban environments [9]. Mostly, radio soundings collected at rural sites nearby cities are used to assist in the analysis of the urban boundary layer [10,11]. Particularly, spatial differences in daytime BLH have been reported due to various surface properties with different surface heat flux characteristics [12,13]. Thus, the BLH estimated from radio soundings launched at rural sites is not representative of urban areas. Moreover, in the frame of the World Meteorological Organization (WMO), the radio soundings are launched only twice per day. In few other locations, i.e., during field campaigns and densified observation, radio-sounding are launched three times a day and therefore cannot capture the diurnal development of the urban boundary layer.
In the last decades, remote sensing techniques and network of instruments became the most valuable and popular approach for detecting the nature of the atmosphere [14], particularly in the urban layer, because of their higher temporal resolutions, such as ceilometers [15,16], sodars [17,18], and lidars [1,11]. The Doppler wind lidar is an optical active remote sensing instrument that became popular after the discovery of the CO 2 laser [19,20]. Compared to other lidar types (e.g., elastic and Raman), Doppler wind lidar is the most suited for the investigation of the vertical wind structures, the development of the urban boundary layer [3,[19][20][21][22][23][24][25], and for validating numerical forecast models [26] as it has the capability to measure the 3D wind vector field with the maximum measurement range exceeding 10 km. For instance, the vertical velocity variance (σ 2 w ), obtained from the 3D wind measurements, can be used to describe the turbulence intensity by the turbulence method, and the BLH can be defined. The urban mixing height was defined by using a threshold method (the so-called threshold method) as the height up to which σ 2 w > 0.1 m 2 s −2 in London, UK [21]. A previous study pointed out that the threshold method provides a reliable height retrieval for convective boundary layer (CBL) with strong turbulence and unstable atmospheric stratification but is not recommended to determine the nocturnal boundary layer (NBL) height, when usually the turbulence is weak and often has σ 2 w < 0.1 m 2 s −2 [27]. For another instance, a Doppler wind lidar was set in urban Beijing during the summer of 2015. Our previous study also evaluated a new improved version of the threshold method [28], to determine the NBL height, which defined the NBL top as the height at which σ 2 w decreases to 10% of its near-surface maximum after subtracting a background value [27]. Based on a case study under heavily polluted conditions (1)(2)(3)(4) December 2016), we found that these two threshold methods in London and Beijing were both inaccurate in determining the NBL height [29]. The former was also proven wrong for severely polluted urban boundary layer conditions with weak turbulence during daytime in winter in Beijing. This might have been because of the lower urban BLH dropping below the minimum observable height (100 m) during the time of the NBL. Additionally, the winter NBL over Beijing is mostly steady owing to surface radiative cooling [30,31], which does not satisfy the near-neutral assumption for the method developed by [27]. As is well known, the NBL under fine weather is the most common example of a stable boundary layer (SBL). The NBL often has a complex structure, especially over polluted urban canopies [32], which makes it far more difficult to retrieve the SBL height by using a solo Doppler lidar. For severely polluted cities, it is therefore necessary to provide measurements of the near-surface layer properties under SBL conditions additional (or auxiliary) to those by Doppler lidar to be able to reliably estimate the diurnal cycle of the BLH.
A tall meteorological tower in Beijing, officially known as "the Beijing 325-m meteorological tower", is located in the city center, providing 15 levels of measurements of WS and temperature. The main limitation of the tower platform is the observational height. The advantages and disadvantages of the Doppler lidar and tower instruments suggest that urban BLH detection might be improved by using the two instruments in synergy. In this work, by combining the retrievals from the co-located Doppler wind lidar and the 325-m meteorological tower, we aim to estimate the urban BLH and investigate the development of urban BLH in Beijing, especially during heavily polluted conditions, and characterize its full diurnal cycle in winter conditions.

Experimental Sites and Instrumentation
The 325-meteorological tower and the Doppler wind lidar (Windcube200, Leosphere, Orsay, France) are co-located at the Institute of Atmospheric Physics (IAP), Chinese Academy Sciences in Beijing city center (39.  tower", is located in the city center, providing 15 levels of measurements of WS and temperature. The main limitation of the tower platform is the observational height. The advantages and disadvantages of the Doppler lidar and tower instruments suggest that urban BLH detection might be improved by using the two instruments in synergy. In this work, by combining the retrievals from the co-located Doppler wind lidar and the 325-m meteorological tower, we aim to estimate the urban BLH and investigate the development of urban BLH in Beijing, especially during heavily polluted conditions, and characterize its full diurnal cycle in winter conditions.

Experimental Sites and Instrumentation
The 325-meteorological tower and the Doppler wind lidar (Windcube200, Leosphere, Orsay, France) are co-located at the Institute of Atmospheric Physics (IAP), Chinese Academy Sciences in Beijing city center (39.97°N, 116.37°E) (See Figure 1a).
The measurements are taken at 15 levels (i.e., 8 m, 15 m, 32 m, 47 m, 65 m, 80 m, 100 m, 120 m, 140 m, 160 m, 180 m, 200 m, 240 m, 280 m, and 320 m above ground level) by a suite of meteorological instruments (010C cup anemometers, 020C wind vanes, Metone, USA, HC2-S3, Rotronic, Switzerland) deployed on the 325-m meteorological tower, which provides observations of the WS, wind direction, relative humidity (RH), and air temperature (T) measurements. In addition, turbulence measurements are taken using sonic anemometers (Model Windmaster Pro, Gill, UK) and four-component radiation using downward-pointing and upward-pointing pyrgeometers, and pyranometers (CNR1, Kipp & Zonen) are maintained at the same heights (i.e., 47 m, 140 m, and 280 m) on this tower. The detectable range of the lidar spans from 100 m to several kilometers above the ground (depending on the aerosol concentration), with a vertical range resolution of 50 m.  In addition, the PM 2.5 mass concentrations measured at the Beijing Olympic Sports Center (Aoti surface station), located approximately 2 km north of IAP station, were obtained from the China National Environmental Monitoring Center website (http://106. 37.208.233:20035/). Cloud fractions observations were from the Nanjiao weather station, about 20 km from the IAP station.
As shown in Figure 1b, the IAP station is surrounded by a highly dense area of tall buildings. Around this station, there are about 92.1% impervious built-up surfaces and commercial and residential tall buildings (10-60 m), within a radius of 5 km of the tower [33,34]. By using the 47 m level wind measurements, we can speculate that strong northerly winds prevailed in winter, while southerly winds are lighter with a speed lower than 3 m s −1 .

Data Processing
The original 3D sonic anemometer records (10 Hz) were processed prior to analysis using the methods of double rotation (i.e., yaw and pitch rotations) and linear detrending. A density effect correction [35], planar fit method [36], and an averaging period of one hour [37] was applied in the analysis. A criterion of threshold carrier-to-noise ratio of −20 dB [21,27] was chosen to reduce the effects of invalid data in this dataset of the Doppler wind lidar. Based on the edited data, the σ 2 w over 30-min segments were calculated [21,27].

Estimation of Urban BLH
Based on radiosonde data, the Ri b method is widely used, as it is a reliable method to determine both the CBL and SBL [38]. Compared with the methodology based on detecting the altitude where there is a sharp change in the temperature and RH profiles, the Ri b method has been found to be more effective, as shown in [39]. The WS and T at 15 levels provided by the 325-m tower make this method available in our study under the condition of BLH < 300 m. The definition of Ri b was given as [40]: where z is height (m), g is the gravity acceleration constant (m s −2 ), θ v is the virtual potential temperature (K), u and v are the horizontal components of WS (m s −1 ), and u * is the surface friction velocity (m s −1 ). In addition, z s is the height of the lower boundary (m), and b is an empirical coefficient. According to [40], z s can be set to 20 m or 40 m, and b can be set to 100 m. The Ri b method defines the BLH as the height at which the Ri b reaches the critical threshold value of 0.25 (the vertical change of Ri b is shown in Figure S1). It is noted that a strong correlation [correlation coefficient (r), r = 0.83] between the BLH derived with z s = 20 m and z s = 40 m. In this paper, we choose z s = 20 m. In many previous studies, b was set to 0 because of the unknown value of u * , as in the experiments [41][42][43]. In our study, u * can be calculated as u w 2 + v w Zhang et al. [44] also found that the surface friction effect had a significant impact on the calculation of BLH, especially for the height of the SBL. Therefore, b = 100 was used in our study.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 19 layer. Zhang et al. [44] also found that the surface friction effect had a significant impact on the calculation of BLH, especially for the height of the SBL. Therefore, b = 100 was used in our study. The standard deviation ( ) of the 30-min vertical velocity was calculated by using the measurements from Doppler wind lidar. The threshold method of > 0.1 m 2 s −2 derived by [13], as mentioned in the Introduction, was used to calculate urban BLH based on the lidar.

Calculation of Heat Flux
From the sonic anemometer measurements, the turbulence kinetic energy (TKE) per unit mass of flow (TKE/m, m 2 s −2 ) and sensible heat flux ( , W m −2 ) were calculated as: where , , and are the air density (kg m −3 ), the specific heat capacity at a constant pressure (J kg −1 K −1 ), and the virtual air temperature (K), respectively.

Urban BLH
The threshold method, based on the lidar data, was not able to yield the NBL height at all times throughout the measurement campaign, because the urban BLH (averaged in one hour) was too low (<100 m, not shown). This method was either available for the urban boundary layer retrievals under daytime heavy pollution conditions. Considering as a case study a red-alert pollution event (the air quality index reaches more than 300, which will last for at least three consecutive days) from 16 to 22 December (Figure 3), rarely the threshold method worked in determining the urban BLH during daytime with a small , especially under heavily polluted conditions (daytime mean concentrations of PM2.5 > 200 μg m −3 on 18-21 December). In contrast, the method successfully retrieved the NBL height and even the height of the polluted urban boundary layer during daytime, which shows reasonable results by reference to the change in . Therefore, in our study, the NBL height was estimated by using the method. Considering that the blind overlap lidar region was 100 m, the results of urban BLH from threshold method, which was lower than 150 m, were not considered. Note that the two methods are different, while both of them are based on similar physical parameters related to turbulence processes. The standard deviation (σ w ) of the 30-min vertical velocity was calculated by using the measurements from Doppler wind lidar. The threshold method of σ 2 w > 0.1 m 2 s −2 derived by [13], as mentioned in the Introduction, was used to calculate urban BLH based on the lidar.

Calculation of Heat Flux
From the sonic anemometer measurements, the turbulence kinetic energy (TKE) per unit mass of flow (TKE/m, m 2 s −2 ) and sensible heat flux (H, W m −2 ) were calculated as: where ρ, C P , and T v are the air density (kg m −3 ), the specific heat capacity at a constant pressure (J kg −1 K −1 ), and the virtual air temperature (K), respectively.

Urban BLH
The threshold method, based on the lidar data, was not able to yield the NBL height at all times throughout the measurement campaign, because the urban BLH (averaged in one hour) was too low (<100 m, not shown). This method was either available for the urban boundary layer retrievals under daytime heavy pollution conditions. Considering as a case study a red-alert pollution event (the air quality index reaches more than 300, which will last for at least three consecutive days) from 16 to 22 December (Figure 3), rarely the threshold method worked in determining the urban BLH during daytime with a small σ 2 w , especially under heavily polluted conditions (daytime mean concentrations of PM 2.5 > 200 µg m −3 on 18-21 December). In contrast, the Ri b method successfully retrieved the NBL height and even the height of the polluted urban boundary layer during daytime, which shows reasonable results by reference to the change in σ 2 w . Therefore, in our study, the NBL height was estimated by using the Ri b method. Considering that the blind overlap lidar region was 100 m, the results of urban BLH from threshold method, which was lower than 150 m, were not considered. Note that the two methods are different, while both of them are based on similar physical parameters related to turbulence processes. In case that only one method is available to retrieve the boundary layer height, that retrieved value is taken as reference. Instead, when the boundary layer height could be retrieved with the two methods, the intercomparison is shown in Figure 4. Combined with Figure 3, it was found that mostly the boundary layer height retrieved by the Ri b method is higher with respect to the retrievals obtained with the threshold method, and the urban BLHs from threshold method is rather insensitive for urban BLH below about 200 m. The correlation is essentially nonexistent (r 2 = 0.16). Considering that the threshold method may underestimate the depth of mixing layer/the urban boundary layer under weak turbulence context [21,23,27], we use the threshold (200 m) of the definition of the onset of CBL [45] and the height of the meteorology tower (325 m) as the references to choose the different results as the real values. When the height could be retrieved from both the two methods, the urban BLHs retrieved from the Ri b method was considered as reference (4% of the total values) if the urban BLHs from threshold method were lower than 200 m, and the averaged values were chosen if the urban BLHs from threshold method were larger than or equal to 200 m, and lower than 325 m (3% of the total values) (see Figure 4). About 1% of the urban BLHs retrievals from threshold method were higher than 325 m and also than the values from Ri b method, which happened during the clean noon time with relatively larger turbulence. Thus, in this condition, the values from threshold method were considered credible and chosen as the real ones.
In addition, through comparison against simulated urban BLH derived by three different BL parameterization schemes (Table S1 in supporting information) in Weather Research and Forecast (WRF) model, it is shown that similar diurnal patterns of urban BLH were captured by our present observation (See Figure S2 in supporting information), implying that the physical process of boundary layer evolution based on our observation is reasonable. While the biases are mainly related to different BLH definitions causing the divergence among different schemes [46]. More details of BLH simulations are available in the study of [46]. Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 19 7  . Considering that was found to account for a much larger amount of the net radiation (40%) than the latent heat flux (10%) during this study period over Beijing, only was shown to represent the heat transfer intensity in this paper. Generally, during this observational period, the urban BLH was positively linked to the near surface WS at all the 15 levels (e.g., r =0.3and r =0.2 for 8 m and 320 m, respectively), and the largest correlation coefficient, r =0.4, was occurred at 47 m and 65 m, and then decreasing with the increasing height. The Urban BLH was also positively linked to DSR (r=0.6), (r=0.5), * (r=0.5), and TKE/m (r=0.5), but negatively linked to the PM2.5 concentrations (r =−0.4). A high urban BLH occurred on strong, windy days with a large amount of solar radiation (e.g., 5 December and 23 December). Note that, for very clean windy days with very low aerosols concentration and strong wind as shown in Figure 5b (e.g., 5 December, 8 December, 15 December, and 22 December), the lidar signals showed a low signal-to-noise, resulting in not available data, because the heterodyne Doppler wind lidar technique using an infrared laser source relies on the presence of aerosols [21,27], leading to discontinuous observations as shown in Figure  5a. The maximum urban BLH during the whole winter month was 1000 m, accompanied by the maximum DSR, approximately 480.7 W m −2 , on 23 December, which is half that reported on clear days during summer in Beijing by using the threshold method [27]. For a heavily polluted urban boundary layer with weak wind, e.g., 4 December, 18-21 December and 31 December, the maximum urban BLH ranged from 260-300 m, and that of DSR from 240-390 W m −2 , during daytime. For the NBL, the heights were mostly around or lower than 200 m, and even lower than 100 m during the heavily polluted conditions.  Considering that H was found to account for a much larger amount of the net radiation (40%) than the latent heat flux (10%) during this study period over Beijing, only H was shown to represent the heat transfer intensity in this paper. Generally, during this observational period, the urban BLH was positively linked to the near surface WS at all the 15 levels (e.g., r = 0.3 and r = 0.2 for 8 m and 320 m, respectively), and the largest correlation coefficient, r =0.4, was occurred at 47 m and 65 m, and then decreasing with the increasing height. The Urban BLH was also positively linked to DSR (r = 0.6), H (r = 0.5), u * (r=0.5), and TKE/m (r = 0.5), but negatively linked to the PM 2.5 concentrations (r = −0.4). A high urban BLH occurred on strong, windy days with a large amount of solar radiation (e.g., 5 December and 23 December). Note that, for very clean windy days with very low aerosols concentration and strong wind as shown in Figure 5b (e.g., 5 December, 8 December, 15 December, and 22 December), the lidar signals showed a low signal-to-noise, resulting in not available data, because the heterodyne Doppler wind lidar technique using an infrared laser source relies on the presence of aerosols [21,27], leading to discontinuous observations as shown in Figure 5a. The maximum urban BLH during the whole winter month was 1000 m, accompanied by the maximum DSR, approximately 480.7 W m −2 , on 23 December, which is half that reported on clear days during summer in Beijing by using the threshold method [27]. For a heavily polluted urban boundary layer with weak wind, e.g., 4 December, 18-21 December and 31 December, the maximum urban BLH ranged from 260-300 m, and that of DSR from 240-390 W m −2 , during daytime. For the NBL, the heights were mostly around or lower than 200 m, and even lower than 100 m during the heavily polluted conditions. The cloud fraction is marked as red crosses (cloudy) and black crosses (clear).

Based on the retrieved vertical distribution of the aerosol extinction coefficient from a Mie-elastic backscatter polarization lidar (AGHJ-1-LIDAR) developed by Anhui Institute of Optics and Fine
Mechanics (AIOFM) also located at the IAP site, the maximum retrieved urban BLH during daytime on 19 and 20 December was 350 m and 300 m in [31], respectively, which is similar to our results (300 m and 278 m). Note, however, that the urban BLH during the early morning was 500 m for these two days but 250 m and 300 m at night in [31], and about 200 m during the early morning on 19 December and 400 m at night on 20 December in [11] by using the same method as [31]. However, urban BLHs were mostly lower than 100 m and exhibited little change during the whole of nighttime (i.e., stable NBL, as shown in Figure S1) in our study. It was found that the stable NBL usually was accompanied by the temperature inversion near the surface during heavily polluted episodes, e.g., during 2-5 December, 16-22 December in [47]. It is important to note that a stable NBL is often coupled with a residual layer (RL) [48]. The aerosol scattering coefficient method probably defines the RL as the stable NBL [49][50][51]. Thus, attention should be given when comparing aerosol-and turbulence-derived urban BLHs [52], especially when evaluating results from model computations. Nonetheless, the aerosol scattering method provides a useful reference for the RL and NBL, which can help provide an indication of the change of the whole ABL.
To reduce the impact of the diurnal change in meteorological parameters on the urban BLH, the values were averaged for daytime and nighttime. The daily mean urban BLH was in the range of 130-536 m and 64-268 m during daytime and nighttime, with monthly averaged values of 278 m and 152 m, respectively. During the heavily polluted event, the daytime mean urban BLH decreased to <200 m. The lowest BLH during daytime occurred on 18 December, with the lowest WS (0.46 m s −1 ) at the 8-m level and notably poor turbulent activity (H = 15.16 W m −2 , u * = 0.12 m s -1 and TKE/m = 0.11 m 2 s −2 at the 47-m level), at which the corresponding daily mean urban BLH during nighttime was only 65 m. The correlation coefficients were calculated between the urban BLH and WS, H and turbulence parameters (u * and TKE/m), which were 0.6, 0.6, 0.6, and 0.5 during daytime, and 0.6, 0.2, 0.6, and 0.5 at nighttime, respectively. The correlation coefficient between BLH and DSR was 0.6 during daytime. All these above results indicate WS greatly contributed to the change in urban BLH during both daytime and nighttime, because of its direct positive effects on H, u * and TKE. For instance, on 17 December 2016, weak high-pressure system is to the southeast of Beijing, resulting in southwesterly prevailing weak winds over Beijing ( Figure S3). This will induce a weak mixing with lower BLH ( Figure 5). In contrast, on 22 December 2016, under the confluence of the high pressure system and the low pressure system ( Figure S3), there is a northwest jet in the boundary layer over Beijing, leading to a higher wind speed and strong mixing with higher BLH (Figure 5).
Previous studies demonstrated that pollutant aerosol concentration is directly linked to the urban BLH in Beijing [53,54]. In our study, the correlation coefficient between urban BLH and PM 2.5 concentration was −0.7 for daytime (at 0.05 significance level) but only −0.5 for nighttime (at 0.05 level significance). This relatively weak correlation for nighttime may have resulted from the small variation range of the BLH during most nights. In addition, as shown in Figure 6, a shallow urban BLH under polluted conditions was usually accompanied by high RH, suggesting that such a humid and shallow BLH favored the accumulation of PM 2.5 . This fact agrees with previous study results in which the generation and transformation of secondary aerosols were found to be significant during humid conditions, leading to an increase of PM 2.5 concentrations and the resultant occurrence of severe air pollution episodes [55,56]. The daytime urban boundary layer was clean (PM 2.5 < 50 µg m −3 ) once the daily mean urban BLH reached 300 m or above. For nighttime, when the urban BLH was higher than 200 m, the concentration of PM 2.5 was mostly lower than 100 µg m −3 . Meanwhile, it should be noted that for some shallow nighttime urban BLHs (e.g., 9-10 December), defined as urban BLH < 150 m, concentrations of PM 2.5 were lower than for others, partly explainable by the lower RH, which implies the level of pollution can be attributed to a humid or dry urban boundary layer. Additionally, compared to daytime, aerosol pollutants accumulated more easily in the urban boundary layer at nighttime in Beijing.

Diurnal Development of the Urban Boundary Layer
The results presented in the previous section suggest that the near surface WS has a significant effect on urban boundary layer variations in Beijing. Avoiding tall building disturbances around the 325-m tower, the daily average WS measured at the 47-m level was used to classify the different wind conditions. Accordingly, eight, fifteen, and eight days (marked in Figure 7) were categorized as low, medium, and high windy days based on the 25th, 50th, and 75th percentiles of the WS, corresponding to 1.29 m s −1 , 1.54 m s −1 , and 2.23 m s −1 . For low and moderate wind classes, the diurnal changes of WS were not obvious, which indicates the thermal forcing dominates the evolution of the urban boundary layer. Notably, the growth of the boundary layer also influences the wind evolution in boundary layer. Using the measurements at 15 levels from the 325-m tower in 20 years, it was found that the maximum WS often appeared at daytime at lower levels, while at higher levels the maximum WS usually occurred at nighttime [33], which was resulted from the surface stability of the urban boundary layer. Based on the case study (1-4 December 2016), it was also found that the stagnating urban boundary layer during heavy aerosol polluted daytime in winter seemed to act like an umbrella, preventing the air flow entrainment from the upper level to the lower level [29]. As shown in Figure 7a, different from the strong wind class, the WS showed a downward trend between sunrise and noon for these two wind classes, especially for the weak wind class, which may partly result from the evolution of urban boundary layer.
The diurnal variations of the urban BLH were consistent with those of thermal forcing factors, e.g., solar radiation and sensible heat flux (Figure 7). Especially for weak wind class, with unapparent diurnal change of WS and TKE, the diurnal pattern of urban BLH was mainly modulated by thermal forcing (DSR and H). While under strong wind condition, the diurnal pattern of urban BLH was largely modulated both by dynamical forcing (WS and TKE) and thermal forcing (DSR and H). As shown in Figure 7b1, the DSR decreased with the decreasing wind speed level, which can be attributed to a combination of aerosol pollutants and cloud. These low or moderate wind days were mostly accompanied by high PM 2.5 concentrations, e.g., 200 µg m −3 and 100 µg m −3 median values, as illustrated in Figure 7b3. For the polluted days, particularly the heavily polluted low-wind days, they were usually covered with cloud, as shown in Figure 3, according to the cloud fraction data. However, it is hard to acquire accurate cloud information during heavily polluted periods. In addition, the existing clouds could be the result of polluted aerosols acting as nuclei for cloud formation, during heavy polluted days, more clouds were indeed found over Beijing in case studies [57]. Here, we considered the effects of both aerosol pollutants and cloud on the DSR. On this basis, reduced DSR and weak wind led to weak thermal and dynamic forcing, resulting in a smaller diurnal range of the urban boundary layer on low-wind days than on high-wind days. In other words, compared to strong wind classes as shown in Figure 7f1-f3, a reduction of approximately 400 m and 500 m to the maximum (median values) of the urban BLH for low and moderate classes, respectively, resulted from the combined changes in the wind, the aerosol pollutant concentration, and the cloud fraction. Therefore, it was considerably challenging to quantitatively analyze the radiative feedback effects of aerosol pollutants on the urban boundary layer in Beijing by only using the experimental measurements. As such, model studies should be conducted in the future to help address this issue.
It can also be seen from Figure 7f that the top of the 325-m tower at most times was above the NBL for all the observational period, and the boundary layer during daytime for some cases with low wind in winter, suggesting that the measurements at some of the top levels may have been decoupled from the surface. As mentioned above, the time of CBL onset is defined at the time when the BLH reaches 200 m [39]. Based on this definition, the time of CBL onset got later with weakening wind, and the peak times showed a similar tendency (Figure 7). The growth rate of the median value of urban BLH also differed among the three classes; for example, from 09:00 China Standard Time (CST) to the peak time, the rates were 105 m h −1 , 48 m h −1 , and 19 m h −1 . Moreover, for the strong wind class, some NBL heights were able to reach 200 m in the early morning before sunrise with larger TKE as shown in Figure 7d3, indicating that wind shear production may contribute greatly to the larger BLH during nighttime than for the other wind classes.  a1,b1,c1,d1,e1,f1), moderate (a2,b2,c2,d2,e2,f2), and strong (a3,b3,c3,d3,e3,f3) wind classes. The shaded areas, blue solid line, red dashed line, and pink lines correspond to the interquartile range (25%-75%), median, 200 m (definition of onset the convective boundary layer), and the slopes of growth rate of the median value of urban BLH from 09:00 CST to the peak time, respectively.

Durations
The CBL starts to grow slowly after sunrise, and a transition feature of both the CBL and SBL, named "morning transition", is the first stage of the development of the ABL. Examination of the transition periods (e.g., passing from night to day or day to night) has important implications for air quality [58,59]. In previous studies, different algorithms were used to define the morning transition (MT) from observations [57,60,61]. It was defined as the period from the time when the surface buoyancy flux first becomes positive (crossover) to the time at which the BLH reaches about 200 m [57] based on the observations from meteorological tower. Considering the positive heat flux also occurred in nighttime and the urban BLH was often greater than 200 m for urban site during summertime, the beginning of the MT was marked as the sunrise, and the end of MT was defined as the rapid onset of the rapid-growth phase in [58] by using the measurements from a tall tower in London. Given that the positive H occasionally occurred during nighttime in our study period, we defined the sunrise as the starting time of MT. We chose the time when the urban BLH reaches about 200 m as the end of the MT after the sunrise in this wintertime period. The period between the maximum H at midday and H < 0 around sunset is usually defined as the afternoon transition (AT) [8,62].
In this paper, we adopted these two methods to estimate the durations of MT and AT during wintertime over horizontally homogeneous Beijing urban surface. Note that the fixed height (200 m) method for MT employed in the present study may not be robust to investigate precise absolute change characteristics of the urban boundary layer. However, at least, it can help us to study the relative change characteristics of the morning development of the urban boundary layer under different conditions (e.g., wind, solar radiation, pollutants) quantitatively. Similarly, AT is also employed to show relative change characteristics of convective urban boundary layer at daytime.
Excluding strong-wind nights (e.g., 5, 15, and 22 December), the MT duration (D MT ) ranged from 1 to 5 h, and the average value was 2.4 h. The maximum D MT occurred on the heavily polluted days, such as 4 and 18 December, with weak morning-averaged winds of 0.6 m s −1 and 0.8 m s −1 , respectively. The variation of D MT was negatively correlated (r = −0.4 for both the 8-m level and 47-m level) with the change in WS, which is consistent with previous studies [57,58,63], and implies that weak wind prolongs the D MT . However, this correlation was weaker than that in [58] (r = −0.7), possibly because of the complexity of the aerosol loading. A positive correlation (r = 0.5) was found between D MT and PM 2.5 concentration, which suggests that the significant loading of the aerosol pollutants had large effects on D MT by reducing the net radiation reaching the surface. This can be proved by model results from a large-eddy simulation in which aerosols delayed the CBL onset in the morning by up to 4.5 h [64]. Note that weak wind is usually associated with high concentrations of PM 2.5 over Beijing; and besides, many other contributing factors, such as stability, advection, the state of the previous day's RL [8], and the entrainment [65], can affect the D MT . The delay in our case was likely caused by a combination of multiple factors.
The AT started at around 11:00-14:00 CST and mostly ended about 1 h before sunset (16:50 CST, monthly averaged), meaning it could last 3-6 h. Interestingly, the AT duration (D AT ) showed no obvious correlation with WS (r < −0.1), while the ratio of H decay (r ∆H ) during D AT showed a strong positive correlation with the WS averaged during D AT (r = 0.5 for the 8-m level and r = 0.6 for the 47-m level). In addition, r ∆H had a negative correlation (r =−0.4) with the PM 2.5 concentration averaged during D AT , due to the aerosol cooing effect. In contrast, as compared with a new stable boundary-layer development near the ground in the morning, the boundary layer in the evening was already developed after a growth in the afternoon [66]. This difference will also modulate the magnitude of duration, due to different thermal-dynamical and aerosol context between morning and evening.
Notably, it has been proven that there are interactions between the urban boundary layer and aerosol pollutants [31,67], from which it can be concluded that the long D MT and rapid r ∆H both support a suitable condition for the accumulation of aerosol pollutants, and in turn, the high concentrations of aerosol pollutants, with their radiation effects, can prolong the D MT and accelerate the r ∆H . Therefore, the transition process is very complex [50], and aerosol pollutants add another level of complexity.
Based on the above findings, a schematic illustration of the diurnal cycle of the ABL on a polluted weak and windy day (taking 19 December as an example) is shown in Figure 8. The RL height (250 m) was taken from [31]. As can be seen, the CBL onset is at 11:00 CST, with a 4-h D MT , which is larger than the mean value (2.4 h). This delayed development is caused by the weak wind (0.7 m s −1 at the 8-m level, averaged during D MT ) and cooling effect of aerosols (PM 2.5 concentration = 267 µg m −3 , averaged during D MT ). Note that the CBL height tends to decrease obviously from about 4 h after the peak time of H, which implies that the response of the urban BLH variations lag behind the thermal forcing. This result is consistent with previous results derived from observational measurements [8] and a large-eddy model [50]. Moreover, the crossover of H is at 18:00 CST, 1 h after sunset, which can be explained by the weak wind and the high thermal conductivity of the concrete in built-up areas (large heat storage of the urban canopy). The existence of pollutants in the RL, formed on the previous polluted day, plays an important role in the changes of pollutants at the surface in the morning, which can reach the ground the following morning through convection [29,68]. This mechanism of pollutant vertical mixing is very complex, and can be further investigated via model simulations.

WS and Vertical Variance Profiles.
Composite vertical profiles of WS and averaged in the three wind classes during the MT and AT periods are illustrated in Figure 9. For the median WS profiles, WS was mostly below 5 m s −1 and 7 m s −1 for weak and moderate wind classes, respectively, while it could reach 17 m s −1 at about 1400 m for the strong wind class. Obviously, during MT, low-level jet (LLJ) was detected in all three wind classes, which was mainly attributed to large-and local-scale forcings, e.g., synoptic pattern, local topographic flows, buoyancy, friction, and entrainment processes [69][70][71]. For the weak wind class, an LLJ occurred at about 400-600 m. Previous studies have revealed that the LLJ can affect the transition of aerosol pollutants via advection flow and vertical mixing [72,73]. A southwesterly LLJ was found in [29], with the jet core at 300-500 m, which might have transferred polluted air from the south via advection to Beijing on 2 December; plus, this LLJ also generated significant wind shear, which could have mixed transported and local pollutants via vertical mixing. In our study, the vertical changes in the wind profile were smaller for weak and moderate wind classes in both the MT and AT periods.
For , its value was smaller in the MT period than in the AT period (Figure 9), suggesting that the vertical mixing intensity for MT was obviously weaker than for AT for all wind classes. Specifically, the median values of were < 0.1 m 2 s −2 for MT, while at AT they could reach 0.35 m 2 s −2 . During the AT period, there was an increasing between 150 m and 300 m, which can be explained by the strong wind shear as mentioned above. Combined with Figure 7e, it can be clearly seen that weak vertical mixing accompanied high PM2.5 concentrations.

WS and Vertical Variance Profiles
Composite vertical profiles of WS and σ 2 w averaged in the three wind classes during the MT and AT periods are illustrated in Figure 9. For the median WS profiles, WS was mostly below 5 m s −1 and 7 m s −1 for weak and moderate wind classes, respectively, while it could reach 17 m s −1 at about 1400 m for the strong wind class. Obviously, during MT, low-level jet (LLJ) was detected in all three wind classes, which was mainly attributed to large-and local-scale forcings, e.g., synoptic pattern, local topographic flows, buoyancy, friction, and entrainment processes [69][70][71]. For the weak wind class, an LLJ occurred at about 400-600 m. Previous studies have revealed that the LLJ can affect the transition of aerosol pollutants via advection flow and vertical mixing [72,73]. A southwesterly LLJ was found in [29], with the jet core at 300-500 m, which might have transferred polluted air from the south via advection to Beijing on 2 December; plus, this LLJ also generated significant wind shear, which could have mixed transported and local pollutants via vertical mixing. In our study, the vertical changes in the wind profile were smaller for weak and moderate wind classes in both the MT and AT periods.
For σ 2 w , its value was smaller in the MT period than in the AT period (Figure 9), suggesting that the vertical mixing intensity for MT was obviously weaker than for AT for all wind classes. Specifically, the median values of σ 2 w were < 0.1 m 2 s −2 for MT, while at AT they could reach 0.35 m 2 s −2 . During the AT period, there was an increasing σ 2 w between 150 m and 300 m, which can be explained by the strong wind shear as mentioned above. Combined with Figure 7e, it can be clearly seen that weak vertical mixing accompanied high PM 2.5 concentrations. Figure 9. Mean diurnal cycles of WS, (a1,a2,a3) and vertical variance vertical profiles (b1,b2,b3), for weak, moderate, and strong wind classes, respectively. The shaded areas correspond to the interquartile range (25%-75%) during morning (blue) and afternoon transitions (red).

Discussions
Note that the aerosol-urban boundary layer interaction is a two-way feedback process [73][74][75]: (1) the shallow urban boundary layer can weaken the diffusion and dilution condition, accumulating the pollutants, and the resultant aerosol concentration increases at the near-surface. (2) The aerosols can reduce surface solar radiation, suppressing boundary layer development due to the weakened surface heat flux and the weakened thermal forming of underlying surface [34,56,[74][75][76][77]. In addition, the aerosols in the atmosphere absorbing the shortwave solar radiation could change the vertical temperature stratification (i.e., heating upper air and cooling surface) and further suppress the development of BL [29,75,78]. Such suppressed boundary layer further enhances near-surface aerosol pollution in megacities. In the present study, based on the observational measurements, we mainly focused on the characteristics of the development of the urban boundary layer and concentrations of PM 2.5 for different wind classes. Some of these characteristics are resulted from the aerosol-urban boundary layer interaction.

Conclusions
In the present study, based on integrated observations from a 325-m meteorological tower and a Doppler lidar in urban Beijing, the entire diurnal cycle of the urban BLH in December 2016 was characterized. The highest daytime urban BLH observed in the measurement campaign was as high as 1000 m; however, under heavily polluted conditions with weak wind, the maximum daytime urban BLH reached up to 300 m only. The Doppler lidar is well suited to retrieve and monitor CBL evolution, while it trudges in monitoring the SBL process. A 325-m meteorological tower provided an important supplement on SBL observation from Doppler lidar, especially under NBL and heavily polluted conditions (urban BLH around or lower than 200 m, and even lower than 100 m). The variations of DSR, WS, H, u * , and TKE correlated positively with the BLH, while there was a negative correlation between BLH and PM 2.5 .
Based on the full image of the boundary layer diurnal cycle, we have improved the knowledge of morning and afternoon transitions of the urban boundary layer under different weather and pollution conditions. That is, the D MT varied between 1 and 5 h, with the largest value occurring under conditions of weak wind and high PM 2.5 concentration, while D AT ranged from 3 to 6 h. In addition, D AT and r ∆H presented positive correlations with WS but negative correlations with PM 2.5 concentration. Our work highlights that weak WS (weak dynamic motion) and high levels of aerosol pollutants (weak thermal forcing due to the cooling effect) can dramatically affect the development of the boundary layer. For the diurnal cycle, under weak WS condition, the pattern of urban BLH was largely modulated by thermal forcing of solar radiation and might partly by WS; while for the strong wind condition, the pattern of urban BLH was largely modulated by dynamical forcing (TKE and WS), and as well by thermal forcing.
In future work, many more examples need to be collected to study the impact of thermal and dynamic forcing on the diurnal development of the wintertime boundary layer and its association with air pollution in urban Beijing by using joint meteorological tower and Doppler lidar observations, especially under calm and heavily polluted conditions. The results reported also have important implications for other urban regions that frequently experience air pollution and may aid in designing effective regional mitigation strategies for the environment and human health at the diurnal scale. For instance, according to the different diurnal stages of boundary layer development, we can carry out staggered production strategy to modulate emissions under various urban boundary layer diffusion conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/23/3935/s1, Figure S1: Height dependence of Ri b (calculated by Equation (1)