Interaction between aerosol and thermodynamic stability within the planetary boundary layer during wintertime over the North China Plain: aircraft observation and WRF-Chem simulation

Aerosol–planetary boundary layer (PBL) interaction has been proposed as a key mechanism for stabilizing the atmosphere and exacerbating surface air pollution. Although the understanding of this process has progressed enormously, its magnitude and impact remain uncertain and vary widely concerning aerosol types, vertical distributions, synoptic conditions, etc. In this study, our primary interest is to distinguish the aerosol–PBL interaction of absorbing and scattering aerosols under contrasting synoptic patterns and aerosol vertical distributions. Detailed in situ aircraft (King Air 350) measurements and online coupled model Weather Research and Forecasting with Chemistry (WRF-Chem) simulations are explored over the North China Plain (NCP). Furthermore, a long-term PBL stability trend from 1980 to 2020 over the NCP is also investigated. The aircraft measurements and surface observations show that the surface air pollution over the city of Baoding on 3 January is heavier than that on 4 January 2020. In addition, the aerosols are restricted to the low layer on 3 January, whereas the aerosols mix more homogeneously in an upwards direction on 4 January. Thus, we focus on the 2 d with distinct synoptic circumstances, PBL stability, and aerosol vertical distributions over the NCP. According to the WRF-Chem modeling, the synoptic pattern over Baoding differs between the 2 d. The prevailing wind direction is opposite with a southwesterly wind on 3 January and a northeasterly wind on 4 January. The results indicate that the synoptic condition may affect the PBL thermal structure, thus affecting the aerosol vertical distribution. Additionally, the sensitive numerical experiments reveal that the light-absorbing and light-scattering aerosols have different effects on altering the PBL thermal structure. The inhibition effect of scattering aerosols on the PBL appears to be independent of the aerosol height distribution and solely depends on its concentration. However, the aerosol–PBL feedback of absorbing aerosols is highly dependent on its vertical distribution. Besides the 2 d case investigation, the analysis of the modeling results for nearly 1 month from 3 to 30 January 2020 in Baoding yields a more robust and representative conclusion. Our analysis highlights that we should principally concentrate on controlling the emissions of scattering aerosols under the stable stratification, while cooperating to control the emissions of scattering and absorbing aerosols in an unstable stratification. Moreover, the long-term interannual variation in the PBL stability shows a strong correlation with the East Asian winter monsoon, which seems to be valuable in determining which pollutants to target in different monsoon years and attaining more precise air pollution control. Based on the numerical simulations and observational constraints, a concept scheme description has been concluded to deepen our recognition of the interactions between thermodynamic stability and aerosols within the PBL over the NCP region.

1 Introduction

Ambient air pollution has been one of the major environmental issues in China, particularly in highly populated and industrialized areas, such as the North China Plain (NCP; Chan and Yao, 2008; Sun et al., 2014; Q. Zhang et al., 2019; Zheng et al., 2018; Fan et al., 2020; Luo et al., 2021). Severe and persistent air pollution episodes commonly break out in the presence of high aerosol emissions and unfavorable synoptic circumstances (e.g., low wind speed, high moisture, and stable stratification), posing great threats to human health (Che et al., 2019; X. Zhang et al., 2019). The interaction of aerosols with the planetary boundary layer (PBL), which is regarded to be a critical process for stabilizing the atmosphere and worsening surface air pollution, has been widely explored in the context of aerosol weather and climate effects (Z. Li et al., 2017; Su et al., 2020; Hung et al., 2021).

Aerosols contribute to considerable uncertainties in interpreting and quantifying the Earth's radiative budget and hydrological cycles through aerosol–radiation interactions (ARIs) and aerosol–cloud interactions (ACIs; Rosenfeld et al., 2014; Luo et al., 2019; Zhao et al., 2019; Letu et al., 2021). Despite the great progress in observational and numerical studies of ARIs and ACIs over recent decades, correctly quantifying the aerosol radiative effect (ARE) on the weather and climate systems remains a challenge. The principal reason for this challenge is the inadequate understanding of the strong variations in aerosol types, loadings, and vertical distributions, as well as the complex mechanisms among large-scale synoptic patterns, local-scale planetary boundary layer (PBL) structures, and AREs (H. Wang et al., 2015; Z. Li et al., 2017; Huang et al., 2018; Su et al., 2020).

The thermodynamic stability of the PBL dictates the planetary boundary layer height (PBLH), thereby dominating the vertical dissipation of surface pollutants to some extent (Q. Zhang et al., 2009; X. Ma et al., 2009; W. Zhang et al., 2018; Su et al., 2020). Aerosols, in turn, have vital feedback on the stability of the PBL, depending on their properties, particularly those of light-absorbing aerosols (e.g., black, brown, and organic carbon; Huang et al., 2018; Menon et al., 2002; Z. Wang et al., 2021). The feedback mechanisms between aerosol and PBL are critical in modulating air pollution outbreaks (Ding et al., 2013, 2016; Huang et al., 2018; Ma et al., 2020). This feedback process can be divided into two aspects, relating to surface cooling and atmospheric heating through ARE. On the one hand, the suspended aerosols in the atmosphere attenuate the incident solar radiation, leading to an overall cooling of the surface, diminishing surface sensible heat flux, and impeding the PBL development, which is known as the umbrella effect (Ma et al., 2020; Shen and Zhao, 2020). On the other hand, light-absorbing aerosols can trap solar energy and heat the atmosphere, which has different effects, depending on the aerosol vertical distribution. First, the absorbing aerosols in the low layer of the atmosphere strongly heat the near-surface layer to form a well-mixed PBL, which is recognized as the stove effect (Ma et al., 2020). Second, the absorbing aerosols in the upper layer of the atmosphere strengthen the PBL inversion intensity and weaken its development, which is identified as the dome effect (Ding et al., 2016).

In addition, the large-scale synoptic patterns disturb the sensitivity of the aerosol–PBL feedback mechanism. PBLH is not only influenced by the local-scale surface properties (surface cover types, sensible heat fluxes, latent heat fluxes, etc) but is also regulated by large-scale synoptic forcing (Hu et al., 2016; Miao et al., 2020). Synoptic conditions can alter the PBL thermodynamic stability through cold or warm advection, thus affecting the aerosol vertical dispersion (J. Zhang et al., 2012; Ye et al., 2016; Miao et al., 2019). Furthermore, due to the prevailing southerly winds, aerosols emitted in the southern regions can be carried to NCP via cross-region transport, and the increased aerosols eventually modulate the PBL thermodynamic stability (Du et al., 2020; Ma et al., 2020). On the interannual scale, previous studies reported that the wintertime air quality is closely connected with the East Asian winter monsoon (EAWM), and the strengthening (weakening) of the EAWM is typically associated with changes in the representative circulation patterns, which can improve (worsen) regional air quality (Niu et al., 2010; Y. Zhang et al., 2016).

According to the aforementioned efforts, most of them are primarily concentrated on the aerosol–PBL interactions and neglect the simultaneous influences of the synoptic condition, the aerosol type, and the vertical distribution. With numerical simulations and observational constraints, the roles of synoptic pattern, PBLH, aerosol type, and vertical distribution in aerosol–PBL interactions warrant further investigation. To address this issue, this study attempts to analyze the effects of synoptic forcing and aerosol types on the PBL thermodynamic stability, aerosol vertical distribution, and aerosol–PBL interaction by using in situ aircraft measurements, surface observations, and online coupled model Weather Research and Forecasting with Chemistry (WRF-Chem) simulations. On 3 and 4 January 2020, 2 d with distinct synoptic conditions, PBL stability and aerosol vertical distributions over NCP are analyzed. A total of five parallel numerical experiments are conducted to investigate the AREs of scattering and absorbing aerosols under different aerosol vertical distributions. Furthermore, to compensate for the 2 d case investigation, an analysis of the modeling results for nearly 1 month in Baoding from 3 to 30 January 2020 gives a more robust and representative conclusion. Moreover, the long-term interannual variation in the PBL stability from 1980 to 2020 over the NCP region is examined, and its driving factors are figured out. The conclusions appear to be beneficial in determining which pollutants to target in different weather conditions and achieving more precise controls of air pollution.

The remainder of this paper is organized as follows. Section 2 describes the data and methods used in this study. Section 3 presents and discusses the aerosol–PBL interactions, as well as the factors of synoptic forcing and aerosol types. The conclusions and summary are provided in Sect. 4.

2 Data and methods

2.1 Flight experiments and observational data

2.1.1 Intensity aircraft observation experiments

To study the vertical distributions of aerosol concentrations and meteorological conditions, the intensity aircraft observation experiments (aircraft model King Air 350, which is displayed in Fig. 1e) were carried out to measure the vertical profiles on 3 and 4 January 2020. The terrain height and three-dimensional flight routes are shown in Fig. 1 during the study periods. The flight routes on the 2 d were almost identical, which is conducive to the comparative analysis. The King Air 350 took off and ascended to about 3000 m over the Shahe airport (40.1^∘ N, 116.3^∘ E), which is located in Beijing. Afterward, the aircraft proceeded towards the southwest for about 150 km, arriving in Baoding and performing vertical measurements ranging from around 650 to 3000 m.

Figure 1 WRF-Chem domains with the terrain height and flight routes. (a) The outermost domain (d01) and the two embedded domains (d02 and d03). (b) The innermost domain (d03). (c–d) The three-dimensional flight routes on 3 and 4 January 2020, respectively.(e) The King Air 350 aircraft. The black dots and brown boxes in panel (b) show the locations of the air quality monitoring sites (AQMSs) and automatic weather stations (AWSs), respectively. The red and blue lines (overlain) in panel (b) represent the horizontal flight routes on 3 and 4 January 2020, respectively.

In this study, the aerosol number concentrations and size distributions were measured by a passive cavity aerosol spectrometer probe (PCASP) mounted on the wing tip of the King Air 350 aircraft. PCASP detected aerosol size distributions ranging from 0.1 to 3.0 µm at 1 Hz. Given the substantial noise of the first channel, we considered all the spectral bands except the first one, i.e., 0.11 to 3.0 µm, in this study (Twohy et al., 2005). A seven-wavelength (370, 470, 520, 590, 660, 880, and 950 nm) Aethalometer (AE-33) was used to detect aerosol absorption. Because black carbon (BC) is the predominant absorber at 880 nm (Ran et al., 2016), measurements at this wavelength were taken to represent BC concentration. In addition, an integrated meteorological measurement system (AIMMS) was employed to measure the aircraft location and ambient meteorological conditions (temperature, wind, etc.). The detailed information of the aircraft observation is summarized in Table 1.

Table 1 Information of the King Air 350 aircraft flight times and instruments. LST is local sidereal time.

2.1.2 Ground-based observations

Surface meteorological observations at three ground-based automatic weather stations (AWSs) in Beijing were collected for model meteorological validation, which includes the hourly datasets of surface pressure, 2 m temperature, 10 m wind speed, etc. In addition, air quality observations, including hourly ground-based PM_2.5 concentration measurements at six air quality monitoring sites (AQMSs) in Baoding, were used in this study to validate the model performance. The observational air pollutant data were collected from the China National Environmental Monitoring Centre (http://www.cnemc.cn/, last access: 20 November 2020). Ground-based observations on 3 and 4 January 2020 were applied in this study, and the locations of the AWSs and AQMSs are shown in Fig. 1b.

2.2 Reanalysis data

The data of isobaric temperature taken from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis dataset were used to summarize the long-term thermal structure of the PBL in the North China Plain. The temperature data of seven isobaric surfaces between 1000 and 850 hPa were considered to be within the PBL. The height between 1000 and 850 hPa is approximately equivalent to the height below 1.5 km. The sea level pressure from the ERA5 reanalysis product was used to calculate and describe the intensity of the EAWM and Siberian High (SiH), which are demonstrated to be the aspects that vigorously govern the PBL thermal structure over the NCP in this study. Here, the monthly average data during the wintertime (December, January, and February) from 1980 to 2020, with a grid of 0.25^∘ × 0.25^∘, was adopted.

2.3 WRF-Chem model

In this study, we simulated the aerosol concentrations and meteorological conditions by using the WRF-Chem model, version 3.9, which is an online, coupled, three-dimensional Eulerian chemical transport model considering complex physical and chemical processes (Grell et al., 2005). The three nesting domains with horizontal resolutions of 27, 9, and 3 km, respectively, are shown in Fig. 1a. There were 29 vertical layers extended from the ground to the top pressure of 50 hPa in the model, with more than 15 layers located below 3 km to fully describe the vertical structure of the PBL. The simulations were conducted from 12:00 UTC (universal coordinated time) on 25 December 2019 to 00:00 UTC on 31 January 2020, with the first 8.5 d as spin-up time for chemistry. The model was run with an 84 h model cycle, with the first 12 h discarded as spin-up time, and the last 72 h results used for the final analysis. The chemical outputs from previous runs were used as the initial conditions for the subsequent overlapping 84 h simulation. The simulations of the case study were carried out from 00:00 UTC on 2 January to 18:00  UTC on 4 January 2020, with the first 16 h as the model spin-up time, and the chemical outputs from the previous run were used as the initial conditions.

The 6 h National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) Final (FNL) reanalysis fields, with a spatial resolution of 1^∘ × 1^∘, were input for the model initial and lateral boundary meteorological conditions. Anthropogenic emissions were adopted from the Multi-resolution Emission Inventory for China (MEIC) in 2017, as developed by Tsinghua University (http://meicmodel.org, last access: 20 November 2020), whose emission sources were classified into the following five sectors: industrial process, power plants, residential combustion, on-road mobile sources, and agricultural activities (M. Li et al., 2017; Zheng et al., 2018). The Model of Emissions of Gases and Aerosols from Nature (MEGAN) was used online in the simulation to calculate the biogenic emissions (Guenther et al., 2006).

The model physical configurations we used in the simulation included the Morrison double-moment microphysics scheme (Morrison et al., 2005), the RRTMG (rapid radiative transfer model for general circulation models) radiation scheme (Iacono et al., 2008), the Yonsei University (YSU) boundary layer scheme (Noh et al., 2003), and the Noah land surface scheme (Ek et al., 2003). PBLH in the YSU scheme is determined from the Richardson bulk number method. The Model for Simulating Aerosol Interactions and Chemistry (MOSAIC) was applied as the sectional aerosol scheme, with aerosols specified as eight size sections (bins) ranging from 39 nm to 10 µm (Zaveri et al., 2008). The extinction, single-scattering albedo, and asymmetry factor of aerosols were computed as a function of wavelength and three-dimensional positions. Each chemical constituent of the aerosol was linked to complex indices of refraction. The refractive indices of the aerosols were calculated using volume averaging for each size bin, and the Mie theory was used to derive the extinction efficiency, the scattering efficiency, and the intermediate asymmetry factor. Aerosol optical properties were then determined by summarizing all size bins (Fast et al., 2006). The refractive indices of various aerosol components were reported in Barnard et al. (2010).

To investigate the interactions between aerosol and PBL thermodynamic stability over the NCP, five parallel numerical experiments were conducted, as presented in Table 2. The lapse rate below 1.5 km was used to quantify the PBL thermodynamic stability. In the control experiment (denoted as the EXP_Ctrl), the anthropogenic emissions were performed at the normal level, and the aerosol optical properties were calculated at each time step and then coupled with the short- and long-wave radiative transfer model. In the EXP_WoF, the emissions were the same as the EXP_Ctrl, except that ARE was not included. In the other three experiments, ARE was taken into consideration, but the emissions varied. Specifically, the BC emission was subtracted in the EXP_WFexBC, the BC emission was 20 times higher in the EXP_WF20BC, and the total aerosol emissions were strengthened to 20 times in the EXP_WF20Aer. To eliminate the influence of cloud condensation nuclei under various emission scenarios, the aerosol indirect effect was turned off in the modeling.

Table 2 Descriptions of five parallel WRF-Chem numerical experiments.

2.4 Statistical methods for model validation

Correlation coefficient (R), mean bias (MB), normalized mean bias (NMB), mean error (ME), and root mean square error (RMSE) were applied to assess the WRF-Chem model veracity in simulating meteorological parameters and air pollutants against the ground-based observations with the following equations (Itahashi et al., 2015; Granella et al., 2021):

$$\begin{array}{}\text{(1)}& {\displaystyle }R={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(S_i-\stackrel{\mathrm{‾}}{S})(O_i-\stackrel{\mathrm{‾}}{O})}{\sqrt{{\sum }_{i=\mathrm{1}}^n(S_i-\stackrel{\mathrm{‾}}{S}{)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^n(O_i-\stackrel{\mathrm{‾}}{O}{)}^{\mathrm{2}}}}},\text{(2)}& {\displaystyle }\mathrm{MB}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(S_i-O_i),\text{(3)}& {\displaystyle }\mathrm{NMB}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(S_i-O_i)}{{\sum }_{i=\mathrm{1}}^n\left(O_i\right)}},\text{(4)}& {\displaystyle }\mathrm{ME}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left|S_i-O_i\right|,\text{(5)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(S_i-O_i{)}^{\mathrm{2}}},\end{array}$$

where S_i and O_i are the simulated and observed parameters, respectively, n is the total number of the values used for validation, and $\stackrel{\mathrm{‾}}{S}$ and $\stackrel{\mathrm{‾}}{O}$ are the averages of the simulation and observation, respectively.

3 Results and discussion

3.1 Validations of meteorological parameters and aerosol concentrations

Considering that the weather conditions play important roles in pollutant transport, aerosol dispersion, and chemical reactions in the atmosphere, it is crucial to validate the simulated meteorological parameters with the observations. Figure 2 shows the validation of the meteorological parameters between the modeled and observed results, and model evaluations suggest that the WRF-Chem model is able to simulate the weather characteristics in the NCP region. The simulated temporal series and magnitudes of temperature, pressure and wind speed generally agree well with the ground-based and vertical aircraft observations. The accuracy of the simulation during the daytime is superior to that at night, which is instrumental in the intention of discussing the PBL thermodynamic stability and ARE in the daytime. In addition, validations of aerosol concentrations between the modeling and in situ observations are shown in Fig. 3. Both the simulation and observation display a high level of air pollution on 3 January and good air quality on 4 January. The vertical profiles of BC concentration suggest a good simulation performance, which can characterize the vertical variations and daily differences (Fig. 3a). The modeled surface PM_2.5 mass concentrations in Baoding compare well with the ground-based measurements, especially during the daytime (Fig. 3b).

Figure 2 Validation of meteorological parameters. (a–c) Time series of 2 m temperature, surface pressure, and 10 m wind speed in Beijing between the modeling (EXP_Ctrl) and ground-based observation, respectively. (d–e) Vertical profiles of temperature and wind speed in Baoding between the modeling (EXP_Ctrl) and aircraft observation, respectively. The shaded areas indicate the error bars (standard deviation).

Figure 3 Validation of aerosol concentration between the modeling (EXP_Ctrl) and in situ observations. (a) Aircraft-measured BC concentration vertical distributions. (b) Ground-based observed PM_2.5 concentration. The shaded areas indicate the error bars (standard deviation).

The statistical evaluations of model performance between in situ observations and EXP_Ctrl are concluded in Table 3. Underestimation of temperature exists in both the surface level and the vertical profile, though the R between the ground-based (aircraft) observation and EXP_Ctrl is 0.87 (0.98), indicating that the model describes the temperature variation well. The surface pressure is well simulated, with an MB of 0.81 and a small NMB of 0.08 %. A slight overestimation of 10 m wind speed (0.40 m s⁻¹) in the simulation might be due to the unresolved topography in the WRF model, which is also demonstrated in previous studies (Jiménez et al., 2013; M. Li et al., 2021). The wind speed profile suggests an underestimation, which could be attributable to flaws in aircraft wind speed measurement. The statistical validations of BC concentration vertical profiles show an R of 0.67 and a total MB of −0.18 µg m⁻³. The statistical validations of PM_2.5 mass concentration indicate an R of 0.79, a total MB of −4.91 µg m⁻³, and an NMB of 4.69 % during the daytime (Table 3). Therefore, in this study, we consider the WRF-Chem simulation to be in line with the observation, and that it can capture the weather characteristics and the general distributions and variations in air pollutants.

Table 3 Statistical evaluations of the WRF-Chem performance between in situ observations and EXP_Ctrl.

^a R – correlation coefficient; MB – mean bias; NMB – normalized mean bias; ME – mean error; RMSE – root mean square error. ^b Ground-based observations are from 10:00 to 18:00 LST on 3 and 4 January 2020. ^c Aircraft measurements are from 650 to 3000 m on 3 and 4 January 2020.

3.2 Aerosol vertical characteristics and PBL thermal structures

Despite the large-scale weather condition, the local-scale PBL thermal structure is also reported as being a critical aspect that substantially affects the surface aerosol loading (Q. Zhang et al., 2009; X. Ma et al., 2009). The PBL thermodynamic stability dictates the PBLH, thereby dominating the vertical dissipation of surface pollutants to some extent. Therefore, the detailed synoptic condition, PBL thermal structure, and their influences on the aerosol vertical characteristics are discussed in the following subsections.

3.2.1 Vertical thermal structure and synoptic condition

The temporal evolution of the EXP_Ctrl simulated temperature profile is presented in Fig. 4, which reflects the variation in the thermal structure and stability during the study period. The results indicate that the thermal structures differ significantly between the 2 d. It is clear that the temperature on 3 January is warmer than that on 4 January, suggesting different wintertime synoptic characteristics. Furthermore, the average lapse rate below 1.5 km is lower on 3 January, with lapse rate values of 3.62 and 5.38 ^∘C km⁻¹ during the aircraft observation periods on 3 and 4 January, respectively.

Figure 4 Temporal evolution of the EXP_Ctrl simulated temperature profile over Baoding. The pink and blue dashed boxes correspond to the aircraft observation periods on 3 and 4 January, respectively. The black line indicates the planetary boundary layer height (PBLH).

The thermal structure of the atmosphere is closely related to and highly sensitive to the synoptic condition over the NCP, especially during the wintertime (Yang et al., 2016; Miao et al., 2019). Figure 5a and b display the surface synoptic patterns during the aircraft observation periods on 3 and 4 January, respectively. In the adjacent regions of Baoding, the wind fields are practically opposite in these two periods, and the contrast is also intuitively shown in Fig. 5c, with the vertical profile difference in the meridional wind speed. The prevailing surface wind direction on 3 January is southwesterly, while on 4 January is northeasterly. The southerly wind at the low layer warms the atmosphere below 1.5 km on 3 January, whereas the northerly wind leads to a cooling effect on 4 January. Thus, there appears a large temperature difference manifested in Fig. 4. Furthermore, the southerly wind warms the atmosphere, particularly the aloft layer at about 1 km altitude, leading to a stable thermal stratification below 1 km and a low PBLH on 3 January.

Figure 5 The EXP_Ctrl simulated synoptic patterns during the study period. (a–b) The surface wind, temperature, and pressure field during the aircraft observation periods on 3 and 4 January, respectively. (c) The temporal evolution of the meridional wind speed over Baoding. The green dots in panels (a)–(b) denote the location of Baoding. The pink and blue dashed boxes in panel (c) correspond to the aircraft observation periods on 3 and 4 January, respectively.

3.2.2 Observed and simulated aerosol vertical distributions

As demonstrated in Fig. 6, contrasting aerosol vertical distributions are observed with aircraft during the afternoon of 3 and 4 January 2020, and the specific times are shown in Table 1. The accumulation mode particles (0.1–1 µm) dominate the aerosol number concentrations during the study period, which is primarily emitted by anthropogenic pollutants, such as BC aerosol loadings (Y. Wu et al., 2017). On 3 January, larger amounts of aerosol particles are constrained to the low layer of the atmosphere (below 650 m), while the pollutants decrease sharply above 650 m and gradually tend to a stable low value. On 4 January, the aerosol number concentrations below 650 m are lower than that on 3 January but higher when the height is above 1000 m. This contrast is attributed to the differences in atmospheric stability (Su et al., 2020). The aloft aerosol layer between 1.0 and 1.7 km altitude on 4 January is related to the PBL vertical mixing transport. The great difference in the lapse rate shown in Fig. 4 leads to the disparity in atmospheric stability and aerosol dispersion ability. On 4 January, the more unstable thermal stratification promotes the dispersion of the near-surface pollutants and, meanwhile, exacerbates the pollution on the elevated layer.

Figure 6 Aircraft observations of aerosol vertical distributions in the afternoon during the flight. (a) The aerosol number concentrations on 3 and 4 January. The dashed and solid lines denote the observations over Beijing and Baoding, respectively. (b–c) The aerosol size distributions on 3 and 4 January, respectively. The shaded areas in panel (a) indicate the error bars (standard deviation). The horizontal black dashed lines are the reference altitude of 650 m.

The temporal evolution of the EXP_Ctrl simulated vertical profiles of total aerosol number concentrations and BC mass concentrations are shown in Fig. 7. Both profiles have similar features, indicating that anthropogenic fine particles, like BC, are the dominant pollutants, which is consistent with the particle size observations exhibited in Fig. 6b and c. The diurnal variation in the aerosol profile is strongly driven by the evolution of PBLH, which is well examined in the previous studies (Qu et al., 2017; Huang et al., 2018). The pollutants are well constrained to the near-surface by the stable PBL at nighttime, while the vertical distribution of aerosol particles is more homogenous by turbulent mixing in the afternoon. In addition, on 4 January, the vertical mixing of aerosols is more homogenous, transporting aerosols to a higher layer, which implies a more unstable thermal stratification and a higher PBLH.

Figure 7 Temporal evolution of the EXP_Ctrl simulated aerosol vertical distributions over Baoding. (a) The aerosol number concentration (particle diameter of 0.15–2.5 µm). (b) The BC concentration. The pink and blue dashed boxes correspond to the aircraft observation periods on 3 and 4 January, respectively.

The contrasting aerosol vertical distributions between the 2 d are interpreted from two aspects. From the perspective of aerosol regional transport, the air quality is remarkably altered by the synoptic condition (Y. Zhang et al., 2016; An et al., 2019). The accumulation of pollutants over NCP is contributed by regional transport from the polluted southwestern region via the northward movement of the air mass, whereas the cessation of pollution over NCP is caused by prevailing northern air mass, and then the pollution hotspot moves to the south (Luo et al., 2021). Therefore, the observed and simulated high surface pollution loading on 3 January is externally affected by the unfavorable synoptic condition, which is consistent with the southerly wind demonstrated in Fig. 5. From the perspective of aerosol vertical distribution, the surface air pollution loading is efficiently driven by the thermal structure and stability of the PBL (Miao et al., 2020; Q. Li et al., 2021). On 3 January, the prevailing surface southerly wind leads to a stable stratification, restraining aerosol vertical dispersion and, finally, exacerbating the surface air pollution. Besides, the northerly wind condition on 4 January has the opposite effect, which reduces the surface aerosol loading.

In conclusion, the stable stratification related by the prevailing surface southerly wind limits the PBL development and exacerbates the surface air pollution by restraining the aerosols to the low layer. Conversely, the unstable stratification generated by surface northerly wind promotes the PBL turbulent mixing and lightens the surface aerosol concentration via vertical diffusion.

3.3 Aerosol and PBL thermodynamic stability interactions

3.3.1 ARE on PBL thermal structure

PBL structure modulates the aerosol vertical distribution by turbulent mixing, and inversely, the suspended aerosols may also modify the PBL thermal structure by ARE to some extent. Figure 8 presents the temperature variation affected by ARE, and it is noticeable that light-absorbing aerosols heat the atmosphere, while light-scattering aerosols contribute to a cooling effect. On 3 January, aerosols are constrained to the low layer during the daytime, and the overall ARE results in a warming effect below 1 km. The conspicuous heating suggests the leading role of absorptive BC aerosols in the variation in PBL thermal structure when coexisting with other types of scattering aerosols on 3 January. On 4 January, due to the strong turbulence mixing (as demonstrated by the temperature profile in Fig. 4), aerosols are carried to the aloft layer. The overall ARE results in a cooling effect below 0.6 km and a warming effect between 0.6 and 1.5 km. When the aerosols are transported to a high altitude, the absorbing and scattering particles inhibit the incident solar radiation from reaching the low layer and exhibit a cooling effect, while the absorbing aerosols heat the upper layer. The contrasting aerosol vertical distribution caused by the varied PBL thermal structure between the two daytime periods results in different AREs; therefore, it is quite crucial to examine the aerosol vertical structure when evaluating the ARE.

Figure 8 Temporal evolution of temperature profile variation influenced by the aerosol radiative effect (ARE). (a) ARE by all the aerosols (EXP_Ctrl – EXP_WoF). (b) ARE by BC (EXP_Ctrl – EXP_WFexBC). (c) ARE by other aerosols (EXP_WFexBC – EXP_WoF). The pink and blue dashed boxes correspond to the aircraft observation periods on 3 and 4 January, respectively.

The simulated ARE modify the thermal structure with a temperature variation of less than 0.5 ^∘C (shown in Fig. 8). To amplify the signal of ARE when discussing its interactions with PBL, we further performed two experiments, i.e., one with strong overall aerosol loading (EXP_WF20Aer) and the other only with strong BC aerosol loading (EXP_WF20BC). The simulated thermal structures are presented in Fig. 9. The amplified simulations reveal a temperature variation of less than 2.5 ^∘C, which is consistent with the altered magnitude of aerosol emissions. ARE by only absorptive BC aerosols heat the low layer more remarkably when compared with the role of overall aerosols, suggesting that the coexistence of absorbing and scattering aerosols partially conceal the signal of ARE. This result is owing to the opposite optical roles of the two aerosol types in the low layer and the strong inhibition of the incident solar radiation by the aerosols in the upper layer.

Figure 9 Temporal evolution of temperature profile with (a) high aerosol emissions (EXP_WF20Aer) and (c) high BC emissions (EXP_WF20Aer). The temporal evolution of the temperature variation profile influenced by (b) high scattering aerosol loadings (EXP_WF20Aer – EXP_WF20BC) and (d) high absorbing aerosol loadings (EXP_WF20BC – EXP_Ctrl). The pink and blue dashed boxes correspond to the aircraft observation periods on 3 and 4 January, respectively.

3.3.2 ARE on the PBLH, lapse rate, and surface pollution

Since ARE can modulate the PBL thermal structure, the related impacts caused by the thermal variation are discussed here. Figure 10 shows the time series of the simulated lapse rate, PBLH, and surface carbon monoxide (CO) concentration under five parallel experiments.

Generally, on 3 January, the PBLH is lower, and the surface air pollution is heavier than that on 4 January. Influenced by the southwesterly winds, the warmer air mass can be guided to the NCP (Fig. 5a), enhancing the thermal stability and restraining the growth of the PBL (Fig. 10b). On the one hand, from the local-scale view, the lower PBLH constrains the air pollution to the near-surface; on the other hand, from the regional transport view, the southwesterly prevailing winds can transport the pollutants emitted from upstream plain regions to the NCP, which further worsen the surface air quality (Fig. 10c). However, the relatively higher PBLH and the northeasterly winds on 4 January are able to reduce the surface air pollution and inversely carry the pollutants to the higher altitude with the strong turbulent mixing. The disparate synoptic conditions regulate the vertical aerosol structures and further affect the ARE and other related variables.

The distinct aerosol types contribute to varied AREs when different aerosol vertical distributions are dominant. The lapse rate and PBLH are closely associated with the vertical thermal structure. A warmer low layer and a colder upper layer induce a higher lapse rate and PBLH, while a colder low layer and a warmer upper layer decrease the lapse rate and PBLH. On 3 January, the pollutants are restricted to the low layer, but the absorbing and scattering aerosols play the opposite roles. Ma et al. (2020) have confirmed that the radiative heating of high concentration absorbing aerosols near the surface will offset and exceed the reduction in the surface sensible heat flux, resulting in near-surface heating (stove effect). The low layer heating, together with the surface buoyancy flux, invigorates the turbulent mixing and leads to a higher lapse rate and PBLH (Ma et al., 2020). On the contrary, the scattering aerosols cool the low layer and inhibit the growth of the PBL (umbrella effect), causing a lower lapse rate and PBLH (Fig. 10a and b). On 4 January, the aerosol vertical profile differs from that on 3 January, with a huge number of pollutants being carried to the high altitude by the unstable stratification. The absorbing and scattering aerosols both prevent the solar radiation to the low layer (umbrella effect), while the absorbing aerosols also heat the upper layer (dome effect), generating a lower lapse rate and PBLH (Fig. 10a and b).

The variations in lapse rate and PBLH further alter the surface air pollution by the thermodynamic and dynamic effect, and the high PBLH ordinarily lightens surface air pollution. Due to CO being a relatively stable pollutant that hardly interacts with other pollutants, and due to the main sources of CO being primary emissions from anthropogenic pollution and open biomass burning (Chi et al., 2013), we use the CO concentration to evaluate the variation in surface air pollution caused by the ARE in this study. The surface CO concentration in the EXP_WF20Aer is divided by a factor of 20 to match other experiments, since the emissions are magnified 20 times. The results presented in Fig. 10c indicate that the surface air pollution variation is closely related to the PBLH variation and driven by the aerosol vertical dispersion. When the stable stratification confines the absorbing aerosols to the low layer, the stove effect destroys the stability, enhancing the thermal convection and decreasing the surface air pollution. The absorbing aerosols in the upper layer or the scattering aerosols in both the low and upper layers stabilize the atmosphere, deteriorating the surface air quality.

Figure 10 Temporal evolution of the (a) lapse rate within 1.5 km, (b) the planetary boundary layer height (PBLH), and (c) the surface CO concentration simulated by the five parallel experiments. The upward or downward arrows indicate an increase or decrease compared to the EXP_Ctrl.

3.3.3 Statistical properties of the PBL and AREs under different synoptic conditions

It is noticeable that different aerosol vertical distributions between the 2 d contribute to distinct AREs due to the synoptic condition and PBL thermal stability differences from the measurements and simulations. In particular, the absorptive BC aerosols have both stove and dome effects, which affect the PBL thermal structure. Here, we further analyze the modeling results for nearly 1 month, from 3 to 30 January 2020, in Baoding to give a more significant and representative conclusion.

Figure 11 shows the correlations between the daily average 10 m meridional wind speed, lapse rate within 1.5 km, and PBLH. The negative correlation between the 10 m meridional wind speed and the lapse rate within 1.5 km (Fig. 11a) suggests that the increased southerly wind stabilizes the PBL, whereas the strong northerly wind destabilizes the PBL. The variation in the lapse rate has a direct impact on the development of the PBL, as evidenced by the PBLH modification shown in Fig. 11b. Figure 12 compares the distinct vertical distributions of aerosols caused by northerly and southerly winds. Samples with a daily average wind speed within ± 0.05 m s⁻¹ are discarded to avoid the north–south reversal of wind direction in a day. Eventually, 16 d with the prevailing northerly wind and 8 d with the prevailing southerly wind are used. The result indicates that the synoptic condition influences the PBL thermal structure, thereby affecting the vertical dispersion of aerosols. The warm and polluted air is carried to the NCP by the southerly winds, which stabilize the PBL, exacerbating the surface air pollution. The cold and clean air is carried to the NCP by the northerly winds, forming an unstable stratification and transporting pollutants to the upper layer. The previous study has also verified the conclusion that the wind field and PBLH tend to be correlated, which means meteorologically favorable horizontal and vertical dispersion conditions are likely to occur together (Su et al., 2018).

Figure 11 Scatterplots of the correlations between the (a) 10 m meridional wind speed (positive – southerly wind; negative – northerly wind) and lapse rate within 1.5 km and (b) the lapse rate within 1.5 km and PBLH. The data are daily averages from 3 to 30 January 2020.

Figure 12 Vertical distributions of the aerosol number concentrations (particle diameter of 0.15–2.5 µm) under the prevailing southerly wind and northerly wind, respectively. The shaded areas indicate the error bars (standard deviation).

When evaluating the AREs of light-absorbing and light-scattering aerosols, the temperature profile variations show various patterns due to differences in aerosol concentration and vertical distribution caused by synoptic conditions and, particularly, the wind direction. Figure 13 demonstrates that light-absorbing aerosols heat the atmosphere, while light-scattering aerosols contribute to a cooling effect. Aerosols are constrained to the low layer under southerly wind conditions, and the BC aerosols result in a warming effect below 1 km (stove effect), while the scattering aerosols cool the layer below 0.6 km (umbrella effect). In contrast, the PBLs exhibit strong turbulence mixing when influenced by the northerly winds, and the aerosols are carried to the aloft layer. The aloft scattering particles prevent incident solar radiation from reaching the low layer, resulting in cooling effects below 1 km (umbrella effect), whereas the aloft absorbing aerosols heat the upper layer between 0.5 and 1.5 km (dome effect). The remarkable aerosol effects under southerly winds attribute to the accumulation of aerosols under adverse weather conditions. The contrasting aerosol vertical distributions caused by the varying synoptic conditions lead to different AREs, which is consistent with the results obtained on 3 and 4 January 2020.

Figure 13 Temporal evolution of the temperature profile variation influenced by aerosol radiative effect (ARE). (a) ARE by other aerosols (EXP_WFexBC – EXP_WoF) and (b) ARE by BC (EXP_Ctrl – EXP_WFexBC) under the prevailing southerly wind. (c) ARE by other aerosols and (d) ARE by BC under the prevailing northerly wind.

Furthermore, based on the nearly 1-month simulations, we quantify the variations in the lapse rate within 1.5 km and PBLH under different synoptic conditions, respectively, as caused by absorptive BC aerosols and other light-scattering aerosols. The results in Fig. 14 reveal that the BC stove effect induces a 0.04 ^∘C km⁻¹ increase in the lapse rate within 1.5 km and a 3 m increase in PBLH under the stable stratifications with the prevailing southerly winds. However, the BC dome effect causes a 0.085 ^∘C km⁻¹ decrease in lapse rate within 1.5 km and a 3 m decrease in PBLH under the unstable stratifications with the prevailing northerly winds. The umbrella effect of scattering aerosols in both stable and unstable conditions reduces the lapse rate by about 0.15 ^∘C km⁻¹ and reduces the PBLH by about 3.5–4 m. The vertical distribution of absorbing aerosols has a significant impact on their aerosol–PBL feedback. The absorbing aerosols concentrated in the low layer that have a strong radiative heating effect on the atmosphere develop in the PBL in the case of stable weather patterns under the influence of the southerly wind. The absorbing aerosols in the upper layer heat the atmosphere and inhibit the development in the PBL in the case of unstable weather patterns under the influence of the northerly wind. The inhibition effect of scattering aerosols on the PBL, on the other hand, is independent of the aerosol height distributions and is solely dependent on the aerosol concentrations.

Figure 14 Box plots of the variations in lapse rate within 1.5 km and PBLH influenced by the BC and other aerosols under (a) the prevailing southerly wind and (b) the prevailing northerly wind, respectively. The squares represent the mean values, the horizontal lines inside the boxes are the medians, and the bottom and top sides of the boxes represent the first and third quartiles. The whiskers are the minimum (maximum) values within 1.5 interquartile ranges of the lower (upper) quartile. The asterisks indicate the minimum (maximum) values.

3.4 Long-term variation in PBL thermodynamic stability over NCP

On a synoptic scale, the PBL stability is influenced by the meridional wind direction to a certain extent, which has been evidenced by the 2 d detailed aircraft measurements and 1-month model simulations. Furthermore, the link between the synoptic pattern (meridional wind intensity) and PBL thermal structure is analyzed from a climatological view, which is critical for the examination of aerosol vertical distribution and ARE. The potential impact of the interannual climate variability on the PBL thermodynamic stability is explored in this study over the NCP during the wintertime (December, January, and February) of 1980–2013. In fact, the interannual variabilities in the PBL are affected by numerous factors, and we only discuss the influence of meridional wind intensity in this study, which is a climatological extension of the case study. Following the method developed by Wu and Wang (2002), the index of SiH and the index of EAWM are used to analyze the possible factors that influence the PBL thermal structure over the NCP. The index of SiH and EAWM are the two variables describing the meridional wind intensity over the coast of East Asia, which are crucial to the PBL thermal structure. The PBL lapse rate is calculated in domain d03 (depicted in Fig. 1b) between 1000 and 850 hPa.

Figure 15 suggests that the relationships between SiH index, EAWM index, and PBL lapse rate are highly correlated. The results indicate that the interannual variabilities in the EAWM and SiH, which are indicators of the meridional wind intensity, are closely related to the PBL thermodynamic stability over the NCP region. Weak EAWM and SiH, in conjunction with the prevailing south wind, decrease the lapse rate between 1000 and 850 hPa and stabilize the PBL, while the strong EAWM and SiH play the opposite role. According to the above synoptic analysis, it is concluded that the changes in the PBL properties will further affect the surface air quality, aerosol vertical distribution, and ARE. As a result, the EAWM and SiH indirectly influence the aerosol vertical distribution by modulating the PBL thermodynamic structure.

4 Conclusions and summary

In this study, the complex relationships among the large-scale synoptic patterns, local PBL thermal structures, aerosol vertical distributions, and AREs of different aerosol types are investigated by combining aircraft observations, surface measurements, reanalysis data, and WRF-Chem simulations. The in situ aircraft observations were carried out over Beijing and Baoding during the daytime on 3 and 4 January 2020, and the numerical simulations were performed over NCP from 3 to 30 January. By and large, the meteorological variables and air pollutants are well simulated when validated by the observations.

Observations show that the surface air pollution over Baoding on 3 January is heavier than that on 4 January, and the aerosols are constrained to the low layer on 3 January, while the aerosols mix more homogeneously and vertically on 4 January. The numerical simulations suggest that the synoptic pattern over Baoding differs between the 2 d, and the prevailing wind direction is opposite, with a southwesterly wind on 3 January and a northeasterly wind on 4 January. Through the 2 d detailed measurement and a long-term simulation, we conclude that the synoptic condition may affect the PBL thermal structure, thus affecting the vertical dispersion of aerosols. The southerly winds flow the warm and polluted air to the NCP and stabilize the PBL, hence exacerbating the surface air pollution. The northerly winds flow the cold and clean air to the NCP, establishing an unstable stratification and transporting the pollutants to the upper layer.

The suspended aerosols may also modify the PBL thermodynamic stability through ARE to some extent, and the synoptic condition can modulate the sensitivity of PBL–aerosol feedback by influencing the PBL thermal structure and aerosol vertical distribution. The sensitive numerical experiments reveal that the absorbing and scattering aerosols have different effects when modifying the PBL thermal structure, and the AREs are also influenced by the aerosol vertical distribution. The aerosol–PBL feedback of absorbing aerosols is highly dependent on its vertical distribution. In the case of stable weather patterns under the influence of the southerly wind, the absorbing aerosols concentrated in the low layer have a strong radiative heating effect on the atmosphere (stove effect). This stove effect disturbs the stratification stability of the low layer, promotes the development of the PBL in coordination with the surface sensible heat flux, and transports heat and pollutants to the upper layer through turbulence. In the case of unstable weather patterns under the influence of the northerly wind, the absorbing aerosols in the upper atmosphere heat the atmosphere (dome effect) and inhibit the development of the PBL, which further leads to the accumulation of surface pollutants. The inhibition effect of scattering aerosols on the PBL is independent of the aerosol height distribution but only depends on its concentration. Scattering aerosols reflect solar radiation back to the space (umbrella effect), thus weakening the surface sensible heat flux and upward heat transfer, inhibiting the development of the PBL and further accumulating surface pollutants. In summary, the umbrella, stove, and dome effects of the scattering and absorbing aerosols under different synoptic patterns are illustrated with four scenarios in Fig. 16, which correspond to the conclusion indicated in the four subgraphs of Fig. 13. The conclusions indicate that we should mainly control the emissions of scattering aerosols under the stable stratification and cooperate to control the emissions of scattering and absorbing aerosols in an unstable stratification, especially to prevent the accumulation of surface pollutants caused by the dome effect of absorbing aerosols.

Figure 15 (a) Time series of the standardized anomaly of the wintertime boundary layer lapse rate (LR) in domain d03 between 1000 and 850 hPa, the index of the Siberian High (SiH), and the index of East Asian winter monsoon (EAWM) from 1980 to 2020. The scatterplots of the correlations between (b) the standardized anomaly of the EAWM and LR and between (c) the standardized anomaly of the SiH and LR. The standardized anomaly is calculated by dividing anomalies by the climatological standard deviation.

Figure 16 The illustration of the umbrella, stove, and dome effects of the scattering and absorbing aerosols under different synoptic patterns. The red plus and blue minus signs indicate the aerosol heating and dimming effect, respectively. The upward and downward arrows denote the increase and decrease in PBLH caused by ARE, respectively.

In addition, the factors affecting the wintertime PBL stability over the NCP area are further analyzed from the climatological perspective, which is mainly controlled by the intensity of EAWM and SiH. A strong winter monsoon creates an unstable PBL stratification, whereas a weak winter monsoon stabilizes the PBL. Combined with the umbrella, stove, and dome effects of scattering and absorbing aerosols under different synoptic conditions, this finding aids us in determining which pollutants to target in different monsoon years and achieving more precise air pollution control.

It should be noted that only 2 d of pollution episodes are measured in this study. While the model simulations produce evidence that is compatible with the observations and substantiate the observational analysis, additional measurements should be carried out in the near future.