Lidar vertical observation network and data assimilation reveal key processes driving the 3-D dynamic evolution of PM_2.5 concentrations over the North China Plain

China has made great efforts to monitor and control air pollution in the past decade. Comprehensive characterization and understanding of pollutants in three-dimensions are, however, still lacking. Here, we used data from an observation network consisting of 13 aerosol lidars and more than 1000 ground observation stations combined with a data assimilation technique to conduct a comprehensive analysis of extreme heavy aerosol pollution (HAP) over the North China Plain (NCP) from November–December 2017. During the studied period, the maximum hourly mass concentration of surface PM_2.5 reached ∼390 µg m⁻³. After assimilation, the correlation between model results and the independent observation sub-dataset was ∼50 % higher than that without the assimilation, and the root mean square error was reduced by ∼40 %. From pollution development to dissipation, we divided the HAP in the NCP (especially in Beijing) into four phases: an early phase (EP), a transport phase (TP), an accumulation phase (AP), and a removal phase (RP). We then analyzed the evolutionary characteristics of PM_2.5 concentration during different phases on the surface and in 3-D space. We found that the particles were mainly transported from south to north at a height of 1–2 km (during EP and RP) and near the surface (during TP and AP). The amounts of PM_2.5 advected into Beijing with the maximum transport flux intensity (TFI) were through the pathways in the relative order of the southwest > southeast > east pathways. The dissipation of PM_2.5 in the RP stage (with negative TFI) was mainly from north to south with an average transport height of ∼1 km above the surface. Our results quantified the multi-dimensional distribution and evolution of PM_2.5 concentration over the NCP, which may help policymakers develop efficient air pollution control strategies.

1 Introduction

Frequent heavy air pollution has exerted significant impacts on air visibility, climate, human health, and other environmental concerns (J. Li et al., 2017; Gao et al., 2017; Pokharel et al., 2019; Su et al., 2020). As a developing country with the largest population in the world, China's air quality has exhibited an obvious improvement trend in recent years (Sun et al., 2015; Cao et al., 2017). Regional air pollution in China is still serious, however, especially the heavy aerosol pollution (HAP) caused by fine particulate matter (PM_2.5) in winter, which has attracted attention worldwide (Zheng et al., 2015, 2019; X. Li et al., 2017; Cheng et al., 2016). Therefore, providing a reliable distribution of the PM_2.5 concentration of HAP, especially at any time and at any height in a given region, is particularly important in the quest for avoiding public health problems and providing government policymakers with help in designing effective controls (Hu et al., 2015).

Compared with other air pollutants (e.g., ozone and nitrogen dioxide), PM_2.5 has a longer atmospheric lifetime (3–5 d), during which it can be transported vertically to great heights and horizontally hundreds of kilometers (Wang et al., 2017; Zhang et al., 2014) depending on the meteorological conditions (e.g., relative humidity and precipitation) and chemical composition (Yang et al., 2017). Previous studies demonstrated that regional transport plays an important role for pollution formation in the major cities of China, e.g., transport contributes over 50 % of the PM_2.5 mass concentration in Beijing city, Shanghai city, Hangzhou city, Guangzhou city, Hong Kong, and Chengdu city during the relatively polluted period (Sun et al., 2017). Annually from 2005–2010, about 35.5 % (32.8 µg m⁻³) of the PM_2.5 in Beijing was attributed to regional transport from the North China Plain (NCP), within which up to 60.4 % (64.3 µg m⁻³) was from southerly and westerly air flows (Wang et al., 2015). Since the 2013 implementation of the most stringent clean air policy in China, the control of local pollution sources has led to the rapid reduction of total PM_2.5 concentration (J. Wang et al., 2019). It should be noted, however, that the local contributions, intra-regional transport, and inter-regional transport accounted for 47 % (12.7 µg m⁻³), 25 % (6.6 µg m⁻³), and 28 % (7.6 µg m⁻³), respectively, of the total PM_2.5 for the Beijing-Tianjin-Hebei (BTH) region from 2014–2017, with the 2017 contribution of regional transport to the BTH concentration rate ranging from 32.5 %–68.4 % (Dong et al., 2020).

Previous studies have shown that it is difficult to use surface observations to characterize the impact of upper-level pollutants in the atmosphere (X. Huang et al., 2018), which is affected by local emissions, regional transport, meteorological conditions, geographical factors, etc. (Tao et al., 2020; Che et al., 2019a). Therefore, understanding the key processes that drive the dynamic temporal and spatial evolutionary characteristics of pollutants on the NCP is essential for revealing the source and transport of aerosols, which has different radiative forcing at different heights (Kumar et al., 2017; Che et al., 2019b). Actually, stereo-monitoring devices and technologies, such as lidar (Sheng et al., 2019; Fan et al., 2019; Z. Chen et al., 2019), MAX-DOAS (Hong et al., 2018; Zhang et al., 2020), and satellite remote sensing (Pang et al., 2018; Schwartz et al., 2012; C. X. Zhang et al., 2019) can reveal the vertical distribution of pollutants at different heights (Tian et al., 2017; Heese et al., 2017). Due to the limited spatial and temporal observations, however, it is impossible to provide physical and chemical properties in the atmosphere at any time period and on any path, which makes it difficult to directly reveal the formation and source of pollution.

On the other hand, although the distribution of pollutants can be simulated by air quality models (M. Huang et al., 2018; Zhang et al., 2008), large uncertainties remain, mainly from the influence of emission inventory, meteorological fields, and some hypothetical conditions (Xu et al., 2016; Zhang et al., 2017; Huang et al., 2016). Fortunately, the above observed data and the results of the model can be fused using data assimilation techniques, which can correct the model simulation results via the observed data (Wang et al., 2013; Ma et al., 2019). Research has shown that mainstream data assimilation (DA) technologies, including 3DVAR (Jiang et al., 2013; Ma et al., 2018), 4DVAR (Yumimoto et al., 2008), and EnKF (D. Chen et al., 2019), can be used to assimilate observation data from the surface, remote sensing data (such as aerosol optical depth – AOD) from satellites, and vertical profile data from lidar, all of which can be used to improve the performance of the model, including the simulation of PM_2.5 and PM₁₀.

In this study, we analyzed the observation data from a vertical observation network consisting of 13 lidars and surface observation stations during an extreme pollution event in eastern China, especially in the NCP. Next, all of the data were utilized by the Gridpoint Statistical Interpolation (GSI) three-dimensional (3-D) variational (3DVAR) data assimilation system to revise the PM_2.5 results from the WRF-Chem simulation (Pagowski et al., 2014). Finally, the multi-dimensional evolutionary characteristics of PM_2.5 at the surface and in the vertical layer, as well as the 3-D distribution, were analyzed in detail. Although data assimilation has been applied in China using surface observation network data (Gao et al., 2017), AOD (Saide et al., 2013, 2014; Schwartz et al., 2012; Liu et al., 2011), and lidar data (Cheng et al., 2019), to our knowledge, this is the first attempt in China to apply lidar network data to assimilation technology, from which the high-precision 3-D distribution of pollutants can be provided, thus supplying effective data support for clarifying the formation mechanism of pollutants (Zheng et al., 2017).

2 Measurements and methods

2.1 Lidar observation network

The vertical aerosol observation network of the NCP was composed of 13 aerosol lidar monitoring stations (Fig. 1) covering four main transport channels of Beijing pollutants, including the southwestern transport path of Baoding City (BD), Shijiazhuang City (SJZ), Xingtai City (XT), Handan City (HD), Xinxiang City (XX), and Yangquan City (YQ); the southern transport path of Dezhou City (DZ) and Jining City (JN); the southeastern transport path of Langfang City (LF), Cangzhou City (CZ), and Zibo City (ZB); the eastern transport path of Tianjin City (TJ); and a lidar in the urban area of Beijing (BJ).

Figure 1 © Google Earth map of (a) China with the studied cities and (b) the North China Plain with all the lidar stations. The data in brackets are the maximum PM_2.5 concentrations (µg m⁻³) at the surface during the observation period. The black arrows in (b) from left to right show that the main pollution pathways of Beijing come from the four directions of southwest, south, southeast, and east.

The lidar system was developed by the Anhui Institute of Optics and Fine Mechanics (AIOFM), Chinese Academy of Sciences (CAS), and was used for the long-term continuous observation of aerosol vertical distribution. The lidar system adopted the Nd:YAG laser, which emits a 532 nm wavelength with 30 mJ single-pulse energy and 10–30 Hz pulse repetition frequency. The vertical resolution is 7.5 m with the original time resolution of 3–10 min. The detection blind area is 0.1 km, and more specific technical details can be found in Xiang et al. (2019). The vertical distribution of the aerosol extinction coefficient was retrieved using the Fernald method (Fernald, 1984), which is more suitable for vertical detection and more accurate than the Collis (Collis et al., 1964) and Klett (Klett, 1981) methods (Schafer et al., 2015). Furthermore, combining the extinction coefficient with the PM_2.5 in situ surface observations, the vertical distribution of the PM_2.5 mass concentration in the boundary layer was obtained using the empirical formula fitting method, which has proven to be reliable and highly accurate; the specific technical details can be found in the literature (L. Lv et al., 2017; Tao et al., 2016; L. H. Lv et al., 2017). In addition, an image recognition algorithm was used to evaluate the height of the atmospheric boundary layer (Xiang et al., 2019; Barrera et al., 2019).

2.2 WRF-Chem model configurations

The WRF-Chem chemical transport model (version 3.8.1) was used to investigate the particulate concentrations and meteorological parameters in the study area and was configured with nested domains consisting of 100×100 (36 km) and 103×103 (12 km) grids (Fig. S1 in the Supplement). The domain had 41 vertical layers from the surface to 50 hPa. To better simulate the conditions within the boundary layer, the resolution of the boundary layer was increased and 20 layers were set in the range of 0–2 km. The initial and boundary meteorological conditions were derived from the 6 h National Centers for Environmental Prediction Final Analysis data with $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution. The inventory of anthropogenic emissions for 2016 was obtained from the Multi-resolution Emission Inventory for China (MEIC) data with $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ resolution (Zhou et al., 2017). Terrestrial biogenic emissions were estimated using the Model of Emissions of Gases and Aerosols from Nature (MEGAN) model (Chatani et al., 2011). The gas-phase chemistry module CBM-Z and the Model for Simulating Aerosol Interactions and Chemistry (MOSAIC) aerosol module were used in this simulation. Detailed information concerning the model configuration is provided in Table S1 in the Supplement. The model runs from 20 November–9 December 2017, and the results from 25 November–9 December 2017 were used for the analysis in Sect. 3.

Figure 2 PM_2.5 mass concentration comparison results from lidar at different heights (b–d, f–h) and surface observations (a, e) with non-assimilation simulations (a–d) and assimilation simulations (e–h).

2.3 GSI 3DVAR DA system

The GSI DA (Gridpoint Statistical Interpolation Data Assimilation) system provides 3DVAR analysis by minimizing the cost function as shown below (Gao et al., 2017):

$$\begin{array}{}\text{(1)}& J\left(\mathit{x}\right)={\left(\mathit{x}-{\mathit{x}}_{\mathrm{b}}\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{b}}\right)+(\mathit{y}-H(\mathit{x}\left){)}^T{\mathbf{R}}^{-\mathrm{1}}\right(\mathit{y}-H\left(\mathit{x}\right)).\end{array}$$

In this equation, x is the analysis vector, x_b denotes the background vector, y is an observation vector, B represents the background error covariance matrix, R represents the observation error covariance matrix, and H is the observation operator used to transform model grid point values to observed variables, which was performed via interpolation in our research. The background error covariance matrix was calculated using the National Meteorological Center (NMC) method (Parrish and Derber, 1992; Saide et al., 2013), which simulated the difference of results at the same time (25 November 2017) with two different starting times (20 and 21 November 2017, respectively). The 1 h assimilated window data included 13 groups (see Fig. 1 for site distribution) of PM_2.5 vertical profiles retrieved from lidar and the surface PM_2.5 data from hundreds of surface monitoring stations (see Fig. 5 for site distribution) from the China Environmental Monitoring Center. The observation errors of ground PM_2.5 and its vertical distribution (through the ground PM_2.5 fitting method in Sect. 2.1) originated from measurement errors and representative errors. The measurement error was computed using ${\mathit{\epsilon }}_{\mathrm{0}}=\mathrm{1.5}+\mathrm{0.0075}\cdot \mathrm{obs}$ (Pagowski et al., 2014), where “obs” indicates observed values. The representative error was computed using ${\mathit{\epsilon }}_{\mathrm{r}}=\mathit{\gamma }{\mathit{\epsilon }}_{\mathrm{0}}\sqrt{\mathrm{\Delta }x/L}$ (Elbern et al., 2007), where γ is the adjustable scale factor (we used the value of 0.5 recommended by the GSI system), Δx is the model grid resolution (we selected 12 km of domain 2), and L is the influencing radius (we used 60 km).

3 Results and discussion

3.1 Evaluation of assimilation performance using vertical PM_2.5 data

In order to evaluate the improvement of model simulation performance from data assimilation using lidar vertical profile data and surface station data, considering the sharp decline of PM_2.5 value at 1 km height (Fig. 6), only the non-assimilation and assimilation results at the surface, 0.2, 0.5, and 1 km were compared, as shown in Fig. 2. These data were selected from five of the most polluted stations, including the cities of TJ, LF, BD, SJZ, and XT. It should be noted that these observation data were not assimilated, which means that the following comparisons are independent (Bocquet et al., 2015). Obviously, the data assimilation used can greatly improve the simulation accuracy. Compared with the observation data at different heights (ground, 0.2, 0.5, and 1 km), the simulation results of PM_2.5 levels under the condition of non-assimilation were higher (Fig. 2a–d), the root mean square error (RMSE) was 52.14±20.27 µg m⁻³, and the correlation coefficient was only 0.56±0.15. Correspondingly, the results of PM_2.5 simulated with assimilation were closer to the observed values (Fig. 2e–h) and the RMSE was 33.07±14.69 µg m⁻³, which represents a reduction of about 40 % in simulation error after assimilation. The correlation coefficient was 0.81±0.10, demonstrating that the simulation accuracy was improved by about 50 % after assimilation.

In addition, compared with the simulation with assimilation (Fig. 5 in Sect. 3.3), the results without assimilation were significantly higher than the observed values (Fig. S2 in the Supplement), especially during the pollution period (Fig. S2d and e), which may be due to the simulation error caused by the model (Zhang et al., 2016). Meanwhile, the comparison of the 3-D results (Fig. 7 in Sect. 3.5 and Fig. S3 in the Supplement) further reveals that the simulation results of upper air PM_2.5 may also overestimate the actual values, which demonstrates the importance of data assimilation in capturing the 3-D structure of pollution.

3.2 The four phases from aerosol pollution development to dissipation

Joint observations and analyses have been widely performed in an effort to reveal the heavy aerosol pollution (HAP) in the NCP region (Li et al., 2016; Zhang et al., 2018). The key processes of a HAP event, from aerosol pollution development to dissipation, usually include an early phase (EP), a transport phase (TP), an accumulation phase (AP), and a removal phase (RP) (Yuan et al., 2019; Zhong et al., 2017), classifications that are based on the increase and decrease of PM_2.5 mass concentration in Beijing (BJ) caused by changes in meteorological conditions. Here, the curves in Fig. 3 show the temporal evolution of PM_2.5 mass concentration monitored at the surface in different cities on the NCP from 25 November–9 December 2017, while the superimposed colors represent the time-varying profiles of the simulated wind fields in BJ, Baoding (BD), Dezhou (DZ), and Langfang (LF), respectively. Overall, PM_2.5 with high concentrations was usually associated with pronounced southerly winds (S in Fig. 3) or southwesterly winds (SW in Fig. 3), while the PM_2.5 concentrations decreased significantly under the prevailing northerly winds (including the wind directions of N, NW, and NE in Fig. 3).

Figure 3 Surface PM_2.5 observations from different cities: (a) Beijing (including Tianjin) and its (b) southwest cities, (c) southeast cities, and (d) east cities for the period 25 November–9 December 2017. Superimposed colors represent the time-varying profiles of the simulated wind fields in Beijing, Baoding, Dezhou, and Langfang, respectively.

Furthermore, in order to characterize the evolution of PM_2.5 during different pollution phases, the period from 29 November–5 December was selected as a typical extreme HAP event covering the four pollution phases. This extreme pollution event lasted more than 4 d and featured a regional transport process. During the EP (29 November–noon 30 November; episode 1 in Fig. 3), the air quality in BJ and its surrounding areas such as Tianjin (TJ) was relatively good, with an average PM_2.5 value of ∼15 µg m⁻³, while slight pollution occurred to the southwest of BJ, including BD, Shijiazhuang (SJZ), Xingtai (XT), and Handan (HD), with an average value of ∼50 µg m⁻³.

During the TP (approximately the morning of 2 December; episode 2 in Fig. 3), the variation of PM_2.5 concentration was more sensitive and responded rapidly to the wind shift from northerly to southerly, causing the PM_2.5 concentration in Beijing to increase quickly from ∼30 to ∼50 µg m⁻³, while southwest of Beijing (e.g., BD, SJZ, XT, and HD) the PM_2.5 concentration increased rapidly to ∼200 µg m⁻³. Research has revealed that the pollutant transport south of Beijing, especially in the southwest areas (the Taihang Mountains), is the most important contribution source to Beijing pollutants (Zhao et al., 2020). During the AP (approximately 3 December; episode 3 in Fig. 3), diffusion of the pollutants was difficult due to the occurrence of a surface temperature inversion in Beijing (Fig. 4) (L. L. Wang et al., 2019), which caused the maximum concentration of PM_2.5 in Beijing to reach ∼250 µg m⁻³. Meanwhile, the PM_2.5 concentrations in TJ, LF, BD, and SJZ reached maximum values of ∼270, 250, 320, and 390 µg m⁻³, respectively. Conversely, the pollution levels in Shanghai (SH), Hefei (HF), and Wuhan (WH) in the southernmost section of the NCP were relatively low with average values < ∼60 µg m⁻³.

Figure 4 Time series of vertical temperatures (a, c, e, g) and temperature gradients (b, d, f, h) from Beijing (a, b), Baoding (c, d), Dezhou (e, f), and Langfang (g, h) simulated by the WRF-Chem model.

Figure 5 Spatial distribution of PM_2.5 in the surface layer during different phases after assimilation. The black arrows indicate the wind direction. The circles represent the in situ surface observations.

Figure 6 Time series of vertical distributions of the aerosol extinction coefficient observed on the North China Plain from 29 November–5 December 2017. The white dashed lines represent the approximate atmospheric boundary layer height. Missing datasets are plotted in white.

During the RP (approximately 5 December; episode 4 in Fig. 3), the wind direction shifted from southwest to north, transporting the relatively clean air in the north to the south and thereby causing the pollutant concentrations in Beijing to decrease rapidly. In just 9 h, the air quality improved from heavy pollution to excellent, and the PM_2.5 concentrations in the NCP also decreased significantly. Finally, by noon on 4 December, the pollutant concentrations in the NCP had reached a low level, with an average value of ∼40 µg m⁻³. In contrast, due to the continuous southward advection of pollutants, serious pollution occurred in SH, HF, and WH, where the PM_2.5 concentrations reached maximum values of ∼210, 310, and 280 µg m⁻³, respectively. These findings are also consistent with the results of previous studies on the regional transport of regional pollutants to the Yangtze River Delta (Hua et al., 2015), which showed them to be due to the continuous southward flow of northwest and northeast winds.

3.3 Spatial distribution of PM_2.5 concentration in the surface layer

Additionally, in order to analyze the pollution characteristics of the NCP, the spatial distribution results of PM_2.5 after data assimilation were plotted in Fig. 5 for all phases. The high concentrations of PM_2.5 in BJ were recorded during the TP, AP, and beginning of the RP, while the PM_2.5 concentrations at other times were lower. Moreover, during the EP, only the eastern cities of Shanxi (SX) Province experienced moderate pollution levels (Fig. 5a). During the TP, the pollutants in the south-central NCP were transported to the north of the NCP (Fig. 5b and c) as a result of the southwesterly wind field, and under the superposition of the local pollutant emissions from each city (J. Li et al., 2017), the cities on the windward side of the Taihang Mountains (e.g., HD, SJZ, and BD) quickly developed varying levels of heavy pollution. In addition, during the AP, due to the large-scale inversion (Fig. 4b, d, f, h) caused by the rapid temperature rise (Fig. 4a, c, e, g) of the NCP region at upper levels, the atmospheric stratification was stable, causing pollutant loading on the NCP (including BJ, BD, SJZ, HD, LF, CZ, and elsewhere) to increase (Fig. 5d), nearly reaching their pollution maxima (Fig. 3). Meanwhile, during the RP, affected by the cold air at upper levels (Fig. 4a, c, e, g) from the northwest and the shift in wind direction over the NCP from southwest to north, the pollution severity gradually eased from north to south (Fig. 5e), with the air quality in the northern part of the region improving significantly (Fig. 5f).

Figure 7 3-D distribution of PM_2.5 during different phases after assimilation. Colors within the boxes depict the PM_2.5 concentrations. The color-coded arrows represent the wind direction and speed at 1 km. On the tops of the boxes, the spatial distributions of wind speed at 10 m are plotted.

3.4 Vertical distribution of aerosols observed by the lidar network

In order to quantify the characteristic vertical distribution of aerosols, the observed aerosol extinction coefficients from the 13 lidar stations in the NCP were plotted, as shown in Fig. 6. These results revealed that on 29 November the aerosol concentration at the surface was relatively low, although pollutant transport at heights of 1–2 km (see Fig. 8a and e) occurred at six stations (BD, SJZ, YQ, XT, HD, and XX) on the windward side of the Taihang Mountains. Figure S4 in the Supplement demonstrates that these pollutants in the upper air come from the local emissions on the ground, which is due to the updraft lifting to 1–2 km above the ground on the night of 28 November. The upper air transport of pollutants continued until 1 December, at which point it merged with the surface flow. Contrary to this, the pollutant transport from north to south occurred at a height of 1 km during the RP (e.g., Fig. 6b and d–g). In addition, the atmospheric boundary layer height (ABLH) reached its highest value of the observation period from 29 to 30 November, averaging more than 1.5 km. The ABLH began to decrease on 1 December, averaging approximately 1 km on that day. The lowest value of the ABLH occurred on 2–3 December, when its average dropped to less than 0.5 km, making it difficult for pollutants to diffuse and causing heavy pollution in the NCP (Z. Q. Li et al., 2017). Fortunately, on 4 December, the atmospheric boundary layer gradually lifted, which was conducive to the diffusion of pollutants.

3.5 Dynamic 3-D evolution of the PM_2.5 concentrations

Figure 7 presents the 3-D distribution of PM_2.5 after assimilation, which clearly shows the generation, dissipation, transport, and diffusion characteristics of pollutants in the atmosphere. The tops of the boxes in the figure depict the wind speeds 10 m above the surface. During the EP, the high-concentration pollutants only occurred in the upper air within ∼1 km of the surface in SX Province (e.g., YQ). During the TP, the high-concentration pollutants were mainly found on the windward side of the Taihang Mountains (southwest pathway) and the loading height of PM_2.5 was <1 km, which is illustrated in Fig. 8. During the AP, the average concentration of pollutants > 200 µg m⁻³ mainly occurred near the surface. Meanwhile, the pollutants with low concentrations at upper levels could be transported to the Bohai Sea. During the RP, high-concentration pollutants > 100 µg m⁻³ simultaneously occurred over the Bohai Sea and the Yellow Sea.

Figure 8 Vertical profiles of PM_2.5 cross sections with wind vectors along the transport pathways, including southwest (a, e, i, m), south (b, f, j, n), southeast (c, g, k, o), and east (d, h, l, p). The first row (a–d) (00:00 LT 2 December 2017) represents the transport stage, the second row (e–h) (10:00 LT 2 December 2017) represents the accumulation stage, and the third row (i–l) (00:00 LT, 3 December 2017) and fourth row (m–p) (14:00 LT, 4 December 2017) represent the removal stage.

3.6 Quantification of regional transport of PM_2.5

To evaluate the variation of pollutants along different transport pathways at different stages, we plotted the vertical profile of the PM_2.5 cross section along the main pollution pathways of Beijing coming from the four directions of southwest, south, southeast, and east (see Fig. 1b and d). As shown in Fig. 8, at XX and XT (located at the start of the southwest transport pathway; Fig. 8a), the PM_2.5 concentration is more than 200 µg m⁻³ at a height of 1 km (Fig. 8a) and the surface PM_2.5 concentration at JN (located in the south pathway) also exceeds 200 µg m⁻³ (Fig. 8b). These high concentrations of pollutants were transported to SJZ, BD, LF, BJ, and other cities via southwest winds (Fig. 8e–g). At the same time, vertical downdrafts reduced the height of loading of aerosol layer to ∼0.6 km (Fig. 8e). Different from the southern (including southwest, south, and southeast) transport pathways, the pollutants in TJ were mainly from BJ outflow in all stages of the eastern transport pathways (Fig. 8d, h, l, p). In addition, wind direction inconsistencies at the origin (XX, JN, and ZB) and target location (Beijing) of the transport pathways occurred at the beginning of the removal phase (Fig. 8i–k), which may have been due to the southward delay of the northerly air flow.

To investigate the vertical variation of PM_2.5 inflow or outflow at different heights and determine the height at which the main transport occurred (H. Zhang et al., 2019), we plotted the vertical distribution of PM_2.5 transport flux in different directions (Fig. 9). Here the PM_2.5 transport flux is defined as the product of PM_2.5 mass concentration (µg m⁻³), wind speed (m s⁻¹), and wind direction projection in the current pathway (Xiang et al., 2020). The southwest, southeast, and east pathways in Fig. 9 were represented by BD, LF, and TJ, respectively, which are the three lidar stations closest to BJ (Fig. 1). TF > 0 indicates that the pollutants were imported to Beijing, while TF < 0 indicates that the pollutants were exported from Beijing. The results revealed that below the height of 3 km, the order of the maximum values of imported pollutants in the to Beijing direction was southwest pathway (1122.8 µg m⁻² s⁻¹) > southeast pathway (423.6 µg m⁻² s⁻¹) > east pathway (278.3 µg m⁻² s⁻¹), while the exported pollutants in the from Beijing direction was southwest pathway (−1571.4 µg m⁻² s⁻¹) > east pathway (−877.7 µg m⁻² s⁻¹) > southeast pathway (−772.4 µg m⁻² s⁻¹). Compared with the PM_2.5 transport flux on the ground surface, the relatively high value (∼200 µg m⁻² s⁻¹) in the southwest pathway (Fig. 9a) occurred on 29 November and early morning on 4 December, while the relatively extreme value ($\sim -\mathrm{400}$ µg m⁻² s⁻¹) on the east pathway (Fig. 9c) was recorded on the night of 2 December.

Figure 9 Time series of PM_2.5 transport flux from different transport pathways. The corresponding directions of the southwest, southeast, and east pathways are shown in Fig. 1.

To further obtain insights into the total transport characteristics in the target area (BJ) and its surrounding area (BD, LF, and TJ) during different evolutionary stages, the time series of the PM_2.5 transport flux intensity (TFI) was shown in Fig. 10. The TFI was calculated by integrating the PM_2.5 transport flux from the ground to a certain height, and the height was selected as 1.5 km, which is consistent with the main transport height of pollutants (Fig. 8) and the height of the boundary layer (Fig. 7). The TFI of PM_2.5 further reveals that pollutants imported into the Beijing area with a maximum PM_2.5 TFI of $\sim \mathrm{4.6}\times {\mathrm{10}}^{\mathrm{5}}$ µg m⁻¹ s⁻¹ were transported mainly via the southwest pathway during the TP, while the extreme TFI of pollutants exported from Beijing via the east pathway was approximately $-\mathrm{1.4}\times {\mathrm{10}}^{\mathrm{5}}$ µg m⁻¹ s⁻¹. In addition, during the RP, the pollutants from Beijing were exported to the southwest and southeast, with extreme values of approximately $-\mathrm{1.03}\times {\mathrm{10}}^{\mathrm{6}}$ and $-\mathrm{4.3}\times {\mathrm{10}}^{\mathrm{5}}$ µg  m⁻¹ s⁻¹, respectively. In contrast, the absolute value of TFI on the southwest pathway was $<\sim \mathrm{1.0}\times {\mathrm{10}}^{\mathrm{4}}$ µg m⁻¹ s⁻¹ during the EP (Fig. 10), which indicates that there was no significant inflow or outflow of pollutants. However, this reason was mainly due to the offsetting of the inflow of pollutants in the upper air and the outflow of pollutants near the ground (Fig. 9a). This special phenomenon also demonstrates that the study of vertical distribution of pollutants has great significance, which can better explain the transport characteristics (H. Zhang et al., 2019).

Figure 10 Time series of PM_2.5 transport flux intensity from different transport pathways. The corresponding directions of the southwest, southeast, and east pathways are shown in Fig. 1.

4 Conclusions

Accurate quantification of the distribution of particulate matter in the atmosphere is a key requirement for predicting air quality and estimating atmospheric environmental capacity from atmospheric observations. We utilized a vertical observation network composed of 13 aerosol lidars combined with data assimilation technology to improve the simulation accuracy of PM_2.5 and further analyzed the multi-dimensional evolutionary characteristics of pollutants in the surface layer, vertical layer, and 3-D space, thereby providing effective data support for clarifying the spatial transport characteristics of heavy pollution.

We found that the average height of the atmospheric boundary layer was <0.5 km during the HAP period. We further demonstrated that the transport of pollutants in the NCP region was mainly via three pathways: southwest, southeast, and east. During the TP, the PM_2.5 advected into Beijing with a maximum transport flux intensity (TFI) of $\sim \mathrm{4.6}\times {\mathrm{10}}^{\mathrm{5}}$ µg m⁻¹ s⁻¹ and was mainly via the southwest pathway, while the polluted air mass in the RP dissipated from Beijing via the southwest and southeast pathways, with extreme PM_2.5 TFI values of approximately $-\mathrm{1.03}\times {\mathrm{10}}^{\mathrm{6}}$ and $-\mathrm{4.3}\times {\mathrm{10}}^{\mathrm{5}}$ µg m⁻¹ s⁻¹, respectively. In addition, the transport of regional pollutants to the Yangtze River Delta was due to the continuous southward flow of northwest and northeast winds. Our results directly revealed that pollutants in the North China Plain can be transported to the Yellow Sea and the Bohai Sea, providing a dataset for a further in-depth study of the mechanism of air pollution in the coastal areas of eastern China. This study also captured the regional transport of air pollutants stretching over 1000 km, proving the necessity and importance of the joint prevention and control of regional air pollution.