Dynamics-based estimates of decline trend with fine temporal variations in China's PM_2.5 emissions

Timely, continuous, and dynamics-based estimates of PM_2.5 emissions with a high temporal resolution can be objectively and optimally obtained by assimilating observed surface PM_2.5 concentrations using flow-dependent error statistics. The annual dynamics-based estimates of PM_2.5 emissions averaged over mainland China for the years 2016–2020 without biomass burning emissions are 7.66, 7.40, 7.02, 6.62, and 6.38 Tg, respectively, which are very closed to the values of the Multi-resolution Emission Inventory (MEIC). Annual PM_2.5 emissions in China have consistently decreased by approximately 3 % to 5 % from 2017 to 2020. Significant PM_2.5 emission reductions occurred frequently in regions with large PM_2.5 emissions. COVID-19 could cause a significant reduction of PM_2.5 emissions in the North China Plain and northeast of China in 2020. The magnitudes of PM_2.5 emissions were greater in the winter than in the summer. PM_2.5 emissions show an obvious diurnal variation that varies significantly with the season and urban population. Compared to the diurnal variations of PM_2.5 emission fractions estimated based on diurnal variation profiles from the US and EU, the estimated PM_2.5 emission fractions are 1.25 % larger during the evening, the morning peak is 0.57 % smaller in winter and 1.05 % larger in summer, and the evening peak is 0.83 % smaller. Improved representations of PM_2.5 emissions across timescales can benefit emission inventory, regulation policy and emission trading schemes, particularly for especially for high-temporal-resolution air quality forecasting and policy response to severe haze pollution or rare human events with significant socioeconomic impacts.

1 Introduction

Anthropogenic emissions have imposed essential influences on the Earth system, from hourly air quality and human health to long-term climate and environment. To reduce anthropogenic emissions, the Chinese government has enforced the Clean Air Action (2013) since 2013 (China State Council, 2013). Studies to date that evaluated the emission controls and understood the climate responses to emission reductions have often used either a fixed meteorology with emission changes or vice versa (Li et al., 2019, 2021; Zhai et al., 2021). Estimated emissions from empirical extrapolation were commonly applied to analyse the meteorological–chemical mechanisms and associated social–economic impacts from occasional events like the 2015 China Victory Day Parade and Coronavirus Disease 2019 (COVID-19) pandemic (Wang et al., 2017; Liu et al., 2020; Huang et al., 2020; Zhu et al., 2021). But to better understand both long-term and short-term influences from emission changes, the continuous, up-to-date, and high-temporal- and high-spatial-resolution emission estimates with coherent interactions of meteorology and emission changes are needed.

The complex contributions from energy production, industrial processes, transportation, and residential consumptions have imposed great challenges on accurately estimating the emissions. The emission inventories created by the traditional bottom-up techniques were typically outdated compared to the present day due to the lack of accurate and timely statistics, and they often used coarse temporal resolutions, from monthly to annual (Zhang et al., 2009; Li et al., 2014; Janssens-Maenhout et al., 2015b; Zheng et al., 2018). Alternatively, up-to-date emission estimates with high temporal and spatial resolutions could be provided by top-down techniques (Miyazaki et al., 2017), but most emissions estimated by top-down techniques were intermittent and analysed at monthly scales or longer (Zhang et al., 2016; Jiang et al., 2017; Qu et al., 2017; Cao et al., 2018; Müller et al., 2018; Chen et al., 2019; Li and Wang, 2019; Miyazaki et al., 2020). Moreover, emissions updated by the top-down techniques based on satellite observations could be insufficient to capture realistic near-surface characteristics (Li and Wang, 2019; Liu et al., 2011; Choi et al., 2020).

Given the development of observation networks and advanced data assimilation strategies, timely and dynamics-based emission estimates with high temporal resolutions can be achieved by harmonically constraining the atmospheric–chemical model with dense observations of trace gas compounds through an optimal assimilation methodology. The ensemble Kalman smoother (EnKS) (Whitaker et al., 2002; Peters et al., 2007; Peng et al., 2015), as a four-dimensional (4D) assimilation algorithm, makes use of chemical observations from the past to the future to provide an optimal estimate of source emissions, and it can capture the “error of the day” and construct fine emission characteristics with high temporal and spatial resolutions by using short-term ensemble forecasts (Kalnay, 2002). Since 2013, the fine particulate matter pollution (PM_2.5, particles smaller than 2.5 µm in diameter) as the most urgent threat to public health has persistently decreased, and ground-based observations of PM_2.5 have progressively increased (Huang et al., 2018). Thus, by harmonically assimilating dense surface PM_2.5 observations into an atmospheric–chemical model through an EnKS, hourly estimates of PM_2.5 emissions that were continuously cycled for the years 2016–2020 are presented in this study.

The timely estimated emissions can provide guidance for emission inventories that usually have time lags and emission trading schemes that often require up-to-date source emissions. Based on the dynamics-based estimated emissions with harmonic combination of the model and observations, a better evaluation of the emission controls and a more comprehensive understanding of the consequent climate responses can be obtained. The high-temporal-resolution estimated emissions can reveal features of emissions that are absent from the traditional ones with coarse temporal resolutions. Moreover, the timely and dynamics-based emission estimates with high temporal resolutions are essential for regional air quality modelling, especially for the occurrence of severe haze pollution, associated with timely evaluations of the impact on public health (Attri et al., 2001; Wang et al., 2014; Ji et al., 2018; Wang et al., 2020; Liu et al., 2021) and events that lead to large changes in emissions and significant socioeconomic impacts such as the COVID-19 pandemic (Huang et al., 2020; Le et al., 2020).

2 Data assimilation and experimental design

The estimate of PM_2.5 emissions can be successfully constrained by the PM_2.5 concentration observations through an ensemble Kalman filter (EnKF; Peng et al., 2017, 2018, 2020). For a retrospective reanalysis mode, here, all available PM_2.5 concentration observations, including those data collected after the analysis time, can be used. Thus, an EnKS, a direct generalization of the EnKF, is applied to incorporate PM_2.5 concentration observations both before and after the analysis time, aiming to provide an optimal estimate of the PM_2.5 emission. In simple words, the emissions are updated by current and future observations though EnKS, while the concentrations are updated by current observations though EnKF. Detailed procedures of the EnKS are described in Sect. 2.1.

2.1 An ensemble Kalman smoother to update the source emission

The ensemble priors of source emissions e^f is created by multiplying a scaling factor λ^f by the prescribed emission e^p (Peng et al., 2017, 2018, 2020), where the superscript f denotes priors. Given a time-invariant e^p, the update of e^f is equivalent to the update of λ^f. Due to a time lag, the prior scaling factor at time t−1 (${\mathit{\lambda }}_{t-\mathrm{1}}^{\mathrm{f}}$) is updated by chemical observations at time t (${\mathbf{y}}_t^{\mathrm{c}}$). At time t−1, the prior scaling factor for the ith member is written as

$$\begin{array}{}\text{(1)}& {\mathit{\lambda }}_{i,t-\mathrm{1}}^{\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{M}}\left[\left(\mathit{\beta }{\displaystyle \frac{{\mathbf{c}}_{i,t-\mathrm{1}}^{\mathrm{f}}}{\stackrel{\mathrm{‾}}{c_{t-\mathrm{1}}^{\mathrm{f}}}}}+\mathrm{1}-\mathit{\beta }\right)+\sum _{j=t-M}^{t-\mathrm{2}}{\mathit{\lambda }}_{i,j\mathrm{|}j+\mathrm{1}:t-\mathrm{1}}^a\right].\end{array}$$

The first term is the concentration ratio given by the prior of the chemical fields (${\mathbf{c}}_{i,t-\mathrm{1}}^{\mathrm{f}}$) normalized by the ensemble mean ($\stackrel{\mathrm{‾}}{c_{t-\mathrm{1}}^{\mathrm{f}}}$), where β is an inflation factor used to compensate for the insufficient ensemble spread (Peng et al., 2017). Through using the concentration ratio, each ensemble member of the source emissions naturally has the spatial correlations given by the chemical fields. The second term is the mean of the posterior scaling factors at previous assimilation cycles, where the superscript a denotes posteriors, M is the length of smoothing, and the subscript $j+\mathrm{1}:t-\mathrm{1}$ indicates that the scaling factor at time j is updated by future observations from j+1 to t−1. The assimilation of future observations will be described below.

The ensemble square-root filter (EnSRF) (Peng et al., 2017) is used to update ${\mathit{\lambda }}_{t-\mathrm{1}}^{\mathrm{f}}$ by assimilating ${\mathbf{y}}_t^{\mathrm{c}}$. For the scaling factor at time t−1, the posterior ensemble mean is given by

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{‾}}{{\mathit{\lambda }}_{t-\mathrm{1}}^a}=\stackrel{\mathrm{‾}}{{\mathit{\lambda }}_{t-\mathrm{1}}^{\mathrm{f}}}+\mathit{\rho }\circ {\mathbf{P}}_{t-\mathrm{1},t}^{\mathrm{ec}}{\mathbf{H}}_t^{c\mathbf{T}}{\left({\mathbf{H}}_t^{\mathrm{c}}{\mathbf{P}}_t^{\mathrm{c}}{\mathbf{H}}_t^{c\mathbf{T}}+{\mathbf{R}}_t^{\mathrm{c}}\right)}^{-\mathrm{1}}({\mathbf{y}}_t^{\mathrm{c}}-H_t^{\mathrm{c}}\stackrel{\mathrm{‾}}{c_t^{\mathrm{f}}}),\end{array}$$

and the posterior ensemble perturbations are given by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{{\mathit{\lambda }}^{\prime }}_{i,t-\mathrm{1}}^a& ={{\mathit{\lambda }}^{\prime }}_{i,t-\mathrm{1}}^{\mathrm{f}}-\mathit{\rho }\circ {\mathbf{P}}_{t-\mathrm{1},t}^{\mathrm{ec}}{\mathbf{H}}_t^{c\mathbf{T}}\\ & {\left[{\left(\sqrt{{\mathbf{H}}_t^{\mathrm{c}}{\mathbf{P}}_t^{\mathrm{c}}{\mathbf{H}}_t^{c\mathbf{T}}+{\mathbf{R}}_t^{\mathrm{c}}}\right)}^{-\mathrm{1}}\right]}^{\mathrm{T}}\\ & {\left[\sqrt{{\mathbf{H}}_t^{\mathrm{c}}{\mathbf{P}}_t^{\mathrm{c}}{\mathbf{H}}_t^{c\mathbf{T}}+{\mathbf{R}}_t^{\mathrm{c}}}+{\sqrt{\mathbf{R}}}_t^{\mathrm{c}}\right]}^{-\mathrm{1}}{\mathbf{H}}_t^{\mathrm{c}}{{\mathit{\lambda }}^{\prime }}_{i,t-\mathrm{1}}^{\mathrm{f}},\end{array}\end{array}$$

where ${\mathbf{P}}_{t-\mathrm{1},t}^{\mathrm{ec}}$ denotes the background error covariance matrix of ${\mathit{\lambda }}_{t-\mathrm{1}}^{\mathrm{f}}$ and ${\mathbf{c}}_t^{\mathrm{f}}$; ${\mathbf{P}}_t^{\mathrm{c}}$ indicates the background error covariance matrix of ${\mathbf{c}}_t^{\mathrm{f}}$; $H_t^{\mathrm{c}}$, ${\mathbf{H}}_t^{\mathrm{c}}$, and ${\mathbf{R}}_t^{\mathrm{c}}$ are the observation forward operator, Jacobian matrix, and observation error covariance matrix of the chemical fields at time t; ρ is the localization matrix; and ∘ denotes the Schur (element-wise) product.

By applying the ensemble Kalman smoother (EnKS) (Whitaker et al., 2002; Peters et al., 2007), the chemical observation ${\mathbf{y}}_t^{\mathrm{c}}$ is also assimilated to update the posterior scaling factor at previous assimilation cycles $j(j=t-K,\mathrm{\cdots },t-\mathrm{2})$. After assimilating the future chemical observation at time t, the posterior ensemble mean of the scaling factor at j is given by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{{\mathit{\lambda }}_{j\mathrm{|}j+\mathrm{1}:t}^a}& =\stackrel{\mathrm{‾}}{{\mathit{\lambda }}_{j\mathrm{|}j+\mathrm{1}:t-\mathrm{1}}^a}+\mathit{\rho }\circ {\mathbf{P}}_{j\mathrm{|}j+\mathrm{1}:t-\mathrm{1},t}^{\mathrm{ec}}{\mathbf{H}}_t^{c\mathbf{T}}\\ & {\left({\mathbf{H}}_t^{\mathrm{c}}{\mathbf{P}}_t^{\mathrm{c}}{\mathbf{H}}_t^{c\mathbf{T}}+{\mathbf{R}}_t^{\mathrm{c}}\right)}^{-\mathrm{1}}({\mathbf{y}}_t^{\mathrm{c}}-H_t^{\mathrm{c}}\stackrel{\mathrm{‾}}{c_t^{\mathrm{f}}}),\end{array}\end{array}$$

and the posterior ensemble perturbations are given by

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{{\mathit{\lambda }}^{\prime }}_{i,j\mathrm{|}j+\mathrm{1}:t}^a& ={{\mathit{\lambda }}^{\prime }}_{i,j\mathrm{|}j+\mathrm{1}:t-\mathrm{1}}^a-\mathit{\rho }\circ {\mathbf{P}}_{j\mathrm{|}j+\mathrm{1}:t-\mathrm{1},t}^{\mathrm{ec}}{\mathbf{H}}_t^{c\mathbf{T}}\\ & {\left[{\left(\sqrt{{\mathbf{H}}_t^{\mathrm{c}}{\mathbf{P}}_t^{\mathrm{c}}{\mathbf{H}}_t^{c\mathbf{T}}+{\mathbf{R}}_t^{\mathrm{c}}}\right)}^{-\mathrm{1}}\right]}^{\mathrm{T}}\\ & {\left[\sqrt{{\mathbf{H}}_t^{\mathrm{c}}{\mathbf{P}}_t^{\mathrm{c}}{\mathbf{H}}_t^{c\mathbf{T}}+{\mathbf{R}}_t^{\mathrm{c}}}+{\sqrt{\mathbf{R}}}_t^{\mathrm{c}}\right]}^{-\mathrm{1}}{\mathbf{H}}_t^{\mathrm{c}}{{\mathit{\lambda }}^{\prime }}_{i,t-\mathrm{1}}^{\mathrm{f}},\end{array}\end{array}$$

where ${\mathbf{P}}_{j\mathrm{|}j+\mathrm{1}:t-\mathrm{1},t}^{\mathrm{ec}}$ denotes the background error covariance matrix of ${\mathit{\lambda }}_{j\mathrm{|}j+\mathrm{1}:t-\mathrm{1}}^a$ and ${\mathbf{c}}_t^{\mathrm{f}}$. After Eqs. (2)–(5), the updated ${\mathit{\lambda }}_{j\mathrm{|}j+\mathrm{1}:t}^a,j(j=t-M+\mathrm{1},\mathrm{\cdots },t-\mathrm{1})$ will be used to construct the prior scaling factor at the next time t+1.

As a Monte Carlo approach, the EnKS uses the forecast analysis error covariances based on ensemble forecasts and/or analyses to compute the Kalman gain matrix with time lags in order to incorporate observations from the past to the future. The first iteration of the EnKS is equivalent to the EnKF that assimilates observations up to the analysis time. The following iterations of EnKS assimilate observations in the future to update the state at the analysis time. The hourly forecasts of PM_2.5 concentration from the cycling assimilation experiment matched the independent observed quantities (Fig. 1). Therefore, the ability of EnKS to retrieve the source emissions has been demonstrated. Previous studies also showed that simulations forced by the posterior emissions could produce improved forecasts for PM_2.5, SO₂, and NO₂ compared to those produced with a priori emissions (Peng et al., 2020).

Figure 1 Times series of hourly PM_2.5 concentration biases (µg m⁻³). The ensemble mean priors compared to the observed quantities for December of the years 2016–2020 (grey and black) and the mean biases of the years 2016–2020 (blue).

2.2 WRF-Chem model, observations, and emissions

To simulate the transport of aerosol and chemical species, the WRF-Chem model version 3.6.1 (Grell et al., 2005) that has the meteorological and chemical components fully coupled is used. The model parameterization schemes follow Peng et al. (2017). Figure 2 shows the model domain that covers most east Asian regions. The horizontal grid spacing is 45 km, with 57 vertical levels and the model top at 10 hPa.

Figure 2 Model domain and observation sites for cycling assimilation. Red and blue dots denote the assimilated and unassimilated observational sites, respectively.

Experiments are conducted for each year from 2016 to 2020 separately. The 6 h meteorological observations, including all in situ observations and cloud motion vectors from the National Centers for Environmental Prediction (NCEP) Global Data Assimilation System (GDAS; http://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_2.htm, last access: 6 November 2023), are assimilated every 6 h. The hourly observed chemical quantities, which contain PM₁₀, PM_2.5, SO₂, NO₂, O₃, and CO from the Ministry of Ecology and Environment of China (https://aqicn.org/map/china/cn/, last access: 6 November 2023), are assimilated every hour. Figure 2 shows the assimilated chemical observation network, which has 560 randomly chosen stations from 1576 stations in total. The thinning of observations is applied to avoid correlated errors of observations. The spatial autocorrelation of the thinning of observations is close to the original observations (Peng et al., 2017). The observation priors are computed by the observer portion of the Grid-point Statistical Interpolation system (GSI) (Kleist et al., 2009).

The hourly and time-invariantly prescribed anthropogenic emissions are obtained from the EDGAR-HTAP (Emission Database for Global Atmospheric Research for Hemispheric Transport of Air Pollution) v2.2 inventory (Janssens-Maenhout et al., 2015b), in which the Chinese emissions are derived from the Multi-resolution Emission Inventory (MEIC) in 2010 (Lei et al., 2011; Li et al., 2014). Natural emissions, including the biogenic (Guenther et al., 1995), dust (Ginoux et al., 2001), dimethyl sulfide, and sea salt emissions (Chin et al., 2000), are computed online.

2.3 Assimilation and ensemble configurations

The PM_2.5 emission directly gives the primary PM_2.5 and then the primary PM_2.5 along with other precursor emissions that could contribute to the secondary PM_2.5. The observations of PM_2.5 concentrations that contain both primary and secondary PM_2.5 are used to constrain the PM_2.5 emission through data assimilation. Thus, the correlations between the concentration observations and source emissions might be contaminated by the secondary PM_2.5. Since the secondary formation process can be captured by the WRF-Chem model, the impact of the secondary PM_2.5 is indirectly considered. The detailed updated state variables with the according observations follow Peng et al. (2018). The concentrations and emissions of PM_2.5 and PM_2.5 precursors (SO₂ and NO) that have observations are updated by the observed quantities. Besides, NH₃ concentrations and emissions are constrained by PM_2.5 observations; however, the VOCs that are also PM_2.5 precursors are not updated due to the lack of direct and limited observations. One possible way to untangle the impact of secondary PM_2.5 on the estimates of PM_2.5 emissions is to jointly estimate the source emission and the primary and secondary PM_2.5 given the concentration observations.

Figure 3 (a) Dynamics-based monthly PM_2.5 emission estimates (Tg d⁻¹) summed over mainland China for each year from 2016 to 2020 (coloured) and the estimated PM_2.5 emission from MEIC (grey). (b) Ratio of PM_2.5 emission changes between two adjacent years from year 2016 to 2020 normalized by the PM_2.5 emissions of the year 2016. (c) Monthly fractions of dynamics-based PM_2.5 emission estimates for years 2016–2020 (light blue); the 5-year mean fractions of dynamics-based monthly PM_2.5 emission estimates, with bars denoting 1 standard deviation of the 5-year variations (dark blue); and the monthly fractions of estimated PM_2.5 emissions from MEIC (grey).

Figure 4 (a) Spatial distribution of dynamics-based PM_2.5 emission estimates (µg m⁻² s⁻¹) for the year 2016 and spatial distributions of dynamics-based PM_2.5 emission changes for the years (b) 2017, (c) 2018, (d) 2019, and (e) 2020 compared to those of the year 2016.

The National Oceanic and Atmospheric Administration (NOAA) operational EnKF system (https://dtcenter.ucar.edu/com-GSI/users/docs/users_guide/GSIUserGuide_v3.7.pdf, last access: 6 November 2023), which is an EnSRF and is modified with the EnKS feature, is used to assimilate the observations. The ensemble size is set to 50. To combat the sampling error resulting from a limited ensemble size, covariance localization and inflation are applied. The Gaspari and Cohn (GC) (1999) function with a length scale of 675 km is used to localize the impact of observations and mitigate the spurious error correlations between observations and state variables. The constant multiplicative posterior inflation (Whitaker and Hamill, 2012), with a coefficient of 1.12 for all meteorological and chemical variables, is applied to enlarge the ensemble spread. The inflation β for advancing the scale factor is 1.2. The smoothing length M for source emissions is 4, and the EnKS lagged length K is 6. The larger the K value, the more future observations are assimilated to constrain the current emission estimate. But the sample estimated temporal correlations could be contaminated by sampling errors and model errors, especially with increased lagged times. Thus, there is a tradeoff between the amount of future observations and the accuracy of the sample estimated temporal correlations. The choice of K (=6) is determined by sensitivity experiments.

At 00:00 UTC on 26 December of the previous year, ensemble initial conditions (ICs) of the meteorological fields are generated by adding random perturbations that sample the static background error covariances (Barker et al., 2012) on the NCEP FNL (final) analyses (Torn et al., 2006). Ensemble ICs of the chemical fields are 0, and source emissions of each ensemble member are adopted from the EDGAR-HTAP v2.2 inventory with random perturbations of mean 0 and variances of 10 % of the emission values. Hourly ensemble lateral boundary conditions (LBCs) are generated using the same fixed-covariance perturbation technique as the ensemble ICs. After a 6 d spin up, ensemble data assimilation experiments start cycling for each year.

Figure 5 (a) The differences in dynamics-based PM_2.5 emission estimates between the years 2017–2020 and 2016 and (b) the differences normalized by the of year 2016.

3 PM_2.5 emission for years 2016–2020

Starting from the time-invariant source emission PR2010 (Janssens-Maenhout et al., 2015b), the dynamics-based estimates of the PM_2.5 emissions are obtained, which include the contributions of both the anthropogenic and biomass burning emissions. The mean annual PM_2.5 emissions from biomass burning in China (2003–2017) amounted to 0.51 Tg (Yin et al., 2019). The annual dynamics-based estimates of PM_2.5 emissions (DEPEs) averaged over mainland China for the years 2016–2020 without biomass burning emissions are 7.66, 7.40, 7.02, 6.62, and 6.38 Tg, respectively. The values from the Multi-resolution Emission Inventory (MEIC; Zheng et al., 2018), which does not consider the contributions of biomass burning emissions, are 8.10, 7.60, 6.70, 6.38, and 6.04 Tg. Thus the annual DEPEs are very close to the values of MEIC. From 2017 to 2020, the estimated annual PM_2.5 emissions are reduced by 3.4 %, 8.4 %, 13.6 %, and 16.7 % compared to that of year 2016. There has been a 3 %–5 % persistent reduction in annual PM_2.5 emissions from the year 2017 to 2020, which demonstrates the effectiveness of China's Clean Air Action (China State Council, 2013), implemented since 2013, and China's Blue Sky Defense War Plan (2018), enforced since 2018 with strengthened industrial emission standards, phased out outdated industrial capacities, promoted clean fuels in residential sector, and so on (Zhang et al., 2019).

Figure 6 Spatial distributions of dynamics-based PM_2.5 emission changes in (a) December 2017 compared to November 2017 and (b) December 2018 compared to November 2018.

Table 1 Dynamics-based PM_2.5 emission estimates of the year 2016 for each province whose value is larger than 0.01 µg m⁻² s⁻¹ are shown in the second column. Ratios of PM_2.5 emission changes for the years 2017–2020 compared to 2016 are shown from the third to the sixth column, with negative (roman)/positive (italic) values indicating a decrease/increase in PM_2.5 emissions.

Table 2 The 5-year mean diurnal fractions (%) of the dynamics-based PM_2.5 emission estimates over mainland China in local solar time (LST) for each month.

The monthly DEPEs show reductions in PM_2.5 emissions nearly every month from 2016 to 2020 (Fig. 3a), which further demonstrates the effectiveness of China's national plan. Compared to 2016, both the reduction amount and the reduction ratio of PM_2.5 emission are more prominent for February, March, June–September, and November than the other months (Fig. 3b). Given larger magnitudes of PM_2.5 emissions in winter than in summer, emission controls with a focus on October to May should be considered in the design of future clean-air actions in China since total PM_2.5 emissions during period account for approximately 75 % of the annual amount. Spatial distributions of the changes in PM_2.5 emissions from 2017 to 2020 compared to 2016 show that significant decreases occurred in the Beijing–Tianjin–Hebei region (BTH), the Yangtze River Delta region (YRD), the Pearl River Delta region (PRD), and the Sichuan–Chongqing Region (SCR), especially for the years 2019–2020 (Fig. 4). From 2016 to 2020, BTH, YRD, and SRC had larger reductions in PM_2.5 emissions than PRD, but SCR had a larger reduction ratio compared to 2016 than BTH and YRD (Fig. 5). Therefore, BTH and YRD have more potential for PM_2.5 emission controls than PRD and SCR, which can give guidance for future clean-air actions. More specifically, most provinces show PM_2.5 emission reductions from 2016 to 2020, and the reduction ratios generally increase from year 2017 to 2020 (Table 1), which confirms continuous and effective emission controls from Clean Air Action and the Blue Sky Defense War Plan in China. The monthly DEPE also demonstrates the effectiveness of strict implementations of emission reduction policies in China, such as the coal ban for residential heating since the 2017–2018 winter. There was a sharp change in PM_2.5 emissions, from an increase in 2017 to a decrease in 2018. As shown by Fig. 6, spatial distributions of the changes in PM_2.5 emissions in December compared to November in 2017 show obvious increases in most of China. However, the changes in 2018 show significant decreases in the areas of the Beijing, Tianjin, Hebei, Shanxi, Henan, and Anhui provinces due to the implementation of the coal ban.

Despite the trend in PM_2.5 emissions from 2016 to 2020, the DEPE of the year 2016 has similar monthly distributions compared to MEIC2016–2020 in general (Fig. 3a). MEIC has a pan-shaped monthly distribution, with nearly time-invariant PM_2.5 emissions from April to October. This seasonal dependence of emissions is mainly contributed by the variations in residential energy use, which are empirically dependent on coarse monthly mean temperature intervals and thus cannot reflect the realistic monthly variations (Streets et al., 2003; Li et al., 2017). The centralized heating system in northern China has fixed dates for being turned on and turned off during each heating season. Therefore, a sudden rise in emissions from October to November and a sudden drop in emissions from March to April are shown. But the turning-on and turning-off dates are variable in different regions, which imposes a smoothing effect on the emissions. However, the DEPE shows a V-shaped monthly distribution, with the minimum occurring in August. The estimated PM_2.5 emission is 11.8 % higher than that of MEIC2016 in April but 12.1 % lower than that of MEIC2016 in August, and these different monthly distributions can influence the consequent climate responses, including the radiative forcing and energy budget (Yang et al., 2020), and can also impact on health issues (Liu et al., 2018). Moreover, monthly fractions of the DEPE are consistent across years (Fig. 3c). The absence of interannual variations of monthly PM_2.5 emission fractions provides a basis for previous studies that follow the same monthly changes in source emissions from different years (Zhang et al., 2009; Zheng et al., 2020, 2021). Monthly allocations of PM_2.5 emissions can be directly and objectively obtained given an estimated total annual amount based on the estimated monthly fractions of DEPE, which is valuable for emission inventories; air quality simulations; and, potentially, applications for future scenarios due to more accurate month fractions of DEPE. Since the hourly priors of PM_2.5 concentrations from the cycling assimilation for optimally estimating PM_2.5 emissions fit to the observed PM_2.5 quantities (Fig. 1), the monthly DEPE provides more realistic monthly fluctuations than the empirical estimate.

4 Diurnal variations of PM_2.5 emissions

The DEPEs with high temporal resolutions given the time-invariant prior PR2010 can reveal features that are unable to be represented in the commonly used emission estimates. Although the prior PR2010 has no diurnal variations, hourly posteriors of PM_2.5 emissions provide the first objectively estimated diurnal variations for different seasons for the years 2016–2020. However, these estimated diurnal variations include the contributions of the time-varying boundary layer. An observing system simulation experiment (OSSE) is performed to investigate the effects of the boundary layer from 00:00 UTC on 29 December 2015 to 00:06 UTC on 1 February 2016. Details of this OSSE are presented in the Supplement. The results indicate that the magnitude of posterior PM_2.5 emissions from the OSSE is closer to the true emission than the prior. Since we have hourly assimilated observations to simultaneously update the chemical concentrations and source emissions, the impacts of the time-varying boundary layer on the posterior PM_2.5 emissions are limited (Fig. S1 in the Supplement). Slightly larger estimated PM_2.5 emission fractions occurred in the morning, and smaller estimated PM_2.5 emission fractions occurred in the afternoon in comparison to the time-invariant true emission. Nevertheless, the influences of the time-varying boundary layer are still important to PM_2.5 emission estimates. To statistically present the diurnal variations, the fractions of hourly PM_2.5 emissions divided by the daily amount are averaged over different years and regions after excluding the impacts of the time-varying boundary layer based on the short-term period simulation, although the influences of the boundary layer could strongly vary with seasons or years (Figs. 7 and 8 and Table 2). The diurnal variations of PM_2.5 emissions are critical for understanding the mechanisms of PM_2.5 formation and evolution and are also essential for PM_2.5 simulations and forecasts.

Figure 7 The 5-year mean diurnal variations of dynamics-based PM_2.5 emission fractions averaged over (a) mainland China, (b) megacities with an urban population ≥ 5 million, and (c) non-megacities with an urban population < 5 million.

Figure 8 Diurnal variations of dynamics-based PM_2.5 emission fractions for the years 2016–2020 (light blue), and the 5-year mean fractions, with bars denoting 1 standard deviation of the 5-year variations (dark blue), are averaged over mainland China for (a) January, (b) April, (c) July, and (d) October.

The 5-year mean diurnal variations of the estimated PM_2.5 emission fraction for mainland China show that, despite the monthly variations of PM_2.5 emissions, the diurnal-variation fractions for November, December, January, and February are similar, while those for June, July, and August are similar (Fig. 7a). There are stronger diurnal variations of PM_2.5 emissions in summer than in winter, which are represented by larger PM_2.5 emission fractions during the morning and fewer PM_2.5 emission fractions during evening. The diurnal variations of PM_2.5 emissions from March to May gradually transform from the patterns of winter to those of summer and vice versa for the diurnal variations of PM_2.5 emissions from September to November. The monthly changes in diurnal variations of PM_2.5 emissions are consistent with the seasonal dependence since monthly variations of PM_2.5 emissions are mainly related to the variations of residential consumptions (Li et al., 2017) in which the space heating has nearly no diurnal variations, and then larger PM_2.5 emissions during winter lead to reduced diurnal variations compared to in summer. Similarly to the monthly fractions of estimated PM_2.5 emissions for mainland China, diurnal variations of PM_2.5 emission fractions are consistent across years for a given month (Fig. 8). Table 2 gives 5-year mean diurnal variations of the estimated PM_2.5 emission fractions for each month. Based on these high-resolution diurnal-variation fractions, hourly estimates of PM_2.5 emissions can be objectively obtained for a given monthly estimated PM_2.5 emission.

Figure 9 Diurnal variations of PM_2.5 emission fractions for each month based on diurnal variation profiles from the US and EU (Wang et al., 2010).

Figure 10 Hourly (light red and blue) and daily (dark red and blue) dynamics-based PM_2.5 emission estimates (kg h⁻¹) summed over mainland China from January to March for the years 2019 and 2020.

Despite the high temporal resolution, the DEPE also has the ability to analyse diurnal variations for specific cities. The monthly changes in the diurnal variations of PM_2.5 emissions estimated for megacities with urban populations larger than 5 million and non-megacities with urban populations smaller than 5 million (China State Council, 2013) are consistent with those estimated from mainland China (Fig. 7). Compared to the diurnal variations of PM_2.5 emissions estimated for mainland China, the megacities have stronger diurnal variations, while the non-megacities have weaker diurnal variations. These detailed descriptions of PM_2.5 emissions that are usually absent in common emission estimates can be essential for PM_2.5 simulation, especially for providing timely and realistic guidance for severe haze events.

Figure 11 Spatial distributions of dynamics-based PM_2.5 emission estimates (µg m⁻² s⁻¹) in (a) February and (c) March of the year 2019 and spatial distributions of dynamics-based PM_2.5 emission reductions of the year 2020 compared to the year 2019 for (b) February and (d) March.

There has been a lack of local measurements for diurnal variations and widely adopted diurnal variation profiles of PM_2.5 emissions in China. Compared to the diurnal variations of PM_2.5 emission fractions estimated based on diurnal variation profiles from the US and EU (Wang et al., 2010; Du et al., 2020), the estimated PM_2.5 emission fractions are 1.25 % larger during the evening, which greatly changes the diurnal variations of DEPE. The noon and evening peaks estimated from DEPE have smaller PM_2.5 emission fractions, with mean underestimations of PM_2.5 emission fractions of 0.40 % and 0.83 % for noon peak and evening peak, respectively (Figs. 7a and 9). In fact, the smaller evening peaks of Wang et al. (2010) occurred in November, December, January, February, and March, while they are almost indistinct from April to October, similarly to those from DEPE. The morning peak of Wang et al. (2010) is similar to that of DEPE for spring and fall, but the former overestimates PM_2.5 emission fractions by 0.57 % for winter, while it underestimates PM_2.5 emission fractions by 1.05 % for summer. Due to the overestimated peaks, diurnal variations of Wang et al. (2010) have a sharper appearance rate for the morning peak and a disappearance rate for the evening peak. Compared to the diurnal variations based on diurnal variation profiles from the US and EU (Wang et al., 2010), the diurnal variations of the DEPE are constrained by the atmospheric–chemical model and observed PM_2.5 concentrations, which can objectively determine the diurnal variations of PM_2.5 emissions for specific regions and seasons.

Figure 12 Mean spatial distributions of PM_2.5 emission differences (µg m⁻² s⁻¹) between the years 2020 and 2019 for 9 weeks starting on 9 January 2020. Negative (positive) values indicate that the PM_2.5 emissions of the year 2020 are smaller (larger) than those of the year 2019. The numbers in (a) denote provinces as follows: (1) Heilongjiang, (2) Neimenggu, (3) Xinjiang, (4) Jilin, (5) Liaoning, (6) Gansu, (7) Hebei, (8) Beijing, (9) Shanxi, (10) Tianjin, (11) Shanxi, (12) Ningxia, (13) Qinghai, (14) Shandong, (15) Xizang, (16) Henan, (17) Jiangsu, (18) Anhui, (19) Sichuan, (20) Hubei, (21) Chongqing, (22) Shanghai, (23) Zhejiang, (24) Hunan, (25) Jiangxi, (26) Yunnan, (27) Guizhou, (28) Fujian, (29) Guangxi, (30) Guangdong, (31) Taiwan, (32) Hong Kong, (33) Macao, (34) Hainan.

5 Impact of COVID-19 on PM_2.5 emissions

The abrupt changes in PM_2.5 emissions during the initial stage of COVID-19 in China provide a natural case study to validate the ability of the dynamic-based data assimilation method to obtain high-temporal-resolution PM_2.5 emission estimates. The abrupt outbreak of the COVID-19 pandemic has produced dramatic socioeconomic impacts in China. To prevent the virus from spreading, a lockdown was first implemented on 23 January 2020 in Wuhan, Hubei Province, and subsequently, the national lockdown was enforced in China (Liu et al., 2020; Huang et al., 2020; Zhu et al., 2021). Consequently, the total PM_2.5 emissions of February 2020 for China show an obvious decrease compared to those of previous years (Fig. 3). The high=temporal-resolution DEPE reveals the detailed changes of PM_2.5 emissions with time (Fig. 10). The PM_2.5 emissions started to decrease right around the COVID outbreak and had been smaller than those of the year 2019 till early March. The emissions in the following months of 2020 were similar to those of 2019 due to the epidemic prevention and control policies enforced by the Chinese government. During February 2020, the DEPE showed significant reductions at the North China Plain and northeast of China, where prominent PM_2.5 emissions occurred, while spotted PM_2.5 emission differences with small magnitudes were shown in the other regions (Fig. 11a–b). Along with recovery from COVID-19, the estimated PM_2.5 emissions rebounded in March (Figs. 3a, 10, 11c–d), which is ascribed to the national work resumption. Thus, the DEPE is able to timeously reflect the dynamic response of PM_2.5 emissions to COVID-19. Although similar emission reductions and emission trends are obtained from the bottom-up technique (Zheng et al., 2021), the reduction amount and ratio from the bottom-up technique are larger than those estimated from DEPE (Fig. 10 and Table 1). This is possibly due to significant reductions in PM_2.5 emissions from the residential sector, as in the bottom-up technique (Zheng et al., 2021); however, PM_2.5 emissions from the residential sector might not have significantly changed around the COVID outbreak.

To avoid fluctuations due to diurnal variations and monthly changes in PM_2.5 emissions, 7 d averaged PM_2.5 emission differences between the years 2020 and 2019 are used to analyse the dynamic impact of COVID-19 on PM_2.5 emissions (Fig. 12). Before the lockdown, there were slight PM_2.5 emission differences over several provinces (Fig. 12a–b). During the first week of lockdown, PM_2.5 emission reductions larger than $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{2}}$ (µg m⁻² s⁻¹) – that is about 60 %–70 % emission reductions – occurred in the Hubei, Hunan, Guangdong, Anhui, and Zhejiang provinces (Fig. 12c). The PM_2.5 emission reduction extended to BTH and the Shandong Province during the second week of lockdown (Fig. 12d) and continuously spread to the three northeast provinces of China during the third week of lockdown (Fig. 12e). During the third week of lockdown, the increased PM_2.5 emissions for BTH and SCR were possibly caused by the long national vacation of the spring holiday of the year 2019 (Ji et al., 2018). The inhomogeneous spatial variations of PM_2.5 emissions possibly relate to different traditions and policy enforcements for different provinces. The PM_2.5 emission reduction had been maintained over the central and northern parts of China till early March when the lockdown was lifted (Fig. 12f–i). Though it is hard to see continuous and consistent signals of lockdown for the whole China, the timely DEPE can provide up-to-date guidance for quantifying the socioeconomic impacts of rare events with large emission changes such as the COVID-19 pandemic.

6 Discussion

High-temporal-resolution and dynamics-based estimations of PM_2.5 emissions can be objectively and optimally obtained by assimilating past and future observed surface PM_2.5 concentrations through flow-dependent error statistics. This advanced assimilation strategy can be applied for emission estimates of other chemical species when corresponding observations are available and can extend to observation types besides the surface concentrations, like the aerosol optical depth (Liu et al., 2011; Choi et al., 2020). Moreover, current estimates of PM_2.5 emissions are lacking in terms of explicit representations of primary and secondary PM_2.5, which could be resolved by joint estimation of the source emissions and primary and secondary PM_2.5 given the concentration observations. Another deficiency of this top-down technique is that it cannot directly determine dynamics-based PM_2.5 emissions for different sectors and contributions from different policies, although the bottom-up technique has the potential to untangle the different contributions from different policies and quantify the different impacts on different sectors. However, this top-down technique can be integrated into the bottom-up technique to retain the advantages of both methods. One future work is to integrate the top-down technique with the bottom-up one, through which the emission estimates for different sectors and polices could be quantified. The annual emission estimates from the bottom-up technique can be further downscaled to hourly estimates by first distributing the annual amount to each month through the monthly allocations estimated from the top-down technique, then assuming evenly daily distribution and finally applying the fractions of diurnal variation estimated from the top-down technique. The information collected by the bottom-up technique is retained, while the common drawback of a coarse temporal resolution for the bottom-up technique is remedied. The integrated bottom-up and top-down techniques can improve spatiotemporal representations of source emissions across timescales and sectors, which is beneficial for emission inventories, air quality forecasts, regulation policies, and emission trading schemes.