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

: 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 2 ) using the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) and a three-dimensional variational (3DVAR) system. The observed hourly surface SO 2 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 2 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 2 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 coefﬁcients increased by 50.0%. The improvements in Southern China were more signiﬁcant 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 coefﬁcients increased by 23.5–53.4%. This SO 2 emission optimization methodology is computationally cost-effective.


Introduction
Sulfur dioxide (SO 2 ) 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 2 emissions come from power plants and industries, accounting for more than 70% of the total SO 2 emissions [4]. It is challenging to predict SO 2 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 2 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 2 to sulfate has been commonly underestimated in current global and regional models, leading to an overestimation of SO 2 concentrations [12][13][14][15]. The conversion of SO 2 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 2 , 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 2 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 2 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 2 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 2 emissions with two satellite instruments (SCIAMACHY and OMI). Koukouli et al. [45] updated SO 2 emissions monthly from 2005 to 2015, assuming a linear relationship between SO 2 emissions and the satellite SO 2 column. The result showed that the SO 2 emissions decreased by 28% from 2010 to 2015. Fioletov et al. [46] developed a new mass balance method, which considered the SO 2 lifetimes, to estimate SO 2 emissions from large SO 2 sources by the OMI SO 2 column. Based on the four-dimensional variational (4DVAR) method and GEOS-Chem adjoint model, several studies developed inverse models to update monthly SO 2 emissions. Wang et al. [47,48] developed an inverse model to update monthly SO 2 emissions by assimilating the OMI SO 2 satellite measurements. Li et al. [49] optimized SO 2 emissions based on OMI SO 2 data and then quantitively evaluated the improvement of SO 2 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 2 emissions. Chen et al. [12] used the ensemble square root filter (EnSRF) system to evaluate SO 2 emission changes in January from 2010 to 2015. The results showed that the SO 2 emissions decreased by 9.9-22.9% in Southern China, while the SO 2 emissions increased by 12.7-72.0% in Northern China and Western China. Dai et al. [50] quantitatively estimated the spatial changes of SO 2  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 2 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 2 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 2 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 2 simulations, using updated emissions, were also verified against observations to show the emission improvement. A discussion and conclusions are presented in Section 4.

Observational Data and a Priori Emission Data
Hourly SO 2 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 −3 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].  [51].
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: 30 December 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: 30 December 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 SO2 emissions in the simulation domain.

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 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 2 emissions in the simulation domain.

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, Remote Sens. 2022, 14, 220 5 of 22 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) 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 2 concentration to improve the emission inventory. The SO 2 observations were assimilated using SO 2 concentrations as a state variable. The incremental cost function for 3DVAR is as follows: here, δx is the N-vector, known as the incremental state variable, which is defined as δx = x − x b , where x b is the forecast or background state of WRF-Chem. B is the N × N-matrix, denoting the background error covariance associated with x b . The M-vector d = y − Hx b is known as the observation innovation vector, where y is an observation vector, and the M × 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 2 processes and can be used to calculate the background error covariance of SO 2 concentrations [62]. In this study, the forecast differences in SO 2 concentrations were used to calculate the BEC.
The horizontal length scale was used to determine the magnitude of SO 2 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 (x1−x) 2 2L 2 s . 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 2 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 2 concentrations decreased with increasing height in the BEC ( Figure 2c). 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 SO2 processes and can be used to calculate the background error covariance of SO2 concentrations [62]. In this study, the forecast differences in SO2 concentrations were used to calculate the BEC.
The horizontal length scale was used to determine the magnitude of SO2 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 ( 1− ) 2 2 2 . 1 and are two points, and is the horizontal length scale. According to Zang et al. [63], when the intersection of the decline curve reaches 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 SO2 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 SO2 concentrations decreased with increasing height in the BEC (Figure 2c).

The Assumptions and Procedure to Optimize the Emissions
As stated in the Introduction, forecast errors associated with SO2 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 SO2 emissions for a single time step is shown in Figure 3. First, the background field at time t0 ( 0 ) was assimilated using 3VDAR with the observation data to obtain the analyzed concentration ( 0 a). Second, 0 was used as the initial field to forecast SO2 at t1 ( 1 ) using the WRF-Chem model. Third, similar  As stated in the Introduction, forecast errors associated with SO 2 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 2 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 2 at t1 (x 1 ) using the WRF-Chem model. Third, similar to the first step, the forecasted SO 2 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 2 concentrations (δx 1 ). Fifth, for each box over a grid area, we assumed that δx 1 is mainly caused by the emission inventory uncertainty at the former time (δE 0 ), that is δx 1 can be expressed as a function of f δE 0 . Some factors contribute to the SO 2 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 2 is small for the long-life-cycle of SO 2 . 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 2 forecast error. Sixth, the emission error (δE 0 ≈ f −1 ( δx 1 ) was calculated by converting the forecast error of the concentration. Finally, the new emission is obtained by adding the emission error (δE 0 ) to the prior emission source (E 0 ).
to the first step, the forecasted SO2 concentrations at t1 ( 1 ), as the background field, were assimilated using 3DVAR with the observation data to obtain the analyzed field ( 1 ). Fourth, the difference between 1 and 1 was calculated to represent the forecast error for the SO2 concentrations ( 1 ). Fifth, for each box over a grid area, we assumed that 1 is mainly caused by the emission inventory uncertainty at the former time ( 0 ), that is 1 can be expressed as a function of ( 0 ). Some factors contribute to the SO2 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 SO2 is small for the longlife-cycle of SO2. 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 SO2 forecast error. Sixth, the emission error ( 0 ≈ −1 ( 1 )) was calculated by converting the forecast error of the concentration. Finally, the new emission is obtained by adding the emission error ( 0 ) to the prior emission source ( 0 ).

Figure 3.
Flowchart of the procedure to optimize SO2 emissions in a single time step from t to t1. The main steps are highlighted below. The blue text means that the steps were produced by using WRF-Chem.

Conversion from Forecast Error to Emission Error
The relationship between the SO2 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) 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 SO2 concentrations is then converted from the emission flux error using the function below ( Observed concentrations at t1 analyses concentrations by using 3DVAR at t1 ( 1 ) analyses concentrations by using 3DVAR at t0 ( 0 ) Observed concentrations at t0 Forecast concentrations at The main steps are highlighted below. The blue text means that the steps were produced by using WRF-Chem.

Conversion from Forecast Error to Emission Error
The relationship between the SO 2 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) δ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 2 concentrations is then converted from the emission flux error using the function below ( Figure 3): step. The emission error can be represented by the forecast error for the SO 2 concentration in the grid box ( Figure 3): in the grid box ( Figure 3): The forecast error of SO2 concentration is defined as = − , where is the background of the model forecast after 1 h and is the analyzed field obtained by assimilating the observation with the background . The − in the assimilation process is also known as the increment field of the SO2 concentration.
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 SO2 increment fields were generated for each hour (00-23 UTC) and then used for emission optimization.

The Operational Consideration for Emission Optimization
To improve the representativeness of at each hour (00-23 UTC), the median value of the SO2 forecast error at each hour in each grid box ( , , ) was obtained from all samples of the SO2 forecast error ( , , , ) over the entire simulation period (the month of October in this study): 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. is the number of simulated days, known as the sample number of that hour, with a maximum of 30 in this study. The median SO2 concentrations' increment over one month of samples was used to improve the representativeness of the hourly SO2 increment field ( , , ) . Once the median increment field is obtained, the updated emission can also be calculated. ( , , ) in Equation (4) is used to replace in Equation (3): SO2 concentrations are closely related to meteorological conditions [65,66]. Wind can promote the horizontal diffusion of SO2 [67,68] and affect the lifetime of SO2. The forecast error of SO 2 concentration is defined as is the background of the model forecast after 1 h and x 1 a is the analyzed field obtained by assimilating the observation with the background in the assimilation process is also known as the increment field of the SO 2 concentration. δ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 2 increment fields were generated for each hour (00-23 UTC) and then used for emission optimization.

The Operational Consideration for Emission Optimization
To improve the representativeness of δx 1 at each hour (00-23 UTC), the median value of the SO 2 forecast error at each hour in each grid box δx 1 (i,j,t) was obtained from all samples of the SO 2 forecast error δx 1 (i,j,t,n) over the entire simulation period (the month of October in this study): 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 2 concentrations' increment over one month of samples was used to improve the representativeness of the hourly SO 2 increment field δx 1 (i,j,t) . Once the median increment field is obtained, the updated emission can also be calculated. δx 1 (i,j,t) in Equation (4) is used to replace δx 1 in Equation (3): SO 2 concentrations are closely related to meteorological conditions [65,66]. Wind can promote the horizontal diffusion of SO 2 [67,68] and affect the lifetime of SO 2 . Fioletov et al. [43] found the correlation between SO 2 emissions and SO 2 concentrations to be 0.93 when the SO 2 life-time was between a few to more than ten hours. In addition, SO 2 is converted effectively into sulfate through aqueous reactions in clouds or fog droplets [14], and precipitation affects the wet removal rate of SO 2 [69]. To reduce the SO 2 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 2 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 −1 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 −4 s −1 were also removed, since SO 2 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 −1 were not removed because in this condition, the SO 2 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 2 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.

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 2 concentrations, and WRF-Chem simulated the "real" SO 2 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 −2 h −1 , 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 −2 h −1 and 24.3 mol km −2 h −1 , respectively, while the correlation coefficient (CORR) increased from 0.2 to 0.9, suggesting that the 3DVAR-CCE method can remarkably improve the SO 2 emission source by assimilating the surface SO 2 concentration.

Experimental Design
Three experiments were performed to evaluate the effect of SO2 emission optimization on SO2 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 −1 in 50% of the area of China within

Experimental Design
Three experiments were performed to evaluate the effect of SO 2 emission optimization on SO 2 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 −1 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 2 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 2 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 2 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).  Figure 6 shows the averaged diurnal variations in the SO 2 concentration increment (δx 1 ), emission (δE 0 in Figure 2) increment, and SO 2 observation concentrations in four cities (Beijing, Shanghai, Chongqing, and Guangzhou) for October 2015. The hourly SO 2 emission increment correlated well with the SO 2 concentration increment, and the CORRs were 0.7, 0.5, 0.8, and 0.7 for Beijing, Shanghai, Chongqing, and Guangzhou, respectively. The SO 2 increment in the Beijing and Chongqing was overall negative during 24 h, implying that the forecast SO 2 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 2 emissions in the Beijing region decreased −0.44 ± 0.11Tg mo −1 from 2010 to 2015. The emission errors during 10-12 UTC were negative in four cities, and the second peak of SO 2 observations weakened. Chen et al. [12] also showed the second peak of SO 2 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 2 emission errors and concentrations were different between different cities, since the implementation of the emission reduction policies was different in different regions [12].

Increments of SO 2 Concentration and Optimized Emissions
were 0.7, 0.5, 0.8, and 0.7 for Beijing, Shanghai, Chongqing, and Guangzhou, respectively. The SO2 increment in the Beijing and Chongqing was overall negative during 24 h, implying that the forecast SO2 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 SO2 emissions in the Beijing region decreased −0.44 ± 0.11Tg mo −1 from 2010 to 2015. The emission errors during 10-12 UTC were negative in four cities, and the second peak of SO2 observations weakened. Chen et al. [12] also showed the second peak of SO2 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 SO2 emission errors and concentrations were different between different cities, since the implementation of the emission reduction policies was different in different regions [12].  Figure 1).

Improvement in Hourly SO2 Simulations
To further investigate the improvement in SO2 simulations over different regions, the modeled (Ctrl, UC_DA, N_EM experiment) diurnal pattern of SO2 concentration in Octo-  Figure 1).

Improvement in Hourly SO 2 Simulations
To further investigate the improvement in SO 2 simulations over different regions, the modeled (Ctrl, UC_DA, N_EM experiment) diurnal pattern of SO 2 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 2 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 2 concentration in the three regions were gradually delayed, indicating the effect of time zone differences. 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 8a). The bias and RMSEs in the UC_DA experiment were smaller than 3 μg m −3 and 20 μg m −3 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 SO2 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.  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 2 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 −3 and 20 µg m −3 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 2 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. 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 SO2 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).  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).
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).  Figure 10 shows the monthly average observed and modeled SO2 concentrations during October 2015. The observed SO2 concentrations showed that the most polluted area was located in Northern China, especially in Shandong, Hebei, and Shanxi, and the observed SO2 concentrations were generally lower than 30 μg m −3 in southern China. Compared with the observations, the SO2 concentrations in the Ctrl experiment were overestimated in Southern China, especially in Shanghai, Hubei, and Chongqing. In addition, the SO2 concentrations were underestimated in Northern China and Western China, especially in Jilin, Gansu, Qinghai, and Xizang. Figure 10c,d shows the monthly average SO2 concentrations for the UC_DA and N_EM experiments, respectively. The SO2 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 SO2 forecast skill by using optimized emissions. The SO2 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 SO2 concentration forecasts, the mean concentration, bias, RMSE, and CORR in the seven regions are compared in Table 3 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 SO2 concentration. The average SO2 concentration of the observations from 449 sites in China was 23.4 μg m −3 during October 2015, while the average SO2 concentration in China during October 2015 was 29.0 μg m −3 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 SO2 from 2010 to 2015. Koukouli et al. [45] found that the SO2 emissions decreased by 28% from 2010 to 2015. Zheng et al. [4] also showed that the SO2 emissions decreased by 62% from 2010 to 2017. The RMSE in the Ctrl experiment was 19.6 μg m −3 , 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. To evaluate the effects of the DA update cycle and the optimization of emissions on the SO 2 concentration forecasts, the mean concentration, bias, RMSE, and CORR in the seven regions are compared in Table 3 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 2 concentration. The average SO 2 concentration of the observations from 449 sites in China was 23.4 µg m −3 during October 2015, while the average SO 2 concentration in China during October 2015 was 29.0 µg m −3 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 2 from 2010 to 2015. Koukouli et al. [45] found that the SO 2 emissions decreased by 28% from 2010 to 2015. Zheng et al. [4] also showed that the SO 2 emissions decreased by 62% from 2010 to 2017. The RMSE in the Ctrl experiment was 19.6 µg m −3 , 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.  In addition, the RMSEs and CORRs in the N_EM experiment were all improved in sevens regions, suggesting the improvement of the SO 2 forecast skill by using the optimized emissions.

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 2 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 2 emission error converted by the SO 2 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 −1 and divergence > 10 −4 s −1 within a 1 h window. For a stable atmosphere with a lower PBLH, the converting estimation will be accurate, since the total SO 2 diffuses within the box.
However, for an unstable atmosphere with a high PBLH, the estimation of SO 2 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 2 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. pared with the control experiment were mainly during 09-23 UTC (17-07* local time). Th PBLH during this period was relatively low (Figure 12), which indicates that the PBL wa 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 general overestimated the concentration during 09-23 UTC, and the N_EM experiment could e fectively 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 th bias is negative in the control experiment, the bias may be lower in the N_EM experimen such as in the Energy Golden Triangle. In addition, the Ctrl experiment underestimate the SO2 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 in creased from 2010 to 2015 in Northeastern China and Xinjiang. Chen et al. [14] also foun that the SO2 emissions increased by 49.4% and 72% for Northeastern China and Xinjian from 2010 to 2015. The experiment for October 2015 can also illustrate the influence of the PBL. In Figure 7, the significant improvements of the SO 2 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 2 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 2 emissions increased by 49

Conclusions
In this study, we developed a 3DVAR-CCE method to optimize and evaluate the SO2 emission inventory using WRF-Chem and a 3DVAR data assimilation system. The hourly observed SO2 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 −2 h −1 was set as the a priori emission. The optimized emission was obtained by using the 273 sites with hourly SO2 "observations". The results showed that the RMSE of the optimized emission decreased from 43.2 9 mol km −2 h −1 to 18.9 mol km −2 h −1 , 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 SO2 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%.

Supplementary Materials:
The following supporting information can be downloaded at: www.mdpi.com/xxx/s1, Figure S1: The location of (a) the real emissions (EM_obs) and (b) the background emission (EM_old) in the model domain. Units are mol km −2 h −1 ; Figure S2

Conclusions
In this study, we developed a 3DVAR-CCE method to optimize and evaluate the SO 2 emission inventory using WRF-Chem and a 3DVAR data assimilation system. The hourly observed SO 2 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 −2 h −1 was set as the a priori emission. The optimized emission was obtained by using the 273 sites with hourly SO 2 "observations". The results showed that the RMSE of the optimized emission decreased from 43.2 9 mol km −2 h −1 to 18.9 mol km −2 h −1 , 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 2 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%.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14010220/s1, Figure S1: The location of (a) the real emissions