Optimization and Evaluation of SO₂ Emissions Based on WRF-Chem and 3DVAR Data Assimilation

Abstract

Emission inventories are important for modeling studies and policy-making, but the traditional “bottom-up” emission inventories are often outdated with a time lag, mainly due to the lack of accurate and timely statistics. In this study, we developed a “top-down” approach to optimize the emission inventory of sulfur dioxide (SO₂) using the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) and a three-dimensional variational (3DVAR) system. The observed hourly surface SO₂ concentrations from the China National Environmental Monitoring Center were assimilated and used to estimate the gridded concentration forecast errors of WRF-Chem. The concentration forecast errors were then converted to the emission errors by assuming a linear response from SO₂ emission to concentration by grids. To eliminate the effects of modelling errors from aspects other than emissions, a strict data-screening process was conducted. Using the Multi-Resolution Emission Inventory for China (MEIC) 2010 as the a priori emission, the emission inventory for October 2015 over Mainland China was optimized. Two forecast experiments were conducted to evaluate the performance of the SO₂ forecast by using the a priori (control experiment) and optimized emissions (optimized emission experiment). The results showed that the forecasts with optimized emissions typically outperformed the forecasts with 2010 a priori emissions in terms of the accuracy of the spatial and temporal distributions. Compared with the control experiment, the bias and root-mean-squared error (RMSE) of the optimized emission experiment decreased by 71.2% and 25.9%, and the correlation coefficients increased by 50.0%. The improvements in Southern China were more significant than those in Northern China. For the Sichuan Basin, Yangtze River Delta, and Pearl River Delta, the bias and RMSEs decreased by 76.4–94.2% and 29.0–45.7%, respectively, and the correlation coefficients increased by 23.5–53.4%. This SO₂ emission optimization methodology is computationally cost-effective.

1. Introduction

Sulfur dioxide (SO₂) is a major air pollutant that contributes to poor air quality. It has a significant impact on the environment and climate [1,2,3]. A major proportion of anthropogenic SO₂ emissions come from power plants and industries, accounting for more than 70% of the total SO₂ emissions [4]. It is challenging to predict SO₂ concentrations using regional air quality models because of many factors. Firstly, the uncertainty of an emission inventory can affect the accuracy of model forecasts [5,6]. Secondly, the concentrations of air pollutants are affected by meteorological conditions (e.g., wind, precipitation), and the uncertainties associated with the physical parameterization scheme could affect the accuracies of SO₂ concentration forecasts [7,8,9]. Thirdly, inaccurate chemistry initial conditions (ICs) can reduce the forecast accuracy. It has been acknowledged that improving the accuracy of ICs by data assimilation (DA) is an effective strategy to improve air pollution forecasting [10,11]. In addition, the conversion of SO₂ to sulfate has been commonly underestimated in current global and regional models, leading to an overestimation of SO₂ concentrations [12,13,14,15]. The conversion of SO₂ to sulfate involves various chemical processes, including gas-phase oxidation, aqueous reactions, and heterogeneous reactions, as well as physical parameters (e.g., relative humidity and cloud fraction). Among the aforementioned factors, emission inventory uncertainty is considered to be the primary factor that influences the forecast accuracy of SO₂, especially for complex air pollution analyses in China, where there are varying, multi-source emissions [16,17,18,19]. Improving the accuracy of the emission inventory in terms of the total amount of emissions and the spatial–temporal distribution would facilitate accurate air quality predictions.

Generally, the “bottom-up” approach is used to estimate SO₂ emissions, which requires statistical information on the activities and emission factors from all possible sources [4,18,20]. It is difficult to obtain an accurate emission estimate associated with the “bottom-up” approach. The accuracy if emission estimation depends on the qualities and uncertainties of the underlying information from the energy activities at the national, provincial, and county levels [21,22,23]. The difference among some global SO₂ emission inventories with the “bottom-up” approach is 42% [24]. In China, the emission estimate is usually updated only once every few years, but the SO₂ emissions have decreased by about 62% on average from 2010 to 2017 [4]. Thus, it is difficult for the emission estimate to reflect the change of emissions in real time [12]. In addition, most of the “bottom-up” emission estimates are annual/monthly total amounts, which need to be spatially and temporally allocated into hourly grid emissions for air quality model application [14,25,26,27,28,29]. The grid emission estimates could be inaccurate due to the spatial allocation, casting uncertainty on the model forecast. For the temporal allocation of emission inventories in air quality models, some methods evenly distribute the total 24 h emission amount by virtue of experience [30,31,32]. The hourly factor has often been obtained from power plants’ reports, human activity characteristics, and previous studies [31,32,33,34,35].

To reduce the uncertainties of “bottom-up” emission inventories, some data assimilation (DA) systems were developed to optimize the emission inventories by using observations as constraints [36,37,38,39,40,41,42,43,44]. Lee et al. [40] applied the mass balance method to update SO₂ emissions with two satellite instruments (SCIAMACHY and OMI). Koukouli et al. [45] updated SO₂ emissions monthly from 2005 to 2015, assuming a linear relationship between SO₂ emissions and the satellite SO₂ column. The result showed that the SO₂ emissions decreased by 28% from 2010 to 2015. Fioletov et al. [46] developed a new mass balance method, which considered the SO₂ lifetimes, to estimate SO₂ emissions from large SO₂ sources by the OMI SO₂ column. Based on the four-dimensional variational (4DVAR) method and GEOS-Chem adjoint model, several studies developed inverse models to update monthly SO₂ emissions. Wang et al. [47,48] developed an inverse model to update monthly SO₂ emissions by assimilating the OMI SO₂ satellite measurements. Li et al. [49] optimized SO₂ emissions based on OMI SO₂ data and then quantitively evaluated the improvement of SO₂ emissions by different models, suggesting the efficacy of the optimized emission estimate. The ensemble Kalman filter (EnKF) method is one of the most popular methods used to optimize SO₂ emissions. Chen et al. [12] used the ensemble square root filter (EnSRF) system to evaluate SO₂ emission changes in January from 2010 to 2015. The results showed that the SO₂ emissions decreased by 9.9–22.9% in Southern China, while the SO₂ emissions increased by 12.7–72.0% in Northern China and Western China. Dai et al. [50] quantitatively estimated the spatial changes of SO₂ emissions in China based on the four-dimensional local ensemble transform Kalman filter method (4D-LETKF). It was shown that the SO₂ emissions decreased by 49.4% in November from 2010 to 2016 by assimilating surface SO₂ observations.

To date, a three-dimensional variational (3DVAR) assimilation algorithm has been commonly used to improve ICs, but has not been used to update emission inventories within the context of improving air quality forecasts. In this study, we attempted to develop a 3DVAR concentrations converted to emissions (3DVAR-CCE) method in order to improve the SO₂ emission inventory based on WRF-Chem and a 3DVAR system. The emission inventory for October 2015 over Mainland China was optimized, and the hourly observed surface SO₂ concentrations were assimilated. This study period was chosen because the meteorological conditions during this period were relatively stable and there was relatively little precipitation.

The remainder of this paper proceeds as follows. The methodology is described in Section 2, including the WRF-Chem and 3DVAR DA system configurations and the newly designed SO₂ emission optimization procedure. The observational data used for evaluating the model and experimental design are also provided. In Section 3, the reliability of the optimized emission estimates is verified by analyzing recent control policies and possible realistic emission changes. The SO₂ simulations, using updated emissions, were also verified against observations to show the emission improvement. A discussion and conclusions are presented in Section 4.

2. Data and Methodology

2.1. Observational Data and a Priori Emission Data

Hourly SO₂ data from the China National Environmental Monitoring Center (CNEMC) (http://www.cnemc.cn, last access: 23 September 2021) were used for assimilation and evaluation. There were altogether 1497 national control measurement sites over China in 2015. Most observational sites are in Central and Eastern China, whereas the sites in the west are relatively sparse. As shown in Figure 1A, only 1048 stations were selected (randomly) to be assimilated, and the data of the remaining 499 were used to verify the improvement of using optimized emissions. To ensure data quality, the values exceeding 650 µg m⁻³ were deemed unrealistic and were not assimilated in the 3DVAR system, following Chen et al. [12]. The dataset is widely used to analyze the air quality in China [12,19,21,23,30,42].

Hourly observational data from national automatic meteorological observation stations in China were used in this study. The ground meteorological observation data were collected from the China Meteorological Data Service Centre (http://data.cma.cn/en, last access: 23 September 2021). This dataset is high-quality and widely used in China [52,53]. There were 1623 observational sites in the year 2015. Meteorological variables including surface temperature, relative humidity, precipitation, and wind speed, and direction were included (Figure 1A). The observational data were rigorously quality controlled to reduce the influence of meteorological data outliers on emission optimization by eliminating outliers (e.g., climate extremes) and verifying the temporal and spatial consistencies of the observational data (e.g., wind speed and precipitation data).

A priori emission data from anthropogenic emissions were obtained from the MEIC, developed by Tsinghua University (http://meicmodel.org/?lang=en, last access: 23 September 2021). The resolution is 0.25° × 0.25°, and the base year was 2010. MEIC is a bottom-up emission inventory covering 31 provinces in Mainland China and including eight major chemical species [23]. The MEIC reports anthropogenic emissions from sources in five sectors (power, industry, residential, transportation, and agriculture). The details of the technology-based approach and source classifications can be found in Zhang et al. [23]. The original emission inventory (0.25° × 0.25°) was pre-processed to match the model grid spacing (27 km). Seven different regions are illustrated to address spatial trends, including the North China Plain, Northeast China, Energy Golden Triangle, Xinjiang, Sichuan Basin, Yangzi River Delta, and Pearl River Delta. Figure 1B shows the spatial distribution of a priori SO₂ emissions in the simulation domain.

2.2. WRF-Chem Forecast Model and the 3DVAR DA System

The Weather Research and Forecasting (WRF) model coupled with online chemistry (WRF-Chem) Version 3.9.1 was used to simulate and predict air pollutant emissions in this study. The WRF-Chem system is a fully coupled online air quality model that considers several physical and chemical processes, such as transport, deposition, emission, chemical transformation, photolysis, and radiation. The domain (Figure 1) is centered at 101.5° E, 37.5° N and covers all of China, with a 27 km horizontal resolution (169 × 211 grids). There are 40 vertical layers extending from the surface to 50 hPa, with a finer resolution near the surface. Initial and meteorological boundary conditions are provided by the 1° × 1° NCEP Global Final Analysis (FNL) data with a 6 h frequency. The WRF-Chem configuration was set using the physics schemes (

Table 1. WRF-Chem model configuration.

Physical or Chemical Process	Option
Microphysics	Lin microphysics scheme [54]
Long-wave radiation	Rapid Radiative Transfer Model [55]
Shortwave radiation	Goddard Space Flight Center shortwave radiation scheme [56]
Boundary layer scheme	Yonsei University [57]
Land surface model	Noah land surface model [58]
Cumulus parameterization	Grell 3-D scheme [59]
Aerosol scheme	Model for Simulating Aerosol Interactions and Chemistry (MOSAIC-4 bin) [60]
Gas scheme	Carbon Bond Mechanism-Z [61]
Initial condition for chemical species	10 d spin-up

) such as the WRF Lin microphysics scheme [54], the Rapid Radiative Transfer Model longwave [55] and Goddard shortwave radiation schemes [56], the Yonsei University (YSU) boundary layer scheme [57], the Noah land surface model [58], and the Grell-3D cumulus parameterization [59]. The Model for Simulating Aerosol Interactions and Chemistry (MOSAIC-4 bin) and the Carbon Bond Mechanism-Z were chosen for the aerosol and gas-phase chemistry scheme, respectively [60,61]. A nudging scheme was used in the WRF-Chem simulation, and the FNL data were assimilated every 6 h. Anthropogenic emissions from the Multi-Resolution Emission Inventory for China (MEIC) in 2010 were used as the a priori emission input.

The framework of the 3DVAR DA system is very similar to that described by Li et al. [62] and Zang et al. [63,64]. However, previous studies focused on the analysis of ICs for PM_2.5 to improve PM_2.5 forecasting, and this study focused on the SO₂ concentration to improve the emission inventory. The SO₂ observations were assimilated using SO₂ concentrations as a state variable. The incremental cost function for 3DVAR is as follows:

$$J\left(\delta x\right)=\frac{1}{2}\delta x^T{B}^{-1}\delta x+\frac{1}{2}{\left(H\delta x-d\right)}^T{R}^{-1}\left(H\delta x-d\right)$$

(1)

here, $\delta x$ is the N-vector, known as the incremental state variable, which is defined as $\delta x=x-x^b$, where $x^b$ is the forecast or background state of WRF-Chem. $B$ is the $N\times N$-matrix, denoting the background error covariance associated with $x^b$. The $M$-vector $d=y-H{x}^b$ is known as the observation innovation vector, where $y$ is an observation vector, and the $M\times M$-matrix $R$ is the observation error covariance. The $H$ is the observation operator that computes the modeled observation estimates from state variables. With the given observations, the background error covariance (BEC) is also important for determining the performance of a 3DVAR scheme to a substantial degree. In this study, we followed the process outlined by Li et al. [62].

The National Meteorological Center (NMC) method has been used for estimating the air pollutant concentrations’ background error covariance [62,63]. Our previous studies, Li et al. [62] and Zang et al. [63], used 1 mo of the 48 h and 24 h forecast differences as the background error and improved forecasts of PM_2.5 concentrations for more than 16 h. The differences between 48 h and 24 h forecasts were generated from 10 October 2015 to November 11 2015. The first initial chemical field at 00 UTC on 10 October 2015 was obtained from a 10 d forecast in consideration of spin-up. The subsequent initial chemical fields were derived from the former forecast one day prior. For a pair of experiments with 24 h and 48 h forecasts, the initial chemical fields were alternatives from different forecasts, resulting in different forecasts for the next cycle. It was assumed that these differences are representative of the short-term forecast errors in transport-related SO₂ processes and can be used to calculate the background error covariance of SO₂ concentrations [62]. In this study, the forecast differences in SO₂ concentrations were used to calculate the BEC.

The horizontal length scale was used to determine the magnitude of SO₂ variance in the horizontal direction. This scale can be estimated by the curve of the horizontal correlation with distances, and the horizontal correlation is approximately expressed by a Gaussian function $e^{\frac{{\left(x1-x\right)}^2}{2L_s^2}}$. $x1$ and $x$ are two points, and $L_s$ is the horizontal length scale. According to Zang et al. [63], when the intersection of the decline curve reaches $e^{1/2}$, the distance can be approximately as the horizontal length scale in Figure 2a. The horizontal length scale was 81 km, which is approximately three-times larger than the scale used in this study. The vertical variance of SO₂ concentrations was considered by vertical correlations in the BEC. A strong relationship was observed in the boundary layer (approximately below the 20th model layer) in the vertical direction (Figure 2b). The standard deviation demonstrates the reliability of the forecasting model, and the standard deviation for the vertical distribution of SO₂ concentrations decreased with increasing height in the BEC (Figure 2c).

2.3. The Methodology Used for Optimizing the SO₂ Emission Inventory

2.3.1. The Assumptions and Procedure to Optimize the Emissions

As stated in the Introduction, forecast errors associated with SO₂ simulation concentrations could be induced by many factors, including uncertainties in emission inventories, physical and chemical mechanisms, and ICs. Thus, three important steps to optimize the emissions are to (1) generate the “best” concentration fields to obtain the forecast error, (2) eliminate errors from aspects other than emissions, and (3) convert the concentration forecast error into an emission error. The foremost goal is to generate the best concentration fields by DA because it is the best initial field for forecasting and the increment field for emission optimization.

A flowchart of the procedure for optimizing SO₂ emissions for a single time step is shown in Figure 3. First, the background field at time t0 ($x^0$) was assimilated using 3VDAR with the observation data to obtain the analyzed concentration ($x^{0a}$). Second, $x^0$ was used as the initial field to forecast SO₂ at t1 ($x^1$) using the WRF-Chem model. Third, similar to the first step, the forecasted SO₂ concentrations at t1 ($x^1$), as the background field, were assimilated using 3DVAR with the observation data to obtain the analyzed field ($x^{1a}$). Fourth, the difference between $x^{1a}$ and $x^1$ was calculated to represent the forecast error for the SO₂ concentrations ($\delta x^1$). Fifth, for each box over a grid area, we assumed that $\delta x^1$ is mainly caused by the emission inventory uncertainty at the former time ($\delta E^0$), that is $\delta x^1$ can be expressed as a function of $f\left(\delta E^0\right)$. Some factors contribute to the SO₂ forecast error, including the initial field, chemical reactions, and meteorological conditions. The contribution of these factors is relatively small compared with that of the uncertainties in the emission inventory for short-time forecasts (e.g., 1 h), because the initial field is accurate after data assimilation and the chemical reaction error for SO₂ is small for the long-life-cycle of SO₂. In particular, under the meteorological conditions of clear sky and stationary winds, dry/wet deposition and the diffusion processes are not expected to affect the SO₂ forecast error. Sixth, the emission error ($\delta E^0\approx f^{-1}(\delta x^1))$ was calculated by converting the forecast error of the concentration. Finally, the new emission is obtained by adding the emission error ($\delta E^0$) to the prior emission source ($E^0$).

2.3.2. Conversion from Forecast Error to Emission Error

The relationship between the SO₂ emission source and the concentrations in a single grid box is shown in Figure 4. It is indicated that the forecast error of the concentrations after 1 h (e.g., 01 UTC) $\delta x^1$ is mainly due to the uncertainty in the emission source 1 h earlier (e.g., 00 UTC) in that grid box. The forecast error of the SO₂ concentrations is then converted from the emission flux error using the function below (Figure 3):

$$\delta x^1\approx f\left(\delta E^0\right)=V_m\times \frac{{\rho }_{air}}{\rho }\times \frac{\Delta S}{\Delta z}\times \delta E^0\times \Delta t$$

(2)

here, $V_m ({22.4\times 10}^{-3} {\mathrm{m}}^{-3})$ is the gas molar volume, $\rho$ is the air density of the actual atmosphere (${\mathrm{kg} \mathrm{m}}^{-3}$), ${\rho }_{air }\left(1.29 \mathrm{kg} {\mathrm{m}}^{-3}\right)$ is the air density in a standard state indicating the gas molar volume, $\Delta S$ is the grid area, $\Delta z$ is the mode layer height, and $\Delta t \left(1 \mathrm{h}\right)$ is the time step. The emission error can be represented by the forecast error for the SO₂ concentration in the grid box (Figure 3):

$$\delta E^0\approx f^{-1}\left(\delta x^1\right)=\frac{1}{V_m}\times \frac{\rho }{{\rho }_{air}}\times \frac{\Delta z}{\Delta S}\times \frac{1}{\Delta t}\times \delta x^1$$

(3)

The forecast error of SO₂ concentration is defined as $\delta x^1=x_a^1-x_b^1$, where $x_b^1$ is the background of the model forecast after 1 h and $x_a^1$ is the analyzed field obtained by assimilating the observation with the background $x_b^1$. The $x_a^1-x_b^1$ in the assimilation process is also known as the increment field of the SO₂ concentration.

$\delta x^1$ is a three-dimensional array, known as the increment field obtained from the background of the model forecast at a 1 h frequency by assimilating the observational data. For each grid, the SO₂ increment fields were generated for each hour (00–23 UTC) and then used for emission optimization.

2.3.3. The Operational Consideration for Emission Optimization

To improve the representativeness of $\delta x^1$ at each hour (00–23 UTC), the median value of the SO₂ forecast error at each hour in each grid box $\delta x_{\left(i,j,t\right)}^1$ was obtained from all samples of the SO₂ forecast error $\delta x_{\left(i,j,t,n\right)}^1$ over the entire simulation period (the month of October in this study):

$$\delta x_{\left(i,j,t\right)}^1=mid\left[\delta x_{\left(i,j,t,n\right)}^1,n\right]$$

(4)

here, i and j are the grid numbers in the east–west and north–south directions, respectively, and t is the hour of a day. $n$ is the number of simulated days, known as the sample number of that hour, with a maximum of 30 in this study. The median SO₂ concentrations’ increment over one month of samples was used to improve the representativeness of the hourly SO₂ increment field $\delta x_{\left(i,j,t\right)}^1$. Once the median increment field is obtained, the updated emission can also be calculated. $\delta x_{\left(i,j,t\right)}^1$ in Equation (4) is used to replace $\delta x^1$ in Equation (3):

$$\delta E_{\left(i,j,t\right)}^0=\frac{1}{V_m}\times \frac{{\rho }_{\left(i,j,t\right)}}{{\rho }_{air}}\times \delta z_{\left(i,j\right)}\times \frac{1}{\delta t}\times mid\left[\delta x_{\left(i,j,t,n\right)}^1,n\right]$$

(5)

SO₂ concentrations are closely related to meteorological conditions [65,66]. Wind can promote the horizontal diffusion of SO₂ [67,68] and affect the lifetime of SO₂. Fioletov et al. [43] found the correlation between SO₂ emissions and SO₂ concentrations to be 0.93 when the SO₂ life-time was between a few to more than ten hours. In addition, SO₂ is converted effectively into sulfate through aqueous reactions in clouds or fog droplets [14], and precipitation affects the wet removal rate of SO₂ [69]. To reduce the SO₂ forecast error caused by horizontal diffusion, chemical reaction, and wet deposition, we chose only the samples that satisfied several strict conditions to optimize the SO₂ emission. Three meteorological variables (precipitation, wind, and divergence) were used as the diagnostic criteria to screen valid samples. The samples with observed hourly precipitation > 0 mm or wind speed > 4 m s⁻¹ were removed to eliminate the effects of wet deposition or heavy diffusion. Since the locations of the observed data are usually not on the model grid, all the samples of the four grids around the observed station were removed. The samples with divergence > 10⁻⁴ s⁻¹ were also removed, since SO₂ concentrations may diffuse out of the grid box under this condition, resulting in an incorrect estimation. Note that the divergence calculation uses the wind speed of the model grid and not the observed wind speed. Samples with divergence < 0 s⁻¹ were not removed because in this condition, the SO₂ concentrations are still in the same box, diffusing upward. The planetary boundary layer height (PBLH) was used to consider the impact of upward diffusion. The SO₂ forecast errors for all layers below the PBLH were converted into emission source errors by Equation (3) and then accumulated as the total emission source error.

2.4. Observing Systems Simulation Experiment

An observing systems simulation experiment (OSSE) was designed to evaluate the effect of the inverse model system. The real emission included 273 sources (Figure S1a). These sources were divided into 13 arrays and 9 columns. The data were randomized and placed on the grid as the real emissions. The hourly emission factors of the real emissions were based on industry and power reports. Then, the real emissions and WRF-Chem model were applied to simulate the 10 d SO₂ concentrations, and WRF-Chem simulated the “real” SO₂ concentrations. The background emission included 273 sources (Figure S1b), which was applied as the a priori emissions. The average value of the background emissions was 50 mol km⁻² h⁻¹, and the hour factors of the background emissions were the same.

The OSSE was performed (Figure 5) with the background emissions as the a priori emissions. The hourly real concentrations, which were simulated by using WRF-Chem model with the real emissions, were assimilated to obtain the optimized emissions. The optimized emissions were obtained according to the steps in Figure 3 and Equation (5). In these 10 d DA cycle simulations, approximately 45% of the samples were removed due to precipitation, wind speed, and divergence. Figure 5 shows the scatter plot between real emissions, a priori emissions, and optimized emissions. It was found that the 3DVAR-CCE method could effectively improve the accuracy of the emission source. Compared with the background emissions, the averaged bias and root-mean-squared error (RMSE) were decreased by 4.4 mol km⁻² h⁻¹ and 24.3 mol km⁻² h⁻¹, respectively, while the correlation coefficient (CORR) increased from 0.2 to 0.9, suggesting that the 3DVAR-CCE method can remarkably improve the SO₂ emission source by assimilating the surface SO₂ concentration.

2.5. Experimental Design

Three experiments were performed to evaluate the effect of SO₂ emission optimization on SO₂ concentration forecasts over China, i.e., control experiment, update cycle DA (UC_DA) experiment, and new emission (N_E) experiment, with the simulation period from 10 October 2015 to 10 November 2015. The meteorological conditions during this period were relatively stable, and there was relatively little precipitation. There were more than 10 d with the average wind speed less than 4 m s⁻¹ in 50% of the area of China within 1 mo. For most areas of China, the rainy days were less than 5 d. The emission source from October 10th to November 10th (October, the same as below) was used for emission optimization due to large variations in emission sources (transportation and fireworks) during 1–7 October, coinciding with China’s National Day. The simulated products were produced at hourly intervals. The FNL data with a 6 h frequency were interpolated to a 1 h frequency by running the WRF Pre-processing System (WPS), and the interpolated data were used to update the meteorological IC and boundary conditions (BCs) to prevent meteorology simulation drifting by nudging.

A control experiment was performed without applying DA, initialized with a former forecast result, starting from 00 UTC on October 1 2015, and the emission inventory was MEIC_2010. The UC_DA experiment was started at 00 UTC on October 10 with a 3DVAR analysis field by assimilating the SO₂ concentrations and restarted every hour with an updated 3DVAR analysis field as the initial chemical field. For UC_DA experiment, not all observation data were assimilated. In a model grid, if there were more than 2 sites, we randomly selected 1 observation site for verification and other for assimilation purposes. In this study, approximately 1048 national control sites were applied to improve the SO₂ ICs, and approximately 449 sites were used to evaluate the model performance (Figure 1A). In the UC_DA experiment, MECI_2010 was used as an a priori emission inventory to simulate the SO₂ concentration in October 2015. Then, the updated optimized emissions for 2015 were obtained. To assess the optimized emissions obtained from the UC_DA experiments, an N_EM experiment, starting from 00 UTC on October 10 for the same period, was conducted and the results were compared with those from the control experiment. The initial chemical conditions and meteorological IC/BC values were consistent with those of the control experiment. The updated optimized emissions were used in the N_EM experiment (

Table 2. Details of different experiments.

Experiment Name	Emissions	Data for DA	Cycle DA
Control (Ctrl)	MEIC_2010	/	/
Update Cycle DA (UC_DA)	MEIC_2010	Surface observations of SO2 concentration	hourly
New Emission (N_EM)	2015 optimized emissions	/	/

).

3. Results

3.1. Increments of SO₂ Concentration and Optimized Emissions

Figure 6 shows the averaged diurnal variations in the SO₂ concentration increment ($\delta x^1$), emission ($\delta E^0$ in Figure 2) increment, and SO₂ observation concentrations in four cities (Beijing, Shanghai, Chongqing, and Guangzhou) for October 2015. The hourly SO₂ emission increment correlated well with the SO₂ concentration increment, and the CORRs were 0.7, 0.5, 0.8, and 0.7 for Beijing, Shanghai, Chongqing, and Guangzhou, respectively. The SO₂ increment in the Beijing and Chongqing was overall negative during 24 h, implying that the forecast SO₂ concentrations was overestimated by WRF-Chem with the a priori emissions. The large overestimation was probably induced by the strict policies of emissions reduction implemented during recent years. Koukouli et al. [45] also found that the SO₂ emissions in the Beijing region decreased −0.44 ± 0.11Tg mo⁻¹ from 2010 to 2015. The emission errors during 10–12 UTC were negative in four cities, and the second peak of SO₂ observations weakened. Chen et al. [12] also showed the second peak of SO₂ emissions in the northern and western regions was much lower than the a priori emissions and was obscure in southern regions. In addition, the trends of hourly SO₂ emission errors and concentrations were different between different cities, since the implementation of the emission reduction policies was different in different regions [12].

3.2. Forecast Performance

3.2.1. Improvement in Hourly SO₂ Simulations

To further investigate the improvement in SO₂ simulations over different regions, the modeled (Ctrl, UC_DA, N_EM experiment) diurnal pattern of SO₂ concentration in October 2015 was compared with the observational data (Figure 7). The observational data showed large regional variations, reflecting the differences in regional energy consumption/emission changes. The diurnal pattern of SO₂ concentration showed a single peak distribution in Southern China and North China Plain. For Northeastern China, the Energy Golden Triangle, and Xinjiang, two peaks were observed at approximately 01/12 UTC, 02/13 UTC, and 03/14 UTC, respectively. The peaks in SO₂ concentration in the three regions were gradually delayed, indicating the effect of time zone differences.

Figure 7 shows that the result of the UC_DA experiment was closest to the observation data. Compared with the Ctrl experiment, the result of the N_EM experiment was significantly improved in China. The largest improvement was at 11 UTC, and the bias decreased from 25.8 μg m⁻³ to 5.7 μg m⁻³. The improvements in the N_EM experiment were different in seven regions. For Southern China (Sichuan Basin, Yangtze River Delta, and Pearl River Delta), the Ctrl experiment significantly overestimated the SO₂ concentrations, especially during 09–23 UTC. Compared with the Ctrl experiment, the results of N_EM experiment were closer to the observation data, and the average bias decreased from 10.6–19.2 μg m⁻³ to 0.6–4.5 μg m⁻³. For North China Plain and the Energy Golden Triangle, the Ctrl experiment underestimated the SO₂ concentrations during 00–08 UTC and overestimated during 09–23 UTC. The results in the N_EM experiment were improved, indicating that the optimized emissions effectively reduced the uncertainty of a priori emissions. For Xinjiang, the improvement of the N_EM experiment was relatively small, because the sites were smallest among the seven regions and the sites were relatively sparse. The moment when the results of the Ctrl and N_EM experiments passed the significance test is marked by the red hollow circle on the X-label to reflect the improvements between the N_EM and Ctrl experiments. It showed that the improvement in the N_EM experiment was significant at most moments, especially in Sichuan Basin.

Figure 8 shows the scatter plots and spatial distributions of both the simulations and observations for the first peaks (at 02 UTC). It showed that the SO₂ concentrations in the UC_DA experiment were similar to the observations (Figure 8a). The bias and RMSEs in the UC_DA experiment were smaller than 3 μg m⁻³ and 20 μg m⁻³ for all hours, respectively. The correlations (CORR) in the UC_DA experiment were higher than 0.6 at all hours except 09–14 UTC. Therefore, the UC_DA experiment was considered as a suitable 3DVAR analysis field (i.e., closest to the observational data) and can be used to verify the model simulations. In Figure 8a, it was found that the bias and CORR in N_EM experiment increased by 35.6% and 164.7% and the RMSE decreased by 27.8%, compared with the Ctrl experiment. The results of the Ctrl experiment overestimated the SO₂ concentrations in Sichuan Basin, Yangtze River Delta, and Pearl River Delta and underestimated in Northeastern China and the Energy Golden Triangle (Figure 8b,c). The positive and negative deviations in the N_EM experiment are both smaller than those in the Ctrl experiment.

Figure 9 is similar to Figure 8a, except at 10 UTC. The CORR in the UC_DA experiment at 10 UTC was worse compared to the other hours. This was mainly because there were more localized characteristics at 10 UTC in most regions. The localized characteristics caused a sharp local gradient of concentrations in the process of DA, which negatively impacted the improvement of the simulation of the UC_DA experiment. It was difficult to describe the local increment field because the BEC was assumed to be static with a large horizontal correlation coefficient of 81 km in the process of 3DVAR. In addition, the overestimation in the control experiment was more significant at 10 UTC than at 02 UTC because of the overestimation of a priori emission at this hour. However, the average bias and RMSE for the N_EM experiments significantly decreased, and the positive biases were smaller in the N_EM than those in the control experiment (Figure 9c,d).

3.2.2. Spatial Distribution of SO₂ Concentrations

Figure 10 shows the monthly average observed and modeled SO₂ concentrations during October 2015. The observed SO₂ concentrations showed that the most polluted area was located in Northern China, especially in Shandong, Hebei, and Shanxi, and the observed SO₂ concentrations were generally lower than 30 μg m⁻³ in southern China. Compared with the observations, the SO₂ concentrations in the Ctrl experiment were overestimated in Southern China, especially in Shanghai, Hubei, and Chongqing. In addition, the SO₂ concentrations were underestimated in Northern China and Western China, especially in Jilin, Gansu, Qinghai, and Xizang. Figure 10c,d shows the monthly average SO₂ concentrations for the UC_DA and N_EM experiments, respectively. The SO₂ concentrations in the UC_DA experiment were close to the observations because of the 3DVAR hourly cycling DA. The simulations of the N_EM experiment performed better than those of the control experiment, owing to the application of optimized emissions. The difference between the Ctrl and N_EM experiments showed the significant improvement of the SO₂ forecast skill by using optimized emissions. The SO₂ concentrations in the N_EM experiment were improved in most regions of China. However, there were few adjustments in Northeastern China and Southwestern China (Yunnan, Guangxi, Tibet, etc.). The main reason for the limited adjustments in the optimized emissions was that there were relatively fewer sites in these areas.

To evaluate the effects of the DA update cycle and the optimization of emissions on the SO₂ concentration forecasts, the mean concentration, bias, RMSE, and CORR in the seven regions are compared in

Table 3. Statistics of the Ctrl, UC_DA, and N_EM experiments in 7 regions (Units: μg m⁻³ for bias and RMSE). NCP: North China Plain; NEC: Northeastern China; EGT: Energy Golden Triangle; XJ: Xinjiang; SB: Sichuan Basin; YRD: Yangtze River Delta; PRD: Pearl River Delta.

	N Data	Mean Concentration				Bias			RMSE			CORR
		Obs	Ctrl	UC_DA	N_EM	Ctrl	UC_DA	N_EM	Ctrl	UC_DA	N_EM	Ctrl	UC_DA	N_EM
China	241,023	23.41	28.95	22.96	23.16	5.54	−0.45	−0.25	33.81	19.56	25.04	0.18	0.51	0.27
NCP	36,866	33.81	38.42	32.19	36.78	4.61	−1.62	2.97	39.10	23.91	32.96	0.24	0.70	0.40
NEC	29,939	24.77	18.12	24.24	19.91	−6.65	−0.53	−4.86	28.67	22.11	25.51	0.17	0.55	0.31
EGT	43,683	30.84	30.09	25.83	27.34	−0.74	−5.01	−3.49	36.42	24.42	30.64	0.14	0.56	0.28
XJ	16,819	11.59	6.99	13.93	8.21	−4.60	2.34	−3.37	17.09	13.25	14.87	0.16	0.37	0.21
SB	58,012	21.79	38.21	21.56	23.64	16.42	−0.23	1.85	41.80	18.49	25.10	0.16	0.42	0.22
YRD	42,964	25.03	35.60	23.82	25.64	10.56	−1.21	0.61	35.22	19.24	25.01	0.25	0.52	0.31
PRD	12,740	16.05	35.23	19.14	20.56	19.17	3.08	4.51	38.39	15.52	21.23	0.13	0.43	0.20

comprehensively. The forecast differences between the control and N_EM experiments reflected the performance and improvement of the optimized emissions. Compared with the control experiment, the N_EM experiment overall decreased the average SO₂ concentration. The average SO₂ concentration of the observations from 449 sites in China was 23.4 μg m⁻³ during October 2015, while the average SO₂ concentration in China during October 2015 was 29.0 μg m⁻³ in the control experiment, being lower than that of the observation concentration. This reduction of concentration implied that there was a significant decrease in emissions of SO₂ from 2010 to 2015. Koukouli et al. [45] found that the SO₂ emissions decreased by 28% from 2010 to 2015. Zheng et al. [4] also showed that the SO₂ emissions decreased by 62% from 2010 to 2017. The RMSE in the Ctrl experiment was 19.6 μg m⁻³, and the CORR was 0.18. Compared with the Ctrl experiment, the RMSE of the N_EM experiment decreased by 25.9%, while the CORR increased by 50.0%, indicating that the optimized emissions could significantly improve the forecast skill.

At the regional scale, compared with the control experiment, the improvements in the N_EM experiment were much more significant in Southern China than those in Northern China, with the bias and RMSE decreasing by 76.4–94.2% and 29.0–45.7%, respectively. For North China Plain, the average SO₂ observation was 33.8 μg m⁻³, while the average SO₂ concentrations in the Ctrl experiment were overestimated. The bias and RMSE in N_EM experiment decreased by 35.6% and 15.7%, and the CORR increased 66.9%. For Northeastern China and Xinjiang, the Ctrl experiment underestimated the SO₂ concentrations, and the biases were −6.7 μg m⁻³ and −4.6 μg m⁻³, respectively. Compared with the Ctrl experiment, the bias in the N_EM experiment increased to −4.9 μg m⁻³ and −3.4 μg m⁻³, reflecting the emissions increased from 2010 to 2015 in Northeastern China and Xinjiang. In addition, the RMSEs and CORRs in the N_EM experiment were all improved in sevens regions, suggesting the improvement of the SO₂ forecast skill by using the optimized emissions.

4. Discussion

The 3DVAR methodology is commonly used to assimilate meteorology and the chemical initial field to improve model forecasts [10,11]. The meteorology forecasts with the assimilated meteorology IC performed better than those without the DA initial field, and the improvement can be maintained for relatively longer time forecasts. The chemical forecasts with the optimization of the chemical IC showed improvements usually in only a few hours. The reason is that emissions are one of the most important driving factors among all the processes in an air quality model, and emissions affect or even dominate the modeled chemical concentrations for long-term forecasts in polluted regions, instead of the initial conditions. Thus, it is important to reduce the uncertainties of emissions to improve the air pollutants’ forecasts. The “top-down” approach is used to optimize emissions, such as 4DVAR and the EnKF method by assimilation observations. 4DVAR could optimize emissions by using the adjoint of a model. The EnKF uses the flow-dependent covariance generated by an ensemble of model outputs to convert concentration observation data into emissions. However, both methods are computationally expensive.

In this study, we developed a “top-down” approach named the 3DVAR-CCE method to optimize and evaluate the SO₂ emission inventory. This 3DVAR-CCE method is computationally cost-effective compared to other complex systems (e.g., EnKF). The basic idea of the method is to estimate the SO₂ emission error converted by the SO₂ forecast concentration error in each grid box. This method significantly relies on the meteorological condition. We assumed the convention within the grid box in the condition of wind speed < 4 m s⁻¹ and divergence > 10⁻⁴ s⁻¹ within a 1 h window. For a stable atmosphere with a lower PBLH, the converting estimation will be accurate, since the total SO₂ diffuses within the box. However, for an unstable atmosphere with a high PBLH, the estimation of SO₂ emissions will be less accurate, since the height of the PBLH in the WRF-Chem model is inaccurate for the unstable PBL, and the upper wind of the PBL is generally stronger, which may cause the SO₂ too diffuse or advect out of the box. This will result in an underestimation of the emissions. Figure 11 shows the OSSE result of the scatter plot among real emissions, a priori emissions, and optimized emissions at (a) 07 UTC and (b) 23 UTC, respectively. The 07 UTC (15 local time) represents the unstable boundary PBL, while the 23 UTC (07 local time) represents the stable PBL. It is obvious that the skill score of optimized emission under a stable PBL is better than that under an unstable PBL.

The experiment for October 2015 can also illustrate the influence of the PBL. In Figure 7, the significant improvements of the SO₂ concentration for the N_EM experiment compared with the control experiment were mainly during 09–23 UTC (17–07* local time). The PBLH during this period was relatively low (Figure 12), which indicates that the PBL was stable. However, the improvements were small during 03–08 UTC (11–16 local time), while the PBL was unstable with a high PBLH. Note that the control experiment generally overestimated the concentration during 09–23 UTC, and the N_EM experiment could effectively decrease the concentration of the overestimation. Therefore, the bias of the N_EM experiment generally reduces in Table 3, depending on the reduction of 09–23 UTC. If the bias is negative in the control experiment, the bias may be lower in the N_EM experiment, such as in the Energy Golden Triangle. In addition, the Ctrl experiment underestimated the SO₂ concentrations for Northeastern China and Xinjiang, and the bias in the N_EM experiment increased by 26.7% and 26.9%, respectively, reflecting that the emissions increased from 2010 to 2015 in Northeastern China and Xinjiang. Chen et al. [14] also found that the SO₂ emissions increased by 49.4% and 72% for Northeastern China and Xinjiang from 2010 to 2015.

5. Conclusions

In this study, we developed a 3DVAR-CCE method to optimize and evaluate the SO₂ emission inventory using WRF-Chem and a 3DVAR data assimilation system. The hourly observed SO₂ concentrations were assimilated by the 3DVAR system to obtain the hourly analysis field. An hourly forecast was performed with the hourly analysis field using WRF-Chem. There were hourly differences between the hourly forecasts and the analysis fields at the same time. The hourly concentration forecast error was used to convert into the emission error, by neglecting the effect of advection and chemistry reaction within an hour.

An OSSE experiment was designed to evaluate the effect of the 3DVAR-CCE method. The average value of 50 mol km⁻² h⁻¹ was set as the a priori emission. The optimized emission was obtained by using the 273 sites with hourly SO₂ “observations”. The results showed that the RMSE of the optimized emission decreased from 43.2 9 mol km⁻² h⁻¹ to 18.9 mol km⁻² h⁻¹, and the CORR increased from 0.2 to 0.9, compared with the a priori emission.

The 3DVAR-CCE method was also used to optimize and evaluate the emissions of October 2015. The a priori emission was from MEIC_2010, which is a bottom-up emission. The optimized emission was obtained by using 1048 hourly observations sites. A control and N_EM experiment were conducted using the a priori emission and optimized emission, respectively, and the independent data were used to evaluate the performance of the two experiments. The average bias and RMSE of the N_EM experiment decreased by 71.2% and 25.9%, and the CORR increased by 50.0%, compared with the control experiment. This indicated that using the optimized emissions can effectively improve the SO₂ forecasting. Significant improvements were found during 09–12 UTC, while the bias and RMSE of the N_EM experiment decreased by 117.4–165.0% and 30.3–58.3%, respectively. For the spatial distribution, the improvements in Southern China were more significant than those in Northern China. For Sichuan Basin, Yangtze River Delta, and Pearl River Delta, the bias and RMSEs decreased by 76.4–94.2% and 29.0–45.7%, respectively, and the CORRs increased by 23.5–53.4%. For Northern China, the bias and RMSEs decreased by 26.9–35.6% and 11.0–15.7%.