Effects of Boundary Layer Height on the Model of Ground-Level PM2.5 Concentrations from AOD: Comparison of Stable and Convective Boundary Layer Heights from Different Methods

The aerosol optical depth (AOD) from satellites or ground-based sun photometer spectral observations has been widely used to estimate ground-level PM2.5 concentrations by regression methods. The boundary layer height (BLH) is a popular factor in the regression model of AOD and PM2.5, but its effect is often uncertain. This may result from the structures between the stable and convective BLHs and from the calculation methods of the BLH. In this study, the boundary layer is divided into two types of stable and convective boundary layer, and the BLH is calculated using different methods from radiosonde data and National Centers for Environmental Prediction (NCEP) reanalysis data for the station in Beijing, China during 2014–2015. The BLH values from these methods show significant differences for both the stable and convective boundary layer. Then, these BLHs were introduced into the regression model of AOD-PM2.5 to seek the respective optimal BLH for the two types of boundary layer. It was found that the optimal BLH for the stable boundary layer is determined using the method of surface-based inversion, and the optimal BLH for the convective layer is determined using the method of elevated inversion. Finally, the optimal BLH and other meteorological parameters were combined to predict the PM2.5 concentrations using the stepwise regression method. The results indicate that for the stable boundary layer, the optimal stepwise regression model includes the factors of surface relative humidity, BLH, and surface temperature. These three factors can significantly enhance the prediction accuracy of ground-level PM2.5 concentrations, with an increase of determination coefficient from 0.50 to 0.68. For the convective boundary layer, however, the optimal stepwise regression model includes the factors of BLH and surface wind speed. These two factors improve the determination coefficient, with a relatively low increase from 0.65 to 0.70. It is found that the regression coefficients of the BLH are positive and negative in the stable and convective regression models, respectively. Moreover, the effects of meteorological factors are indeed related to the types of BLHs.


Introduction
Air pollution seriously affects the environment and endangers human health.Smaller particles with aerodynamic diameters less than 2.5 µm (also known as PM 2.5 ) are dominant in industrial and urban pollution and can increase the incidence of heart disease, cardiovascular disease and Atmosphere 2017, 8, 104; doi:10.3390/atmos8060104www.mdpi.com/journal/atmospherelung cancer [1][2][3][4].Most aerosols are emitted into the atmospheric boundary layer and result in serious pollution near the ground.However, the space coverage of PM 2.5 monitoring stations is sparse, especially in the suburban environment.Aerosol optical depth (AOD) data from satellites or ground-based sun photometer spectral observations can play an auxiliary role in air quality monitoring.AOD is the integral of the aerosol extinction due to scattering and absorption in the vertical plane, and surface PM 2.5 concentrations are related to the aerosol extinction.There is a wide body of literature that shows that there are strong correlations between AOD and surface PM 2.5 concentrations [5][6][7].
The regression method is widely employed for estimating ground-level PM 2.5 concentrations.Some meteorological parameters such as temperature, relative humidity, and wind are introduced into the regression model of AOD-PM 2.5 to improve the accuracy of estimating PM 2.5 concentrations [8][9][10].Boundary layer height (BLH) is a popular parameter in the regression model, because it can indicate the height of turbulent diffusion.A higher BLH implies stronger turbulent diffusion, and usually results in a lower surface PM 2.5 concentration [11,12].
However, some studies have suggested that the effect of BLH is uncertain or insignificant to the regression model of AOD-PM 2.5 [8,13,14].Tian and Chen [8] investigated the effects of the BLH and meteorological factors on the regression model of MODIS AOD in southern Ontario, Canada during 2004.They found that the P value of the BLH is 0.71, which is far larger than 0.01.It is suggested that the BLH is an ineffective factor for the regression model.Liu et al. [13] reported that the BLH is not significant at the α = 0.05 level in the model of Multiangle Imaging Spector Radiometer (MISR) AOD in the St Louis, Missouri, USA.One possible reason for the failure of the BLH in the model is that there are considerable calculation errors, especially for the BLH from the model product.Even using the radiosonde data, it is difficult to estimate an actual BLH, because the structure of the actual boundary layer is complex and varies with different definition and calculation methods [15].
Generally, the boundary layer is classified into the stable boundary layer and the convective boundary layer [16,17].The stable boundary layer height is commonly defined as the height where the negative buoyancy flux at the surface damps the turbulence.It is always associated with a surface-based temperature inversion, and can be estimated by the top of the inversion using the temperature profile [18][19][20].The convective boundary layer height is commonly defined as the height where the positive buoyancy flux at the surface creates a thermal instability.There are several methods to estimate convective BLH using the temperature, humidity, and refractivity profiles of sounding data, and the values of BLH calculated from these methods are usually different [17,21,22].In addition, BLHs can also be estimated from a model, but the values depend on the model and its boundary layer parameterization.These differences may affect the application and assessment of BLH in the regression model of AOD-PM 2. 5 .
In this study, the BLHs calculated from radiosonde data and NCEP reanalysis data were compared for the station in Beijing, China during 2014-2015.Then, these BLHs were respectively introduced into the regression model to seek the optimal stable and convective BLHs.Finally, we used the optimal BLH and meteorological parameters to improve the prediction of ground-level PM 2.5 concentrations.In this paper, the materials and methodology are described in Sections 2 and 3, respectively.A comparison of BLHs and their effect on the AOD-PM 2.5 regression equation are analyzed and discussed in Section 4. Several major conclusions and potential future improvements to the model of AOD-PM 2.5 are summarized in Section 5.

Materials
To assess the effects of BLHs in the regression model of AOD-PM 2.5 , four types of data sets were collected, including hourly PM 2.5 concentrations, Aerosol Robotic Network (AERONET) AOD, NCEP reanalysis data, and radiosonde data.These data sets were collected in the Beijing area.However, the specific sites of stations for these data sets vary (Figure 1).The distances between the PM 2.5 , AERONET, and radiosonde stations are less than 10 km, and the center resides at the AERONET station.The NCEP grid data were also interpolated with respect to the AERONET station.Sections 2.1-2.4 describe each data set in more detail.

Radiosonde Data
The radiosonde measurements are obtained from the National Climate Data Center's (NCDC) Integrated Global Radiosonde Archive (IGRA) [17,23].The latest version of the NCDC IGRA provides quality-assured daily radiosonde records at mandatory pressure levels, additionally required levels, and thermodynamically significant levels [24,25].The records of these levels can describe the detailed construction of the troposphere, and are widely used to calculate the BLH [18,20,26].The radiosonde records consist of temperature, pressure, relative humidity, wind direction, and wind speed.Typically, radiosondes are launched two times per day at 00:00 UTC and 12:00 UTC.In this study, the radiosonde at 00:00 UTC from 1 January, 2014 to 31 December, 2015 were used and matched with other data.

Ground-Based PM2.5 Concentration Data
The PM2.5 concentration data are collected from the Ministry of Environmental Protection (MEP) of China.The concentration is measured hourly by a Tapered Element Oscillating Microbalance (TEOM) with an accuracy of 1.5 µg/m 3 .There are 12 PM2.5 monitoring stations in the Beijing area.The Guan-Yuan station is the nearest PM2.5 monitoring station to the radiosonde station, which is approximately 10 km away (Figure 1).Since the radiosonde is launched once per day at 00:00 UTC, the hourly PM2.5 concentration at 00:00 UTC is selected.Invalid data (values recorded as NAN) because of equipment breakdown are removed.

AERONET AOD Data
The AERONET is a globally distributed network of sun photometers that provides multiwavelength (340, 380, 440, 500, 675, 870, 940, and 1020 nm) AOD measurements [27][28][29].The surface sun photometer employed by the AERONET has a very narrow field of view, and is therefore rarely affected by surface reflectance and aerosol forward scattering [30].The AOD-retrieval uncertainty of AERONET AOD is only 0.01-0.02,and is widely used to validate satellite AOD values [31][32][33].There are five AERONET stations in the Beijing area.The Beijing-CAMS AERONET station is located 5.0 km from the closest radiosonde station (Figure 1) and can provide a complete data record during

Radiosonde Data
The radiosonde measurements are obtained from the National Climate Data Center's (NCDC) Integrated Global Radiosonde Archive (IGRA) [17,23].The latest version of the NCDC IGRA provides quality-assured daily radiosonde records at mandatory pressure levels, additionally required levels, and thermodynamically significant levels [24,25].The records of these levels can describe the detailed construction of the troposphere, and are widely used to calculate the BLH [18,20,26].The radiosonde records consist of temperature, pressure, relative humidity, wind direction, and wind speed.Typically, radiosondes are launched two times per day at 00:00 UTC and 12:00 UTC.In this study, the radiosonde at 00:00 UTC from 1 January 2014 to 31 December 2015 were used and matched with other data.

Ground-Based PM 2.5 Concentration Data
The PM 2.5 concentration data are collected from the Ministry of Environmental Protection (MEP) of China.The concentration is measured hourly by a Tapered Element Oscillating Microbalance (TEOM) with an accuracy of 1.5 µg/m 3 .There are 12 PM 2.5 monitoring stations in the Beijing area.The Guan-Yuan station is the nearest PM 2.5 monitoring station to the radiosonde station, which is approximately 10 km away (Figure 1).Since the radiosonde is launched once per day at 00:00 UTC, the hourly PM 2.5 concentration at 00:00 UTC is selected.Invalid data (values recorded as NAN) because of equipment breakdown are removed.

AERONET AOD Data
The AERONET is a globally distributed network of sun photometers that provides multi-wavelength (340, 380, 440, 500, 675, 870, 940, and 1020 nm) AOD measurements [27][28][29].The surface sun photometer employed by the AERONET has a very narrow field of view, and is therefore rarely affected by surface reflectance and aerosol forward scattering [30].The AOD-retrieval uncertainty of AERONET AOD is only 0.01-0.02,and is widely used to validate satellite AOD values [31][32][33].There are five AERONET stations in the Beijing area.The Beijing-CAMS AERONET station is located 5.0 km from the closest radiosonde station (Figure 1) and can provide a complete data record during 2014-2015.AERONET AOD data are computed for three data quality levels: level 1.0 (unscreened), level 1.5 (cloud-screened), and level 2.0 (cloud-screened and quality-assured).Data processing, cloud-screening algorithm, and inversion techniques are described by Holben et al. (1998Holben et al. ( , 2001)).For the Beijing-CAMS station, only level 1.0 and level 1.5 data were available for the dates of interest.The data of level 1.5 are used widely in the previous studies [34][35][36].In this study, we used the AERONET Level 1.5 (cloud-screened) aerosol product from the Beijing-CAMS station at a wavelength of 500 nm.The radiosonde is at 00:00 UTC, but the AERONET observations are usually not taken exactly at 00:00 UTC.There is a small interval of a few minutes near 00:00 UTC.The averaged AOD from 23:30 to 00:30 UTC was matched with the radiosonde at 00:00 UTC.In addition, if the standard deviation of the average AOD was greater than 0.5, the average AOD was removed to reduce the possibility of spurious AOD values [13,37].

BLH Data From NCEP
In addition to the BLH calculated from the sounding profile, we also used the BLH from the NCEP's global reanalysis data, with a horizontal resolution of 1 • × 1 • and a time step of 6 h (00:00, 06:00, 12:00, and 18:00 UTC).The NCEP uses the bulk Richardson number (Ri) approach to iteratively estimate BLH starting from the ground upward [38,39].This approach is suitable for both stable and convective boundary layers [26].More detailed information of this reanalysis product can be found on the NOAA website [40].In this study, only the 00:00 UTC BLH during 2014-2015 are collected.Note that the NCEP reanalysis data are on a regular latitude/longitude grid.The grid-cell BLHs were interpolated to the AERONET station using the inverse distance weighting method.

The Method Used to Estimate the BLH
The BLH is divided into two types: the stable BLH and convective BLH.If there is a surface-based inversion in the sounding profile, the BLH is classified as the stable BLH (BLH Sta ) [18,19].It is usually estimated by the top of the surface-based inversion (SBI), written as BLH Sta  SBI .If there is not a surface-based inversion in the sounding profile, the BLH is classified as the convective BLH (BLH Con ).It can be calculated by the following five traditional methods: (1) The "parcel method" evaluates the BLH using the profile of virtual potential temperature (θ v ).The BLH is the height of the value of surface θ v equaling the value of aloft θ v [17,41].The θ v is calculated by where T is the atmospheric temperature, q is the specific humidity, P is the atmospheric pressure, R d is the gas constant for dry air, and c pd is the constant pressure specific heat for dry air.The BLH evaluated by the θ v is written as BLH Con θ v .(2) The height of the maximum vertical gradient of potential temperature (θ) [42,43], written as BLH Con θ .The θ is calculated by where T, P, R d , and c pd are the same as the variables included in the definition of θ v .
(3) The base of an elevated inversion (EI) [18].The height is estimated by the elevated inversion of the temperature profile and written as BLH Con  EI .(4) The height of the minimum vertical gradient of relative humidity (RH) [44].The BLH is estimated using the RH profile and written as BLH Con  RH .
(5) The height of the minimum vertical gradient of refractivity (N) [21,45], written as BLH Con N .The refractivity is calculated by [46,47] where N is the refractivity, P is the atmospheric pressure, T is the atmospheric temperature, and e is water vapor pressure.
Table 1 provides the basic information of these calculation methods from the sounding profile data.For all these methods, we restricted the available data of all the sounding records to 4000 m to avoid mistaking free tropospheric features for the top of the boundary layer [18,47].Moreover, to avoid spurious estimates of large vertical gradients resulting from horizontal (or vertical) separation of the surface instrument shelter from the radiosonde launch site, the height of the first level in the sounding was taken as the surface level [20].
The BLH from the NCEP reanalysis data is written as BLH RE .Because the BLH RE is based on the boundary layer parameterization of the model, it is not assigned to the stable or convective types in the model.We classified the BLH RE into stable or convective types also according to the surface-based inversion from the sounding profile data.Namely, if there is a surface-based inversion from sounding profile data, the BLH RE is classified as stable BLH, written as BLH Sta RE .If not, the BLH RE is classified as convective BLH, written as BLH Con RE .The effect of the BLH was assessed using the regression model to predict PM 2.5 concentrations based on AOD.First, a simple lognormal regression model (termed as M-I) was developed to predict the surface PM 2.5 concentrations using AERONET AOD as the only factor.
where PM 2.5 is the surface PM 2.5 concentration (µg/m 3 ) and AOD is the AERONET AOD (unitless).a 1 is the intercept and a 2 is the slope for the M-I.
Then, the BLH from different methods was added into the M-I to form a new lognormal model (termed as M-II). ln The regression coefficients b 2 and b 3 are associated with predictor factors of AOD and BLH (m), respectively.By using the M-II, we expect to find the optimal BLHs out of the two stable BLHs and six convective BLHs noted in Table 1.
Finally, a semi-empirical model (termed as M-III) with more surface meteorological factors was proposed for improving the prediction of surface PM 2.5 concentrations [13,48].A stepwise regression method was used in the M-III, and the factors of AOD, optimal BLH, and surface meteorological parameters were added into the model one by one.If all factors are introduced into the model, the regression model is given by the following form: where the regression coefficients c 2 -c 6 are associated with AERONET AOD, BLH(m), surface temperature ( • C), surface relative humidity (%) and surface wind speed (m/s), respectively.Certainly, the optimal stepwise regression may only include parts of the previously mentioned factors.
The predicted PM 2.5 concentrations were fitted against the observed values to evaluate the performance of the regression model.In addition, the determination coefficient (R 2 ) and the Root Mean Square Error (RMSE) were calculated to evaluate the degree of fit between predicted and observed PM 2.5 concentrations.

Seasonal Difference of Stable and Convective Boundary Layer
The samples of sounding profiles at 00 UTC during 2014-2015 are divided into stable or convective boundary layers according to whether or not there is a surface-based inversion.There are 314 samples for the stable boundary layer, and 416 samples for the convective boundary layers (Table 2).These samples are further analyzed for four seasons.The four seasons are defined as follows.Spring includes March, April and May; summer includes June, July and August; autumn includes September, October and November; and winter includes December, January and February.The percentage of the stable boundary layer is highest for winter with 61.67% that is almost two times as large as the percentage of the convective boundary layer (38.33%).If the performed time of the sounding profiles is at 00:00 UTC that corresponds to 08:00 local time in the morning.It can be inferred that the sunrise time is late in winter, and the surface temperature is still low.Thus, most stable boundary layers that formed during the night have not been destroyed.On the contrary, the percentage of stable boundary layer for the summer is only 21.74%, since most stable boundary layers are destroyed at 08:00 local time by the quickly increasing surface temperature after the sunrise.Figure 2 shows the average vertical temperature profiles of stable boundary layer and convective boundary layer for different seasons.There is significant surface-based inversion for the stable boundary layer (red line) in the spring, autumn and winter.Moreover, the height of the inversion layer is lower in the winter than that in other seasons.This may be an important reason for the serious aerosol pollution in the winter, since the aerosol is inhibited in the lower inversion layer.However, for the summer, the temperature profile (red line) displays a weak and shallow inversion layer.For the convective layer (blue line), the tendencies of temperature all decrease with height in four seasons.
There is a weak isothermal layer between 0.25 and 0.5 km for the spring (Figure 2a).This probably resulted from elevated inversion layers.However, since the height of elevated inversions is various, there are is not a significant elevated inversion for the average profile.
Atmosphere 2017, 8, 104 7 of 17 the convective layer (blue line), the tendencies of temperature all decrease with height in four seasons.
There is a weak isothermal layer between 0.25 and 0.5 km for the spring (Figure 2a).This probably resulted from elevated inversion layers.However, since the height of elevated inversions is various, there are is not a significant elevated inversion for the average profile.are in summer and spring, since the sunrise is early in the morning for these two seasons.The surface temperature and θ quickly increase after sunrise, and the boundary layer will transform from stable to convective.However, the boundary layer is more likely to be neutral in the seasons of autumn and winter, and there is not a θ on the aloft profile equal to the surface θ for these samples with neutral stratification.We can increase the sample sizes of BLH by adding an excess δθ on the surface θ [17].However, it is difficult to determine the δθ using the radiosonde data for its coarse vertical levels.In addition, the sample size of BLH is also less than the sample sizes of other convective BLHs, because there are a few samples without elevated inversion in the total samples of the convective boundary layer [18].θ v are in summer and spring, since the sunrise is early in the morning for these two seasons.The surface temperature and θ v quickly increase after sunrise, and the boundary layer will transform from stable to convective.However, the boundary layer is more likely to be neutral in the seasons of autumn and winter, and there is not a θ v on the aloft profile equal to the surface θ v for these samples with neutral stratification.We can increase the sample sizes of BLH Con θ v by adding an excess δθ v on the surface θ v [17].However, it is difficult to determine the δθ v using the radiosonde data for its coarse vertical levels.In addition, the sample size of BLH Con EI is also less than the sample sizes of other convective BLHs, because there are a few samples without elevated inversion in the total samples of the convective boundary layer [18].In Figure 3, the stable type of BLH (two bars on the left) is significantly lower than the convective type of BLH (six bars on the right).This is because the stable boundary layer with negative buoyancy flux suppresses turbulence and results in a lower BLH, whereas the convective boundary layer with positive buoyancy flux enhances turbulence and results in a higher BLH [15,18].Although the BLHs of the convective type are higher than the BLHs of the stable type, the BLHs within the same type also show significant differences.For the stable type of BLH, the BLH is lower than the BLH by about 117 m.Moreover, there is a negative correlation coefficient (R = −0.29) between these two BLHs (Figure 4).

Statistical Comparison among Different BLH Methods
Because the BLH is estimated from the sounding profile data, it should more closely approximate to the actual situation of BLH.Thus, it is concluded that the BLH derived from the reanalysis products may be inaccurate, due to the horizontal and temporal resolution of the model, as well as the difference in BLH definition in the model [49][50][51].It should be noted that the method of estimating BLH in the model is not based on the temperature inversion, which is different than the method of estimating BLH .Thus, the correlation between BLH and BLH is low.For the convective type of BLH, the correlations between the six BLHs vary greatly, but most of them are low (Figure 5).Except for the squares along the diagonal line, which indicate autocorrelation of 1.0, the highest correlation is between BLH and BLH at about 0.66.A possible explanation for this is that these two BLHs are estimated using the temperature profile.The second highest correlation is between BLH and BLH at 0.63, and the correlation between BLH and BLH also reaches 0.35.It seems that the BLH is relatively reliable for the convective boundary layer compared to BLH for the stable boundary layer.A possible reason for this is that the In Figure 3, the stable type of BLH (two bars on the left) is significantly lower than the convective type of BLH (six bars on the right).This is because the stable boundary layer with negative buoyancy flux suppresses turbulence and results in a lower BLH, whereas the convective boundary layer with positive buoyancy flux enhances turbulence and results in a higher BLH [15,18].Although the BLHs of the convective type are higher than the BLHs of the stable type, the BLHs within the same type also show significant differences.For the stable type of BLH, the BLH Sta RE is lower than the BLH Sta SBI by about 117 m.Moreover, there is a negative correlation coefficient (R = −0.29) between these two BLHs (Figure 4).
Because the BLH Sta SBI is estimated from the sounding profile data, it should more closely approximate to the actual situation of BLH.Thus, it is concluded that the BLH Sta RE derived from the reanalysis products may be inaccurate, due to the horizontal and temporal resolution of the model, as well as the difference in BLH definition in the model [49][50][51].It should be noted that the method of estimating BLH in the model is not based on the temperature inversion, which is different than the method of estimating BLH Sta SBI .Thus, the correlation between BLH Sta SBI and BLH Sta RE is low.
Atmosphere 2017, 8, 104 8 of 17 In Figure 3, the stable type of BLH (two bars on the left) is significantly lower than the convective type of BLH (six bars on the right).This is because the stable boundary layer with negative buoyancy flux suppresses turbulence and results in a lower BLH, whereas the convective boundary layer with positive buoyancy flux enhances turbulence and results in a higher BLH [15,18].Although the BLHs of the convective type are higher than the BLHs of the stable type, the BLHs within the same type also show significant differences.For the stable type of BLH, the BLH is lower than the BLH by about 117 m.Moreover, there is a negative correlation coefficient (R = −0.29) between these two BLHs (Figure 4).
Because the BLH is estimated from the sounding profile data, it should more closely approximate to the actual situation of BLH.Thus, it is concluded that the BLH derived from the reanalysis products may be inaccurate, due to the horizontal and temporal resolution of the model, as well as the difference in BLH definition in the model [49][50][51].It should be noted that the method of estimating BLH in the model is not based on the temperature inversion, which is different than the method of estimating BLH .Thus, the correlation between BLH and BLH is low.For the convective type of BLH, the correlations between the six BLHs vary greatly, but most of them are low (Figure 5).Except for the squares along the diagonal line, which indicate autocorrelation of 1.0, the highest correlation is between BLH and BLH at about 0.66.A possible explanation for this is that these two BLHs are estimated using the temperature profile.The second highest correlation is between BLH and BLH at 0.63, and the correlation between BLH and BLH also reaches 0.35.It seems that the BLH is relatively reliable for the convective boundary layer compared to BLH for the stable boundary layer.A possible reason for this is that the For the convective type of BLH, the correlations between the six BLHs vary greatly, but most of them are low (Figure 5).Except for the squares along the diagonal line, which indicate auto-correlation of 1.0, the highest correlation is between BLH Con θ and BLH Con EI at about 0.66.A possible explanation for this is that these two BLHs are estimated using the temperature profile.The second highest correlation is between BLH Con θ v and BLH Con RE at 0.63, and the correlation between BLH Con RE and BLH Con EI also reaches Atmosphere 2017, 8, 104 9 of 18 0.35.It seems that the BLH Con RE is relatively reliable for the convective boundary layer compared to BLH Sta RE for the stable boundary layer.A possible reason for this is that the variation within the convective BLH is large, and the tendency of BLH from the model usually agrees with that from the sounding profiles.In Figure 5, the lowest correlation is between BLH Con θ v and BLH Con N at just about 0.01, which indicates that these two BLHs are approximately independent.In addition, the correlations between BLH Con RH and BLH Con RE , BLH Con N and BLH Con RE , BLH Con θ v and BLH Con θ are also less than 0.1.Thus, these BLHs will result in different effects on the regression model of AOD-PM 2.5 .
Atmosphere 2017, 8, 104 9 of 17 variation within the convective BLH is large, and the tendency of BLH from the model usually agrees with that from the sounding profiles.In Figure 5, the lowest correlation is between BLH and BLH at just about 0.01, which indicates that these two BLHs are approximately independent.In addition, the correlations between BLH and BLH , BLH and BLH , BLH and BLH are also less than 0.1.Thus, these BLHs will result in different effects on the regression model of AOD-PM2.5.  Figure 6a shows the scatter plot of observed vs. predicted surface PM2.5 concentrations derived from the M-I-Sta.This model can explain 50% (R 2 = 0.50) of the variability in the corresponding PM2.5 concentrations.The RMSE of this regression model is 44.30µg/m 3 .For the M-I-Con (Figure 6b), the R 2 is 0.65 and the RMSE is 31.76µg/m 3 , which is superior to the results of the M-I-Sta.It is found that there are some points with a large observed concentration but small predicted concentration that decrease the accuracy of the M-I-Sta (Figure 6a).We reasoned that while the observed PM2.5 may be large, it may exhibit a low AOD in the M-I-Sta, because the surface-based inversion can inhibit the turbulent diffusion of aerosols by acting like a lid [52].However, aerosols are relatively uniformly mixed in the convective boundary layer.Thus, the prediction accuracy of the M-I-Con is higher than that of the M-I-Sta.

Performances of Simple Lognormal Models
A simple lognormal model (M-I) was developed to describe the relationship between surface PM 2.5 concentrations and AOD.To analyze the performance of the model for the two types of boundary layers, the M-I was divided into two types for the stable boundary layer (M-I-Sta) and the convective boundary layer (M-I-Con).There are 342 samples after matching the AOD and PM 2.5 data during 2014-2015.These samples are also divided into two parts for the two types of models.There are 156 and 186 samples for the M-I-Sta and M-I-Con after matching, respectively.
Figure 6a shows the scatter plot of observed vs. predicted surface PM 2.5 concentrations derived from the M-I-Sta.This model can explain 50% (R 2 = 0.50) of the variability in the corresponding PM 2.5 concentrations.The RMSE of this regression model is 44.30µg/m 3 .For the M-I-Con (Figure 6b), the R 2 is 0.65 and the RMSE is 31.76µg/m 3 , which is superior to the results of the M-I-Sta.It is found that there are some points with a large observed concentration but small predicted concentration that decrease the accuracy of the M-I-Sta (Figure 6a).We reasoned that while the observed PM 2.5 may be large, it may exhibit a low AOD in the M-I-Sta, because the surface-based inversion can inhibit the turbulent diffusion of aerosols by acting like a lid [52].However, aerosols are relatively uniformly mixed in the convective boundary layer.Thus, the prediction accuracy of the M-I-Con is higher than that of the M-I-Sta.

Ability of Different Blhs to Estimate Surface PM2.5 Concentrations
In this section, the BLH values from different methods were added into the M-I to develop a new model, M-II, according to Equation ( 5).The M-II was also developed for the stable BLH (M-II-Sta) and convective BLH (M-II-Con), respectively.We expect to find the optimal BLHs out of the two stable BLHs and six convective BLHs.The samples of the M-II-Sta and the M-II-Con should be matched with the BLH based on the samples of the M-I-Sta and the M-I-Con, respectively.There are 156 samples for the M-II-Sta.For the M-II-Con, however, the sample sizes vary with the different convective BLHs.There are 60 and 149 samples after matching with BLH and BLH , respectively.There are 186 samples after matching with BLH , BLH , BLH and BLH .Figure 7 shows the scatter plots of the observed versus the predicted surface PM2.5 concentrations for the results of the M-II-Sta.The R 2 of the M-II-Sta with BLH (R 2 = 0.59, Figure 7a) is significantly higher than the R 2 of the M-I-Sta (R 2 = 0.50, Figure 5a), but the R 2 of the M-II-Sta with BLH (R 2 = 0.51, Figure 7b) is similar to the R 2 of the M-I-Sta.The results indicate that BLH but not BLH can effectively improve the model performance for the prediction of PM2.5.Moreover, the regression coefficient associated with BLH is 0.38 (Figure 7a), but the regression coefficient associated with

Ability of Different Blhs to Estimate Surface PM 2.5 Concentrations
In this section, the BLH values from different methods were added into the M-I to develop a new model, M-II, according to Equation ( 5).The M-II was also developed for the stable BLH (M-II-Sta) and convective BLH (M-II-Con), respectively.We expect to find the optimal BLHs out of the two stable BLHs and six convective BLHs.The samples of the M-II-Sta and the M-II-Con should be matched with the BLH based on the samples of the M-I-Sta and the M-I-Con, respectively.There are 156 samples for the M-II-Sta.For the M-II-Con, however, the sample sizes vary with the different convective BLHs.There are 60 and 149 samples after matching with BLH Con θ v and BLH Con EI , respectively.There are 186 samples after matching with BLH Con θ , BLH Con RH , BLH Con N and BLH Con RE .Figure 7 shows the scatter plots of the observed versus the predicted surface PM 2.5 concentrations for the results of the M-II-Sta.The R 2 of the M-II-Sta with BLH Sta SBI (R 2 = 0.59, Figure 7a) is significantly higher than the R 2 of the M-I-Sta (R 2 = 0.50, Figure 5a), but the R 2 of the M-II-Sta with BLH Sta RE (R 2 = 0.51, Figure 7b) is similar to the R 2 of the M-I-Sta.The results indicate that BLH Sta SBI but not BLH Sta RE can effectively improve the model performance for the prediction of PM 2.5 .Moreover, the regression coefficient associated with BLH Sta SBI is 0.38 (Figure 7a), but the regression coefficient associated with BLH Sta RE is −0.12 (Figure 7b).This is consistent with the negative correlation between BLH Sta SBI and BLH Sta RE shown in Figure 3.We suggest that the positive regression coefficient is reasonable, since a higher BLH Sta SBI means a deeper inversion that strongly inhibits the diffusion of surface pollution and increases the PM 2.5 concentration [53][54][55][56][57].
Sta) and convective BLH (M-II-Con), respectively.We expect to find the optimal BLHs out of the two stable BLHs and six convective BLHs.The samples of the M-II-Sta and the M-II-Con should be matched with the BLH based on the samples of the M-I-Sta and the M-I-Con, respectively.There are 156 samples for the M-II-Sta.For the M-II-Con, however, the sample sizes vary with the different convective BLHs.There are 60 and 149 samples after matching with BLH and BLH , respectively.There are 186 samples after matching with BLH , BLH , BLH and BLH .Figure 7 shows the scatter plots of the observed versus the predicted surface PM2.5 concentrations for the results of the M-II-Sta.The R 2 of the M-II-Sta with BLH (R 2 = 0.59, Figure 7a) is significantly higher than the R 2 of the M-I-Sta (R 2 = 0.50, Figure 5a), but the R 2 of the M-II-Sta with BLH (R 2 = 0.51, Figure 7b) is similar to the R 2 of the M-I-Sta.The results indicate that BLH but not BLH can effectively improve the model performance for the prediction of PM2.5.Moreover, the regression coefficient associated with BLH is 0.38 (Figure 7a), but the regression coefficient associated with Figure 8 shows the results of the M-II-Con with the BLH Con θ v , BLH Con EI , BLH Con θ , BLH Con RH , BLH Con N , and BLH Con RE , respectively.The performance of the M-II-Con with BLH Con EI is optimal, with R 2 = 0.69 (Figure 8b), and the performance of the M-II-Con with BLH Con RE is the worst, with R 2 = 0.66 and RMSE = 31.66µg/m 3 (Figure 8f).It was found that the performance differences were not significant among these models of BLHs.Note that the regression coefficients associated with convective BLHs are all negative, and the largest is −0.28 (Figure 8b).It can be inferred that a higher convective BLH means a higher bottom of inversion.As such, the negative regression coefficient for convective BLH indicates that PM 2.5 concentrations from the surface are diluted in the boundary layer as convective BLH increases.However, the regression coefficient associated with the BLH Sta SBI is positive in Figure 6a, since there is not a lifting of the bottom of the surface-based inversion but rather a stronger inhibition of aerosol diffusion.Thus, the effects of the BLH on the model of AOD-PM 2.5 are contrary for the stable and convective boundary layers.

Optimal Model to Estimate Surface PM2.5 Concentrations
The BLH shows optimal performance in the M-II-Sta, and the BLH shows optimal performance in the M-II-Con model.Thus, these two BLHs are still used as factors in the stepwise

Optimal Model to Estimate Surface PM 2.5 Concentrations
The BLH Sta SBI shows optimal performance in the M-II-Sta, and the BLH Con EI shows optimal performance in the M-II-Con model.Thus, these two BLHs are still used as factors in the stepwise regression model of the M-III-Sta and the M-III-Con, respectively.The sample size of the M-III-Sta is 156, the same as the sample size of the M-II-Sta.For the M-III-Con, the samples size for the M-III-Con is 149, the same as the sample size of BLH Con EI in M-II-Con.Table 3 shows the results of the stepwise regression of the M-III-Sta.The optimal regression is the fourth equation, which includes the factors of surface relative humidity (RH sur ), BLH Sta SBI , and surface temperature (T sur ).The AOD, RH sur , and BLH Sta SBI are all significant at the 99% confidence level (α = 0.01), and T sur is significant at the 95% confidence level (α = 0.05).By introducing these factors, the R 2 increases from 0.50 to 0.68, and the RMSE decreases from 44.30 to 35.44 µg/m 3 .The first-order factor's influence on the two-variate regression (AOD-PM 2.5 ) is RH sur , which improves the R 2 from 0.50 to 0.62, and the second-order factor is BLH Sta SBI , which improves the R 2 from 0.62 to 0.65.Note that there is a more significant improvement if the factor of BLH Sta SBI is introduced first (Figure 6a), such that the R 2 increases from 0.50 to 0.59.This is because of the complex correlation between BLH Sta SBI , RH sur , and PM 2.5 .The influence of BLH Sta SBI on PM 2.5 is partly represented by RH sur in the second or third equation of Table 3, which decreases the improvement of the BLH Sta SBI .It could be inferred that a higher BLH Sta SBI means a stronger surface-based inversion layer that inhibits the diffusion of water vapor and aerosols, subsequently raising the RH sur and PM 2.5 .Thus, the regression coefficients of BLH Sta SBI and RH sur are all positive in the equation for predicting PM 2.5 .In the fourth equation of Table 3, the regression coefficient of T sur is negative (−0.01).It may be speculated that a higher T sur means a weaker surface-based inversion that decreases the surface PM 2.5 .In a word, the influences of these factors on the PM 2.5 are all reasonable based on the correlation between surface-based inversion and PM 2.5 .Table 4 shows the results of the stepwise regression of the M-III-Con.The optimal regression is the third equation, which includes the factors of AOD, BLH Con EI , and surface wind speed (WS sur ).By introducing these factors, the R 2 increases from 0.65 to 0.70 and the RMSE decreases from 31.76 to 29.89 µg/m 3 .The first-order factor's influence on the two-variate regression (AOD-PM 2.5 ) is BLH Con EI , which improves the R 2 from 0.65 to 0.69.The second-order factor is WS sur , which slightly improves the R 2 from 0.69 to 0.70; no other significant factors are introduced into the stepwise regression.It is suggested that the prediction accuracy is difficult to improve for the model of the convective type since the aerosol is uniformly mixed in this situation.Moreover, in the situation of the convective boundary layer, RH sur is also uniformly mixed and its value is usually less than 80% near the ground.There are only 14 samples with the RH sur being greater than 80% out of all 186 samples.The influence of humidity on the correlation of AOD-PM 2.5 is small because of the slight aerosol hygroscopic growth in the condition of RH sur less than 80% [58,59].Unlike RH sur , the factor of WS sur may vary greatly in the convective boundary layer.A higher WS sur means a stronger convective process and a higher BLH Con  EI , but a lower PM 2.5 .As such, the regression coefficient of WS sur is negative (−0.22) in the optimal stepwise equation of Table 4.However, the WS sur is not introduced in the stepwise equation of Table 3, since it is generally small in the stable boundary layer [60][61][62].Thus, the number, order and effect of factors in the stepwise equation depend on the boundary layer types.

Conclusions
To evaluate the effect of the BLH on the regression model of AOD-PM 2.5 , we divided the BLH into the stable and convective boundary types.These two types of boundary layer obviously exhibit seasonal difference.The percentage of the stable boundary layer is highest for winter.In contrast, the percentage of the stable boundary layer is lowest for summer.In addition, these two types of BLH are estimated by several different calculation methods, respectively.The stable BLH is estimated by the method of the top of surface-based inversion layer from the radiosonde profile (BLH Sta SBI ) and by the method of interpolation from NCEP reanalysis (BLH Sta RE ).The convective BLH is estimated by the methods using the profiles of virtual potential temperature (BLH Con θ v ), potential temperature (BLH Con θ ), elevated temperature inversion (BLH Con EI ), relative humidity (BLH Con RH ), refractivity (BLH Con N ), and by another method from NCEP reanalysis (BLH Con RE ).It is found that there are significant differences among the BLHs determined by these methods.The correlation between the two stable BLH of BLH Sta SBI and BLH Sta RE is even negative (−0.29).Moreover, the correlations between the convective BLHs are all lower than 0.70 with the lowest being 0.01.The differences between these BLHs could cause different regression coefficients and prediction accuracies in the regression model of AOD-PM 2.5 .
The regression model of AOD-PM 2.5 is also divided into stable and convective types according the type of the BLH.For the model with AOD as the only factor, it is found that the prediction accuracy of the model for the stable type (M-I-Sta) is significantly lower than that of the model for the convective type (M-I-Con).After the BLH was introduced into the regression models (M-II-Sta and M-II-Con), the improvement of BLH for the M-II-Sta is more significant than that for the M-II-Con.For the M-II-Sta, the optimal BLH is BLH Sta SBI , which improved the R 2 from 0.50 to 0.59.For the M-II-Con, the optimal BLH is BLH Con EI , which improved the R 2 from 0.65 to 0.69.It is found that the regression coefficient of BLH Sta  SBI in the M-II-Sta is positive, but the regression coefficient of BLH Con EI in the M-II-Con is negative.This result suggests that that the effects of BLH are contrary in the types of stable and convective boundary layers.
Based on the M-II-Sta and the M-II-Con models with the optimal BLH, other surface meteorological parameters were introduced into the regression model (M-III-Sta and M-III-Con) using the stepwise regression method.For the M-III-Sta, the optimal regression improved the R 2 from 0.50 to 0.68, which included the factors of surface relative humidity, BLH Sta SBI , and surface temperature.The surface relative humidity is the first-order factor with a positive coefficient.For the M-III-Con, the optimal regression included the factors of BLH Con EI and surface wind speed.It is found that the number and the order of factors are different in the two types of models, and the effects of the meteorological factors on the model also closely related to the types of BLHs.
Most notably, this is the first study to our knowledge to investigate the effectiveness of BLH on the regression model of PM 2.5 -AOD with respect to the stable and convective boundary layers.Our results provide compelling evidence for varying effects and mechanism of these two types of BLHs.However, our results are encouraging and should be validated for more radiosonde stations.In addition, methods to improve the accuracy of the PM 2.5 -AOD model for areas without radiosonde data and only with reanalysis data should also be studied in future work.
NCEP grid data were also interpolated with respect to the AERONET station.Sections 2.1-2.4 describe each data set in more detail.

Figure. 1 .
Figure. 1.The spatial distribution of the PM2.5 station, radiosonde station and Aerosol Robotic Network (AERONET) station used in this study.

Figure 1 .
Figure 1.The spatial distribution of the PM 2.5 station, radiosonde station and Aerosol Robotic Network (AERONET) station used in this study.

3 . 2 .
The Method of Estimating PM 2.5

Figure 2 .
Figure 2. The average vertical temperature profiles of stable BLH and convective BLH for spring (a), summer (b), autumn (c) and winter (d).

Figure 3
Figure 3 shows the averaged BLHs of Beijing during 2014-2015 calculated using the methods described in Section 3.1.The left two bars belong to the stable type of BLH, and the right bars belong to the convective type of BLH.The sample sizes are presented above the bars.The total sample size of sound profiles is 730 during 2014-2015.They are divided into two types based on the existence of a surface-based inversion layer.The sample size of stable BLH is 314, and that of the other is 416.Note that the sample size of BLH is only 104, which is much less than the sample size of the other convective BLHs.Most samples with the BLHare in summer and spring, since the sunrise is early in the morning for these two seasons.The surface temperature and θ quickly increase after sunrise, and the boundary layer will transform from stable to convective.However, the boundary layer is more likely to be neutral in the seasons of autumn and winter, and there is not a θ on the aloft profile equal to the surface θ for these samples with neutral stratification.We can increase the sample sizes of BLH by adding an excess δθ on the surface θ[17].However, it is difficult to determine the δθ using the radiosonde data for its coarse vertical levels.In addition, the sample size of BLH is also less than the sample sizes of other convective BLHs, because there are a few samples without elevated inversion in the total samples of the convective boundary layer[18].

Figure 2 .
Figure 2. The average vertical temperature profiles of stable BLH and convective BLH for spring (a), summer (b), autumn (c) and winter (d).

Figure 3
Figure 3 shows the averaged BLHs of Beijing during 2014-2015 calculated using the methods described in Section 3.1.The left two bars belong to the stable type of BLH, and the right bars belong to the convective type of BLH.The sample sizes are presented above the bars.The total sample size of sound profiles is 730 during 2014-2015.They are divided into two types based on the existence of a surface-based inversion layer.The sample size of stable BLH is 314, and that of the other is 416.Note that the sample size of BLH Con θ v is only 104, which is much less than the sample size of the other convective BLHs.Most samples with the BLH Con θ v are in summer and spring, since the sunrise is early in the morning for these two seasons.The surface temperature and θ v quickly increase after sunrise, and the boundary layer will transform from stable to convective.However, the boundary layer is more likely to be neutral in the seasons of autumn and winter, and there is not a θ v on the aloft profile equal to the surface θ v for these samples with neutral stratification.We can increase the sample sizes of BLH Con θ v by adding an excess δθ v on the surface θ v[17].However, it is difficult to determine the δθ v using the radiosonde data for its coarse vertical levels.In addition, the sample size of BLH Con EI is also less than the sample sizes of other convective BLHs, because there are a few samples without elevated inversion in the total samples of the convective boundary layer[18].

Figure 3 .
Figure 3. Average BLHs of stable boundary layer (a) and convective boundary layer (b) in Beijing during 2014-2015 calculated by different methods.

Figure 3 .
Figure 3. Average BLHs of stable boundary layer (a) and convective boundary layer (b) in Beijing during 2014-2015 calculated by different methods.

Figure 3 .
Figure 3. Average BLHs of stable boundary layer (a) and convective boundary layer (b) in Beijing during 2014-2015 calculated by different methods.

A
simple lognormal model (M-I) was developed to describe the relationship between surface PM2.5 concentrations and AOD.To analyze the performance of the model for the two types of boundary layers, the M-I was divided into two types for the stable boundary layer (M-I-Sta) and the convective boundary layer (M-I-Con).There are 342 samples after matching the AOD and PM2.5 data during 2014-2015.These samples are also divided into two parts for the two types of models.There are 156 and 186 samples for the M-I-Sta and M-I-Con after matching, respectively.

Figure 6 .
Figure 6.Scatter plots of the observed vs. the predicted surface PM2.5 concentrations for the M-I-Sta (a) and the M-I-Con (b).

Figure 7 .
Figure 7. Scatter plots of the observed vs. the predicted surface PM2.5 concentrations from the M-II-Sta with BLH (a) and BLH (b).

Figure 6 .
Figure 6.Scatter plots of the observed vs. the predicted surface PM 2.5 concentrations for the M-I-Sta (a) and the M-I-Con (b).

Figure 7 .
Figure 7. Scatter plots of the observed vs. the predicted surface PM2.5 concentrations from the M-II-Sta with BLH (a) and BLH (b).

Figure 7 .
Figure 7. Scatter plots of the observed vs. the predicted surface PM 2.5 concentrations from the M-II-Sta with BLH Sta SBI (a) and BLH Sta RE (b).

Figure 8 .
Figure 8. Scatter plots of the observed vs. the predicted surface PM 2.5 concentrations from the M-II-Con with BLH Con θ v (a), BLH Con EI (b), BLH Con θ (c), BLH Con RH (d), BLH Con N (e) and BLH Con RE (f).

Table 1 .
Seven methods for estimating boundary layer height (BLH).

Table 2 .
Number and percentage of samples of stable and convective boundary layers for different seasons during 2014-2015.

Table 3 .
Results of the stepwise regression of M-III-Sta (sample size n = 156).
Note: the numbers in parentheses represent the regression coefficients of stepwise regression model; * and ** at the upper right corner of number represent statistically significant at the 0.05 and 0.01 level, respectively.

Table 4 .
Results of the stepwise regression of the M-III-Con (sample size n = 149).: the numbers in parentheses represent the regression coefficients of stepwise regression model; ** at the upper right corner of number represent statistically significant at the 0.05 and 0.01 level, respectively. Note