Using Lidar technology to assess regional air pollution and improve estimates of PM_2.5 transport in the North China Plain

The overall structure of the Mie Lidar system consists of an emitting system, receiving system, data acquisition system, control system, and GPS system (Lv et al 2017b). The emitting system comprises a laser and an emission optical unit. The receiving system is designed based on the principle of the Cassegrain optical telescope, which has a 1-mard receiving field of view with a vertical resolution of 7.5 m. The data acquisition and control system comprise one of the key components used to determine the data quality of the Lidar. In addition, the response speed, bandwidth, dynamic range, and signal-to-noise ratio of the measured electronic system directly determine the spatial and temporal resolutions of the Lidar system and the analytical ability of weak signals. The main technical parameters of the Lidar system are listed in table S1 (available online at stacks.iop.org/ERL/15/094071/mmedia). Additional details regarding the ground-based Lidar information can be found in Xiang et al (2019).

In this study, we established two vehicle-mounted Lidars and five ground-based Lidar networks, covering the Beijing-Tianjin-Hebei region. The two vehicle-mounted Lidars performed closed observations along the 6th Ring Road of Beijing at 1:40 p.m. local time on 27 November 2017. One vehicle started north of the 6th Ring Road and traveled clockwise to the starting point, while the other vehicle started south of the 6th Ring Road. The observation route of the 6th Ring Road in Beijing is depicted by the red curve in figure 1. The 6th Ring Road in Beijing is a circular expressway with four lanes in two directions and a total length of approximately 187.6 km. The vehicle speed was set to 70 km h⁻¹ in order to maximize the detection precision, detection time, and other factors. The two vehicles returned to the starting position at 4:20 p.m. on the same day. Meanwhile, the blue stars in figure 1 represent the five ground-based Lidar networks, which were used for long-term continuous observations of the Beijing-Tianjin-Hebei region. The observation data from 26–29 November 2017, were selected to analyze the pollution status of the Beijing-Tianjin-Hebei region.

Zoom In Zoom Out Reset image size

2.3.1. Particle extinction coefficient

The Lidar equation for the retrieved aerosol Lidar signal can be written as

$$\begin{align}P\left( {z,\lambda } \right) & = \frac{{K{P_0}\left( \lambda \right)A\Delta Z}}{{{Z^2}}}\left[ {{\beta _m}\left( {z,\lambda } \right) + {\beta _\alpha }\left( {z,\lambda } \right)} \right]\nonumber\\ & \quad \exp \left\{ { - 2\mathop \smallint \limits_{{z_0}}^z \left[ {{\alpha _\alpha }\left( {z,\lambda } \right) + {\alpha _m}\left( {z,\lambda } \right)} \right]dz} \right\}\end{align} \tag{ 1 }$$

where $P\left( {z,\lambda } \right)$ is the received aerosol Lidar signal at the height of z, $\lambda$ is the wavelength of laser, K is the correction constant related to the Lidar system, ${P_0}\left( \lambda \right)$ is the power to launch a laser beam, A is the optical receiving area of the receiving telescope, and ${z_0}{\text{ }}$ is the height of the Lidar system. ${\beta _\alpha }\left( {{\text{z}},\lambda } \right)$ and ${\beta _m}\left( {{\text{z}},\lambda } \right)$ are the aerosol backscatter coefficient and molecular backscatter coefficient, respectively. ${\alpha _\alpha }\left( {{\text{z}},\lambda } \right)$ and ${\alpha _m}\left( {{\text{z}},\lambda } \right)$ the aerosol extinction coefficient and molecular extinction coefficient, respectively.

Previous studies have shown that compared with the Collis method (Collis and Russell 1976) and the Klett method (Klett 1981), the Fernald method (Fernald 1984) was the most stable, mature, and accurate inversion method. Therefore, the extinction coefficient of aerosols was retrieved using the Fernald method. The Fernald method assumes that the scattering of particles is proportional to the extinction coefficient for both aerosols and molecules. As shown in the following equations (2) and (equation (3)), where ${\text S_\text a}$ is the ratio of extinction coefficient and backscatter coefficient in atmospheric aerosol (usually also referred to as the Lidar ratio), which was set as 50 Sr in this study. Moreover, ${\text S_\text m}$ is the ratio of extinction coefficient and backscatter coefficient in atmospheric molecules, which was calculated as $8\pi /3$ according to the atmospheric molecular model

$$\begin{equation}{{\text{S}}_{\text{a}}} = \frac{{{{{\alpha }}_{{\alpha }}}\left( {{{\text{R}}_{\text{c}}}} \right)}}{{{{{\beta }}_{{\alpha }}}\left( {{{\text{R}}_{\text{c}}}} \right)}}\end{equation} \tag{ 2 }$$

$$\begin{equation}{{\text{S}}_{\text{m}}} = \frac{{{{{\alpha }}_{\text{m}}}\left( {{{\text{R}}_{\text{c}}}} \right)}}{{{{{\beta }}_{\text{m}}}\left( {{{\text{R}}_{\text{c}}}} \right)}}.\end{equation} \tag{ 3 }$$

Therefore, the equation (equation (1)) for the aerosol extinction coefficient can be written as:

$$\begin{equation}{{{\alpha }}_{\text{a}}}\left( {\text{z}} \right) = - \frac{{{{\text{S}}_{\text{a}}}}}{{{{\text{S}}_{\text{m}}}}} \cdot {{{\alpha }}_{\text{m}}}\left( {\text{z}} \right) + \frac{{{\text{P}}\left( {\text{z}} \right){{\text{Z}}^2} \cdot {\text{exp}}\left[ {2\left( {\frac{{{{\text{S}}_{\text{a}}}}}{{{{\text{S}}_{\text{m}}}}} - 1} \right)\mathop \smallint \nolimits_{\text{z}}^{{{\text{z}}_{\text{c}}}} {{{\alpha }}_{\text{m}}}\left( {{\text{z}}\prime} \right){\text{dz}}\prime} \right]}}{{\frac{{{\text{P}}\left( {\text{z}} \right){{\text{Z}}^2}}}{{{{{\alpha }}_{\text{a}}}\left( {{{\text{z}}_{\text{c}}}} \right) + \frac{{{{\text{S}}_{\text{a}}}}}{{{{\text{S}}_{\text{m}}}}}{{{\alpha }}_{\text{m}}}\left( {{{\text{z}}_{\text{c}}}} \right)}} + 2\mathop \smallint \nolimits_{\text{z}}^{{{\text{z}}_{\text{c}}}} {\text{P}}\left( {{\text{z}}\prime} \right){{\text{z}}^{{{\prime}}2}} \cdot {\text{exp}}\left[ {2\left( {\frac{{{{\text{S}}_{\text{a}}}}}{{{{\text{S}}_{\text{m}}}}} - 1} \right)\mathop \smallint \nolimits_{\text{z}}^{{{\text{z}}_{\text{c}}}} {{{\alpha }}_{\text{m}}}\left( {{\text{z}}\prime\prime} \right){{dz}}\prime\prime} \right]{\text{dz}\prime}}}\end{equation} \tag{ 4 }$$

where ${{\text{z}}_{\text{c}}}$ is the reference point, which height is the minimum aerosol concentration and has a high enough signal-to-noise ratio.

2.3.2. PM_2.5 mass concentration profile.

Previous studies have shown that the extinction coefficient of aerosols is closely related to the concentration of particulate matter (Li and Han 2016, Tao et al 2016). Based on single-scattering theory, the formula for converting the aerosol extinction coefficient (${{\alpha }}$) profile to particle mass concentration (PM_2.5) can be expressed as follows (Tao et al 2016, Lv et al 2017b):

$$\begin{equation}\text P{\text M_{2.5}} = \frac{{4\rho {r_{eff}}\eta }}{{3{Q_{ext}}}}\alpha = \delta *\alpha,\end{equation} \tag{ 5 }$$

where $\delta = \frac{{4\rho {r_{eff}}\eta }}{{3{Q_{ext}}}}$ is the proportion coefficient determined by the aerosol size distribution, components, and optical refractive index. In addition, $\rho$ in equation (5) represents the aerosol mass density, ${r_{eff}} = \frac{{\mathop \smallint \nolimits_0^\infty {r^3}n\left( r \right)dr}}{{\mathop \smallint \nolimits_0^\infty {r^2}n\left( r \right)dr}}$ depicts the effective radius of the particle, and n(r) is the aerosol size distribution,${\text{ }}{Q_{ext}}$ represents the extinction efficiency factor, and $\eta = \frac{{\text P{\text M_{2.5}}}}{{\text P{\text M_{total}}}}$ is the ratio of the PM_2.5 mass concentration to the total mass concentration.

Considering the unstable factors in the measurement (e.g. system deviation), a constant C was added to equation (5). Thus, the final formula for the PM_2.5 mass concentration in the linear model can be expressed as

$$\begin{equation}{\text{ }}\text P{\text M_{2.5}} = {{\delta *\alpha }} + {\text{C}}{\text{.}}\end{equation} \tag{ 6 }$$

Therefore, in the linear fitting of the near ground extinction coefficient (α) provided by the Lidar and the ground surface PM_2.5 mass concentration observed synchronously in situ, the slope and intercept were obtained by the fitting corresponding to the values of δ and C in equation (6).

In order to evaluate the accuracy of the model, the data from the two most polluted observation stations (Shijiazhuang City and Baoding City) were selected. The fitting results between the extinction coefficient at 0.1 km and the measured surface PM_2.5 concentration are presented in figure 2, revealing that the correlations were 0.89 and 0.87, respectively, and the corresponding root-mean-square errors (RMSEs) were 24.82 and 20.57 μg m⁻³, indicating that the model is reliable.

Zoom In Zoom Out Reset image size

2.4.1. Modeling of regional PM_2.5 production

2.4.2. Data assimilation system

The GSI DA system provides a 3DVAR analysis by minimizing the cost function as shown below (Gao et al 2017):

$$\begin{align}J\left( x \right) & = {\left( {x - {x_b}} \right)^T}{B^{ - 1}}\left( {x - {x_b}} \right) \nonumber\\ & \quad + {\left( {y - H\left( x \right)} \right)^T}{R^{ - 1}}\left( {y - H\left( x \right)} \right).\end{align} \tag{ 7 }$$

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 quantities, which is achieved via interpolation.

The hourly profile PM_2.5 concentrations of the two sets of Lidar data within 1 h windows were assimilated. The observation errors of the Lidar data originated from measurement errors and representative errors. The measurement error was computed using ${\varepsilon _0} = 1.5 + 0.0075{\text{*obs}}$ (Schwartz et al 2012, Jiang et al 2013, Gao et al 2017), where obs indicates the observed values. The representative error was computed using ${\varepsilon _r} = \gamma \sqrt {{{\Delta }}x/L}$ (Pagowski et al 2014), where $\gamma$ is the adjustable scale factor, ${{\Delta }}x$ is the model grid resolution (we selected 4 km for the domain), and L is the influencing radius (we used 60 km for the Lidar).

Figure 3 shows a comparison of the non-DA, DA, and observed values, in which the observed values of PM_2.5 were provided by the China National Environmental Monitoring Center (http://106.37.208.233:20035). These results clearly indicate that the Lidar data were a significant improvement over the simulation results. The correlation values of the non-DA and DA data were 0.74 and 0.86, respectively, representing an increase of 16.2%. Additionally, the RMSEs of the non-DA and DA data were 24.62 and 15.60 µg m⁻³, respectively, representing a decrease of 36.6%.

Zoom In Zoom Out Reset image size

2.4.3. Retrieval method for TF and TFI

Equation (8) was used to calculate the TF of the PM_2.5 based on the revised PM_2.5 concentration after assimilation and the wind field data from the simulation:

$$\begin{equation}{\text{Flux}} = {\text{ }}\text P{\text M_{2.5}}{\text{*}}WS{\text{*}}Cos\left( {\left( {WD - dir} \right){\text{*}}\frac{{{\pi }}}{{180}}} \right).{\text{ }}\end{equation} \tag{ 8 }$$

In this equation, PM_2.5 (µg m⁻³) is the concentration of particulate matter at the current height; WS (m/s) and WD are the wind speed and wind direction at the same position, respectively; $dir$ represents the angle between the current position and the direction of the transport pathway; and Flux (µg m⁻² s⁻¹) represents the TF of the PM_2.5 at the current position. Thus, positive values of Flux signify input fluxes, while negative values signify output fluxes. The TFI of the PM_2.5 was calculated by integrating the TF of the PM_2.5.

Figure S2 shows the wind field distribution at 4:00 p.m. on 26, 27, and 28 November 2017, which had respective winds from the southeast, southwest, and northwest, and these wind directions at the surface were consistent with those of the upper air. On the North China Plain, the average surface wind speeds for the 3 d were 3.3, 5.0, and 7.6 m s⁻¹, respectively, while the average wind speeds 1 km above the surface were 5.1, 10.7, and 10.9 m s⁻¹, respectively. For Beijing, the average surface wind speeds on the 26th and 27th were < 2 m s⁻¹, causing pollutants to accumulate (Slater et al 2020). Correspondingly, the average upper air wind speeds were > 6 m s⁻¹, which was conducive to the transport of pollutants at upper levels. Meanwhile, the overall wind speeds at both the surface and upper levels in Beijing were relatively high on the 28th, which was conducive to the diffusion of pollutants by the northwest winds.

The Fernald method was used to invert the Lidar data, thereby obtaining the extinction coefficient profile (Lv et al 2016). As shown in figure 4, the particles were mainly concentrated below 1.8 km during the observation period, appearing in decoupled layers. Additionally, the data obtained by the two Lidar systems exhibited good consistency. The maximum extinction coefficient was 1.2 km⁻¹, which occurred at 0.8 km. The stratification of the upper-air and surface pollutants had disappeared by 5:00 p.m.

Zoom In Zoom Out Reset image size

In order to further examine the paths of pollutant transport, figure 5 shows the observation results of the extinction coefficient for the five ground-based Lidars on the North China Plain from 26–29 November 2017. We found that pollutants began to accumulate in Shijiazhuang, Baoding, and other areas in the early morning of the 26th. Meanwhile, under the influence of the southwest wind, pollutants began to travel in the direction of Shijiazhuang, Baoding, and Beijing, with the transport height increasing from 0.6 km in the morning to 1.6 km at night, resulting in high concentrations of pollution in Beijing, Langfang, and other cities. The height of the pollutants began to decrease during the day on the 27th, with the pollutants reaching their peak values on the night of the 27th and in the early morning of the 28th, and the maximum extinction coefficient exceeding 1.5 km⁻¹. During the day on the 28th, under the action of northwest winds, pollutant advection began from north to south. The pollutants in Beijing, Tianjin, Langfang, and Baoding were successively reduced, while the pollutant diffusion conditions in Shijiazhuang were relatively poor due to the low wind speed.

Zoom In Zoom Out Reset image size

Figure 6 presents the spatial distribution of the PM_2.5 concentration after DA using the two sets of Lidar data. At 2:00 p.m., the PM_2.5 appeared in decoupled layers, of which the heavily polluted layer was 200–700 m. At the same time, a layer of low-concentration pollution appeared, ranging from 700–1700 m. These findings are consistent with the Lidar observation results, confirming that the results of the Lidar assimilation system were correct. At 2:00 p.m., the PM_2.5 concentration in northwest Beijing was high, with an upper air concentration of 135 µg m⁻³; the surface concentration, however, was only about 30 µg m⁻³. When combined with the meteorological results from the WRF-Chem simulation, an updraft near the ground was revealed, which caused the pollutants to spread out into the ambient atmosphere. Hence, the concentration of pollutants near the ground was low. Under the influence of the southwest wind flow, pollutants ahead of the Taihang Mountains were transported to Beijing, causing the PM_2.5 concentration in southwest Beijing to reach its maximum at 5:00 p.m. At the same time, the pollutant height decreased, causing the PM_2.5 concentration at the surface to reach 126 µg m⁻³. Meanwhile, in order to further estimate the possible sources and transport pathways of the aerosols, the backward and forward trajectory analyses from the HYSPLIT model are presented in figure S3.

Figure 7 shows the TF results of PM_2.5 along the 6th Ring Road in Beijing. A positive TF indicates that PM_2.5 was transported to Beijing; otherwise, the pollutants were exported from Beijing. Both the southern portion and the western portion of the 6th Ring Road primarily received PM_2.5 during the observation period, with the transmission occurring at 1 km. At the same time, the transmission intensity along the southern part of 6th Ring Road was higher than that along the western part of 6th Ring Road. In the northeast region of Beijing, PM_2.5 was mainly transmitted to other cities, and the transmission height was similar to that of PM_2.5 transported to Beijing in the southwest region, although the output intensity was significantly lower than the input intensity. The maximum PM_2.5 transmission values of 931 and − 701 µg m⁻² s⁻¹ in the southern and eastern parts of the 6th Ring Road, respectively, both occurred at 4:00 p.m. Moreover, figure S4 shows the TFI results of PM_2.5 along the 6th Ring Road in Beijing. The TFI values of PM_2.5 in the eastern, southern, western, and northern parts were −1133, 847, 1060, and −457 µg km⁻¹ s⁻¹, respectively, indicating that pollutants were transported from the southwest to Beijing and diffused from Beijing to the northeast.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size