2.1 AQM

The WRF–Chem model version 3.9.1 was selected as the AQM. The model has 40 vertical layers between the surface and 50 hPa, with the resolution gradually decreasing from the bottom up. The model domains are double-nested, and the second domain (D02) is centered at 114.57^∘ E and 37.98^∘ N and has 175×166 grid points with a grid interval of 9 km. D02 covers the central and eastern regions of China (most of North China, northern Central China, northern East China, and eastern Northwest China) (Fig. 1). The MOSAIC_4bin aerosol scheme was adopted for the simulations. This scheme, which will be described in Sect. 2.4, can be used to predict the profiles of eight aerosol types. For each aerosol type, there are four particle-size bins (4bins). The following summarizes the other physical and chemical schemes used in this study: the carbon-bond mechanism version Z (CBMZ) chemical reaction mechanism, the fast-J photolysis calculation scheme, the rapid radiative transfer model for general circulation model (RRTMG) shortwave radiation scheme, the RRTMG longwave radiation scheme, the WRF single-moment 5-class microphysical scheme, the unified Noah land-surface parameterization scheme, the Grell 3D ensemble cumulus parameterization scheme, the Yonsei University planetary boundary layer scheme, and the revised MM5 Monin–Obukhov near-surface layer scheme.

2.2 Data

The AEC profiles used in this study were derived from data captured by five conventional Mie scattering lidars (positioned in Beijing, Shijiazhuang, Taiyuan, Xuzhou, and Wuhu; Fig. 1) at a wavelength of 532 nm between 00:00 and 12:00 Coordinated Universal Time (UTC) on 13 November 2018 (Chen et al., 2019; Zhang et al., 2020). The temporal resolution of the data measured by the lidars in Shijiazhuang, Taiyuan, Xuzhou, and Wuhu was 1 min; that is, data were captured, and a vertical AEC profile was derived every minute. The vertical resolution of these data was 7.5 m; that is, one AEC was determined in one profile 7.5 m away from the next one. The blind zone of these lidars was 100 m; that is, these systems could not effectively capture AEC data between the surface and the height of 100 m. The temporal and vertical resolutions of the AEC profiles captured by the lidar in Beijing were 1 h and 15 m, respectively, and the blind zone of this lidar was 210 m. The relative standard deviation of the aerosol parameter profiles captured by the lidar over Beijing was 20.4 % in the height range of 1–2 km. This lidar was calibrated via comparative observation of several lidars (Chen et al., 2019). The precision of the AEC profiles released by the other four lidars was below the quality margins (25 % of the typical AEC observed in the planetary boundary layer or ±0.01 km⁻¹), as defined by Matthias et al. (2004). However, the relative standard deviation of the aerosol parameter profiles in the height range of 2–5 km released by lidar over Beijing was 35.9 %. To improve the effectiveness of the DA, it was necessary to first perform quality control on and preprocess the original AEC profiles. This ensured that the lidar data matched the numerical model in terms of temporal and spatial resolution. Quality control involved four steps:

1. Entire AEC profiles passing through low clouds and AEC measurements in mid- and high-cloud regions were eliminated. Clouds were defined as regions in which the AEC was higher than $\mathrm{5000}\times {\mathrm{10}}^{-\mathrm{6}}$ m⁻¹ (assuming the AEC in the near-surface layer (below 150 m) was lower than $\mathrm{3000}\times {\mathrm{10}}^{-\mathrm{6}}$ m⁻¹).

2. AEC profile data were subjected to maximum and minimum control. AEC measurements higher than $\mathrm{3000}\times {\mathrm{10}}^{-\mathrm{6}}$ m⁻¹ were each reassigned with a value of $\mathrm{3000}\times {\mathrm{10}}^{-\mathrm{6}}$ m⁻¹. AEC measurements lower than $\mathrm{20}\times {\mathrm{10}}^{-\mathrm{6}}$ m⁻¹ were eliminated.

3. For spatial continuity, data were required to be continuous within a vertical space L_con, which was set to be 90 m in this study. Specifically, two metrics were used to examine the spatial continuity of the data. First, the profile with vertical resolution L_res was examined. After the first two steps of quality control, the remaining number of data points (N_remain) within the L_con could not be less than one-third the total number of data points within the L_con ($N_{\text{total}}=L_{\text{con}}/L_{\text{res}}$). Otherwise, no valid data would be available for the center of the L_con. Second, the deviation of the valid data from the mean value of the data within the L_con could not exceed 3 times the standard deviation (SD).

4. Data within the blind zone of a lidar were eliminated. In addition, because lidar signals are relatively weak and AMCs are extremely low above 5000 m, data for the region above 5000 m were also eliminated in this study. After the quality control process, 84.32 % of the original AEC data from the lidar over Beijing were accepted as valid data, and 88.75 %, 54.10 %, 26.74 %, and 10.95 % of the data from the Taiyuan, Wuhu, Shijiazhuang, and Xuzhou lidars, respectively, were valid.

Preprocessing of quality control-treated AEC profiles involved two steps:

1. Temporal and spatial smoothing. Profiles were subjected to moving averaging over 30 m in the vertical direction. Temporally, the AEC profiles were averaged over the previous hour.

2. Data thinning. If there were multiple data points between two adjacent model layers in the vertical direction, only one was selected for assimilation. In this study, the nearest data point below each model layer was selected for assimilation. After processing, the number of assimilated AEC measurements per profile did not exceed 25, as there were no more than 25 model layers between the top of the lidar blind zone and the height of 5000 m.

PM_2.5 and PM₁₀ data (hereinafter referred to as PM data) used in this study, including 1 h MC data collected at more than 1500 ground environmental monitoring stations, originated from the China National Environmental Monitoring Center. Most of the monitoring stations were distributed in cities in economically developed regions, including the Yangtze River Delta, the Beijing–Tianjin–Hebei region, and the Pearl River Delta. Of these monitoring stations, more than 790 were located within the D02 region (Fig. 1). The assimilated PM data were collected between 00:00 and 12:00 UTC on 13 November 2018. After assimilation, forecasts for PM_2.5 from 12:00 UTC on 13 November 2018 to 12:00 UTC on 14 November 2018 were produced. In addition, the effects of DA on the forecast performance of the model were evaluated based on surface PM_2.5 measurements. To improve the DA performance and the representativeness of the evaluation metrics, the original PM data were subjected to quality-control and preprocessing treatments. Quality control involved two main steps:

1. Anomalous elimination. Measurements that remained unchanged over a continuous period of 24 h were considered anomalous and removed.

2. Maximum and minimum control. PM_2.5 MC measurements higher than 600 µg m⁻³, PM₁₀ MC measurements higher than 1200 µg m⁻³, and PM MC measurements less than 0 were considered anomalies and were removed.

During the DA and verification processes, there could be multiple PM MC measurements for one grid cell. To allow the measurements to represent the average PM MC within a certain area, the PM data used for DA and verification were subjected to grid-cell averaging. The PM data used for assimilation were averaged within 5×5 grid cells. Specifically, the PM data within the same 5×5 grid cell area were first examined to determine their spatial consistency. Data greater than twice the SD were removed. Next, the arithmetic mean of the data within the area was calculated and assimilated. The PM_2.5 MC measurements used for verification and model forecasts were averaged within 1×1 grid cells. Specifically, model forecasts were first interpolated to the location of each ground environmental monitoring station. Next, the arithmetic mean of the measured and forecasted values within the same grid cell was calculated and used as a sample for quantifying the evaluation metrics. The processed PM MC data for the D01 and D02 regions were assimilated, while only the PM_2.5 MC data for the D02 region were used to evaluate the effects of the DA. After the grid-cell averaging treatment, approximately 190 data points in the D02 region were assimilated each time.

2.3 Basic theoretical DA model

To mathematically achieve 3-DVAR DA, it is necessary to establish an objective function to transform the DA problem to a problem of finding the extreme values of the function. By calculating the extreme values of the function using the variational method, an “optimal” analysis field is obtained. The following shows the mathematical form of such a function:

This function describes the sum of the distance between the analysis field (x) and the background field (x^b) and the distance between the analysis field (x) and the observation field (y), with the background error covariance B and the observation error covariance R as weights, respectively. In Eq. (1), x is the control variable in the DA system, which is a one-dimensional (1D) vector composed of aerosol variables at all the 3D grid cells in the DA analysis field; x^b is the background value (or best guess) of the control variable (as the forecast level of AQM increases, model forecasts are generally used as background fields); B is the background error covariance; y is the observation variable, which is a 1D vector composed of all the measurements; H is the observation operator, which maps the control variable to the observation space to ensure that the observation data can provide observation information for the control variable even if they are not direct measurements of the control variable; and R is the observation error covariance. For simultaneous assimilation of two or more types of observation data, the second term on the right side of Eq. (1) can be expanded to multiple terms, each of which corresponds to one type of observation data. This will facilitate the simultaneous assimilation of observational data from various sources.

2.4 Control variables and B

The MOSAIC_4bins aerosol scheme adopted in this study accommodated eight aerosol types, namely, black ∕ elemental carbon (EC ∕ BC), organic carbon (OC), sulfates (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), nitrates (NO${}_{\mathrm{3}}^{-}$), ammonium salts (NH${}_{\mathrm{4}}^{+}$), chlorides (Cl⁻), sodium salts (Na⁺), and other unclassified inorganic compounds (OIN). There were four particle-size bins (4bins) for each aerosol type, namely, 0.039–0.1, 0.1–1.0, 1.0–2.5, and 2.5–10 µm. Thus, there were 32 model variables that represented the various aerosols. However, limitations in computer memory and computational capacity necessitated a reduction in the total number of control variables. In addition, the AECs of fine (PM_2.5) and coarse (PM_2.5−10) particles differed significantly. Thus, two control variables for each aerosol type were designed – one corresponding to fine particles (formed by combining the first three particle-size bins) and one corresponding to coarse particles (the fourth particle-size bin). Thus, there were 16 control variables in the DA scheme, namely, EC_2.5, EC_2.5−10, OC_2.5, OC_2.5−10, SO${}_{{\mathrm{4}}_{\mathrm{2.5}}}$, SO4_2.5−10, NO${}_{{\mathrm{3}}_{\mathrm{2.5}}}$, NO3_2.5−10, NH${}_{{\mathrm{4}}_{\mathrm{2.5}}}$, NH4_2.5−10, CL_2.5, CL_2.5−10, NA_2.5, NA_2.5−10, OIN_2.5, and OIN_2.5−10.

There were two problems associated with calculations involving B:

1. In this scheme, B contained 3.5×10¹⁴ (square of 16 (number of control variables) $\times \mathrm{175}\times \mathrm{166}\times \mathrm{40}$ (number of grid cells)) elements. Thus, it was necessary to mathematically treat and simplify B to facilitate numerical calculations. Following the method used by Li et al. (2013) and Zang et al. (2016), B was decomposed into a background-error standard deviation (BESD) matrix, a background-error horizontal correlation coefficient (BEHCC) matrix, and a background-error vertical correlation coefficient (BEVCC) matrix for calculations.

2. As the true value of B was unknown, it was necessary to develop a reasonable statistical method to estimate it. The National Meteorology Center (NMC) method (Parrish and Derber, 1992) was employed in this study to statistically estimate B. Specifically, the differences between the 48 and 24 h forecasts of the control variables were assumed to be a proxy of the background error. Next, B was estimated based on the covariance of the difference field, which was obtained by producing continuous 24 and 48 h forecasts for a month using the WRF–Chem model.

2.5 Observation operator and its adjoint

Obtaining the observation operator involved two calculations:

1. The control variables at each grid cell were mapped to the observation space; that is, the control variables were mapped to the AEC values (or PM_2.5 and PM₁₀ MCs).

2. The mapped values at the eight vertices of the model grid cell associated with the observation data were interpolated using the inverse distance-weighted method to the observation location. Here, we only describe the first step of the derivation of the observation operators, which are different for different observation data.

The AEC observation operator was based on the IMPROVE equation. The following shows the specific form of the IMPROVE equation:

The left side of Eq. (2) is the AEC value Ext (unit: 10⁻⁶ m⁻¹). The variables in the brackets on the right side of Eq. (2) are combinations of the 16 control variables (unit: µg m⁻³). The coefficient variables f_s(RH), f_l(RH), and f_ss(RH) reflect the effects of hygroscopicity of fine, coarse, and sea-salt aerosols, respectively, under various relative humidity (RH) conditions. The values of the parameters given by Gordon et al. (2018) were used in this study. The variables (in square brackets) at each grid cell were obtained by combining the 16 control variables using the following method:

The principle for determining α involved preferentially allocating NH${}_{{\mathrm{4}}_{\mathrm{2.5}}}$ to SO${}_{{\mathrm{4}}_{\mathrm{2.5}}}$. The remaining NH${}_{{\mathrm{4}}_{\mathrm{2.5}}}$ was allocated to NO${}_{{\mathrm{3}}_{\mathrm{2.5}}}$.

The observation operators for PM_2.5 and PM₁₀ were the sums of control variables in the corresponding particle-size bin, that is,

The corresponding adjoint process on the operators for PM and AEC were developed and passed the adjoint sensitivity test. For the adjoint test method, please refer to Zou et al. (1997).

Table 1Numerical experiment schemes. n/a – not applicable.

Download Print Version | Download XLSX

2.6 DA and forecast experimental design and verification analysis method

To analyze the effects of DA on aerosol analysis and forecasts, one control experiment and three DA experiments were designed for a pollution event that occurred from 13 to 14 November 2018 (Table 1). In the control experiment, no chemical observation data were assimilated. Forecasts were produced for a 36 h period, starting at 00:00 UTC on 13 November 2018. In the DA experiments, aerosol data were assimilated every hour for the DA period of 00:00–12:00 UTC on 13 November 2018. Next, with the analysis field obtained from the DA as the initial chemical field, forecasts were performed for a 24 h period starting at 12:00 UTC on 13 November 2018. For the first DA cycle in each of the three DA experiments, the initial field of the control experiment was used as the background field, the observation data for 00:00 UTC on 13 November 2018 were assimilated, and a DA analysis field was generated. With this DA analysis field as the initial field at 00:00 UTC on 13 November 2018 in the DA experiment, 1 h forecasts were produced. The forecasts produced for 01:00 UTC on 13 November 2018 were used as the background field for the second DA cycle. The process was repeated for 13 assimilation cycles. Thus, a DA analysis field for 12:00 UTC on 13 November 2018 was generated. The effects of DA on forecast performance during the forecast comparison period from 12:00 UTC on 13 November 2018 to 12:00 UTC on 14 November 2018 were analyzed by comparing the forecast performance of the DA and control experiments. In the first DA experiment, PM data alone were assimilated (DA_PM). In the second DA experiment, the lidar data alone were assimilated (DA_Ext). In the third DA experiment, PM and lidar data were assimilated simultaneously (DA_PM_Ext). Furthermore, ${\mathrm{0.25}}^{\circ }\times {\mathrm{0.25}}^{\circ }$ 6 h reanalysis data provided by the US National Centers for Environmental Prediction (NCEP) were used as the meteorological field of the model.

Two metrics, the regional mean and root-mean-square error (RMSE), were used to evaluate simulation and forecast accuracy of the PM_2.5 MC in the experiments. The closer the mean of the simulated values to the mean of the measurements and the smaller the RMSE, the higher the performance. Let M_i, O_i, N, $\stackrel{\mathrm{‾}}{M}$, and $\stackrel{\mathrm{‾}}{O}$ be the simulated value sample, the measured value sample, the number of samples, the mean of simulated values, and the mean of the measurements, respectively. The following summarizes the equations for calculating the metrics:

Figure 2Vertical BESD profiles of the 16 control variables.

Download

Figure 3BEVCCs of six control variables.

Download

3.1 BESD and BEVCC

Under the same conditions, the larger the BESD, the larger the DA increment field (the difference between the “optimal” analysis field and the background field). Therefore, the structural pattern of the BESD significantly affected the distribution pattern of the DA increment field. The vertical BESD profiles of the 16 control variables are shown in Fig. 2. The BESD differed significantly among the control variables. The seven control variables with the largest BESDs below the height of 1000 m (corresponding to the 22nd layer of the model) in descending order of BESD were OIN_2.5−10, NO${}_{{\mathrm{3}}_{\mathrm{2.5}}}$, OIN_2.5, NH${}_{{\mathrm{4}}_{\mathrm{2.5}}}$, SO${}_{{\mathrm{4}}_{\mathrm{2.5}}}$, OC_2.5, and EC_2.5. As height increased, the BESD of each control variable decreased. The rates of decrease were the highest above the boundary layers at heights of 1000–2000 m (corresponding to the 20th–25th layers of the model). The BEVCC matrix can spread the observation information contained in measurements around one model layer to nearby vertical layers. Therefore, even if the PM data are only available at the surface, there will still be increments of PM near the surface (in-air) after DA. Furthermore, even though the lidar AEC data are not available at the surface, assimilating lidar data can still correct the surface PM_2.5 MC distribution. Figure 3 shows the BEVCC matrices of six control variables with relatively large BESDs (OIN_2.5−10, NO${}_{{\mathrm{3}}_{\mathrm{2.5}}}$, OIN_2.5, NH${}_{{\mathrm{4}}_{\mathrm{2.5}}}$, SO${}_{{\mathrm{4}}_{\mathrm{2.5}}}$, and OC_2.5). The BEVCCs of the six control variables share certain common characteristics. The correlation decreases as the interlayer spacing of the model increases. Each in-air layer is positively correlated with the surface layer, although the correlation decreases as height increases. For OIN_2.5−10, the correlation coefficient between the surface and 10th layers is 0.34, compared with 0.49–0.51 for other variables. This indicates that OIN_2.5−10 has a significantly weaker vertical correlation, and hence DA increments of these particles settle more rapidly than the other variables do. This is mainly because coarse particles settle faster vertically than fine particles and are concentrated near the surface in larger quantities.

Figure 4Surface PM_2.5 MC measurements in the D02 region and NCEP reanalysis wind field for the period from 00:00 UTC on 13 November 2018 to 12:00 UTC on 14 November 2018 (Bu-O: Buyant-Ovoo; Ba-O: Bayan-Ovoo; CY: Chaoyang; BJ: Beijing; SJZ: Shijiazhuang; TY: Taiyuan; ZB: Zibo; X'A: Xi'an; ZZ: Zhengzhou; NJ: Nanjing).

Figure 5AEC profiles measurements (black lines), the AEC profiles in the analysis fields of the control (blue lines), DA_PM (green lines), DA_Ext (purple lines), and DA_PM_Ext (red lines) experiments and the simulated RH profiles (orange lines) at four lidar stations at 00:00 UTC on 13 November 2018 (BJ: Beijing; SJZ: Shijiazhuang; TY: Taiyuan; WH: Wuhu).

Download

3.2 Analysis of the pollution process

Figure 4 shows the evolutionary process of the surface PM_2.5 MC and the NCEP reanalysis surface wind field in the D02 region for the period from 00:00 UTC on 13 November 2018 to 12:00 UTC on 14 November 2018 (the time interval between Fig. 4a–d is 12 h). At 00:00 UTC on 13 November 2018, the D02 region was predominantly controlled by a high-pressure circulation centered over Zibo. There was a clockwise wind field around the high-pressure center. Therefore, the northerlies (easterlies) east (south) of the high-pressure center brought clean air over the sea land-ward. As a result, the PM_2.5 MCs over East China were relatively low. For example, the mean PM_2.5 MC measured at the ground environmental monitoring stations in Nanjing was 41.8 µg m⁻³. There were relatively slow southerlies west and northwest of the high-pressure center, which led to favorable conditions for pollutant accumulation east of the Taihang Mountains and south of the Yan Mountains. As a result, North China was heavily polluted by PM_2.5. For example, the mean PM_2.5 MCs in Beijing and Shijiazhuang were 122.7 and 149.3 µg m⁻³, respectively. In addition, within the D02 region, there was a northeast–southwest-trending cold front near Buyant-Ovoo–Bayan-Ovoo in Mongolia. As time passed (Fig. 4b–d), the high-pressure center gradually moved northeastward and reached near the eastern boundary of the D02 region by 12:00 UTC on 14 November 2018 (Fig. 4d). The cold front gradually moved southeastward and reached the Chaoyang–Beijing–Taiyuan–Xi'an line by 12:00 UTC on 14 November 2018 (Fig. 4d). As the high-pressure center and the cold front moved, the level of pollution in North China continued to rise, and pollution gradually expanded northeastward to Chaoyang, southward to Zhengzhou, and westward to Taiyuan. The level of pollution gradually increased in the Wei and Yellow river valleys east of Xi'an due to the dual action of advection by the easterlies and the narrow terrain, while the PM_2.5 MCs decreased considerably with the passing of the cold front due to the good dispersion conditions. There were no significant changes in the PM_2.5 MCs in East China due to the continuous impact of sea winds.

3.3 Analysis of the direct effects of DA

Figure 5 shows the AEC profile measurements, the AEC profiles in the analysis fields of the control and DA experiments, and the simulated RH profiles at four lidar stations at 00:00 UTC on 13 November 2018, when the first DA cycle was performed. The results of the control experiment were used as the background field in the three DA experiments. Figure 5a–d show the results for Beijing, Shijiazhuang, Taiyuan, and Wuhu, respectively. As the in-air RH profile (brown lines) below 1 km was basically consistent with that of the surface RH, the vertical changes in the AEC values in this region were only slightly affected by the RH. Thus, the AEC profiles were used to study the vertical changes in the PM_2.5 MC. For Beijing, the simulated AEC results from the control experiment (blue lines) agreed with the lidar AEC measurements well (Fig. 5a – black lines). However, for Shijiazhuang and Taiyuan, the simulation underestimated the empirical results (Fig. 5b and c, respectively), particularly near the height of 100 m (the lowest height of valid lidar data), while for Wuhu, it overestimated them (Fig. 5d).

The DA increments of AEC values from the DA_PM, that is, the AEC values obtained from the DA_PM experiment (green lines) minus those from the control experiment (blue lines), were negative for Beijing (Fig. 5a), Taiyuan (Fig. 5c), and Wuhu (Fig. 5d) at the surface. They were also negative from the near-surface to a height of about 1000 m, although their absolute values were smaller than those at the surface. This is because the BEVCCs between each in-air layer and the surface layer were positive and decreased with height (Fig. 3) so that the information contained in the surface PM MC measurements was spread to the air. However, the results of the adjustment of the AEC profiles were not always positive, because the aerosol bias of the control experiment at the surface was not always the same as it was in the atmosphere. Thus, they were overall positive for Beijing and Wuhu but negative for Taiyuan, reflecting the fact that the PM DA did not effectively account for the vertical aerosol distribution adjustment.

Compared to those from the DA_PM experiments, the AEC values from the DA_Ext experiments (purple lines) for Taiyuan (Fig. 5c) at heights of approximately 100 and 700 m were significantly larger than those from the DA_PM experiment and were consistent with the measurements (black line), and those for Wuhu (Fig. 5d) were very close to the measurements across the entire profile. This suggests that the AEC observation operator whose design was based on the IMPROVE equation effectively facilitated 3D variational assimilation of lidar AEC data. In addition, although lidar data were not available at the surface, the DA_Ext adjusted the surface PM MCs, corrected the overestimation of surface PM_2.5 MCs in Beijing and Wuhu but increased the overestimation of surface PM_2.5 MCs in Taiyuan. This is because the information contained in the in-air AEC was spread to the surface, while the aerosol bias of the control experiment in the air did not always match that at the surface.

The in-air AEC profiles obtained from the DA_PM_Ext experiment (red lines) for the four cities almost coincided with those from the DA_Ext experiments above 400 m. The near-surface AEC values obtained from the DA_PM_Ext experiment for Beijing (Fig. 5a), Taiyuan (Fig. 5c), and Wuhu (Fig. 5d) almost coincided with those from the DA_PM experiment, were between those from the DA_PM and DA_Ext experiments, and were smaller than those from both the DA_PM and DA_Ext experiments. This suggests that simultaneously assimilating the two types of data can fully integrate their observation information and reflect their respective advantages, thereby generating the most accurate analysis field.

Figure 6AEC profiles measurements (solid black lines), the AEC profiles in the control experiment (solid blue lines), in the background field of the DA_PM (dotted green lines), DA_Ext (dotted purple lines), and DA_PM_Ext (dotted red lines) experiments, and in the analysis fields of the DA_PM (solid green lines), DA_Ext (solid purple lines), and DA_PM_Ext (solid red lines) experiments at four lidar stations at 12:00 UTC on 13 November 2018 (BJ: Beijing; SJZ: Shijiazhuang; TY: Taiyuan; WH: Wuhu).

Download

Figure 7Surface PM_2.5 MC measurements (a), surface PM_2.5 MCs in the initial field of control experiment (b) and its bias (c), the inverse DA increments of PM_2.5 MC of DA experiments, that is, the PM_2.5 MCs obtained from the control experiment minus that from the DA experiments (d–f) at 12:00 UTC on 13 November 2018 (black triangles signify the locations of the lidar stations, and green triangles mark the locations of the two cities without lidar) (CY: Chaoyang; BJ: Beijing; SJZ: Shijiazhuang; TY: Taiyuan; ZZ: Zhengzhou; XZ: Xuzhou; WH: Wuhu).

Figure 6 shows the AEC profiles measured, simulated by the control experiment, in the background fields and the analysis fields of the DA experiments at four lidar stations at 12:00 UTC on 13 November 2018. The time of 12:00 UTC on 13 November 2018 was the last time point of the DA period, the starting time point of the forecast period, and the time point at which 13 DA cycles had elapsed. The background field for each of the three DA experiments was generated during the continuous DA period, whereas the results of the control experiment were obtained by a 12 h forecast starting at 00:00 UTC on 13 November 2018. As a result, there was a significant difference between the background fields of the three DA experiments and those of the control experiment.

The DA increments of the AEC values from the DA_PM experiment were significant below 1000 m (green lines). These adjustments corrected the near-surface overestimation of the AEC values for the four cities in the control experiment; however, they increased the underestimation for Taiyuan at heights of 120–400 m (Fig. 6c) and overestimation for Wuhu above 400 m (Fig. 6d). Additionally, it is worth noting that there were small direct DA increments generated in the DA_PM experiment at this time point. This means that for the surface PM DA, a DA period of 11 h or less was sufficient to effectively adjust aerosol distribution in this experiment. This may be because aerosols were primarily concentrated near the surface and surface PM data covered a wide area and had a high spatial resolution. Thus, surface PM data measured at a few time points contained the main aerosol distribution information for the whole region.

Compared to the DA_PM experiment, the DA_Ext experiment (purple lines) reflected the advantages of adjusting the vertical aerosol distribution. The overestimations for Beijing above 300 m (Fig. 6a), Taiyuan above 600 m (Fig. 6c), and Wuhu below 400 m (Fig. 6d) in the control experiment were effectively corrected. The rapid decrease in the AEC from the surface to a height of 1000 m over Beijing (Fig. 6a) and the maximum-AEC layer at a height of 1300 m over Wuhu (Fig. 6d) were accurately reproduced by the DA_Ext experiment. However, the near-surface overestimation for Taiyuan (Fig. 6c) increased. Moreover, the direct DA increments generated in the DA_Ext experiment at this time point remained notable. This suggests that the background field errors at each lidar station at 12:00 UTC remained relatively large, even after the continuous DA period. To improve the effects of the DA, it was necessary to increase the length of the continuous DA period. This may have been due to the limited number of lidars and the fact that the lidars were relatively far apart from one another. Thus, the simulation error for the region upstream of a lidar was difficult to correct through DA and affected the lidar location due to the effects of advection at the next time point. In addition, because the 12:00 UTC (20:00 LST) was only 2–3 h after sunset, large changes of PM concentration profile may occur due to large changes in the planetary boundary layer height after sunset.

Figure 7 shows the surface PM_2.5 MC measurements, the surface PM_2.5 MCs of the initial field of the control experiment and their biases, and the inverse DA increments of PM_2.5 MCs from the DA experiments, that is, the PM_2.5 MCs obtained from the control experiment minus those from the DA experiments at 12:00 UTC on 13 November 2018. The measurements (Fig. 7a) showed that the PM_2.5 MCs were relatively high in North China, with a heavily polluted zone in the Beijing–Shijiazhuang–Zhengzhou region, while the PM_2.5 MCs were relatively low surrounding North China. The control experiment (Fig. 7b) successfully simulated regions with relatively high and low PM_2.5 MCs. However, the PM_2.5 MCs were overestimated for most stations in D02 (Fig. 7c), especially in the Beijing–Shijiazhuang–Zhengzhou region, and underestimated for stations near Chaoyang.

Figure 8Variation of the regional mean PM_2.5 MC over time measured and simulated by the four experiments (the vertical orange line separates the DA and forecast periods; the black line signifies measurements; the blue line signifies that obtained from the control experiment; the green, purple, and red lines signify that obtained from the DA_PM, DA_Ext, and DA_PM_Ext experiments, respectively).

Download

The inverse DA increments of the PM_2.5 MCs of the DA_PM experiment (Fig. 7d) were relatively consistent with the bias of the control experiment (Fig. 7c), indicating that the overestimation for most regions and the underestimation for some regions in the initial field of the control experiment were corrected by the PM DA. The inverse DA increments of the PM_2.5 MCs of DA_Ext (Fig. 7e) were significant in the regions surrounding and downstream of the five lidar stations. In addition, certain DA increments were also present in regions far away from the lidar stations. This indicates that long-term continuous lidar AEC DA can affect a relatively large area. Overall, the DA_Ext corrects the overestimation for most stations and underestimation for a few stations in the control experiment. However, the DA_Ext increments were smaller than the DA_PM increments in terms of horizontal spatial range and absolute values. This is mainly because there are relatively few lidars, and these lidars cover a limited spatial area. It is worth noting that DA_Ext yields a negative effect for northern Beijing and the region around Taiyuan, a result which will be discussed later in Sect. 4. The inverse DA increments of PM_2.5 MCs of DA_PM_Ext (Fig. 7f) were relatively consistent with those of the DA_PM (Fig. 7c). This is mainly because the quantity and spatial coverage of the PM data were larger and more complete than those of the lidar data. As a result, the DA increments of the surface PM_2.5 MCs originated primarily from the observation information contained in the PM data. Because the AEC profiles of the DA_PM_Ext almost coincided with those of the DA_Ext above 400 m (Fig. 5), the DA_PM_Ext reflected the 3D spatial distribution pattern of the aerosols most accurately.

3.4 Effects of DA on the forecast performance for surface PM_2.5 MCs

In this section, the forecast performances of the DAs for surface PM_2.5 are evaluated based on measurements that cover most of the D02 region.

Figure 8 shows the variation of the regional mean of the PM_2.5 MC over time from the four experiments. The regional mean of the PM_2.5 MC (black line) exhibited a notable diurnal pattern. Two notable minimum PM_2.5 MC values (69.1 and 77.9 µg m⁻³) appeared at 08:00 UTC (16:00 local time) on 13 and 14 November 2018, respectively. High PM_2.5 MCs appeared between 13:00 UTC on 13 November 2018 and 02:00 UTC on 14 November 2018 (from night to morning), with a maximum PM_2.5 MC of 96.0 µg m⁻³. Meanwhile, there was a relative minimum PM_2.5 MC (87.0 µg m⁻³) appearing at 22:00 UTC on 13 November 2018 (around dawn local time) during the high-PM_2.5-MC period.

The control experiment (blue line) simulated the periodic variation pattern of the mean PM_2.5 MC but significantly overestimated the value of this parameter during the entire forecast period. The mean PM_2.5 MC of the control experiment at the initial time for the forecast period (12:00 UTC on 13 November 2018) was 128.6 µg m⁻³, which is 36.3 µg m⁻³ (39.3 %) larger than that of the measurements (92.3 µg m⁻³). The DA_PM (green line, which almost coincides with the red line) significantly reduced the overestimation of the control experiment, with a mean PM_2.5 MC of 91.4 µg m⁻³ that is 0.9 µg m⁻³ (1.0 %) lower than the measurement. As a result of the decrease in the MC levels in the initial field, the PM_2.5 MC forecasts of the DA_PM were significantly lower than those of the control experiment during the entire forecast period. This suggests that the overestimation of the initial field is the primary cause of the overestimated forecasts of the control experiment. The overestimation of the control experiment at the initial time point was reduced by the DA_Ext (purple line) from 36.3 µg m⁻³ (39.3 %) to 20.5 µg m⁻³ (22.2 %), which improved the forecast performance significantly (even though there were only five lidars within the region). There was no significant difference between the results of the DA_PM_Ext (red line) and DA_PM (green line) at the surface. This suggests that in these experiments, after DA of surface PM data, the DA of lidar data did not significantly affect the surface PM_2.5 MC levels. There are two reasons for this. The PM data set was far larger than the lidar data set in terms of quantity and spatial coverage. In addition, after surface PM DA, lidar DA mainly directly adjusted the AMC values not at surface but in-air and hence affected the surface AMC forecasts only indirectly, via processes such as settling. However, in this simulation process, the surface AMC levels remained relatively high, while the vertical air movement was weak due to the relatively stable meteorological conditions, particularly in the heavily polluted zone. Therefore, the effects of the lidar DA on the surface PM_2.5 MCs are far smaller after the surface PM DA.

Figure 9Variation in the RMSE of surface PM_2.5 MC forecasts over time (the vertical orange line separates the DA and forecast periods; the blue line signifies that obtained from the control experiment; the green, purple, and red lines signify that obtained from the DA_PM, DA_Ext, and DA_PM_Ext experiments, respectively).

Download

Figure 9 shows the variation in the RMSE of surface PM_2.5 MC forecasts over time. The RMSEs for simulations and forecasts were relatively large (small) when the mean PM_2.5 MCs were relatively high (low) (Fig. 8). The RMSE in the control experiment was 59.6 µg m⁻³ at the initial time for the forecast period (12:00 UTC on 13 November 2018) and fluctuated between 44.5 and 67.1 µg m⁻³ instead of linearly increasing or decreasing throughout the forecast period. The RMSEs in the DA_PM (green line), DA_Ext (purple line), and DA_PM_Ext (red line) experiments at the initial time point were 21.0, 49.1, and 21.2 µg m⁻³, respectively, which were 38.6 µg m⁻³ (64.8 %), 10.5 µg m⁻³ (17.6 %), and 38.4 µg m⁻³ (64.4 %) lower than that of the control experiment. Owing to the optimized initial field, the RMSE of the forecasts of each of the DA experiments was lower than that of the control experiment during the forecast period. For the 24th forecast hour, the RMSEs of the forecasts of the Da_PM, Da_Ext, and DA_PM_Ext were 6.1 µg m⁻³ (11.8 %), 1.5 µg m⁻³ (2.9 %), and 6.5 µg m⁻³ (12.6 %) smaller than that of the control experiment, respectively. This suggests that the optimization of the initial field has a lasting (more than 24 h in all cases) positive effect on model forecasts. It is worth noting that while there are very few lidar stations, the results of the DA_Ext experiment were still better than those of the control experiment, and the results of the DA_PM_Ext experiment were also slightly better than those of the DA_PM experiment. This indicates that even in relatively low quantities, lidar data still improve the forecast performance of the model. As lidar data become increasingly rich and provide more vertical and horizontal aerosol distribution information in the future, lidar DA will further improve PM_2.5 MC forecasts.