Comparison of four different types of planetary boundary layer heights during a haze episode in Beijing

Comparison of four different types of planetary boundary layer heights during a haze episode in Beijing

Journal Pre-proofs Comparison of four different types of planetary boundary layer heights during a haze episode in Beijing Yu Shi, Fei Hu, Zhisheng Xi...

3MB Sizes 0 Downloads 72 Views

Journal Pre-proofs Comparison of four different types of planetary boundary layer heights during a haze episode in Beijing Yu Shi, Fei Hu, Zhisheng Xiao, Guangqiang Fan, Zhe Zhang PII: DOI: Reference:

S0048-9697(19)34920-4 https://doi.org/10.1016/j.scitotenv.2019.134928 STOTEN 134928

To appear in:

Science of the Total Environment

Received Date: Revised Date: Accepted Date:

14 August 2019 8 October 2019 10 October 2019

Please cite this article as: Y. Shi, F. Hu, Z. Xiao, G. Fan, Z. Zhang, Comparison of four different types of planetary boundary layer heights during a haze episode in Beijing, Science of the Total Environment (2019), doi: https://doi.org/10.1016/j.scitotenv.2019.134928

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 The Author(s). Published by Elsevier B.V.

1

Comparison of four different types of planetary boundary layer heights during a haze episode in Beijing Yu Shia, Fei Hua,b,*, Zhisheng Xiaoc, Guangqiang Fand, and Zhe Zhanga,b aState Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China bUniversity of Chinese Academy of Sciences, Beijing 100049, China cChinese Academy of Meteorological Sciences, Beijing 100081, China dKey Laboratory of Environmental Optics and Technology, Anhui Institute of Optics and Fine Mechanics, Chinese Academy of Sciences, Hefei 230031, China *Corresponding author. Email address: [email protected] ABSTRACT The planetary boundary layer (PBL) height mainly determines the environmental capacity for the diffusion of atmospheric pollutants, and has always been a hot issue in the study of air pollution. However, there still remains great uncertainty, partly because different PBL heights definitions and the PBL heights are obtained by various measurement instruments. Pollutants are the substances emitted, different from the atmospheric background physical properties such as wind, temperature and turbulence flux that always exist even without pollution. It is very important to distinguish PBL heights obtained from wind, temperature, turbulence quantities and the concentration of pollutants. In this paper, we express the PBL heights determined on the above four parameters as Hu, Hθ, Ht and Hc respectively, and compare them during a heave haze pollution process in Beijing using observation data and simulation results. The comparison results show that: (1) Hθ, namely the inversion layer height, decreased from approximately 1250 m to 450 m from 26 to 30 December, resulting in deteriorating pollution situation. Hc, calculated by lidar and characterizes the maximum depth of vertical diffusion of particulates, also dropped below 500 m, and on the whole, the values of Hc estimated by gradient method and Hθ were in good agreement; (2) Generally, Hc was relatively lower than Hθ and Hu, despite a high bias caused by the existence of the residual layer, multilayer aerosol structure, or lower inversion; (3) Ht estimated from turbulence quantities simulated by WRF model mainly approximated Hu, Hθ and Hc in the daytime during haze pollution, however for the nocturnal boundary layer height in the winter, Ht was seriously underestimated. The averaged PBL heights according to the pollution level showed that Hc, Hθ, Hu and Ht differed greatly on clean days, and the maximum PBL height Hu exceeded 1400 m. On clean days, the inversion intensities observed were lower, so the blocking effect of the inversion layer to pollutant diffusion was not strong enough, Hθ (886m) deviated from Hc (1111m). However, Hc and Ht were very close, approximately 1100 m. The decrease of PBL height led to heavy pollution, Hc, Hθ and Ht were almost

2 700 m. Hu was slightly higher and reduced by about 450 m during heavy pollution. The detailed analyses and comparisons of the PBL height from different variables can help improve the rational application of different methods in the determination of PBL height. Key words: PBL height, Lidar, Radiosonde, WRF, Haze 1. Introduction The planetary boundary layer (PBL) is about 1-2 km away from the ground, adjacent to the Earth’s surface (Stull, 1988), and also is an important place for human activities and the major occurrence place for air pollution (Brook et al., 2004; Sun et al., 2014). The PBL height is a very important parameter in the atmospheric boundary layer and determines the atmospheric environmental capacity for the mixing, diffusion and transport of pollutants emitted into the PBL (Tegen et al., 2000; Dawson et al., 2007; Kleeman, 2008). The PBL height is also of great scientific and practical significance for weather, climate, and air pollution prediction. (Suarez et al., 1983; Holtslag and Nieuwstadt, 1986; Seibert et al., 2000). The most common definition of PBL height is the height up to which the influence of the presence of the lower surface is detectable (Lenschow, 1986). The air flows within the PBL are mainly turbulent (Caughey et al., 1979; Stull, 1988), so many definitions of PBL height are based on the level where the turbulence quantities basically decrease significantly to a proportion of their surface values (Mahrt and Vickers, 2006; Balsley et al., 2006) (e.g., when the Turbulent Kinetic Energy, TKE, decreases to 5 percent of its ground value). This is turbulent boundary layer height. However, direct measurements of turbulence fluctuations are very difficult, rare and expensive in practical routine observations. Other methods for determining PBL heights are based on the vertical distribution of the temperature, wind, water vapor or even the pollutants in the PBL. However, previous studies have seldom mentioned that the concentration of conservative pollutants in the atmosphere is essentially different from physical properties of the atmosphere itself. Despite the strong control strategies exerted by China government, haze events still often occur in Beijing, the capital of China, especially in the winter when the diffusion conditions are unfavorable. Aerosol lidar has been widely used to monitor pollutant concentrations. Lidar can collect the backscattering signals of aerosol particles and then to retrieve extinction coefficient. The pollutant concentration in the PBL is much higher compared to free atmosphere, so that there is often a sudden change in the extinction coefficient signal from the boundary layer top to the free atmosphere (Emeis et al., 2008; Summa et al., 2013). And this phenomenon has been used to obtain PBL heights, and the calculation methods are also improving continually (Balsley et al., 2006; Yang et al., 2017). This is material boundary layer height representing the depth where pollutants concentrated. Many studies have shown that the PBL heights obtained by lidar and radiosonde matched well under unstable and stable conditions (Davis et al., 2000).

3 Moreover, the correlation coefficient was higher under unstable conditions (Eresmaa et al., 2005; Emeis and Schäfer, 2006; Martucci et al., 2007), and the reason for a lower correlation coefficient at night was the effect of aerosol multilayer structure (Münkel et al., 2007). The daily averaged PBL heights from three remote sensing observations wind profiler radar, microwave radiometer and lidar demonstrated that lidar overestimated the PBL height in the morning and at night (Quan et al., 2013). The PBL height estimated from radiosonde detection is generally used as a reference for evaluating remote sensing techniques (Beyrich, 1997; Seidel et al., 2012). The most common methods for determining PBL heights by the radiosonde are based on the potential temperature profiles observed, so that the PBL heights determined by potential temperature are thermal boundary layer heights representing the atmospheric thermal properties. The shortcoming lies in its poor continuity, so that the complete diurnal variation of the PBL height cannot be obtained. The radiosonde launch times are 08:00 and 20:00 (local station time, LST, the same below), the transitional period of the boundary layer for Beijing. The PBL heights calculated by radiosonde were found to be higher in spring and summer than those in fall and winter (Guo et al., 2016; Si et al., 2018). Wind is the basic information vector in the atmosphere, and its variations are highly random. Unlike the rare observations of lidar, wind observations are very widely and routinely recorded; thus, wind observations have played an important role in the prediction of pollution situations and the development of wind energy resources. Generally, the height where boundary layer wind converts to the geostrophic wind in free atmosphere is used to estimate PBL height (Pichugina and Banta, 2010; Beyrich, 1997), and this is dynamical boundary layer height. Previous studies have consistently focused on the features of PBL height; however, it is noticeable that the variables used to determine PBL heights are different, in particular, due to frequent occurrences of haze pollution in China, aerosol lidar has been widely applied to estimate PBL heights. The distinctions between conservative pollutant matters and traditional atmospheric physical parameters to obtain PBL height have not yet been clearly figured out. Based on the above considerations, the PBL heights during heavy haze pollution from 26 to 31 December (Dec) in Beijing were classified and studied with a combination of aerosol lidar, radiosonde and WRF model. The WRF model application is because the turbulence profiles measured are very rare. Although there are towers, the observation height is limited. The boundary layer parameterization schemes in the WRF model can output the PBL heights through turbulent kinetic energy (TKE) or bulk Richardson number (Hu et al., 2010; LeMone et al., 2014; Banks et al., 2016), so the turbulent boundary layer height in this paper is simulated at present. For convenience, in this paper the PBL heights were classified based on four different variables are listed as follows: Hθ is estimated from the potential temperature; Hu is calculated by the wind profiles; Ht is determined by turbulence fluctuations such as turbulence kinetic energy (TKE) or the

4 bulk Richardson number. Hθ, Hu and Ht represent the height of thermal boundary layer, dynamic boundary layer and turbulent boundary layer respectively (Smedman, 1988). Hc is material boundary layer height, proposed in this paper, reflecting the depth where pollutants accumulate. Hc is determined by conservative atmospheric pollutant concentration,and the aerosol concentration detected by aerosol lidar is used in this paper. 2. Data and method 2.1 Observation data description All data measured come from the heavy haze pollution lasting about one week in Beijing from 26 to 31 Dec 2017. Beijing, the capital of China, is a typical urban megacity with a population exceeding 21 million. Beijing is surrounded by Yanshan Mountains and Taihang Mountains, and Bohai Bay is in the east (Fig. 1a). This kind of topography is not conducive to the diffusion of air pollutants (Sun et al. 2014; Wang et al. 2018). The aerosol lidar (AGHJ-I-Lidar, China) used in this paper is located at the Institute of Atmospheric Physics (IAP), Chinese Academy of Sciences (39.96ºN, 116.36ºE, Fig. 1b), between the north third and fourth ring roads of Beijing. Aerosol lidar can provide high-resolution aerosol extinction coefficients at wavelengths of 532 nm and 355 nm at a vertical resolution 7.5 m and a temporal resolution 5-10 min. Surface visibility and other meteorological variables were measured by the observation station of Beijing Meteorological Bureau (the position is shown in Fig. 1b, and the international code is ZBAA) from Wyoming website (http://weather.uwyo.edu/surface/meteorogram/ index.shtml). Twice daily (08:00 and 20:00) routine radiosonde data at ZBAA were obtained from the Wyoming website http://weather.uwyo.edu, and the vertical resolution below 2500 m is 4-5 m, and the resolution of radiosonde data affects the estimation results for PBL heights (Seidel et al., 2012). 1-h averaged surface air pollutants data at the Olympic Sports Central Stadium (OSCS, Fig. 1b) environmental monitoring station, about 2 km to the northeast of IAP, was downloaded from the official website of the Beijing Environmental Protection Agency (http://beijingair.sinaapp.com).

5 Fig. 1. (a) Local topography of Beijing and its surrounding areas. (b) The locations of the observation sites in Beijing; red circle: IAP (aerosol lidar), blue circle: ZBAA observation station, cyan circle: pollution observation station OSCS positioned approximately 2 km northeast of the IAP. Both maps are from Google maps. 2.2 PBL height calculation methods As a remote sensing technique that has recently been widely used in the atmospheric environment, aerosol lidar can continuously monitor the extinction coefficient vertical profiles (Seibert et al. 2000; Shimizu et al., 2016). Three popular algorithms are used in this paper to extract the PBL heights based on the extinction coefficient (designated as EX), namely the gradient method (Lidar_gra) (Flamant et al., 1997), standard deviation method (Lidar_std) (Hooper and Eloranta, 1986) and wavelet covariance method (Lidar_wav) (Davis et al., 2000; Brooks, 2003). The calculation formulas for the gradient method and standard deviation method are as follows:

 EX ) z

Hc(Lidar _ gra)  max(

(1)



n

Hc(Lidar _ std )  max(

i 1

(E(zi )  E(z))2 n

)

(2)

The wavelet covariance transform technique was also used to derive the PBL heights from the extinction coefficient profiles. The mother wavelet function (Haar wavelet) function h(x) is defined as follows:

a  1, b   ≤z≤b  2  z b  a )  1, b<z≤b  h( 2 a  0, elsewhere  

(3)

where a is called the dilation of the function, b is the translation of the Haar wavelet, and z is the altitude. The wavelet covariance transform Wf (a,b) is expressed as: 1

Wf (a, b)  a

zb

 EX (z)h(

za

z b )dz a

(4)

where EX(z) is the extinction coefficient at z of the lidar, zb represents the lower limit and zt is the higher limit of the profile. The PBL height calculated with Lidar_wav means that:

Hc(lidar _ wav)  max[Wf (a, b)]

(5)

The PBL height estimated from the potential temperature (θ) is calculated with the formula ( means a discrete difference):

6

 ) z

H =max(

(6)

There are many ways to determine the PBL heights from wind profiles (Pichugina and Banta, 2010; Beyrich, 1997). Because of the strong randomness of wind, it is difficult to utilize a uniform standard method. The PBL heights Hu can be determined by a significant change in the wind direction or the wind shear. 2.3 Brief introduction of WRF numerical simulation model The haze pollution was also simulated by the WRF model (version 3.9) to obtain the turbulent PBL heights and weather circulations. The simulated area is divided into two domains, the grid resolution of the first domain is 36 km, with 210×200 grid points, the grid resolution of the nest domain is 12 km, and the number of grid points is 130×112. The initial and boundary condition data for starting up the model are based on the European Centre for Medium-Range Weather Forecasts (ECMWF, ERA-5) reanalysis 0.75×0.75. The parameterization schemes used in this paper are listed in Table. 1.

Table 1. Parameterization schemes of the WRF model used in this paper. Parameterization scheme Scheme name Microphysical WSM3 Shor wave radiation Dudhia Long wave radiation RRTM Boundary layer YSU/MYJ Land surface Noah Surface layer MM5

3. Results and discussion 3.1 Variation in the surface pollutant concentration and meteorological conditions Before studying the PBL height, some basic characteristics of this pollution process were analyzed. Fig. 2 shows the temporal evolution of PM2.5, SO2, NO2, O3 concentrations, visibility and other various meteorological elements observed on the ground. The surface

7 PM2.5 concentration measured at the OSCS increased sharply, ranging from 4 g m-3 to 247 g m-3 from 26 to 30 Dec 2017. Additionally, the NO2 concentration monitored on the ground in Beijing was relatively higher. Visibility, serving as a representative index of air quality (Singh et al. 2016), continued to deteriorate during this period, dropping from 10 km during clean days to about 1 km in polluted episode. The diurnal variations in temperature and RH were not obvious when the pollution situation was serious from 29 to 30 Dec and the ground wind directions belonged to the southeast winds, generally bringing abundant water vapor from Bohai Bay to the Beijing area. After 12:00 30 Dec, relatively high wind speeds from the northwest were observed on the ground, and the PM2.5 concentration dropped rapidly, the visibility increased and the air quality improved significantly.

Fig. 2. Time series of PM2.5 (a), SO2, NO2, and O3 concentrations at the surface (b), visibility (c), pressure (d), RH and temperature (e), wind direction and wind speed (f) during 26 to 31 Dec 2017; the units of all variables are shown in the figure.

3.2 Material boundary layer height Hc calculated by aerosol lidar Fig. 3 shows the temporal and spatial evolutions of the extinction coefficient observed by lidar in the pollution process. The white triangles in the figure represent the PBL heights H obtained by the potential temperature profiles observed by radiosonde. The height where the potential temperature gradient reaches its maximum is defined as the PBL height Hθ, namely the inversion layer height. The red line (Lidar_gra), grey line (Lidar_std) and purple line (Lidar_wav) are the extracted PBL heights Hc by the gradient method, standard deviation method and wavelet method using extinction coefficient, respectively. The

8 calculation formulas of these three methods have been explained in Section 2.2. Additionally, the PBL heights calculated based on lidar are smoothed.

Fig. 3. Temporal and spatial variations of extinction coefficient (shaded, unit: km-1) from 26 to 31 Dec 2017, and the PBL heights (m) determined by different methods. The red line, grey line, and purple line represent Hc determined by lidar using gradient method, standard deviation method, and wavelet method, respectively. The white triangles are H determined by potential temperature gradient. As shown in Fig. 3 the PBL heights determined by the gradient method and the wavelet method were very similar, but the values calculated by standard deviation method departed slightly. Moreover, Hc determined by the three methods was basically below the inversion layer cap from 26 to 30 Dec, revealing that most aerosols were mainly concentrated below Hθ. If the H observed by the radiosonde was taken as a reference, Hc obtained by the gradient method was more accurate. Hθ was generally less than 1500 m, mostly concentrated in the vicinity of 1000 m. The inversion layer decreased to about 450 m at midnight on 30 Dec when the PM2.5 reaching peak. The minimum Hc was basically consistent with the PBL heights reported in previous studies during pollution days in the winter over the Beijing area (Zhang et al. 2015). In the early morning on 30 Dec 2017, the downward trend of the inversion layer has seriously hindered the diffusion of particles, therefore the extinction coefficient increased significantly, reaching beyond 3 km-1. It is noteworthy that there were some moments when Hc was higher than Hθ; this phenomenon was discussed in detail in Section 3.3. 3.3 Thermal boundary layer height H estimated from the temperature profiles

9 The pollutants distribution is closely related to the inversion layer, so the potential temperature profiles observed by radiosonde were analyzed during this pollution process. The formula (Stull, 1988) to calculate the potential temperature is as follows:

1000 0.286 ) p

 =T (

p is air pressure and T is temperature. The potential temperature , RH, and corresponding extinction coefficient profiles at some moments were selected to analyze the inversion layer distribution and its influence on Hc. As shown in Fig. 4, the heights at which the potential temperature gradient achieved its maximums basically decreased from 20:00 26 Dec to 08:00 30 Dec. At 20:00 26 Dec, Hθ approximated 1100 m with an inversion intensity of about 0.05 C m-1, and the RH from the ground to the upper air was around 40%. At 08:00 27 Dec, H appeared at 750 m, and the residual layer characteristics were maintained from 750 m to about 1400 m with an RH nearing 80%. Combined with the wind profile, the whole layer was southwest wind at this time, revealing that the southwest wind could indeed bring water vapor to Beijing area. Moreover, there were two peaks in the extinction coefficient profile, demonstrating that a multilayer aerosol structure, and the extinction coefficient was larger at approximately 500 m. The extinction coefficient in the residual layer was still high, indicating that there were still some pollutants remaining in the residual layer. The residual layer or multilayer aerosol structure can affect the determination of Hc, causing Hc a little higher than Hθ. By 20:00 27 Dec, H was approximately 1236 m and H agreed well with Hc, but the potential temperature at 1500 m increased with the inversion intensity reaching nearly 0.14 C m-1 due to the effect of warm advection. The RH at 500~1000 m increased to 90 % where there was still southeast wind, and the extinction coefficient within this layer was obviously higher. The maximum extinction coefficient below the inversion layer exceeded 3 km-1 and had only one peak value at 20:00 28 Dec. The inversion intensity was maintained at 0.08 C m-1, and Hc was a little lower than Hθ. Until 08:00 30 Dec, the inversion layer eventually decreased to approximately 450 m, and meanwhile, the potential temperature and the RH above 500 m were greatly reduced. The maximum extinction coefficient near the surface was close to 5 km-1. The potential temperature profile at 20:00 30 Dec showed almost no obvious change as illustrated in Fig. 4d, and the whole layer was nearly neutral; in addition, at 08:00 31 Dec, except for the surface inversion layer, the PBL height Hθ at this time was approximately 625 m, but the inversion intensity was obviously lower, reducing to 0.04 C m-1. The extinction coefficient decreased dramatically after 20:00 30 Dec at night on 31 Dec, the inversion layer began to form from the ground and there were still some pollutants accumulating under the inversion layer; however, above the inversion layer, the pollutants with a relatively high mixing degree in the daytime had not completely settled under the inversion layer forming form the ground. When the inversion intensity was not low, the

10 blocking effect of the inversion layer on the pollutants was not obvious causing inconsistency between Hc and Hθ.

Fig. 4. Vertical profiles of the potential temperature (, unit: K) and RH (unit: %) observed by the radiosonde at the ZBAA station and extinction coefficients (unit: km-1) measured by lidar. 3.4 Effects of warm advection on the PBL heights In order to further understand the continuous decrease of the PBL heights extracted by lidar and the radiosonde, we also simulated this process using WRF mesoscale numerical model, showing the distribution of potential temperature, water vapor and wind vectors at 850 hPa (approximately 1500 m from the ground). At 08:00 on 26 Dec, the northwest wind prevailed on 850 hPa in the Beijing area (which coincided with the wind direction observed by the radiosonde, as shown in Fig. 6c). There was cold advection, meaning that the wind moved from a cold area to a warm area. At this time, H was approximately 1200 m, the PM2.5 concentration was less than 10 g m-3, and the visibility was good. By 20:00 on 27 Dec, the weak warm advection appeared at 850 hPa in North China, and the intersection angle between the wind direction and the temperature gradient

11 became even smaller at 08:00 on 28 Dec, illustrating stronger warm advection. At 20:00 on 28 Dec, the northwest wind was observed in Beijing, and the warm advection disappeared. The extinction coefficient of this period was basically below 2 km-1. Until 08:00 on 29 Dec, because of the persistent warm advection, the inversion layer height H decreased from more than 1000 m to about 750 m. The continuous decline of the inversion layer seriously compressed the environmental capacity for the vertical diffusion of pollutants, so Hc decreased to nearly 500 m. After 12:00 on 29 Dec, the concentration of PM2.5 increased significantly, and the corresponding extinction coefficient increased rapidly to more than 3 km-1. By 08:00 on 30 Dec, the inversion layer was the lowest at 459 m. The pollution situation after 12:00 on 30 Dec improved significantly. At 20:00 on 30 Dec, the inversion height increased to 1500 m, and Hc also increased to approximately 1200 m. There was northwest wind at 850 hPa, and the warm advection disappeared. (a)

(b)

(c)

(d)

12 (e)

(f)

Fig. 5. Potential temperature (shaded, unit: K) and wind vectors (black arrows, unit: m simulated by WRF at 850hPa: (a) 08:00 on 26 Dec, (b) 20:00 on 27 Dec, (c) 08:00 on 28 Dec, (d) 20:00 on 28 Dec, (e) 08:00 on 29 Dec, and (f) 20:00 on 30 Dec.

s-1)

3.5 Dynamical boundary layer Hu determined from the wind profiles 3.5.1 PBL heights Hu estimated from the wind direction As shown in Fig. 6a, there was an obvious change in the wind direction at 20:00 on 26 Dec; specifically, the southwest wind maintained from the ground to approximately 1200 m, and above 1200 m, the wind direction began to rotate counterclockwise and turned into a northwest wind, implying that it had become a geostrophic wind at the top of the boundary layer. At this time, both Hθ and Hc were near 1160m and were in good agreement with the height where the wind direction began to convert. That is, although Hθ, Hc and Hu are based on different elements, sometimes they are highly consistent, and pollutants usually accumulate below the height at which the wind direction begins to change distinctly. 3.5.2 PBL heights Hu estimated from wind shear Wind shear can generally be used to measure the degree of wind change. When a distinct low-level wind maximum exists, the wind maximum generally corresponds to the place where the wind shear is zero. In this case, for example, 20:00 on 29 Dec (see Fig. 6b), the height of the maximum wind speed or the zero wind shear can determine the PBL height Hu. The wind speed at approximately 1100 m increased beyond 15m s-1, manifesting a typical low-level jet structure. The low-level jet zone usually corresponds to the PBL height. The value of Hc extracted by lidar coincided with Hθ very well, and both values were less than the height of the low-level jet, illustrating that aerosols were also still under the inversion layer or the low-level jet zone. The inversion intensity increased to 0.07 C m-1; wind speeds often affect the existence and maintenance of the inversion layer (Andreas et al., 2000). In addition, the inversion layer prevents momentum exchange between high and low altitudes, and it is easy to trigger a low-level jet above the inversion layer.

13 -0.04 2500

-0.02

0

0.02

0.04

-0.04 2500

2500 WD

U

-0.02

U/ z

0.02

0.04

2500 WD

U/ z

H

H

Hc

2000

0

U

Hc

2000

2000

1500

1500

1500

1500

1000

1000

1000

1000

500

500

500

500

0

0

5

10

0

15

(a)

S

W

N

-0.04 2500

-0.02

E

0

0

S

0.02

0.04

0

2000

5

10

15

0

S

W

N

E

S

(b) 2500 WD

U U/ z H Hc

2000

2000

1500

1500

1000

1000

500

500

0

0

5

10

0

15

S

W

N

E

S

(c)

Fig. 6. Wind profile (U, unit: m s-1), wind shear ( U/z, unit: s-1) and wind direction (WD, unit: degree) observed by the radiosonde at (a) 20:00 on 26 Dec 2017, (b) 20:00 on 29 Dec 2017, (c) 08:00 on 26 Dec 2017. The blue dotted lines (H) represent the PBL heights calculated by the maximum potential temperature gradients, and the green dotted lines (Hc) represent the PBL heights extracted by lidar. Because the wind speed at the ground is zero and the ground can exert friction on the air flow, the wind in the boundary layer increases gradually with height (Stull, 1988). In particular, the convective boundary layer, due to the mixing effect within the boundary layer, often exhibits a sudden change in the wind speed profile at the top of the boundary layer (Huang et al., 2017). The wind profile distribution was characterized by the transformation zone of the wind speed at 08:00 on 26 Dec (Fig. 6c), although this moment did belong to convective boundary layer. The PM2.5 concentration was not very high, H corresponded to the height of the maximum wind shear, and Hc deviated slightly more. Three typical cases to determine Hu have been presented, and the determination of Hu should be based on the specific wind profile features.

14 3.6 Turbulent boundary layer height Ht simulated by WRF model Above all, we have discussed the PBL heights obtained by different methods based on the data observed, and then we analyzed the simulation results by the WRF model to study the characteristics of the turbulent boundary layer height Ht. WRF model can output the PBL heights calculated from the turbulence quantities, the TKE or the bulk Richardson number, and YSU scheme calculating the PBL height based on the bulk Richardson number was used in this paper. In addition, the TKE vertical profiles outputted by MYJ scheme were also given. As shown in Fig. 7e, the PBL height Ht outputted by the YSU scheme was seriously low at 08:00 on 30 Dec. At this time, TKE simulated maintained background value. The wind speed simulated by WRF basically matched with the radiosonde observation (see Fig. 7a) but was a little overestimated below approximately 200 m. The higher wind speed simulated in the low layer may be due to the failure to couple the urban canopy effects. Equally, the wind direction deviated from the measurement in the lower level (see Fig. 7b). At this time, Hu determined by WRF and the radiosonde using the zero wind shear method were basically close to 700 m (see Fig. 7c). The simulated inversion intensity was slightly lower than that of the observation (see Fig. 7d), and generally speaking, the simulated temperature was colder. The inversion layer appeared at approximately 500 m; at this time, Hc (see Fig. 7f) and H were very close,2017/12/30/08:00 and both were below Hu. 2500 (a)

WD(WRF) WD(Radio)

(b)

U(WRF) U(Radio)

2000

(c)

U/ z(WRF) U/ z(Radio)

1500 1000 500 0

0

5

15 S

10

W

N

2500 (WRF) (Radio)

(d)

2000

S -0.05

E

TKE(WRF) Ht(WRF)

(e)

0

(f)

0.05

0.1

EX(Lidar) Hc(Lidar)

1500 1000 500 0 270

280

290

300 0

0.1

0.2

0.3

0.4 0

2

4

6

15 Fig. 7. Vertical profiles of the wind speed (U, unit: m s-1), wind direction (unit: degree), wind shear (U/z, unit: s-1), and potential temperature (, unit: K) of the WRF simulation and the radiosonde observation. The vertical profile of the TKE (unit: m2 s-2) simulated by the WRF model and the vertical profile of the extinction coefficient (EX, unit: km-1) observed by lidar. Ht: the PBL height outputted from the WRF model. Hu, Hc, and H represent the PBL heights determined from the wind, extinction coefficient and potential temperature, respectively. The times of these profiles are at 08:00 on 30 Dec. Hu determined by the maximum wind shear at 02:00 on 29 Dec was 668 m (see Fig. 8c) according to the simulation results, corresponded to the inversion layer height (see Fig. 8d). Ht simulated at night was obviously underestimated by the model, and the TKE maintained at the background value (see Fig. 8e). Hc calculated by gradient method was close to 1 km, and the wavelet transform method (615 m) seemed more reasonable at this moment (see Fig. 8f). 2017/12/29/02:00

2500

2500

2500 U(WRF)

WD(WRF)

(b)

(a)

2000

2000

2000

1500

1500

1500

1000

1000

1000

500

500

500

0

0

5

15

10

2500 2000

0

W

S

N

E

S

2500

0 -0.05

(e)

2000

1500

1000

1000

1000

500

500

500

290

300

0

0

0.2

EX(Lidar) Hc(Lidar)

2000

1500

280

0.05

0

(f)

TKE(WRF) Ht(WRF)

1500

0 270

(c)

2500

(d)

(WRF) H (WRF)

U/ z(WRF) Hu (WRF)

0.4

0

0

2

4

6

Fig. 8. Vertical profiles of the wind speed (U, unit: m s-1), wind direction (unit: degree), wind shear (U/z, unit: s-1), potential temperature (, unit: K), and TKE (unit: m2 s-2) simulated by WRF and extinction coefficient vertical profile (EX, unit: km-1) observed by lidar at 02:00 on 29 Dec. Ht: PBL height outputted from the WRF model using YSU scheme, Hu, Hc and H represent the PBL height determined from the wind, extinction coefficient and potential temperature, respectively.

16

At 14:00 on 29 Dec (see Fig. 9), due to relatively strong surface radiation, Ht outputted from the WRF model was about 635 m; moreover, Ht given by the YSU scheme was basically the same as the minimum height where the TKE decreased to its background value of 0.1 m2 s-2 (LeMone et al., 2014) (see Fig. 9e). The simulation showed that from the ground to approximately 600 m, there was basically a small wind layer, and the value of Hu determined by the maximum wind shear was nearly 668 m (see Fig. 9c). Hu showed no obvious diurnal variation compared with the value at 02:00 on 29 Dec. Additionally, Hu was also basically consistent with the height where the wind direction obviously changed. Because of the strong mixing of the air below the boundary layer during the day, Hc was very consistent with Ht and a little higher than H, revealing that the pollutants were affected not only by the inversion layer but also by the turbulent diffusion motions. Although these elements were different, Hu, H, Ht and Hc were very similar at this 2017/12/29/14:00 moment. 2500

2500

2500 U(WRF)

WD(WRF)

(b)

(a)

(c)

2000

2000

2000

1500

1500

1500

1000

1000

1000

500

500

500

0

0

5

10

15

2500 2000

0

S

W

N

E

S

2500

0 -0.05

(e)

(f)

TKE(WRF) Ht(WRF)

2000 1500

1500

1000

1000

1000

500

500

500

280

290

300

0

0.05

0

0.2

EX(Lidar) Hc(Lidar)

2000

1500

0 270

0

2500

(d)

(WRF) H (WRF)

U/ z(WRF) Hu (WRF)

0.4

0

0

2

4

6

Fig. 9. Similar to Fig. 8 but for 14:00 on 29 Dec. According to the pollution level, the averaged PBL heights Hc, H, Hu and Ht are listed in Table. 2. Because the nocturnal boundary layer heights are seriously underestimated by WRF model, this paper has averaged the four PBL heights from 11:00 to 16:00 in the daytime when the mixing process strong. Hc was calculated based on aerosol liar, H and Hu were estimated according radiosonde data. Hc, H, Hu and Ht differed greatly on clean days without pollution. Hu was the highest, exceeding 1400 m, but the

17 inversion layer height H was only 886 m. The inversion intensities measured on clean days were lower, and the pollutants emitted from the surface were dispersed by the turbulent vertical mixing process, so that Hc and Ht were more similar, approximately 1100 m. The averaged values of Hc, H, Hu and Ht decreased obviously under the heavy polluted condition; Hc, H and Ht were very similar, approximately 700 m on heavily polluted days. Hu also decreased by approximately 450 m, reducing to 900 m. Table 2. Averaged PBL heights of different levels of pollution. good (0-7 5g m-3), slightly polluted (75-115g m-3), moderately polluted (115-150 g m-3), heavily polluted (150-250 g m-3), the values in the square brackets represent the number of data points averaged.

Quality level Good Slightly polluted Moderately polluted Heavily polluted

H

Hc

Hu

Ht

(m)

(m)

(m)

(m)

886[5] 928[3] 714[3]

1111[74] 1070[22] 925[10] 799[35]

1362[5] 1464[3] 914[3]

1057[20] 682[8] 581[3] 784[4]

4. Conclusions There are a variety of comparison work on the PBL height, but it is worth noting that the elements used to obtain PBL heights are different. Pollutants in the atmosphere are substances, different from the physical properties of the atmosphere such as wind, temperature and turbulence that always exist even without pollution. This paper classified four types PBL heights based on four parameters turbulent fluctuations, temperature, wind and aerosol concentration during a haze pollution process in Beijing from 26 to 31 Dec 2017. Four types PBL heights are expressed by Ht, H, Hu and Hc, representing the height of turbulent boundary layer, thermal boundary layer, dynamical boundary layer and material boundary layer respectively. H , also representing the height of the inversion layer, decreased from about 1250m to 450 m, causing an increase in PM2.5 concentration and a rapid deterioration in the visibility. The peak inversion intensity reached nearly 0.14 ºC m-1 due to the persistent warm advection on 850 hPa. Hc decreased from about 1000 m in the clean stage to 500 m in the heavy polluted episode. Hc was basically below H and was in good agreement with H, but the residual layer, multilayer aerosol structure or the low inversion intensity can cause higher Hc than H. Hu estimated from the wind profiles observed by the radiosonde were also analyzed. The methods for determining Hu were based on the wind direction or wind shear. Hu

18 determined from the wind direction approximated H and Hc, meaning that pollutants also accumulated below the height at which the wind direction began to obviously change or began to turn into the prevailing winds. When there were obvious wind maxima, such as a low-level maximum wind or low-level jet, H and Hc were generally below Hu. Strong wind can affect the stable existence of the inversion layer and the height where pollutants accumulate. The turbulent boundary layer height Ht outputted by the WRF model were very close to Hu and Hc in the daytime at 14:00 on 29 Dec, basically near 600 m. Hc and Ht were more consistent, and both were a litter higher than Hc. However, Ht simulated by WRF at night in the winter was seriously underestimated. At this time, the PBL height simulated by the model can be estimated from the vertical profiles of the wind and temperature. The averaged Hc, H, Hu and Ht calculated according to the pollution level showed that Hc, H, Hu and Ht differed greatly during clean days and Hu was highest reaching approximately 1400 m, but the inversion layer height H was only 886 m. Hc and Ht were very similar on clean days, approximately 1100 m, revealing that the diffusion of the pollutants in the boundary layer was more affected by turbulent motions in clean days with low inversion intensity. When the PBL height decreased, the pollution level increased significantly; Hc, H and Ht approximated 700 m on heavily polluted days. Hu was slightly higher, but also decreased by about 450 m compared with clean days. This study helps us to figure out the PBL height between material and physical elements. The physical meaning of four types PBL heights classified in this paper is clear. Detailed analyses and discussions of the PBL height can improve the rational applications of different method in the determination of PBL height. Future work is still necessary to characterize the PBL heights from four different variables using more data and further to obtain more statistical results. Acknowledgments. All the authors have declared that no competing interests exist. This work was supported by the National Key Research and Development Program of China (2017YFC0209605) the National Natural Science Foundation of China (Grant No. 41975018). REFERENCES Andreas, E. L., Claffy, K. J., Makshtas, A. P., 2000. Low-level atmospheric jets and inversions over the western weddell sea. Boundary-Layer Meteorology, 97(3),459-86. doi:10.1023/a:1002793831076. Balsley, B. B., Frehlich, R. G., Jensen, M. L., Meillier, Y., 2006. High-resolution in situ profiling through the stable boundary layer: Examination of the sbl top in terms of minimum shear, maximum stratification, and turbulence decrease. Journal of the Atmospheric Sciences, 63(4):1291-307. doi: 10.1175/JAS3671.1. Banks, R. F., Tiana, A. J., Baldasano, J. M., Rocadenbosch, F., Papayannis, A., Solomos, S., Tzanis, C. G., 2016. Sensitivity of boundary-layer variables to pbl schemes in the

19 wrf model based on surface meteorological observations, lidar, and radiosondes during the hygra-cd campaign. Atmospheric Research, 176:185-201. doi:10.1016/j.atmosres.2016.02.024. Beyrich, F., 1997. Mixing height estimation from sodar data: A critical discussion. Atmospheric Environment, 31(23):3941-53. doi: 10.1016/S1352-2310(97)00231-8. Brook, R. D., Franklin, B., Cascio, W., Hong, Y., Howard, G., Lipsett, M., Luepker, R., Mittleman, M., Samet, J., Smith, S. C., Tager, I., 2004. Air pollution and cardiovascular disease: A statement for healthcare professionals from the expert panel on population and prevention science of the American heart association. Circulation. 109(21):2655-71. doi:10.1161/01.CIR.0000128587.30041.C8. Brooks, I. M., 2003. Finding boundary layer top: Application of a wavelet covariance transform to lidar backscatter profiles. Journal of Atmospheric and Oceanic Technology, 20(8):1092-105. doi:10.1175/1520-0426(2003)020<1092: FBLTAO>2.0.CO;2. Caughey, S. J., Wyngaard, J. C., Kaimal, J. C., 1979. Turbulence in the evolving stable boundary layer. Journal of the Atmospheric Sciences, doi:10.1175/1520-0469(1979) 0362.0.CO;2. Davis, K. J., Gamage, N., Hagelberg, C. R., Kiemle, C., Lenschow, D. H., Sullivan, P. P., 2000. An objective method for deriving atmospheric structure from airborne lidar observations. Journal of Atmospheric and Oceanic Technology,17(11):1455-68. doi:10.1175/1520-0426(2000)017<1455:AOMFDA>2.0.CO;2. Dawson, J. P., Adams, P. J., Pandis, S. N. 2007. Sensitivity of pm 2.5 to climate in the eastern us: a modeling case study. Atmospheric Chemistry and Physics, 7(16):4295309. doi:10.5194/acp-7-4295-2007. Emeis, S., Schäfer, K., 2006. Remote sensing methods to investigate boundary-layer structures relevant to air pollution in cities. Boundary-Layer Meteorology,121(2): 377-85. doi:10.1007/s10546-006-9068-2. Emeis, S., Schäfer, K., Münkel, C., 2008. Surface-based remote sensing of the mixing -layer height: a review. Meteorologische Zeitschrift, 17(5):621-30. doi: 10. 1127/0941-2948/2008/0312. Eresmaa, N., Karppinen, A., Joffre, S. M., Räsänen, J., Talvitie, H., 2005. Mixing height deter-mination by ceilometer. Atmospheric Chemistry and Physics, 6(6):1485-93. doi:10.5194/acp-6-1485-2006. Flamant, C., Pelon, J., Flamant, P. H., Durand, P., 1997. Lidar determination of the entrainment zone thickness at the top of the unstable marine atmospheric boundary layer. Boundary-Layer Meteorology, 83(2):247-84. doi: 10.1023/A: 1000258318944. Guo, J. P., Miao, Y. C., Zhang, Y., Liu, H., Li, Z., Zhang, W., He, J., Lou, M., Yan, Y.,

20 Bian, L., Zhai, P., 2016. The climatology of planetary boundary layer height in china derived from radiosonde and reanalysis data. Atmospheric Chemistry and Physics,16(20):13309-19. doi:10.5194/acp-16-13309-2016. Holtslag, A. A. M., and F. T. M. Nieuwstadt, 1986: Scaling the atmospheric boundary layer. Boundary-Layer Meteorology,36:201-9. doi: 10.1007/BF00117468. Hooper, W. P., and E. W. Eloranta, 1986: Lidar measurements of wind in the planetary boundary layer: The method, accuracy and results from joint measurements with radiosonde and kytoon. Journal of Applied Meteorology, 25(25):990-1001. doi:10.1175/1520-0450(1986)025<0990:LMOWIT>2.0.CO;2. Hu, X. M., Nielsen-Gammon, J. W., Zhang, F., 2010. Evaluation of three planetary boundary layer schemes in the wrf model. Journal of Applied Meteorology and Climatology,49(9):1831-44. doi:10.1175/2010JAMC2432.1. Huang, M., Gao, Z., Miao, S., Chen, F., Lemone, M. A., Li, J., Hu, F., Wang, L., 2017. Estimate of boundary-layer depth over beijing, china, using doppler lidar data during surf-2015. Boundary-Layer Meteorology,162(3):503-22. doi: 10.1007/s10546-0160205-2. Kleeman, M. J., 2008. A preliminary assessment of the sensitivity of air quality in california to global change. Climatic Change, 87(1):273-92. doi: 10.1007/s10584-0079351-3. LeMone, M. A., Tewari, M., Chen, F., Dudhia, J., 2014. Objectively determined fairweather nbl features in arw-wrf and their comparison to cases-97 observations. Monthly Weather Review,142(8):2709-32. doi: 10.1175/MWR-D-13-00358.1. Lenschow, H., 1986. Probing the atmospheric boundary layer. American Meteorological Society. Mahrt, L., Vickers, D., 2006. Extremely weak mixing in stable conditions. BoundaryLayer Meteorology,119(1):19-39. doi: 10.1007/s10546-005-9017-5. Martucci, G., Matthey, R., Mitev, V., Richner, H., 2007. Comparison between backscatter lidar and radiosonde measurements of the diurnal and nocturnal stratification in the lower troposphere. Journal of Atmospheric and Oceanic Technology, 24(7):1231-44. doi:10.1175/JTECH2036.1. Münkel, C., Eresmaa, N., Räsänen, J., Karppinen, A., 2007. Retrieval of mixing height and dust concentration with lidar ceilometer. Boundary-Layer Meteorology, 124(1): 117-28. doi:10.1007/s10546-006-9103-3. Pichugina, Y. L., Banta, R. M., 2010. Stable boundary layer depth from high-resolution measurements of the mean wind profile. Journal of Applied Meteorology and Climatology, 2010, 49(1):20-35.,49(1):20-35. doi:10.1175/2009JAMC2168.1. Quan, J., Gao, Y., Zhang, Q., Tie, X., Cao, J., Han, S., Meng, J., Chen, P., Zhao, D., 2013. Evolution of planetary boundary layer under different weather conditions, and its impact on aerosol concentrations. Particuology, 11(1):34-40. doi: 10.1016/j.partic. 2012.04.005.

21 Seibert, P., Beyrich, F., Gryning, S. E., Joffre, S., Rasmussen, A., Tercier, P., 2000. Review and intercomparison of operational methods for the determination of the mixing height. Atmospheric Environment, 34(7): 1001-27.doi: 10. 1016/S1352-2310(99)00349-0. Seidel, D. J., Zhang, Y., Beljaars, A., Golaz, J., Jacobson, A. R., Medeiros, B., 2012. Climatology of the planetary boundary layer over the continental united states and Europe. Journal of Geophysical Research, 117. doi:10.1029/2012JD018143. Shimizu, A., Nishizawa, T., Jin, Y., Kim, S. W., Wang, Z., Batdorj, D., Sugimoto, N., 2016. Evolution of a lidar network for tropospheric aerosol detection in east asia. Optical Engineering, 56(3):31219-. doi:10.1117/1.OE.56.3.031219. Si, Y., Li, S., Chen, L., Yu, C., Wang, Z., Wang, Y., Wang, H., 2018. Validation and spatiotemporal distribution of geos-5based planetary boundary layer height and relative humidity in china. Advances in Atmospheric Sciences, 35(4): 479-92. doi:10.1007/s00376-017-6275-3. Singh, A., Bloss, W. J., Pope, F. D., 2016. 60 years of uk visibility measurements: impact of meteorology and atmospheric pollutants on visibility. Atmospheric Chemistry and Physics, 17(3): 2085-101. doi:10.5194/acp-17-2085-2017. Smedman, A. S., 1988. Observations of a multi-level turbulence structure in a very stable atmospheric boundary layer. Boundary-Layer Meteorology, 44(3): 23-53. Stull, R. B., 1988. An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, Dordrecht. Suarez, M. J., Arakawa, A., Randall, D. A., 1983. The parameterization of the planetary boundary layer in the ucla general circulation model: Formulation and results. Monthly Weather Review, 111(11):2224-43. doi: 10.1175/ 1520-0493(1983)111<2224:TPOTPB>2.0.CO;2. Summa, D., Girolamo, P. D., Stelitano, D., Cacciani, M., 2013. Characterization of the planetary boundary layer height and structure by raman lidar: comparison of different approaches. Atmospheric Measurement Techniques, 6(12): 3515-25 doi:10.5194/amt-6-3515-2013. Sun, Y., Jiang, Q., Wang, Z., Fu, P., Li, J., Yang, T., Yin, Y., 2014. Investigation of the sources and evolution processes of severe haze pollution in beijing in January 2013. Journal of Geophysical Research,119(7): 4380-98.doi:10.1002/ 2014JD021641. Tegen, I., Koch, D., Lacis, A. A., Sato, M., 2000. Trends in tropospheric aerosol loads and corresponding impact on direct radiative forcing between 1950 and 1990: A model study. Journal of Geophysical Research, 105(22): 26971-89. doi: 10. 1029/2000JD900280. Wang, Q., Sun, Y., Xu, W., Du, W., Zhou, L., Tang, G., Chen, C., Cheng, X., Zhao, X., Ji, D., 2018. Vertically resolved characteristics of air pollution during two severe winter

22 haze episodes in urban beijing china. Atmospheric Chemistry and Physics, 18(4):128.doi:10.5194/acp-18-2495-2018. Yang, T., Wang, Z., Zhang, W., Gbaguidi, A., Sugimoto, N., Wang, X., Matsui, I., Sun, Y., 2017. Technical note: Boundary layer height determination from lidar for improving air pollution episode modeling: development of new algorithm and evaluation. Atmospheric Chemistry and Physics, 17(10): 6215-25. doi:10.5194/acp-176215-2017. Zhang, Q., Quan, J., Tie, X., Li, X., Liu, Q., Gao, Y., Zhao, D., 2015. Effects of meteorology and secondary particle formation on visibility during heavy haze events in beijing, china. Science of The Total Environment, 502:578-84. doi: 10.1016/j.scitotenv.2014.09.079.

23

Dear editors: All the authors have declared that no competing interests exist. Yours sincerely, Yu Shi, Fei Hu, Zhisheng Xiao, Guangqiang Fan, Zhe Zhang Fund Rf: Name of Funder:The National Key Research and Development Program of China,Grant Agreement No. / Award No.:2017YFC0209605 Name of Funder:National Natural Science Foundation of China Grant Agreement No. / Award No.:41975018

24 Highlights

(1) The concentration of pollutants used to estimate PBL height are classified from the atmospheric background physical properties, and Hc was proposed; (2) Hc was relatively lower than Hθ and Hu, despite a high bias caused by the existence of the residual layer, multilayer aerosol structure, or lower inversion; (3) Hc, Hθ, Hu and Ht differed greatly on clean days. During the haze pollution process, Hc, Hθ and Ht were almost 700 m.

25 Graphical abstract