Immission levels and identification of sulfur dioxide sources in La Oroya city, Peruvian Andes.

La Oroya is a city in the Peruvian Andes that has suffered a serious deterioration in its air quality, especially due to the high rate of sulfur dioxide (SO2) emissions, which underlines the importance of knowing its sources of contamination and variation over the years. In this sense, this study aimed to evaluate the immission levels and determine the sources of SO2 contamination in La Oroya. This analysis was performed using the hourly concentration data of SO2, and meteorological variables (wind speed and direction), which were analyzed for a period of three years (2018-2020). Graphs of time series, wind and pollutant roses, bivariate polar graphs, clustering k-means, nonparametric statistical tests, and the application of the conditional bivariate probability function were performed to analyze the data and identify the emission sources. The mean concentration of SO2 was 264.2 μg m-3 for the study period, where 55.66 and 2.37% of the evaluated days exceeded the guideline values recommended by the World Health Organization and the Peruvian Environmental Quality Standard for air for 24 h, respectively. The results showed a defined pattern for the daily and monthly variations, with peaks in the morning hours (0900-1000 h LT) and at the end of the year (December), respectively. The main sources of SO2 emissions identified were light and heavy vehicles that travel through the Central Highway, the La Oroya Metallurgical Complex, the transit of vehicles within the city, and the diesel-electric locomotives that provide cargo transportation services and tourism passenger transportation. The article attempts to contribute to the development of adequate air quality management policies.

Introduction

Exposure to air pollution represents a potential risk to human health, which has caused 300,000 deaths a year only in the Americas, an area where nine out of ten people breathe polluted air (CEPAL, 2020). Air pollution originates from natural and anthropogenic sources, the latter being the largest emitter of atmospheric pollutants (Lang et al., 2016; Yang et al., 2016).

Sulfur dioxide (SO2) is a criteria pollutant that, when emitted in large quantities into the atmosphere, not only pollutes water and soil through the generation of acid rain (Sudalma-Purwanto et al., 2015; Mosquera et al., 2018; Sagan et al., 2018; García et al., 2013) but also affects human health, harming the respiratory system and altering lung function (Kim et al., 2018; Yang et al., 2016).

Different studies have shown that short periods of exposure to SO2 (5–10 min) decrease lung function and increase respiratory symptoms in adults with asthma and children, the latter being the most vulnerable, due to their higher inhalation rate per unit of body mass (WHO, 2021; Pal, 2018). The study of SO2 is of great importance because it is a precursor to other air pollutants, such as PM2.5, where the SO4−2 is part of its composition (García et al., 2013; Lang et al., 2016).

Among the anthropogenic sources of SO2 are mineral smelters, combustion processes in the vehicle fleet, electric power generation plants, and industrial processes (Bari et al., 2020; Carn et al., 2007; Mosquera et al., 2018; Orellano et al., 2021; Sagan et al., 2018). Muñoz (2020) and Carn et al. (2007) reported that, in the processes of smelters and refineries, significant amounts of SO2 are released; an example is the extraction of copper from chalcopyrite (CuFeS2), which in its processing can release up to 2 tons of SO2 per ton of copper produced. Also, González et al. (2017) found that motorized vehicles in the city of Manizales in Colombia that use gasoline, diesel, and natural gas emit 26 t yr−1 of SOx, while industrial activity (metallurgy, manufacturing, and food production) generates 113.5 t yr−1 of SOx; as a consequence, these are considered important sources of SO2 emissions.

La Oroya is considered one of the cities with the highest level of pollution worldwide, caused mainly by high concentrations of SO2 (Arellano, 2019; Muñoz, 2020; Sagan et al., 2018). The city is located in central Peru, whose geomorphology is characterized by steep mountains that induce the formation of thermal inversions, which avoids the dispersion of atmospheric pollutants, favoring the increase of pollution levels in the city (Spieler, 2010). One of the economic activities that deteriorated the air quality of this city was the processing of minerals in the La Oroya Metallurgical Complex (CMLO by its acronym in Spanish), which is made up of a set of smelters and refineries, that transform polymetallic minerals into different metals and by-products, releasing considerable amounts of SO2 (Plasencia & Cabrera 2009).

SO2 has been monitored for three decades in La Oroya city. Neuman (2016), in the year 2000, found that the concentration levels of SO2 registered were almost twice the safe value established by the WHO. On the other hand, the General Directorate of Environmental Health and Food Safety of Peru (DIGESA by its acronym in Spanish) has monitored concentrations of SO2 in the city since 1999, obtaining that from 2009 to 2012 (years that coincide with the cessation of operations from the metallurgical complex), concentrations were below the Environmental Quality Standard for Air (ECA by its acronym in Spanish) of Peru (DIGESA, 2013). It should be noted that in 2012 the CMLO partially reactivated its operations; despite the increase in SO2 concentration levels, it remained below the ECA. Similarly, the Peruvian Environmental Assessment and Enforcement Agency (OEFA by its acronym in Spanish) monitors SO2 levels in La Oroya city; thus, in 2017 concentrations of SO2 were reported to exceed the ECA value (250 μg m−3) for 24 h (OEFA, 2017). The aforementioned institutions also identified some emission sources of SO2, such as the traffic of heavy and light vehicles, and sporadic emissions from the combustion of stoves and ovens of some establishments’ chimneys (DIGESA, 2013; OEFA, 2017).

The air pollutants concentrations are strongly influenced by transport, chemical transformations, and dry–wet deposition (Bartnicki et al., 2018; De Simone et al., 2014). To evaluate these variations, air quality models are used, which are carried out to estimate the distribution of pollutants released into the atmosphere (Demirarslan & Zeybek, 2021). Numerous approaches for mathematical modeling of air quality can be found in the scientific literature, one of these being the receptor model. This model uses statistical techniques to estimate the amount of emissions at local and long distances. Receptor models are often used to identify, characterize, and monitor the distribution of sources of air pollutants (Demirarslan & Zeybek, 2021; Hopke, 2016; Uria-Tellaetxe & Carslaw, 2014).

Among these models, we have the bivariate polar graphs, the k-means clustering, and the conditional bivariate probability function (Demirarslan & Zeybek, 2021), which were implemented in the openair package of the R statistical software, designed for atmospheric data analysis (Carslaw, 2021; Uria-Tellaetxe & Carslaw, 2014). Likewise, the polar annulus analysis allows the evaluation of the variation of the contributions from a specific source as a function of the wind direction and time (Yeganeh et al., 2021). The use of distribution graphs with two or more variables is highly reliable for determining and interpreting the sources of air pollution in a particular area (Hui et al., 2021; Jain et al., 2020; Wei et al., 2020). Recent studies have used the receptor models mentioned above to identify sources of air pollution (Hama et al., 2020; Jorquera & Villalobos, 2020; Prabhu et al., 2020). Hama et al. (2020) used bivariate polar plots and k-means clustering in the Delhi region, India, to identify local sources of PM2.5 and PM10, finding that road transport and domestic heating were the main sources of pollution. Jorquera and Villalobos (2020) employed k-means clustering to identify the main sources of PM2.5 and PM10 pollution in three urban areas of Chile, with particular meteorological characteristics and diverse sources of pollution. The study found vehicular traffic, residential wood burning, and windblown dust from the nearby desert environment as the main sources of these pollutants. Prabhu et al. (2020) used the conditional bivariate probability function in the Doon Valley, in the foothills of the Himalayas, to identify local sources of black carbon; the study found that local biomass burning activities for heating and cooking are the main sources of black carbon emissions. Furthermore, Yeganeh et al. (2021) employed polar annulus analysis to visualize the temporal aspects of black carbon concentrations at a spatial scale by wind direction. Due to the aforementioned, these graphs are very useful because they allow the observation of the behavior of pollutants based on wind direction and speed (Ali-Taleshi et al., 2021; Nguyen et al., 2022). Similarly, a comparison of pollutant time series data is essential because it reveals information on contributions from local sources and meteorological conditions (Eunhwa-Woogon et al., 2016).

Therefore, due to the problems that La Oroya has been facing concerning SO2 (Caycho, 2018; Muñoz, 2020) and the absence of research to date that has identified sources of SO2 in cities with a complex topography and variable meteorological conditions (Demirarslan & Zeybek, 2021; Hama et al., 2020), it is pertinent to develop air quality management plans. For this reason, the present study aims to identify emission sources and evaluate SO2 immission levels in La Oroya city through receptor models (bivariate polar plots, the k-means clustering, and the conditional bivariate probability function) during the 2018–2020 period. This research could be a useful tool in air quality management since it will provide valuable information for the formulation of strategies in favor of the local air quality improvement.

Materials and methods

Study area

La Oroya city is located in the Mantaro river basin, Yauli Province, Junín Department, Peru, between geographical coordinates 11° 30' and 12° 00' south latitude and between 75° 30' and 76° 00' west longitude, occupying a total area of 388.42 km2 of the Peruvian territory (Caycho, 2018; INGEMMET, 2021). It has an approximate population of 13,595 inhabitants (INEI, 2021), and presents a varied rugged topography with slopes between 2 and 6% in the upper area and 1 to 10% in the lower area (Ureta, 2013). The city is characterized by having a cold and dry climate, with rainfall that can reach up to 800 mm, and is surrounded by mountainous regions, a typical morphological characteristic of the Peruvian Andes (CONAM, 2006). Owing to its geographical and strategic position, it is a forced path for land transport by road and rail, due to its location along both sides of the Central Highway of Peru (Plasencia & Cabrera 2009). Figure 1 presents the study area map and the air quality monitoring station location.

Fig. 1

Map of the study area and location of the monitoring station

SO2 concentration data and meteorological parameters

The hourly SO2 concentration data for the 2018–2020 period was collected from the only air quality monitoring station managed by the OEFA, which is part of the National Air Quality Surveillance Network (OEFA, 2016; OEFA, 2018a). The station is located on the roof of the House of Culture of the Provincial Municipality of Yauli (OEFA, 2021). This station has a Thermo model 43i automatic analyzer and a Campbell model CR1000 automatic meteorological station. The first one uses ultraviolet fluorescence method analysis for the quantification of SO2 concentrations (MINAM, 2019a), while the second measures meteorological variables, such as wind speed (WS) and wind direction (WD) (OEFA, 2016; OEFA, 2018a).

The analysis method to determine the pollutant concentration, the criteria for the establishment of the monitoring station, and the reliability of the recorded data were carried out following the guidelines defined in the National Air Quality Monitoring Protocols of the years 2005 and 2019 (DIGESA, 2005a; MINAM, 2019a). Figure 1 shows the location of the monitoring station in the study area, and Table 1 presents its general description.

Table 1

General description of the air quality monitoring station

Code	Location	Latitude	Longitude
CA-CC-01	Comandante Zarate street, block N ° 1—La Oroya, on the roof of the house of culture of the Yauli Provincial Municipality. approximately 700 m from the La Oroya metallurgical complex	11°31′12.44"S	75°54′1.18"O

The station registers hourly SO2 concentration data continuously; then, these data are validated and published by the same institution (OEFA, 2021). In that sense, due to the easy access, continuous recording, robustness, and reliability of the data, the present research used these station’s data from the 2018–2020 period, to accomplish the stated objectives.

Although using a single monitoring station may limit the air quality monitoring of the entire study area (Jiang et al., 2021), it is necessary to mention that La Oroya has only been studied for its deterioration of air quality caused by the CMLO. Due to this, it was decided to implement a continuous air quality monitoring station following the guidelines of the European Union Commission Directive 2015/1480, which defines the number of monitoring stations based on the number of inhabitants within a city (EU, 2015; OEFA, 2016). However, evaluating the air quality of La Oroya from other perspectives could open the possibility of expanding the monitoring network in the case that more potential sources of pollution are found, which is the main objective of this study.

Methods

Time series

The temporal variation of a pollutant allows to identify the time when the minimum or maximum concentrations occur. Likewise, it shows the behavior pattern, which helps identify the days or months when the concentrations tend to increase or decrease and infer possible sources of contamination (Agustine et al., 2017; Althuwaynee et al., 2020).

In order to study the temporal variation, four-timescale graphs (hour-day, hourly, weekly, and monthly) were made. To compare the concentrations of SO2 with both the Environmental Quality Standard for Air (ECA) and the guideline value from the World Health Organization (WHO) a time series graph with daily averages was plotted, using the timeVariation and timePlot functions from the R software version 4.0.3 openair package, respectively (R Core Team, 2021).

Wind rose and pollutant rose

A wind rose shows the frequency distribution of the wind direction and speed in different intervals as well as its predominance and maximum value, respectively, within a time interval (Agustine et al., 2017). On the contrary, the pollutant rose applies the same logic but substitutes the wind speed for another set of data, for instance the concentration of a pollutant (Carslaw, 2021; Puentes, 2019). For the elaboration of the wind rose and pollutants rose, hourly records of wind direction, wind speed, and hourly concentrations of SO2 were used. The windRose and pollutionRose functions of the openair package of R software version 4.0.3 (R Core Team, 2021), were used to obtain the graphs mentioned above, respectively.

Bivariate polar plot

Bivariate polar plots are useful tools for the identification of potential emission sources of an atmospheric pollutant, through a graphical analysis (Carslaw & Ropkins, 2012; Habeebullah, 2013; Rojano et al., 2018). These plots are presented in polar coordinates, which show the variation of a pollutant concentration together with the speed and direction of the wind, allowing the identification of possible emission sources and their dispersion directions, as well as the dependence of the pollutant concentration on the wind speed (Khan et al., 2016). To identify possible SO2 emission sources in La Oroya city, bivariate polar plots were developed using hourly records of wind direction, wind speed, and hourly SO2 concentrations. The bivariate polar plots were obtained using the polarPlot function of the openair R software package (Carslaw, 2021). To facilitate the diagram interpretation, this function uses the generalized additive model (gam) for smoothing, which is carried out using the mgcv package (Carlaw, 2021; Carslaw & Beevers, 2013).

K-means clustering of bivariate polar plot

Although bivariate polar plots provide helpful information based only on a graphical analysis determined by the direction and wind speed ranges, they have several limitations in subsequent analyses, since they do not allow the characterization of the different ranges of pollutant concentrations (Carslaw & Beevers, 2013), limitations that the k-means grouping can overcome (Carslaw, 2021). K-means clustering is a method in which similar characteristics of the bivariate polar plot can be identified, selected, and grouped, to improve the understanding of the characteristics of potential sources (Althuwaynee et al., 2020; Carslaw & Beevers, 2013). This grouping was carried out using the polar Cluster function of the openair R package. It is a function that applies grouping techniques to the surfaces generated by the polarPlot; thus, it helps to identify potential characteristics for a more detailed analysis of SO2 sources (Carslaw, 2021).

Nonparametric statistical tests

The nonparametric statistical tests Kruskal–Wallis and U of Mann–Whitney were applied, with a significance level of 5%, to identify significant differences between the clusters formed by the k-means analysis. The first test reveals the significant differences between the defined groups by comparing the medians of the SO2 concentration (Chávez, 2020). The second test was used post hoc to compare the pairing of each group, where a Bonferroni adjustment was necessary to increase the level of significance required (Boso et al., 2019; Dinno, 2015). The analyses of both tests were carried out with the kruskal.test and pairwise.wilcox.test functions of the statistical software R version 4.0.3, respectively (R Core Team, 2021).

Conditional bivariate probability function

The conditional probability function (CPF) is used to identify pollution sources and estimates the probability come from a defined wind direction (Δθ) (Begum et al., 2010). In addition, conventional CPF estimates the probability of a measured concentration exceeding a fixed threshold. The CPF is mathematically defined as (Uria-Tellaetxe & Carslaw, 2014):where m∆ is the number of samples in the wind sector Δθ whose concentration C is greater than or equal to the reference value x, and n∆ is the total number of samples in the wind sector Δθ.

On the other hand, the conditional bivariate probability function (CBPF) couples the CPF with the wind speed, being this the third parameter, subject to be replaced by another variable that allows the distinction between the types of sources based on dispersion (Carslaw & Beevers, 2013). This function designates the analyzed contaminant concentration to the cells of direction and wind speed ranges, providing greater reliability than conventional methods that only establish the concentrations of pollutants in the sectors determined by the wind direction (Carslaw, 2021). The incorporation of a third variable provides more information regarding the type of source of interest. In that sense, the CBPF is appropriate for the identification of pollution sources in areas characterized by a high complexity of emission sources (Althuwaynee et al., 2020; Uria-Tellaetxe & Carslaw, 2014). The CBPF can be defined as (Carslaw, 2021):where m∆ is the number of samples in the wind sector Δθ with a wind speed interval Δu, whose concentration C is greater than a reference value x, and n∆ is the total number of samples in that wind speed–direction interval. The CBPF was used to validate the most significant sources of SO2. The implementation of this technique was carried out using the polarPlot function with the “cpf” statistic from the openair package of the R software version 4.0.3 (R Core Team, 2021).

Results and discussions

The descriptive statistics of SO2 in the study area are presented in Table 2. It is observed that the hourly concentrations of SO2 ranged from 6.5 to 3,315.1 μg m−3, 7.9 to 2,964.8 μg m−3 and 4.5 to 2,894.2 μg m−3, with averages of 78.74 μg m−3, 57.45 μg m−3 and 65.41 μg m−3 for the years 2018, 2019 and 2020, respectively. These variations may be associated with the influence that meteorological parameters (Althuwaynee et al., 2020) and terrain topography have over the behavior of air pollutants (Alván, 2010; Arrieta, 2016).

Table 2

Descriptive statistics of SO2 levels in La Oroya city

	2018	2019	2020
Minimum	6.5	7.9	4.5
Maximum	3,315.10	2,964.80	2,894.20
1. Quartile	11.3	10.7	10.7
3. Quartile	35.6	27.5	16.5
Mean	78.74 ± 2.27	57.45 ± 1.58	65.41 ± 2.33
Median	13.9	14.7	12.3
Stdev	207.1	144.3	211.1
Skewness	5.9	6.1	6.4
Kurtosis	48.6	55.4	52
Missing data %	4.9	5.3	6.7

On the other hand, it is observed that the asymmetry and kurtosis values varied in the range of 5.9 to 6.4 and 48.66 to 55.4, respectively. These results demonstrate the asymmetric behavior of the SO2 concentration. According to O'Leary et al. (2016), Ottosen and Kumar (2019), an asymmetric behavior of the data occurs possibly due to the influence of outliers above the 75 percentile.

Box-and-whisker plots

Due to the significant presence of outliers, the box-and-whisker diagram is presented without drawing the outliers in order to achieve a better graphical representation (Fig. 2) (R core team, 2021). In this research, outliers were not eliminated, because their presence is indicative of the diversity and a high number of SO2 pollution sources in La Oroya city (Caycho, 2018; Estevan et al., 2019, 2020; Plasencia & Cabrera, 2009). Likewise, it is common that in atmospheric sciences, data sets with outliers are presented, so their exclusion, systematic elimination, or transformation could hide important atmospheric environmental processes and lead to incorrect conclusions that can difficult the decision making, or even lead to the adoption of inadequate management programs (Benhadi-Marín, 2018; O'Leary et al., 2016). Finally, it is essential to mention that the hourly data concentrations of SO2 are validated by OEFA, and the monitoring station is subject to periodic evaluations carried out by this institution (OEFA, 2021).

Fig. 2

Box-and-whisker plot of SO2 concentrations

The graph shows that the interquartile range is comparatively higher in 2018 and lower in 2020. Therefore, it is inferred that there was higher variability of SO2 concentrations in 2018 and intermediate variability in the year 2020.

The temporal variation in four timescales (hour-day, hourly, weekly, and monthly) of SO2 concentrations during 2018, 2019, and 2020 (Carslaw, 2021) is shown in Fig. 3. It is observed that during the day there is a peak between 0900 and 1000 h local time (LT) with an average concentration of 264.2 μg m-3 that exceeds the guideline values, of 40 μg m−3 and 250 μg m−3, established by WHO and the national ECA, respectively. At a monthly level, SO2 concentrations do not follow a homogeneous trend in the three years period, but they do register a minimum average value in June (26.4 μg m−3) and a maximum average in December (126.1 μg m−3). At the weekly level, the variation was more notable, without defined patterns.

Fig. 3

Temporal variation of SO2 concentrations

Sagan et al. (2018) and Zarauz et al. (2010) also recorded variations in SO2 concentrations in La Oroya and found that the main factors influencing the dispersion of pollutants were the oscillation of the planetary boundary layer, topography, thermal inversions, and changes in wind direction and speed. This is due to the fact that the mountainous areas of the Peruvian Andes form a natural physical barrier for wind circulation and make it impossible for the polluted air to be removed from the area. In addition, the dry climate in these places increases the frequency of thermal inversions (García et al., 2013). The complex topography of La Oroya city determines the speed and direction of airflows, hampering the circulation during the contaminant transport process, thus preventing dispersion (Cheng et al., 2022; Caycho, 2018; Ibañez et al., 2021; Sagan et al., 2018; Arrieta, 2016). This factor also influences the amount of solar radiation received by the earth's surface throughout the day, which generates atmospheric pressure differences that produce air mass movements (Cheng et al., 2013, 2022; Zardi & Whiteman, 2013), and give rise to mountain–valley wind systems, which are characteristic of the city (Mamani, 2013; Vidal & Pérez, 2017).

Bei et al. (2018) found that the atmospheric contaminants formed in the afternoon are transported by the valley winds from an urban area to the mountains, being stored in the latter, while at night, due to the action of mountain winds, they are put back to the urban area, causing greater pollution. Cao et al. (2020) reported that the mountain–valley winds were stronger on sunny days than on cloudy days in the city of Sichuan, China; they also mentioned that these systems originate the transport of pollutants between the city and the mountains, increasing the level of air pollution in the city. The phenomenon described above has also been evidenced in the city of La Oroya due to its topographical characteristics, as reported by Mamani (2013), who states that the valley–mountain wind system of the city influences the wind flows and therefore the dispersion of SO2 concentrations.

On the other hand, due to its geographical, meteorological, and topographical characteristics, the city of La Oroya is exposed to temperature inversions that prevent adequate diffusion of gaseous pollutants, settling in the city for long periods of time (Alván, 2010; Cerdestav & Barandiarán, 2002; Mendiola et al., 2017; Tello, 2014). The planetary boundary layer of the city is between 700 and 800 m; however, this vertical structure varies during the day and night, and between wet and dry periods.

Below the boundary layer, the thermal inversion phenomenon occurs, which originates and varies according to the aforementioned periods, presenting a greater negative impact on air quality during the dry period. During this period, the thermal inversion generally becomes frequent at the end of the day between 1800 h LT and 1900 h LT, until 0900 h LT and 1000 h LT the next day, so there are high levels of gas concentration in the mornings, while in the afternoon the intensity of the winds increases and the atmosphere becomes unstable, which favors the dispersion of gases and reduces the impact on air quality (CONAM, 2004; CONAM, 2006; Córdova, 2008). According to Giovannini et al. (2020), this phenomenon in mountainous regions can even extend to periods ranging from several days to several weeks, since the low visibility factor of the sky limits the penetration of solar radiation in the valleys, and the night drainage winds favor the accumulation of cold air at low levels. The temperature decrease during the mornings generates thermal inversion and, in combination with the low humidity of the environment, causes contaminants to be trapped in the bottom and prolong its duration (García et al., 2013).

The hourly peak in the mornings is possibly related to the greater intensity of vehicular traffic on the Central Highway, where heavy and light vehicles circulate, and to which the emissions from fugitive sources from the CMLO would be added (Alarcón, 2016; Integral Consulting, 2005; Red Muqui, 2020; OEFA, 2018b,c–d). In Peru, it has been reported that 90 and 58.1% of cargo and passenger transport vehicles, respectively, use petroleum (diesel), a fossil fuel with high concentrations of sulfur and an important source of SO2 (MINEM, 2017; Fei et al., 2018). The patterns identified between vehicle fleet and SO2 concentrations are explained because SO2 is generated as a product of the incomplete combustion of sulfurized hydrocarbon fuels (USEPA, 2020).

It is important to mention that the proximity of the monitoring station to the main highway may have favored the recording of peak concentrations from heavy traffic at the hourly level (Mohtar et al., 2018). In their study, De Souza et al. (2019) reported a maximum peak of high concentrations of SO2 at 1100 h LT and related it to the increase of vehicular traffic in Campo Grande, Brazil, in 2016. The influence of transport on SO2 concentration was also reported by Yang et al. (2018), as they determined that transportation contributed to 10% of China's total emissions during 2002.

In addition, the industrial sector has been recorded as a large contributor to SO2 emissions, according to various studies. For example, Min et al. (2021) found that industrial processes and industrial combustion in Sichuan, China were the main contributors of SO2 emissions with 46% and 29%, respectively, in 2017. Akyuz and Kaynak (2019) reported that in Turkey, the energy sector was the main source of SO2 emission because 70% of the total electricity was generated with a 50% share of coal-fired power plants (CPP) using sub-bituminous and lignite coal as fuel, which contain a large amount of sulfur. Li et al. (2017) calculated China’s emissions inventory for the year 2010 by sectors and found that the industrial sector contributed 57.39% to total SO2 emissions.

However, in Peru, according to the Lima-Callao Air Quality Management Diagnosis of 2019, mobile sources were the main contributors of SO2 emissions with a contribution of 54.79%, while point and area sources only contributed 33.02 and 12.18%, respectively. Within the point sources, the industrial activity of zinc refining was the main contributor to SO2 emissions (MINAM, 2019b). On the other hand, although in La Oroya, the main contribution of SO2 emissions (380,075 t-SO2 yr−1) in 2005 came from the CMLO's metal smelting and refining processes (representing 99% of total SO2 emissions) (CONAM, 2004; CONAM, 2006; DIGESA, 2005b), fugitive sources from it have also been responsible for most of the local impacts, especially in La Oroya Antigua and La Oroya Nueva (DIGESA, 2005b; Integral Consulting, 2005).

Finally, as described in the previous paragraph, in Lima-Callao and La Oroya, SO2 contributions from metal refining activities have been considerable, but their contribution has been decreasing. In the case of La Oroya, these emissions have been decreasing with the partial cessation of the CMLO's activities. Therefore, currently other sources, mainly mobile sources, could equal or exceed the contribution of the polymetallic smelter and refining.

The reduction of SO2 emissions from the CMLO affects the immission levels of La Oroya. Morales (2018) evaluated the concentrations of SO2 during the operation and shutdown of the CMLO, and found that the concentrations of SO2 were higher than 2,500 μg m−3 in the year 2008, while, after the operations stopped, the SO2 immission level was below the minimum value established by DIGESA to consider La Oroya city in a state of alert (arithmetic average of 24 h > 500 μg m−3 of SO2).

Comparison with ECA and WHO guideline values.

The temporal variation of SO2 in La Oroya with a daily resolution (24-h average) is presented in Fig. 4. The daily levels exceeded the WHO guideline value (40 μg m−3) on 610 days, representing 55.66% of the days of the study period (2018–2020). The national ECA for SO2 (250 μg m−3) was exceeded during 26 days, representing 2.37% of days of the study period. The days with the highest concentrations that exceeded the national regulations were recorded at the end of 2020, which was possibly associated with the partial lifting of health restrictions in the framework of the COVID-19 pandemic (Rojas et al., 2021; Toro et al., 2021; Morales et al., 2021). Several studies worldwide identified an increase in air pollutant concentrations as the pandemic restrictive measures were lifted (Anil & Alagha, 2021; He et al., 2021; Wang et al., 2021). According to the WHO, high levels of exposure to this gas can affect the respiratory system, lung function, cause eye irritation and aggravate asthma and chronic bronchitis; it also increases hospital admissions for heart disease and mortality rates (WHO, 2006; WHO, 2021).

Fig. 4

Daily temporal variation of SO2 in La Oroya city

Source identification of SO2 by wind rose and pollutant rose

The wind rose for the study period is shown in Fig. 5. The predominant direction of the wind is from the west (W) throughout the four seasons. Likewise, variable directions are evident in all seasons. Zarauz (2011) observed similar wind direction behavior at the CMLO main chimney, located approximately 700 m SE from the monitoring station. The mean wind speed ranged from 1.177 m s−1 (autumn) to 1.4197 m s−1 (spring). The lower wind speeds recorded in March, April, and May, added to the presence of thermal inversions, cause stagnant air masses and a consequent increase of SO2 concentrations in La Oroya (Espinoza & Alderete, 2021; Sagan et al., 2018; Silva et al., 2017).

Fig. 5

Wind rose by seasons of the year during the 2018–2020 period

The inverse relationship between wind speed and SO2 concentration is evidenced in the pollutant rose (Fig. 6). At lower wind speeds (0 to 0.8 m s−1), there is a higher average concentration (93.732 μg m−3) of the pollutant, while at higher wind speeds (1.7 to 5.6 m s−1), there is a lower average concentration (32.166 μg m−3) of SO2. Mid-concentrations of the pollutant were recorded in the speed ranges from 0.8 to 1.1 m s−1 (72.739 μg m−3) and 1.1 to 1.7 m s−1 (64.593 μg m−3). Other studies also reported that at lower wind speeds, higher concentrations of pollutants are recorded, due to its low dispersion (Sooktawee et al., 2020; Althuwaynee et al., 2020; Perišić, 2020; Jindamanee et al., 2020; Mohtar et al., 2018; Jang et al., 2017).

Fig. 6

Pollutant roses for SO2 hourly concentrations by wind speed intervals during the 2018–2020 period

When analyzing the yearly seasonal correlation between wind speed and SO2 concentration, it was determined that during the summer, the correlation was higher (|− 0.10|≤ r ≤|− 0.48|), followed by spring (|− 0.12| ≤ r ≤|− 0.33|), autumn (|− 0.04|≤ r ≤|− 0.34|), and winter (|− 0.02|≤ r ≤|− 0.09|). These correlations improve when analyzing the data set by month. It was observed that all the correlations between both parameters are negative, which confirmed the inverse relationship between wind speed and SO2 concentrations, these being stronger in the summer and spring seasons. Dandotiya et al. (2020) and Chao et al. (2021) also found an inverse correlation (r = − 0.38 and r = − 0.34, respectively) between both parameters. Li and Xie (2016) explained that the inverse relationship between wind speeds and SO2 concentrations is because lower wind speeds do not promote the diffusion of the pollutant, which results in higher monthly concentrations of SO2, while, in conditions with higher wind speeds, the diffusion of SO2 is favored, so monthly concentration ​​are lower.

The inverse relationship found in La Oroya could also be related to the proximity between the emission source and the monitoring station, due to its direct impact on the readings of the monitoring equipment, but this would not limit the registration of lower concentrations coming from possible distant sources (Althuwaynee et al., 2021; Moufarrej et al., 2020; Sooktawee et al., 2020).

Pollution peaks within the range of 150 to 3,315.1 μg m−3 occurred more frequently at low velocities and in the south (S), southeast (SE), and east (E) directions, where one of Peru's main highways, the CMLO and the railroad with the Huancayo–La Oroya–Lima route are located. Between 4,000 and 4,300 vehicles per day travel along the Central Highway, including heavy and light vehicles, most of which use petroleum (diesel) as their main fuel (Alarcón, 2016; MINEM, 2017). The CMLO is classified as one of the largest refineries in Peru and has been declared a major source of SO2. However, in recent years it has been operating at a limited capacity and a non-constant period (Requejo, 2020; TFA, 2021; MINEM, 2020; Red Muqui, 2020). Trains with diesel–electric locomotives run on the railroad, mainly for freight transportation and, to a lesser extent, for tourist passenger service with scheduled departures along the Huancayo–La Oroya–Lima route (OSITRAN, 2020; FVCA, 2006).

Additionally, high concentrations were observed at low wind speeds coming from the west (W) of the monitoring station, which may be associated with vehicle traffic within La Oroya city (Concepción & Rodríguez, 2014). In addition, to the west–south–west (WSW) and approximately one kilometer away, there is a critical point of vehicular congestion due to the convergence of three highways from Cerro de Pasco, Huancayo, and Lima, where light and heavy vehicles circulate (CONAM, 2004). Finally, to the west (W) and southwest (SW) there is the railway that connects the cities of Cerro de Pasco, Junín, La Oroya, and Lima, where freight trains from important mining centers frequently pass through, being the main transported products: copper concentrate, zinc concentrate, lead concentrate, and industrial sulfuric acid (OSITRAN, 2020; MINCUL, 2019). The aforementioned sources explain the high concentrations recorded in the west (W) direction.

Given the presence of mobile sources in all wind directions and their significant contribution to SO2 concentrations in La Oroya, to complement the study, SO2 emissions from this source were estimated for the 2018–2020 period. The calculation of emissions followed the updated Tier 1 methodology described in the Air pollutant emission inventory guidebook of the European Environment Agency (EEA, 2021). Vehicle fleet data estimated by the Ministry of Transport and Communications (MTC by its acronym in Spanish) (MTC, 2021) were used and it was considered that the fuels used by these vehicles were Gasohol, diesel, and LPG (OSINERGMIN, 2021). It was estimated that for the years 2018, 2019, and 2020, SO2 emissions were 9,852, 10,103, and 8,704 t-SO2 yr−1, respectively. These values ​​far exceeded the emissions from mobile sources of SO2 reported in the latest emissions inventory of La Oroya (83 t-SO2 yr−1) (DIGESA, 2005b), which clearly shows the progressive growth of the vehicle fleet throughout the years (MTC, 2021), bringing with it, the increase in the concentration of SO2 in the city of La Oroya. These SO2 emissions lead to an increase in the concentration of the pollutant. In the study conducted by Yang et al. (2016) in China, it was found that the emission of 10,000 tons of SO2 increases the environmental concentration of the pollutant by 0.463 μg m−3. This estimation confirms the importance of La Oroya’s vehicle fleet in the contributions of SO2 concentrations in recent years. The high rate of emission by this source can be alarming. A deeper study and more detailed analysis of the sulfur dioxide emissions in this city or the department of Junín and its contribution to the concentrations of sulfur dioxide can be a relevant subject for future research.

On the other hand, the contamination range from 0 to 50 μg m−3 of SO2 occurs in all wind directions and wind speeds, which can be attributed to the sources close to the monitoring point (Sooktawee et al., 2020), in addition to those already mentioned. The OEFA (2018e, f, g, h, i), during the execution of its air quality surveillance plan, observed sporadic emissions from chimneys at 30 m from the monitoring point, which were generated as a result of the combustion of fuels in ovens, restaurant kitchens, and other establishments. According to the air quality baseline diagnosis reported in 2004, these emissions would contribute minimally to SO2 concentrations, since the use of oil and kerosene in some establishments, mainly bakeries and restaurants, emits low amounts of SO2 (CONAM, 2004).

Source identification of SO2 by bivariate polar plots and k-means clustering

The bivariate polar graphs of SO2 concentrations by seasons for the 2018–2020 period in La Oroya city are shown in Fig. 7. The color scale in the figure shows the SO2 concentration in μg m−3, and the radial scale represents the wind speed in m s−1 which increases radially from the center of the diagram. High concentrations of SO2 are observed mainly in the southwest (SW), south (S), and southeast (SE) directions throughout the four seasons of the year, being particularly high during the spring and summer of the southern hemisphere. The sources of pollution that might contribute to the high concentrations of SO2 in these directions could be the Central Highway, the CMLO, and the railways through which diesel freight and passenger trains travel (Requejo, 2020; TFA, 2021; MINEM, 2020; Red Muqui, 2020). In the air quality surveillance monitoring reports from the years 2016, 2017, and 2018, the presence of emissions from the CMLO, from the automotive fleet (cars, trucks, and trains), and to a lesser extent from some commercial and service establishments (restaurants and bakeries) was evidenced (OEFA 2016, OEFA, 2017, OEFA, 2018g, h–i).

Fig. 7

Bivariate polar graphs from the concentrations of SO2, WS, and WD

Likewise, a possible source of contamination in the SW direction, with less contribution than the previous ones and that is remarkable during the summer and autumn seasons in the wind speed range from 2 to 3 m s−1, would also be associated with the transit of cargo and passenger trains that follow the Huancayo–La Oroya–Lima and Unishp (Cerro de Pasco)–La Oroya–Lima routes, that circulate along the railway located at a distance of approximately 500 m from the monitoring station (OSITRAN, 2020; FVCA, 2006). Finally, the graphs also showed that in the rest of the wind directions, SO2 concentrations are minimal, which could indicate that there is no significant source of pollution in those directions.

The k-means clustering of the bivariate polar graphs for the SO2 concentrations is shown in Fig. 8. The diagram was divided into 3 clusters, to characterize the SO2 emission sources in the Oroya city. The red (CL1), blue (CL2), and green (CL3) clusters represent high, medium, and low concentrations of SO2, respectively. The clusters with the highest concentration in the SE, S and SW directions are mainly influenced by vehicular emissions from the Central Highway, produced by the oil or diesel incomplete combustion process, where sulfur as part of the composition of these fuels is burned and emitted in the form of SO2 (Khobragade et al., 2018). On the other hand, CMLO emissions might also contribute to the high concentrations in the red cluster (CL1). The CMLO mainly processes copper, lead, and zinc concentrate, which contain high concentrations of sulfur, an element that is released in high amounts as SO2 during the heating process of the concentrate (Carn et al., 2007; Della et al., 2017; Fioletov et al., 2016). Zhang et al. (2021) identified metal smelters as one of the important sources of SO2 in China, contributing to 20% of total SO2 emissions in 2012. Mohtar et al. (2018) reported that the main source of SO2 in power plants and industrial facilities is the burning of fossil fuels. Finally, the trains that travel on the railroad also contribute to SO2 emissions since they use diesel–electric locomotives (OSITRAN, 2020). The high emissions of SO2 generated during the incomplete combustion of diesel or gasoline could be reduced by lowering the sulfur content in these fuels (Miranda et al., 2016; Mohtar et al., 2018; Rodríguez, 2017).

Fig. 8

Clustering of k-means on the bivariate polar graph

Figure 8 reaffirmed that the highest concentrations of SO2 come from the SE, S, and SW directions, being associated with three possible sources of contamination such as the traffic of vehicles on the Central Highway, the CMLO, and the trains.

The multiannual temporal variation of SO2 concentrations for each defined cluster is shown in Fig. 9. Temporal analysis shows average concentrations with a 95% confidence interval (Carslaw, 2021). The cluster with the highest concentration (CL1) reflects the characteristic behavior of vehicular flow in La Oroya city (Alarcón, 2016; Concepción & Rodríguez, 2014). In general, there are two daily peaks, one in the morning and the other in the afternoon, which is related to the traffic congestion caused in those periods (Alarcón, 2016). In addition, the decrease in the concentration of SO2 is observed during noon, which may be related to the turbulent conditions of the atmosphere at this time, which enhances the dispersion of pollutants (Espinoza & Alderete, 2021; Habeebullah, 2013; Sagan et al., 2018; Zarauz et al., 2010). Likewise, it was identified that the concentrations of SO2 associated with CL1 were different from those of CL2 and CL3, but, between the latter, there were no significant differences.

Fig. 9

Multiyear temporal variation of the SO2 clusters in La Oroya city

Monthly there is a marked and differentiated behavior between the three clusters. The CL2 and CL3 show an almost constant average in all months of the year, which may indicate that the concentrations come from fixed sources, while the monthly variation of CL1 follows a variable behavior in which lower values are evidenced in June and higher in the last months of the year. The characteristic monthly behavior of CL1 could be mainly associated with the development of economic activities, particularly the production, transportation, and processing of minerals. According to records from the Ministry of Energy and Mines (MINEM, by its acronym in Spanish) and the Central Reserve Bank of Peru (BCRP, by its acronym in Spanish), the average production of metals such as Cu, Zn, and Pb in the Junín region for the study period was higher at the end of the year (September, October, November, and December) and less in the middle of the year (April, May, June, and July). The production was exceptionally high at the end of 2020, a period that coincided with the progressive lifting of sanitary restrictions due to COVID-19 (MINEM, 2021; BCRP, 2021; Rojas et al., 2021). In addition, Alarcón (2016) reported that the interprovincial transport service along the Central Highway also followed a defined pattern, registering the highest number of trips at the beginning and the end of the year and the lowest at midyear. In this sense, the processing of the minerals and the number of trips made in trains and vehicles throughout the year would explain the monthly behavior of the cluster that registered the highest concentrations of SO2.

Table 3 presents the results of the Kruskal–Wallis and Mann–Whitney U tests for the clusters of SO2 concentrations. The nonparametric Kruskal–Wallis test identified the existence of differences between the groups formed by the k-means cluster analysis at a significance level of 5% (p < 0.05), which further reinforces the findings. The Mann–Whitney U test identified the group that differs from the other groups formed, through pairwise multiple comparisons (Calazans et al., 2018). The results of the multiple comparison test confirmed that the red cluster (CL1) contributed more to the concentrations of SO2, which is associated with emissions from heavy and light vehicle traffic, with contributions from the CMLO and with emissions from train traffic. Clusters 2 and 3 also showed significant differences at 5%; however, their average values and behavior are similar, which may indicate predominant emission sources with similarity for both clusters or the absence of important sources in the directions they represent (Al-Harbi et al., 2020). The Mann–Whitney U test proves to be an effective tool for determining significant differences between clusters (Chávez, 2020).

Table 3

Results of the Kruskal–Wallis and Mann–Whitney U tests in the SO2 clusters

Cluster	Mean	Chi-squared	df	p-value	Pairwise	Mann–Whitney
CL1	238.0	4823	2	< 2.2e–16	CL1–CL2	< 2e–16
CL2	28.1				CL1–CL3	< 2e–16
CL3	30.6				CL2–CL3	< 2e–16

Source identification of SO2 by conditional bivariate probability function

Finally, the graphs generated with the CBPF for concentrations of SO2 are shown in Fig. 10. CBPF plots show relevant hidden information about other pollution sources that are not distinguishable in bivariate polar plots, so their use is appropriate for places with high emission source complexity (Althuwaynee et al., 2020, 2021). A 75th percentile was established to represent the sources with the highest contribution to SO2 concentrations (Uria-Tellaetxe & Carslaw, 2014). The analysis validated the main sources identified. Thus, significant concentrations of SO2 (> 75th percentile—25 μg m−3) were identified with a probability of occurrence of 90% or more in the SE, S, and SW directions, during the four seasons of the year and throughout the speed range. This is mainly related to the emissions from mobile sources that travel along the Central Highway and the contributions from the CMLO. A greater probability of the occurrence of significant concentrations greater than the 75th percentile was also observed in the spring and summer seasons compared to the autumn and winter seasons, which agrees with the findings from the other methods applied in this study. A higher concentration of SO2 in the spring and summer seasons, and a lower one in the autumn and winter seasons, could be related to the climatic conditions of the place under study, considering that areas with dry climates such as La Oroya, can contribute to severe haze episodes with high concentrations of SO2 (CONAM, 2006; García et al., 2013).

Fig. 10

Conditional bivariate probability function at the 75th and 90th percentiles for SO2 concentrations

The CBPF at the 90th percentile showed similar behavior to the 75th percentile. However, it helped to identify an additional significant source of SO2 in the south–southwest direction (SSW), whose predominance occurred in March, April, and May and in the spring and summer seasons. That behavior might be associated with the movement of freight and passenger trains at different times of the day, which constitute another of the mobile sources of air quality pollution in La Oroya (DIGESA, 2013).

Conclusions and recommendations

In this study, the immission levels were evaluated and the potential sources of SO2 emissions were identified from hourly records of SO2, direction, and wind speed, in La Oroya city. The mean SO2 concentration for the study period was 264.2 μg m−3. Likewise, it was found that the daily average concentrations of this pollutant exceeded the guideline value established by the WHO and the ECA air during 610 and 26 days, respectively. The highest number of days that exceeded these values was recorded in 2020, which coincided with the progressive lifting of sanitary restrictions due to COVID-19 in Peru. The daily variation of the pollutant concentrations presented a daily peak, recorded in the morning hours (0900–1000 h LT), which coincides with the characteristic behavior of the vehicular flow of the study area, added to the complex relief topographical and meteorological conditions of the place. At a monthly level, the highest monthly averages of SO2 were recorded in December; this was associated with the highest production of minerals and interprovincial passenger transport recorded in this month. The predominant winds were from the West (W) during the four seasons, with an average speed between 1.177 m s−1 (autumn) and 1.4197 m s−1 (spring), as shown in the wind roses. While the analysis with the pollutant roses allowed to determine that at lower wind speeds (0–0.8 m s−1), higher concentrations (93.732 μg m−3) of the pollutant occur in the S, SE, and E directions of the monitoring station. The bivariate polar graphs and their subsequent grouping using the k-means algorithm showed that the highest concentrations of SO2 mainly occur in the S, SE, and SW directions throughout the year, which are mainly influenced by the emissions from the fleet of light and heavy vehicles that travel daily on the Central Highway, possible emissions from the La Oroya Metallurgical Complex and the traffic of diesel–electric locomotives that travel the railways surrounding the city. The analysis using the conditional bivariate probability function with 75th and 90th percentiles made it possible to validate the sources of contamination identified with the other methods used and to estimate the concentrations most influenced by the main sources of contamination. Finally, the need arises to continue working on the development of more advanced models that incorporate other atmospheric and topographic variables, thus allowing more accurate characterization of temporal variations and the influence of these variables on the dispersion of the pollutant.