Simulating PM2.5 Concentrations during New Year in Cuenca, Ecuador: Effects of Advancing the Time of Burning Activities

Fine particulate matter (PM2.5) is dangerous to human health. At midnight on 31 December, in Ecuadorian cities, people burn puppets and fireworks, emitting high amounts of PM2.5. On 1 January 2022, concentrations between 27.3 and 40.6 µg m−3 (maximum mean over 24 h) were measured in Cuenca, an Andean city located in southern Ecuador; these are higher than 15 µg m−3, the current World Health Organization guideline. We estimated the corresponding PM2.5 emissions and used them as an input to the Weather Research and Forecasting with Chemistry (WRF-Chem 3.2) model to simulate the change in PM2.5 concentrations, assuming these emissions started at 18:00 LT or 21:00 LT on 31 December 2021. On average, PM2.5 concentrations decreased by 51.4% and 33.2%. Similar modeling exercises were completed for 2016 to 2021, providing mean decreases between 21.4% and 61.0% if emissions started at 18:00 LT. Lower mean reductions, between 2.3% and 40.7%, or even local increases, were computed for emissions beginning at 21:00 LT. Reductions occurred through better atmospheric conditions to disperse PM2.5 compared to midnight. Advancing the burning time can help reduce the health effects of PM2.5 emissions on 31 December.


Introduction
Fireworks to welcome in the new year are customary across different parts of the world [1,2]. Around midnight on 31 December in Ecuador, it is also a tradition to burn the "años viejos", cloth dolls filled with sawdust or paper. Both the fireworks and the burning of puppets take place in a festive environment amid social and family connections [3]. The use of fireworks results in raising short-term particulate matter concentrations [4], which affects public health and visibility [5], reaching levels that typically exceed both national regulations and the guidelines of the World Health Organization (WHO). Most of the particles aerosolized by burning fireworks consist of metals and compounds in colorful firework displays [2]. PM 2.5 produces health effects after both short-term and long-term exposure. The WHO guideline does not guarantee complete protection against PM 2.5 effects [6], requiring the lowest concentrations to be achieved. Short-term exposure promotes cardiovascular and respiratory effects [7]. The International Agency for Research on Cancer (IARC) classified particulate matter and outdoor air pollution as carcinogenic to humans [8,9]. Moreover, in terms of effects on the brain, particulate matter appears to be the most concerning air pollutant [10]. PM 2.5 accumulates in the reproductive organs, disrupting hormone levels and affecting fertility [11]. Moreover, pyrotechnic fireworks introduce significant quantities of toxic metals into the atmosphere, affecting mammalian cells and lungs [12]. PM 2.5 affects red blood cell counts, monocyte counts, and hemoglobin concentration [13]. Fireworks have increased health problems, particularly in infants, women (including pregnant women),

•
Are there reductions in PM 2.5 concentrations if the emission activities associated with the arrival of a new year in Cuenca occur some hours before midnight? • What is the magnitude of these changes?

Primary Sources and Past PM 2.5 Concentrations in Cuenca
Cuenca (2500 m.a.s.l.), an Andean city located in southern Ecuador (Figure 1), has a complex topography, with altitudes ranging from 1000 to 4000 m.a.s.l. This city has approximately 640,000 inhabitants, representing approximately 3.6% of the Ecuadorian population [33]. This territory has diverse land-use categories. The primary sources of PM 2.5 are diesel vehicles, industries at the northeast zone of the city, and a power facility located in the northeast of the urban area ( Figure 1d). Moreover, irregular PM 2.5 emissions come from artisanal brick production in the northwest of the urban region, due to the combustion of biomass [34]. More information about the magnitude of the emission sources in Cuenca are described in Parra (2018) [35] and Parra and Espinoza (2020) [36].

Primary Sources and Past PM2.5 Concentrations in Cuenca
Cuenca (2500 m.a.s.l.), an Andean city located in southern Ecuador (Figure 1), has a complex topography, with altitudes ranging from 1000 to 4000 m.a.s.l. This city has approximately 640,000 inhabitants, representing approximately 3.6% of the Ecuadorian population [33]. This territory has diverse land-use categories. The primary sources of PM2.5 are diesel vehicles, industries at the northeast zone of the city, and a power facility located in the northeast of the urban area ( Figure 1d). Moreover, irregular PM2.5 emissions come from artisanal brick production in the northwest of the urban region, due to the combustion of biomass [34]. More information about the magnitude of the emission sources in Cuenca are described in Parra (2018) [35] and Parra and Espinoza (2020) [36]. During the first hours of 2016 to 2021, the air quality network in the MUN (Municipio) station measured maximum hourly PM2.5 concentrations between 55.8 and 182.1 µg m −3 . The MUN station has been operating since 2012, measuring both meteorological parameters and air quality. Some years later, the Colegio Carlos Arízaga (CCA) and Escuela Ignacio Escandón (EIE) stations were installed to monitor air quality in the zones of the industrial park and in the southwest of the city (Figure 1d). In the last weeks of 2021, three During the first hours of 2016 to 2021, the air quality network in the MUN (Municipio) station measured maximum hourly PM 2.5 concentrations between 55.8 and 182.1 µg m −3 . The MUN station has been operating since 2012, measuring both meteorological parameters and air quality. Some years later, the Colegio Carlos Arízaga (CCA) and Escuela Ignacio Escandón (EIE) stations were installed to monitor air quality in the zones of the industrial park and in the southwest of the city (Figure 1d). In the last weeks of 2021, three new stations were installed (Condamine, CON; Terminal Terrestre, TER; Cebollar, CEB), equipped with sensors to record PM 2.5 and PM 1 . The air quality network is operated by the Municipality of Cuenca, the entity accredited by the National Environmental Authority.

Estimation of PM 2.5 Emissions through Burning Activities Associated with the Arrival of 2022
In general, emission inventories have high levels of uncertainty. This feature is especially critical when estimating New Year's emissions. After exploring different approaches, we estimated the total PM 2.5 emissions at the beginning of 1 January 2022 as the contribution of the use of fireworks and the burning of puppets, through the basic emission model: Activity: amount of fireworks or puppets burned during the New Year's festivities. Emission factor: amount of PM 2.5 emitted by each firework or puppet burned.

Emissions Due to the Use of Fireworks
Firstly, we estimated the amount of fireworks used in Cuenca, through the following approach.
Based on information about the amount of money spent for the importation of fireworks in 2021 [37,38], we estimated that in Cuenca, approximately USD 112,000 was allocated for the purchase of fireworks in December. We identified the pyrotechnic cake with 30 shots, which has approximately 500 g of net explosive content (NEC), as the most consumed product. Using an average price of USD 24 per cake, we estimated that there was 3497 kg NEC in the amount of fireworks legally sold. Moreover, we considered the additional amount of fireworks from illegal sales by non-controlled producers or importers. Based on information about illegal fireworks [39] and an assumed percentage of noncontrolled sales, we estimated a 3279 kg NEC as a result of these sales. Therefore, the total amount of fireworks were estimated to equal 6776 kg NEC.
Then, we selected the average emission factor of 200 g PM 2.5 NEC −1 proposed by Keller and Schragen (2021) [40] for common pyrotechnic articles. Applying Equation (1), the corresponding emission was 1.36 t PM 2.5 .

Emissions Due to the Burning of Puppets
The amount of puppets burned was estimated through the following approach. The Municipality of Cuenca determined 369 authorized points for selling puppets, with an average amount of 25 units sold per day [41]. Considering the number of sales over two days, an average weight of 2 kg per unit, and an increase of 25% due to domestic production and non-authorized sales, we estimated 46.1 t as the amount of biomass contained in the puppets. Although in recent years authorities have disincentivized the use of sawdust, it is still used for filling puppets. Therefore, we selected the average emission factor of 8.10 g PM 2.5 per kg of biomass, proposed for wood waste in the European emission factor database [42]. Therefore, the estimated emission due to the burning of puppets was 0.37 t PM 2.5 .
The total estimated emission due to the burning of fireworks and puppets was 1.73 t PM 2.5 , corresponding to an average emission factor of 2.7 g PM 2.5 per inhabitant. We distributed the total emission over 1 January 2022 at 00:00 LT, 1 January 2022 at 01:00 LT, and 1 January 2022 at 02:00 LT, assuming percentages of 40%, 32%, and 28%, respectively. Then, the hourly emissions were spatially distributed in the modeling grid (Section 2.3), based on the population density map and index of unsatisfied needs [33]. The black component of Figure 2a indicates the total PM 2.5 hourly emissions associated with the use of fireworks and puppets burning, starting on 1 January at 00:00 LT. Figure 2b indicates the PM 2.5 gridded emissions on 1 January 2022 at 00:00 LT. Figure 2a also indicates the hourly emissions from on-road traffic, industries, and the power facility, as background PM 2.5 sources for modeling. Figure 2b,c depict the gridded emissions on 1 January 2022 at 00:00 LT from background sources and from New Year's burning activities, respectively.

Approach for Modeling the Dispersion of PM 2.5
We used the 3D Eulerian Weather Research and Forecasting with Chemistry (WRF-Chem 3.2) [43] to model meteorology and PM 2.5 dispersion in Cuenca on 31 December 2021 and 1 January 2022. WRF-Chem is a last-generation non-hydrostatic model used for numerical modeling and solving equations of atmospheric motion from global to local scales. Meteorological simulations were performed through a master domain of 70 × 70 cells (27 km each) and 3 nested subdomains. The cells of the third domain (100 × 82) have 1 km of resolution and cover the region of Cuenca (Figure 1c). For the third subdomain, hourly PM 2.5 emissions were applied to WRF-Chem. As background sources, we used the PM 2.5 emissions from on-road traffic, industries, and the power facility in the northeast of the city corresponding to a festival day, the hourly emissions of which are indicated in Figure 2a. The PM 2.5 emissions on 1 January 2022 totaled 2.8 t d −1 , corresponding with the 1.7 t d −1 (Section 2.2), 0.6 t d −1 , 0.3 t d −1 , and 0.1 t d −1 of the New Year's combustion activities, on-road traffic, power facility, and industries, respectively. Background emissions were obtained from the last emission inventory from Cuenca. Initial and boundary conditions were generated from the final NCEP FNL operational global analysis data [44]. Table 1 summarizes the options selected for modeling.
First, we assessed the performance for modeling the hourly temperature and wind speeds at the surface through the following metrics: GE: gross error; RMSE: root-mean-square error; N: number of values; Pm: mean modeled value; Om: mean observed value; Pi: modeled value; Oi: observed value. Table 2 indicates the benchmark values for these indicators. According to the WHO guidelines, we obtained the maximum mean PM 2.5 concentrations during 24 consecutive hours from 31 December 2021 to 1 January 2022, using the records provided by each of the six stations currently measuring PM 2.5 ( Figure 1d). Then, the modeling values at the location of each station were compared with the corresponding records. We considered that the model captured the observed records if the difference was less than 50%.   Figure 5 compares the hourly PM 2.5 records with the corresponding modeling levels at the six stations. In general, the modeled values captured the trend of records. Figure 6 compares the PM 2.5 (maximum 24 h mean) records and the corresponding computed values at each station. Differences between records and modeled values ranged between 0.3% and 13.9%. The linear fit reached a value of 0.89 for the coefficient of determination, highlighting the strength of the relationship between records and modeled values [54]. These metrics indicated that the simulation of PM 2.5 dispersion was properly performed.
Although it was a component with a high level of uncertainty, the results indicated that the total PM 2.5 emissions were acceptable, both in magnitude and in their temporal and spatial distribution. Therefore, it is feasible that these emissions were used for the numerical experiment to simulate the effects of PM 2.5 levels, in case the emissions started before midnight on 31 December 2021, and to estimate the changes in PM 2.5 concentrations for the meteorological conditions corresponding to 31 December and 1 January of the previous years.       Applying the decrease of 0.65% in mortality per 10 µg m −3 decrease in PM2.5 concentrations [53,55], the emissions on 31 December at 21:00 implied reductions between 0.17%       Applying the decrease of 0.65% in mortality per 10 µg m −3 decrease in PM 2.5 concentrations [53,55], the emissions on 31 December at 21:00 implied reductions between 0.17% and 2.15% in all non-accidental mortality. The emissions on 31 December at 18:00 resulted in decreases between 0.55% and 3.25%.

Modeling PM2.5 Dispersion on 31 December and 1 January
The results of 31 December 2020 (2020-2021) with emissions at 21:00 show an average decrease of 2.6 µg m −3 . At the EIE station, the model indicated an average decrease of 14.0 µg m −3 , although an increase of 11.6 µg m −3 was obtained at the CCA station (Table 3). These results deserve specific analysis. Figure 8  For the reference scenario, on 1 January 2021 at 00:00 LT, wind was mainly generated from the northeast, transporting the pollutants to the southwest. In addition, in most of the urban area, PBL heights between 100 and 500 m were modeled, with values greater than 200 m in the historic center, an area in which concentrations of up to 300 µg m −3 were obtained. In the northern part of the urban area, concentrations between 10 and 100 µg m −3 were modeled.
For the scenario with New Year's emissions starting at 21:00 LT, the wind moved from the south in the urban area, transporting emissions to the north. The PBL presented values greater than 200 m in a larger area than the values on 1 January 2021 at 00:00 LT. However, concentrations of up to 200 µg m −3 were obtained in the northern part of the urban area. This increase was modeled in the coverage areas of the CCA and CEB stations. The numerical experiment indicates that wind transports the emissions to the north of the historic center. Although the greater heights of the PBL at the CCA and CEB stations allow for more efficient dispersion, the transport of PM 2.5 from the historic center has a more significant effect and, consequently, higher concentrations than the reference scenario. This particular result highlights the complexity of the interaction between emissions and meteorological conditions, and indicates that wind direction can promote local increases in PM 2.5 levels, despite the greater height of the PBL. Although the model's performance in capturing the temperature and wind speed for 31 December 2021 and 1 January 2022 was evaluated in this study, this result highlights the need to include the evaluation of the modeled wind direction in the future. It also highlights the importance of carrying out specific evaluations for other regions or cities due to their topographic conditions, geographic location, and the temporal/spatial configuration of their emissions. In future, the competent authority could start with the generation and delivery of the forecasted PM 2.5 concentrations, considering that emissions will occur on 1 January at 00:00 LT. This information will act as a warning and recommendation for preventive measures to reduce exposure, especially for vulnerable populations or people with pre-existing diseases. Suppose the forecast indicates atmospheric conditions with high stability and, therefore, high PM 2.5 concentrations. In that case, the competent authority could request that citizens avoid or reduce emissions arising from the use of fireworks, firecrackers, and burning puppets and limit the time they spend outdoors.
Based on forecasted PM 2.5 levels, the competent authority could also inform citizens about the benefits of advancing emission activities before midnight on 31 December. In this sense, this information can promote the collaboration of at least part of the population. Implementing a forecasting system based on the approach described in this study is challenging. For this purpose, it is necessary to have trained personnel and sufficient computational resources to have the modeled results sufficiently in advance.
In future, the performance of numerical models to capture unusual conditions, such as situations of strong atmospheric stability, should be investigated. The influence of the initial conditions also deserves dedicated research, because the performance of numerical models for meteorological forecasting tends to gradually degrade after a specific time scale (time after the initial conditions) [21]. Weather forecasting at small scales, as in the case of Cuenca, deteriorates more rapidly than at larger scales. In addition, the first hours of a weather forecast are relatively useless while the numerical model adjusts to imbalances in the initial conditions.

Limitations of This Study and Future Activities
Our estimation of the New Year's emissions (1.73 t PM 2.5 ) represents approximately 0.2% of the total PM 2.5 emissions (907 t per year) reported in the emission inventory of Cuenca for 2014 [34]. This value is consistent with the criteria of the Environment Agency of Austria [56], an entity which indicated that on New Year's Eve, fireworks account for less than 0.6% of the annual PM 2.5 emissions. However, we highlight the uncertainty of the emissions, and the limitation or lack of information in Ecuador to characterize the corresponding emission activities through a bottom-up approach. Likewise, the use of emission factors taken from the literature contribute to uncertainty, as values could differ for the altitude of Cuenca (2500 m.a.s.l.), in which the availability of oxygen in the atmosphere is 25% lower compared to sea levels [57].
We assumed that the emissions generated by the use of fireworks and the burning of puppets occurred during the first three hours of 1 January 2022. However, it was observed that the use of fireworks mainly occurs at 00:00 LT, and although the burning of puppets also mostly begins at 00:00, its emissions can last for several hours. The estimation of PM 2.5 emissions must be improved by collecting information that allows for the better characterization of activity levels and emission factors. This is a challenging task because emissions can differ year to year, both in magnitude and even location. Particulate material composition must also be characterized in terms of heavy metals, black, and elemental carbon, and monitoring the levels of particulate matter in the number of particles per unit volume. In addition, it is necessary to study the influence of relevant PM 2.5 sources, as industries with permanent and significant emissions are simultaneous with emissions associated with the arrival of a new year.
The Andean Ecuadorian region has poor coverage of radiosonde observations. This limitation implies an incomplete description of the atmospheric conditions for meteorological and air quality parameters. Although in Cuenca atmospheric monitoring has improved, the air quality network was predominantly designed to satisfy national regulation requirements, covering only surface measurements. Vertical observations need to be promoted to improve our understanding of atmospheric conditions and their interaction with air pollutants. This information will allow for more complete studies about modeling performance and improved results for air quality management.
We explored a novel and promising approach for controlling the effects of PM 2.5 emissions due to New Year's burning activities, in an Andean city with complex topography. Other ways of reducing PM 2.5 emissions, such as limiting or controlling the sale of fireworks and the burning of puppets, should be assessed. These options have been analyzed in other regions (e.g., [31,58]). However, we did not find studies based on the approach used in this study.
Animated .gif files of the modeled PM 2.5 dispersion for emissions on 1 January at 00:00 LT ( Figure S1), and for the alternatives on 31 December 2021 at 21:00 LT ( Figure S2) and 31 December 2021 at 18:00 LT ( Figure S3), are available in the Supplementary Materials. These files were created using the integrated data viewer tool [59].

Conclusions and Summary
We modeled the dispersion of PM 2.5 emitted into the atmosphere, due to the use of fireworks and other combustion activities associated with the arrival of a new year in an Andean Ecuadorian city. Our numerical experiments indicated that one way to reduce the effects on air quality is advancing the time of emissions to hours before midnight on 31 December. Our modeling approach applied to the last six years and provided mean decreases between 2.3% and 40.7% in the 24 h mean PM 2.5 concentrations, if the emissions begin on 31 December at 21:00 LT. The decrease ranged between 21.4% and 61.0% if emissions began on 31 December at 18:00 LT. The reduction in concentrations is produced by better atmospheric conditions to disperse PM 2.5 . Higher benefits can be obtained if the emissions take place at times of higher planetary boundary layer depths, when atmospheric conditions promote the pollutant's dispersion. This approach can be applied to assess the benefits of advancing New Year's emissions in other regions, based on their topographic conditions, geographic location, and the temporal/spatial configuration of their emissions.
In the future, the performance of numerical models to capture unusual conditions, such as situations of strong atmospheric stability, should be investigated. Moreover, the influence of the initial conditions deserves a dedicated study. The particulate material composition must be characterized in terms of heavy metals, black and elemental carbon contents. Monitoring should cover the levels of particulate matter in the number of particles per unit volume. In addition, it is necessary to have a dedicated study investigating the influence of relevant PM 2.5 sources in industries with permanent and significant emissions, with the emissions associated with the arrival of a new year. Vertical observations need to be promoted in the Andean Ecuadorian region, to improve understanding of atmospheric conditions and their interaction with air pollutants. The estimation of New Year's PM 2.5 emissions must be improved by collecting information that allows for better characterization of the activity levels and emission factors. Author Contributions: Conceptualization, R.P.; methodology, R.P. and C.S.; monitoring, C.E.; emissions, C.S. and R.P.; formal analysis and modeling, R.P.; data curation, C.S., C.E. and R.P.; writingreview and editing, R.P.; visualization, C.S. and R.P. All authors have read and agreed to the published version of the manuscript. Acknowledgments: This research is part of the project "Emisiones y Contaminación Atmosférica en el Ecuador 2021-2022". Simulations were done at the High Performance Computing system at the USFQ. Publication of this article was funded by the Universidad San Francisco de Quito USFQ Research Publication Fund.

Conflicts of Interest:
The authors declare no conflict of interest. Figure A1 compares the observed and modeled temperatures at the MUN station. The modeled values captured the trend of records, which indicate higher levels (approximately 22 • C) two hours after midday, and the minimum (11 • C) on 1 January 2022 at 06:00 LT. The values of the GE and BIAS metrics were 1.0 and 0.18 • C, respectively; both values were in the corresponding benchmark ranges (Table 2). Figure 4 compares the hourly records and modeled values of wind speed. The highest records (3 m s −1 ) were measured at midday. The minimum record (0.4 m s −1 ) was measured on 1 January 2022 at 01:00 LT. The RMSE and BIAS metrics were 0.8 and 0.4 m s −1 , respectively, which are in the corresponding benchmark ranges. Therefore, these meteorological parameters were captured by model-ing, implying the dynamic of PBL depths and atmospheric stability was also acceptably modeled. mately 22 °C) two hours after midday, and the minimum (11 °C) on 1 January 2022 at 06:00 LT. The values of the GE and BIAS metrics were 1.0 and 0.18 °C, respectively; both values were in the corresponding benchmark ranges ( Table 2). Figure 4 compares the hourly records and modeled values of wind speed. The highest records (3 m s −1 ) were measured at midday. The minimum record (0.4 m s −1 ) was measured on 1 January 2022 at 01:00 LT. The RMSE and BIAS metrics were 0.8 and 0.4 m s −1 , respectively, which are in the corresponding benchmark ranges. Therefore, these meteorological parameters were captured by modeling, implying the dynamic of PBL depths and atmospheric stability was also acceptably modeled.