Impact of the ‘Coal-to-Natural Gas’ Policy on Criteria Air Pollutants in Northern China

During the last decades, China had issued a series of stringent control measures, resulting in a large decline in air pollutant concentrations. To quantify the net change in air pollutant concentrations driven by emissions, we developed an approach of determining the closed interval of the deweathered percentage change (DPC) in the concentration of air pollutants on an annual scale, as well as the closed intervals of cumulative DPC in a year compared with that in the base year. Thus, the hourly mean mass concentrations of criteria air pollutants to determine their interannual variations and the closed intervals of their DPCs during the heating seasons from 2013 to 2019 in Qingdao (a coastal megacity) were analyzed. The seasonal mean SO2 concentration decreased from 2013 to 2019. The seasonal mean CO, NO2, and PM2.5 concentrations also generally decreased from 2013 to 2017, but increased unexpectedly in 2018 (from 0.9 mg m−3 (CO), 42 μg m−3 (NO2), and 51 μg m−3 (PM2.5) in 2017 to 1.1 mg m−3, 48 μg m−3, and 64 μg m−3 in 2018, respectively). The closed intervals of DPC in concentrations of CO, NO2, and PM2.5 from the 2017 heating season (2017/2018) to the 2018 heating season (2018/2019) were obtained at (27%, 30%), (15%, 18%), and (30%, 33%), respectively. Such high positive endpoint values of the closed intervals, in contrast to their small interval lengths, indicate increased emissions of these pollutants and/or their precursors in 2018/2019 compared with 2017/2018, by minimizing the meteorological influences. The rebounds of CO, NO2, and PM2.5 in 2018/2019 were likely associated with a doubled increase in natural gas (NG) consumption implemented by the “coal-to-NG” project, as the total energy consumption showed little difference. Our results suggested an important role of the “coal-to-NG” project in driving concentrations of air pollutant increases in China in 2018/2019, which need integrated assessments.


Introduction
China has experienced rapid economic growth and urban expansion in the last four decades, leading to excessive emissions of multiple air pollutants [1][2][3][4][5]. Focus has been placed on severe haze pollution in winter, which is characterized by high mass concentrations of atmospheric particles with an aerodynamic diameter of ≤2.5 µm (PM 2.5 ) or ≤10 µm (PM 10 ) [1,3,6]. PM exposure is well-known to cause several adverse health outcomes, such as acute and chronic respiratory illnesses and cardiovascular diseases, or even death, and ambient PM 2.5 human-health effects have been widely investigated in China [7][8][9]. For example, it was revealed that a 2.9% increase in household healthcare expenditure was estimated to be attributed to a 1% increase in yearly PM 2.5 exposure in China [10].
To relieve PM 2.5 pollution, China has issued a number of stringent control measures since 2013 under the "Action Plan on Prevention and Control of Air Pollution" [11], including the enforcement of flue gas desulfurization and denitrification, elimination of in-use vehicles exhibiting high emissions, implementation of ultralow air pollutant emission standards, and the implementation of the "coal-to-natural gas (NG)" policy in 2018 across northern China [12,13]. With the implementation of the emission control policies, the overall air quality in China has been significantly improved [2]. The concentrations of PM 2.5 , PM 10 , and SO 2 were noticeably decreased, but the still high-level PM and stagnant NO 2 present further challenges [2,11]. Ambient air quality standards in China (GB3095-2012) include six criteria air pollutants: PM 2.5 , PM 10  NG is commonly considered as a clean fossil fuel, resulting in fewer emissions of air pollutants including primary carbonaceous particles, SO 2 , volatile and semi-volatile organic species, etc., than burning coal to produce an equal amount of energy [14,15]. The implementation of the "coal-to-NG" policy is expected to directly reduce air pollutant emissions, and thereby relieve PM 2.5 air pollution in northern China. Replacing coal with NG for industrial facilities would reduce emissions of SO 2 and PM 2.5 ; however, its impact on NO x and CO emissions reduction requires further investigation. PM 2.5 concentrations fluctuated largely with the banned use of solid fuels in Krakow, Poland during the last decades, although a long-term decreasing trend was there [16]. The impact of the "coal-to-NG" program on air quality, especially on criteria pollutants, is contentious; for example, increased concentrations of NO x and particulate nitrate were reported in Beijing and Tianjin in China in 2018 [12,13]. The impact of this policy on air pollution has not been thoroughly investigated in other parts of northern China because of the lack of air pollutant emission data associated with the program. However, the big challenges associated with the "coal-to-NG" program include: (1) no catalyst available for properly operating for vehicles powered by gasoline and NG; and (2) lack of denitrification facilitates after the replacement. Moreover, the net change in pollution levels after excluding perturbations caused by varying weather conditions is still unknown.
To explore the impact of increased NG consumption on ambient mass concentrations of criteria air pollutants, the influence of varying weather conditions on their interannual variations must first be minimized [17][18][19][20][21][22][23]. In the literature, several approaches have been documented to decouple the influence of varying weather conditions from emission-driven changes on the interannual variations in concentrations of air pollutants [19,[24][25][26][27]. In these approaches, either meteorological and other related parameters or random noises were introduced to estimate the deweathered concentrations of criteria pollutants. However, the uncertainties of the deweathered concentrations were poorly understood. Yao and Zhang proposed an approach to convert the time series of data to concentration series of data so that the perturbation from varying weather conditions can be minimized in assessing the effects of emission reduction on the wet deposition trends of several inorganic ions [21]. Building on this study, we further developed an approach to determine the closed interval of the deweathered percentage change (DPC) in ambient concentrations of criteria pollutants on an annual scale, with the aim of quantifying mitigation effects after minimizing weather impacts on air pollution.
Qingdao is a coastal megacity in northern China that pioneered in the implementation of the "coal-to-NG" program starting in 2018. The annual statistical reports of Qingdao showed that NG consumption almost doubled from 2017 to 2018-2019, and its contribution to the total energy consumption significantly increased from 5% in 2017 to 8-9% in 2018-2019, even though the total energy consumption showed little change in 2017-2019 [28]. Thus, it would be an ideal location for studying the impact of the "coal-to-NG" program on air quality by comparing it with the impact of implemented policies in 2013-2017.
In this study, we focused on interannual changes in PM 2.5 , PM 2.5-10 (calculated as PM 10 minus PM 2.5 ), CO, NO 2 , and SO 2 ; however, the use of NO x instead of NO 2 can better reflect NO x emissions. We applied the closed intervals of the DPC method to observations in Qingdao during the heating season (from November to the following March) as moderate and heavy air pollution can sometimes occur [17,18,22]. We then rationalized the closed intervals of DPC in concentrations of criteria pollutants in terms of the mitigation measures implemented in different years, with particular attention on the significant increase in NG consumption since 2018. Overall, this study aims at evaluating the effects of the "coalto-NG" program on air quality by minimizing the perturbation from varying weather conditions. In addition, PM 2.5-10 was considered as a coarse particle in this study. Dust storms, dust from unsealed roads, dust from construction and demolition, road traffic, and marine aerosols were usually identified as the major emission sources of coarse particles. Except for dust storms, PM 2.5-10 was more related to local emissions rather than long-range transport because of its relatively short residence time in the atmosphere.

Study Location
Qingdao is located approximately 500 km from Beijing and Tianjin in northern China ( Figure 1). Fourteen air quality monitoring sites, which have been operating since 2013 or earlier and are distributed mainly in urban and semi-urban areas in Qingdao (Figure 1), were chosen to analyze the annual means of air pollutants. All the hourly mean air pollutant concentrations are reported online, such as on IQAir [29], and have been updated every hour in Qingdao since 2013. Hourly mean air pollutant concentrations at the city level were calculated as the mean values of those simultaneously measured at the individual stations. For any hour with missing data in five or more stations, the city-level hourly concentration was also treated as missing data. A total of 3624 (or 3648 in a leap year) hourly data points for each pollutant would be available for one heating season if no data were missing. Approximately 93-94% of the hourly data of the criteria air pollutants for each heating season were available; for example, the number of hourly PM 2.5 data ranged from 3372 to 3407 during the heating season from 2013 to 2019.
NOx emissions. We applied the closed intervals of the DPC method to observat Qingdao during the heating season (from November to the following March) as m and heavy air pollution can sometimes occur [17,18,22]. We then rationalized the intervals of DPC in concentrations of criteria pollutants in terms of the mit measures implemented in different years, with particular attention on the signifi crease in NG consumption since 2018. Overall, this study aims at evaluating the ef the "coal-to-NG" program on air quality by minimizing the perturbation from v weather conditions. In addition, PM2.5-10 was considered as a coarse particle in this Dust storms, dust from unsealed roads, dust from construction and demolition, ro fic, and marine aerosols were usually identified as the major emission sources of particles. Except for dust storms, PM2.5-10 was more related to local emissions rath long-range transport because of its relatively short residence time in the atmosphe

Study Location
Qingdao is located approximately 500 km from Beijing and Tianjin in northern ( Figure 1). Fourteen air quality monitoring sites, which have been operating since earlier and are distributed mainly in urban and semi-urban areas in Qingdao (Fig were chosen to analyze the annual means of air pollutants. All the hourly mean ai tant concentrations are reported online, such as on IQAir [29], and have been u every hour in Qingdao since 2013. Hourly mean air pollutant concentrations at level were calculated as the mean values of those simultaneously measured at th vidual stations. For any hour with missing data in five or more stations, the ci hourly concentration was also treated as missing data. A total of 3624 (or 3648 in year) hourly data points for each pollutant would be available for one heating se no data were missing. Approximately 93-94% of the hourly data of the criteria ai tants for each heating season were available; for example, the number of hourly PM ranged from 3372 to 3407 during the heating season from 2013 to 2019.  In addition, air quality in Qingdao has notably improved each year in response to the implementation of various mitigation measures since 2013, which occurred prior to the large increase in NG consumption in 2018. For example, the annual mean mass concentrations of PM 2.5 , PM 10 , and SO 2 were 44%, 29%, and 74% lower, respectively, in 2017 than in 2013 [30]. However, the annual mean mass concentrations of PM 2.5 , PM 10 , and NO 2 were 9%, 3%, and 3% higher in 2019, respectively, than in 2018 [30].

Study Period
This study focused on the heating season, which includes November to the following

Energy Consumption during Study Period
The consumption of coal, oil, and NG in Qingdao from 2013 to 2020 is shown in Figure S1 in Supplementary Materials, and the data were downloaded from the Qingdao Bureau of Statistics [28]. NG consumption sharply increased from 3.7-4.5 × 10 12 kcal in 2013-2017 to 8.5-9.3 × 10 12 kcal in 2018-2019 ( Figure S1 in Supplementary Materials). Relative to the levels in 2017, coal consumption decreased from 6.9 × 10 13 kcal to 6.6-6.8 × 10 13 kcal and oil consumption decreased from 3.1 × 10 13 kcal to 2.4-2.7 × 10 13 kcal in 2018-2019. The total energy consumption narrowly fluctuated around 1 × 10 14 kcal between 2014 and 2019, but the contribution of NG consumption to the total energy consumption increased from 4-5% in 2014-2017 to 8-9% in 2018-2019. The total energy consumption in Qingdao decreased to 9.6 × 10 13 kcal in 2020 because of the COVID-19 pandemic. During the COVID-19 lockdown period in Qingdao, the total energy consumption and power plant energy consumption decreased by 8-9% and 19-23%, respectively, compared to those of 2019 [28]. The total global energy consumption and air pollutant emissions were reportedly led to a decrease during the COVID-19 lockdown period, based on daily global CO 2 emissions estimates, satellite data, and thousands of ground-level air pollutant observations [31][32][33]. Thus, we assume that energy consumption in areas neighboring and upwind of Qingdao also decreased in early 2020.
Based on our incomplete survey (no official data in the public domain, personal communication with people managing facilities that adopted "coal-to-NG"), the large increase in NG consumption was mainly attributed to two large chemical facilities (over 80%) to support new production capacity. Less than 20% of the NG consumption increase was attributed to small facilities that were serviced for old production capacities. The national NG consumption in China only increased by approximately 1% per year in 2017, 2018, and 2019 [34].

Climate Anomalies during Study Period
Meteorological data in Qingdao were downloaded from http://data.cma.cn/ (assessed on 30 November 2021) with permission after registration. A total of 1200 (or 1208 in a leap year) hourly data for precipitation, wind speed (WS), air temperature (T), and relative humidity (RH) were for one heating season if no data were missing. More than 95% of the data were available for each heating season. The anomalies of climate factors during the  [18,20,22], disfavoring the dispersion of air pollutants. However, a large increase in precipitation during 2019/2020 may have overridden the impact of warm winters. The ambient RH anomaly from the decade seasonal mean of 62% varied from −3.1% to 4.4% ( Figure S2d in Supplementary Materials). The anomalies were −1.1% and +4.4% in 2018/2019 and 2019/2020, respectively. Increased ambient RH may accelerate the heterogeneous formation of secondary PM 2.5 [3,6], while increased precipitation can minimize this impact to some extent.

Method of Closed Interval of DPCs in Concentrations of Air Pollutants
Unlike the estimation of the deweathered concentrations of criteria air pollutants adopted in the literature [19,[24][25][26][27], this study alternatively determined the closed interval of the DPC of criteria air pollutant mass concentrations between any two years. To achieve the target, we developed a five-step approach for data processing and analysis because of criteria air pollutant mass concentrations in different years with varying data sizes.
Step 1 reconstructed each year's dataset while retaining the original statistical metrics so that each year's dataset was the same size.
Step 2 conducted a correlation analysis of the observational data between two consecutive years using the newly constructed datasets from Step 1. The infection point is visibly identified from the regression curve as the first guess and all data points having a value larger than the inflection point are considered as outliers caused by varying weather conditions or extreme events.
Step 3 statistically screened out more outliers by repeating the correlation analysis using the datasets that excluded the outliers identified in Step 2.
Step 4 calculated the range of the true DPC between two years, and Step 5 evaluated the residual perturbation, assuming that the impact of varying weather conditions may not be fully excluded from the first four steps. Detailed descriptions of the five steps are provided in the Supplementary Materials. Here, we summarize the key points of each step.
Throughout Steps 1-3, 91-98% of the observational data were left to reconstruct the final arrays for a given pollutant between any two consecutive years of the same size. In Step 4, the linear regression (LR) analysis was repeated for each pair of the final arrays with zero intercept. The LR slope was defined as the primary DPC (DPC primary ). Moreover, the LR with a non-zero intercept was conducted using the data in each pair of the final arrays. The LR slope was defined as the secondary DPC (DPC secondary ). The R 2 values of the LR analysis for calculating both DPC secondary and DPC primary or calculating DPC secondary were mostly larger than 0.99, except for a few values larger than 0.96. The true DPC should be theoretically between DPC primary and DPC secondary .
Step 5 evaluated the residual perturbation by varying weather conditions. Theoretically, the intercept in the LR analysis of Step 4 should infinitely approach zero when the Atmosphere 2022, 13, 945 6 of 14 perturbation by varying weather conditions is negligible. This can be clearly identified in approximately 1/3 of the calculated pairs of DPC primary and DPC secondary . Moreover, higher consistency in each pair of DPC primary and DPC secondary was found for PM 2.5 , the concentrations of which were affected mainly by regional sources [35]. The reverse was generally true for PM 2.5-10 , SO 2 , CO, and NO 2 , whose concentrations were determined largely by local emissions. If DPC secondary is smaller than DPC primary , the varying weather conditions incompletely removed by Step 3 likely disfavor the accumulation and/or formation of the pollutant. The reverse would be true if DPC secondary was larger than DPC primary . Thus, DPC primary and DPC secondary were used to construct the closed interval of DPC in concentrations of criteria pollutants, i.e., (DPC primary , DPC secondary ).
The  Table 1. The DPC can be considered as the cumulative effect of mitigation measures in multiple years and the corresponding interval was hereafter referred to as the closed interval of cumulative DPC in this study. There are four advantages of the above-described approach: (1) allowing statistical identification and exclusion of outliers in estimating two endpoints of the closed interval of DPC; (2) confirming the accuracy of DPC when the interval length of DPC is sufficiently small; (3) identifying the large perturbation from varying weather conditions on DPC when the interval length is large; and (4) avoiding direct estimation of deweathered concentrations of air pollutants in which the uncertainties are difficult to accurately evaluate. When the third is the case, the closed intervals can be narrowed simply by combining results extracted from the approaches reported in the literature. Moreover, the links between emissions and ambient concentrations of air pollutants are generally nonlinear. The perturbation from varying weather conditions on ambient concentrations cannot be fully excluded. Thus, the true DPC should always be a closed interval, but not a single value. The closed interval of DPC was thereby obtained, i.e., (−19%, −20%), with the interval length as low as 1%. Thus, the impact of perturbations caused by varying weather conditions on the estimated DPC were negligible. The differences between DPC primary and DPC secondary were within 2% for 2014/2015, 2016/2017, 2017/2018, and 2019/2020 compared with those in each previous year, and the values were negative (Figure 3a); however, the DPC primary and DPC secondary were positive and their differences were slightly larger for 2015/2016 (3% and 8%) and 2018/2019 (30% and 33%), respectively. Note that the difference in PM 2.5 concentrations between two consecutive years was generally significant (p < 0.01), except for 2015/2016 vs. 2014/2015 (Figure 3a). of PM2.5 were determined and further discussed.
The hourly mean data of PM2.5 concentration reconstructed by Step 1 during the ing season between two consecutive years were plotted in Figure 2a  Step 1 (a-f) (blue markers and regression curves use all data points, and yellow ones use the se data by excluding data points that suffered from severe perturbations from the anomalies). Step 1 (a-f) (blue markers and regression curves use all data points, and yellow ones use the selected data by excluding data points that suffered from severe perturbations from the anomalies).   [31,33] Clearly, the interval lengths of the cumulative DPCs were larger than any of the aforementioned values for two consecutive years. Investigating the closed intervals of both cumulative DPCs and DPCs in any two consecutive years can solidify the emission-driven increase or decrease in any particular year.

Closed Intervals of DPCs in PM2.5-10
The invasion of dust storms, which can be easily identified in satellite images and supported by PM2.5-10 concentrations of >150 µg m −3 combined with a mass PM2.5-10 to PM2.5 ratio of >2, usually occurs during the heating season and early spring in northern China [3,17]. Data points with PM2.5-10 >150 µg m −3 generally deviate from the regression curve to a large extent (Figure 4a-f); these data points were therefore removed from the calculation of the closed interval of DPC because they are unrelated to the mitigation of air pollutants.  [31,33]. Clearly, the interval lengths of the cumulative DPCs were larger than any of the aforementioned values for two consecutive years. Investigating the closed intervals of both cumulative DPCs and DPCs in any two consecutive years can solidify the emission-driven increase or decrease in any particular year.

Closed Intervals of DPCs in PM 2.5-10
The invasion of dust storms, which can be easily identified in satellite images and supported by PM 2.5-10 concentrations of >150 µg m −3 combined with a mass PM 2.5-10 to PM 2.5 ratio of >2, usually occurs during the heating season and early spring in northern China [3,17]. Data points with PM 2.5-10 >150 µg m −3 generally deviate from the regression curve to a large extent (Figure 4a-f); these data points were therefore removed from the calculation of the closed interval of DPC because they are unrelated to the mitigation of air pollutants.
Following the same approach used for PM 2.5 , the estimated endpoint values of the PM 2.5-10 closed interval of DPC for each year before 2016/2017 were mostly negative compared with those of each previous year (Figure 3b). The closed interval of cumulative DPC from 2013/2014 to 2016/2017 were obtained at (−29%, −31%), respectively (Table 1). There were somewhat large interval lengths of DPC between two consecutive years (e.g.,  Step 1 (a-f) (blue markers and regression curves use all data points, and yellow ones use the select data by excluding data points that suffered from severe perturbations from the anomalies).
Following the same approach used for PM2.5, the estimated endpoint values of t PM2.5-10 closed interval of DPC for each year before 2016/2017 were mostly negative com pared with those of each previous year (Figure 3b). The closed interval of cumulative DP from 2013/2014 to 2016/2017 were obtained at (−29%, −31%), respectively ( Table 1) Step 1 (a-f) (blue markers and regression curves use all data points, and yellow ones use the selected data by excluding data points that suffered from severe perturbations from the anomalies).  (Table 1).

Closed Intervals of DPCs in SO2, CO, and NO2
The estimated DPC in concentrations of SO2 were always negative compared with those in each previous year, with the largest two endpoint values of −16% and −17% from  (Table 1).

Discussion
The observed SO 2 in Qingdao was generally derived from local sources associated with industrial sectors, while the SO 2 from long-range transport is typically diluted or oxidized to a large extent [17,18,33]. The significant negative DPCs in SO 2 throughout the study period are likely due to the local adoption of continuous mitigation measures. NG consumption is commonly regarded as a negligible source of SO 2 , and NG replacing coal should therefore decrease SO 2 emissions to some extent. During the COVID-19 lockdown, reduced industrial activities also lowered SO 2 emissions.
The concentrations of CO and NO 2 usually exhibit large spatial variations in Qingdao [29], indicating overwhelming contributions from local sources. The closed interval of cumulative DPC in concentrations of CO from 2013/2014 to 2017/2018 was obtained at (−25%, −18%) (Table 1) . The rebound in 2018/2019 was consistent with a recent report that found that the increased NG consumption in Beijing and Tianjin in northern China may have caused increased concentrations of NO x and particulate nitrate [13].
Although the DPCs in concentrations of SO 2 and NO 2 were both highly negative (<−16%) from 2014/2015 to 2015/2016, the corresponding DPC in concentration of PM 2.5 was positive at (3%, 8%). This suggests that, in addition to local sources, long-range transport also contributed significantly to the PM 2.5 concentration in this city [18,38], particularly because of the longer residence time of PM 2.5 relative to its gaseous precursors. The exact causes of the unexpected positive DPC of PM 2.5 in 2015/2016 require further investigation in terms of its chemical composition (e.g., inorganic/organic fractions) and dominant air mass origins, etc. The highly positive DPC of PM 2.5 from 2017/2018 to 2018/2019 at (30%, 33%) indicates a largely emission-driven rebound in PM 2.5 , similar to the cases of CO and NO 2 , which raises the possibility of increased secondary aerosol formation from increased gaseous precursors at the local scale. Although direct PM 2.5 emissions from NG combustion should be negligible, incomplete NG combustion may emit organic species with various volatilities [39], which can act as important gaseous precursors of secondary organic aerosols. Further investigation is therefore needed to understand NG emission-related secondary aerosol formation [40]. As newly formed secondary aerosols are mostly fine particles [3,41,42], the above discussion can also partially explain the positive DPC from 2017/2018 to 2018/2019. The DPCs of PM 2.5-10 were generally negative during 2014/2015-2016/2017, although the interval lengths were relatively large (Figure 3b). The negative DPCs in the first several years were likely due to the overall effects associated with the implemented measures (Tables S1-S3 in Supplementary Materials). In contrast, the positive DPC values in the subsequent two years were likely due to accelerated construction activities (Table S3 in Supplementary Materials), as well as increased local dust emissions. However, the reduced construction activities, on-road vehicle use, and industrial activities because of the COVID-19 lockdown in 2019/2020 changed the DPC to be negative.

Conclusions
In this study, we developed a new approach to calculate the closed intervals of DPC in air pollutant concentrations between two consecutive years, as well as the closed intervals of cumulative DPC in a year compared with that in the base year. Thus, the changes in concentrations of air pollutants caused by varying weather conditions were excluded and then the emission-driven changes in concentrations of air pollutants could be studied. The obtained closed intervals indicated that the DPCs extracted before 2018/2019 were generally negative in Qingdao, mainly due to the implementation of mitigation measures. However, the positive DPCs, to be distributed at (27%, 30%) for CO, (15%, 18%) for NO 2 , and (30%, 33%) for PM 2.5 with small interval lengths, from 2017/2018 to 2018/2019 implied emissiondriven rebounds of these pollutants rather than changes in meteorological conditions. The comparison between the closed intervals of cumulative DPCs from 2013/2014 to 2018/2019 and those from 2013/2014 to 2017/2018 also supported the analyses. The rebound was likely associated with a doubled increase in NG consumption in 2018, as the total energy consumption showed little difference. By combining our results for Qingdao with previously reported results for the Beijing-Tianjin-Hebei region, we conclude that the "coal-to-NG" program likely caused the increase in NO 2 concentrations on a large scale across northern China. With more official statistics going public, integrated assessments of "coal-to-NG" could be provided.
The obtained closed intervals of DPCs suffered from a broad range in many cases. The interval length was generally even broader for the calculated cumulative DPC. The broad interval length implied that the perturbations caused by varying weather conditions cannot be eliminated. In the cases, a combination of multiple meteorological normalization techniques reported in the literature may be needed to narrow the interval.
Supplementary Materials: The following supporting information can be downloaded at https://www. mdpi.com/article/10.3390/atmos13060945/s1: Table S1: Mitigation measures implemented to reduce emissions of air pollutants from coal combustion in Qingdao from 2014 to 2019; Table S2: Mitigation measures implemented to reduce emissions of air pollutants from road dust in Qingdao from 2014 to 2019; Table S3: Mitigation measures implemented to reduce emissions of VOC in Qingdao from 2014 to 2019; Table S4: Statistical comparison between raw array (A 2 ) and reconstructed array (B 2 ) for hourly average concentrations of PM 2.5 in 2014 heating season; Figure S1: Annual coal, oil and natural gas consumption (a) and annual total energy consumption and the percentage of NG in total energy consumption (b) in Qingdao from 2013 to 2020; and Figure S2