The Effects of Fireworks Discharge on Atmospheric PM2.5 Concentration in the Chinese Lunar New Year

Discharging fireworks during the Chinese Lunar New Year celebrations is a deep-rooted custom in China. In this paper, we analyze the effect of this cultural activity on PM2.5 concentration using both ground observations and satellite data. By combining remote sensing data, the problem of uneven spatial distribution of ground monitoring has been compensated, and the research time span has been expanded. The results show that the extensive firework displays on New Year’s Eve lead to a remarkable increase in nationwide PM2.5 concentration, which were 159~223% of the average level, indicating the instantaneous effect far exceeds that of any other factor over the whole year. However, the averaged PM2.5 concentrations of the celebration period were 0.99~16.32 μg/m3 lower compared to the average values of the corresponding pre-celebration period and post-celebration period, indicating the sustained effect is not very significant. The implementation of firework prohibition policies can greatly reduce the instantaneous PM2.5 increase, but no obvious air quality improvement is observed over the entire celebration period. Combining these findings and the cultural significance of this activity, we recommend that this custom is actively maintained, using new technologies and scientific governance programs to minimize the negative effects.


Introduction
Discharging fireworks during the Chinese Lunar New Year is a custom that has continued for thousands of years. As early as the Northern and Southern Dynasties (420-589 AD), a book called "The Chronicle of Jingchu" had already recorded this custom [1]. According to the legend that a cannibal named Nian always attacked villages in winter. In the fight against Nian, people unexpectedly discovered that Nian did not like the color red and was afraid of flames and explosions, so they threw bamboo into the fire, and the fire plus the crackling sound successfully scared the monster away. Subsequently, this act of repelling the monster became a symbolic part of the New Year celebration. The word "Nian", which has the same pronunciation as the word "year" in Chinese, has also become a synonym for "Chinese New Year". Over time, with the development of gunpowder, fireworks and firecrackers gradually replaced bamboo and became the core element to the Chinese New Year celebrations.
However, despite the profound cultural significance of this traditional custom, the atmospheric pollution it causes has become a topic of national interest [2][3][4]. In recent years, China has faced

Study Region and Data
The study region is China ( Figure S2), and the research period is mainly from 18 January 2013, to 6 November 2016. Three kinds of data were used in this study: (1) Ground PM 2.5 observations. In recent years, the Chinese government has accelerated the construction of the PM 2.5 ground monitoring network. By the end of 2016, the number of online sites had reached 1600. The collection and processing of pollutant data is undertaken in accordance with the Technical Regulation on Ambient Air Quality Index (AQI) and national quality control guidelines [50,51]. (2) Moderate Resolution Imaging Spectroradiometer (MODIS) aerosol optical depth (AOD) products.
Compared with station observations, satellite data enjoy a wider coverage and longer observation times, which means that satellite data show significant advantages in air quality related research [52][53][54][55][56][57][58][59]. In this study, MODIS Level 2 Collection 6 10-km AOD products were adopted for the national PM 2.5 retrieval. The quality of MODIS AOD products has been proven to meet the requirements of atmospheric-related research [60][61][62]. The newly released 10-km AOD products in collection 6 have further enhanced the retrieval capacity in the highlighted regions by combining the dark target algorithm and deep blue algorithm in AOD retrieval. The overall correlation between Collection 6 AOD and AERONET-observed AOD over land is R = 0.86, and 69.4% of Collection 6 AOD fall within expected uncertainty of ±(0.05 + 12%) [63]. The Terra and Aqua are polar-orbiting satellites and cross the equator at around 10:30 a.m. (descending orbit) and 1:30 p.m. (ascending orbit) local sun times, respectively. Since the crossing time of these two sensors is different, the data of Terra and Aqua should be combined to better represent daily AOD. For each pixel, if there is no data for both Terra and Aqua products, it was identified as "null"; if both values are identified, the average value was used; if there is only one valid value, this value was taken. (3) Meteorological data. Meteorological conditions can greatly influence the distribution of atmospheric pollutants. The NCEP GDAS/FNL 0.25 Degree Global Tropospheric Analyses and Forecast Grids product is used in this research. These data have a 6-hr temporal resolution and fully cover the research area since 1991. We extracted relative humidity (RH, %), temperature (TEMP, K), wind speed, (WS, m/s), surface pressure (SP, Pa), and height of planetary boundary layer (HPBL, m) (Table S1) and resampled them into 0.1 degree (10 km) for historical PM 2.5 retrieval. Because the spatial resolution of the reanalysis data was relatively coarse, the ground monitoring data was a better choice for the detailed analysis. The daily meteorological ground observations were also used to supplement the analysis of the firework policy effects.

Historical PM 2.5 Retrieval
Ground observations are extensively used in surface air pollution. However, long-term research of China's PM 2.5 pollution has been severely restricted because of the late construction of the ground PM monitoring network. Besides, the distribution of ground stations is uneven. Thus, the remote sensing products have become popular in ground PM 2.5 monitoring, because they can not only extend the research time, but also obtain the spatially continuous PM 2.5 data [56,59,64].
Several attempts have been made to describe the AOD-PM relationship, and most of them can be classified into simulation-based models or observation-based models [65]. Because the requirement of datasets is more comprehensive for simulation-based models (especially emission inventories can be hard to acquire sometimes) [64,66], the observation-based models can be a good compromise [67].
According to previous researches [68][69][70][71], the AOD-PM relationship has gradually been regarded as a nonlinear problem of multiple variables. Therefore, artificial neural networks, which can better present complex nonlinear relationships, have been used to estimate PM 2.5 concentrations [67,72,73].
A three-layer back-propagation neural network (BPNN) model was constructed to estimate historical daily PM 2.5 [67].
It takes a total of nine factors as input, including spatiotemporal information, AOD, and meteorological elements, and gives PM 2.5 as the only output, the model structure is shown in Figure 1. The optimal performance node number of the hidden layer is usually between 2 + µ and 2n + 1 (where n and µ represent the number of nodes in the input layer and the output layer, respectively) [74]. In previous work, 18 nodes were proved to have the best performance [75]. The model performance was evaluated using the correlation coefficient (R) and root-meansquare error (RMSE) between the measured PM2.5 and the retrieved PM2.5, and 10-fold crossvalidation was introduced to assess the possibility of model overfitting [76]. The formulas of R and RMSE are defined as follows: where n describes the sample size; and are the observation and prediction values of the PM2.5 concentrations, respectively; * is the average value of PM2.5 observations.
The R describes the fitting degree of the constructed model, and usually the higher value indicates the better model reliability. The RMSE describes the prediction error of the model, and usually the lower value indicated the better model accuracy. As for the 10-fold cross-validation method, the sample-based cross-validation is adopted to evaluate the overall predictive performance of satellite-based PM2.5 retrieval models [77]. In sample-based cross-validation, all the samples are almost equally divided into 10 folds randomly, and the complete validation needs 10 rounds. Each round, nine folds will be used for model fitting and the remaining one will be used for model validation.
In this work, the main purpose of remote sensing data retrieval is to obtain historical ground PM2.5 data, so as to explore the influence of the Chinese New Year fireworks on PM2.5 concentrations in a longer time span. Thus, we used the data from 1 March 2014 to 29 February 2016 to construct the model, to estimate the data from 1 March 2013 to 28 February 2014. The model fitting R is 0.7581 and RMSE is 27.1316 μg/m³, the model cross-validation R is 0.7523 and RMSE is 27.4554 μg/m³ (see Figure  2a). These values indicated that the model is not seriously overfitting. The model prediction R is 0.7243 and RMSE is 39.9400 μg/m³ (see Figure 2b), and the model 17-day prediction R is 0.7930 and RMSE is 26.9968 μg/m³ (see Figure 2c). These values indicated that the model can be considered robust for historical PM2.5 retrieval, especially for the latter historical 17-day fusion research. Furthermore, we retrieved PM2.5 data from 2002 to 2012 for the evaluation of fireworks' long-term The model performance was evaluated using the correlation coefficient (R) and root-mean-square error (RMSE) between the measured PM 2.5 and the retrieved PM 2.5 , and 10-fold cross-validation was introduced to assess the possibility of model overfitting [76]. The formulas of R and RMSE are defined as follows: where n describes the sample size; y i and y i are the observation and prediction values of the PM 2.5 concentrations, respectively; y * is the average value of PM 2.5 observations. The R describes the fitting degree of the constructed model, and usually the higher value indicates the better model reliability. The RMSE describes the prediction error of the model, and usually the lower value indicated the better model accuracy. As for the 10-fold cross-validation method, the sample-based cross-validation is adopted to evaluate the overall predictive performance of satellite-based PM 2.5 retrieval models [77]. In sample-based cross-validation, all the samples are almost equally divided into 10 folds randomly, and the complete validation needs 10 rounds. Each round, nine folds will be used for model fitting and the remaining one will be used for model validation.
In this work, the main purpose of remote sensing data retrieval is to obtain historical ground PM 2.5 data, so as to explore the influence of the Chinese New Year fireworks on PM 2.5 concentrations in a longer time span. Thus, we used the data from 1 March 2014 to 29 February 2016 to construct the model, to estimate the data from 1 March 2013 to 28 February 2014. The model fitting R is 0.7581 and RMSE is 27.1316 µg/m 3 , the model cross-validation R is 0.7523 and RMSE is 27.4554 µg/m 3 (see Figure 2a). These values indicated that the model is not seriously overfitting. The model prediction R is 0.7243 and RMSE is 39.9400 µg/m 3 (see Figure 2b), and the model 17-day prediction R is 0.7930 and RMSE is 26.9968 µg/m 3 (see Figure 2c). These values indicated that the model can be considered robust for historical PM 2.5 retrieval, especially for the latter historical 17-day fusion research. Furthermore, we retrieved PM 2.5 data from 2002 to 2012 for the evaluation of fireworks' long-term effect on PM 2.5 concentrations.

The 17-Day Time Division Schema
In this research, we consider the entire Chinese New Year celebration period as a whole and denote its time span (17 days) as Δt. Using Δt as division unit, and taking each celebration period as an anchor point, we were able to divide the time before and after each celebration period. As the result, the time from 2002 to 2016 can be circularly divided into twenty-one periods in sequence (the schematic diagram is shown in Figure 3a). The order of periods is: the 10th pre-celebration period~ the 1st pre-celebration period, the celebration period, and the 1st post-celebration period ~ the 10th post-celebration period.
Since the date of the Chinese New Year was different every year (see Table S2), after divided the 9th periods, the number of days left in each year would be slightly different. These days were mainly in the summer and had little effect on the study of the Chinese New Year celebration periods in winter. Therefore, we still divided the 10th post-celebration period by Δt and put the rest of the days into the 10th pre-celebration period. Thus, each twenty-one period is composed of twenty 17-day periods and one residual period (the 10th pre-celebration period).

The 17-Day Time Division Schema
In this research, we consider the entire Chinese New Year celebration period as a whole and denote its time span (17 days) as ∆t. Using ∆t as division unit, and taking each celebration period as an anchor point, we were able to divide the time before and after each celebration period. As the result, the time from 2002 to 2016 can be circularly divided into twenty-one periods in sequence (the schematic diagram is shown in Figure 3a). The order of periods is: the 10th pre-celebration periodt he 1st pre-celebration period, the celebration period, and the 1st post-celebration period~the 10th post-celebration period.
Since the date of the Chinese New Year was different every year (see Table S2), after divided the 9th periods, the number of days left in each year would be slightly different. These days were mainly in the summer and had little effect on the study of the Chinese New Year celebration periods in winter. Therefore, we still divided the 10th post-celebration period by ∆t and put the rest of the days into the 10th pre-celebration period. Thus, each twenty-one period is composed of twenty 17-day periods and one residual period (the 10th pre-celebration period).

The Historical PM2.5 Data Fusion Schema
As mentioned above, ground observations may have some deficiencies in related researches due to the construction time and location constraints of the stations. At present, the satellite retrieved PM2.5 data are generally considered to be good supplementary data. However, constrained by imaging conditions, such as bright ground, clouds, and fog, most satellite products are troubled by data missing. Thus, the retrieved PM2.5 data usually cannot achieve full daily coverage. Therefore, it is necessary to adopt appropriate fusion strategies to improve data coverage before further analysis.
Through the 17-day time division schema, the division results can be illustrated as Figure 3b. To compare the PM2.5 levels between different time periods across these years, the retrieved data of the same period from fifteen years were fused. For example, the data of every Chinese New Year celebration period in 2002-2016 were fused into one distribution map of PM2.5 concentrations (illustrated in Figure 3b with red letters). The data of other periods were fused with the same strategy.

Instantaneous Effect of Fireworks Discharge on PM2.5
Based on all the available ground monitoring data from January 18, 2013, to November 6, 2016 (including the four celebration periods, see Table S3), we obtained continuous hourly average PM2.5 concentrations for mainland China.
Firstly, we recorded the daily maximum PM2.5 ( Figure 4a). It can be seen that the daily maximum PM2.5 reaches its lowest point in summer and its highest point in winter, with a clear pattern of seasonal variation [55,78]. The severe PM2.5 pollution in winter is likely the result of winter heating [27,28] and stable meteorological conditions [79,80].

The Historical PM 2.5 Data Fusion Schema
As mentioned above, ground observations may have some deficiencies in related researches due to the construction time and location constraints of the stations. At present, the satellite retrieved PM 2.5 data are generally considered to be good supplementary data. However, constrained by imaging conditions, such as bright ground, clouds, and fog, most satellite products are troubled by data missing. Thus, the retrieved PM 2.5 data usually cannot achieve full daily coverage. Therefore, it is necessary to adopt appropriate fusion strategies to improve data coverage before further analysis.
Through the 17-day time division schema, the division results can be illustrated as Figure 3b. To compare the PM 2.5 levels between different time periods across these years, the retrieved data of the same period from fifteen years were fused. For example, the data of every Chinese New Year celebration period in 2002-2016 were fused into one distribution map of PM 2.5 concentrations (illustrated in Figure 3b with red letters). The data of other periods were fused with the same strategy.

Instantaneous Effect of Fireworks Discharge on PM 2.5
Based on all the available ground monitoring data from January 18, 2013, to November 6, 2016 (including the four celebration periods, see Table S3), we obtained continuous hourly average PM 2.5 concentrations for mainland China.
Firstly, we recorded the daily maximum PM 2.5 (Figure 4a). It can be seen that the daily maximum PM 2.5 reaches its lowest point in summer and its highest point in winter, with a clear pattern of seasonal variation [55,78]. The severe PM 2.5 pollution in winter is likely the result of winter heating [27,28] and stable meteorological conditions [79,80]. It should be noted that these are the highest values of the corresponding years, and these peaks all appeared at 2 a.m. on Lunar New Year's Day. So, we recorded the time when the national average PM2.5 reached the maximum each day to observe its distribution. According to Figure 4b, daily PM2.5 usually reached a maximum around 7:00-11:00 or 21:00-24:00 (00:00). Of the 1389 days considered in this study, only 12 maximum daily PM2.5 concentrations appeared at 2 a.m., and the PM2.5 concentrations of the other eight days were far below the PM2.5 level of each Lunar New Year (i.e., 94 μg/m³, 61 μg/m³, 62 μg/m³, 40 μg/m³, 99 μg/m³, 72 μg/m³, 69 μg/m³, and 38 μg/m³, respectively).
Lunar New Year's Eve and Lunar New Year's Day are the time when most intensive fireworks discharge takes place. In terms of time, these New Year peaks are most likely caused by fireworks. Meanwhile, the closure of many factories and enterprises has greatly suppressed the PM2.5 pollutants emitted by economic activities [45,81]. Under such favorable conditions, the daily maximum PM2.5 still reaches extremely high levels on the Lunar New Year, indicating that the Chinese Lunar New Year fireworks discharge has the strongest transient influence in the whole year. These conclusions confirmed some city-level conclusions on a national scale.
Although we had confirmed that the effect of fireworks discharge is serious nationwide, we still wanted to check if this effect was equally significant in individual cities. We therefore decided to conduct an analysis in 31 provincial capital cities, including the capital Beijing, because ground data are usually more widely available in these cities. Because fireworks seem to have the greatest impact on PM2.5 in the morning of New Year's Day ( Figure S3), and PM2.5 tends to reach its highest level at 2 a.m., we decided to explore the response of PM2.5 to fireworks based on the time window of 00:00-04:00. During each celebration period, each city was categorized into one category (significant response, obvious response, or weak response), according to the response of PM2.5 to fireworks discharge.
The classification was conducted as follows: (1) If the maximum value of PM2.5 in a city appeared during 0:00-04:00 on New Year's Day, it was considered that the city's PM2.5 responded significantly to the New Year fireworks. (2) If the maximum value of PM2.5 in a city did not appear during 0:00-04:00 on New Year's Day, but both the increase and the increase rate were the highest compared with those during the Meanwhile, the closure of many factories and enterprises has greatly suppressed the PM 2.5 pollutants emitted by economic activities [45,81]. Under such favorable conditions, the daily maximum PM 2.5 still reaches extremely high levels on the Lunar New Year, indicating that the Chinese Lunar New Year fireworks discharge has the strongest transient influence in the whole year. These conclusions confirmed some city-level conclusions on a national scale.
Although we had confirmed that the effect of fireworks discharge is serious nationwide, we still wanted to check if this effect was equally significant in individual cities. We therefore decided to conduct an analysis in 31 provincial capital cities, including the capital Beijing, because ground data are usually more widely available in these cities. Because fireworks seem to have the greatest impact on PM 2.5 in the morning of New Year's Day ( Figure S3), and PM 2.5 tends to reach its highest level at 2 a.m., we decided to explore the response of PM 2.5 to fireworks based on the time window of 00:00-04:00.
During each celebration period, each city was categorized into one category (significant response, obvious response, or weak response), according to the response of PM 2.5 to fireworks discharge.
The classification was conducted as follows: (1) If the maximum value of PM 2.5 in a city appeared during 0:00-04:00 on New Year's Day, it was considered that the city's PM 2.5 responded significantly to the New Year fireworks. (2) If the maximum value of PM 2.5 in a city did not appear during 0:00-04:00 on New Year's Day, but both the increase and the increase rate were the highest compared with those during the 0:00-04:00 period of the remaining days, it was considered that the city's PM 2.5 responded obviously to the New Year fireworks. (3) If the variation of PM 2.5 in a city did not match the above two conditions, it was considered that the city's PM 2.5 responded weakly to the New Year fireworks.
To avoid any possible interference, the data from the Lantern Festival period were excluded from this classification. In addition, according to the Tibetan calendar (Table S4) and related customs, the data from the Tibetan New Year firework period would not affect the classification results.
According to the classification map shown in Figure 5 and the corresponding statistical results (Table S5), it is unsurprising to see that the cites where the PM 2.5 showed a significant response to fireworks were in the majority during the 2013-2016 celebration periods in the Chinese Lunar New Year, amounting to 74.2%, 61.3%, 51.6%, and 58.1% of the total, respectively. In contrast, the proportion of cities with a weak response was almost the lowest, amounting to 6.5%, 9.7%, 19.4%, and 25.8%, respectively. Notably, only two cities fell into the weak response category in 2013. Even taking different diffusion conditions into account, it is still reasonable to say that the discharge of Chinese Lunar New Year fireworks is causing significant instantaneous PM 2.5 increases. More importantly, such a phenomenon does not only occur in individual cities, but in 31 cities that are across the whole country. 0:00-04:00 period of the remaining days, it was considered that the city's PM2.5 responded obviously to the New Year fireworks. (3) If the variation of PM2.5 in a city did not match the above two conditions, it was considered that the city's PM2.5 responded weakly to the New Year fireworks.
To avoid any possible interference, the data from the Lantern Festival period were excluded from this classification. In addition, according to the Tibetan calendar (Table S4) and related customs, the data from the Tibetan New Year firework period would not affect the classification results.
According to the classification map shown in Figure 5 and the corresponding statistical results (Table S5), it is unsurprising to see that the cites where the PM2.5 showed a significant response to fireworks were in the majority during the 2013-2016 celebration periods in the Chinese Lunar New Year, amounting to 74.2%, 61.3%, 51.6%, and 58.1% of the total, respectively. In contrast, the proportion of cities with a weak response was almost the lowest, amounting to 6.5%, 9.7%, 19.4%, and 25.8%, respectively. Notably, only two cities fell into the weak response category in 2013. Even taking different diffusion conditions into account, it is still reasonable to say that the discharge of Chinese Lunar New Year fireworks is causing significant instantaneous PM2.5 increases. More importantly, such a phenomenon does not only occur in individual cities, but in 31 cities that are across the whole country.

Sustained Effect of Fireworks Discharge on PM2.5
Since the transient effect is so strong, we could not help but wonder how strong the sustained effect is. As can be seen in Figure 4a and Figure S3, although the PM2.5 on New Year's Day showed a fairly sharp increase, the PM2.5 concentrations dropped back to the usual level over the next few days.

Sustained Effect of Fireworks Discharge on PM 2.5
Since the transient effect is so strong, we could not help but wonder how strong the sustained effect is. As can be seen in Figure 4a and Figure S3, although the PM 2.5 on New Year's Day showed a fairly sharp increase, the PM 2.5 concentrations dropped back to the usual level over the next few days. In fact, there is usually a great variety of spontaneous celebrations and commemorations over the entire 17-day celebration period, which often involve the discharge of fireworks. Therefore, it was necessary to explore the continuous influence of New Year's fireworks on PM 2.5 throughout the entire celebration period. To undertake this analysis, each Chinese lunar year was divided into 21 periods by the same time span as the celebration period through the 17-day interval division scheme. The average national PM 2.5 monitoring values during each period in 2013-2016 were calculated and are shown in Figure 6. celebration period. To undertake this analysis, each Chinese lunar year was divided into 21 periods by the same time span as the celebration period through the 17-day interval division scheme. The average national PM2.5 monitoring values during each period in 2013-2016 were calculated and are shown in Figure 6. As can be seen, the phased concentrations of PM2.5 is in accordance with the seasonal variation appeared in Figure 4, and interestingly, even under the dramatic transient effect of the fireworks, the celebration period PM2.5 concentration does not break the seasonal variation pattern. The PM2.5 of four celebration periods were 87.43 μg/m³, 83.37 μg/m³, 51.13 μg/m³, and 56.68 μg/m³, respectively. The corresponding PM2.5 averages of pre-and post-celebration periods were 94.71 μg/m³, 96.46 μg/m³, 67.45 μg/m³, and 57.68 μg/m³, respectively. The PM2.5 in four celebrations were all below the average level of pre-and post-celebration periods, and the PM2.5 level was even lower than that of both the pre-and post-celebration periods in 2014, 2015, and 2016. In addition, although the PM2.5 level during the 2013 celebration period was slightly higher than that of the post-celebration period, it was still much lower than that of the pre-celebration period. It therefore seems that the sustained effect of the New Year fireworks is not as strong as the instantaneous effect.  As can be seen, the phased concentrations of PM 2.5 is in accordance with the seasonal variation appeared in Figure 4, and interestingly, even under the dramatic transient effect of the fireworks, the celebration period PM 2.5 concentration does not break the seasonal variation pattern. The PM 2.5 of four celebration periods were 87.43 µg/m 3 , 83.37 µg/m 3 , 51.13 µg/m 3 , and 56.68 µg/m 3 , respectively. The corresponding PM 2.5 averages of pre-and post-celebration periods were 94.71 µg/m 3 , 96.46 µg/m 3 , 67.45 µg/m 3 , and 57.68 µg/m 3 , respectively. The PM 2.5 in four celebrations were all below the average level of pre-and post-celebration periods, and the PM 2.5 level was even lower than that of both the preand post-celebration periods in 2014, 2015, and 2016. In addition, although the PM 2.5 level during the 2013 celebration period was slightly higher than that of the post-celebration period, it was still much lower than that of the pre-celebration period. It therefore seems that the sustained effect of the New Year fireworks is not as strong as the instantaneous effect.
To further establish the sustained effect of New Year's fireworks over a longer period of time, Terra and Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data were employed. We used the 10-km AOD products to retrieve the near-ground PM 2.5 data from 2002 to 2016, and further combined this with the 17-day interval division scheme to obtain fusion maps of PM 2.   Figure 7, we can see the retrieved PM 2.5 showed similar spatial distributions to the ground observations. According to previous studies [82], the PM 2.5 level in northern China was generally higher, and the PM 2.5 level in the southern area was relatively lower, which is mainly affected by population distribution and urbanization progress.
As is known to us all, the celebration period is usually in the middle of winter when PM 2.5 pollution is at its worst. As in Figure 7 we can see that the PM 2.5 does show a significant pattern of seasonal variation, it decreases during the spring (from the 1st post-celebration period to the 5th post-celebration period), and reaches the lowest point in the summer (from the 6th post-celebration period to the 10th pre-celebration period), then gradually increases through the autumn (from the 9th pre-celebration period to the 5th pre-celebration period), and finally reaches the highest point in the winter (from the 4th pre-celebration period to the celebration period). The regional aggregations indicated that there are many other factors that play important roles in PM 2.5 accumulation. For example, the PM pollution in Xinjiang is possibly due to the existence of deserts, while the PM 2.5 hot spots in Sichuan and Beijing-Tianjin-Hebei region are not only affected by population density, but also because of the surrounding mountains make it difficult for the pollutants to spread [83].
As for smaller time window, compared with the pre-celebration period, PM 2.5 in the celebration period shows a remarkable decrease, especially in areas with higher PM 2.5 concentrations, such as the Beijing-Tianjin-Hebei region and the Sichuan Basin. In the post-celebration period, the pollution is further alleviated, which is more prominent in central and eastern China. Although the station-based overall PM 2.5 level of the celebration period shows a gradual decreasing trend which is consistent with its seasonal variation, it is still unexpected to see that the PM 2.5 of the entire celebration period shows no obvious increase across the country.

Effect of Firework Prohibition Policies
Since the debate is still fierce over whether to ban fireworks discharge to improve air quality, some local governments have attempted to implement such policies in China as attempts [84]. As the statistical results showed in Table S5, the number of cities with a weak response showed a certain increase throughout 2013-2016, which is very likely to be a policy-oriented phenomenon. It is therefore necessary to analyze and evaluate the effects of these policies.
Considering many factors, we decided to take the cities of Wuhan, Chengdu, and Zhengzhou as examples. These cities implemented firework prohibition policies in 2014, 2015, and 2016, respectively, which makes the data of these cities very suitable for comparing the changes in PM 2.5 levels before and after the policies' implementation.
The 34-day mean PM 2.5 concentration during the pre-celebration period (17 days) and the post-celebration period (17 days) were used to represent the normal PM 2.5 level, and the difference between the hourly PM 2.5 concentrations and the corresponding averages was then calculated to indicate the policy efficiency.
From Figure 8, we can see when focusing on the variation of PM 2.5 pollution in the Chinese Lunar New Year's Eve and the early morning of New Year's Day, it is clear that the implementation of the firework prohibition policies significantly reduced the instantaneous PM 2.5 increase. For example, PM 2.5 concentration sharply increased in Wuhan on New Year's Day in 2013, but after the implementation of the prohibition policy in 2014, the PM 2.5 concentration decreased significantly in the next three celebration periods. For Chengdu, the PM 2.5 concentrations in the last two years, under the prohibition policy, were also lower than those of the previous two years, despite the remaining visible peaks. As for Zhengzhou, the PM 2.5 concentrations were all very high in 2013-2015, but the contamination level was greatly reduced in 2016 when the prohibition policy took effect. continual PM2.5 increase, and the reduction of pollutants by firework discharge regulations were not effective enough. Similarly, the effect of daily pollution can also be observed in Wuhan and Chengdu from 2013 to 2016. For most of the time after the seventh day of the Lunar New Year, the PM2.5 was lower than the normal level. It is especially apparent for Chengdu in 2013, where the PM2.5 level was almost always lower than usual, except for the early hours of New Year's Day.

Discussion
The results of this study indicate that the intense discharge of fireworks at the Chinese New Year can dramatically increase PM2.5 over a short period of time, and this effect exceeds that of any other factor nationwide. However, the sustained effect of the fireworks is not strong. On the one hand, the discharge of fireworks at New Year is intense, but the duration is short, and when combined with certain meteorological conditions [86], it allows the PM2.5 to diffuse and settle quite quickly in the relatively less-polluted surroundings. On the other hand, during the celebration period, some major However, when analyzing the whole of the Chinese Lunar New Year celebration period, it is difficult to say if the prohibition policy is still that effective. This is likely because other factors such as meteorological conditions and daily emissions can also profoundly influence the pollution level of PM 2.5 . For example, the purification effect of heavy precipitation made the air much cleaner in 2015 than in 2014 and 2016 in Wuhan, although prohibition policies were in effect in all these three years (see Table S6). In Chengdu in 2013, although the fireworks discharge significantly increased the instantaneous PM 2.5 concentrations, the favorable weather conditions still made the PM 2.5 level in the few days around the New Year's Day lower than those in the years when the prohibition policies were implemented (see Table S7). In contrast, in 2014, the stable meteorological conditions in Chengdu greatly exacerbated the accumulation of pollutants, resulting in this city suffering from much more severe PM 2.5 pollution [85]. Another example about industrial and municipal discharges is in 2016, the New Year peak of PM 2.5 in Zhengzhou was noticeably eased under the implementation of the prohibition policy, but under the same weather conditions (see Table S8), the PM 2.5 pollution situation began to deteriorate from the second day of the new year, and then the serious pollution lasted for four days. This phenomenon was probably because daily emissions were still causing a continual PM 2.5 increase, and the reduction of pollutants by firework discharge regulations were not effective enough. Similarly, the effect of daily pollution can also be observed in Wuhan and Chengdu from 2013 to 2016. For most of the time after the seventh day of the Lunar New Year, the PM 2.5 was lower than the normal level. It is especially apparent for Chengdu in 2013, where the PM 2.5 level was almost always lower than usual, except for the early hours of New Year's Day.

Discussion
The results of this study indicate that the intense discharge of fireworks at the Chinese New Year can dramatically increase PM 2.5 over a short period of time, and this effect exceeds that of any other factor nationwide. However, the sustained effect of the fireworks is not strong. On the one hand, the discharge of fireworks at New Year is intense, but the duration is short, and when combined with certain meteorological conditions [86], it allows the PM 2.5 to diffuse and settle quite quickly in the relatively less-polluted surroundings. On the other hand, during the celebration period, some major sources of daily PM 2.5 pollution, such as factory and traffic emissions, are greatly reduced due to factory shutdown, festive population migration, and other related factors. Thus, the "holiday effect" may also play an important role in curbing PM 2.5 pollution [45,81]. Therefore, it is rational to believe that the aggravation of PM 2.5 pollution caused by the Chinese Lunar New Year fireworks discharge is no stronger than the mitigation brought by the natural cleaning processes and the "holiday effect".
As for the effectiveness of firework prohibition policies, it is clear that these policies can significantly alleviate the PM 2.5 bursts caused by fireworks discharge in the early hours of Lunar New Year's Day. However, we have to admit that its influence has limitation over the whole celebration period, which is probably because there are still many other strong influence factors, such as meteorological conditions and daily emissions, that can offset or even overshadow the effect of policies. Additionally, it is important to take the cultural and symbolic meaning of Chinese Lunar New Year fireworks into account, in that people may still want to celebrate the Chinese Lunar New Year in the traditional way, no matter what [87]. During 2014 and 2015 in the Chinese Lunar New Year celebration periods in Chengdu, the PM 2.5 level still showed an obvious increase in the early morning of New Year's Day, despite the prohibition policy being in effect. Additionally, during the 2014 celebration period in Wuhan, even if they may get punished, people still celebrated the New Year in the traditional way [88].
Firework prohibition policies can indeed alleviate the transient pollution, but their effect on the sustained PM 2.5 concentration during the entire celebration period of the Chinese Lunar New Year is not obvious. As a Chinese culture passed down for thousands of years, the custom of discharging fireworks during the Chinese Lunar New Year is deeply rooted. It therefore seems overly simplistic and unreasonable to completely prohibit this celebration activity only for the air pollution consideration. More appropriate approaches should be adopted in a more humanized and scientific way. For example, the government could consider guiding the public to rationally discharge fireworks at appropriate times and places, according to diffusion conditions and safety factors. Alternatively, we could encourage the development of low-pollution substitutes, including replacing the traditional high-polluting gunpowder with environmentally friendly desulfurized explosives, and adopting electronic fireworks products.

Conclusions
To date, there have been some studies of PM 2.5 variation with firework displays, but most of the works have focused on individual cities [40,42,45,48,49,89]. However, considering that China has a large span of latitude and longitude, the climatic characteristics and geographical conditions are quite different, which makes the researches based on a single city have certain limitations. Therefore, a comprehensive analysis of the effect of Chinese Lunar New Year fireworks discharge on PM 2.5 concentration in China was necessary.
In this study, we jointly used station monitoring data and remote sensing products to comprehensively analyze the effects of the Chinese Lunar New Year fireworks discharge on PM 2.5 concentration within the entire mainland China. By using remote sensing data, the problem of uneven spatial distribution of ground monitoring has been compensated, and the research time span has been expanded from four years (ground observation time) to fifteen years (satellite observation time).
The results indicate that the intensive Chinese Lunar New Year fireworks discharge can contribute to a sharp increase of PM 2.5 within a short period of time. Compared with the concentrations before and after each Lunar New Year, the PM 2.5 of the four Lunar New Years increased by 109%, 59%, 123%, and 106%, respectively. In most cities (ranging from 51.6~74.2%), PM 2.5 variations during the Lunar New Year showed strong correlations with the discharge of fireworks. However, due to the suitable diffusion conditions and the reduction of factory and traffic emissions, the overall PM 2.5 concentration during the 17-day celebration period of the Chinese Lunar New Year is no higher, and even lower than the normal level (which were 7.28 µg/m 3 , 13.09 µg/m 3 , 16.32 µg/m 3 , and 0.99 µg/m 3 lower compared to normal level, respectively). Therefore, the sustained effect of fireworks discharge on PM 2.5 concentration is not as serious as many people may once have thought. As for the firework prohibition policies, its implementation has been proved to have a significant effect on the reduction of the instantaneous PM 2.5 peaks. The PM 2.5 increase during the concentrated fireworks discharge hours (00:00-04:00 on Lunar New Year's Day) has been significantly reduced. However, its impact on the entire celebration period is limited. The PM 2.5 levels did not show any significant changes during the rest of the celebration period before and after the implementation of the policies. At the same time, other factors may also have a certain impact on the effect of policy implementation.
The custom of setting off fireworks during the Chinese Lunar New Year has been passed down for thousands of years, and it seems too arbitrary to completely abolish this custom only for air pollution consideration. We believe that this significant custom should be continued with minimized negative effects, with the assistance of new technology and scientific management.