Revealing the sulfur dioxide emission reductions in China by assimilating surface observations in WRF-Chem

The anthropogenic emission of sulfur dioxide (SO₂) over China has significantly declined as a consequence of the clean air actions. In this study, we have developed a new emission inversion system based on a four-dimensional local ensemble transform Kalman filter (4D-LETKF) and the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) to dynamically update the SO₂ emission grid by grid over China by assimilating the ground-based hourly SO₂ observations. Sensitivity tests for the assimilation system have been conducted firstly to tune four system parameters: ensemble size, horizontal and temporal localization lengths, and perturbation size. Our results reveal that the same random perturbation factors used throughout the whole model grids with assimilating observations within about 180 km can efficiently optimize the SO₂ emission, whereas the ensemble size has only little effect. The temporal localization by assimilating only the subsequent hourly observations can reveal the diurnal variation of the SO₂ emission, which is better than updating the magnitude of SO₂ emission every 12 h by assimilating all the observations within the 12 h window. The inverted SO₂ emission over China in November 2016 has declined by an average of 49.4 % since 2010, which is well in agreement with the bottom-up estimation of 48.0 %. Larger reductions of SO₂ emission are found over the a priori higher source regions such as the Yangtze River Delta (YRD). The simulated SO₂ surface mass concentrations using two distinguished chemical reaction mechanisms are both much more comparable to the observations with the newly inverted SO₂ emission than those with the a priori emission. These indicate that the newly developed emission inversion system can efficiently update the SO₂ emissions based on the routine surface SO₂ observations. The reduced SO₂ emission induces the sulfate and PM_2.5 surface concentrations to decrease by up to 10 µg m⁻³ over central China.

1 Introduction

China and India are the top two emitters of anthropogenic sulfur dioxide (SO₂) in the world (C. Li et al., 2017). SO₂ is a toxic air pollutant and the precursor of sulfate aerosol which leads to the acidification of the atmosphere and the current heavy haze problem in China (G. Wang et al., 2016; Huang et al., 2014; Yao et al., 2018). Sulfate aerosol can further perturb the radiative energy budget on Earth through directly scattering solar radiation (Goto et al., 2011) and the hydrological cycle by aerosol–cloud interactions (Ramanathan et al., 2001; Sato et al., 2018; Rosenfeld et al., 2019). Sulfate coating on dust leads to a shorter lifetime of dust by increasing the deliquescence of the mixed dust, inducing a great impact on radiative properties and climate modelling (Zhang et al., 2003; Bauer et al., 2007; Fu et al., 2009; Wang et al., 2013; Qi et al., 2013; Penner, 2019). Hydrophilic polluted continental aerosols such as sulfate and mixed dust serve as cloud condensation nuclei (CCN) and thus have a substantial effect on cloud properties and the initiation of precipitation (Rosenfeld et al., 2008). The liquid and ice water paths of dust-contaminated clouds were found to be significantly smaller than those of dust-free conditions over East Asia (Huang et al., 2006a, b). Asian dust altering cloud microphysics and precipitation was revealed by observations and model simulations (Liu et al., 2019a, b, 2020). This, in turn, plays a key role in the climate system. To mitigate climate change and control air quality, the emission control policies, especially for SO₂ implemented by China since 2006, cover all the major source sectors and have become increasingly stringent over time (Zhang et al., 2012). Consequently, the decreasing trends of SO₂ loading over China have been revealed by satellite observations, demonstrating that SO₂ emissions in China declined by 75 % during 2007–2016 (Wang et al., 2018; C. Li et al., 2017). The relative change rate of SO₂ emission in China during 2010–2017 is also estimated to be −62 % using a bottom-up emission inventory (Zheng et al., 2018).

Timely precise emission inventories such as for SO₂ are the primary inputs to models for air quality prediction and mitigation. All atmospheric chemistry and aerosol models rely on their descriptions of the emissions virtually, which are mostly from bottom-up emission inventories. Bottom-up emission inventories are compiled based on indirect information such as activity data and emission factors (Zhang et al., 2009; Kurokawa et al., 2013; Zheng et al., 2018). Due to the uncertainties of the activity rates and emission factors, large discrepancies of global and regional emissions are identified among different emission inventories (Li et al., 2018; Granier et al., 2011). This demonstrates that there is still no consensus on the best estimates for the emissions of atmospheric compounds. Moreover, bottom-up anthropogenic emission inventories often lag several years behind the present and may quickly become outdated (Zheng et al., 2018), leaving the model without up-to-date emission inventories.

The emission inversion approach can feed historical and near-real-time observations into the models, providing a top-down approach to estimate and update the primary emissions of air pollutants in a timely way (Streets et al., 2013). Generally speaking, variational and ensemble data assimilation approaches are the two most widely used methodologies to estimate the emission fluxes of gases (such as NO_x, CO, VOCs) (Tang et al., 2011; Qu et al., 2017; Wu et al., 2020; Miyazaki et al., 2012b; Cheng et al., 2010; Feng et al., 2020a) and/or aerosols (Dai et al., 2019a; Cohen and Wang, 2014; Peng et al., 2017; Yumimoto et al., 2008). NO_x emission changes over China during the COVID-19 epidemic were inferred from surface NO₂ observations based on an ensemble data assimilation approach (Feng et al., 2020b). The emission reductions during the 2015 China Victory Day Parade were successfully detected with an ensemble data assimilation system (Chu et al., 2018). SO₂ emission inventories over China were updated on monthly or seasonal timescales, assuming a linear relationship between SO₂ emissions and satellite-observed SO₂ column amounts (Koukouli et al., 2018; Lee et al., 2011), known as the mass balance approach (Martin, 2003), although the sulfur chemistry, especially in polluted areas as well as in the interactions of clouds, should be non-linear (Goto et al., 2011; Liao et al., 2003). Fioletov et al. (2015) described a new mass balance approach to simultaneously estimate the SO₂ lifetimes and emissions from large SO₂ point sources using satellite measurements. Based on the variational data assimilation approach in the framework of the GEOS-Chem adjoint model, Y. Wang et al. (2016) developed a new sophisticated inverse modelling (IM) method to update monthly anthropogenic SO₂ emissions in a timely way by assimilating the Ozone Monitoring Instrument (OMI) SO₂ satellite measurements. The non-linear full sulfur chemistry and life cycle in the atmosphere were accounted for for the first time to conduct the top-down estimation of the anthropogenic SO₂ emissions from the GEOS-Chem adjoint model (Y. Wang et al., 2016). However, a great limitation to the application of the variational data assimilation approach is the requirement of developing the complicated adjoint model (Henze et al., 2007; Liang et al., 2020). The ensemble data assimilation approach requires neither linearization of the observation operator and nor an adjoint model; therefore it is much more easily implemented and flexible (Evensen, 2003). Additionally, the ensemble data assimilation and the variational data assimilation use the flow-dependent and precalculated model error covariances respectively (Descombes et al., 2015; Zang et al., 2016). Based on the ensemble square root filter (EnSRF) approach (D. Chen et al., 2019), the recent SO₂ emission changes from the year 2010 in China were successfully updated to improve the model forecast skill. An ensemble Kalman filter data assimilation system was developed to simultaneously optimize the chemical initial conditions and emissions including SO₂ with multi-species chemical observations (Peng et al., 2018). The effects of meteorological assimilation on SO₂ emission inversions were also studied recently (Peng et al., 2020).

Retrievals of SO₂ from satellite-based spectrometers are often contaminated by factors such as interference between ozone and SO₂, and there are significant regional differences between different satellite instruments (Fioletov et al., 2013). This subsequently induces the inconsistency of the inversed regional emissions by assimilating different satellite observations (Lee et al., 2011). Meanwhile, satellite observations are usually assimilated on the monthly timescale due to data availability. Compared with satellite observations, the surface SO₂ observations have higher accuracy and temporal frequency. Therefore, the assimilation of intensive direct surface SO₂ observations can provide more spatial–temporal characteristics of emissions (D. Chen et al., 2019). The China National Environmental Monitoring Centre (CNEMC) started to monitor hourly concentrations of PM_2.5 (particulate matter with diameter ≤ 2.5 µm), PM₁₀, SO₂, NO₂, CO and O₃ in 2012, and it had included 1436 monitoring sites from 369 cities by March 2017 (Wu et al., 2018). Those important direct intensive surface SO₂ observations provide a new chance to estimate the more spatial–temporal characteristics of the SO₂ emission in China (D. Chen et al., 2019).

Due to the limited ensemble members, the ensemble Kalman filter (EnKF) generally has a spurious long-distance correlation problem (Houtekamer and Mitchell, 2001; Miyazaki et al., 2012a). Compared with the EnKF, the local ensemble transform Kalman filter (LETKF) can assimilate measurements simultaneously over different model grids in the parallel architecture (Miyoshi et al., 2007; Hunt et al., 2007). Generally speaking, the LETKF computational time is robust with increasing observations, while that of most other ensemble Kalman filters is essentially proportional to the number of observations (Miyoshi et al., 2007). Moreover, the global analysis consists of linear combinations of ensemble members in different regions, which is not confined to the limited ensemble members and provides better results in many cases (Ott et al., 2004). A four-dimensional LETKF (4D-LETKF) was recently developed to assimilate hourly aerosol optical properties observed by satellite, which can avoid frequent switching between the assimilation and the ensemble aerosol forecasting to significantly reduce computational load (Dai et al., 2019b). In the current study, we implement a 4D-LETKF in the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem). Our major objectives are to investigate whether 4D-LETKF together with the intensive CNEMC SO₂ observations can be applied to quantitatively estimate the spatially resolved changes of SO₂ emissions in China and how sensitive the estimated SO₂ emissions are to the system parameters of the 4D-LETKF.

The remainder of the paper is organized as follows. In Sect. 2, the methodology of our emission inversion system is described in detail. Section 3 presents our experimental design and purposes. The emission inversion results and validations are provided in Sect. 4, before concluding in Sect. 5.

2 Methodology

In order to optimize the SO₂ emissions in this study, we need to formally minimize a scalar cost function J in a Bayesian framework (Hunt et al., 2007; Huneeus et al., 2012). J can be formulated as the sum of the departures of a potential gridded SO₂ emissions x and the corresponding simulated SO₂ surface mass concentrations to the a priori SO₂ emissions x^f and the CNEMC-observed surface SO₂ concentrations y^o:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}J\left(x\right)& =\mathrm{1}/\mathrm{2}{\left(x-x^{\mathrm{f}}\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(x-x^{\mathrm{f}}\right)\\ & +\mathrm{1}/\mathrm{2}{\left(H\left(x\right)-y^{\mathrm{o}}\right)}^T{\mathbf{R}}^{-\mathrm{1}}\left(H\left(x\right)-y^{\mathrm{o}}\right)\end{array},\end{array}$$

where H is the observation operator that forwards the SO₂ emissions to the simulated CNEMC measurements; and B and R are the covariance matrix of the error statistics of the a priori SO₂ emissions and CNEMC observations.

2.1 Forward model and observation operator

The relationship between the emission and the surface concentration of the short-lived reactive gas SO₂ is mainly determined by the atmospheric chemical reactions, transport, and deposition. The fully coupled “online” Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) version 4.1.2 (Grell et al., 2005) serves as the forward model to relate the SO₂ emissions to the simulated observations of surface mass concentration in the current study, which can reflect the complex non-linear relationship between atmospheric chemical concentrations and emissions. Our primary aim is to understand how sensitive the estimated SO₂ emissions are to the parameters of the assimilation system, which requires huge computing resources for sensitivity experiments as described later. Therefore, the model is configured with a domain covering most of China as shown in Fig. 1, with a relatively low horizontal resolution of 50 km and 32 vertical levels (Snyder et al., 2015). A state-of-the-art and highly non-linear gas-phase chemical mechanism named the second- generation Regional Acid Deposition Model (RADM2) (Stockwell et al., 1990) coupled with the Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) aerosol model (Chin et al., 2000, 2002) (i.e. chem_opt = 303) is adopted to simulate the atmospheric sulfur cycle. The RRTMG radiation scheme with prognostic aerosols is selected to consider the aerosol direct effect on atmospheric radiation and photolysis calculations (Iacono et al., 2008). The other main selected physics are identical to those of Dai et al. (2019a). The initial and lateral boundary meteorological conditions are from the NCEP Final (FNL) analysis. To reduce the uncertainties associated with the meteorological fields and facilitate a more straightforward comparison of simulations and observations, the predicted wind (u, v), temperature (t), and specific humidity (q) by the WRF dynamical core are also nudged to the NCEP FNL analysis every 6 h (Dai et al., 2018). The meteorological fields in the planetary boundary layer (PBL) are not nudged. The WRF-Chem simulated surface gridded SO₂ volume mixing ratios in the unit of parts per million (ppmv) are firstly converted to micrograms per cubic metre (µg m⁻³) for comparison to the observations (D. Chen et al., 2019) and then linearly interpolated to the CNEMC site locations.

Figure 1 WRF-Chem model computational domain with the topography. The locations of the assimilated and independent verification observation sites of the China National Environmental Monitoring Centre (CNEMC) are shown by the black circles and red squares, respectively. The three magenta boxes mark the North China Plain (NCP), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD) subregions, where relatively dense observation sites are available.

2.2 SO₂ observations and uncertainties of CNEMC

The quality-assured and quality-controlled measurements of hourly SO₂ surface mass concentration from the CNEMC, which is partly purposefully built for assimilation (Wu et al., 2018), are used to minimize the cost function J. There are a total of 1424 sites in November 2016, and those sites span most of central and eastern China and are primarily located in urban and suburban areas (Peng et al., 2017). Due to unresolved emission variations between urban and suburban areas, the model may have large representativeness errors. To overcome the spatial-scale gaps and to produce more representative observations, super-observations are adopted to average all observations located within a model grid cell (Miyazaki et al., 2012a). Altogether, 463 of 7221 model grid cells are covered by the super-observations (Fig. 1). The locations of the super-observations are assumed to be the locations of the covered model grid cells. To independently verify the assimilation results, we further randomly eliminate the super-observations located in 155 of the 463 grid cells to be assimilated. In other words, the assimilated and independent verification observation sites are randomly decided. The observation error covariance matrix R is assumed diagonal. In other words, the observational error covariance is assumed uncorrelated. The observation error of CNEMC is calculated the same as in D. Chen et al. (2019), which contains both the measurement and representativeness errors. In the assimilation data quality control process, an SO₂ observation leading to absolute innovation exceeding 3 times of the prior total spread is considered to be an outlier and discarded. The innovation is calculated as the observation minus the model-simulated ensemble mean observation determined from the first-guess field, and the prior total spread is the square root of the sum of the background ensemble variance and the observational error variance (D. Chen et al., 2019; Rubin et al., 2016).

2.3 4D-LETKF

The 4D-LETKF assimilation approach generalizes a flow-dependent B from ensemble simulation and finds the minimum of the cost function J through the following five formulas (Cheng et al., 2019):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{x}}^{\mathrm{a}}={\stackrel{\mathrm{‾}}{x}}^{\mathrm{f}}+{\mathbf{X}}^{\mathrm{f}}{\stackrel{\mathrm{‾}}{\mathbf{w}}}^{\mathrm{a}}\text{(3)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{w}}}^{\mathrm{a}}={\stackrel{\mathrm{̃}}{P}}^a{\left({\mathbf{Y}}^{\mathrm{f}}\right)}^T{\mathbf{R}}^{-\mathrm{1}}f\left(r\right)\left(y^{\mathrm{o}}-{\stackrel{\mathrm{‾}}{\mathit{y}}}^{\mathrm{f}}\right)\text{(4)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{P}}^{\mathrm{a}}={\left[(k-\mathrm{1})\mathbf{I}/\mathit{\rho }+{\left({\mathbf{Y}}^{\mathrm{f}}\right)}^T{\mathbf{R}}^{-\mathrm{1}}f\left(r\right){\mathbf{Y}}^{\mathrm{f}}\right]}^{-\mathrm{1}}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathbf{X}}^{\mathrm{a}}={\mathbf{X}}^{\mathrm{f}}{W}^{\mathrm{a}}\text{(6)}& {\displaystyle }& {\displaystyle }{W}^{\mathrm{a}}={\left[(k-\mathrm{1}){\stackrel{\mathrm{̃}}{P}}^{\mathrm{a}}\right]}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{x}}^{\mathrm{f}}$ and ${\stackrel{\mathrm{‾}}{x}}^{\mathrm{a}}$ represent the ensemble mean of the first-guess (a priori) and analysis (a posteriori) SO₂ emissions in this study; the ensemble perturbation matrix X is calculated as $x\left(i\right)-\stackrel{\mathrm{‾}}{x}$, {i=1, 2, …, k}, which k represents the ensemble size; the matrix ${\stackrel{\mathrm{‾}}{\mathbf{w}}}^{\mathrm{a}}$ is the Kalman gain, which specifies the increment between the first guess and the analysis; the vector ${\stackrel{\mathrm{‾}}{\mathit{y}}}^{\mathrm{f}}$ represents the first-guess SO₂ surface concentrations averaged over the ensemble members; the matrix Y^f is calculated as $y^{\mathrm{f}}\left(i\right)-{\stackrel{\mathrm{‾}}{\mathit{y}}}^{\mathrm{f}}$, {i=1, 2, …, k}; I represents the identity matrix. The ensemble analyses are calculated as the sum of the ${\stackrel{\mathrm{‾}}{x}}^{\mathrm{a}}$ and each of the columns of X^a, which serves as part of a priori emission information for the next analysis as described later. The multiplicative inflation factor ρ is used to avoid the filter divergence, which is fixed at 1.1 to inflate the analysis covariance the same as in our previous studies (Dai et al., 2019b; Cheng et al., 2019). In our implementation of the 4D-LETKF, the temporal and spatial localizations are achieved by multiplying the R⁻¹ by a factor f(r) as described in Sect. 3, which makes the effect of an observation on the analysis decay smoothly to zero as the time and physical distance of the observation from the analysis grid point increase (Hunt et al., 2007).

Figure 2 Flow chart of the WRF-Chem/4D-LETKF SO₂ emission inversion system by assimilating the SO₂ observations.

As shown in Fig. 2, each assimilation cycle with 4D-LETKF includes two steps: a first guess and a state analysis. In our implementation, the first guess is the WRF-Chem ensemble forecasting for 12 h with hourly model output. The state analysis optimizes the SO₂ emissions in the past 12 h. The advantages of 4D-LETKF used here are threefold: (1) each member of the ensemble WRF-Chem simulations is continuously integrated for 12 h; therefore, this avoids frequent switching between the ensemble WRF-Chem forecasts and the assimilation (Peng et al., 2017; D. Chen et al., 2019). (2) The asynchronous observations can be assimilated to optimize the current state (Hunt et al., 2007; Dai et al., 2019b). (3) The assimilation time window of 12 h could avoid filter convergence and divergence by finite ensemble samples, since more frequent assimilation forces the experiments more closer together, inducing the underestimation of the model spread and overconfidence in the first-guess state estimate (Schutgens et al., 2010; Miyazaki et al., 2012a; Hunt et al., 2007).

2.4 State variable and forecast model for emission

In this study, the state variable to be optimized is the SO₂ emission. A forecast model for emission is required to propagate observation information and determine the first guess for the next assimilation cycle (Miyazaki et al., 2012a). We adopt the same forecast model for SO₂ emission proposed by D. Chen et al. (2019). The forecast model for SO₂ emission weights 75 % and 25 % toward the SO₂ emission ensemble $E_{t_n}^{\mathrm{a}}$ from the previous analysis and the static initial prior ensemble $E_{t_{\mathrm{0}}}$ as in the following formula:

$$\begin{array}{}\text{(7)}& E_{t_{n+\mathrm{1}}}^{\mathrm{f}}=\mathrm{0.75}\times \mathbf{M}{E}_{t_n}^{\mathrm{a}}{\mathbf{M}}^T+\mathrm{0.25}\times E_{t_{\mathrm{0}}},\end{array}$$

where M is the identity matrix. The optimized SO₂ emission ensemble $E_{t_n}^{\mathrm{a}}$ has SO₂ emissions at 12-hourly time slots, which are used to calculate the first-guess SO₂ emission ensemble $E_{t_{n+\mathrm{1}}}^{\mathrm{f}}$ in sequence for the next assimilation cycle. The SO₂ emission inversion depends on the forecast model; therefore, sensitivity experiments for various different emission forecasts are conducted to tune the assimilation system as given in Table 1. The detailed settings of the sensitivity experiments will be described in the next section. As shown in Figs. S1 and S2 in the Supplement, the temporal and spatial distributions of the ensemble spread of the forecast emissions $E_{t_{n+\mathrm{1}}}^{\mathrm{f}}$ are significantly sensitive to the assimilation system parameters. The initial prior ensemble of SO₂ emission is generated by perturbing the freely available MIX Asian inventory S for November 2010 (M. Li et al., 2017). For example, the SO₂ emission for ensemble member i at a given location (x,y) is calculated as $f_i(x,y)S(x,y)$ (Rubin et al., 2016), and the perturbation f_i(x,y), {i=1, 2, …, k}, follows a log-normal distribution in the k-dimensional space. The mean and the variance of the perturbations f(x,y) are equal to 1 and the MIX SO₂ uncertainty (i.e. 35 %) (M. Li et al., 2017). The horizontal perfectly correlated and random uncorrelated perturbations are both created to generate the initial prior ensemble $E_{t_{\mathrm{0}}}$ and the associated first-guess SO₂ emission ensemble $E_{t_{n+\mathrm{1}}}^{\mathrm{f}}$ as described later. The spatial distribution of the ensemble spread of the $E_{t_{\mathrm{0}}}$ with either horizontal perfectly correlated or random uncorrelated perturbations has a similar pattern to the MIX Asian inventory S, which is generally equal to 35 % multiplied by S. In the MIX inventory, anthropogenic emissions are aggregated into five sectors: power, industry, residential, transportation, and agriculture. However, only the combined total emission is used in the model and updated in the analysis. This aims to decrease the degree of freedom in the analysis (Miyazaki et al., 2012a). A total of 10 chemical species including both gaseous and aerosol species are included in the MIX inventory (M. Li et al., 2017). The original monthly MIX anthropogenic emissions with a horizontal resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ are remapped to the model resolution of 50 km. The residential, transportation, and agriculture emissions are allocated in the lowest model layer, whereas the power and industry emissions are allocated in the lowest seven model layers with the vertical profiles of the emission factors from the Model Inter-Comparison Study for Asia (MICS-Asia) phase III (L. Chen et al., 2019). An improved speciation framework for mapping Asian anthropogenic emissions of non-methane volatile organic compounds (NMVOCs) to multiple chemical mechanisms (Li et al., 2014) is adopted to prepare the initial hourly anthropogenic emissions every 12 h with two separated emission files (i.e. io_style_emissions = 2). We do not apply any diurnal variation for the MIX emissions. Therefore, the initial a priori emissions are identical throughout the 24 h. The emissions of aerosol species for WRF-Chem are prepared according to the study of Peng et al. (2017). Notably, only the SO₂ emission is perturbed and optimized by CNEMC SO₂ observations in this study.

Table 1 Experimental design in this study.

The chemical initial conditions (i.e. atmospheric SO₂ concentrations) for the next forward simulation of the WRF-Chem ensemble also need to be updated with the optimal emission ensemble from the previous analysis (Peng et al., 2015; Peters et al., 2005), and this is achieved by recalculation of the WRF-Chem ensemble with the optimized emissions (Fig. 2). In other words, the WRF-Chem ensemble is performed twice in one assimilation cycle. Theoretically, the uncertainties of the forecast SO₂ concentrations by recalculation of the WRF-Chem ensemble are dependent on the optimized emissions. Lower uncertainties of the initial SO₂ conditions for the next assimilation cycle should be found with higher accurate optimized SO₂ emissions, which in turn makes the SO₂ emission inversion more reasonable. Sensitivity experiments for the SO₂ emission inversions as described in the next section are performed to choose the best assimilation system parameters.

3 Experimental design

The effectiveness of 4D-LETKF is highly dependent on having sufficient spread in the ensemble members in order for the observations to impact the first guess (Rubin et al., 2016; Dai et al., 2019b; Hunt et al., 2007). The ensembles represent the uncertainty in the model first guess; therefore, the method for generating the ensemble is an important consideration for an optimal top-down emission inversion. Meanwhile, 4D-LETKF allows a flexible choice of observations to be assimilated for a specific grid point through horizontal, vertical, and temporal observation localizations (Miyoshi et al., 2007; Dai et al., 2019b; Cheng et al., 2019). The observation localization gradually reduces the effect of an observation as the increasing departure from the analysis grid. In this study, the horizontal localization factor is calculated as the Gaussian function (Miyoshi et al., 2007):

$$\begin{array}{}\text{(8)}& f\left(r\right)=\exp \left(-r^{\mathrm{2}}/\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\right)\end{array}$$

where σ is the localization length, and r is defined as the physical distance between the observation and the analysis grid, and we force the localization factor to zero at 3.65 times the localization length (Zhao et al., 2015). In other words, we ignore observations beyond the cut-off distance. The tuneable horizontal and temporal localization lengths are defined in physical distance (km) and time (h), respectively. The vertical localization is not applied for the SO₂ emission inversion in this study. In other words, we trust the vertical profiles of the emission factors from the Model Inter-Comparison Study for Asia (MICS-Asia) phase III (L. Chen et al., 2019).

A correct choice of the assimilation system parameters such as the ensemble size and correlation length is important to improve the data assimilation performance (Miyazaki et al., 2012b). A series of sensitivity experiments are performed to tune the assimilation system as listed in Table 1. A control experiment assuming the same emissions in November 2016 as in November 2010 (i.e. the standard MIX emissions) is conducted as the deterministic simulation to assess the influence of data assimilation. Considering the GOCART aerosol scheme uses a simple representation of the aerosol chemistry for reducing the computational load, we also conduct another deterministic simulation using a more sophisticated aerosol chemical scheme named Model for Simulating Aerosol Interactions and Chemistry (MOSAIC), coupled with the “lumped-structure” Carbon Bond Mechanism (CBMZ) (Zaveri et al., 2008) (i.e. chem_opt = 9) to investigate the effects of different chemistry and aerosol schemes on SO₂ oxidation. The data assimilation experiments are divided into three groups. In the first group, the same random perturbation factor throughout the whole domain emission grids including vertical and temporal spaces per member is applied to the MIX SO₂ emission to generate 10 ensemble members for the WRF-Chem ensemble forward simulation. The spatial correlation coefficients among the initial prior ensemble of SO₂ emissions over every two model grids are equal to 1, and this makes the spatial correlations among the grid points of the forecast emissions also equal to 1. The same random perturbation factor generates a perfect correlation of emission in both the spatial and temporal spaces; however, this should not be seen as overly restrictive (Schutgens et al., 2010). Firstly, the bottom-up SO₂ emission inventories are to a large extent based on the used activity rates and emission factors (Li et al., 2018). Therefore, with the same random perturbation factors, we effectively create an ensemble of inventories derived with different activity rates and emission factors. Secondly, the same emission standards for SO₂ emission mitigation are implemented in China (Zheng et al., 2018), and this causes the SO₂ reductions to be correlated to a certain extent in both spatial and temporal spaces. Thirdly, the analysis is conducted locally in 4D-LETKF, and the analysis at two grids separated by a distance over about 7.3 times the localization length is mostly independent (Schutgens et al., 2010). In this group, the strongest temporal localization is applied to assimilate only the observations within 1 h of the analysis time. In other words, the hourly SO₂ emission is optimized using only the CNEMC SO₂ observation within the subsequent 1 h, making the inverted SO₂ emission variable hour by hour. The difference of the experiments in this group is only the horizontal localization length, which is assumed to be 10, 30, 50, and 100 km respectively. The purpose of the experiments in this group is to investigate the effects of horizontal localization length on SO₂ emission inversion. Based on the results in the first group as described later in Sect. 4, the second group of experiments by fixing the horizontal localization length of 50 km is subsequently performed with 10, 20, and 40 ensemble members to investigate the effects of ensemble size on SO₂ assimilation. In this group, we remove the temporal localization to investigate the effects of temporal localization on the SO₂ emission inversion. In other words, the hourly SO₂ emission is optimized using all the CNEMC SO₂ observations within the 12 h assimilation window, making the inverted SO₂ emission constant within every 12 h. In the third group, the experiments are performed the same as those of the second group except that the ensembles are generated by independently perturbing the emission in horizontal space but dependently in vertical and temporal spaces. These last two groups of experiments are used to investigate the effects of the ensemble size and perturbation factor on SO₂ emission inversion. The sensitivity data assimilation experiments are all performed for 10 d over the period of 00:00 UTC 8 November to 00:00 UTC 18 November 2016. The global model MOZART-4/GEOS-5 provides the initial and lateral boundary conditions used in this study (https://www.acom.ucar.edu/wrf-chem/mozart.shtml, last access: 10 August 2020). Since we do not know the uncertainties of the global model MOZART-4/GEOS-5, the initial and lateral boundary chemical fields are not perturbed in this study. The first 3 d is used as the spin-up of the data assimilation system, and the subsequent simulation results for 1 week are analysed in the next section. Based on the sensitivity tests of the SO₂ emission inversion system, the experiment H50kmT1hE10Ps, which generally performed better than other experiments, is extended to 00:00 UTC on 1 December 2016. This provides a longer period of 20 d to further validate the assimilation system. We also perform a recalculation experiment with the sophisticated CBMZ/MOSAIC scheme and the updated SO₂ emissions to verify the new SO₂ emission and the associated effects of SO₂ emission reduction.

Figure 3 Spatial distributions of the MIX SO₂ emission in November 2010 (a) and the MEIC SO₂ emission in November 2016 (b) in the model's lowest layer. Spatial distributions of the inverted SO₂ emissions in November 2016 in various data assimilation experiments (c–l).

Figure 4 Comparisons of the MIX SO₂ emissions in November 2010 (a) and the inverted SO₂ emission in November 2016 in various data assimilation experiments (b–k) to the MEIC SO₂ emissions in November 2016.

4 Results

4.1 Sensitivity of the inverted SO₂ emission to the assimilation parameters

The spatial distribution of the MIX SO₂ emission in November 2010 at the model lowest layer is shown in Fig. 3a, which serves as the base of the initial a priori SO₂ emission for our experiments in November 2016. The hotspots of anthropogenic SO₂ emission are found over the economically developed areas such as the North China Plain (NCP), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD). The Multi-resolution Emission Inventory for China (MEIC; http://www.meicmodel.org, last access: 15 February 2021) developed by Tsinghua University can provide the updated SO₂ emission for November 2016 (Fig. 3b), which is used as the independent bottom-up SO₂ emission to validate our inverted SO₂ emission. It is apparent that significant negative changes of SO₂ emission are found over the a priori higher source regions such as the NCP, YRD, and PRD between the years 2010 and 2016, which is in agreement with the changes of the column SO₂ concentrations observed by satellites (Wang et al., 2018). As a consequence of the clean air actions (Zheng et al., 2018), the SO₂ emissions over most areas of China showed a systematic decline from the year 2010 to 2016 (Fig. 4a). Can we reveal the reductions of the SO₂ emission by assimilating the CNEMC-observed surface SO₂ concentration?

As shown in Figs. 3c–f and 4b–e, both spatial distribution and magnitude of the inverted SO₂ emission in November 2016 firstly become closer to the independent MEIC ones but get worse subsequently as the horizontal localization length of the assimilation system increases. The inverted SO₂ emissions of each assimilation experiment are obtained by averaging the ones over the ensemble members. The spatial distributions of the mean differences of the MIX and inverted SO₂ emissions minus the MEIC ones are shown in Fig. S3, and the spatial distributions of the mean ratios between the inverted SO₂ emissions and the MIX ones are shown in Fig. S4. The time series of the hourly SO₂ emissions averaged over China of the initial MIX prior, the forecast, and the analysis of the assimilation experiment H50kmT1hE10Ps from 00:00 UTC 8 November to 23:00 UTC 17 November 2016 are also shown in Fig. S5, which illustrates the adjustment of SO₂ emissions with data assimilation. The experiment with the smallest horizontal localization length (i.e. 10 km) only optimizes the SO₂ emission over the specific grids where there are observations to be assimilated. In such a case, the significant reductions of the SO₂ emission over the grids with no observation sites are unable to be revealed, such as Shandong province in the NCP. With a larger localization length, an observation can constrain the emissions in more grids surrounding the observation, and the observation error more gently increases as the distance from the observation location increases (Hunt et al., 2007). It is obvious that the systemic SO₂ emission reductions especially over the Shandong province are detected by increasing the horizontal localization length. However, the perfect correlations of the emission perturbations over the domains with too large a horizontal localization length cause spurious error covariance, causing the more local emission changes to be undetectable. This is demonstrated by the inverted SO₂ emissions with a localization length of 100 km tending to be lower than the independent MEIC ones, with a mean bias of −0.44 mol km⁻² h⁻¹. Generally speaking, the inverted SO₂ emissions with a horizontal localization length of 50 km are best in agreement with the MEIC ones, with a mean bias of −0.15 mol km⁻² h⁻¹ and root mean square error (RMSE) of 5.34 mol km⁻² h⁻¹.

With a horizontal localization length of 50 km, the spatial distribution of the inverted SO₂ emission by removing the temporal localization is shown in Fig. 3g. It is clearly found that the inverted SO₂ emissions over the Shandong province, YRD, and PRD without temporal localization are lower than those with temporal localization, inducing a larger negative bias and RMSE (Fig. 4f). This demonstrates that it is important to reveal the diurnal variations of the SO₂ emission (Wang et al., 2010). The experiment with temporal localization can reveal the hourly variation of the SO₂ emission by assimilating only the subsequent hourly observations, whereas the experiment without temporal localization only adjusts the magnitude of SO₂ emission every 12 h by assimilating all the observations within the 12 h window.

As shown in Figs. 3g–i and 4f–h, there are no significant differences of the horizontal distribution and magnitude of the inverted SO₂ emission between 10, 20, and 40 ensemble members. This indicates that the ensemble size has little effect on the SO₂ emission inversion when randomly correlated perturbing the emissions. The ensemble forecast with 10 members seems feasible to reveal the SO₂ reductions in China, although the inverted emissions have not converged properly. This in turn significantly reduces the required computational resources and time for the forward calculation of the ensemble model, making the dynamical update of air pollutant emissions affordable when assimilating near-real-time observations.

The inverted SO₂ emission with horizontal random uncorrelated perturbations becomes closer to the independent MEIC one as the number of the ensemble members is increased (Figs. 3j–l and 4i–k). However, the performances of the horizontal distribution and magnitude of the inverted SO₂ emission using 40 ensemble members with horizontal random uncorrelated perturbations are clearly even worse than those using 10 ensemble members with horizontal correlated perturbations. This demonstrates that the independent emission perturbations over each model grid tend to underestimate the model spread due to the current limited ensemble members and the cancellation of neighbouring cells (Pagowski and Grell, 2012; Schutgens et al., 2010).

The mean bias and RMSE of SO₂ emission over China using the MIX SO₂ emission in November 2010 for November 2016 are 2.70 and 9.78 mol km⁻² h⁻¹, respectively (Fig. 4a). For the inverted SO₂ emission by data assimilation, the bias and RMSE reduction rates (Miyazaki et al., 2012b) are estimated as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle \frac{\mathrm{2.70}-\mathrm{|}{B}_{\mathrm{DA}}\mathrm{|}}{\mathrm{2.70}}}\times \mathrm{100}\text{(10)}& {\displaystyle }& {\displaystyle \frac{\mathrm{9.78}-\mathrm{|}{\mathrm{RMSE}}_{\mathrm{DA}}\mathrm{|}}{\mathrm{9.78}}}\times \mathrm{100},\end{array}$$

where B_DA and RMSE_DA are the mean bias and RMSE between the inverted SO₂ emission and the MEIC SO₂ emission in November 2016. As shown in Fig. 5, it is found that (1) the inverted SO₂ emission in every assimilation experiment can reduce both the bias and the RMSE; (2) the randomly correlated perturbation factor is superior to the randomly uncorrelated perturbation factor in reducing the bias and RMSE, and it is generally unaffected by the ensemble size; (3) the experiment H50kmT1hE10Ps shows the best performance in reducing both the bias and the RMSE, decreasing the bias and RMSE by 94.5 % and 45.4 % respectively.

Figure 5 Reductions of the bias and root mean square error (RMSE) between the inverted SO₂ emissions in various data assimilation experiments and the MEIC ones referring to those between the MIX and MEIC SO₂ emissions.

Figure 6 Spatial distributions of the mean biases between the simulated surface SO₂ concentrations in various experiments and the CNEMC-observed ones over both the assimilated and independent sites. The locations of the assimilated and independent verification observation sites of the CNEMC are shown by the circles and squares, respectively.

Figure 7 Same as Fig. 6 but for the RMSEs.

4.2 Sensitivity of the surface SO₂ concentration to the emission

Figures 6 and 7 show the horizontal distributions of the biases and RMSEs between the surface SO₂ concentrations simulated in various experiments and the CNEMC-observed ones over both the assimilated and independent sites. The SO₂ concentrations in each assimilation experiment are obtained by averaging the ones over the WRF-Chem ensemble recalculations with the optimized emissions. The spatial distributions of the mean SO₂ concentrations simulated with the original MIX emissions and the updates of the simulated SO₂ concentrations with the inverted SO₂ emissions are shown in Fig. S6. The spatial distributions of the mean differences of the SO₂ concentrations simulated in the FR and FR_CM experiments are also shown in Fig. S6. It is apparent that significant RMSEs and positive biases are found over the a priori SO₂ emission hotspot regions such as the NCP, YRD, and PRD in both the two free-run experiments, whereas slight RMSEs and negative biases are both found over northwestern China. Furthermore, the horizontal distributions of both the biases and RMSEs of the two free-run experiments are generally similar. As given in Table 2, the relative differences of the RMSEs in the FR and FR_CM experiments are both less than 1 % over the assimilated and independent sites, although the mean biases in the FR_CM experiment both tend to be slightly smaller than those in the FR experiment. This demonstrates that the biases and RMSEs between the simulated and observed surface SO₂ concentrations are not induced by the uncertainties of the different chemical reaction mechanisms but are due to the uncertainties of the used SO₂ emissions. The simulated ensemble mean surface SO₂ concentrations by recalculating the WRF-Chem with the inverted SO₂ emissions in all assimilation experiments are shown to be more comparable to the observations, and the performances of the simulated SO₂ surface concentrations are clearly affected by the inputs of the different inverted SO₂ emissions due to assimilation system parameters. This indicates that the uncertainties of the different chemical reaction mechanisms in simulating SO₂ concentrations are much smaller than those of the SO₂ emissions. In the first group of data assimilation experiments, the largest biases and RMSEs of the simulated and observed SO₂ surface concentrations over both the assimilated and independent sites are found in the H10kmT1hE10Ps experiment. This indicates that the SO₂ emission changes existing grid correlations, and the SO₂ emission inversions over only the grids with available assimilated sites are not sufficient to reveal the real SO₂ emission changes in the grids without observation sites. In addition, the largest biases and RMSEs over both the assimilated and independent sites are still found in the third group of data assimilation experiments, although the biases and RMSEs decrease with the increase of the number of ensemble members. This further illustrates there are correlations of the gridded SO₂ emission changes, and the random perfectly correlated emission perturbation factors over the model grids are superior to the random uncorrelated emission perturbations for current emission inversions. The latter is probably due to the currently limited number of ensemble members for reducing the computational resources. However, the sophisticated random uncorrelated emission perturbations should have better performances with large or unlimited ensemble members. Similar to the inverted emissions, the experiments in the second group show the ensemble size has little effect on the biases and RMSEs of the SO₂ surface concentrations over both the assimilated and independent sites when the ensemble members are generated by perturbing the emissions perfectly correlated over the domain grids. The reductions of the biases of the SO₂ surface concentrations in both the assimilated and independent sites benefit from the temporal localization, although the RMSEs are slightly increased. It is interesting that the smallest RMSE of the SO₂ surface concentrations over the independent sites is also found in the H50kmT1hE10Ps experiment, with a value of 36.20, which shows that the inverted SO₂ emissions are also best in agreement with the independent MEIC ones. This further indicates that the assimilation system parameters used in this experiment are suitable for the SO₂ emission inversion, decreasing the biases of SO₂ surface concentrations over assimilated and independent sites by 87.2 % and 88.9 % respectively. The underestimation of the surface SO₂ concentration with the original MIX emission over northwestern China such as the Gansu province is potentially attributable to the increasing SO₂ emissions due to energy industry expansion and relocation over northwestern China (Ling et al., 2017). The SO₂ emissions and surface concentrations over the Gansu province are increased to reduce the negative biases in the assimilation experiments as shown in Figs. S4 and S6, indicating that our emission inversion system also works well when the prior emissions are underestimated. However, the simulated surface SO₂ concentrations with the inverted emissions are still underestimated over the Gansu province. The reason for the underestimation is twofold: (1) there are limited observations to be assimilated over northwestern China because the observation sites are sparse; (2) the initial a priori MIX SO₂ emission over northwestern China is small and underestimated, causing the model uncertainty to be small relative to the observation one. This translates to a reduced impact of the observation on the a priori emission.

Table 2 The mean biases and root mean square errors (RMSEs) of the simulated SO₂ surface concentrations in various experiments and the CNEMC-observed ones over all assimilated and independent sites.

Figure 8 Frequency distributions of the deviations of the simulated SO₂ surface concentrations in various experiments minus the observed ones.

Figure 8 illustrates the frequency distributions of the deviations of the simulated SO₂ surface concentrations in various experiments minus the observed ones. It is expected that the distributions of the SO₂ surface concentrations deviations for the two free-run experiments in China and the three subregions are all positively biased due to the known overestimation of the SO₂ emissions. The distributions of the SO₂ surface concentration deviations with the updated SO₂ emissions in all the data assimilation experiments show reduced biases over both the assimilated and independent sites. However, the distributions of the deviations with the updated SO₂ emissions in the third group of experiments and the H10kmT1hE10Ps experiment are still positively biased, whereas slightly negative biased distributions are found in the second group of experiments and the H100kmT1hE10Ps experiment. The distributions of the SO₂ concentration deviations with the updated SO₂ emissions in the H50kmT1hE10Ps experiment, as expected, show the best performance, with a peak closer to 0 in both the assimilated and independent sites.

Figure 9 SO₂ emission reductions in November over the period 2010 to 2016 in China and the three subregions estimated by the bottom-up and top-down approaches.

4.3 SO₂ reduction in China and associated effects

Based on the aforementioned sensitivity tests of the SO₂ emission inversion system, the experiment H50kmT1hE10Ps is extended to 00:00 UTC on 1 December 2016. This provides a longer period for 20 d to estimate the reduction of the SO₂ emission over China in November over the period 2010–2016. The bottom-up and top-down estimations of the SO₂ emission reduction from 2010 to 2016 are calculated by comparing the MEIC and inverted SO₂ emissions by data assimilation in November 2016 to the MIX SO₂ emission in November 2010. As shown in Fig. 9, the top-down estimation of the SO₂ emission reduction over China is 49.4 %, which is well in agreement with the bottom-up estimation of 48.0 %. In addition, larger SO₂ emission reductions over the three subregions estimated by the bottom-up approach are correctly revealed by the emission inversion system. The top-down and bottom-up estimations of the SO₂ emission reduction over NCP are generally comparable, with values of 56.0% and 52.4 % respectively. The largest SO₂ emission reductions both with the top-down and bottom-up approaches are found over the YRD region, with values of 73.1 % and 61.8 % respectively. The SO₂ emission reduction using the top-down approach is 10 % higher than that using the bottom-up approach over the PRD region. To validate the inverted SO₂ emissions and explore the possible reasons of the overestimation of SO₂ emission reduction over the YRD and PRD using the top-down approach, the time series of the simulated SO₂ surface concentrations in various experiments and the observed ones are shown in Fig. 10. The simulated SO₂ surface concentrations, especially over the YRD subregion in the two free-run experiments, show significant positive biases over all the period, revealing the drawback of the prescribed SO₂ emissions in November 2016 being the same as in November 2010. The simulated SO₂ surface concentrations with the inverted SO₂ emissions using both the RADM2/GOCART and CBMZ/MOSAIC chemical reaction mechanisms are much closer to the observations in both the assimilated and independent sites over all the period. This demonstrates that the WRF-Chem/4D-LETKF emission inversion system can continuously and dynamically update the SO₂ emissions by assimilating the newly available observations as shown Fig. S5. The SO₂ surface concentrations simulated by the FR_CM experiment are sometime lower than those in the FR experiment, especially over the YRD and PRD subregions, indicating the overestimations of the SO₂ emission reduction using the top-down approach over the YRD and PRD are probably due to the simple aerosol chemistry schemes used in RADM2/GOCART (Chin et al., 2000). This is proved as the simulated SO₂ surface concentrations in YRD with the RADM2/GOCART scheme and the inverted SO₂ emissions over the period 18–22 November 2016 are generally comparable to the observed ones, whereas the simulated SO₂ surface concentrations with the sophisticated CBMZ/MOSAIC scheme and the inverted SO₂ emissions are lower than the observed ones. The simulated SO₂ surface concentrations at all sites with the inverted emission in both the FR_CM and assimilation recalculation are generally underestimated. This is due to the inverted emission being sufficient to reduce the overestimations of SO₂ concentration over the a priori SO₂ emission hotspot regions but insufficient to eliminate the underestimations over northwestern China.

Figure 10 Time series of the simulated SO₂ surface concentrations in various experiments and the CNEMC observations.

Based on the inverted SO₂ emissions from 11 November to 1 December 2016, the daily and diurnal variations of the SO₂ emission reductions over China and the NCP subregion are estimated as shown in Fig. 11a and b respectively, and the diurnal variations of the inverted SO₂ emissions over China and the NCP subregion are also shown in Fig. 11c. Generally speaking, the daily variation of the SO₂ emission reduction over China is not so significant. Larger SO₂ emission reductions over the period 17 to 19 November in the NCP induced by the first orange alert for heavy winter air pollution in 2016 are clearly detected by the inverted emissions (Shi et al., 2019). Lower SO₂ emission reductions over China and NCP from 21 to 22 November are probably contaminated by the strong cold wind from the northwestern direction, inducing the low SO₂ concentrations and underestimating the associated ensemble spread. The latter causes the inverted emission to be overconfident in the background emission (Hunt et al., 2007). Since the emissions are constant over time in the a priori MIX inventory, the diurnal variations of the SO₂ emission reduction over China and NCP both reveal higher emission reductions in the night-time, causing the SO₂ emissions in the night-time to be lower than those in the daytime (Fig. 11c). This is generally reasonable as fewer human and economic activities happen in the night-time (L. Chen et al., 2019).

Figure 11 Daily (a) and diurnal (b) variations of the SO₂ emission reductions over China and the NCP subregion based on the inverted emissions. Diurnal variations of the inverted SO₂ emissions over China and the NCP subregion (c).

Figure 12 shows the spatial distributions of the averaged surface concentrations of the sulfate, ammonium, nitrate, and PM_2.5 over 11 November to 1 December 2016 simulated with the CBMZ/MOSAIC mechanism and the original MIX emissions and the absolute and relative changes of the associated aerosol surface concentrations with the newly inverted emissions by data assimilation. It is found that the SO₂ emission reductions cause the sulfate surface concentrations to be reduced by up to 10 µg m⁻³ (50 % over central China), and this is due to the sulfate aerosols being dominated by the occurrence of in-cloud oxidation (Chin et al., 2000; Goto et al., 2015) and more clouds being found over central China (Li et al., 2015; Ma et al., 2014). The nitrate surface concentrations are found slightly increased in central China with the reduction of sulfate aerosols, and this is due to the emissions of the nitrate precursors (i.e. NO and NO₂) not being updated in this study and because NH₄NO₃ is formed only in sulfate-poor aerosols (Zaveri et al., 2008; Chen et al., 2016). The synergy effects of sulfate–nitrate–ammonium induce slight reductions of ammonium surface concentrations, decreasing the PM_2.5 surface concentrations by about 10 µg m⁻³ (10 %) over central China.

Figure 12 Spatial distributions of the averaged surface concentrations of the sulfate, ammonium, nitrate, and PM_2.5 over 11 November to 1 December 2016, simulated with the CBMZ/MOSAIC mechanism and the MIX emissions in November 2010, and the absolute and relative changes of the associated aerosol surface concentrations with the updated emissions by data assimilation.

5 Conclusions

Timely precise emission inventories are crucial to air quality prediction and mitigation. To dynamically update the emissions of air pollutants, we have developed a new emission inversion system based on the 4D-LETKF and the fully coupled model named WRF-Chem. Our emission inversion system considers the complex non-linear relationship between atmospheric chemical concentrations and emissions by ensemble forecasting with perturbed emissions. The emission inversion system is examined to update the outdated MIX SO₂ emissions in November 2010 by assimilating the quality-assured and quality-controlled observations of SO₂ surface concentration from the CNEMC in November 2016. The inverted SO₂ emissions over China by data assimilation for November 2016 are validated with the independent MEIC emissions in November 2016.

Sensitivity tests for the emission inversion system demonstrate that the assumption of the covariance error matrix of the a priori SO₂ emissions has the largest effect on the inverted emissions. The random perfectly correlated emission perturbations throughout the whole model grids with a horizontal localization length of 50 km can best reproduce the independent MEIC SO₂ emissions, decreasing the MIX emission bias and RMSE by 94.5 % and 45.4 % respectively. The independent emission perturbations over each model grid tend to underestimate the model spread due to the currently limited number of ensemble members and the cancellation of neighbouring cells. With the random perfectly correlated emission perturbations, the ensemble size only has little effect on the inverted SO₂ emissions, and the ensemble forecast with 10 members seems feasible to reveal the SO₂ reductions in China. The temporal localization by assimilating only the subsequent hourly observations can reveal the diurnal variation of the SO₂ emission, which is better than updating the magnitude of SO₂ emission every 12 h by assimilating all the observations within the 12 h window.

The known overestimates of the prescribed SO₂ emissions in November 2016 being the same as in November 2010 are successfully detected as the simulated SO₂ surface concentrations, especially over the SO₂ emission hotspot subregions, with two distinguished chemical reaction mechanisms both being significantly positive biased. The simulated SO₂ surface concentrations with the inverted SO₂ emissions in all assimilation experiments are shown to be more comparable to the observations, and the performances of the simulated SO₂ surface concentrations are clearly affected by the inputs of the different inverted SO₂ emissions due to assimilation system parameters. This indicates that the uncertainties of the different chemical reaction mechanisms in simulating SO₂ concentrations are much smaller than those of the SO₂ emissions. The smallest RMSE of the simulated and observed SO₂ surface concentrations over the independent verification sites is also found in the experiment in which the inverted SO₂ emissions are best in agreement with the independent MEIC ones, decreasing the biases of SO₂ surface concentrations by 88.9 %.

The SO₂ emission reduction over China in November over the period 2010 to 2016 is estimated as 49.4 % by assimilating the observations of surface SO₂ concentrations, which is well in agreement with the bottom-up estimation of 48.0 %. In addition, larger SO₂ emission reductions over the NCP, YRD, and PRD estimated by the bottom-up approach are correctly revealed by the emission inversion system. The largest SO₂ emission reductions both with the top-down and bottom-up approaches are found over the YRD region, with values of 73.1 % and 61.8 % respectively, and the simple parameterizations of the aerosol chemistry in the GOCART scheme may cause the overestimates of the SO₂ emission reductions by about 10 %. The SO₂ emission reductions cause the sulfate and PM_2.5 surface concentrations to decrease by up to 10 µg m⁻³ over central China.