Air pollution during New Year’s fireworks and daily mortality in the Netherlands

Short-term exposure to air pollution has been associated with cardiovascular and respiratory mortality and morbidity. Little is known about associations between air pollution caused by firework events and daily mortality. We investigated whether particulate matter from fireworks during New Year’s celebrations was associated with daily mortality. We analyzed the celebrations of the years 1995–2012. PM10 concentrations increased dramatically during the firework events. Countrywide, the daily average PM10 concentrations from 27–30 December was 29 μg/m3 and increased during the first hour of the New Year by 277 μg/m3. In the more densely populated areas of the Netherlands the increase was even steeper, 598 μg/m3 in the first hour of the New Year. No consistent associations were found using linear regression models between PM10 concentrations during the first six hours of 1 January and daily mortality in the general population. Yet, using a case-crossover analysis firework-days and PM10 concentrations were associated with daily mortality. Therefore, in light of the contradictory results obtained with the different statistical analyses, we recommend further epidemiological research on the health effects of exposure to firework emissions.

inventory of diagnoses made on patients admitted to a hospital in Philadelphia during the week of July 4, revealed two cases of asthma exacerbation, one fatal and one near-fatal, following exposure to elevated PM concentrations from fireworks 32 . Smith et al. performed a small-scale study (n = 9) around New Year's festivities in which they found a decrease in pulmonary function in 2 volunteers with a history of respiratory disease, while the 7 subjects without a history of respiratory disease showed no significant change following exposure to firework emissions 33 .
An increase in emergency room visits following a fireworks episode was described by Bach et al. 34 . Furthermore, Beig et al. 4 estimated an increase in mortality and morbidity attributed to population exposure to PM 2.5 and PM 10 mass concentrations within areas of 2 kilometers radii from the fireworks displays. Godri et al. 9 described in an in vitro study a relationship between the oxidative potential of PM and trace metals associated with fireworks, suggesting a potential negative impact of fireworks emissions on health. In this study, we evaluated the association between hourly PM 10 concentrations, observed during New Year's Eve fireworks, and daily mortality in the Netherlands. During the research period setting off fireworks was permitted from 31 December 10 am until 1 January 2 am.

Results
Peak exposures were found between midnight and 1 am.
Mortality. Summary statistics of daily mortality, PM 10 concentrations, and daily temperatures are presented in Table 1.
Average non-accidental daily mortality decreased from 478 in 1995 to 383 in 2011 (Fig. 1). The decline was mainly found for individuals over 65 years of age (395 in 1995 to 309 in 2011) and for cardiorespiratory mortality (260 in 1995 to 165 in 2011). Daily non-accidental and cardiorespiratory mortality was slightly higher in January as compared to December (Table 1).
Air quality. The average PM 10 concentration from 27-30 December was 29 μg/m 3 (Table 1) and slightly decreased (Fig. 2) over the years. PM 10 concentrations increased during the first hour of the New Year on average in the less densely populated Dutch municipalities by 143 μg/m 3 (range, 35-255 μg/m 3 ), and in the densely populated municipalities by 598 μg/m 3 (range, 335-1132 μg/m 3 ). PM 10 concentrations during the first hour of 1 January showed no clear pattern over the years (Fig. 2). Figure 2 does not show the PM 10 concentrations of the 'densely populated' and the 'less densely populated' municipalities in December separately, because these differences are so small that they would be indiscernible in this figure.  www.nature.com/scientificreports www.nature.com/scientificreports/ Associations between daily mortality and PM 10 from fireworks. Although non-accidental mortality slightly increased following exposure to increased PM 10 concentrations on exposure days and lag days, none of the associations between the increase in PM 10 concentrations during the first 6 hours of 1 January and mortality were statistically significant ( Table 2). Cardiorespiratory mortality showed some increases and decreases and as with non-accidental mortality none of these associations were statistically significant ( Table 2). Results for the increase in PM 10 concentrations during the first hour (0-1) and the first four hours (0-4) were comparable (Additional file 1: Table S1 and Table S2).
Results for the different age groups are shown in Figs 3 and 4. Results for the densely and the less densely populated regions separately are presented in the Supplemental material (Additional file 1: Table S3 and Table S4). In the densely populated regions no associations between PM 10 concentrations during the first hours of 1 January and daily non-accidental and cardiorespiratory mortality were found. In the less densely populated areas a 10 μg/m 3 increase in PM 10 concentration from 0-6 am 1 January, increased the risk of non-accidental mortality for 0-65 years on 1 January by 1.68% (95%CI: 0.64%, 2.72%). Adjusted for daily temperature the increase was 1.53% (95%CI: 0.36%, 2.70%) as presented in the Supplemental material (Additional file 1: Table S5 and Table S6). No other associations with and without adjustment for daily temperature were found for other age groups.
The case-crossover analysis showed both positive associations between firework-events and mortality, and between PM 10 concentrations and mortality ( Table 3). The same pattern of associations was found amongst the age group over 65 years (Additional file 1: Tables S8 and S9). However, if we put firework-events and PM 10 levels in the model together, associations strongly decreased (Additional file 1: Table S10).

Sensitivity analyses of the linear regression analysis.
No associations between PM 10 concentrations and mortality were found in the sensitivity analyses for monitoring stations that had data for at least 15 years (results not shown). The results of the sensitivity analyses in which we excluded specific years are given in the supplement (Additional file 1: Table S7). The results of these sensitivity analyses did not show consistent positive associations between PM 10 and non-accidental and cardiorespiratory mortality.  Table 2. Mean percent change in daily mortality associated with 10 µg/m 3 difference in PM 10 concentration on 1 January 0-6 hours compared to the pre-firework concentration.

Discussion
In general, linear regression analysis revealed no statistically significant associations between PM 10 concentrations during the first hours of 1 January and daily non-accidental and cardiorespiratory mortality on the first 4 days of January. Additionally, no associations were found in the age group over 65 years. Yet, in a symmetric bi-directional case-crossover analysis we did find positive associations between PM 10 concentrations and mortality, and between firework-events and mortality using PM 10 concentrations and firework-events separately. Additionally, the same pattern was found in the age group over 65 years. However, associations considerably weakened when both variables were put in the model together.    Table 3. Case-crossover analysis in The Netherlands: ORs and 95% CIs for daily mortality associated with fireworks adjusted for temperature and daily mortality associated with 10 µg/m 3 PM 10 concentration on hour 0-6 adjusted for daily temperature. www.nature.com/scientificreports www.nature.com/scientificreports/ Based on the results of the linear regression analysis, effect estimates for non-accidental mortality increase per 10 µg/m 3 PM 10 in our study were between 0.1% and 0.2% whereas the non-accidental increases reported earlier in the Netherlands were between 0.3 and 0.8% 24,35 . The increases reported in international meta-analyses were between 0.36% and 0.48% 36,37 . Effect estimates are also lower than the increase (0.49%) associated with PM 10 from wildfire exposures 23 .
The absence of associations in the linear regression analysis of this study that contrast with the associations found in the abovementioned studies 23,24,35,37 , are potentially due to the short exposure duration to these peak concentrations. Consequently, adverse health effects might be limited, as has been suggested earlier 10 . The peak PM 10 concentrations lasted for a very short time, found on 1 January from 0-1 hours after which the concentrations quickly leveled off during the following hours as supported by other studies 10,19 . Therefore, daily averaged concentrations PM 10 on 1 January were lower than the reported 6-h concentrations. Because the results of the analyses with the peak concentrations from 0-1 am and the concentrations of the 6-h time period did not differ, we reported mainly the analyses of the 6-h period to allow for a better comparison to the literature. On the other hand, using the case-crossover approach we did find associations between PM 10 concentrations and mortality. Because PM 10 concentrations are strongly associated with New Year's celebrations, the found associations might also be caused by other putative determinants that are typical for New Year's celebrations, such as emotional stress, changes in diet and alcohol consumption, and so forth 38 . Therefore, we additionally performed a case-crossover analysis using the Christmas days during our study period as case days. Christmas has many of the abovementioned determinants in common with New Year's celebrations, with the exception of elevated PM 10 concentrations. No associations were found between Christmas days and mortality, while we did find an association with New Year's celebrations. A putative explanation of the difference between Christmas and New Year might be that air pollution due to fireworks caused increased daily mortality.
There are potential alternative explanations for the negative results of the linear regression analysis in our study. First, the data set was limited to 17 periods (i.e. 17 data points) because air pollution by fireworks was restricted to one episode per year (1996-2012). Furthermore, around 75% of the Dutch inhabitants live in the less densely populated municipalities, in which the air quality is relatively less affected by the fireworks than in the densely populated municipalities. However, in the densely populated areas no significant positive associations were found between PM 10 and daily non-accidental mortality. Therefore, our study may have had insufficient power. Second, susceptible groups such as people with asthma and the elderly probably reduce the exposure to fireworks emissions by staying inside during the peak of firework events 39,40 . Third, the Christmas season, which includes New Year's Day, is an annually recurring period characterized by days off, eating, drinking, celebrating and stress. The influence of these factors on daily mortality might obscure potentially existing associations with air quality. However, this is not supported by the results from the case-crossover analyses which suggest that factors such as eating, drinking and celebrating might be insufficient to totally obscure the impact of particulate matter on health. The case-crossover design has the advantage of controlling for potential confounders caused by time-invariant individual characteristics (e.g. age, sex, body mass index, and comorbidity) [41][42][43][44] . The results of a simulation study by Bateson and Schwartz showed that the symmetric case-crossover design performed best in terms of bias 41 . Unlike other studies, the results in our study obtained with a case-crossover analysis were not similar to those obtained with a different statistical analysis 44,45 . At present, we cannot give insight in what may have caused the difference between the results obtained with the case-crossover analysis and the linear regression analysis, respectively.
A limitation of the study was that we only adjusted for daily temperature. Given the low number of data points we did not adjust for relative humidity and influenza. Furthermore, measurements did not discriminate between fireworks and other sources of particulate matter. However, because of the steep increase of PM 10 within a short period of time, it is most likely that the peak concentrations on 1 January consisted largely of particulate matter from fireworks. Moreover, there are no indications that other human activities, which emit particulate matter, increase during setting of fireworks 46 .
In this study, PM 10 was chosen as a proxy for air pollution by fireworks because of the relative abundance of the monitoring data. Although PM 10 is just one of the pollutants emitted by setting off fireworks, others have shown that the changes in PM 10 levels correspond with changes in PM 2.5 , PM 1 , carbon monoxide, sulphur dioxide and other pollutants [2][3][4]13,17,46,47 . Therefore, we think that the use of PM 10 instead of PM 2.5 or other pollutants will not have led to exposure misclassification.
To date only a handful of small scale studies have addressed the potential adverse health effects stemming from exposure to firework emissions. In 2016 Lin reviewed evidence related to ambient particulate matter during firework periods and associated human health risks 48 . According to Lin, among the 49 reviewed articles, only 7 have reported health risk evaluations directly related to particulate matter from fireworks. Smith and Dinh described a 26% decrease in maximal midexpiratory flow rate (FEF 25-75% ) following exposure to air pollution from fireworks in two subjects with a history of chronic respiratory disease. No significant decrease was found in the FEF 25-75% of seven healthy subjects 33 . Do et al. found that post-firework particulate matter was more toxic to human bronchial epithelial BEAS-2B cells than pre-firework particulate matter 49 . The other reviewed studies described quantifications of adverse health impacts caused by exposure to particulate matter 4,50 or associated metals 16,18,51 . Furthermore, two case studies reported a case of acute eosinophilic pneumonia and a fatal, and a near-fatal asthma exacerbation of two asthmatic children following exposure to elevated PM concentration from fireworks 31,32 . Bach et al. described an increase in emergency room visits following a fireworks episode 34 . Godri et al. 9 described in an in vitro study a relationship between the oxidative potential of PM and trace metals associated with fireworks, suggesting a potential negative impact of fireworks emissions on health. (2019) 9:5735 | https://doi.org/10.1038/s41598-019-42080-6 www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusion
As far as we know, this is the first large-scale observational epidemiological study on the association between exposure to air pollution by firework events and mortality. Although positive associations were absent in the linear regression analysis, the case-crossover analysis showed some positive associations between firework-events, PM 10 concentrations and daily mortality. As for now, we do not know why the different analysis methods resulted in dissimilar results. Therefore, in light of the contradictory results obtained with the different statistical analyses, we recommend further epidemiological research on the consequences of exposure to firework emissions.
Mortality data were stratified in age groups 0-65 years and older than 65 years. Furthermore, data were stratified according to population density ('densely populated': municipalities with ≥2500 inhabitants/km 2 and 'less densely populated': all other municipalities in the Netherlands).
Air quality. PM 10 data were obtained from the National Institute for Public Health and the Environment (RIVM), which operates the Ambient Air Quality Monitoring Network in the Netherlands and used as a proxy for air pollution by fireworks. The network had 46 monitoring stations for PM 10 in both densely and less densely populated areas during the period 1995-2012 53 .
For the entire period, concentration differences between 1-h average concentrations on 1 January and the mean of 1-h average concentrations over 27-30 December of the preceding year were calculated for each monitoring station. Besides the 1-h average concentrations from 0-1 am, the 4-hours averaged concentrations from 0-4 am and the 6-hours averaged concentrations from 0-6 am on 1 January were used.
Imputation of one or more missing values between 1 am and 6 am on 1 January was based on linear interpolation of preceding and following 1-h average concentrations per monitoring station. Missing values between midnight and 1 am were not imputed, because maximum levels were found from 0-1 am.
Meteorology. Mean daily temperature was obtained from the Royal Dutch Meteorological Institute (KNMI) for a weather station at a central rural location (de Bilt, 05° 17′ 70″, 52° 10′ 10″) of the national meteorological network 54 .
Data analysis. Linear regression analysis. SPSS version 20.0 statistical software was used (SPSS Inc., Chicago, IL, USA). Associations between exposure variables and mortality with and without adjustment for daily temperature were calculated using a linear regression analysis. The level of statistical significance was set at p < 0.05.
Daily mortality differences between 1 January, 2 January, 3 January, 4 January, the average over 1-4 January and the average daily mortality over 27-30 December of the preceding year were analyzed in association with PM 10 concentration differences (per 10 µg/m 3 ) between the 1-h average concentrations during which the peak concentration occurs on 1 January (0-1 am and 0-6 am) and the 1-h average concentrations over 27-30 December of the preceding year. Regression coefficients were recalculated to mean percent increase (95% CI) in daily mortality per 10 µg/m 3 PM10 concentration on 1 January 0-6 hours. Therefore, data were restricted to 17 periods of four days before (27)(28)(29)(30) and four days after (1-4 January) New Year. Stratified data for different age groups and population densities were analyzed separately. For the analyses stratified by population densities we used PM 10 data for stations in densely populated and less densely populated areas separately. For each analysis, we calculated the percent change in mortality per 10 µg/m 3 increase in PM 10 by dividing the regression coefficient for PM 10 by the baseline daily mortality (i.e. the age-group-specific and cause-specific average daily mortality on 27-30 December during the entire study period).
Case-crossover analysis. Additionally, a case-crossover analysis 55,56 was used in which each deceased individual served as her or his own control. We used the symmetric bi-directional approach 41 , in which the 1 January was selected as case day (see first column) and 4 control days were chosen per case day (the same days of the week in the 2 preceding weeks and in the 2 successive weeks). PM 10 concentrations were calculated as the mean of the first six hours per day. Using conditional logistic regression, we estimated the OR in three models for mortality associated with 1) firework-days adjusted for temperature, 2) 10 µg/m 3 PM 10 concentration adjusted for temperature, and 3) 10 µg/m 3 PM 10 concentration adjusted for temperature and firework-day (yes/no). For both the case day and PM 10 concentrations four lags were evaluated (lags 0, 1, 2, and 3 days).
Sensitivity analysis of the linear regression analysis. Sensitivity analyses were performed by analyzing only monitoring stations that had data available during the entire measurement period (n = 8), and by analyzing stations that had data available for at least 15 of the 17 years (n = 16). We also performed sensitivity analyses to investigate if the results changed when we excluded: a. the end-of-year-period with the highest level of PM 10 ; b. the end-of-year-period with the lowest level of PM 10 and c. the two end-of-year-periods with the highest level of PM 10 (there were 2 periods with much higher levels than the other years (see Fig. 2)).