Quantification of SO₂ Emission Variations and the Corresponding Prediction Improvements Made by Assimilating Ground-Based Observations

Abstract

In this research, a new time-resolved emission inversion system was developed to investigate variations in SO₂ emission in China during the COVID-19 (Corona Virus Disease 2019) lockdown period based on a four-dimensional variational (4DVar) inversion method to dynamically optimize the SO₂ inventory by assimilating the ground-based hourly observation data. The inversion results obtained were validated in the North China Plain (NCP). Two sets of experiments were carried out based on the original and optimized inventories during the pre-lockdown and lockdown period to quantify the SO₂ emission variations and the corresponding prediction improvement. The SO₂ emission changes due to the lockdown in the NCP were quantified by the differences in the averaged optimized inventories between the pre-lockdown and lockdown period. As a response to the lockdown control, the SO₂ emissions were reduced by 20.1% on average in the NCP, with ratios of 20.7% in Beijing, 20.2% in Tianjin, 26.1% in Hebei, 18.3% in Shanxi, 19.1% in Shandong, and 25.9% in Henan, respectively. These were mainly attributed to the changes caused by the heavy industry lockdown in these areas. Compared to the model performance based on the original inventory, the optimized daily SO₂ emission inventory significantly improved the model SO₂ predictions during the lockdown period, with the correlation coefficient (R) value increasing from 0.28 to 0.79 and the root-mean-square error (RMSE) being reduced by more than 30%. Correspondingly, the performance of PM_2.5 was slightly improved, with R-value increasing from 0.67 to 0.74 and the RMSE being reduced by 8% in the meantime. These statistics indicate the good optimization ability of the time-resolved emission inversion system.

1. Introduction

Over the last two years, the whole world has suffered from the COVID-19 pandemic, which changed the way people lived through the spread of the mutant Delta and Omicron variants. Lockdowns proved to be the most efficient control measure for cutting down the transmission of the virus. Thus, investigating the effects of lockdowns on chemical mechanisms and emissions has become a new hot topic in the atmospheric chemistry field across the whole world. The outbreak of COVID-19 has greatly influenced industry and transportation, leading to emission reductions and corresponding changes in the levels of pollutants, which have been captured by both satellite retrievals and ground observations [1,2,3,4].

The lockdown measures used to combat COVID-19 in China after the virus spread from Wuhan from 23 January to 17 February 2020, were extremely strict, and many researchers have focused on this period when investigating changes in emissions [4,5]. Huang, et al. [6] updated emission data based on dynamic economic and industrial activity levels for individual provinces in China, showing different reduction levels for gases and particulate matter. Jia, et al. [7] conducted black carbon (BC) inversion based on Bayes’ theorem, with 70% and 48% BC emission reductions seen in eastern China and northern China during the lockdown. Feng, et al. [8] inferred daily NO_x emissions using the EnKF DA (Data Assimilation) system, obtaining results showing that NOx emissions fell by more than 60% in many large cities and by 36% across the whole of mainland China due to the lockdown.

These changes in emissions bring into question the accuracy of the use of atmospheric chemical transport models (CTMs) to forecast air quality, as emission inventory is an important component of CTMs [3,6,9,10]. The traditional “bottom-up” inventory is based on the statistical data obtained for different sectors, such as the Task Force on Hemispheric Transport of Air Pollution (HTAP) and Multi-resolution Emission Inventory for China (MEIC) [11], which failed to account for these changes. This kind of inventory has generally moved away from real time for several years, since their formation is systematic and complicated. Considering the uncertainty of the emission factors and activity levels involved, the biases between different gridded emission inventories are relatively high [12,13,14]. In order to reduce the deviation, the “top-down” method was used to invert the emission inventory by adding observation data to update the original emission inventory, including data assimilation, multi-source data fusion, mass balance, and linear constraint [4,15,16,17,18].

In these inverse methods, data assimilation is relatively accurate and widely used [19,20,21,22]. There are two main methods used to assimilate emission inventory, 4DVar (four-dimensional variational assimilation) and EnKF (Ensemble Kalman Filter) [23]. Feng, Jiang, Wu, Wang, Ju and Wang [23] constructed EnKF air pollution data assimilation system based on the Community Multiscale Air Quality Modelling System coupled with the Weather Research and Forecasting (WRF/CMAQ); their results showed that CO emissions were reduced from December 2013 to 2017 over the mainland of China. Dai, et al. [24] dynamically updated the SO₂ emission data obtained for China based on a four-dimensional local ensemble transform Kalman filter (4D-LETKF) and the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem). As for 4DVar emission inversion, the adjoint model of the Goddard Earth Observing System–Chemistry (GEOS-CHEM-ADJ) has been widely used around the world in different territories and for different species. In China, Zhang, et al. [25] assimilated satellite observations of the NH₃ column concentration using the Tropospheric Emission Spectrometer (TES) based on GEOS-CHEM-ADJ to optimize Chinese anthropogenic NH₃ emissions. Zhang, Henze, Grell, Carmichael, Bousserez, Zhang, Torres, Ahn, Lu, Cao and Mao [17] used the same system to invert the black carbon (BC) emission inventory in East Asia.

However, so far there has been no research focused on emission inversion using the 4DVar method in China during the lockdown period. Furthermore, Chen, Mao, Hong, Zang, Chen, Zhang, Gan, Gong and Xu [18] pointed out that the benefits of emission adjustment for model improvement lasted longer than the initial condition assimilation period during the lockdown. Compared with other inversion methods, the backward adjoint model was better able to calculate the sensitivity of different species, including the chemical relationship. Thus, it is meaningful to evaluate the SO₂ emission changes in China using the 4DVar method. SO₂ emissions in China are mainly attributed to power plants and heavy industry based on coal consumption. In normal times, it is difficult to obtain the consumption level of coal and its variation. The lockdown presented an opportunity to investigate these energy changes and evaluate the inversion results, since the industry indicators and carbon emissions could be obtained from open channels. Industry-related emissions are closely related to carbon emissions, and COVID-19 temporarily reduced China’s CO₂ emissions by a quarter. This decrease in emissions was estimated to be 27% based on the operating rate of the steel and chemical industry in several cities in Shandong and Hebei provinces, which are located in NCP [26,27]. In addition, the inversion inventory for this period has a more accurate temporal and spatial distribution and will support studies focusing on the mechanism of formation of secondary inorganic aerosol.

This paper describes a new study carried out using a 4DVar time-resolved emission inversion system to quantitatively estimate the emission changes in SO₂ by using ambient observations made during the COVID-19 lockdown and to update the original emission inventory. The structure of this paper is as follows. Section 2 describes the methodology used for the inversion method, including the descriptions of the China Meteorological Administration Unified Atmospheric Chemistry Environment (CUACE) 4DVar time-resolved emission inversion system and the model setup and observation data. Section 3 presents the evaluation of the different inventories obtained for the lockdown and pre-lockdown period, the effects of the lockdown on SO₂ emissions and concentrations, as well as sulfate concentrations in the NCP region. Section 4 gives our conclusions and a summary of this research.

2. Methodology

2.1. The Priori Inventory and Observation Data

The anthropogenic emissions were derived from the MEIC emission inventory (http://www.meicmodel.org/, accessed on 3 February 2022), which is one of the most frequently used inventories in China, including information on CH₄, CO, SO₂, NO_x, NMVOC, NH₃, PM₁₀, PM_2.5, BC, and OC with a 0.25° × 0.25° grid resolution [11,28]. In this research, the latest MEIC inventory taken in January 2017 was adopted as the a priori inventory to estimate the posterior gridded SO₂ emission inventory for the lockdown and pre-lockdown periods.

The ground-based SO₂ observation data from the network of environmental monitoring stations which were operated under the Ministry of Ecology and Environment (MEE) of the People’s Republic of China (https://air.cnemc.cn:18007, accessed on 3 February 2022) were assimilated by the CUACE-4DVar (China Meteorological Administration Unified Atmospheric Chemistry Environment—four-dimensional variational assimilation system) time-resolved emission inversion system to optimize the original emission inventory. There were up to 245 monitor stations and over 50 grids in the assimilation domain. For the monitors located in the same model grid, the SO₂ concentrations were averaged with an hourly time resolution. The interpolation locations of the observation data to assimilate in the simulation domain are marked as red dots in Figure 1.

2.2. CUACE-4DVar Time-Resolved Emission Inversion System

WRF-CUACE (the Weather Research and Forecasting model coupled with China Meteorological Administration Unified Atmospheric Chemistry Environment) is a new-generation air quality model with air quality and meteorological components fully coupled in an “online” approach [28,29]. There are two gas mechanisms in version 1.0 of WRF-CUACE, including the second generation of the Regional Acid Deposition Model (RADM2) and Carbon Bond Mechanism Z (CBM-Z) [30,31]. In this research, the modified Carbon Bond Mechanism IV (CBM-IV) from CMAQ [32] was newly added to the updated version of the WRF-CUACE to match the adjoint model. The adjoint model of WRF-CUACE (WRF-CUACE-ADJ) consisted of a gas adjoint module (CBM-IV-ADJ) and an aerosol adjoint module (CAM-ADJ), which were migrated from CMAQ-ADJ (the adjoint model of CMAQ) and GRAPES-CUACE-ADJ (the adjoint model of Global/Regional Assimilation and Prediction System-CUACE), which have been widely used and verified in many studies [33,34,35,36]. The reverse wind was used for the calculation of the meteorological module in the backward integration [37]. In other words, the adjoint sensitivity includes the backward trajectory and the chemistry relationship at the same time. The flowchart, code structure, and other details of the updated WRF-CUACE and its adjoint model are presented in the Supplementary Materials (Figures S1 and S2). Considering the consistency, the backward steps are the same as the steps in the forward model, since the adjoint theory performs calculation in the same state. Thus, the time steps are kept the same for 5 min without variation in either the forward or backward model. The relationship of each time step between the forward simulation and backward simulation is shown directly in Figure S2.

The purpose of data assimilation is to minimize the cost function. As for the emission assimilation, the definition is presented in Equation (1). The 4DVar method and adjoint model are combined by the adjoint operator K^T, which is also called the adjoint sensitivity and represents the sensitivity of the concentration to related emissions. Thus, the adjoint model is the foundation of the 4DVar data assimilation system, which is the output of the adjoint model and the input of the optimization algorithm and influences the inversion performance. In addition, the observation data are another input to the assimilation process. Adjoint forcing is the difference between the observed and modeled concentrations, and the adjoint operator K^T is calculated by the adjoint code. Both of these drive the WRF-CUACE-ADJ to calculate the gradient of the cost function. After the continuous iterative adjustment of the optimization algorithm through observation data and adjoint sensitivity, the estimated inversion inventory is obtained to approximate the real situation.

Based on the WRF-CUACE-ADJ model, the CUACE-4DVar time-resolved emission inversion system was developed. The flowchart of its data assimilation process is shown in Figure 2. The blue rectangle represents a single DA window and was set to one day (24 h) to optimize the daily emission inventory with a resolution of 30 × 30 km. The framework includes three kinds of simulation: forward run, restart run, and backward run. First, a forward WRF-CUACE run with a 10-day spin-up was carried out to form the daily restart files at the positions of the red circles. Second, the restart run represented by the green arrow marked as 1 was operated at the positions of the blue circles after the restart files were revised by observation data. In this operation time, the state variables were stored at every dynamic time step to drive the adjoint run. Third, the backward run represented by yellow arrow 2 was carried out to calculate the gradient of the cost function (J), defined as the difference between hourly observation data and modeling concentration (Φ) multiplied by the adjoint operator (K^T). Last, an L-BFGS optimization algorithm was applied to minimize the J values by 1 + 2 forward and backward iteration loops. The iteration reached convergence when the adjacent delta J was less than 1%; in the meantime, the optimized emission inventory at this DA window was obtained. This system was applied to SO₂ in the current paper, while other components are still under development, since the assimilation parameters introduced in the following part and the chemical characteristics are different for different species and will include NO_x, VOC, and BC in the future.

Generally, the cost function J is defined by the following equation, including observation and background terms.

$$J\left(\mathit{x}\right)=\frac{1}{2}{\left({\mathit{y}}_{\mathit{m}}-{\mathit{y}}_{\mathit{o}\mathit{b}\mathit{s}}\right)}^T{S}_o^{-1}\left({\mathit{y}}_{\mathit{m}}-{\mathit{y}}_{\mathit{o}\mathit{b}\mathit{s}}\right)+\frac{1}{2}{\gamma }_r{\left(\mathit{x}-x_a\right)}^T{S}_a^{-1}\left(\mathit{x}-x_a\right)$$

(1)

where ${\mathit{y}}_{\mathit{m}}$ is the concentration of the forward model, ${\mathit{y}}_{\mathit{o}\mathit{b}\mathit{s}}$ is the observation data, $\mathit{x}$ is the posterior inventory, $x_a$ is the priori inventory. $S_o^{-1}$ and $S_a^{-1}$ are the observational and background error covariance, both of which are diagonal matrices, and ${\gamma }_r$ is the penalty term used to adjust the weight between observation and background. S_a is the uncertainty of the priori inventory, defined as 100% for SO₂, since the emission varied a lot during the lockdown. $S_o^{-1}$ is the sum of simulation bias and observation error, followed by the relative residual error (RRE) method [38,39]. The scaling factor σ was chosen to avoid the negative value that might occur in the optimization process, defined as σ = ln (x/x_a) [17,40]. The ${\gamma }_r$ was set to be 1–10 from the L-curve approach [41]. It generally takes 8 to 15 iterations to converge with a 25–40% reduction in the cost function [22,39].

The CUACE-4DVar time-resolved emission inversion system can invert both gas and aerosol emissions simultaneously. This is a great improvement compared with the original GRAPES-CUACE-ADJ and CMAQ-ADJ. The former only included the aerosol inversion of BC and the latter only included the assimilation of gas species. The newly developed inversion system has a better daily temporal resolution than the original monthly prior inventory. In the current situation of the rapid change in emission sources in China, the 4Dvar inversion system can rapidly capture the drastic emission changes, such as those observed for the lockdown period, resulting in a better performance of atmospheric chemistry models. Moreover, it is worth noting that the resolution of the inventory could be defined by the user, leading to the improvement of grid resolution and a better distribution of emissions. Benefitting from the framework of WRF, the CUACE-4DVar inversion system could implement parallel computing, which greatly improves the efficiency of calculation.

2.3. Model Setup and Experiments Design

The model domain used in this study covered central and eastern China, as shown in Figure 1, with a resolution of 30 km, which was sufficiently large to consider the pollutant transboundary transport from outside NCP. There were 32 vertical layers, with the top pressure level being at 100 hPa. Based on previous lockdown studies, two individual weeks were chosen to represent the pre-lockdown and lockdown situations [9,42,43]. The simulation time of the pre-lockdown period was six days from 15 to 20 January in 2020, and the lockdown period was 4–9 February. The lockdown effect on the SO₂ emission changes was quantified by the differences between the pre-lockdown and lockdown period. Due to the high computing resource and storage requirements of 4DVar data assimilation with input and output at every time step, the simulation period did not have a one-month length. However, the features of the lockdown period are quite consistent with those of other studies [6,7,8].

The main physical options used in this research included the Lin microphysics scheme [44], the Rapid Radiative Transfer Model (RRTM) longwave radiation scheme [45], the Goddard shortwave scheme, the YSU (Yonsei University) boundary layer scheme [46], and the unified Noah land surface model [47]. The chemical options are described in Section 2.1. The FNL reanalysis data were provided by the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) with a resolution of 1° × 1°, and they were used as the meteorological initial condition and lateral boundary condition with an update frequency of every 6 h. The biological emissions showed no change during the simulation time, with a spatial resolution of 30 km, as calculated by MEGAN (Model of Emissions of Gases and Aerosols from Nature).

3. Results and Discussion

3.1. Spatial Distribution of the MEIC and Optimized SO₂ Emissions

The spatial distributions of the grid SO₂ emission for MEIC (January 2017) inventory and average optimized inventory during the pre-lockdown and lockdown periods are presented in Figure 3. The optimized results obtained during the pre-lockdown period by this study indicated that the MEIC inventory might underestimate the SO₂ emissions from 15–20 January (Figure 3b), which is mainly attributed to fireworks and sacrifice burning [48,49] during the week before the Spring Festival. The firecrackers and sacrificial offerings could increase the SO₂ emission in this special week, as was also found by [8] in a NOx emission study carried out over these days. The reduction grids for the SO₂ emission in urban locations confirmed this point of view, since fireworks and sacrifice burning were forbidden by the government in these regions.

During the lockdown period, the Chinese government took strict measures to keep people at home. Transportation was kept at an extremely low level, and industries were almost shut down, except for non-stop industries [26,42]. In the meantime, these phenomena led to a sharp decrease in anthropogenic emissions, including SO₂, and coal consumption was also reduced by about 30% [26]. Accordingly, this reduction was captured by the CUACE-4DVar time-resolved emission inversion system, as presented in Figure 3c. The optimized emission levels were mainly decreased in Beijing, Tianjin, Hebei, and parts areas of Shandong and Henan province. This distribution was in agreement with the changes in indicators of industrial activity for China during this period, while the operating rate varied from 75% to 56% for steel bar plants that were mainly located in Tangshan, Hebei province, and the oil refinery utilization in Shandong province varying from 66% to 48% [27]. However, the optimization result of SO₂ emission during the lockdown period showed a similar increment in Shanxi as the pre-lockdown period, which could illustrate that the MEIC 2017 inventory might have underestimated the SO₂ emissions. For Shanxi province, the SO₂ concentration was at a relatively high level compared to other key regions in China during recent years. This could be explained by the energy structure, as Shanxi was the main production area of coal [50], which was also pointed out by Cao and Cui [51] in the research on the Fenwei Plain’s pollution status. These variations in SO₂ emissions over time are also in agreement with the results from satellite retrieval [52].

As a result, the spatial distribution and value of the original MEIC inventory were updated by this study using observation data. Comparing the differences between the MEIC and the optimized inventory, the SO₂ emissions of the MEIC inventory were relatively low compared to the real situation in most regions during the pre-lockdown period. However, the distribution was higher in Beijing, Tianjin, and Hebei provinces during the lockdown period.

3.2. Variation of SO₂ Emission Due to the Lockdown

The variations in SO₂ emissions due to the lockdown measures were calculated based on the differences between the averaged daily amount of the optimized emission inventories during the pre-lockdown and lockdown periods; the results are shown in Figure 4. It can be seen intuitively from Figure 4a that the SO₂ emissions were reduced in the NCP region, except for some individual grids. There was commonly about a 50–300 mol/km²/h reduction with a ratio of 10 to 30% in this area. However, the distribution in central and northern Shandong province seems different from that in other regions in the whole domain, with a maximum increment of 300 mol/km²/h. This was mainly attributed to the enhanced heating requirements during the lockdown period, since Shandong province is also a coal-burning area, as with Shanxi [53]. The heating episode was captured by the Ecological Environment Bureau in Jinan, Shandong province, as the SO₂ concentration reached 18 mg/m³ at the boiler of the local thermal production center (http://jnepb.jinan.gov.cn/art/2020/3/20/art_27565_4292195.html, accessed on 3 February 2022). As presented in Figure 4b, the difference in the SO₂ emissions due to the lockdown measures showed a consistent reduction trend in the NCP. For the changes in some individual grids, the difference was over ±40%. These grids were mainly located in the northwest part of Hebei province and the adjacent regions of the major cities, where the original SO₂ emission levels were relatively low (Figure 3a–c). A small change would cause a great difference, but the actual increase value was not large.

Except for the gridded distribution, the accumulated amount of the SO₂ emissions for different provinces was calculated by the sum of the model grids located in the corresponding regions. The emission amount for individual regions during different periods is presented in the left part of Figure 5, and the reduction ratios are presented in

Table 1. The estimation of SO₂ emission reduction ratios in NCP and different provinces during the lockdown period.

Province	Ratio	Province	Ratio
Beijing	−20.7%	Shandong	−19.1%
Tianjin	−20.2%	Henan	−25.9%
Hebei	−26.1%
Shanxi	−18.3%	NCP	−20.1%

. The average reduction ratio was 20.1% for the whole NCP region, while the values were 20.7% in Beijing, 20.2% in Tianjin, 26.1% in Hebei, 18.3% in Shanxi, 19.1% in Shandong, and 25.9% in Henan province. The total ratios fluctuated around 20%, which was consistent with the reduction data provided by Huang, Ding, Gao, Zheng, Zhou, Qi, Tang, Wang, Ren, Nie, Chi, Xu, Chen, Li, Che, Pang, Wang, Tong, Qin, Cheng, Liu, Fu, Liu, Chai, Davis, Zhang and He [6] based on dynamic economic and industrial factors for the different provinces varying from 16% to 26%. Unlike the multiple sources of NO_x, the emission sources of SO₂ were mainly industries and power plants, and they could be used as indices to represent the industrial level [12]. Though the industry activity decreased, there were still some non-stop industries such as steel, petrochemical factories, and power plants, which were the main reasons for the air pollution that took place in the NCP during the lockdown period [42].

Furthermore, the daily variations in the SO₂ emissions for different regions as well as the corresponding daily SO₂ observations are shown in Figure 5. The trend of the daily emissions was consistent with the variation in the observations for the NCP region, which showed the good performance of the time-resolved capacity of the CUACE-4DVar emission inversion system. For each period in the NCP, the SO₂ emissions decreased first and then increased, while the variations were different for the individual provinces. There was an increment in SO₂ concentration in Beijing from 16 to 17 in January while the SO₂ emissions decreased. However, the emissions in Hebei province increased at the same time, which means that the SO₂ increment in Beijing might have been caused by the transport from Hebei province, as the chemical characteristics of the atmosphere were similar in both regions. This phenomenon was also found at other times for Tianjin, Shanxi, and Henan provinces, which indicates the importance of the pollutant transboundary transports in NCP.

3.3. Surface SO₂ and Sulfate Concentration Changes Due to the Lockdown

Except for the direct SO₂ emission changes, the surface SO₂ and sulfate concentrations simulated by the optimized inventories (Figure 3b,c) are analyzed in this section. The lockdown effect on the surface SO₂ and sulfate concentrations is represented by the concentration differences between the pre-lockdown and lockdown periods shown in Figure 6. As a response to the reduction in the emission amount in the NCP, the SO₂ concentration decreased by 4–8 µg/m³ in the main region. The decreases were most obvious in Datong, Tangshan, and Zibo city, at 10 µg/m³, as marked in Figure 6a. For Datong and Tangshan, the decrease was mainly attributed to the change in the heavy industries because the lockdown control affected the operating rate of the steel and coal consumption factories. As for sulfate, the reduction regions coincided with the SO₂-reduced locations. There was a decrease of 0.5–1.0 µg/m³ in Hebei and Shandong provinces, with the maximum value found in central Shandong. However, the main difference was seen for Shanxi province, where the reduction in sulfate showed a relatively lower value compared to the SO₂ decrease. This might be due to the oxidability of the ambient environment, since mining is the main industry in Shanxi province, while the manufacturing and chemical industries are located in Hebei and Shandong provinces. Thus, SO₂ was more likely to be converted to sulfate by heterogeneous reactions in Beijing, Tianjin, and Hebei provinces [54].

Additionally, we noted an interesting phenomenon that occurred in Zibo, which is located in central Shandong. Though the SO₂ emissions increased in the lockdown period (Figure 4), the concentration decreased at the same time. This phenomenon indicated that the increase in emission in several grids could not affect the total reduction trend in the surrounding regions, since the accumulated SO₂ emissions were reduced in Shandong province. Moreover, this also reflected the importance of the transboundary transport of air pollution. Comparing the average emission factor method, which varied by province during the lockdown, the CUACE-4DVar time-resolved emission inversion system not only had a better daily time resolution but also redistributed the emission more reasonably based on the adjoint sensitivity.

3.4. Simulation Improvement of the Optimized SO₂ Emission Inventory

Unlike the analysis of the historical inventory, the most common factor used for the evaluation of the near-real-time inventory was the CTM model performance. Since the observations of PM_2.5 components are barely related to location and time series, it is better to use PM_2.5 concentration rather than individual sulfate to reflect the model’s performance across the whole domain. Thus, the simulated SO₂ and PM_2.5 concentrations in the monitor locations were verified by the corresponding ground-based observations from the network of MEE introduced in Section 2.1. The root-mean-square error (RMSE) and Pearson correlation coefficient (R) were chosen as the main statistic indicators for evaluating the accuracy of the original MEIC and two optimized inventories. Additionally, mean fractional bias (NMB) and mean bias (BIAS) were calculated. These indices are widely used to validate the performance of the numerical model and inventory [8,24,55]. The calculation results obtained for different scenarios are presented in

Table 2. The statistics of model simulation and observation data of SO₂ and PM_2.5 concentrations for MEIC and optimized inventories during the lockdown and pre-lockdown periods.

Period	Statistics	MEIC	Optimization
SO2 (pre-lockdown)	R	0.46	0.87
	RMSE (µg/m3)	17.0	13.7
	NMB (%)	3.7	−39.0
	BIAS (µg/m3)	0.7	−7.8
SO2 (lockdown)	R	0.28	0.79
	RMSE (µg/m3)	12.7	8.8
	NMB (%)	23.1	−32.3
	BIAS (µg/m3)	2.6	−4.8
PM2.5 (pre-lockdown)	R	0.62	0.67
	RMSE (µg/m3)	61.3	56.2
	NMB (%)	−45.9	−40.1
	BIAS (µg/m3)	−37.8	−33.1
PM2.5 (lockdown)	R	0.67	0.74
	RMSE (µg/m3)	20.4	18.9
	NMB (%)	−31.2	−30.1
	BIAS (µg/m3)	−11.3	−10.9

, while the scatter plots are shown in Figure 7.

Both priori and posterior inventories showed the better relevance of SO₂ simulation during the pre-lockdown time compared to the lockdown period, as the emission changes were much greater during the lockdown period. However, the results obtained for PM_2.5 were the opposite, and all the statistics calculated during the lockdown period were better than those for the pre-lockdown period since the deviation of PM_2.5 was smaller during the lockdown. As for MEIC inventory, the simulation results obtained by WRF-CUACE for PM_2.5 were better than those obtained for SO₂ in terms of R-values but worse in terms of NMB. The SO₂ concentrations were relatively better in the pre-lockdown period than in the lockdown period. The R-value of SO₂ in the pre-lockdown period was 0.46 and passed the 99% significance test, but the quality was relatively poor, at only 0.28 during the lockdown period. However, the R-values for PM_2.5 were both over 0.6 during the whole simulation period.

It can clearly be seen from the scatter plots (Figure 7) and Table 2 that the optimized inventories showed the better performance of SO₂ and PM_2.5 in terms of both R and RMSE for all periods. Compared to the original MEIC inventory, the optimized SO₂ emission inventory improved the R from 0.46 to 0.87 while reducing the RMSE by 19.2% from 17.0 to 13.7 µg/m³ in the pre-lockdown period. There was an approximately 31.3% reduction in RMSE during the lockdown period, with R being improved significantly from 0.28 to 0.79. For PM_2.5, the improvement was small, as the sulfate accounted for a relatively small proportion, with R varying from 0.62 to 0.67 during the pre-lockdown period and from 0.67 to 0.74 during the lockdown period, while the RMSE was reduced by 8%.

Furthermore, both bias indices NMB and BIAS showed that the MEIC 2017 inventory overestimated the SO₂ concentration in all periods while the optimized inventory underestimated it, which was also found for the slopes of the scatter plot. As for PM_2.5, both MEIC and the optimized inventories underestimated the concentration. However, the optimized inventories reduced NMB and BIAS by 12% and 4% in the pre-lockdown and lockdown periods, respectively. These model performance values were comparable to those obtained in other emission inversion studies [4,24,39].

There was a certain degree of underestimation for the optimized inventory, and this might be due to the lack of a heterogeneous adjoint module for SO₂ sensitivity in the current WRF-CUACE-ADJ, since the optimization performance mainly relies on the accuracy of the adjoint sensitivity. More adjoint modules, such as heterogeneous reaction and thermodynamic equilibrium, should be added in the future to complete the inverse adjoint model and decrease the uncertainty.

Due to the high computing resource and storage requirements of 4DVar data assimilation with input and output at every time step, the simulation period did not have a length of one month. The simulation time chosen for the pre-lockdown and lockdown period was two specific weeks, and this might have caused some bias. However, the features and quantification results of this lockdown period are quite consistent with those obtained for other long-term modeling research.

The different combinations of parameterizations could have led to different modeling performances. The results obtained for SO₂ variation in this article were based on the CBM-IV gas mechanism and a distance of 30 km. Furthermore, the grid resolution could have affected the inversion results, since the calculation of the meteorological parameterization scheme, emission amount, and chemistry scheme was influenced by the resolution. As the emission amount was directly influenced by the grid area of the chemistry model in the emission rate equation, the improvement of the resolution will lead to a better spatial distribution of the sensitivity information.

As a result, the inversion inventory significantly improved the performance of the WRF-CUACE of SO₂, especially in the lockdown period, and slightly improved the performance of PM_2.5, since the ratio of the sulfate was relatively small, as shown in Figure 6b. Although the optimized inventory underestimated the SO₂ emissions to a certain extent, the correlation coefficient and simulation error were improved obviously, which indicated the good optimization ability of the CUACE-4DVar time-resolved emission inversion system. The inversion system decreased the lag and uncertainty of the priori “bottom-up” inventory and led to a better temporal and spatial distribution.

4. Conclusions

In this research, a time-resolved emission inversion system was developed based on the adjoint model of WRF-CUACE to optimize the MEIC emission inventory using the hourly surface SO₂ observations from MEE in China. The quantification of the SO₂ emission changes caused by lockdown measures in the NCP was characterized by the variation in the optimization inventories during the pre-lockdown and lockdown periods. We commonly found an approximate 10 to 30% reduction for the whole NCP region with an average ratio of 20.1%. The reduction amounts were 20.7% in Beijing, 20.2% in Tianjin, 26.1% in Hebei, 18.3% in Shanxi, 19.1% in Shandong, and 25.9% in Henan province, respectively, while the emissions in central and northern Shandong were increased due to the enhanced heating requirements during the lockdown period. The trend of the daily emissions was consistent with the variation in the observations in the NCP region, which showed the good performance of the time-resolved capacity of the CUACE-4DVar emission inversion system.

Furthermore, the SO₂ concentration decreased by 4–8 µg/m³ in the main NCP as a response to the emission reduction. The maximum reduction values reached 10 µg/m³ in Datong, Tangshan, and Zibo city, which we mainly attribute to the change in the heavy industry operating rate due to the lockdown controls imposed in these areas. The sulfate concentration changed correspondingly, though at a smaller magnitude. The evaluation results indicated that the daily optimization inventory of SO₂ significantly improved the model performance in terms of SO₂ during the COVID-19 pandemic in January and February 2020, with the R-values increasing from 0.28 to 0.79 and the RMSE being reduced by over 30% or about 4 µg/m³. Meanwhile, the corresponding changes in sulfate slightly improved the PM_2.5 simulation, with the R-values increasing from 0.67 to 0.74 and the RMSE being reduced by 8% or about 1.6 µg/m³. The inversion system might lead to the underestimation of the emissions, as the NMB and BIAS showed negative values.

As a result, the CUACE-4DVar time-resolved emission inversion system reduced the lag and uncertainty of the priori “bottom-up” inventory and had a better temporal and spatial distribution, especially during the lockdown period. Compared with the monthly averaged MEIC inventory, the optimized “top-down” daily inventory not only had a better daily time resolution but also redistributed the SO₂ emissions more reasonably based on the adjoint sensitivity in an arbitrary resolution that could be defined by the user. In the current situation of the rapid change in emission sources under the dual control and carbon policy in China, the 4Dvar inversion method has strong scientific research and practical application value.