CO Emissions Inferred From Surface CO Observations Over China in December 2013 and 2017

China has implemented active clean air policies in recent years, and the spatiotemporal patterns of major pollutant emissions have changed substantially. In this study, we construct a regional air pollution data assimilation system based on the WRF/CMAQ model and ensemble Kalman filter algorithm to quantitatively optimize gridded CO emissions using hourly surface CO measurements over China. The Multi‐resolution Emission Inventory of China CO emission inventories in December 2012 and 2016 are treated as prior emissions, and the CO emissions in December 2013 and 2017 are optimized using the CO observations of December of 2013 and 2017, respectively. The results show that in both periods, assimilation of CO observations significantly improves the CO simulations and emission estimates. Assimilation increases the CO emissions in most areas of mainland China, especially in northern China, and the spatial patterns of the increases in the two periods are similar. Overall, the posterior CO emissions in December 2017 are 17% lower than those in December 2013. Large emission decreases are mainly found in most key urban areas and developed regions, and emission increases are mainly located in their surrounding areas and certain central and western regions, which might reflect the emission migration from developed regions or urban areas to developing regions or surrounding areas. These changes are not found in the prior emissions but are basically consistent with the emission control strategies and industrial transformation and upgrade phenomena in recent years, indicating that our CO assimilation system could successfully capture the temporal and spatial variations.


Introduction
Anthropogenic air pollutant emissions have been increasing rapidly in China and contribute 18-35% of the global air pollutant emissions (Hoesly et al., 2018). In response to severe haze pollution, the Chinese government updated the Ambient Air Quality Standard (GB3095-2012) (Ministry of Environmental Protection of China, 2012) according to the World Health Organization air quality guidelines and launched the Twelfth Five-Year Plan on Air Pollution Prevention and Control in Key Regions. To attain this air quality standard, the State Council of China issued the Air Pollution Prevention and Control Action Plan (Action Plan, 2013 in 2013. Since then, China began to strengthen its emission standards and implement comprehensive measures to reduce air pollutant emissions, including phasing out the outdated industrial capacity, transferring polluted factories, and replacing residential coal with electricity and natural gas. The effective implementation of these measures finally decreased PM 2.5 concentrations by 27.7-39.6% from 2013-2017 over Beijing-Tianjin-Hebei (BTH), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD) (China, 2018). Space-and ground-based observations have confirmed the improvement of the air quality in China (Krotkov et al., 2016;Zhang et al., 2017;B. Zhao et al., 2017).
Carbon monoxide (CO) is one of the main air pollutants and is a product of incomplete combustion and a by-product of the oxidation of hydrocarbons. The reaction of CO and the OH radical not only affects the oxidizing power of the troposphere but also yields ozone (O 3 ) and carbon dioxide (CO 2 ), which results in indirect positive radiative forcing and climate change (Intergovernmental Panel on Climate Change (IPCC), 2013). Therefore, accurate and precise estimates of CO emissions are important. However, the estimates of the CO emissions of China still subject to large uncertainties. Saikawa et al. (2017) compared five emission inventories of anthropogenic air pollutants in China and revealed that the CO emissions have the largest differences among these inventories. Moreover, the uncertainty of emissions at the regional/local scale is even larger than the uncertainty of the total emissions across the country (Andres et al., 2012).
Data assimilation (DA) is a statistical optimization technique that can improve emission estimates by combining observations and background fields (prior emission inventory). Because of the large uncertainties in CO emission inventories, many efforts have been focused on constraining CO emissions by assimilating CO observations, and a series of assimilation studies have demonstrated the feasibility of improving CO emission estimates (Hooghiemstra et al., 2011;Jones et al., 2007;Koohkan & Bocquet, 2012;Yumimoto & Uno, 2006). For example, by assimilating satellite retrieval data, Y. Y. Yin et al. (2015) constrained global CO emissions for the period of , and Z. Jiang et al. (2017 estimated global CO emissions and investigated the dominant reasons for the trends from 2000-2015 using MOPITT retrievals. However, due to the bias drift in MOPITT lower-tropospheric retrievals, the emission estimates in these two studies show opposite trends in India and Southeast Asia. B. Zheng et al. (2018) also used MOPITT retrievals to evaluate the 2005-2016 trends of CO emissions over East Asia and analyzed the underlying drivers of CO changes. To provide extensive global constraints on CO sources, Kopacz et al. (2010) performed a global inversion from May 2004 to April 2005 as constrained by three different satellites (MOPITT, SCIAMACHY, and AIRS). However, there are still challenges in getting consistent inversion results because of the differences across instruments in sensitivity, retrieval techniques, and spatiotemporal sampling. Large differences between the posterior CO emission estimates were also reported by Jones et al. (2009) when using TES and MOPITT satellite retrievals. Most of these satellite data are assimilated at the global scale, and the results of specific regions are then analyzed. However, fine-scale variability (e.g., cities) is difficult to resolve, which may result in large aggregation or opposite errors (Kopacz et al., 2009). Moreover, biases in satellite retrieval can greatly affect the magnitude of estimated emissions (Miyazaki et al., 2015(Miyazaki et al., , 2012. Retrieval-based satellite observations have large differences (Deeter et al., 2015;Warner et al., 2007;, which result in large uncertainties in inversions, whether computed with different satellite retrievals (Silver et al., 2016) or different products (e.g., column data and profile data) of the same satellite (Z. Jiang et al., 2017). Data missing in many regions because of cloud cover, severe air pollution, and so forth are also one of the drawbacks of satellite retrievals. In addition, in order to match with the satellite retrievals (e.g., aerosol optical depth or weighted column average concentration), the simulations should be converted using a radiative transfer model (Z. Liu et al., 2011) or the same weighted column average method (Miyazaki et al., 2015) during constructing an observation operator, which may bring new uncertainties and thereby further increase the inversion uncertainties.
Ground-based observations are also widely used in emission optimization. Yumimoto and Uno (2006) assimilated three surface sites' data to optimize CO emissions in April 2001 in eastern Asia. Hooghiemstra et al. (2011) optimized global CO emission estimates for 2003 and 2004 using 31 surface observations. Although assimilated results showed better agreement with observations, with a limited amount of observations, the posterior emissions largely depended on the prior emissions and error settings. Ground-based observations are more often applied to inversions at the regional scale. For example, Koohkan and Bocquet (2012) estimated CO emissions for 2005 across France using 89 surface sites' data. In China, which has special weather conditions and much higher CO emissions levels, DA studies for CO emission estimates in China are still insufficient. As a pioneering study, X. Tang et al. (2013) assimilated surface CO observations from 25 sites to estimate local CO emissions across BTH with the ensemble Kalman filter (EnKF) method to improve short-term simulations. China has implemented unprecedented measures to reduce emissions. However, there has been no attempt to estimate the change characteristics of CO emissions using the DA method targeting the mainland China, especially after the emission reduction policy implemented in 2013. Since 2013, the Ministry of Ecology and Environment of China has gradually established nationwide, unified calibrated, and well-maintained air pollution monitoring stations. Through present day, hourly CO concentrations have been measured and released nationwide. This widespread coverage of the observation network provides suitable conditions for a CO DA study in China, and it is sufficient to constrain the spatial and temporal scales dictated by the inhomogeneity of emissions and relatively short atmospheric residence time .
In this study, an EnKF-based CO DA system is constructed and used to investigate the characteristics of CO emissions using nationwide surface observations. The advantage of EnKF lies in its easy implementation for complicated systems. It does not require the development of an adjoint model, which is technically difficult and cumbersome for complex chemical transport models, but which is necessary for 4D-VAR, another method commonly used in inversions. The objective of this work is to provide updated CO emission estimates, and the effects of the Action Plan on spatiotemporal patterns of CO emissions spanning 2013 and 2017 are emphatically analyzed. Model errors of transport and chemistry notably influence CO emission inversions (Z. Jiang, Jones, et al., 2013;. Instead of optimizing the CO concentrations and emissions simultaneously, we performed a two-step inversion (see section 3) to mitigate the influences of model errors (i.e., the initial conditions). A description of the DA inversions system and the observation uncertainties are presented in section 2, and the experimental design is described in section 3. The results of the inversion on the emission estimate and the effects of the emission reduction policies on the spatiotemporal patterns of CO emissions are presented and discussed in section 4, and conclusions are summarized in section 5.

Method and Data
2.1. A Regional CO Assimilation System 2.1.1. Procedure of the Assimilation System Figure 1 shows the procedure of the cycling assimilations in this system. Since the lifetime of CO is relatively long, the errors in the simulated concentrations have a notable impact on the inversion result of the next step over a short time scale (Meirink et al., 2006;Peylin et al., 2005). Therefore, before the cycling assimilations, a 5-day spin-up simulation using the WRF/CMAQ model is conducted to generate a relatively good initial field, and then, the 3DVAR CO DA method within the Grid point Statistical Interpolation (GSI) framework of the National Centers for Environmental Prediction is used to further optimize the CO initial field (Kleist et al., 2009;Wu et al., 2002). The GSI CO 3DVAR framework is basically the same as that used by Feng et al. (2018), who constructed a WRF/CMAQ-3DVAR DA system for the simulation of PM 2.5 using surface measurements. During the cycling assimilations, instead of simultaneously optimizing the CO concentrations and emissions (e.g., Fortems-Cheiney et al., 2012;Y. Yin et al., 2015), in which the emission and concentration contributions to the model bias are hard to distinguish (Z. Jiang et al., 2017), in this study, we assume that all the biases between the simulations and observations are from the emissions since we have generated a relatively perfect initial field through the spin-up simulation and 3DVAR DA as mentioned before. In each DA window, a two-step method is applied. In the first step, based on the prior emission (X b ) and initial field, hourly CO concentrations are simulated using the CMAQ model, which are entered into the ensemble square root filter (EnSRF) algorithm along with the corresponding CO measurements, and the optimized CO emissions (X a ) are generated. In the second step, these optimized CO emissions are again entered into the CMAQ model to generate the initial field of the next DA window. It should be noted that in the first DA window, the prior emission is from emission inventory, which will be introduced in the next section, while in the subsequent windows, the prior emissions are from the optimized emission of the previous window. The DA window is set to 1 day because the model needs a longer time to integrate CO emission information into the concentration ensembles (Ma et al., 2019). In addition, due to the complexity of the hourly emissions, it is very difficult to simulate hourly CO concentrations that can match the observations, so daily mean simulations and observations are used in the EnSRF algorithm, and daily emissions are optimized in this study. Following X. Tang et al. (2013), the posterior CO emissions used for analysis in this study are temporally averaged using the optimized emissions of each DA window.

WRF/CMAQ Model
The "offline" regional air quality model WRF/CMAQ is applied in this study to simulate the meteorological fields and atmospheric compositions. WRF is a next-generation mesoscale numerical weather prediction system designed for both atmospheric research and operational forecasting needs (Skamarock & Klemp, 2008), and the CMAQ model is a regional 3-D Eulerian atmospheric chemistry and transport modeling system (Byun & Schere, 2006;Yu et al., 2012). With "one-atmosphere" design, CMAQ could address the complex interactions among multiple pollutants/air quality issues simultaneously. Many studies have shown that WRF/CMAQ has a good performance in extensive applications (N. Wang et al., 2015). In this study, the WRF/CMAQ settings are very similar to those of Feng et al. (2018). WRF simulations are performed across East Asia (169 west-east and 129 south-north cells) with a 36-km horizontal grid spacing ( Figure 2). Vertically, there are 51 sigma levels, with the model top fixed at 100 hPa. The CMAQ model is run with the same domain but with three gird cells removed from each side of the WRF domain. There are 15 layers in the CMAQ vertical coordinate, including approximately 7 layers in the boundary layer and the rests in the free troposphere, which are compressed from 51 WRF layers. In this study, WRF version 3.6 drives the CMAQ model. CMAQ version 4.7.1 (CMAQv4.7.1) optimizes the emissions and evaluates the effects of emission reductions, and CMAQ version 5.0.2 (CMAQv5.0.2) investigates the sensitivity of the emission estimates to the different model versions and parameter settings.
The meteorological initial conditions and lateral boundary conditions are both downscaled from the Final (FNL) Operational Global Analysis data of the National Centers for Environmental Prediction, which has a resolution of 1°× 1°at 6-hr intervals. The anthropogenic emissions over China are from the Multi-resolution Emission Inventory of China (MEIC), while those over the regions outside China are obtained from the mosaic Asian anthropogenic emission inventory (MIX) . The spatial resolutions of both the MEIC and MIX inventories are 0.25°× 0.25°. This combined emission inventory is employed as the prior emission of the first DA window (Figure 1). The biogenic emissions are calculated offline using the Model of Emissions of Gases and Aerosols from Nature (Guenther et al., 2012), which is also driven by the WRF model in this study. The Carbon Bond 05 chemical mechanism (CB05) is chosen as the gas phase chemical mechanism (Yarwood et al., 2005). Detailed physical and chemical configurations are listed in Table 1.

EnKF Assimilation Algorithm
EnKF uses ensembles of random samples to obtain error statistics of the model state variable or parameter. Error statistics can be propagated with linear or nonlinear dynamic models by simply implementing  (2002), is used to constrain the CO emissions in this study. In contrast to the traditional EnKF algorithm, the EnSRF algorithm obviates the need to perturb observations. The implementation process and setup are detailed below.
The uncertainty of the prior emissions is represented by a set of perturbed ensemble samples: where b represents background state and i is the identifier of the perturbed samples; δX b i represents the randomly perturbed samples that are added to the prior emissions X b 0 to produce ensemble samples of the inputs X b i . δX b i is drawn from Gaussian distributions with a mean of zero, and the perturbation magnitude is assumed to be 40% of the prior emission in each grid, which is set according to the work of Miyazaki et al. (2012). The ensemble size N is set as 30, which is a compromise between the computational efficiency and assimilation performance. After obtaining the ensemble of state vectors, CMAQ ensemble runs are conducted to propagate these errors in the model with each ensemble sample of state vectors from equation 1 as inputs. The background (a priori) error covariance P b is calculated based on the forecast ensemble from equation 2: where X b represents the mean of the ensemble samples. Based on the background error covariance, the response of the uncertainty in the simulated concentrations to the uncertainty in emissions is obtained. Combined with observational vector y, the state vector (prior emissions) is updated according to the following equations based on inverse theory and methods (Tarantola, 2004): While employing sequential assimilation and independent observations, where H is the observation operator that maps the state variable (CO concentrations) from the model space to the observation space, K is the Kalman gain matrix of the ensemble mean depending on the background and observation error covariance R, representing the relative contributions to analysis, and e K is the Kalman gain matrix of the ensemble perturbation. Thereafter, emission perturbations after inversions δX a i can be calculated. At the analysis step, the ensemble mean X ā is taken as the best estimate of the CO emissions.  To reduce the computational cost and the influence of representativeness errors, "superobservation" approach is implemented based on the optimal estimation theory (Miyazaki et al., 2012). Generally, the EnSRF works relatively well for sparse observational networks (Houtekamer & Zhang, 2016) or for observations with model-data mismatch errors uncorrelated (Houtekamer & Mitchell, 2001). However, in this study, as shown in Figure 2, the observation sites are in cluster distributions. Due to the low model spatial resolution (36 km), cluster sites (in one city) are usually located in a same model grid, which leads to the observations in one cluster may have significantly correlated or fully consistent model-data mismatch errors, resulting in repeated adjustments in a same direction for the emissions near one cluster sites. A test without "superobservation" method shows that the inversion performance will significantly decrease ( Figure S2 in the supporting information and Table 4). A superobservation is generated by averaging all observations located within a superobservation grid (36 km). Equations 7-9) show the calculation of a superobservation. We assume that the observation errors of the different stations at different times are independent of each other.
where j is the identifier of observations, m is the number of observations within a superobservation grid, y j represents observations, r j is observational error of y j , x ij represents a simulated concentration corresponding to the jth observation y j and modeled using the ith prior sample, and w j ¼ 1 r j 2 is the weighting factor. y new is a superobservation, r new is the corresponding observational error, and x new,i is the corresponding simulations of the ith ensemble sample. The superobservation error decreases as the number of observations used for the superobservation increases. To avoid storing and inverting very large matrices during analysis, superobservations are assimilated into the model in a sequential way with the observations assimilated one after another, which have better performance than assimilated simultaneously as long as the observation errors are uncorrelated (Houtekamer & Mitchell, 2001).
The model-data mismatch error stems not only from the state vector and observations but also from the model error. Neglecting model errors in the EnKF may lead to filter divergence, leading to a small ensemble spread and disregarding observations during analysis. Many studies use covariance inflation to maintain a suitable ensemble spread (H. Liu et al., 2017;Schwartz et al., 2014). However, this method lacks a physical basis and induces a spurious linear increase of background error in areas far from observation sites (X. Tang et al., 2011). In this study, considering compensation of the model error and prior inventory and daily emission uncertainties, we alternatively perturb emissions at each DA window with a same uncertainty of 40% . This method is similar to that used by Miyazaki et al. (2012), whereby the analysis spread is artificially inflated to a predefined value to prevent underestimation of error covariance during the analysis step. Additionally, Ma et al. (2019) also imposed a same spread at each DA window (time-invariant spread) to avoid filter divergence. The advantage of this method is that Gaussianity is maintained all the time. It should be noted that the uncertainty of 40% is not one of the raw inventories (MEIC inventory), which might be significantly greater than 40% . As shown in Figure 1, except for the first DA window, the prior emission used in each window is from the optimized emission of the previous window. Therefore, the emission analysis is generally no longer sensitive to the prior error of original emission inventory after several assimilation cycles (Gurney et al., 2004;Miyazaki et al., 2012). This setting of 40% is the same as the one of Miyazaki et al. (2012). Two additional sensitivity tests with prior uncertainty of 20% and 60% are conducted in this study. The results show that there is only very small differences in the  Figures S5 and S6 and Table S1), which is largely due to the two-step method used in this study that the inversion errors of the previous window could be transferred to the current one for further optimization (Figure 1).
Covariance localization is performed to reduce spurious correlations due to sampling errors caused by the finite ensemble size. Covariance localization determines that only measurements located within a certain distance (cutoff radius) from a grid point will influence the analysis of this grid. By using the Schur product of the ensemble-based covariance with a smooth correlation function, the impact of more noises and discontinuities in the analysis close to cutoff boundary could be reduced. Localization scale is an important parameter in an assimilation system. Generally, large localization scales may produce spurious increment of emissions, while small localization scales may cause insufficient emission adjustment in the distance and overadjusted emissions in the vicinity to compensate for simulation errors. An optimal localization scale is related to the ensemble size, the assimilation window, the dynamic system, and lifetime of chemical species in the atmosphere. CO is rather stable in the atmosphere, with lifetime more than 1 month, especially in winter (Minschwaner et al., 2010;W. Tang et al., 2019). Therefore, the localization scale is mainly depended on the transport of CO within one assimilation window (here, 1 day). According to Table 3 in section 4.1, on average, the simulated mean wind speed near the surface (10 m) is 3.6-3.86 m/s, and the corresponding transport distance is about 311-333 km. Considering the whole boundary layer, the mean wind speed will be greater than the ones near the surface, as we use a localization scale of 360 km (10 grid boxes).

Observation Data and Errors
Hourly As in Z. Jiang, Liu, et al. (2013) and Feng et al. (2018), to ensure data reliability before assimilating the CO data, we implemented two quality control steps in this study, including value range and time-continuity checks. First, only measurements in the range of 0.1~7 mg/m 3 are assimilated. The reason for setting data bounds is that the simulated concentrations represent an average value of one model grid (here, 36 km), which needs observations that must also have a certain spatial representation. Generally, too high or too low observations are considered to be without spatial representations, which may be caused by a sudden emission event, a huge wet deposition by extreme rain, or abnormal observations. In this study, the lower limit (0.1 mg/m 3 ) was set based on the background concentrations (Pani et al., 2019), and the upper bound (7 mg/m 3 ) was set based on experience. We checked the data carefully and found that most observations are lower than 7 mg/m 3 . Then, the time continuity is checked to eliminate gross outliers using the function of | y(t) − y(t ± 1)| ≤ f(t), where y(t) and y(t + 1) represent the observations at times t and t + 1, respectively, and f(t) equals 0.5 + 0.15y(t), which is determined empirically. The reason for applying the continuous check is similar to the value range check. Generally, due to the limitation of model spatial resolution and the assumption of a unify diurnal variation of emissions, the model cannot capture a sudden variation of the concentrations. Figure S1 shows the rejected data of four selected sites due to the value range check, and Figure S2 shows the rejected in one site of Nanjing as an example due to the time continuity check. Overall, the value range check rejects 1.2% and 0.1%, and the time continuity check rejects 4.4% and 2.4% of the total amount of data in 2013 and 2017, respectively. As shown in Figure S2, most data rejected by the time continuity check are real signals; therefore, an additional sensitivity test was conducted for the impact of the time continuity check on the inversion results for December 2017, in which the time continuity check was closed. The results show that the impact is very small ( Figure S3); the estimated CO emissions will only increase by 0.02% across the whole China mainland. After these two steps of quality controls, the sites with variable data less than a half are also removed. In addition, to maintain the inversion consistency, only those sites that have observations in both 2013 and 2017 are assimilated in this study. Finally, 681 measurement sites are assimilated for each year across the mainland China. Figure 2 shows the locations of the CO measurement sites.
The observation error covariance (R) contains both measurement and representation error contributions. To associate higher CO values with larger measurement errors, similar to the work of Schwartz et al. (2012), the measurement error ε 0 is defined as 10.1029/2019JD031808 Journal of Geophysical Research: Atmospheres where ermax is a base error, which is set based on Hooghiemstra et al. (2011) and Peng et al. (2018). Hooghiemstra et al. (2011) used a base error of 1.5 ppb, but they pointed out that a larger error is more robust, and Peng et al. (2018) applied a larger one of 10 ppb. In this study, it is set to be 0.01 mg/m 3 (8 ppb), which is in the range of the previous two studies. Π 0 denotes the observed CO concentration (unit: mg/m 3 ). Representativeness errors arise due to inaccuracies in the observation operator, which is parametrized using the equations of Elbern et al. (2007) defined as where γ is a tunable scaling factor (γ = 0.5 is used), Δl is the grid spacing (36 km), and L indicates the range that an observation can reflect, which is dependent on its location (here, 3 km for simplication). The total CO observation error (r) is defined as r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ε 0 2 þ ε r 2 p .

Evaluation Data and Methods
The hourly meteorological data, including temperature and relative humidity at 2 m (T2 and RH2) and wind speed at 10 m (WS10), are obtained from the National Climate Data Center integrated surface database (http://www.ncdc.noaa.gov/oa/ncdc.html). Meteorological data from a total of 35 sites are downloaded for validations of the meteorological field simulations. The locations of these sites are shown in Figure 2.
As presented in section 2.3, hourly averaged surface CO concentrations from 938 and 1,505 measurement sites are obtained for December 2013 and 2017, respectively. The sites with observations in both years (681 sites) are assimilated, and those with only data in 2017 are used for independent evaluation. Additionally, the CO observations in January 2014 and 2018 for those assimilated sites are also obtained for independent validation. As introduced in section 2.3, before those data are used, two quality controls steps are implemented.
During the evaluation, three basic statistical measures, the mean bias (BIAS), root-mean-square error (RMSE), and correlation coefficient (CORR), are calculated against the surface meteorological and CO observations to validate modeled meteorological fields and posterior emissions, respectively. The BIAS represents the overall model tendency, expressed as where x j and y j denote the modeled and the observational values, respectively, at the jth out of M locations, and the overbars denote averages. The RMSE reflects both the model bias and error variance, which is defined as The CORR quantifies the linear correspondence between the modeled and observational values, expressed as When compared with the observations, the model values are extracted from the model grids of the lowest level where the corresponding observation sites are located.

Experimental Design
Herein, the study periods are December 2013 and 2017.  Table 2).
AEP4 is run fully according to the procedure and settings as descripted in section 2.1, and AEP5 is conducted the same as AEP4 but using CMAQv5.0.2. CEPD and VEPD are directly run with the combined emission inventory as descripted in section 2.1.2 and the posterior emissions of AEP4, respectively. Both experiments use the same initial field at 0000 UTC on 1 December, which is generated through the spin-up simulation runs and GSI 3DVAR DA. These two experiments are conducted to investigate the performance of the assimilation experiments. CEPJ and VEPJ are also directly run with the December origin and posterior emission inventories, respectively, of AEP4. The initial field of these two experiments at 0000 UTC on 1 January is generated using the same method as before. These two experiments are run for an independent validation under the assumption that the emissions for subsequent January are the same as those for December. Table 2 lists all experiments conducted in this study.

Evaluation of the Simulated Meteorological Fields
Meteorological processes notably affect pollutant transport, photochemical reactions, and aerosol processes and determine the estimation of the flow-dependent background error covariance. Therefore, the biases of the meteorological simulations may highly affect the assimilation estimates of the CO emissions. Generally, both higher temperature and lower relative humidity lead to a faster oxidation of CO, and a higher wind speed corresponds to better diffusion and transport of air pollutants, resulting in lower CO concentrations in the atmosphere. Because all the biases between the simulated and observed concentrations are attributed to the emissions during assimilation, the lower concentrations caused by the higher wind speeds and temperatures and the lower relative humidity may lead to overestimation of the emissions.

Evaluation of the Posterior CO Emissions
Direct validation of the optimized emission inventories is impossible, and as a rule, we indirectly validate the posterior emissions by comparing the forward simulated atmospheric concentrations against measurements (e.g., F. Jiang et al., 2014;Jin et al., 2018;Koohkan & Bocquet, 2012;Peters et al., 2007). Based on this method, three groups of CO measurements over China are adopted, and two experiments (VEPD and VEPJ, Table 2) are conducted to evaluate the posterior CO emissions. For comparisons, the atmospheric concentrations simulated using the prior emissions are also evaluated in this section. For the three groups of CO measurements, one group consists of those assimilated during the two study periods, and the two other groups consist of the independent observations, namely, the remaining observations in December 2017 (all data are assimilated in 2013) and the measurements in January 2014 and 2018. Figure 3 shows the spatial distributions of the mean biases between the observed and simulated CO concentrations of the different experiments in the different years. For the CEPs (Figures 3a and 3b), most sites have negative biases (the simulated values are lower than the observed values; the same hereafter) in both years, and generally, the negative biases over northern China are larger than those over southern China. For the differences between these two stages, the biases in 2013 are basically larger than those in 2017, especially in the NCP and PRD. In the NCP, the negative biases at more than half of the sites are larger than 1 mg/m 3 in 2013, while in 2017, the biases at most sites are in the range of approximately −1 to 0.5 mg/ m 3 . After DA, for both years, the biases at most sites are reduced to smaller than 0.5 mg/m 3 . In 2013, the biases at 76.8% of the sites are within 0.5 mg/m 3 , and in 2017, the performance is even better than that of 2013, with biases at 94.1% of the sites within 0.5 mg/m 3 . Less improvement is mainly attained in Northwest and Southwest China, which only have very few sites, especially in 2013. We checked the CO records in these areas in the 2 years and found that in December 2013, many observations are missing, and the data are very discontinuous. In addition, it is found that there are some sites at which the bias is increased, indicating that the emissions in those areas may be overadjusted by DA. This could be explained by sites with opposite biases before DA in adjacent areas, and therefore, the contradictory adjustments partially decreased the site performance, resulting in larger biases. Statistics show that there are different levels improvements at 90.9% and 91.9% of the total 681 assimilation sites for 2013 and 2017, respectively.  adjustments with the monthly factors. It also could be found that in the YRD, in January 2014, the performance of the CO simulations in the VEPJ is even worse than the one in the CEPJ, with the absolute BIAS and RMSE increasing by 9% and 43%, respectively. One possible reason for this situation is that from December 2013 to January 2014, the CO emissions in the YRD are significantly reduced probably due to emergency control measures, since the observed CO concentration in the YRD decreases from 1.50 to 1.15 mg/m 3 , while in other regions, the observed CO concentrations between these 2 months are very similar. In January 2018, the improvements in the YRD are much better than those in 2014, with smaller BIAS, RMSE, and larger CORR. Although the BIAS is still much larger than those in the other regions, the RMSE and CORR in the YRD are close to or even better than those in the other regions. Figure 5 shows the spatial distributions of the biases calculated against the independent observations in December 2017. With the prior emissions, the distributions of the biases at these independent sites are similar to those at the assimilated sites (Figure 3b), and the BIAS, RMSE, and CORR values in all regions are also comparable (Table 4 vs. Table 6). Moreover, with the posterior emissions, in central eastern China, the distributions of the biases at the independent sites are also similar to those at the assimilated sites ( Figure 3d). However, in western China, due to the sparse assimilated sites, the biases at independent sites are much larger than those at the assimilated sites. Compared with the simulations with prior emissions, it could be determined that the improvements are very limited, indicating that the uncertainties of the posterior  (Table 4 vs. Table 6), indicating that the inferred emissions in these areas are reasonable. were 224% and 271%, respectively, in North China, and 600% and 376%, respectively, in Northwest China. The largest underestimation in northern China may be partly attributed to residential coal combustion, which is very common in the rural areas in northern China (S. Chen, Xu, Zhang, et al., 2017). B. Zheng et al. (2018) pointed out that the residential sector is one of the main sources of CO, and it has the largest uncertainties (Y. Zhao et al., 2012). Zhi et al. (2017) conducted a village energy survey and revealed a huge amount of missing rural raw coal for winter heating in northern China, which implies an extreme underestimation of the rural coal consumption by the prior emissions. In the NCP, YRD, and PRD, the relative increases in CO emissions are 174%, 168%, and 210%, respectively, for 2013 and 183%, 126%, and 109%, respectively, for 2017. Underestimation of the prior emissions is partly attributed to the fact that we use the emissions in December 2012 and 2016 as a priori, which are both 1 year behind the two assimilation months. However, considering the differences between the prior emissions in 2012 and 2016 and the differences between the prior and posterior emissions (the former is much lower than the latter), we believe that the underestimation caused by time 1 year ahead of the assimilation is basically negligible, and the CO Note. Statistics are calculated against daily regional-averaged observations (31 pairs), which are from 681, 159, 112, and 63 sites over four regions, respectively.

10.1029/2019JD031808
Journal of Geophysical Research: Atmospheres   The decreases mainly occur in North, East, Central, and Northwest China, with change rates of −6%, −24%, −38%, and −46%, respectively. In Northeast, South, and Southwest China, in general, the emissions are increased. In the NCP, YRD, and PRD, the reductions in CO emissions are −14%, −27%, and −29%, respectively. The total decrease in these three regions is account for 35% of the decrease in mainland China. The significant emission reductions in these three regions may be mainly due to the strict emission control strategies, including moving factory from urban regions to remote regions and changing the energy consumption and energy structure (Bian et al., 2019;D. Chen, Liu, et al., 2019;van der A et al., 2017). In BTH, the emissions are consistently decreased, but in its surrounding areas, including Shanxi, eastern Shandong, the junction area of Hebei, Shandong, and Henan provinces, as well as the rural areas of Tangshan and Qinhuangdao, the emissions are increased. Additionally, in the PRD, the emissions are consistently reduced, while in the eastern and western parts of the PRD, that is, eastern and western Guangdong Province, the emissions are  Journal of Geophysical Research: Atmospheres significant increased. In addition, it could also be observed that in or around most of the 74 key cities, the emissions are significantly reduced, while in the areas surrounding or near those larger cities, the emissions are increased. For example, in the SCB, the emissions are decreased in Chengdu but are significantly increased in northern Chengdu; in Hubei, the emissions are decreased in Wuhan but are increased in the surrounding areas. These results indicate emission migration from key urban areas or developed regions to their surrounding areas or developing regions.
Although the inventories based on the "bottom-up" methodology have large uncertainties, these inventories based on the same method may reflect the emission change characteristics to a certain extent. Figure 8 shows the differences in the prior emissions (MEIC 2016minus MEIC 2012. Compared with the emission changes in prior emissions, the changes in posterior emissions between December 2013 and 2017 are more complex and notable. In the prior emissions, the changes across the different areas consistently decrease in northern China and consistently increase in southern China, with the most significant increase in Guizhou. In addition, no emission migration, that is, significantly reduced emissions in key cities and regions and notable increases in their surrounding or nearby areas, is observed based on the changes in prior emissions. For the whole mainland of China, the changes in prior emissions between the 2 years are comparable to those in posterior emissions (−15% vs. −17%, Table 7). In the three key regions, in the NCP, the changes in prior emissions and the changes in posterior emissions are also comparable (−16% vs. −14%), while in the YRD, the change in posterior emissions is 2 times the change prior emissions, and in the PRD, the prior emissions show a slight increase (5%), while the posterior emissions exhibit a significant decrease (−29%). For the different areas in China, the prior emissions exhibit the most significant reduction in Northwest China, which is consistent with the posterior emissions in this study, but the reduction rate is lower than one half of that of the posterior emissions (−20% vs. −46%). Emission reductions are also observed in North, East, and Central China in both the prior and posterior emissions, but in North China, the reduction in prior emissions is much larger than that in the posterior emissions, while in East and Central China, the opposite occurs. Moreover, in Northeast, South, and Southwest China, the emissions are reduced in the prior results but are increased in the posterior estimates.

Sensitivity Analysis of Model Uncertainty
Model uncertainty plays an important role in the emission estimates because it relates the emissions with the concentrations. Although remarkable advances have been made in air quality modeling, large uncertainties remain. Model intercomparison studies within the framework of MICS-Asia III revealed that the poorly/differently parameterized physical, dynamical, and chemical processes of models greatly contribute to the large CO variability Kong et al., 2019). To provide insights into the impact of different models on the uncertainties of inverted CO emissions, sensitivity tests are performed using the CMAQv4.7.1 and CMAQv5.0.2 models. Because the simulated CO concentrations using the CMAQv5.0.2 model are slightly higher than that using CMAQv4.7.1 model (not shown), the inferred emissions in December 2013

10.1029/2019JD031808
Journal of Geophysical Research: Atmospheres and 2017 are both less than those using the CMAQv4.7.1 model. Compared with the posterior emissions using CMAQv4.7.1, the posterior emissions using CMAQv5.0.2 over mainland China, NCP, YRD, and PRD are decreased by 6.7%, 14.0%, 3.8%, and 8.6%, respectively, for December 2013 and 9.4%, 16.4%, 3.1%, and 7.6%, respectively, for December 2017. The highest sensitivity occurs in the NCP, which may be related to the complex air pollution mechanisms. For the posterior emission changes between December 2013 and 2017, there are no significant differences in the spatial distributions between the results inferred using CMAQv4.7.1 and CMAQv5.0.2. For the whole mainland of China, the emissions are reduced by 19.8%, while the reduction rates in the NCP, YRD, and PRD are 16.0%, 26.6%, and 28.5% using CMAQv5.0.2, respectively. These change rates are all comparable with those inferred using CMAQv4.7.1.  Journal of Geophysical Research: Atmospheres 4.6. Discussion The emission migration from urban or developed regions to the surrounding areas or developing regions found in this study is consistent with the fact that in recent years. With the acceleration of industrial transformation and upgrading, the adjustment of the economic structure of China is reflected not only in the differentiation of industrial structure that developed regions promote industrial structure to gather in low-polluted high-end service industries but also in the relayout and large migration of the industrial chain (Zhao & Fan, 2019). Traditional heavy industries (e.g., coal-fired power, iron, and steel industry) are the main sources of CO emissions (Dickerson et al., 2002;H. Liu et al., 2018;Streets et al., 2006). Therefore, the migration of industry is bound to affect CO emissions significantly.
In North China, due to the impact of capacity removal and production restriction , the area around Beijing has experienced production capacity clearance. A large number of iron and steel industries as a whole have gathered in Hebei and other regions with a higher production efficiency (G. Liu et al., 2016;H. Zheng et al., 2017). The coal industry has undergone similar capacity removal. These industries gather from BTH to Inner Mongolia and Shanxi (e.g., Demonstration Zone of Undertaking Industrial Transfer Cooperation in the East Mongolia and Demonstration Zone for Undertaking Industrial Transfer in the Golden Triangle of the Yellow River), where the coal industry is a pillar of the economic Zhan et al., 2019), leading to an increase in CO emissions in parts of these regions.
In East China, especially in the YRD, it has been reported that some industries have been transferred to Pan-YRD areas (e.g., the Demonstration Zone of Undertaking Industrial Transfer Cooperation in the Anhui River City Belt) and Jiangxi (e.g., the Demonstration Zone of Undertaking Industrial Transfer Cooperation in South Jiangxi) (Chang et al., 2013;H. Zhao et al., 2014). This transfer not only could be reflected in the Development Plan of Urban Agglomeration in the YRD (National Development and Reform Commission, 2016) but also coincides with the increases in CO emissions in central Anhui and Northeast Jiangxi Provinces (Figure 7). In Southern China, the CO emissions decreased in the PRD but increased in eastern and western Guangdong Province and southern Guangxi Province. Guangdong government proposed the macrocontrol strategy of "Bi-transfer" in which the eastern and western wings of Guangdong and the mountainous area of northern Guangdong actively undertake industrial transfer in the PRD to solve the problem of unbalanced regional economic development and promote ecological civil civilization (Bian et al., 2019;Jing et al., 2017;X. Yin et al., 2017). This may lead to increased emissions from the center to the periphery (L. Chen, Xu, & Yang, 2017;Wei et al., 2019;Xiao et al., 2017). Additionally, Guangxi (e.g., the Demonstration Zone of Undertaking Industrial Transfer Cooperation in East Guangxi) also undertakes the industrial transfer from Guangdong.
In the central and western regions, relying on abundant energy sources, specific industrial bases, and increasingly accessible transportation networks, a large number of traditional industries from the YRD and PRD are transferred to the central hinterland (e.g., the Demonstration Zone of Undertaking Industrial Transfer Cooperation in Jingzhou, Hubei). In Central China, many automobile manufacturing and iron and steel industries gather in Hubei (Li, 2014). In Northwest China, under the national strategy for energy development and safety (National Development and Reform Commission, 2012, National Energy Administration, 2014), energy industries (e.g., coal mining and thermal power) have been expanded and relocated toward energy resource-abundant (e.g., coal, oil, and natural gas) areas such as the northern Energy Golden Triangle (intersection of Inner Mongolia, Shaanxi, Gansu, and Ningxia provinces) and northern Xinjiang Province. A mass migration of industries has led to increased emissions in these areas. Ling et al. (2017) also revealed that the CO emissions in these areas had significantly increased based on satellite measurements.
Nevertheless, the emission changes still have large uncertainties. Generally, if these changes are believable, they should be greater than the uncertainties of the posterior emissions. As shown in Figure 3, the overall bias in 2013 for the concentrations simulated with posterior emissions is slightly greater than the one in 2017, indicating that the posterior uncertainties in 2013 are larger than those of 2017. Therefore, the emission changes and the posterior emission uncertainties in 2013 are compared (Figure 9). It could be found that the emission changes on the grid scale are larger than the emission uncertainties in most parts of central and eastern China, especially in the areas with dense observations. Instead, in western China, except for some key cities, the emission changes are lower than or comparable with the emission uncertainties. These indicate that in North China, East China, and South China, the uncertainties of the emission changes may be relatively low. In Southwest China, although the differences in west Sichuan, most Yunnan and Tibet are very small, the emissions in these areas are also very weak, indicating that the uncertainty may be also low in Southwest China. Oppositely, in Central China, Northeast China, and Northwest China, due to the sparse site distribution in the areas where they have significant emissions, the uncertainties of emission changes in these three regions may be quite large. For example, in Central China, the differences between the emission changes and posterior uncertainties are very small in southern Henan, western Hubei, and most Hunan, but the emissions in these areas are significant ( Figure 6). Furthermore, the uncertainties of emission changes in different regions may also be explained using the concentration biases in 2013 and 2017 (Figures 3c and 3d) at a certain degree. Generally, a similar bias in both years may suggest a low uncertainty with the emission changes, while a significantly different bias may suggest a large uncertainty. For example, in central Shaanxi province, the biases of all sites in 2013 are very large, while in 2017, they are rather small, indicating that the emission in this area was significantly overestimated in 2013 but well estimated in 2017, which may result in an overestimation of the emission change. Central Shaanxi is the most economically developed area in Northwest China. Therefore, the largest emission change in Northwest China shown in Table 3 (−46%) is overestimated in this study.

Summary and Conclusions
In this study, we construct a regional air pollution DA system based on the WRF/CMAQ model and the EnKF algorithm to quantitatively optimize gridded CO emissions using hourly surface CO observations across China obtained from the MEP network. The CO emissions in December 2013 and 2017 are optimized. We find that assimilation of the CO observations could significantly improve the CO simulations and emission estimates. Compared with the different types of CO observations, the BIAS and RMSE could be reduced by 66.7-98.5% and 12.1-43.2%, respectively, and the CORR was increased by 26.9-78.4% over mainland China. After assimilation, the uncertainties in the CO emissions in December 2013 and 2017 are comparable.

Journal of Geophysical Research: Atmospheres
China. For the whole mainland of China, NCP, YRD, and PRD, the increase rates are 186%, 174%, 168%, and 210%, respectively, in 2013 and 178%, 183%, 126%, and 109%, respectively, in 2017. Compared with the posterior emissions between December 2013 and 2017, overall, the CO emissions in December 2017 are lower than those in December 2013, with a mean decrease rate for the whole mainland of China of −17%. Large emission decreases mainly occur in key urban areas and developed regions. For the three key regions of the NCP, YRD, and PRD, the CO emissions are decreased by 14%, 27%, and 29%, respectively. Emission increases mainly occur in their surrounding areas and certain central and western regions. These changes are basically consistent with the implemented emission control strategies, and industrial transformation and upgrading have occurred in recent years, probably reflecting the emission migration from developed regions or urban areas to developing regions or surrounding areas. However, these changes are not found in the prior emissions, indicating that the posterior emissions could better reflect the temporal and spatial variations.