Boundary Layer Height as Estimated from Radar Wind Proﬁlers in Four Cities in China: Relative Contributions from Aerosols and Surface Features

: The turbulent mixing and dispersion of air pollutants is strongly dependent on the vertical structure of the wind, which constitutes one of the major challenges a ﬀ ecting the determination of boundary layer height (BLH). Here, an adaptive method is proposed to estimate BLH from measurements of radar wind proﬁlers (RWPs) in Beijing (BJ), Nanjing (NJ), Chongqing (CQ), and Wulumuqi (WQ), China, during the summer of 2019. Validation against simultaneous BLH estimates from radiosondes (RSs) yielded a correlation coe ﬃ cient of 0.66, indicating that the method can be used to derive BLH from RWPs. Diurnal variations of BLH and the ventilation coe ﬃ cient (VC) at four sites were then examined. A distinct diurnal cycle of BLH was observed over all four cities; BLH gradually increased from sunset, reached a maximum in the afternoon, and then dropped sharply after sunset. The maximum hourly average BLH (1.426 ± 0.46 km) occurred in WQ, consistent with the maximum hourly mean VC larger than 5000 m 2 / s observed there. By comparison, the diurnal variation of VC was not strong, with values ranging between 2000 and 3000 m 2 / s, likely owing to the high-humidity environment. Furthermore, surface sensible heat ﬂux, latent heat ﬂux, and dry mass of particulate matter with aerodynamic diameter ≤ 2.5 µ m (PM 2.5 ) concentrations were found to somehow a ﬀ ect the vertical structure of wind and thermodynamic features, leading to a di ﬀ erence between RS and RWP BLH estimates. This indicates that the atmospheric environment can a ﬀ ect BLH estimates using RWP data. The BLH results from RWPs were better in some speciﬁc cases. These ﬁndings show great potential of RWP measurements in air quality research, and will provide key data references for policy-making toward emission reductions.


Introduction
The atmospheric boundary layer (ABL) is the lowest layer of the atmosphere that is readily affected by surface properties such as soil moisture and roughness [1][2][3][4]. The determination of accurate vertical structures of the ABL has great implications for environmental and climate studies [5][6][7]. Boundary layer height (BLH) and wind speed in the ABL are found to significantly affect the accumulation and diffusion of air pollutants [8][9][10][11][12][13]. Therefore, to build and maintain an environmentally friendly society and clean atmospheric environment, the atmospheric ventilation coefficient (VC) and BLH at short time scales should be determined [1,2].
A radar wind profiler (RWP) is an active remote sensing instrument. A number of countries have established nationwide RWP networks that can provide continuous observations of wind and estimation algorithm, the RWP data from four cities of China, each of which has one RWP, were used in this study: Beijing (BJ), Nanjing (NJ), Chongqing (CQ), and Wulumuqi (WQ). More importantly, these four cities also provide atmospheric measurements from RS, which makes it possible to verify BLH estimates from RWPs. The RWP data were used to test the proposed method, and the RS measurement was taken as the reference value for the RWP. Table 1 lists the site details. All observation data were collected from June to August 2019. To eliminate the influence of rain and clouds, only the observation data on clear-sky days were used to retrieve BLH. Rainy and cloudy days were identified when the vertical velocity at 0-3 km altitude was larger than 0.5 m/s [32,33].

L-Band Radiosonde Measurements
The L-band RS is a new-generation sounding system that can provide fine-resolution profiles of temperature, pressure, relative humidity, wind speed, and wind direction twice a day at 08:00 and 20:00 local time (LT). The profiles of meteorological parameters at 14:00 LT, particularly at BJ, were collected during several recent intensive observational periods in summer [34,35]. The RS data were obtained from June to August 2019 at the four stations. The accuracy of the radiosonde temperature profile measurements is within 0.1 K in the troposphere [27]. Humidity profiles contain large uncertainties in the presence of clouds due to the low sensitivity of the hygristor [28]. For this reason, BLH is generally calculated under cloud-free conditions, which ensures that the RS measurements are good enough to derive BLH values. Here, the bulk Richardson number (Ri) method was employed due to its robustness under stable and convective ABL regimes [36,37]. Guo et al. conducted a sensitivity analysis, and suggested that Ri should be set to 0.25 to calculate BLH [37]. Notably, the method of bulk Ri is unable to determine BLH due to a decoupled structure that frequently occurs in the developing stages of the ABL [7]. Most of the ABL throughout China remains coupled with the ground surface at 20:00 LT, especially in the summer. However at 08:00 LT, the Ri method may be not applicable in some cases. In WQ in particular, there is a 2-hour shift in local solar time due to the longitude difference. Therefore, the gradient method was used to calculate BLH from the virtual potential temperature (PT V ) profiles when the Ri method was not valid. In this case, BLH refers to the first point in the PT V profile where its vertical gradient value crosses 0.5 • /∆r, where ∆r is the range resolution of the RS [25]. The BLH from RSs was taken as the reference value for the retrievals from RWPs.

Reanalysis Data
The MERRA-2 dataset contains reanalysis data developed by NASA's Goddard Earth Sciences Data and Information Services Center [38]. In particular, the surface turbulent flux dataset, which is named "MERRA-2 tavg1_2d_flx_Nx: 2d, 1-hourly, time-averaged, single-level, assimilation, surface flux diagnostics V5.12.4," is used [39]. These data are available hourly at a spatial resolution of 0.625 • longitude and 0.5 • latitude. Surface fluxes were computed using the surface-layer turbulence parameterization described by Helfand and Schubert [40].

PM 2.5 Data
The Ministry of Environmental Protection of China has established more than 1000 environmental monitoring stations throughout the country. They routinely measure the dry mass of particulate matter with aerodynamic diameter ≤2.5 µm (PM 2.5 ), which is released to the public in real time with relatively high credibility [31]. PM 2.5 data can be obtained from the official website of the Ministry of Environmental Protection (http://data.cma.cn/en, last access: 30 January 2020).

BLH Determined from RWPs
The atmospheric refractive index structure constant (C 2 n ) and meteorological parameters have evident characteristics at the top of the boundary layer [4,41]. C 2 n reaches a maximum at BLH due to the small-scale buoyancy fluctuations associated with the entrainment process [18,25]. The value of C 2 n is directly proportional to the SNR of RWPs [22]. On this basis, various researchers have employed a wavelet transform or threshold method to invert BLH using SNR profiles and achieved results that correlate well with the RS [18][19][20][21]. However, the accuracy of these methods strongly depends on the set of threshold or initial parameters. By contrast, several researchers have proposed using the average wind profile to retrieve BLH. The boundary layer is the lowest layer of the atmosphere that is directly affected by surface properties such as soil moisture and roughness [1][2][3][4]. Hence, atmospheric motion in the boundary layer can be regarded as a turbulent phenomenon [4]. The friction effects can be neglected above the boundary layer [42]. Thus, wind speed and direction above and below the boundary layers are significantly changed. The mean or instantaneous vertical profiles of wind could be used to invert BLH [12,13,43,44]. Previous studies also indicated that the second-moment variables of wind profiles result in the local minimum or maximum being at the top of the boundary layer. These variables include the standard deviation of vertical velocity [45], standard deviation of wind direction [46], and vertical flux of the component of horizontal momentum along the surface wind direction [47]. However, those methods are used to retrieve the stable boundary layer height at nighttime, and are only verified on the observation data of the meteorological tower. Therefore, to obtain BLH automatically and robustly in multiple radar wind profilers, the normalized standard deviation method (NSDM) was developed to retrieve BLH from RWP observations.

Normalized Standard Deviation Method
The NSDM was developed based on the profile of horizontal wind speed and direction. The horizontal wind speed and direction of turbulent wind are slightly different from those of laminar wind. Moreover, turbulence intensity is related to horizontal wind speed and direction. Therefore, the profiles of horizontal wind speed and direction were used to locate turbulent boundaries. To obtain the mutation information in the vertical direction, the vertical standard deviation of horizontal wind speed and direction were calculated over three adjacent altitude bins. The red lines in Figure 1a,b represent the vertical standard deviation of horizontal wind speed and direction, and the blue circles represent local peaks. SD WS and SD WD were used to evaluate the turbulence intensity difference between adjacent altitude bins. However, the order of magnitude of SD WD was considerably larger than that of SD WS . Therefore, SD WS and SD WD were normalized by the max normalization method. This normalization can be expressed as follows: where SD WS and SD WD represent the standard deviation of horizontal wind speed and direction in the vertical direction, respectively; max(SD WS ) and max(SD WD ) represent the maximum values in SD WS and SD WD , respectively; and NSD WS and NSD WD represent the normalized standard deviation of horizontal wind speed and direction, respectively. In Figure 1c, the profile of NSD WS + NSD WD was used to describe the turbulence intensity difference in the vertical direction. Local peak points (blue circles in Figure 1c) indicate strong turbulence differences between altitude bins. This can be understood as the critical point of strong and weak turbulence or of turbulence and laminar flow. The BLH result (blue line) could be determined based on the local peak point. Figure 1d   To clarify the process of the algorithm, more cases are shown in Figure 2 under different wind vertical structures. Strong wind shear was defined when the deviation of the horizontal wind direction between adjacent sample points exceeded 45°. According to the amount of strong wind shear, vertical wind profiles were divided into three categories: one wind shear, multiple wind shears, and no wind shear. When only one local peak existed in the NSDWS + NSDWD profile, the height of the local peak was output as BLH regardless of the wind structure change. Figures 2a-c illustrate that under the three wind structures, the only peak point was directly regarded as the BLH. By contrast, local peaks were filtered when multiple local peaks formed in the NSDWS + NSDWD profile. Multiple local peaks would be formed in the vertical direction (blue circles in Figure 2d,e,g) due to the turbulent motion pattern of the ABL and the turbulence in the cumulus and free atmosphere, particularly clear-sky turbulence [4]. That is, turbulence within or above the boundary layer resulted in the appearance of multiple local peaks. Therefore, BLH should be determined from multiple local peaks. Fortunately, previous research has shown that the near-surface maximum of radial velocity variance and near-surface minimum curvature of wind speed profiles are good criteria for determining BLH under a stable boundary layer [12,43]. The vertical flux of the component of horizontal momentum and vertical velocity variance can be used to estimate BLH for a weakly stable boundary layer [12,47,48]. Therefore, a filter was designed based on vertical variation of the wind structure. This filter is used to determine whether significant wind shear occurred at the local peak To clarify the process of the algorithm, more cases are shown in Figure 2 under different wind vertical structures. Strong wind shear was defined when the deviation of the horizontal wind direction between adjacent sample points exceeded 45 • . According to the amount of strong wind shear, vertical wind profiles were divided into three categories: one wind shear, multiple wind shears, and no wind shear. When only one local peak existed in the NSD WS + NSD WD profile, the height of the local peak was output as BLH regardless of the wind structure change. Figure 2a-c illustrate that under the three wind structures, the only peak point was directly regarded as the BLH. By contrast, local peaks were filtered when multiple local peaks formed in the NSD WS + NSD WD profile. Multiple local peaks would be formed in the vertical direction (blue circles in Figure 2d,e,g) due to the turbulent motion pattern of the ABL and the turbulence in the cumulus and free atmosphere, particularly clear-sky turbulence [4]. That is, turbulence within or above the boundary layer resulted in the appearance of multiple local peaks. Therefore, BLH should be determined from multiple local peaks. Fortunately, previous research has shown that the near-surface maximum of radial velocity variance and near-surface minimum curvature of wind speed profiles are good criteria for determining BLH under a stable boundary layer [12,43]. The vertical flux of the component of horizontal momentum and vertical Remote Sens. 2020, 12, 1657 6 of 17 velocity variance can be used to estimate BLH for a weakly stable boundary layer [12,47,48]. Therefore, a filter was designed based on vertical variation of the wind structure. This filter is used to determine whether significant wind shear occurred at the local peak point. The height of the local peak point where significant wind shear occurred was output as BLH. Figure 2d reveals a case with multiple local peaks and one wind shear. The height of the local peak point, which corresponded to strong wind shear, was regarded as the BLH. When either no or multiple evident wind shears occurred, vertical velocity variance in the height of the local peak point was compared, and the local peak point close to the near-surface local minimum of vertical velocity variance was output as the BLH. Figure 2e demonstrates the case with multiple local peaks and no corresponding wind shear. By comparison, Figure 2g shows the condition with multiple local peaks and corresponding multiple strong wind shears. In these cases, the standard deviation of vertical velocity should filter the error points. The near-surface local minimum of vertical velocity variance is consistent with the BLH [13,45]. Therefore, the local peak point closest to the near-surface local minimum was indicated as the BLH (Figure 2f,h).
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 17 point. The height of the local peak point where significant wind shear occurred was output as BLH. Figure 2d reveals a case with multiple local peaks and one wind shear. The height of the local peak point, which corresponded to strong wind shear, was regarded as the BLH. When either no or multiple evident wind shears occurred, vertical velocity variance in the height of the local peak point was compared, and the local peak point close to the near-surface local minimum of vertical velocity variance was output as the BLH. Figure 2e demonstrates the case with multiple local peaks and no corresponding wind shear. By comparison, Figure 2g shows the condition with multiple local peaks and corresponding multiple strong wind shears. In these cases, the standard deviation of vertical velocity should filter the error points. The near-surface local minimum of vertical velocity variance is consistent with the BLH [13,45]. Therefore, the local peak point closest to the near-surface local minimum was indicated as the BLH (Figure 2f,h).  Figure 3 shows a flowchart of the standard deviation method. Four key steps of this method were identified: standard deviation of wind profiles, normalization of standard deviation, local peak search, and filtration of the error point. The uncertainty of this algorithm comes from two aspects, data precision and the filter model. First, given the accuracy of wind profile data, the accuracy of the wind profile directly affects the accuracy of BLH. The location of the boundary is based on the wind speed and direction profiles. Thus, for various instruments at different sites, errors of collected data directly lead to errors in the determined BLH. The uncertainty of the horizontal wind speed and direction is approximately 5%-10% [49]. Therefore, the uncertainty of the resulting BLH is  Figure 3 shows a flowchart of the standard deviation method. Four key steps of this method were identified: standard deviation of wind profiles, normalization of standard deviation, local peak search, and filtration of the error point. The uncertainty of this algorithm comes from two aspects, data precision and the filter model. First, given the accuracy of wind profile data, the accuracy of Remote Sens. 2020, 12, 1657 7 of 17 the wind profile directly affects the accuracy of BLH. The location of the boundary is based on the wind speed and direction profiles. Thus, for various instruments at different sites, errors of collected data directly lead to errors in the determined BLH. The uncertainty of the horizontal wind speed and direction is approximately 5%-10% [49]. Therefore, the uncertainty of the resulting BLH is approximately 5%-10%. Second, the filter also introduces errors. For multiple local peaks, an empirical model is used to set the filter and consequently find the BLH. However, the empirical model is suitable for most situations, but not at all cases. If the case fails to conform to the empirical model, then the BLH result of this method is unreliable. Notably, a convective cloud is accompanied by strong turbulence, which results in its boundary being misjudged as BLH. Therefore, this method is unsuitable for convective cloud conditions. approximately 5%-10%. Second, the filter also introduces errors. For multiple local peaks, an empirical model is used to set the filter and consequently find the BLH. However, the empirical model is suitable for most situations, but not at all cases. If the case fails to conform to the empirical model, then the BLH result of this method is unreliable. Notably, a convective cloud is accompanied by strong turbulence, which results in its boundary being misjudged as BLH. Therefore, this method is unsuitable for convective cloud conditions.

Analytical Method
The atmospheric ventilation coefficient can be characterized by the BLH and wind speed in the boundary layer (WSBL). Following Liu et al. [50], the VC was calculated as follows: where WSBL is the average WS within ABL, calculated by Equation (3); BLH is the boundary layer height; and WSi is the WS observed at a certain height.

Results and Discussion
In this section, we test the feasibility and evaluate the performance of the algorithm. Furthermore, we investigate the variation of boundary layer meteorology in different cities during the summer season. Finally, we analyze the influence of surface air pollution on the vertical structure of the wind field and BLH.

Case studies
Figures 4 and 5 present the case studies for BLH obtained from RWPs in BJ and NJ, respectively. Correspondingly, BLH determined from the temperature profiles was used to compare with BLH estimated from RWPs. Figure 4 reports the case study on 11 June 2019 over the BJ site; wind speed is shown in Figure 4a. Evident wind shear was revealed in the profiles of NSDWD + NSDWS at 08:00 and 20:00 LT, where corresponding BLHs were 630 and 990 m, respectively (Figure 4b). Evidently, a clear demarcation zone of wind vector fields between ABL and free troposphere, which roughly

Analytical Method
The atmospheric ventilation coefficient can be characterized by the BLH and wind speed in the boundary layer (WS BL ). Following Liu et al. [50], the VC was calculated as follows: where WS BL is the average WS within ABL, calculated by Equation (3); BLH is the boundary layer height; and WS i is the WS observed at a certain height.

Results and Discussion
In this section, we test the feasibility and evaluate the performance of the algorithm. Furthermore, we investigate the variation of boundary layer meteorology in different cities during the summer season. Finally, we analyze the influence of surface air pollution on the vertical structure of the wind field and BLH.  Figure 4 reports the case study on 11 June 2019 over the BJ site; wind speed is shown in Figure 4a. Evident wind shear was revealed in the profiles of NSD WD + NSD WS at 08:00 and 20:00 LT, where corresponding BLHs were 630 and 990 m, respectively (Figure 4b). Evidently, a clear demarcation zone of wind vector fields between ABL and free troposphere, which roughly corresponds to the ABL height, was observed (Figure 4c). Bulk Ri and virtual potential temperature profiles from RS measurements at 08:00 and 20:00 LT are presented in Figure 4d. The local gradient maximum of virtual potential temperature profiles was consistent with the results of the Ri method. The BLH from RSs at 08:00 LT was 510 m, which was similar to the 630 m determined from RWPs. However, for 08:00 LT, the BLH from RSs (650 m) was less than that determined from RWPs (990 m). The RS observation at 14:00 LT was conducted at the BJ site during summer. The BLH from RSs at 14:00 LT was 1331 m, which is consistent with the 1350 m determined from RWPs. profiles from RS measurements at 08:00 and 20:00 LT are presented in Figure 4d. The local gradient maximum of virtual potential temperature profiles was consistent with the results of the Ri method. The BLH from RSs at 08:00 LT was 510 m, which was similar to the 630 m determined from RWPs. However, for 08:00 LT, the BLH from RSs (650 m) was less than that determined from RWPs (990 m). The RS observation at 14:00 LT was conducted at the BJ site during summer. The BLH from RSs at 14:00 LT was 1331 m, which is consistent with the 1350 m determined from RWPs.

Case studies
The algorithm was also tested at the NJ site. Notably, the vertical resolution of this station was 120 m from 0 to 2 km, but shifted to 240 m above 2 km. The distribution of BLH over the NJ site was mostly less than 2 km [37,51]. Hence, the change in vertical resolution above 2 km did not significantly influence the estimation results. Figure 5 illustrates the case on 3 June 2019 at the NJ site. The BLHs from RWPs at 08:00 and 20:00 LT were 940 and 1540 m, respectively. BLH by RS measurement was calculated from the local gradient of the virtual potential temperature profile instead of bulk Ri. The BLH from RS (1010 m) was consistent with that determined from RWPs at 08:00 LT (940 m). However, for 20:00 LT, the BLH by RS (540 m) was different from the RWP-derived BLH (1540 m) due to the instability of the boundary layer structure at 20:00 LT, resulting in measurement differences. From the vertical distribution of NSDWD + NSDWS at 20:00 LT (Figure 5b), one of the local peaks (blue circle) corresponded to the BLH by RS, but due to the change in the vertical structure of the wind vector, the BLH by RWPs was located at 1540 m.  The algorithm was also tested at the NJ site. Notably, the vertical resolution of this station was 120 m from 0 to 2 km, but shifted to 240 m above 2 km. The distribution of BLH over the NJ site was mostly less than 2 km [37,51]. Hence, the change in vertical resolution above 2 km did not significantly influence the estimation results. Figure 5 illustrates the case on 3 June 2019 at the NJ site. The BLHs from RWPs at 08:00 and 20:00 LT were 940 and 1540 m, respectively. BLH by RS measurement was calculated from the local gradient of the virtual potential temperature profile instead of bulk Ri. The BLH from RS (1010 m) was consistent with that determined from RWPs at 08:00 LT (940 m). However, for 20:00 LT, the BLH by RS (540 m) was different from the RWP-derived BLH (1540 m) due to the instability of the boundary layer structure at 20:00 LT, resulting in measurement differences. From the vertical distribution of NSD WD + NSD WS at 20:00 LT (Figure 5b), one of the local peaks (blue circle) corresponded to the BLH by RS, but due to the change in the vertical structure of the wind vector, the BLH by RWPs was located at 1540 m.

Evaluation of BLH Estimates
To evaluate the BLH estimations, BLH at 08:00, 14:00, and 20:00 LT from RWPs by NSDM were compared with corresponding BLH from RS. A total of 298 sample points were used to conduct the comparison analysis. Figure 6 shows that the correlation coefficients between BLH from RWPs and RSs at 08:00, 14:00, and 20:00 LT were 0.59, 0.61, and 0.54, respectively. For all collected data, the correlation coefficient between BLH determined by NSDM and RSs was 0.66. The results indicate that the NSDM can be applied to these radar wind profilers. Moreover, the local maxima method (LMM) was employed to compare with the NSDM. The LMM defines the largest local maximum of hourly SNR profile as BLH [15]. Figure 7 shows the mean absolute difference of BLH determined using NSDM and LMM relative to BLH determined by RS. At 14:00 LT, the difference of NSDM was similar to that of LMM. At 08:00, 20:00 LT, and all cases, the difference of NSDM was less than that of LMM. Overall, the estimation results determined by NSDM were close to the estimation results of RS, and the performance of NSDM was better than LMM. However, certain errors remained. The difference between BLH determined by NSDM and RSs reached more than 1 km in some cases. RWPs rely on the vertical structure of wind fields to invert the BLH, but the RS is dependent on the thermodynamic structure to invert the BLH [22,46]. Therefore, the difference in measurement principle may be the major cause of the difference.

Evaluation of BLH Estimates
To evaluate the BLH estimations, BLH at 08:00, 14:00, and 20:00 LT from RWPs by NSDM were compared with corresponding BLH from RS. A total of 298 sample points were used to conduct the comparison analysis. Figure 6 shows that the correlation coefficients between BLH from RWPs and RSs at 08:00, 14:00, and 20:00 LT were 0.59, 0.61, and 0.54, respectively. For all collected data, the correlation coefficient between BLH determined by NSDM and RSs was 0.66. The results indicate that the NSDM can be applied to these radar wind profilers. Moreover, the local maxima method (LMM) was employed to compare with the NSDM. The LMM defines the largest local maximum of hourly SNR profile as BLH [15]. Figure 7 shows the mean absolute difference of BLH determined using NSDM and LMM relative to BLH determined by RS. At 14:00 LT, the difference of NSDM was similar to that of LMM. At 08:00, 20:00 LT, and all cases, the difference of NSDM was less than that of LMM. Overall, the estimation results determined by NSDM were close to the estimation results of RS, and the performance of NSDM was better than LMM. However, certain errors remained. The difference between BLH determined by NSDM and RSs reached more than 1 km in some cases. RWPs rely on the vertical structure of wind fields to invert the BLH, but the RS is dependent on the thermodynamic structure to invert the BLH [22,46]. Therefore, the difference in measurement principle may be the major cause of the difference.

Variations of Boundary Layer Meteorology
In this section, we analyze the variations of boundary layer meteorology, such as BLH, VC, WSBL, and surface PM2.5 concentration at the four sites using the developmental algorithm. BLH and WSBL represent the vertical and horizontal atmospheric dilution capability, and VC represents the atmospheric dilution capability. These parameters significantly affect the accumulation and diffusion of air pollutants [51,52]. Figure 8 presents diurnal variations of BLH, VC, WS, and PM2.5 concentrations at the four sites averaged over the whole time period analyzed here. BLH exhibited strong daily variation characteristics, with low values at nighttime that increased rapidly after sunrise (Figure 8a). Peak BLH values commonly appeared in the afternoon (approximately 13:00-17:00 LT). For high-altitude areas such as WQ, the hourly average BLH in the daytime was larger than that at other sites, and the maximum hourly average value reached 1.426 km. For the other sites, diurnal variations of BLH were similar. Figure 8b reports the diurnal variation of VC at four sites. For BJ, NJ, and WQ, the diurnal variation in VC was similar to that of BLH. Hence, the dilution capability was strong before sunset; gradually weakened after sunset, and remained stable at night. The daily maximum VC value of at those sites was larger than 5000 m 2 /s. However, for CQ, the variation of VC was stable and changed

Variations of Boundary Layer Meteorology
In this section, we analyze the variations of boundary layer meteorology, such as BLH, VC, WSBL, and surface PM2.5 concentration at the four sites using the developmental algorithm. BLH and WSBL represent the vertical and horizontal atmospheric dilution capability, and VC represents the atmospheric dilution capability. These parameters significantly affect the accumulation and diffusion of air pollutants [51,52]. Figure 8 presents diurnal variations of BLH, VC, WS, and PM2.5 concentrations at the four sites averaged over the whole time period analyzed here. BLH exhibited strong daily variation characteristics, with low values at nighttime that increased rapidly after sunrise (Figure 8a). Peak BLH values commonly appeared in the afternoon (approximately 13:00-17:00 LT). For high-altitude areas such as WQ, the hourly average BLH in the daytime was larger than that at other sites, and the maximum hourly average value reached 1.426 km. For the other sites, diurnal variations of BLH were similar. Figure 8b reports the diurnal variation of VC at four sites. For BJ, NJ, and WQ, the diurnal variation in VC was similar to that of BLH. Hence, the dilution capability was strong before sunset; gradually weakened after sunset, and remained stable at night. The daily maximum VC value of at those sites was larger than 5000 m 2 /s. However, for CQ, the variation of VC was stable and changed

Variations of Boundary Layer Meteorology
In this section, we analyze the variations of boundary layer meteorology, such as BLH, VC, WS BL , and surface PM 2.5 concentration at the four sites using the developmental algorithm. BLH and WS BL represent the vertical and horizontal atmospheric dilution capability, and VC represents the atmospheric dilution capability. These parameters significantly affect the accumulation and diffusion of air pollutants [51,52]. Figure 8 presents diurnal variations of BLH, VC, WS, and PM 2.5 concentrations at the four sites averaged over the whole time period analyzed here. BLH exhibited strong daily variation characteristics, with low values at nighttime that increased rapidly after sunrise (Figure 8a). Peak BLH values commonly appeared in the afternoon (approximately 13:00-17:00 LT). For high-altitude areas such as WQ, the hourly average BLH in the daytime was larger than that at other sites, and the maximum hourly average value reached 1.426 km. For the other sites, diurnal variations of BLH were similar. Figure 8b reports the diurnal variation of VC at four sites. For BJ, NJ, and WQ, the diurnal variation in VC was similar to that of BLH. Hence, the dilution capability was strong before sunset; gradually weakened after sunset, and remained stable at night. The daily maximum VC value of at those sites was larger than 5000 m 2 /s. However, for CQ, the variation of VC was stable and changed between 2000 and 3000 m 2 /s. Combined with the ABL (Figure 8c), WS could reach up to 5 m/s in plain areas such as BJ and NJ. WS could reach a relatively high value in the daytime even for WQ. However, the WS in CQ was lower than at other sites, at less than 3 m/s throughout the day. This finding was probably because the CQ site is part of a mountainous landform and the terrain blocks the circulation of the surface atmosphere, thereby reducing the atmospheric dilution capability. Overall, the atmospheric dilution capability was strong in plain areas (BJ and NJ) during the daytime, and horizontal dilution was dominant at night. In CQ, vertical dilution was dominant during the daytime, and the atmospheric dilution capability was weak at night. In WQ, vertical dilution during the daytime was the main component. between 2000 and 3000 m 2 /s. Combined with the ABL (Figure 8c), WS could reach up to 5 m/s in plain areas such as BJ and NJ. WS could reach a relatively high value in the daytime even for WQ. However, the WS in CQ was lower than at other sites, at less than 3 m/s throughout the day. This finding was probably because the CQ site is part of a mountainous landform and the terrain blocks the circulation of the surface atmosphere, thereby reducing the atmospheric dilution capability. Overall, the atmospheric dilution capability was strong in plain areas (BJ and NJ) during the daytime, and horizontal dilution was dominant at night. In CQ, vertical dilution was dominant during the daytime, and the atmospheric dilution capability was weak at night. In WQ, vertical dilution during the daytime was the main component.

Effects of Surface Parameters
In this section, the surface SHF, LHF, and PM2.5 concentrations data are collected to study their effects on the estimates of BLH. Figure 9 shows the difference in BLH from RSs and RWPs between all four sites under different surface parameters. Results indicate that the height difference decreased as PM2.5 concentration increased (Figure 9a). Similarly, Figure 9c reveals that the height difference tended to drop when SHF increased. By contrast, the height difference initially decreased and subsequently increased as the LHF increased (Figure 9b). To quantify their effects, the difference of BLH from RSs and RWPs under different conditions was calculated. Figure 9d reveals that the largest difference reached up to 0.35 ± 0.33 km when PM2.5 concentration was less than 20 μg/m 3 . The smallest BLH difference of 0.2 ± 0.16 km occurred when the PM2.5 concentration was greater than 60 μg/m 3 . By comparison, the smallest difference was 0.19 ± 0.17 km when the LHF lay within the range of 50 to 100 W/m 2 (Figure 9e) The difference of BLH became larger than 0.3 when the LHF lay in other value ranges. Interestingly, the difference tended to decline as the SHF increased. In particular, the smallest BLH difference (0.18 ± 0.17 km) occurred when the SHF was larger than 120 W/m 2 . As a side note, BLH from RSs is based on vertical thermodynamic structures, whereas BLH from RWPs mainly relies on vertical structures of wind [7,18]. It is well known that different methods and instruments will inevitably lead to differences in BLH estimates to varying degrees [53]. Given the results shown in Figure 9, the

Effects of Surface Parameters
In this section, the surface SHF, LHF, and PM 2.5 concentrations data are collected to study their effects on the estimates of BLH. Figure 9 shows the difference in BLH from RSs and RWPs between all four sites under different surface parameters. Results indicate that the height difference decreased as PM 2.5 concentration increased (Figure 9a). Similarly, Figure 9c reveals that the height difference tended to drop when SHF increased. By contrast, the height difference initially decreased and subsequently increased as the LHF increased (Figure 9b). To quantify their effects, the difference of BLH from RSs and RWPs under different conditions was calculated. Figure 9d reveals that the largest difference reached up to 0.35 ± 0.33 km when PM 2.5 concentration was less than 20 µg/m 3 . The smallest BLH difference of 0.2 ± 0.16 km occurred when the PM 2.5 concentration was greater than 60 µg/m 3 . By comparison, the smallest difference was 0.19 ± 0.17 km when the LHF lay within the range of 50 to 100 W/m 2 (Figure 9e). The difference of BLH became larger than 0.3 when the LHF lay in other value ranges. Interestingly, the difference tended to decline as the SHF increased. In particular, the smallest BLH difference (0.18 ± 0.17 km) occurred when the SHF was larger than 120 W/m 2 . As a side note, BLH from RSs is based on vertical thermodynamic structures, whereas BLH from RWPs mainly relies on vertical structures of wind [7,18]. It is well known that different methods and instruments will inevitably lead to differences in BLH estimates to varying degrees [53]. Given the results shown in Figure 9, the difference between RSs and RWPs of less than 0.3 km can be considered acceptable. However, under extreme conditions such as low PM 2.5 and SHF, or extreme high or low LHF, the difference of BLH seems to be too large, and the influential mechanism merits further analysis based on model simulation in the future.
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 17 difference between RSs and RWPs of less than 0.3 km can be considered acceptable. However, under extreme conditions such as low PM2.5 and SHF, or extreme high or low LHF, the difference of BLH seems to be too large, and the influential mechanism merits further analysis based on model simulation in the future. To understand the effect of these parameters, the vertical and thermodynamic structures of wind fields under different atmospheric conditions were investigated. Figure 10 demonstrates the vertical distribution of wind fields and virtual potential temperature (PTv) under different surface PM2.5, LHF, and SHF. Figure 10a indicates that the wind shear and PTv gradient profiles became more evident as PM2.5 increased. This result is similar to those in previous studies; polluted days were usually accompanied by significant atmospheric stratification [54]. As shown in Figure 9, the difference of BLH from RSs and RWPs was insignificant when the vertical wind profile exhibited strong wind shear and the PTv profile had a strong gradient. Similarly, Figure 10b shows that the LHF was within the range of 50-100 W/m 2 , wind shear was evident, and certain gradient change in the PTv profile existed. This condition had the smallest mean absolute deviation. A similar situation was observed for SHF (Figure 10c). The appearance of wind shear was consistent with that of the gradient in the PTv profile when the SHF was larger than 80 W/m 2 , thereby leading to a small mean absolute deviation. Therefore, the BLH determined from RWPs was only consistent with that from RSs when wind shear and the PTv profile gradient existed simultaneously. This finding is explained given that the surface parameters affect the vertical and thermodynamic structures of wind fields. The RWP relies on wind shear to invert the BLH. However, the RS is dependent on the thermodynamic structure to invert the BLH. Therefore, the height difference between the BLH from RSs and RWPs is only evident when the vertical field or thermodynamic structures of the wind are unstable. The results of the effects of surface parameters on atmospheric vertical structure provide new insights for the application of RWP BLH observations. To understand the effect of these parameters, the vertical and thermodynamic structures of wind fields under different atmospheric conditions were investigated. Figure 10 demonstrates the vertical distribution of wind fields and virtual potential temperature (PT v ) under different surface PM 2.5 , LHF, and SHF. Figure 10a indicates that the wind shear and PT v gradient profiles became more evident as PM 2.5 increased. This result is similar to those in previous studies; polluted days were usually accompanied by significant atmospheric stratification [54]. As shown in Figure 9, the difference of BLH from RSs and RWPs was insignificant when the vertical wind profile exhibited strong wind shear and the PT v profile had a strong gradient. Similarly, Figure 10b shows that the LHF was within the range of 50-100 W/m 2 , wind shear was evident, and certain gradient change in the PT v profile existed. This condition had the smallest mean absolute deviation. A similar situation was observed for SHF ( Figure 10c). The appearance of wind shear was consistent with that of the gradient in the PT v profile when the SHF was larger than 80 W/m 2 , thereby leading to a small mean absolute deviation. Therefore, the BLH determined from RWPs was only consistent with that from RSs when wind shear and the PT v profile gradient existed simultaneously. This finding is explained given that the surface parameters affect the vertical and thermodynamic structures of wind fields. The RWP relies on wind shear to invert the BLH. However, the RS is dependent on the thermodynamic structure to invert the BLH. Therefore, the height difference between the BLH from RSs and RWPs is only evident when the vertical field or thermodynamic structures of the wind are unstable. The results of the effects of surface parameters on atmospheric vertical structure provide new insights for the application of RWP BLH observations. Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 17 Figure 10. Vertical distribution of wind (black arrows) and virtual potential temperature (PTv, in red) shown as a function of (a) surface PM2.5 concentration, (b) LHF, and (c) SHF. Blue lines represent root mean squared difference (RMSD) of the difference of BLH between RWP and RS measurements at each bin.

Conclusion
Using continuous high-temporal resolution wind measurements from RWPs in BJ, NJ, CQ, and WQ, we developed an adaptive method to estimate BLH from the vertical structure of wind fields. This algorithm is advantageous over previously developed algorithms, since it can calculate BLH without assuming initial parameters such as threshold and window function.
The proposed method was applied to continuous RWP measurements at the BJ, NJ, CQ, and WQ sites during the summer. To assess the accuracy of this method, a comprehensive comparison was conducted between BLH determined from RSs and from RWPs. The results indicate that the correlation coefficient between BLH retrievals from RWPs and from RSs was 0.66. However, in some cases, the difference between BLH from RWPs and from RSs reached more than 1 km. To understand the differences, the effects of surface SHF, LHF, and PM2.5 concentration on the accuracy of BLH from RWPs were subsequently investigated. The results show that their height difference decreased as PM2.5 concentration increased, while the same phenomenon was observed for SHF. The height difference initially decreased and subsequently increased as the LHF increased. Moreover, the diurnal variation of summertime boundary layer meteorology at the four sites was examined. Results Figure 10. Vertical distribution of wind (black arrows) and virtual potential temperature (PTv, in red) shown as a function of (a) surface PM 2.5 concentration, (b) LHF, and (c) SHF. Blue lines represent root mean squared difference (RMSD) of the difference of BLH between RWP and RS measurements at each bin.

Conclusions
Using continuous high-temporal resolution wind measurements from RWPs in BJ, NJ, CQ, and WQ, we developed an adaptive method to estimate BLH from the vertical structure of wind fields. This algorithm is advantageous over previously developed algorithms, since it can calculate BLH without assuming initial parameters such as threshold and window function.
The proposed method was applied to continuous RWP measurements at the BJ, NJ, CQ, and WQ sites during the summer. To assess the accuracy of this method, a comprehensive comparison was conducted between BLH determined from RSs and from RWPs. The results indicate that the correlation coefficient between BLH retrievals from RWPs and from RSs was 0.66. However, in some cases, the difference between BLH from RWPs and from RSs reached more than 1 km. To understand the differences, the effects of surface SHF, LHF, and PM 2.5 concentration on the accuracy of BLH from RWPs were subsequently investigated. The results show that their height difference decreased as PM 2.5 concentration increased, while the same phenomenon was observed for SHF. The height difference initially decreased and subsequently increased as the LHF increased. Moreover, the diurnal variation of summertime boundary layer meteorology at the four sites was examined. Results revealed that the boundary layer exhibited clear diurnal variations, and the atmospheric dilution capability had significant regional differences. For the BJ, NJ, and WQ sites, the dilution capability was strong before sunset, gradually weakened after sunset, and remained stable at night. The maximum hourly mean VC at BJ, NJ, and WQ was larger than 5000 m 2 /s. However, for CQ, the variation of hourly mean VC was stable, changing between 2000 and 3000 m 2 /s.
We propose an adaptive method in an attempt to estimate BLH from high-resolution RWP measurements in BJ, NJ, CQ, and WQ. The proposed algorithm can provide reasonable estimates of BLH from the vertical structure of wind fields. The findings of this study provide key references for ABL-related research using large-scale RWP measurements. However, this algorithm cannot be applied when the vertical structure of the wind field is disordered and when either multiple or no strong wind shear would cause uncertain errors. Therefore, much work will be needed in the future to improve the algorithm. For example, we can try to improve the performance of the algorithm by combining deep learning methods. Afterward, efforts will be devoted to estimating BLH using the vast wealth of height-revolved wind from the well-established RWP network in China.

Author Contributions:
The study was completed with cooperation among all authors. J.G. and B.L. designed this study; B.L., J.G., and W.G. conducted the experiment and wrote the manuscript with input from all coauthors; Y.S. and S.J. checked the experimental results. All authors have read and agreed to the published version of the manuscript.