Abstract

Based on the theory of inverse problem and data assimilation, the adjoint method is applied for the estimation of parameters including the initial condition (IC), the source and sink (SS) in a PM2.5 transport model. To reduce the ill-posedness of the inverse problem, an independent point scheme (IPS) is implemented during the estimation process. In twin experiments, both the prescribed IC and SS can be inverted successfully and better inversion results are obtained when the IPS is used than not, suggesting the feasibility and validity of the PM2.5 transport model as well as the IPS. In practical experiments, several inversion strategies are compared based on the simulation results of PM2.5 concentrations over China. It is found that IC and SS are better estimated with smaller difference between simulated results and observations, when IC and SS are inverted simultaneously than separately. And the simulated results can reproduce the temporal and spatial variation feature of the observed PM2.5 concentrations. On the basis of the numerical results, it is shown that the adjoint method and the IPS are the powerful way to improve the precision of the simulation of the PM2.5 concentrations.

1. Introduction

In China, a combination of rapid industrialization and high population density inevitably deteriorates the air pollution problems [1]. As the 2013 Report on the State of Environment in China shows, the number of hazy days is about 36, reaching a record high since 1961. In some areas the number is even up to 100. And when haze occurs, the PM2.5 (particles with aerodynamic diameters less than 2.5 μm) is usually reported as the principal pollutant. To control the haze pollution effectively and provide increased protection for public health, the Chinese government made the new ambient air quality standards in 2012 to which the PM2.5 standard was added. Naturally PM2.5 pollution has become a major concern in China in recent years, due to its possible impacts on health [26], visibility [7], and climate changes [8].

Numerical model is one of the most important ways to investigate PM2.5. And many studies about PM modeling have been conducted using Eulerian box models (e.g., [9]), Lagrangian trajectory models (e.g., [10]), Lagrangian plume models (e.g., [11]), and 3-dimensional Eulerian models (e.g., [12]). The aforementioned PM models have developed rapidly and become increasingly complex, and details about them are available in several recent review articles (e.g., [1315]). However, there are several indeterminate factors in the simulation, such as primary PM [1618], gas-phase chemistry [19, 20], and dry and wet deposition [21, 22]. McKeen et al. [21] evaluated seven PM2.5 forecast models using data collected during the ICARTT/NEAQS 2004 field campaign and found that the ensemble PM2.5 forecast, created by combining six separate forecasts with equal weighting, yielded the best possible forecast in terms of the statistical measures considered. And the ensemble forecast is also taken as the technique to improve forecasting skills in Zhang et al. [15]. In addition, Zhang et al. [15] indicated that the chemical data assimilation techniques could combine the predicting model with observations to improve the initial and boundary conditions as well as emissions and decrease inaccuracy in the predicting process. Furthermore, the inverse model using data assimilation could estimate the uncertain parameters in the simulation process. Previous studies have demonstrated the effectiveness of various inverting modeling method to obtain spatial distribution of pollution emissions. For large fields, the 4-dimensional variational method (4D-Var) is chosen. Quélo et al. [23] proved its feasibility on real emission factor of NOx. Black carbon emission over East Asia was estimated by Hakami et al. [24] using the adjoint STEM model. Yumimoto and Uno [25] applied 4D-Var to a chemical transport model and estimated carbon monoxide emissions over East Asia. Elbern et al. [26] carried out the large-scale inversions of the precursors of O3, SO2, NH3, and so forth, over Europe. And the adjoint method has been applied to several chemical transport models [2729]. In addition, the adjoint method has been used to study the PM2.5. Henze et al. [30] used the adjoint of GEOS-Chem to do the inverse modeling and mapping US air quality influences of inorganic PM2.5 precursor emissions. The source attribution of PM2.5 pollution over North China was examined by Zhang et al. [31] using the adjoint method. Capps et al. [32] used the GEOS-Chem adjoint to reveal the relative contributions of global emissions to PM2.5 air quality attainment in the US.

But the parameters estimation is often ill-posed and beset by instability and nonuniqueness, particularly if the parameters are distributed in space and time domain [3336]. The parameters in the general chemical transport model are often distributed complexly; even they are not understood fully. So reducing the number of control variables in the process of parameters estimation is a good way to increase the well-posedness of the inverse problem and improve the estimated results. In this paper, an independent point scheme (IPS) is employed in which a subset of the full grid points is chosen as the independent points and the parameters at the other gird points are obtained through interpolation of the values at the independent points.

In this paper, the parameters of a PM2.5 transport model are estimated using the adjoint method (i.e., 4D-Var) with the IPS and the PM2.5 concentrations are simulated over China. In Section 2, the adjoint PM2.5 transport model is constructed and the IPS is described. In Section 3, twin experiments are carried out to test the reasonability and feasibility of this model and verify the effectiveness of the IPS. Section 4 estimates the parameters using this model and simulates the PM2.5 over China. Section 5 summarizes the paper and provides some discussions.

2. Model

2.1. The PM2.5 Transport Model

Generally speaking, PM2.5 varies due to interactions among many processes including emissions (anthropogenic emission and natural dust production), transport (as well as convection-influenced dispersion and dilution), photochemical transformation (new particle speciation and production of secondary PM2.5), and deposition (dry and wet), with meteorology playing an overarching role [37]. Hence the general PM2.5 model includes a rich description of the photochemical oxidant cycle and a chemical mechanism to simulate the formation of secondary PM2.5, and the adjoint model is significantly difficult to build.

Just to verify the estimation ability of adjoint method and the superiority of IPS, we take the primary PM2.5 and secondary PM2.5 as a whole, calling it “source and sink” (SS) without considering the specific details. Considering that the observations are of ground level and the vertical meteorology environment is indeterminate, we establish a two-dimensional PM2.5 transport model in rectangular coordinates as follows:

Here represents the PM2.5 concentration, and are the horizontal wind velocity in -coordinate and -coordinate, respectively, is the horizontal diffusivity coefficient, and is the SS. The model has the initial conditions (ICs) and is subject to constant boundary conditions at the inflow boundary and to no gradient boundary conditions at the outflow boundary .

The finite difference scheme of (1) is as follows:

Here, the upwind scheme is used in the advection term; that is,

And it is similar for the term in -coordinate.

With the IC, the SS, and the background value which is the constant in (3), the PM2.5 transport model could be solved by a sequence of time steps of length in the discrete schemes (i.e., (5)).

2.2. The Adjoint Model

In order to construct the adjoint model, the cost function is defined as

And based on the theory of Lagrangian multiplier method, the Lagrangian function is defined as

Here is the set of the observations, and are the simulated and observed PM2.5 concentrations, respectively, denotes the adjoint variable of , and is the weighting matrix and theoretically should be the inverse of observation error covariance matrix. Assuming the errors of the data are uncorrelated and equally weighted, can be simplified [38], and in this study is 1 if observations are available and 0 otherwise.

To let the cost function reach the minimum, we make the first-order derivates of Lagrangian function with respect to all the variables and parameters be zero, as follows:

Actually (9) is (1). The adjoint equation can be derived from (10), which is in the following form:

Furthermore, the ICs and open boundary conditions are similar to (2)–(4).

The finite difference scheme of (12) is similar to (5), but it should be noted that the difference scheme of the advection term is opposite to (6) because the adjoint model is integrated backward.

From (11) and the difference scheme, the gradients of the cost function with respect to parameters including the horizontal diffusivity coefficient, the SS, and IC could be obtained and displayed as follows:

2.3. Independent Point Scheme

In oceanic and atmospheric studies, the adjoint method (i.e., 4D-Var) has been more and more widely implemented. However, the ill-posedness is normal if the parameters which are estimated are distributed in space or time domain, especially when the parameters are distributed complexly [36]. In this study the IC and the SS spatially vary and the SS also depends on the time. So the dimension of the estimated model parameters is large. Meanwhile, the observations are just obtained at a limited number of locations. Hence the ill-posedness problem is ineluctable.

An IPS is developed in this paper. In this scheme, some grid points are selected uniformly in the space domain as independent points at which the values of parameters are adjusted and the values at other points are obtained by the interpolation of the values at the independent points. It reduces the number of control variables so as to increase the well-posedness and ensure the continuity of variables which accords with the physical properties.

Let be the value of parameters at the grid point and let be the value of parameters at the independent point . In this paper the is the value of IC or the SS. Then we have the relationship between and :where is the coefficient of linear interpolation, as follows:

is the weighted coefficient of the Cressman form [36]; that is,where is the influence radius and is the distance between the two points.

The gradients of the cost function with respect to the independent parameters are as follows:

Note that the partial derivative term in the right side of (19) is computed by (14) or (15).

With the PM2.5 transport model and corresponding adjoint model constructed, the PM2.5 observations could be assimilated to optimize the parameters with IPS and improve the simulated results by repeating the following steps. (1) Run the PM2.5 transport model with the given parameters. (2) The simulated PM2.5 concentrations could be obtained, and the cost function can be calculated by (7) with the observations. (3) The adjoint model is run and the adjoint variable can be obtained. (4) The gradients of the cost functions with respect to the parameters at the independent points could be calculated by (19). (5) The parameters at the independent points are adjusted by the typical steepest descent method with the gradients. (6) Then the parameters at other points are determined by the interpolation of the values at the independent points by (16).

2.4. Model Settings

The model is applied over a domain encompassing China (see Figure 1). The horizontal resolution for the model is . The integral time step is 600 s and the assimilation window is 168 hours. Considering the actual PM2.5 concentration and the air quality standards in China, the background value is set as 15.0 μg/m3 which is the constant boundary condition at the inflow boundary. The independent points are selected uniformly over each area and the influence radius of Cressman interpolation in IPS is .

The records of PM2.5 concentrations every two hours for 82 major Chinese cities which are shown in Figure 1 are used as observations in the present study. Elbern et al. [26] pointed out that the validity of the assimilation results can be shown satisfactorily by independent observations. So in this study, Jinan (JNa), Zhengzhou (ZhZ), Shenyang (SY), Quanzhou (QZ), Hangzhou (HaZ), Kunming (KM), Chengdu (CD), and Urumchi (UMQ) are selected as “checked cities,” in which the observations are not assimilated and used for testing the assimilated results. And the other cities are seen as “assimilated cities,” in which the observations are assimilated. The wind data are derived from the National Centers for Environmental Prediction (NCEP) winds daily averaged for each latitude by longitude region.

3. Twin Experiments

It is well known that the SS and the IC are significant for the simulation of PM2.5. As the SS is taken as a black box in our model, the inversion of the SS and IC is the vital part of the whole model. In this section, the ideal twin experiments (TEs) are carried out to verify the inversion ability of the adjoint model and test the reasonability and availability of the IPS.

3.1. Setup of Twin Experiments

The process of TE is designed as follows. First, a distribution of the SS (or IC) is prescribed and taken as “true values.” Then the PM2.5 transport model is run using the “true values” and the simulated results recorded at grid points of “checked cities” and “assimilated cities” are taken as the “observations”; that is, only the spatiotemporal position of real observation is used and the value of observation is obtained by running the PM2.5 transport model. Considering that the real observations contain noises, artificial random errors are added to the “observations” and the maximum percentage of random errors is 5%. Having obtained the “observations,” a guess value of SS (or IC) that is 0 in this work (or the background value for IC) is assigned to run the PM2.5 transport model. Then the differences between simulated results and “observations” will drive the adjoint model. With the backward integration of the adjoint model and the calculation of the gradient, the SS (or IC) are optimized. With the above procedures repeated, the SS (or IC) will get closer to the “true values” continually and the differences between the simulated results and the “observations” will be decreased simultaneously. The iteration of optimization will terminate when the number of iteration steps is equal to 300 exactly which is sufficient for this work.

Because there is no fixed trend of SS in real life, in TEs the SS varies only in space and keeps invariant at every time step. Considering the distribution and emission of PM2.5 concentrations in China, six different TEs are designed which are described in Table 1.

3.2. Results of Twin Experiments

For inversion ability evaluation, the mean absolute errors (MAEs) between simulated results and “observations” in “checked cities” and “assimilated cities” and between the prescribed “true values” of SS (or IC) and the inverted values are calculated and shown in Table 2. It is calculated that the MAEs at “assimilated cities” decrease by 98.11% in TE1, 89.52% in TE2, 98.07% in TE3, 95.47% in TE4, 89.40% in TE5, and 95.34% in TE6. And the MAEs in “checked cities” decrease by 97.90% in TE1, 86.90% in TE2, 97.78% in TE3, 93.59% in TE4, 90.11% in TE5, and 92.68% in TE6. It is obvious that the “observations” have been assimilated well in all TEs from the MAEs. Comparing the MAEs in TE1 and TE4, it is found that the MAEs in TE1 are less than those in TE4, stating that the IPS could improve the simulated results. And the conclusion is also found from the other TEs.

The prescribed IC, SS, and the inversion results when they are inverted separately are shown in Figure 2. Although the MAEs in TE4 have shown that the PM2.5 observations are assimilated well and the inverted IC is close to the prescribed one, from Figure 2(e) it is found that the inverted IC is inconsecutive and the characteristic of Figure 2(a) which is the prescribed value does not reappear. However, Figure 2(c) which is the inverted one with IPS reproduces the feature of Figure 2(a). In Figure 2(f), the SS cannot reflect the feature of the prescribed one. However, the inverted SS with IPS which is illustrated in Figure 2(d) has the same character of spatial distribution with Figure 2(b). It is concluded that the IPS can ensure the continuity of the inverted IC and the SS which accords with the prescribed one.

Figure 3 illustrates the prescribed IC, SS, and the inverted ones in TE3 and TE6. It is worth noting that the prescribed initial and SS have the overlapped area in space. From Figure 3(e), it is found that the inverted IC without IPS is discontinuous and has different values in Shanxi and Jiangxi provinces. What is worse, the inverted SS without IPS in Figure 2(f) does not have the same pattern with Figure 2(b) and the values in Liaoning province are greater than the prescribed ones. In fact, the pattern of the inverted SS in TE6 is more close to the pattern of the prescribed IC, rather than the prescribed SS, and it is ascribed to the ill-posedness of the inverse problem. But the inverted IC and SS with IPS reproduce the prescribed ones well. It is shown that the IPS can increase the well-posedness of the inverse problem. In addition, the prescribed IC and SS can be inverted successfully with IPS no matter separately or simultaneously.

4. Practical Experiments

4.1. Setup of Practical Experiments

In order to find the best way to invert the IC and the SS in practice, several practical experiments (PEs) are carried out. In PEs, both the position and the values of observations are used, which is different from TEs. It is noted that in this section the observations during 18 to 24 May 2014 are assimilated to invert the SS and the IC, that is, the PM2.5 concentrations at 18 May 0:00.

In the adjoint model, there are several modes to invert the IC and the SS. To get the best simulation results, ten different modes are compared in the PEs, which are described as follows. In PE1, the SS and the IC are inverted simultaneously with 300 optimization iteration steps using IPS. In PE2, only the SS is inverted with 300 optimization iteration steps using IPS. In PE3, only the IC is inverted with 300 optimization iteration steps using IPS. In PE4, the SS is inverted firstly with 150 optimization iteration steps and the IC is inverted subsequently with 150 optimization iteration steps using IPS. In PE5, the IC is inverted firstly with 150 optimization iteration steps and the SS is inverted subsequently with 150 optimization iteration steps using IPS. PE6–PE10 are corresponding to PE1–PE5 without using IPS.

4.2. Results and Analysis of Practical Experiments

Table 3 gives the error statistics for PEs. It could be found that in PE1 the cost function which is normalized by its initial values at the first iteration step reaches the minimum value and the MAE between simulated values and observations in “assimilated cities” is the smallest of the PEs using IPS. At the same time, in PE6 which is the corresponding one without IPS those evaluation indexes are the smallest of PEs without IPS. It is shown that when the IC and the SS are inverted simultaneously with or without IPS, the effect of the data assimilation is the best.

In addition, the normalized cost functions and the MAE between simulated values and observations in “assimilated cities” in the PEs without IPS are less than those in the corresponding ones with IPS, but the MAE between simulated values and observations in “checked cities” is greater than that in the corresponding ones with IPS. It may be stated that, in the PEs without IPS, the simulated results are optimized just in the “assimilated cities” where the observations are assimilated. To verify this, four cities, including Zhanjiang (ZhaJ), Jingzhou (JZ), Hangzhou (HaZ), and Chengdu (CD), are selected to compare the time-varying PM2.5 concentrations. The four cities are representative, because the MAE between simulated results and observations in Zhanjiang is the smallest of the “assimilated cities” and that in Jingzhou is the largest; Hangzhou and Chengdu are major cities of the “checked cities.” The time-varying PM2.5 concentrations of the four cities in the PEs are shown in Figure 4. The results in Zhanjiang and Jingzhou indicate that all PEs could capture the variation trend of PM2.5 concentrations, except PE3 and PE8. In detail, it is found that the simulated results in the PEs without IPS are more close to the observations than those with IPS. However, in Hangzhou and Chengdu, it is obvious that the simulated results in the PEs with IPS except PE3 could show the variation character of PM2.5 observations well, and the simulated results in the PEs without IPS are bad. The results in the four representative cities agree with the aforementioned inference. In fact, the “assimilated cities” are in a minority, and the cities in which there are no observations may gain more attention. So the simulated results in the PEs with IPS are more reasonable and believable. With the MAEs in Table 3 and Figure 4, it is concluded that the best simulated result is gained in PE1 in which the IC and SS are inverted simultaneously using IPS.

The IC is important for the simulation of PM2.5. To compare the IC obtained from different modes, the observations at the initial time are shown in Figure 5. The PM2.5 concentrations are larger than 120 μg/m3 in Chengdu (CD), Wuhu (WHu), Nanjing (NJ), and Yangzhou (YZ), while the value is less than 30 μg/m3 in the other cities whose names are marked in Figure 5. The ICs obtained from PEs in which the IC is inverted are exhibited in Figure 6. It is obvious that the inverted ICs got in the PEs without IPS are discontinuous in several areas, especially in PE6; there are several extrema at the grid points where the observations are assimilated. So it is indicated that in the PEs without IPS just the observations are assimilated fully and the inverted ICs do not accord with the physical properties. In PE3 and PE5 the PM2.5 concentrations in Tianjin and Shandong province are larger than 140 μg/m3 where the observations are about 90 μg/m3 and the simulated results are different obviously with the observations in Nanjing, Yangzhou, and Wuhu. In PE1, the inverted IC conforms to the spatial distribution characteristics of the observations and the MAE between the inversion results and observations is the minimum. So the inverted IC in PE1 is the best one.

From the analysis above, it is obvious that the simulated results and inverted IC in PE1 are the best. So the inverted SS in PE1 will be investigated. The SS at the current time step will affect the PM2.5 concentrations at the next time step and there should be correlation between them. The correlation coefficients between the PM2.5 observations and the corresponding SS are calculated in all cities. The mean values of PM2.5 concentrations during assimilation window and the value of correlation coefficient in all cities are shown in Figure 7. From Figure 7, the following conclusions can be deduced. Firstly, the mean values of PM2.5 concentrations are larger in northwest China which is covered by desert and in middle China where the emission of PM2.5 is enormous. And the mean values are smaller in southeast China where precipitation is rich. Secondly, the correlation coefficient is positive in most cities. With the setting of the PE1, the convection and diffusion would not play roles in just one time step. Therefore, the SS at the current time step will determine the change of the PM2.5 concentrations at the next time step. It is accordant with the fact and proves that the inverted SS in PE1 is roughly reasonable.

For detailed model performance evaluation in PE1, regression statistics along with two measures of bias, the mean bias (MB) and the normalized mean bias (NMB), and two measures of error, the root mean square error (RMSE) and normalized mean error (NME) [14], are calculated. The mean values of MB and RMSE in PE1 are 0.70 μg/m3 and 20.58 μg/m3, respectively, and those for NMB and NME are 1.26% and 23.08%, respectively. Furthermore, the mean value of all the observations is 55.62 μg/m3 and that of the simulated results is 56.32 μg/m3, and the correlation coefficient between the observations and simulated results is 0.85. In detail, the scatterplot of Figure 8(a) indicates that the model captured a majority (95.4%) of observations with a factor of 2 and underestimated the observations in the high PM2.5 concentration range. Since Tapered Element Oscillating Microbalance analyzers (TEOMs, which are commonly used in the EPA-China’s air quality monitoring network) measurements for PM2.5 should be considered as lower limits because of volatilization of soluble organic carbon species in the drying stages of the measurement [12, 37, 39], the undersimulation is likely to be more severe than this evaluation suggests. In order to investigate the model’s performance over time, the values of mean PM2.5 concentrations, NMB, NME, and correlation coefficients are calculated (spatial averages) and plotted as a time series (Figure 8(b)). It is noted that there are no observations at several hours. As can be seen, the mean values of observations and simulated results are almost equal and have the same time-varying trend. In addition, it is found that the mean PM2.5 concentrations vary analogously every day and the values reach the minimum every afternoon. The NMB values range from −14.25% to 18.41%. The NME values are approximate to 23.08% which is the NME value of all the observations and it proved that the model is steady during the time. Spatially, the model estimates observed PM2.5 well in most areas (74.39%) where the NME value is less than 0.25. Moreover, in 7 cities the absolute value of NMB is larger than 0.25 and only Quanzhou (QZ) and Urumchi (UMQ) are the “checked cities.” It indicates that the model can simulate good results in the cities where the observations are not assimilated through the data assimilation. Therefore, whether in terms of time or space, even on the whole, this model simulates the PM2.5 concentrations well through estimating the IC and the SS simultaneously using IPS.

4.3. Discussion

To further indicate that the adjoint method and the IPS can improve the precision of the simulation of the PM2.5 concentrations, the PM2.5 concentrations over China are simulated during 10 to 16 December 2014. It is in different season and different meteorological conditions with the previous PEs and the other conditions are the same. As before, ten same modes numbered as NPE1 to NPE10 are compared. The error statistics for NPEs are shown in Table 4. It is found that when the IPS is used, the effect of the data assimilation and the simulation results are better. As the PM2.5 concentrations in winter are larger and the differences in concentration are steeper than spring, the MAEs are slightly larger than those of PEs aforementioned. In addition, it is proved that the inverted IC and SS are reasonable and the simulated results are the best in NPE1 using the same way aforementioned, which is not shown in this part.

5. Summary and Discussion

To estimate the parameters and simulate the PM2.5 concentrations over China, a PM2.5 adjoint model is set up in this paper. In this model, the IC and the SS are estimated by assimilating the observations with the adjoint method. In addition, an IPS is used in the estimation to decrease the ill-posedness of the inversion problem. Through TEs, the inversion ability of the adjoint model is verified and the effectiveness of the IPS is shown. In PEs, the IC and the SS are estimated in ten different modes. The simulated results indicate that when they are estimated simultaneously with IPS in PE1, the best results are obtained. In PE1 the inverted IC is better than that of other PEs, and the SS is accordant with the fact. In detail, the PM2.5 concentrations in four representative cities have the same trend with the observations. The simulated results in PE1 show that the model is able to capture a majority (95.4%) of PM2.5 observations within a factor of 2, with NMB and NME of 1.26% and 23.08%, respectively. And whether in time or space domain, even on the whole, the PM2.5 concentrations are simulated well in PE1. To further explain the conclusions in PEs, the NPEs are done in winter. And the same conclusions are obtained.

The PM2.5 transport model in this paper is simple and the complex multiphase chemical pathways are lumped into a single term, SS. It is obvious that the SS is circumscribed because the primary PM2.5 is not obtained from direct emission observations and the chemical mechanism that simulates the formation of secondary PM2.5 is not considered. However, the PM2.5 concentrations over China are simulated well by estimating the IC and the SS using the adjoint method with IPS. It is shown that the adjoint method is a powerful way to estimate the parameters in the models and improve the simulated results by assimilating the observations. The adjoint method has a bright prospect in the atmospheric pollutants simulation and forecast model. In addition, the IPS presented in this paper is an effective method to decrease the ill-posedness of inverse problem and further improve the precision of parameters estimation and simulated results. In further studies, the selection of independent points in space will be optimized and the independent points will be even selected over time.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Partial support for this research was provided by the Natural Science Foundation of Shandong Province of China through Grant ZR2014DM017, the State Ministry of Science and Technology of China through Grant 2013AA122803, the National Natural Science Foundation of China through Grants 41371496 and 41206001, and Beijing Environment Foundation for Young Talents.