Aerosol–radiation feedback deteriorates the wintertime haze in the North China Plain

Atmospheric aerosols scatter or absorb a fraction of the incoming solar radiation to cool or warm the atmosphere, decreasing surface temperature and altering atmospheric stability to further affect the dispersion of air pollutants in the planetary boundary layer (PBL). In the present study, simulations during a persistent and heavy haze pollution episode from 5 December 2015 to 4 January 2016 in the North China Plain (NCP) were performed using the Weather Research and Forecasting model with Chemistry (WRF-Chem) to comprehensively quantify contributions of aerosol shortwave radiative feedback (ARF) to near-surface (around 15 m above the ground surface) PM_2.5 mass concentrations. The WRF-Chem model generally performs well in simulating the temporal variations and spatial distributions of air pollutants concentrations compared to observations at ambient monitoring sites in the NCP, and the simulated diurnal variations of aerosol species are also consistent with the measurements in Beijing. Additionally, the model simulates the aerosol radiative properties, the downward shortwave flux, and the PBL height against observations in the NCP well. During the episode, ARF deteriorates the haze pollution, increasing the near-surface PM_2.5 concentrations in the NCP by 10.2 µg m⁻³ or with a contribution of 7.8 % on average. Sensitivity studies have revealed that high loadings of PM_2.5 attenuate the incoming solar radiation reaching the surface to cool the low-level atmosphere, suppressing the development of the PBL, decreasing the surface wind speed, further hindering the PM_2.5 dispersion, and consequently exacerbating the haze pollution in the NCP. Furthermore, when the near-surface PM_2.5 mass concentration increases from around 50 to several hundred µg m⁻³, ARF contributes to the near-surface PM_2.5 by more than 20 % during daytime in the NCP, substantially aggravating the heavy haze formation. However, when the near-surface PM_2.5 concentration is less than around 50 µg m⁻³, ARF generally reduces the near-surface PM_2.5 concentration due to the consequent perturbation of atmospheric dynamic fields.

1 Introduction

Atmospheric aerosols, produced both naturally and anthropogenically, influence the radiative energy budget of the Earth's atmospheric system in many ways. They scatter or absorb a fraction of the incoming solar radiation to cool or warm the atmosphere, decreasing surface temperature and altering atmospheric stability (e.g., Ackerman, 1977; Jacobson, 1998, 2002). Also, they serve as cloud condensation nuclei (CCN) and ice nuclei (IN), thus modifying cloud optical properties and lifetime (e.g., Zhang et al., 2007; Li et al., 2008a, b, 2009). Among those impacts, the scattering and absorption of solar radiation by aerosols and the associated feedbacks (hereafter referred to as aerosol–radiation feedback or ARF) not only constitute one of the main uncertainties in climate prediction (IPCC, 2007), but also substantially affect the atmospheric chemistry by perturbing the temperature profile and moisture, winds, and planetary boundary layer (PBL) stability (Boucher et al., 2013). In particular, as a short-lived pollutant with uneven distribution and physical and chemical heterogeneities in the atmosphere, ARF varies by more than a factor of 10 with location or time of emissions (Penner et al., 2010).

During wildfire with high-loading absorbing aerosols, ARF has been reported to heat the atmosphere and cool the surface, and thence enhance the PBL stability (e.g., Grell et al., 2011; Fu et al., 2012; Wong et al., 2012). In addition, numerous studies have been performed to evaluate impacts of ARF of dust on the regional meteorology and climate (e.g., Perez et al., 2006; D. F. Zhang et al., 2009; Santese et al., 2010). Anthropogenic aerosols, dominated by scattering components, such as organics and sulfate, primarily attenuate the incoming solar radiation down to the surface, cooling the temperature of the low-level atmosphere to suppress the development of the PBL and hinder the aerosol dispersion in the vertical direction (e.g., Fast et al., 2006; Vogel et al., 2009; Zhang et al., 2010). In addition, the temperature profile perturbation caused by ARF also alters cloud formation and development, possibly causing the precipitation delay or decrease (e.g., Zhao et al., 2005; Koch and Del Genio, 2010; Ding et al., 2013).

Rapid industrialization and urbanization in China have significantly elevated the concentrations of aerosols or fine particulate matter (PM_2.5), causing frequent occurrence of haze pollution, particularly during wintertime in North China (e.g., Zhang et al., 2013; Pui et al., 2014). Guo et al. (2014) have elucidated the haze formation mechanism in China, highlighting the efficient aerosol nucleation and growth during haze episodes. Moreover, high-loading aerosols during heavy haze episodes induce efficient ARF, encumbering the PBL development and further deteriorating the haze pollution. It is worth noting that ARF increases precursors for the aerosol nucleation and growth in the PBL, such as sulfuric and organic gases, causing efficient aerosol nucleation and growth (Zhang et al., 2004; Guo et al., 2014). Based on field measurements, recent studies have proposed that the high level of PM_2.5 increases the stability of PBL due to ARF and further decreases the PBL height (PBLH), consequently enhancing PM_2.5 concentrations ([PM_2.5]) (Quan et al., 2013; Petäjä et al., 2016; Yang et al., 2016; Tie et al., 2017; Ding et al., 2017). Online-coupled meteorology and chemistry models have also been used to verify the impact of ARF on the PBLH and near-surface [PM_2.5] during haze episodes in Europe, eastern China, and northern China (Forkel et al., 2012; Z. F. Wang et al., 2014; Wang et al., 2015; Zhang et al., 2015; Gao et al., 2015). However, the ARF impact on near-surface [PM_2.5] varies, depending on the evaluation time and location (Table 1). For example, the two-way coupled WRF–CMAQ (Community Multiscale Air Quality) system has been employed to evaluate the ARF contribution to the haze formation in January 2013 over the North China Plain (NCP), showing that ARF reduces the PBLH by 100 m and enhances near-surface [PM_2.5] by up to 140 µg m⁻³ in Beijing (J. Wang et al., 2014). Therefore, it is still imperative to comprehensively quantify the ARF contribution to near-surface [PM_2.5] under various pollution levels to provide the underlying basis for supporting the design and implementation of emission control strategies.

Table 1 Impact of ARF on near-surface [PM_2.5] in China.

In this study, simulations are performed using the Weather Research and Forecasting model with Chemistry (WRF-Chem) to interpret the relationship between the near-surface [PM_2.5] and the PBLH and further quantify the ARF contribution to near-surface [PM_2.5] under various pollution levels. The model and methodology are described in Sect. 2. Analysis results and discussions are presented in Sect. 3, and summary and conclusions are given in Sect. 4.

2 Model and methodology

2.1 WRF-Chem model and configurations

The WRF-Chem model (Grell et al., 2005) with modifications by Li et al. (2010, 2011a, b, 2012) is applied to evaluate effects of ARF on the wintertime haze formation in the NCP. The model includes a new flexible gas-phase chemical module, which can be used with different chemical mechanisms, such as the Carbon Bond Mechanism (CBIV), Regional Acid Deposition Mechanism 2 (RADM2), and Statewide Air Pollution Research Center mechanism (SAPRC). In the study, the SAPRC99 chemical mechanism is used based on the available emission inventory. For the aerosol simulations, the CMAQ/models3 aerosol module (AERO5) developed by the US EPA has been incorporated into the model (Binkowski and Roselle, 2003). The wet deposition is based on the method in the CMAQ module, and the dry deposition of chemical species follows Wesely (1989). The photolysis rates are calculated using the FTUV (fast radiation transfer model) with the aerosol and cloud effects on photolysis (Li et al., 2005, 2011a).

It is worth noting that the most recent extension of ISORROPIA, known as ISORROPIA II, has incorporated a larger number aerosol species (Ca, Mn, K salts) and is designed to be a superset of ISORROPIA (Fountoukis et al., 2009). However, the ISORROPIA version II uses the exact same routines as ISORROPIA to compute the equilibrium composition, which produces identical results to ISORROPIA when crustal species are not considered. Therefore, the inorganic aerosols in this study are predicted using ISORROPIA version 1.7, calculating the composition and phase state of an ammonium–sulfate–nitrate–water inorganic aerosol in thermodynamic equilibrium with gas-phase precursors (Nenes et al., 1998). In addition, a parameterization of sulfate heterogeneous formation involving aerosol liquid water (ALW) has been developed and implemented into the model, which has successfully reproduced the observed rapid sulfate formation during haze days (Li et al., 2017a). The sulfate heterogeneous formation from SO₂ is parameterized as a first-order irreversible uptake by ALW surfaces, with a reactive uptake coefficient of $\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{4}}$ assuming that there is enough alkalinity to maintain the high iron-catalyzed reaction rate.

The OA module is based on the volatility basis set (VBS) approach, with aging and detailed information found in Li et al. (2011b). The POA components from traffic-related combustion and biomass burning are represented by nine surrogate species with saturation concentrations (C^*) ranging from 10⁻² to 10⁶ µg m⁻³ at room temperature (Shrivastava et al., 2008) and assumed to be semivolatile and photochemically reactive (Robinson et al., 2007). The SOA formation from each anthropogenic or biogenic precursor is calculated using four semivolatile VOCs with effective saturation concentrations of 1, 10, 100, and 1000 µg m⁻³ at 298 K. The SOA formation via the heterogeneous reaction of glyoxal and methylglyoxal is parameterized as a first-order irreversible uptake by aerosol particles with an uptake coefficient of $\mathrm{3.7}\times {\mathrm{10}}^{-\mathrm{3}}$ (Liggio et al., 2005; Zhao et al., 2006; Volkamer et al., 2007).

A persistent air pollution episode from 5 December 2015 to 4 January 2016 in the NCP is simulated using the WRF-Chem model. During the study episode, the average hourly [PM_2.5] in the NCP are approximately 127.9 µg m⁻³, within the fourth grade of National Ambient Air Quality Standards with [PM_2.5] between 115 and 150 µg m⁻³ (moderately polluted, Feng et al., 2016). The persistent and widespread haze pollution episode with high [PM_2.5] in the NCP provides a suitable case for observation analyses and model simulations to investigate the ARF effect on haze pollution. Figure 1a shows the model simulation domain, and detailed model configurations can be found in Table 2.

Figure 1 (a) WRF-Chem simulation domain with topography and (b) the Beijing–Tianjin–Hebei area. In (a), the blue circles represent centers of cities with ambient monitoring sites, and the size of blue circles denotes the number of ambient monitoring sites of cities. In (b), the blue and red filled circles denote the NCNST and IRSDE site, respectively; the red filled rectangle denotes the meteorological site. The red numbers denote the CERN sites with the solar radiation measurement. (1) Beijing urban; (2) Jiaozhouwan; (3) Yucheng; (4) Luancheng.

Table 2 WRF-Chem model configurations.

2.2 Aerosol radiative module

In the present study, the Goddard shortwave module developed by Chou and Suarez (1999, 2001) is employed to take into account the ARF effect on the haze formation. The aerosol radiative module developed by Li et al. (2011b) has been incorporated into the WRF-Chem model to calculate the aerosol optical depth (AOD or τ_a), single scattering albedo (SSA or ω_a), and the asymmetry factor (g_a).

In the CMAQ aerosol module, aerosols are represented by a three-moment approach with a lognormal size distribution:

$$\begin{array}{}\text{(1)}& n\left(\ln D\right)={\displaystyle \frac{N}{\sqrt{\mathrm{2}\mathit{\pi }}\ln {\mathit{\sigma }}_{\mathrm{g}}}}\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\ln D-\ln D_{\mathrm{g}}}{\ln {\mathit{\sigma }}_{\mathrm{g}}}}\right)}^{\mathrm{2}}\right],\end{array}$$

where D is the particle diameter, N is the number distribution of all particles in the distribution, D_g is the geometric mean diameter, and σ_g is the geometric standard deviation. To calculate the aerosol optical properties, the aerosol spectrum is first divided into 48 bins from 0.002 to 2.5 µm, with radius r_i. The aerosols are classified into four types: (1) internally mixed sulfate, nitrate, ammonium, hydrophilic organics and black carbon, and water; (2) hydrophobic organics; (3) hydrophobic black carbon; and (4) other unidentified aerosols. These four kinds of aerosols are assumed to be mixed externally. For the internally mixed aerosols, the complex refractive index at a certain wavelength (λ) is calculated based on the volume-weighted average of the individual refractive index. Given the particle size and complex refractive index, the extinction efficiency (Q_e), ω_a, and g_a are calculated using the Mie theory at a certain wavelength (λ). The lookup tables of Q_e, ω_a, and g_a are established according to particle sizes and refractive indices to avoid multiple Mie scattering calculation. The aerosol optical parameters are interpolated linearly from the lookup tables with the calculated refractive index and particle size in the module.

The aerosol optical depth (AOD or τ_a) at a certain wavelength (λ) in a given atmospheric layer k is determined by the summation over all types of aerosols and all bins:

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{\mathrm{a}}\left(\mathit{\lambda }k\right)=\sum _{i=\mathrm{1}}^{\mathrm{48}}\sum _{j=\mathrm{1}}^{\mathrm{4}}{Q}_{\mathrm{e}}\left(\mathit{\lambda }{r}_{i}jk\right)\mathit{\pi }{r}_i^{\mathrm{2}}n\left(r_{i}jk\right)\mathrm{\Delta }{Z}_k,\end{array}$$

where $n\left(r_i,j,k\right)$ is the number concentration of jth kind of aerosols in the ith bin. ΔZ_k is the depth of an atmospheric layer. The weighted-mean values of σ and g are then calculated by (d'Almeida et al., 1991)

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\omega }}_a\left(\mathit{\lambda }k\right)={\displaystyle \frac{{\displaystyle \begin{array}{c}{\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_{j=\mathrm{1}}^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\\ {\mathit{\omega }}_{\mathrm{a}}(r_i,j,k)\mathrm{\Delta }{Z}_k\end{array}}}{\sum _{i=\mathrm{1}}^{\mathrm{48}}\sum _{j=\mathrm{1}}^{\mathrm{4}}{Q}_e(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\mathrm{\Delta }{Z}_k}},\text{(4)}& {\displaystyle }{g}_{\mathrm{a}}\left(\mathit{\lambda }k\right)={\displaystyle \frac{{\displaystyle \begin{array}{c}{\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_{j=\mathrm{1}}^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\\ {\mathit{\omega }}_{\mathrm{a}}(r_i,j,k)g_{\mathrm{a}}(\mathit{\lambda },r_i,j,k)\mathrm{\Delta }{Z}_k\end{array}}}{\begin{array}{c}{\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_{j=\mathrm{1}}^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\\ {\mathit{\omega }}_{\mathrm{a}}(r_i,j,k)\mathrm{\Delta }{Z}_k\end{array}}}.\end{array}$$

When the wavelength-dependent τ_a, ω_a, and g_a are calculated, they can be used in the Goddard shortwave module to evaluate the ARF. Detailed information can be found in Li et al. (2011b).

2.3 Data and statistical methods for comparisons

The model performance is validated using the available measurements in the NCP, including AOD, SSA, PBLH, downward shortwave flux (SWDOWN), aerosol species, and air pollutants. The daily AOD is retrieved from Terra and Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) level 2 products, with a resolution of 0.1^∘ × 0.1^∘. The hourly SSA is calculated using the measurement of the turbidity meter at the National Center for Nanoscience and Technology (NCNST), Chinese Academy of Sciences (116.33^∘ E, 39.99^∘ N), in Beijing (Fig. 1b). The daily PBLH at 12:00 Beijing time (BJT) is diagnosed from the radiosonde observation at a meteorological site (116.47^∘ E, 39.81^∘ N) in Beijing. The SWDOWN is measured by CM-11 pyranometers at four sites from the Chinese Ecosystem Research Network (CERN) in the NCP (Liu et al., 2016). The hourly measurements of O₃, NO₂, SO₂, CO and PM_2.5 concentrations have been released by China's Ministry of Ecology and Environment (China MEP) since 2013. The hourly submicron sulfate, nitrate, ammonium, and organic aerosols are measured by the Aerodyne Aerosol Chemical Speciation Monitor (ACSM) at NCNST. The primary organic aerosol (POA) and SOA concentrations are obtained from the ACSM measurement analyzed using the positive matrix factorization (PMF). In addition, we have also analyzed the relationship between near-surface [PM_2.5] and the PBLH retrieved from the lidar measurement at the Institute of Remote Sensing and Digital Earth (IRSDE), Chinese Academy of Sciences (116.38^∘ E, 40.00^∘ N), in Beijing (Fig. 1b).

In the present study, the mean bias (MB), root mean square error (RMSE), and the index of agreement (IOA) are used to assess the performance of WRF-Chem model simulations against measurements. The detailed description can be found in the Supplement.

3 Results and discussions

3.1 Model performance

We first define the base simulation in which ARF is considered (hereafter referred to as f_base), and results from f_base are compared to observations in the NCP. Generally, the model simulates the horizontal distributions and temporal variations of PM_2.5, O₃, NO₂, and SO₂ mass concentrations against measurements in the NCP well. Additionally, the model also reproduces the temporal profiles of the aerosol species reasonably well compared to observations in Beijing. The detailed model validation of air pollutants in the NCP and the aerosol species in Beijing can be found in the Supplement.

Figure 2 (a) Scatter plot of the MODIS-retrieved and simulated daily AOD; (b) Taylor diagram (Taylor, 2001) to present the variance, bias, and correlation of the retrieved and simulated daily AOD averaged in the NCP from 5 December 2015 to 4 January 2016.

3.1.1 Aerosol radiative property simulations in the NCP

Aerosol radiative forcing mainly depends on the AOD, SSA, and asymmetry parameter (g). The model validations of AOD and SSA are provided in this study to further evaluate the aerosol radiative effect on the air pollution. The daily AOD at 550 nm, retrieved from Terra and Aqua MODIS level 2 products, is compared with the simulation. Figure 2a shows the scatter plot of the daily retrieved and simulated AOD averaged in the NCP from 5 December 2015 to 4 January 2016. The simulated daily average AOD correlates well with the observation, with a correlation coefficient of 0.86. Generally, the retrieved and simulated AOD increases with deterioration of the haze pollution, but the model considerably underestimates the AOD against the observation. Figure 2b presents the Taylor diagram (Taylor, 2001) to show the variance, bias, and correlation of the simulated and retrieved AOD from 5 December 2015 to 4 January 2016. There exists a good relationship between the simulated and retrieved daily AOD during the study episode, with correlation coefficients generally ranging from 0.5 to 0.9 and standard deviation mostly varying from 0.25 to 1.0. Figure 3 shows the pattern comparison of the retrieved and simulated AOD averaged during the simulation period. The model reasonably reproduces the AOD distribution compared to the observations in the NCP but considerably underestimates the AOD. The simulated and retrieved AOD averaged in the NCP during the simulation period is 0.43 and 0.59, respectively. It is worth noting that the simulated AOD is not only dependent on the column aerosol content and constituent, but is also significantly influenced by the relative humidity (RH) controlling the aerosol hydroscopic growth. Additionally, the satellite-retrieved AOD is subject to contamination by the existence of clouds, and considering the high occurrence frequency of clouds during haze days, the retrieved AOD is generally higher than the simulation (Engstrom and Ekman, 2010; Chand et al., 2012; Grandey et al., 2013).

Figure 3 Spatial distribution of (a) retrieved and (b) simulated AOD averaged from 5 December 2015 to 4 January 2016 in the NCP.

Aerosols are the mixture of absorbing and scattering constituents in the atmosphere. Their radiative effect of cooling or warming the atmosphere relies on many parameters, and SSA is one of the most important ones (Satheesh et al., 2010). Figure 4 depicts the comparison of the measured and simulated diurnal profiles of SSA at NCNST in Beijing during the episodes. The model performs reasonably in simulating the daily variation of SSA in Beijing, with an IOA of 0.69 and a MB of 0.0, but the overestimation or underestimation is rather large. SSA is the ratio of scattering to extinction, which is highly sensitive to the relative distribution of scattering and absorbing aerosol constituents in the atmosphere, as well as the RH determining the hygroscopic growth of aerosols. Therefore, the uncertainties of the simulated SSA probably originated from the model biases of aerosol constituents and the RH.

Figure 4 Comparison of measured (black dots) and predicted (red line) diurnal profiles of SSA in Beijing from 5 December 2015 to 4 January 2016.

3.1.2 Downward solar radiation simulations in the North China Plain

Figure 5 presents the daily profiles of simulated and observed SWDOWN at ground surfaces in Beijing, Jiaozhouwan, Luancheng, and Yuancheng from 5 December 2015 to 4 January 2016. The WRF-Chem model simulates the daily variation of SWDOWN well, especially in Jiaozhouwan, Luancheng, and Yucheng, with IOAs around 0.90. The model is subject to overestimating the SWDOWN against measurements, with MBs ranging from 6.3 to 86.2 W m⁻². The SWDOWN reaching the ground surface is very sensitive to the cloud cover and optical thickness. However, the WRF-Chem model still has difficulties in accurately predicting the cloud cover and optical thickness, which might constitute one of the most important reasons for model biases of the SWDOWN. In addition, the horizontal resolution used in simulations cannot adequately resolve the cumulus clouds, also causing uncertainties in the simulations of the SWDOWN.

Figure 5 Comparison of measured (black dots) and predicted (red line) diurnal profiles of the SWDOWN reaching the ground surface in (a) Beijing, (b) Jiaozhouwan, (c) Luancheng, and (d) Yucheng from 5 December 2015 to 4 January 2016.

3.1.3 PBLH simulations in Beijing

Figure 6 shows the temporal variations of the observed and simulated PBLH at a meteorological site in Beijing from 5 December 2015 to 4 January 2016. The average PBLH at 12:00 BJT during the episode at the meteorological site is 465.2 m, with the minimum of 101.8 m and the maximum of 1017.9 m, showing decreased PBLH during the haze episode. In general, the WRF-Chem model tracks the daily variation of the PBLH in Beijing reasonably, with an IOA of 0.70. However, the model has difficulties in reproducing the observed very low PBLH, e.g., less than 200 m. The PBLH varies substantially with time due to many factors, including large-scale dynamics, cloudiness, convective mixing, and the diurnal cycle of solar radiation (Sivaraman et al., 2013). Therefore, the simulation uncertainties of meteorological conditions constitute the main reason for the simulation bias of PBLH. For example, the overestimation of SWDOWN at 12:00 BJT (Fig. 5a) probably caused the overestimation of PBLH in Beijing.

Figure 6 Comparison of predicted diurnal profile (red line) of PBLH from 5 December 2015 to 4 January 2016 with observations at 12:00 BJT in Beijing.

In general, the simulated variations of SWDOWN, PBLH, aerosol radiative properties, air pollutants (PM_2.5, O₃, NO₂, SO₂, CO), and aerosol species are in good agreement with observations, indicating that the simulations of meteorological conditions, chemical processes, and the emission inventory used in the WRF-Chem model are reasonable, providing a reliable basis for the further investigation.

3.2 Relationship between near-surface [PM_2.5] and PBLH

Figure 7 presents the scatter plot of the lidar-retrieved PBLH at IRSDE and near-surface [PM_2.5] at a monitoring site close to IRSDE during daytime (08:00–17:00 LT) from 8 January to 20 February 2014. The wind speeds (WSPD) at a meteorological site close to IRSDE are shown by the color of the filled circles in Fig. 7. Additionally, near-surface [PM_2.5] during daytime are also subdivided into 20 bins with the interval of 25 µg m⁻³. The PBLH as the bin of near-surface [PM_2.5] is assembled, and an average of PBLH in each bin is calculated (Nakajima et al., 2001; Kawamoto et al., 2006), which is represented by the rectangle in Fig. 7. Generally, on average, when the PBLH decreases from 1500 to around 400 m, the near-surface [PM_2.5] increase from 10 to more than 200 µg m⁻³. When near-surface [PM_2.5] exceed 200 µg m⁻³, the PBLH remains 400–500 m. Previous studies have also reported the nonlinear relationship between the PBLH and near-surface [PM_2.5] and proposed that increasing [PM_2.5] reduce the PBLH or the ARF is attributed to the PBLH decrease (e.g., Petäjä et al., 2016; Tie et al., 2017; Liu et al., 2018).

Figure 7 Scatter plot of the PBLH and near-surface [PM_2.5] at the IRSDE site from 12 January to 20 February 2014. The black rectangle shows the bin average of PBLH. The color of the filled circles denotes the WSPD at the meteorological site close to IRSDE in Fig. 1b.

The PBLH is primarily determined by the wind shear in the vertical direction and the thermal condition of ground surfaces. The occurrence of low near-surface [PM_2.5] generally corresponds to efficient dispersions of PM_2.5 in horizontal and/or vertical directions. The strong horizontal winds in the lower atmosphere not only disperse PM_2.5 emitted or formed efficiently, but also intensify the wind shear in the vertical direction, increasing the PBLH and facilitating the rapid vertical exchange of PM_2.5 in the PBL. When near-surface [PM_2.5] are less than 50 µg m⁻³, the PBLH exceeding 1000 m is observed, which is chiefly determined by strong horizontal winds and less influenced by the ground thermal condition during wintertime, and the observed average WSPD is about 2.4 m s⁻¹. The occurrence of high near-surface [PM_2.5] indicates that the lower atmosphere is stable or stagnant, with weak horizontal winds and inactive convections, hindering the dispersion of PM_2.5 in the horizontal and vertical directions. Additionally, as the horizontal winds become weak or calm, the wind shear in the vertical direction is diminished and the PBLH is dominated by the ground thermal condition. When near-surface [PM_2.5] increase from 50 to around 200 µg m⁻³, the PBLH decreases from around 700 to 400 m, and the average WSPD decreases to 1.8 m s⁻¹. However, the increased PM_2.5 reducing PBLH still cannot be fully attributed to ARF, which is more likely caused by the decrease in winds or the formation of stagnant situations in the low-level atmosphere. When near-surface [PM_2.5] exceed 200 µg m⁻³, the observed PBLH fluctuates between 400 and 500 m with the average WSPD of around 1.0 m s⁻¹ and does not exhibit a continuous decrease with the increasing near-surface [PM_2.5].

Under the stagnant situation with weak winds, the PBLH is more sensitive to the ground thermal condition. Increasing aerosols or PM_2.5 in the low-level atmosphere attenuate the SWDOWN to the ground surface and decrease the surface temperature (TSFC) and turbulence kinetic energy, suppressing the PBL development and further enhancing near-surface [PM_2.5]. Therefore, with near-surface [PM_2.5] exceeding 200 µg m⁻³, the inert PBLH might be caused by the defect of the lidar-retrieved PBLH. The aerosol backscatter signal received by lidar is used to retrieve the PBLH. If the atmosphere is stable, the aerosols near the maximal PBLH are subject to being confined in situ, and the retrieved PBLH is generally the maximal one. Additionally, it is worth noting that the occurrence of the wintertime severe haze pollution in the NCP is often accompanied with the high-level convergence between 500 and 700 hPa, producing a persistent and strong sinking motion in the middle lower troposphere to reduce the PBLH and facilitate accumulation of air pollutants (Wu et al., 2017; Ding et al., 2017). Therefore, a subsidence inversion appears in the lower layer as a result of the air masses sinking in the middle troposphere, restraining the PBL development and determining the maximal PBLH. Hence, it is imperative to evaluate the contribution of ARF to the PBLH and near-surface [PM_2.5].

3.3 Sensitivity studies

The conceptual model about the ARF contribution to the heavy haze formation has been established in previous studies (e.g., Tie et al., 2017; Liu et al., 2018). During wintertime, under stagnant meteorological situations with weak winds and humid air, air pollutants are subject to accumulation in the PBL, facilitating the formation of PM_2.5. Increasing PM_2.5 in the PBL absorbs or scatters the incoming solar radiation to decrease the TSFC and facilitate anomalous temperature inversion, subsequently suppressing the vertical turbulent diffusion and decreasing the PBLH to further trap more air pollutants and water vapor to increase the RH in the PBL. Increasing RH enhances aerosol hygroscopic growth and multiphase reactions and augments the particle size and mass, causing further dimming and decrease in the TSFC and PBLH. The whole process constitutes a positive feedback induced by the aerosol radiation effect to enhance near-surface [PM_2.5], which has been proposed in many studies (Quan et al., 2013; Petäjä et al., 2016; Yang et al., 2016; Tie et al., 2017; Ding et al., 2017; Liu et al., 2018). The noted positive meteorological condition feedback has also been considered as the main reason for the near-surface PM_2.5 explosive growth (Zhong et al., 2018; X. Zhang et al., 2018b).

To comprehensively evaluate the influence of ARF on near-surface [PM_2.5] during the haze episode, a sensitivity study in which ARF is turned off has been conducted (hereafter referred as f_rad0). Therefore, the contribution of ARF to near-surface [PM_2.5] can be determined by the difference between f_base and f_rad0 (f_base−f_rad0). The most polluted area in the NCP is first selected to verify the conceptual model of the ARF contribution to the heavy haze formation, with the average near-surface [PM_2.5] during the haze episode exceeding 150 µg m⁻³. Figure 8 provides the temporal variation of near-surface [PM_2.5], SWDOWN, TSFC, PBLH, and RH averaged in the selected area during the episode in f_base and f_rad0. Apparently, ARF considerably decreases the solar radiation reaching the ground surface and correspondingly lowers the TSFC (Fig. 8b and c). Subsequently, the PBLH is decreased and the surface RH is increased due to decreasing TSFC during daytime (Fig. 8d and e). However, the variation trend of near-surface [PM_2.5], PBLH, TSFC, and RH due to ARF is not similar to that proposed in the conceptual model. During the haze development stage, whether ARF is considered or not, the TSFC and RH exhibit an increasing trend, showing the air mass originated from the south, and the PBLH does not consistently decrease with increasing near-surface [PM_2.5]. Additionally, the ARF contribution to near-surface [PM_2.5] is generally marginal during the haze development stage. During the haze maturation stage, ARF commences to elevate near-surface [PM_2.5] appreciably. It is worth noting that, even if ARF is not considered in f_rad0, the heavy haze pollution still occurs during the episode. For example, from 17 to 20 December 2015, without ARF, near-surface [PM_2.5] still continue to increase from around 30 to 300 µg m⁻³ and fluctuate between 150 and 300 µg m⁻³ until the occurrence of favorable meteorological conditions on 25 December. Hence, according to the variation trend of near-surface [PM_2.5] with and without the ARF contribution, the continuous accumulation of PM_2.5 during the haze episode is not primarily caused by ARF but predominantly induced by the stagnant meteorological conditions as well as the massive air pollutant emissions in the NCP. Figure 9 presents the temporal variation of AOD at 550 nm averaged in the selected area during the episode in f_base and f_rad0 to evaluate the impact of ARF on AOD. Apparently, except from 8 to 11 December, the ARF contribution to AOD is generally marginal, indicating that ARF does not play an important role in the column-integrated aerosol abundance. Additionally, the considerable AOD enhancement from 8 to 11 December is more likely caused by the substantial increase in RH due to ARF, which facilitates aerosol hygroscopic growth to augment particle size and further increases AOD. It is worth noting that the extinction of haze aerosols in the PBL also decreases the photolysis to suppress the photochemistry, further hindering the secondary aerosol formation to offset effects of ARF on near-surface [PM_2.5].

Figure 8 Temporal variations of the average (a) near-surface [PM_2.5], (b) SWDOWN at the ground surface, (c) TSFC, (d) PBLH, and (e) RH in the most polluted area in the NCP with [PM_2.5] of more than 150 µg m⁻³ in f_base (red solid line) and f_rad0 (blue solid line) from 5 December 2015 to 4 January 2016.

Figure 9 Temporal variations of the average AOD at 550 nm in the most polluted area in the NCP with [PM_2.5] of more than 150 µg m⁻³ in f_base (red solid line) and f_rad0 (blue solid line) from 5 December 2015 to 4 January 2016.

Figure 10 Average (a) percentage decrease in SWDOWN at the ground surface, (b) decrease in TSFC, (c) decrease in WSPD, (d) percentage decrease in PBLH, (e) increase of RH, and (f) percentage contribution of near-surface [PM_2.5] caused by ARF, as a function of the near-surface [PM_2.5] in the NCP during daytime from 5 December 2015 to 4 January 2016.

In order to quantitatively evaluate effects of ARF on near-surface [PM_2.5], which cannot be reflected by the temporal variation of near-surface [PM_2.5], TSFC, PBLH, and RH, an ensemble method is used in this study. The daytime near-surface [PM_2.5] in the NCP during the episode in f_base are first subdivided into 30 bins with an interval of 20 µg m⁻³. The SWDOWN, TSFC, PBLH, the near-surface WSPD, RH, and [PM_2.5] in f_base and f_rad0 in the same grid cell are assembled as the bin [PM_2.5], respectively, and an average of these variables in each bin is calculated. Figure 10 shows the decrease in SWDOWN (%), TSFC (^∘C), PBLH (%), WSPD (m s⁻²) and the increase of RH (%, not percentage change) and near-surface [PM_2.5] contribution (%) caused by ARF as a function of bin [PM_2.5]. The SWDOWN reaching the ground surface almost decreases linearly with the enhancement of near-surface [PM_2.5]. When ARF is considered, aerosols in the atmosphere absorb or scatter the incoming solar radiation, directly attenuating the radiation reaching the ground surface. When near-surface [PM_2.5] exceed 200 µg m⁻³, the SWDOWN at ground surfaces decreases by more than 20 % (Fig. 10a). Moreover, the decrease in the SWDOWN correspondingly lowers the TSFC, and the decrease in the TSFC is generally proportional to near-surface [PM_2.5], about 0.35 ^∘C per 100 µg m⁻³ PM_2.5 (Fig. 10b). Interestingly, ARF also decreases near-surface WSPD by about 0.1–0.2 m s⁻¹ with near-surface [PM_2.5] exceeding 80 µg m⁻³ (Fig. 10c). When severe air pollution occurs in the NCP during wintertime, atmospheric convergence occurs in the PBL (Liao et al., 2015; Ding et al., 2017). However, the ARF-induced cooling in the low-level air generates a divergence in the NCP, causing the decrease in near-surface WSPD.

Figure 11 Average decrease in (a) near-surface water vapor content and (c) vertical velocity below 400 m caused by ARF, and (b) average vertical velocity below 400 m in f_rad0 as a function of the near-surface [PM_2.5] in the NCP during daytime from 5 December 2015 to 4 January 2016.

The PBLH is primarily determined by the atmospheric dynamic and thermal condition of ground surfaces. Therefore, the decrease in WSPD and TSFC due to ARF subsequently suppresses the PBL development and diminishes the PBLH (Fig. 10d). When near-surface [PM_2.5] are less than 250 µg m⁻³, the PBLH decreases rapidly with increasing [PM_2.5]. When the near-surface [PM_2.5] are between 250 and 350 µg m⁻³, the decrease in PBLH is around 28 %. With near-surface [PM_2.5] more than 350 µg m⁻³, the decrease in PBLH exceeds 30 %. As for the ARF effect on water vapor in the PBL, the conceptual model has proposed that the decreased PBL induced by ARF weakens the vertical exchange of water vapor or the dispersion of water vapor is constrained by the shallow PBL (Tie et al., 2017; Liu et al., 2018). However, Fig. 11a shows that ARF decreases the near-surface water vapor content slightly, by more than 0.1 g kg⁻¹ with near-surface [PM_2.5] exceeding 100 µg m⁻³. During the haze episode in the NCP, the abundant moisture in the PBL is mainly transported from the south. The divergence due to cooling caused by ARF weakens the prevailing southerly wind and decreases the moisture transport from the south, reducing the water vapor content in the NCP. Considering that the RH is sensitive to the temperature with a constant water vapor content, the ARF-induced cooling still increases the near-surface RH (Fig. 10e). When near-surface [PM_2.5] exceed 300 µg m⁻³, the RH is increased by more than 5 %, so the heavy haze generally causes the air to be more humid.

Figure 12 Near-surface [PM_2.5] contribution caused by ARF, averaged from 5 December 2015 to 4 January 2016 in the NCP.

More PM_2.5 emitted or formed is trapped by a shallow PBL caused by ARF, and increased RH promotes the aerosol hygroscopic growth and further multiphase reactions, progressively enhancing near-surface [PM_2.5] (Fig. 10f). When near-surface [PM_2.5] are more than 50 µg m⁻³, the contribution of ARF to near-surface [PM_2.5] consistently increases with the haze deterioration. When the severe haze occurs, i.e., near-surface [PM_2.5] exceed 250 µg m⁻³, more than 12 % or 30 µg m⁻³ PM_2.5 is contributed by ARF. The simulated ARF effects on near-surface [PM_2.5] are generally comparable to those reported by previous studies. Z. F. Wang et al. (2014) have shown that ARF increases the monthly PM_2.5 concentration by 10 %–30 % in Beijing–Tianjin–Hebei in January 2013. Using the WRF-Chem model, Gao et al. (2015) have indicated that ARF increases the PM_2.5 concentration by 10–50 µg m⁻³ (2 %–30 %) over Beijing, Tianjin, and south Hebei from 10 to 15 January 2013, a period with the simulated maximum hourly surface PM_2.5 concentration of more than 600 µg m⁻³. X. Zhang et al. (2018a) have also quantified the aerosol–meteorology interaction effect on PM_2.5 concentrations in China in 2014 using the WRF-Chem model, showing that the increase of PM_2.5 concentrations associated with ARF is up to 16 % in China. Other previous studies have also confirmed the ARF effect during the heavy haze pollution episode (Wang et al., 2015; Zhang et al., 2015; Gao et al., 2016). However, when near-surface [PM_2.5] are less than 50 µg m⁻³, the contribution of ARF to near-surface [PM_2.5] is negative, although ARF decreases PBLH and increases RH. One of the possible reasons for the negative contribution of ARF is perturbations of wind fields caused by the ARF-induced cooling. Figure 11b presents the average vertical velocity (the net velocity by combining updrafts and downdrafts) below about 400 m in f_rad0 as a function of near-surface [PM_2.5]. Apparently, when ARF is not considered, the area with near-surface [PM_2.5] less than 100 µg m⁻³ is generally controlled by downward airflow, and vice versa, for the area with near-surface [PM_2.5] more than 100 µg m⁻³. The ARF-induced cooling generally causes a downward motion in the PBL (Fig. 11c), which suppresses the upward motion in the area with near-surface [PM_2.5] more than 100 µg m⁻³ to enhance near-surface [PM_2.5], but accelerates the downward motion in the area with near-surface [PM_2.5] less than 100 µg m⁻³ to strengthen the divergence intensity, further decreasing near-surface [PM_2.5]. Countered by the decrease in PBLH and increase of RH, the ARF contribution becomes positive with near-surface [PM_2.5] exceeding 50 µg m⁻³.

Figure 13 TSFC and wind filed variations caused by ARF, averaged from 5 December 2015 to 4 January 2016 in the NCP.

Figure 12 presents spatial distributions of the average near-surface PM_2.5 contribution due to ARF during the episode. The average near-surface PM_2.5 contribution caused by ARF in the NCP is 10.2 µg m⁻³ or 7.8 %, with the maximum exceeding 40 µg m⁻³ in the south of Hebei. On average, the ARF contribution to near-surface [PM_2.5] is the most significant in Tianjin, about 17.6 µg m⁻³ or 10.3 %, followed by Hebei (11.6 µg m⁻³ or 9.3 %), Shandong (11.5 µg m⁻³ or 7.3 %), Henan (11.2 µg m⁻³ or 7.7 %), Anhui (7.7 µg m⁻³ or 7.4 %), Beijing (7.3 µg m⁻³ or 6.9 %), and Jiangsu (7.0 µg m⁻³ or 6.2 %). It is noteworthy that the ARF contribution during the episode in North China is generally positive, but in its surrounding area the contribution becomes negative. At a large scale, when the air pollution occurs during wintertime in North China, the vertical motion over the polluted area generally shows an ascending–descending–ascending distribution from the surface to the middle level of the troposphere, and wind directions present a structure of convergence–divergence–convergence accordingly (Liao et al., 2015; Wu et al., 2017; Ding et al., 2017). ARF cools the low-level atmosphere and induces a downward motion, which suppresses the upward motion in the convergence area in North China to increase near-surface [PM_2.5], but accelerates the downward motion in the divergence area to decrease [PM_2.5].

Furthermore, when ARF is considered, near-surface [PM_2.5] over the East and South China Sea are also increased, with an enhancement less than 5 µg m⁻³ (about 3 % to more than 15 %). Considering the low near-surface [PM_2.5] over sea, the [PM_2.5] enhancement might be caused by the PM_2.5 transport from the continent. Figure 13 shows the spatial distribution of the TSFC and wind field variation caused by ARF averaged during the episode. Apparently, ARF causes a widespread cooling effect in East China, and the cooling is the most significant in the NCP, with the maximum TSFC decrease exceeding 1.5 ^∘C. The cooling effect in the NCP induces a weak northerly wind, decreasing the prevailing southerly wind during the haze episode (Fig. 13). Additionally, the cooling effect over the continent also intensifies the temperature contrast between land and sea, producing a secondary circulation to transport the PM_2.5 from the continent to the East and South China Sea.

4 Conclusions

In the study, a persistent haze pollution episode in the NCP from 5 December 2015 to 4 January 2016 is simulated using the WRF-Chem model to verify the ARF contribution to the haze formation. Generally, the model reproduces the spatial distributions and temporal variations of PM_2.5, O₃, NO₂, SO₂, and CO mass concentrations against observations in the NCP well. The calculated temporal variations of aerosol species are also consistent with the ACSM measurement in Beijing, particularly with regard to the simulation of sulfate, nitrate, and ammonium. Moreover, the model simulates the variation of SWDOWN, PBLH, and aerosol radiative properties during the episode reasonably well, compared to the measurement.

Previous studies have established that a positive feedback induced by ARF causes the heavy haze formation by modulating the PBL and RH. However, model results demonstrate that, during the haze development stage in the NCP, ARF does not dominate the accumulation of near-surface [PM_2.5], while ARF considerably enhances near-surface [PM_2.5] during the haze mature stage.

Ensemble analyses of model results show that, during daytime, ARF attenuates SWDOWN reaching ground surfaces efficiently, and correspondingly the TSFC progressively decreases with increasing near-surface [PM_2.5] in the NCP. The ARF-induced cooling generates a divergence in the low-level atmosphere in the NCP, lowering the near-surface WSPD and decreasing the water vapor transport from the south. The decreased WSPD and TSFC caused by ARF hinder the PBL development and the PBLH decreases rapidly with increasing near-surface [PM_2.5]. Although the water content in the NCP is decreased slightly, the RH is still increased due to the ARF-induced cooling. A shallow PBL and more humid air caused by ARF accelerate the PM_2.5 accumulation and secondary pollutant formation, facilitating heavy haze formation. The contribution of ARF to near-surface [PM_2.5] increases from 12 % to 20 % when near-surface [PM_2.5] increase from 250 to 500 µg m⁻³. However, ARF decreases the PM_2.5 level with near-surface [PM_2.5] less than 50 µg m⁻³.

The average near-surface PM_2.5 contribution of ARF during the episode in the NCP is 10.2 µg m⁻³ or 7.8 %. ARF aggravates the heavy haze formation in North China, but in its surrounding area ARF slightly mitigates the haze pollution. Generally, there is a structure of convergence–divergence–convergence over the polluted area of North China from the surface to the middle level of the troposphere. A downward motion is induced due to the widespread cooling effect of the low-level atmosphere caused by ARF, impeding the upward motion in the convergence area in North China to increase near-surface [PM_2.5] but accelerating the downward motion in the divergence area to decrease [PM_2.5].

Although the model performs generally well in simulating air pollutants, aerosol species and radiative properties, SWDOWN, and PBLH, the uncertainties from meteorological fields and the emission inventory still have the potential to influence ARF evaluation. In particular, further studies need to be conducted to improve the AOD simulations. In this study, ARF only considers the aerosol effect on the solar radiation, and the influence of longwave radiation also needs to be included. It is worth noting that modification of photolysis by aerosol scattering or absorbing solar radiation ultimately alters the atmospheric oxidizing capacity to influence the secondary aerosol formation, which potentially offsets the ARF effect on the haze pollution. Hence, further studies need to be performed to evaluate the effect of aerosol photolysis interaction on the haze pollution. In addition, aerosols play an important role in the cloud process serving as cloud condensation nuclei (CCN) and ice nuclei (IN). Therefore, aerosol–cloud interactions (aerosol indirect effect) modify temperature and moisture profiles and further influence precipitation, leading to potential effects on the atmospheric chemistry (Wang et al., 2011). Future studies should be performed to investigate the feedbacks of the aerosol indirect effect on the air pollutants.