DOI: 10.20937/ATM.53482

Received: March 26, 2025; Accepted: June 10, 2025

Evaluation of air quality in Puebla, Mexico: A wavelet transform and predictive modeling approach

RESUMEN

Este artículo presenta un análisis detallado de los contaminantes del aire en la ciudad de Puebla, México, basado en datos recopilados entre 2016 y 2024. La investigación examina la variación diaria de cinco contaminantes principales: ozono (O₃), material particulado menor de 10 micras (PM₁₀), material particulado menor de 2.5 micras (PM_2.5), dióxido de azufre (SO₂) y dióxido de nitrógeno (NO₂). Para identificar tendencias significativas y patrones estacionales, se aplicaron la prueba de Mann-Kendall, el método de análisis de tendencia innovador y la transformada wavelet. Los resultados indican tendencias ascendentes estadísticamente significativas en las concentraciones de O₃, SO₂ y NO₂, mientras que los niveles de PM₁₀ y PM_2.5 han mostrado una disminución sostenida a lo largo del periodo de estudio. El análisis con escalogramas revela concentraciones de energía estacional en SO₂, posiblemente asociadas con la actividad industrial y las condiciones meteorológicas. Además, se utilizó el modelo de pronóstico Prophet para estimar los niveles de PM_2.5 y PM₁₀ entre 2022 y 2024, observándose un mejor desempeño del modelo en horizontes temporales más largos. Este estudio resulta especialmente relevante ante el crecimiento urbano y la actividad industrial en Puebla, factores que pueden contribuir al deterioro de la calidad del aire y afectar la salud de la población. La identificación de tendencias y patrones en la contaminación del aire es esencial para la implementación de estrategias de mitigación y políticas públicas orientadas a mejorar la calidad del aire en la región.

ABSTRACT

This article presents a detailed analysis of air pollutant dynamics in Puebla City, Mexico, using data collected between 2016 and 2024. The research examines the daily variation of five main pollutants: ozone (O₃), particulate matter smaller than 10 microns (PM₁₀), particulate matter smaller than 2.5 microns (PM_2.5), sulfur dioxide (SO₂), and nitrogen dioxide (NO₂). To identify significant trends and seasonal patterns, the Mann-Kendall test, innovative trend analysis (ITA), and wavelet transform were applied. The results indicate statistically significant upward trends in O₃, SO₂, and NO₂ concentrations, while PM₁₀ and PM_2.5 levels have exhibited a sustained decrease throughout the study period. The scalogram analysis highlights seasonal energy concentrations of SO₂, potentially linked to industrial activity and meteorological conditions. Additionally, the Prophet forecasting model was used to estimate PM_2.5 and PM₁₀ levels from 2022 to 2024, achieving better performance over longer time horizons. This study is particularly relevant given the urban growth and industrial activity in Puebla, factors that can contribute to the deterioration of air quality and affect the health of the population. The identification of trends and patterns in air pollution is essential for the implementation of mitigation strategies and public policies aimed at improving air quality in the region.

Keywords: Puebla City, wavelet transform, Mann-Kendall test, ITA, Prophet forecasting model.

1. Introduction

Air pollution is one of the most significant environmental issues worldwide, impacting human health, the climate, and ecosystems. According to the World Health Organization (WHO), 99% of the world’s population breathes air with pollutant levels higher than recommended. Mexico faces critical levels of air pollution, especially in its large cities. According to data from the Secretaría de Medio Ambiente y Recursos Naturales (Secretariat of Environment and Natural Resources, Semarnat) and the Red Nacional de Monitoreo Atmosférico (National Air Monitoring Network), the cities with the highest levels of air pollution are Mexico City, Monterrey, Guadalajara, Toluca, and Puebla. High levels of air pollution have led to the implementation of measures such as vehicle restriction programs, restrictions on industrial activity, and the promotion of public transport. Puebla is one of the most important and populated cities in Mexico, with growing industrial activity and a strong connection to Mexico City. Studying pollution in this area is key for several reasons: urban and industrial expansion, atmospheric interconnection with Mexico City, meteorological and orographic conditions, impact on public health, and regional mitigation strategies. These factors highlight the importance of studying air pollution in Puebla and nearby urban areas. A better understanding of pollutant behavior is essential to design effective public health interventions, develop regional mitigation strategies, and promote sustainable development across central Mexico.

According to the 2020 Population and Housing Census conducted by the Instituto Nacional de Estadística y Geografía (National Institute of Statistics and Geography, INEGI), the population of the municipality of Puebla was 1 692 181 inhabitants, making it the fourth largest city in the country after Mexico City, Monterrey, and Guadalajara. Puebla is part of the Puebla-Tlaxcala metropolitan area, whose geographic position and relatively flat topology are key for trade throughout the country.

Puebla is located at an average altitude of 2160 masl, and is surrounded by mountain ranges of the Trans-Mexican Volcanic Belt (TMVB). To the west lie the Popocatépetl and Iztaccíhuatl volcanoes, to the north La Malinche, and to the east, the Pico de Orizaba. The TMVB, where Puebla is located, is a complex geological region characterized by diverse volcanism and active crustal deformation (Gómez-Tuena et al., 2007; Suárez et al., 2019). The green area of Amozoc, located on the eastern edge of the urban area, serves as a significant green lung for the city (Fig. 1a).

Fig. 1. (a) Map of Puebla City and the surrounding areas. (b) SINAICA atmospheric monitoring stations in Puebla city. (c) Map of Mexico highlighting the proximity of Mexico City and Puebla.

The predominant climate in Puebla is temperate, with varying humidity levels and an average temperature of 16 ºC, with summer being the rainiest season. The COVID-19 lockdown temporarily reduced the concentration of carbon monoxide (CO), nitrogen dioxide (NO₂), ozone (O₃), sulfur dioxide (SO₂), and particulate matter smaller than 10 microns (PM₁₀). However, Puebla City’s air quality has been affected by a significant increase in pollutants since 2020, according to national air quality reports (INECC, 2024).

Mexican Official Standards (NOMs) that define maximum permissible limits for air pollutants incorporate recent evidence on health impacts and air quality management practices. Two types of standards are used for air quality monitoring: environmental health NOMs, which establish permissible limits for criteria pollutants, and technical NOMs, which define measurement methods for criteria pollutants.

Over recent years, extensive research has demonstrated the impact of air pollution on health. Exposure to pollutants such as particulate matter (PM) and O₃ is linked to increased death rates and hospitalizations due to respiratory and cardiovascular disease. Both short-term studies, examining daily pollution fluctuations, and long-term studies, following groups exposed to pollution over time, have found adverse effects even at very low exposure levels. Furthermore, it has been observed that these adverse effects occur even at very low levels of exposure (Brunekreef and Holgate, 2002).

This work analyzes the temporal dynamics of air pollutants in the City of Puebla. However, given the limited literature on air quality studies in this region, we compared some results with those obtained by other authors in Mexico City, a metropolitan area that has a comparable average altitude above sea level (2160 for Puebla City vs. 2240 m for Mexico City) and is located at a relatively short distance, approximately 100 km (Fig. 1c). Topographic conditions are also comparable, since both cities are surrounded by mountain ranges that influence the dispersion of pollutants (de Foy et al., 2006).

Although the meteorological and geographic conditions in both cities may be comparable, it is important to consider the variations in pollution sources and their management when analyzing their indices (Molina and Molina, 2002). Mexico City’s metropolitan area is among the largest in the world, and due to its strong industrial activity, high transportation traffic, and enormous population, it emits a lot of ozone precursors. However, Puebla City’s economy is closely tied to the automobile and textile sectors, which generate pollutants at a different rate than Mexico City.

Previous studies in Mexico City have highlighted the role of photochemical activity, complex terrain, traffic volume, and meteorology in shaping ozone dynamics (Molina and Molina, 2004). Similar investigations in Guadalajara and Monterrey report persistent ozone and particulate matter levels and identify local emission sources, topographic effects, and chemical regimes as key influencing factors (Limón-Sánchez and Smith, 2011; Hernández-Paniagua et al., 2017; Díaz-Torres et al., 2022; Fajardo-Montiel and Ramírez-Sánchez, 2023). Modeling efforts in northeastern Mexico further emphasize that the effectiveness of emission control strategies depends on the dominant precursor conditions (Sierra et al., 2013). Despite these advances, cities like Puebla, experiencing rapid urbanization and industrial growth, have received comparatively less attention in the scientific literature, underscoring the need for localized studies.

To understand the temporal dynamics of pollutant concentrations, several analytical methods have been developed and applied in recent air quality studies. Due to its non-parametric nature and robustness to non-normal distributions, the Mann-Kendall test has been widely used for trend detection in long-term pollutant data (Güçlü et al., 2019; Aladağ, 2023). The innovative trend analysis (ITA), developed by Şen (2012), offers a graphical approach that can detect linear and nonlinear trends without statistical assumptions and has been applied in megacities such as Istanbul (Güçlü et al., 2019). Wavelet transform methods allow for the identification of periodic components and changes at various temporal scales, and have proven particularly effective in capturing short-term and seasonal behaviors in pollutants such as PM_2.5 and O₃ (Sahu et al., 2019; Aladağ, 2024). For forecasting, the Prophet model developed by Facebook has shown high flexibility and performance in time series with strong seasonality and irregular trends (Hasnain et al., 2022). Each of these methods provides complementary insights, and their integration allows for a more comprehensive characterization of air quality evolution.

This study is original in that it integrates non-parametric statistical trend tests, multiscale signal decomposition, and machine learning-based forecasting into a unified framework for analyzing air pollution. To our knowledge, this is the first study to apply this combined methodology to the case of Puebla, Mexico, a mid-sized urban center with significant industrial activity but limited scientific documentation. The study contributes to a broader understanding of pollution dynamics in Mexican cities beyond the well-studied metropolitan areas, such as Mexico City, Monterrey, or Guadalajara.

The purpose of this study is to conduct a detailed analysis of the main air pollutants identified in the city of Puebla in recent years. To achieve this, we employed statistical and signal processing techniques, as well as a machine learning model for pollutants whose concentration has shown an increasing trend.

2. Data

Air pollutants are monitored hourly by the Sistema Nacional de Información de la Calidad del Aire (National Air Quality Information System, SINAICA), which has five stations covering the metropolitan area. For this study, we selected the central region of Puebla, using data from three monitoring stations: Las Ninfas (NIN), Benemérito Instituto Normal del Estado (BINE), and Universidad Tecnológica de Puebla (UTP), covering the period from January 2016 to April 2024 (Fig. 1b). The pollutants analyzed include O₃, PM₁₀, fine particulate matter (PM_2.5), SO₂, and NO₂. These data were obtained from the public access site https://sinaica.inecc.gob.mx (3043 daily-averaged observations). The percentage of null data was not greater than 5% for each station. Measurement gaps were addressed using the multifractal interpolation method proposed by Diosdado et al. (2013), which uses local multifractal interpolation based on temporal patterns in pollutant series, allowing the imputation of missing values without altering the long-term structure or introducing artificial periodicities. It has been shown to preserve statistical properties such as variance and autocorrelation.

These three monitoring stations were selected due to their data completeness, temporal continuity from 2016 to 2024, and their strategic spatial distribution across the Puebla Metropolitan Area. The NIN station is situated in a residential neighborhood near the city center; BINE is located in a peri-urban area influenced by both commercial and residential activities; and UTP lies within an industrial corridor. This selection provides a comprehensive representation of urban, industrial, and transitional environments within the dataset.

Figure 1 provides geographic and topographic context, including the position of Puebla relative to Mexico City, its surrounding volcanic features, and the location of the monitoring stations.

3. Methodology

Several statistical and analytical methods were applied to assess pollutant trends and better understand their temporal and seasonal variations.

3.1 Mann-Kendall test

The non-parametric Mann-Kendall (MK) test was used to identify significant trends in the pollutant data time series (Gilbert, 1987). This methodology is widely used in trend analysis due to its robustness and simplicity. The test statistic (s) of the series p₁, p₂, p₃,…, p_n is estimated using Eq. (1), where n represents the total observations, and p_k and p_j denote the observed data at time k and j, respectively,

(1)

where the sgn function is defined as:

(2)

The variance of S is estimated as follows:

(3)

where, g is the number of the linked group, and T_j refers to the extent of the j-th linked number. From S and VAR(S), the standardized test measure is calculated using Eq. (4).

(4)

A positive value of Z indicates increasing trends, while a negative value indicates decreasing trends in the time series data.

3.2 Innovate trend analysis

The innovative trend analysis (ITA) method was introduced by Şen (2012) to evaluate trends in different time series. This method does not require normality, serial autocorrelation, or specific data length. Additionally, it is robust to outliers and autocorrelated data. The ITA method has been applied in urban air quality studies, including major metropolitan areas such as Istanbul, where it has proven effective for combining trend diagrams with air quality index evaluations (Güçlü et al., 2019). In this method, the given time series data is divided into two equal subseries, and then each subseries is arranged in ascending order independently. The subseries are plotted against each other, with the first half of the subseries on the x-axis and the second half on the y-axis to obtain a scatter diagram. Finally, the 1:1 straight line is drawn. Data points above this line indicate increasing-positive trends, while those below indicate decreasing trends (Şen, 2017). The ITA indicator (ITA_ind) is computed as follows:

(5)

where p_i and p_k indicate the values of the first and second half, n is the extent of each half, and p̅₁ is the mean of the first half. The slope (m_s) of the series can be calculated using the following equation:

(6)

where p̅₂ is the mean of the second half.

3.3 Wavelet transform

In recent years, the wavelet transform theory has been widely investigated due to its applications in various disciplines, including physics, numerical analysis, signal processing, probability, and statistics, among others. Its usefulness lies in its ability to approximate a function or signal while retaining spatial information using a set of functions known as wavelets. Additionally, the wavelet transform is employed to extract localized features of interest in a signal, resulting in improved data compression. Hybrid approaches combining wavelet analysis with traditional time-series forecasting methods have shown promise in improving prediction accuracy for PM concentrations (Aladağ, 2021).

The continuous wavelet transform of a function f(t) is defined as the convolution of f(t) with an analyzer function ψ(σ). To be considered a wavelet, the function must be localized in time as well as in frequency space; moreover, it must be an integrable function, that is, with zero mean (Farge, 1992). It is also assumed that ψ is normalized, that is, , where ψ* is the complex conjugate. For a scale s, location u, and time t, the variable

(7)

can be seen as a dimensionless time scale. For a given time function f(t), the continuous wavelet transform is obtained by

(8)

where the factor 1/ is necessary to satisfy the normalization condition.

The set of basis functions is derived from a single function called the wavelet mother function, which is the function that undergoes modifications to perform the analysis; it is expanded or compressed, and translated along the signal. These modifications occur through the scaling and displacement parameters. Scaling stretches or compresses the wavelet, which is lengthened or compressed, allowing both detailed and global components of the signal to be observed. Displacement refers to the shifting of the wavelet along the signal.

This paper concentrates on the three different mother wavelets given in Torrence and Compo (1998):

The Morlet wavelet (Fig. 2a):

(9)

The Paul wavelet (Fig. 2b):

(10)

And the derivative of Gaussian (DOG) wavelet (Fig. 2c):

(11)

Fig. 2. (a) Morlet wavelet, (b) Paul wavelet, and (c) derivative of Gaussian (DOG) wavelet in the time domain.

The value of k controls the number of oscillations present in the mother wavelet and will therefore strongly influence the frequency and time resolution of the corresponding wavelet transform. The Morlet wavelet has a reasonably large number of oscillations, which will ensure good frequency resolution. The Paul wavelet has far fewer oscillations but is highly localized in time. This will give it a very fine time resolution and at the same time a reduced frequency resolution. The DOG wavelet has relatively few oscillations over a much larger time domain. Note that both the Morlet and Paul wavelets are complex-valued, while the derivative of the Gaussian wavelet is real-valued. Wavelet analysis was used to extract localized features of interest in the signal, allowing for better data compression and the identification of seasonal cycles and variation patterns. Wavelet transform methods allow for the identification of periodic components and changes at different temporal scales, and have proven particularly effective in capturing short-term and seasonal behaviors in pollutants such as PM_2.5 and O₃ (Sahu et al., 2019).

The three wavelet functions were selected for their complementary analytical capabilities, particularly in the context of scalogram representation, which visualizes how signal energy is distributed across time and frequency. The Morlet wavelet, with its strong frequency localization, excels at identifying long-term periodicities in the data and appears in the scalogram as persistent patterns across broader time intervals and lower frequencies. In contrast, the Paul wavelet provides high temporal resolution, making it suitable for detecting transient or short-duration events, which are represented in the scalogram as sharp, localized bursts of energy at higher frequencies. Lastly, the DOG wavelet is sensitive to sudden transitions and is effective for detecting abrupt changes or edges in the time series, which manifest in the scalogram as sharp contrasts in energy distribution. Together, these wavelets allow for a multiscale decomposition of the pollutant time series, effectively capturing dynamic atmospheric processes under varying environmental conditions.

3.4 Prophet forecasting model

The Prophet forecasting model (PFM) is a time series regression tool developed by the Core Data Science team at Facebook (Taylor and Letham, 2017). It is specifically designed to model daily data with strong seasonal components, such as air pollution concentrations in environmental applications. PFM decomposes the time series into three key components (trend, seasonality, and special events), resulting in a flexible and interpretable model. The additive form of the model is:

(12)

Here, y(t) represents the predicted value obtained from either a linear or logistic equation. The functions g(t) and s(t) capture the seasonality or time series patterns based on yearly, monthly, daily, or other periodic cycles, h(t) accounts for outliers related to holidays, and ϵ(t) represents the random or unexpected error. A distinctive feature of PFM is its ability to incorporate trend change points, which mark periods when the underlying behavior of the time series shifts, often due to interventions, holidays, or exogenous events. These change points can be specified manually or detected automatically. To avoid overfitting, Prophet applies L1 regularization to reduce a large set of potential change points to only the most relevant.

Thanks to its capacity to handle structured seasonality and incorporate external regressors, PFM has been successfully used to forecast air pollutant concentrations (e.g., Shen, 2020; Hasnain et al., 2022). In this study, we trained the model using data from 2016 to 2021 and evaluated it using forecasts for the year 2022. As part of the modeling, we included official Mexican holidays and a COVID-19 restriction regressor based on the epidemiological semaphore (ES) system, to account for the impact of pandemic-related lockdowns and reduced human activity on pollutant levels.

All data processing, statistical analyses, and modeling procedures were conducted using the Python (3.11) programming language with the following libraries: PyWavelets, fbprophet, and statsmodels.

4. Results and discussion

4.1 Descriptive statistics

Table I presents the summary statistics of the pollutants. Positive skewness is exhibited by the five series, indicating a trend to extremely high values of the pollutants rather than a trend showing reduced values. Fat tails are present in the five series as shown in terms of kurtosis and Shapiro-Wilk statistics. Figure 3 exhibits the time series of the five pollutants in the period from January 2016 to June 2024. All pollutants showed an annual cycle. The maximum and minimum of the O₃ cycle are found for the winter-spring and summer-autumn seasons, respectively. The green line represents the average concentration. The time series of O₃ concentration in Puebla City indicates sustained oscillations, suggesting stationarity, i.e., increases and decreases occur in a predictable seasonal cycle pattern. In addition, elevated O₃ levels tend to form during periods of warm temperatures, so health-hazardous levels are common during the summer season. These levels can increase during the warm season due to several interrelated factors. O₃ is mainly formed through chemical reactions involving sunlight, nitrogen oxides (NO_x), and volatile organic compounds (VOCs).

Table I. Summary statistics of the five air pollutants studied.

Pollutant	Minimum (μg m–3)	Mean (μg m–3)	Maximum (μg m–3)	Standard deviation (μg m–3)	Skewness	Kurtosis	Shapiro-Wilktest
O3	0.014	0.044	0.103	0.013	0.956	1.192	0.948(000)
PM10	5.739	41.575	109.680	15.498	0.597	0.372	0.978(000)
PM2.5	1.666	17.122	59.158	7.403	0.995	1.665	0.947(000)
SO2	0.0008	0.003	0.010	0.001	0.987	1.843	0.952(000)
NO2	0.007	0.025	0.059	0.007	0.657	0.446	0.973(000)

 

Fig. 3. Daily variations of five air pollutants in the Puebla City center region. (a) Ozone, (b) particles with aerodynamic diameter smaller than 10 μm, (c) particles with aerodynamic diameter smaller than 2.5 μm, (d) sulfur dioxide, and (e) nitrogen dioxide. The vertical line denotes the mean of the time series.

Evaluating distributional characteristics, such as skewness, kurtosis, and the results of the Shapiro-Wilk normality test, is essential for determining whether the assumptions of classical statistical models are met. In the context of air pollutant data, which often exhibit marked positive skewness, high kurtosis, and deviations from normality, these metrics help reveal underlying asymmetries, tailing behavior, and the presence of outliers. Such deviations can distort parametric analyses that assume normality, leading to misleading results. The Shapiro-Wilk test, in particular, offers a formal assessment of normality. Given the observed non-normality in the data, the use of non-parametric methods such as the Mann-Kendall trend test and ITA is statistically justified, as these approaches are more robust to outliers, skewness, and heavy-tailed distributions, ensuring more reliable inferences in environmental time series analysis.

4.2 Trend analysis

Table II shows the trends obtained from the Mann-Kendall test. The p-value of the test is less than 0.05 for all pollutants, indicating enough evidence to conclude that a trend is present in the. A positive Z value indicates increasing trends for O₃, SO₂ and NO₂ while a negative value of Z indicates decreasing trends for PM₁₀, and PM_2.5. Using the ITA method, we can contrast this fact in the graphs (Fig. 4, left side), where the pollutants with negative trends appear below the 1:1 line, and those with positive trends appear above it. Finally, using Eq. (6), a straight line is obtained that illustrates the trend in the time series data of the pollutants (Fig. 4, right side), reaffirming the results given by the Mann-Kendall test.

Table II. Summary of the Mann-Kendall test for the five air pollutants studied.

Pollutant	Trend	p-value	Z	τ	S	Var(s)
O3	Increasing	0.00	17.59	0.21	984559	3.132395e9
PM10	Decreasing	0.00	–7.79	–0.09	–436067	3.132397e9
PM2.5	Decreasing	0.01	–2.53	–0.03	–141998	3.132386e9
SO2	Increasing	0.00	15.53	0.18	869403	3.132395e9
NO2	Increasing	0.00	17.23	0.20	964650	3.132391e9

 

Fig. 4. Left column: scatter plots from the innovative trend analysis (ITA) method comparing the first and second halves of each pollutant time series. Right column: daily time series of each pollutant with two overlaid trend lines: the red line represents the linear slope obtained from ITA, while the yellow line shows a 10-day moving average, used here as a simple low-pass filter to reduce high-frequency variability and highlight long-term trends in the data.

4.3 Wavelet analysis

4.3.1 Ozone (O₃)

Figure 5 presents the scalograms of the Morlet, Paul, and DOG wavelet spectra of O₃ and PM₁₀. In the case of O₃, the Morlet wavelet is used to generate the scalogram in Figure 5a. This complex wavelet is useful for detecting periodic patterns or dominant frequencies in the data. The x-axis is the period (in days) of the original signal corresponding to the wavelet scale on the y-axis, related to frequency. Wavelets with higher frequencies are plotted on lower scales and vice versa. The intensity or color indicates the magnitude of the transformation at a specific point in time and wavelet scale. The deeper the color or the higher the value, the higher the signal energy at that location. Therefore, the yellow contour encloses regions with confidence greater than 95%. The largest yellow band is observed in the 360 days, indicating that we have a highly correlated signal, i.e., a pattern of strong influence. The hatched regions at each end indicate the “cone of influence” where edge effects become important. It is essential to note that, from March 2020, the city of Puebla experienced a partial or total shutdown of factories, in addition to a lockdown due to the COVID-19 pandemic. However, this analysis showed that the O₃ concentration did not decrease despite the measures imposed in response to the health crisis. A yellow coloration is even observed from 2020 to 2023, in approximately 360 days. In addition, high-energy bands can be visualized at smaller scales. The right side of Figure 5a shows the Morlet wavelet spectrum and the Fourier spectrum. Here, the frequencies of the signals and their amplitude are shown. The peak of the dominant frequency is observed at 360 days. While the Fourier spectrum describes the frequencies present in a signal globally, the wavelet spectrum provides a time- and frequency-focused representation, making it easier to identify signal characteristics at various times and time scales.

Fig. 5. Scalograms of O₃ and PM₁₀ using the (a, d) Morlet, (b, e) Paul, and (c, f) derivative of Gaussian (DOG) wavelets.

Figure 5b shows that the Paul wavelet reflects a high-energy region on a 360-day scale. The cone of influence is narrower, suggesting that local perturbations have a more limited effect on the outcome of the analysis. Being a complex function that takes complex values instead of real values provides a more complete representation of the signal, as it includes more information about amplitude and phase. Here, the Fourier spectrum shows two dominant peaks of energy, the predominant peak at about 360 days and a fainter peak at 130 days.

Figure 5c presents the analysis performed on the time signal using the mother wavelet DOG. The most notable difference is the fine-scale structure using DOG. This might be because DOG is a real function that captures the positive and negative oscillations of the time series as separate peaks in the wavelet power. Morlet and Paul wavelets are complex functions containing more oscillations than DOG, and hence, the wavelet power combines positive and negative peaks into a single broad peak. Overall, the figures shown contain similar power features on the 360-day scale. The DOG wavelet is narrower in space-time but broader in spectral domain than the Morlet wavelet. In Figure 5c, the peaks appear sharper in time and more elongated in scale.

4.3.2 Particles smaller than 10 microns (PM₁₀)

Regarding the Morlet wavelet (Figure 5d), we note that at the 360-day scale, there is a persistent high-energy band throughout the entire period. Additionally, significant autocorrelations are observed at scales shorter than eight days across the entire period. However, the use of this wavelet does not offer relevant information for this pollutant. The Paul wavelet exhibits a high-energy band throughout the entire period, spanning approximately 360 days, which weakens at the beginning of 2020 and again at the beginning of 2022. This band shows a significant autocorrelation peak around 32 days in mid-2019. Furthermore, significant autocorrelations are present at scales below eight days throughout the entire period. Finally, in the scalogram corresponding to the DOG wavelet, we observe a high energy concentration at approximately semi-annual periods, ranging from 200 to 500 days.

4.3.3 Particles smaller than 2.5 microns (PM_2.5)

Regarding the Morlet wavelet (Fig. 6a), we note that on a 360-day scale, there is a narrow high-energy band throughout the period, which decreases its power between 2020 and 2022. Additionally, significant autocorrelations are observed at scales shorter than eight days throughout the entire period. The Paul wavelet exhibits a high-energy band around 360 days that begins to extend to very large scales starting in 2019. In the scalogram corresponding to the DOG wavelet, we observe a high-energy concentration in approximately semiannual periods ranging from 200 to 500 days and expanding to larger scales starting in 2019. Note that in all three scalograms, an important peak in signal energy appears in 2019.

Fig. 6. Scalograms of PM_2.5 and SO₂ using the using the (a, d) Morlet, (b, e) Paul, and (c, f) derivative of Gaussian (DOG) wavelets.

4.3.4 Sulfur dioxide (SO₂)

Figure 6a-c presents the analysis for the SO₂ pollutant with the Morlet, Paul, and DOG wavelets, respectively. With the Morlet wavelet, we can highlight that on a 360-day scale, there is a high-energy band in time from 2019 to 2023; this band extends to the cone of influence, which is the area where the wavelet has a significant influence on the reconstruction of the original signal. Furthermore, it is noteworthy that we observe significant frequencies at both high and low scales. Recalling that higher scales represent lower frequency components and vice versa. The Paul wavelet shows some high energy zones between 2019 to 2020 and between 2021 to 2022. This energy extends from a scale of 64 to 512 days to the cone of influence. It is worth noting that the signal also exhibits high energy on smaller scales. We notice that the 95% confidence band of the spectrum is limited, which may be due to the fact that the Paul wavelet is a real wavelet that is based solely on magnitude and does not consider phase, a crucial aspect for understanding the complete structure of the signal. Finally, the DOG wavelet derived from the difference between two Gaussian functions, is useful for detecting edges and local features. It is relatively simple and effective for identifying abrupt transitions. Here, we see that the SO2 concentration increased over a period of 360 days. It is important to note that this wavelet reveals high energy at specific points in time.

The presence of high-energy zones in SO₂ scalograms, particularly between 2019 and 2023, may be associated with seasonal changes in industrial activity, use of fossil fuels during colder months, or meteorological conditions such as wind stagnation or temperature inversions. However, additional data on emissions inventories would be necessary to confirm these associations.

4.3.5 Nitrogen dioxide (NO₂)

Figure 7 shows the analysis performed for the NO₂ pollutant with the Morlet, Paul, and DOG wavelets. From here, we can highlight that on a scale greater than 256 days, there is a high-energy band throughout the entire period, increasing the power of the signal from 2021 onwards. Figure 7b shows that the Paul wavelet indicates a high concentration of NO₂ at small scales. We can see that there is no autocorrelation at scales less than 1000 days during the periods 2019-2021 and 2022-2023. High concentration levels were present on scales greater than 16 days starting in 2021. For the DOG wavelet, we can observe the lack of autocorrelation at scales of less than 1000 days during the period 2019-2021. A significant peak of autocorrelation can be observed at the beginning of 2022 on all scales less than 512 days.

Fig. 7. Scalograms of NO₂ using the (a, d) Morlet, (b, e) Paul, and (c, f) derivative of Gaussian (DOG) wavelets.

The marked NO₂ peak observed in early 2022 may relate to the post-lockdown economic recovery and increased vehicular activity, particularly after the relaxation of COVID-19 restrictions. This rebound effect has been documented in multiple cities during 2022, where pollutant concentrations rose again after confinement periods due to the resumption of transport and industrial emissions (Dong et al., 2022).

4.4 PM₁₀ and PM_2.5 forecast

It is worth noting that in the scalograms for PM₁₀ and PM_2.5 (shown in the previous section), a prominent peak in energy appears in 2019. This is possibly because the Spearman correlation coefficient r between these two pollutants is 0.8, while the correlation between PM₁₀ and the other pollutants remains close to r = 0. The existence of this peak coincides with the fact that the highest concentrations of PM_2.5 were recorded in May 2019 (Fig. 8b). In Mexico, the Official Mexican Standard NOM-025-SSA1-2014 (SSA, 2014) establishes that the concentration limit for PM_2.5 is 30 μg m^–3, although the WHO recommends 25 μg m^–3 (Fig. 8a). For this reason, in this section our primary objective is to model concentrations of PM_2.5 and PM₁₀ air pollutants.

Fig. 8. (a) PM_2.5 time series. The horizontal dotted red line denotes the maximum daily exposure according to NOM-025-SSA1-2014 guidelines. (b) The days in which the maximum daily exposure was exceeded are shown in red.

To model both short-term and long-term trends in PM concentrations, we incorporated official Mexican holidays and a COVID-19-related regressor into the PFM. This regressor accounts for the impact of lockdown periods enforced via Mexico’s ES system, which categorized regions by color based on the severity of the COVID-19 spread. During red and orange phases, significant reductions in mobility, industrial activity, and emissions were observed, particularly affecting pollutants such as PM₁₀ and PM_2.5. Incorporating this regressor allows the model to isolate the effect of these anomalies and yield more reliable trend estimations.

The PFM was trained using data from 2016 to 2021 and evaluated with forecasts for the period 2022-2024. Figure 9 presents the observed and predicted values for (a) PM_2.5 and (b) PM₁₀. The training data (2016-2021) are shown in black, the test data (2022-2024) in red, and the predicted trend in blue, with a shaded area representing the model’s uncertainty intervals. Including the full-time span in the figure provides a clearer visualization of both models’ fit and forecast behavior.

Fig. 9. Forecast results using the Prophet model for (a) PM_2.5 and (b) PM₁₀. Black dots indicate training data (2016-2021), red dots represent test data (2022-2024), and the blue line corresponds to model predictions with uncertainty intervals in light blue.

The model’s performance was evaluated using the correlation coefficient for 2022 (after training), and the results are presented in Table III for PM_2.5 and PM₁₀. As shown, the PFM exhibits weak performance in the short term, with a very low correlation between actual and predicted values for both PM_2.5 and PM₁₀. In this case, the model cannot accurately capture the short-term fluctuations and dynamics of these pollutants, due to the complex and rapidly changing factors. Over six months, the correlation increases to moderate levels, indicating that the model can capture more of the underlying trends in both PM_2.5 and PM₁₀ for long-term patterns. For PM₁₀, the PFM shows a strong prediction ability over nine to 12 months, with correlations values ranging 0.61 to 0.64. This suggests that the model is more effective at capturing long-term trends for PM₁₀, due to more stable temporal patterns in its behavior. On the other hand, for PM_2.5, the correlation remains moderate (around 0.40) at the nine-month horizon, but decreases to 0.30 over a full year, indicating that the PFM for PM_2.5 is less reliable over extended periods.

Table III. Correlation coefficient R between observed and predicted concentrations of PM_2.5 and PM₁₀ using the Prophet model over increasing forecast horizons.

Time	PM2.5	PM10
30 days	0.03	0.05
90 days	0.13	0.11
180 days	0.40	0.42
270 days	0.40	0.64
360 days	0.30	0.61

 

According to our calculations, the PFM performs better for PM₁₀ than PM_2.5 over longer periods, implying that the utility of the model is more suitable for medium- to long-term forecasting for PM₁₀. Additionally, short-term predictions for both pollutants should be interpreted with caution due to weak correlations.

The improved performance of the PFM over longer time horizons can be attributed to its ability to capture seasonal and trend components more effectively than short-term high-frequency variability. Short-term fluctuations in pollutant concentrations are often driven by transient meteorological events and localized emissions, which are harder to model without high-frequency meteorological covariates. This behavior has been observed in similar air quality forecasting applications using Prophet (Shen et al., 2020).

Finally, we use the PFM to observe the underlying trend. The analysis of air quality data using the PFM revealed significant trends for both PM_2.5 and PM₁₀ concentrations. After isolating the trend components of the model, the trend analysis showed a clear decreasing pattern over the study period. We computed the average rate of change for the test dataset starting from 2022, obtaining –0.0031 and –0.0053 for PM_2.5 and PM₁₀, respectively. In Figure 10, the trend component of the model indicates a consistent decline in PM_2.5 and PM₁₀ concentrations, in accordance with the results of previous sections for PM₁₀ and PM_2.5.

Fig. 10. Decreasing trends for PM_2.5 and PM₁₀ using the Prophet forecasting model.

4.5 Discussion

This study provides a comprehensive temporal characterization of air pollutants in Puebla City by combining non-parametric trend analysis, wavelet-based spectral decomposition, and time-series forecasting through the PFM. The results obtained reveal important spatial and temporal patterns that align with, but also extend beyond, findings in other Mexican urban contexts.

The trend analysis, using the Mann-Kendall and ITA methods, showed significant increases in O₃, SO₂, and NO₂, while PM₁₀ and PM_2.5 levels decreased, consistent with post-pandemic transitions and structural changes in emissions. A similar behavior has been reported in Monterrey and Mexico City (Peralta et al., 2021; Fajardo-Montiel y Ramírez-Sánchez, 2023), where economic reactivation was accompanied by a partial rebound in gaseous pollutants but not in PM, possibly due to long-term emission control strategies.

The wavelet transform analysis confirmed that pollutant behavior in Puebla exhibits seasonal cycles, especially at scales near 360 days, with additional high-frequency activity (e.g., <16 days) detected in pollutants like O₃ and NO₂. These findings suggest that meteorological variability, industrial activity, and vehicular emissions have a strong influence on daily and seasonal fluctuations in pollutant levels. The presence of energy peaks in 2019 across multiple wavelet scales is particularly noteworthy and may reflect regional pollution episodes or meteorological anomalies.

A key contribution of this study is the integration of the PFM. While short-term correlations were weak due to transient fluctuations and lack of high-frequency covariates, the model showed better predictive capacity over longer horizons, especially for PM₁₀. This aligns with results from Shen et al. (2020) and Hasnain et al. (2022), who also noted the PFM’s improved performance in medium and long-term forecasting. In recent studies, hybrid forecasting models based on Empirical Mode Decomposition and ARIMA have also demonstrated potential for improving air quality prediction, particularly for complex pollutant patterns (Aladağ, 2023).

Moreover, the inclusion of COVID-19 restriction phases as external regressors enhanced the interpretability of trends. The pandemic offered a natural experiment to assess how mobility and activity restrictions affect emissions. Despite reduced activity during lockdowns, O₃ levels remained elevated, suggesting the influence of photochemical mechanisms and meteorological constraints, as previously observed in Mexico City (Peralta et al., 2021). Recent studies have used wavelet coherence analysis to explore the dynamic relationship between air pollution and the COVID-19 pandemic, revealing temporal dependencies between lockdown measures and pollutant levels (Aladağ, 2024).

However, the study also presents limitations. The use of only three monitoring stations may limit spatial generalizability, and the absence of meteorological data (e.g., wind, humidity, temperature) constrains the ability to explain variability. The PFM, although useful, does not capture short-term pollution spikes or dynamic interactions between emissions and weather.

Overall, these findings highlight the importance of using multiscale and hybrid methodologies to assess air quality in mid-sized cities like Puebla, where rapid industrial growth and urbanization may outpace traditional monitoring and policy responses. Future studies should incorporate more granular data, advanced dispersion models, and integrate remote sensing and satellite data to improve forecasting and support an evidence-based environmental policy.

5. Conclusion

This study highlights the dynamics of air pollution in Puebla City and its implications for public health and environmental policy. Key findings indicate significant increases in ozone (O₃), sulfur dioxide (SO₂), and nitrogen dioxide (NO₂) levels, while particulate matter (PM₁₀ and PM_2.5) and carbon monoxide (CO) levels have shown a decreasing trend. These patterns reflect a complex interaction between anthropogenic emissions, urban activity, and atmospheric processes. While increases in some pollutants may be related to the resumption of industrial operations and vehicular traffic, the decline in PM concentrations suggests the influence of both temporary and long-term structural factors.

The COVID-19 lockdowns in early 2020 contributed to short-term reductions in air pollutant concentrations, particularly in NO₂ and PM_2.5. However, the continued downward trend in PM₁₀ and PM_2.5 beyond the pandemic period may be indicative of lasting improvements in emissions control, such as stricter regulations, cleaner transportation technologies, and improved industrial practices.

Similar behaviors have been observed in other urban centers, where O₃ levels did not decrease significantly during periods of reduced activity, highlighting the role of photochemical reactions and meteorological conditions in ozone formation. These dynamics are probably also relevant in Puebla, underscoring the need for integrated pollution control strategies that consider local chemical and climatic processes.

The wavelet analysis performed in this study revealed the presence of high-frequency energy bands and seasonal cycles in pollutant concentrations, suggesting that small meteorological variations can significantly influence air pollution levels in Puebla, as has been shown in other major cities.

This study has some limitations, including the use of only three monitoring stations, which may not fully capture the spatial variability of air pollution in the city. Additionally, the absence of meteorological variables such as temperature, wind speed, precipitation, and humidity limits the ability to analyze atmospheric dispersion and chemical transformation processes. Phenomena such as thermal inversions and weekly activity cycles were also not explicitly addressed. Furthermore, while the PFM was effective in capturing long-term trends, it lacked accuracy in short-term forecasts due to the absence of high-frequency covariates. Future studies should incorporate high-resolution meteorological data and dynamic dispersion models to better understand the seasonal and topographic influences on air quality in the region.

To address the challenges associated with poor air quality, it is essential to implement a dynamic and responsive management strategy. This should include precise identification of pollution sources through spatial analysis, comprehensive health impact assessments, and the deployment of smart monitoring systems for real-time data acquisition and response.

Policymakers should use these findings to strengthen environmental regulations, promote sustainable industrial and urban development, and increase public awareness regarding the sources and effects of air pollution. Continued research and policy action are key to effectively mitigating the environmental and health impacts of urban air pollution.

References

Aladağ E. (2021). Forecasting of particulate matter with a hybrid ARIMA model based on wavelet transformation and seasonal adjustment. Urban Climate 39: 100930. https://doi.org/10.1016/j.uclim.2021.100930

Aladağ E. (2023). Prediction of air pollution in Van province using ARIMA model based on empirical mode decomposition. Yuzuncu Yil University Journal of the Institute of Natural and Applied Sciences 28: 495-509. https://doi.org/10.53433/yyufbed.1220578

Aladağ E. (2024). The nexus between air pollution and the COVID-19 pandemic in Turkey: Further insights from wavelet coherence analysis. Aerosol Science and Engineering 8: 108-119. https://doi.org/10.1007/s41810-023-00209-1

Brunekreef B, Holgate ST. 2002. Air pollution and health. The Lancet 360: 1233-1242. https://doi.org/10.1016/S0140-6736(02)11274-8

De Foy B, Varela JR, Molina LT, Molina MJ. 2006. Rapid ventilation of the Mexico City basin and regional fate of the urban plume. Atmospheric Chemistry and Physics 6: 2321-2335. https://doi.org/10.5194/acp-6-2321-2006

Díaz-Torres JJ, Ojeda-Castillo V, Hernández-Mena L, Vergara-Sánchez J, Saldarriaga-Noreña HA, Murillo-Tovar MA. 2022. Long-term analysis of tropospheric ozone in the urban area of Guadalajara, Mexico: A new insight of an alternative criterion. Atmosphere 13: 152. https://doi.org/10.3390/atmos13020152

Diosdado AM, Coyt GG, López JB, del Rio Correa JL. 2013. Multifractal analysis of air pollutants time series. Revista Mexicana de Física 59: 7-13. https://www.redalyc.org/articulo.oa?id=57030970002

Dong X, Zheng X, Wang C, Zeng J, Zhang L. 2022. Air pollution rebound and different recovery modes during the period of easing COVID-19 restrictions. Science of The Total Environment 843: 156942. https://doi.org/10.1016/j.scitotenv.2022.156942

Fajardo-Montiel LA, Ramírez-Sánchez HU. 2023. Evolution of air pollution in the Monterrey metropolitan area, Mexico. Asian Journal of Environment & Ecology 21: 55-75. https://doi.org/10.9734/ajee/2023/v21i1453

Farge M. 1992. Wavelet transforms and their applications to turbulence. Annual Review of Fluid Mechanics 24: 395-458. https://doi.org/10.1146/annurev.fl.24.010192.002143

Gilbert RO. 1987. Statistical methods for environmental pollution monitoring. Wiley, NY, 320 pp.

Gómez-Tuena A, Orozco-Esquivel MT, Ferrari L. 2007. Igneous petrogenesis of the Trans-Mexican Volcanic Belt. In: Geology of Mexico: Celebrating the centenary of the Geological Society of Mexico (Alaniz-Álvarez SA, Nieto-Samaniego AF, Eds.). Geological Society of America. https://doi.org/10.1130/2007.2422(05)

Güçlü YS, Dabanlı İ, Şişman E, Şen Z. (2019). Air quality (AQ) identification by innovative trend diagram and AQ index combinations in Istanbul megacity. Atmospheric Pollution Research 10: 88-96. https://doi.org/10.1016/j.apr.2018.06.011

Hasnain A, Sheng, Y, Hashmi MZ, Bhatti UA, Hussain A, Hameed M, Marjan S, Bazai SU, Hossain MA, Sahabuddin M, Wagan RA, Zha Y. 2022. Time series analysis and forecasting of air pollutants based on Prophet forecasting model in Jiangsu Province, China. Frontiers in Environmental Science 10: 945628. https://doi.org/10.3389/fenvs.2022.945628

Hernández-Paniagua IY, Clemitshaw KC, Mendoza A. 2017. Observed trends in ground-level O3 in Monterrey, Mexico, during 1993-2014: Comparison with Mexico City and Guadalajara. Atmospheric Chemistry and Physics 17: 9163-9185. https://doi.org/10.5194/acp-17-9163-2017

INECC. 2024. Informe nacional de la calidad del aire. Instituto Nacional de Ecología y Cambio Climático, Mexico. https://sinaica.inecc.gob.mx/archivo/informes/Informe2021.pdf

Limón-Sánchez MT, Smith W. 2011. Black carbon in PM_2.5 data from two urban sites in Guadalajara, Mexico during 2008. Atmospheric Pollution Research 2: 358-365. https://doi.org/10.5094/APR.2011.040

Molina LT, Molina M.J. 2002. Air quality in the Mexico megacity: An integrated assessment. Springer Dordrecht. https://doi.org/10.1007/978-94-010-0454-1

Molina LT, Molina MJ. 2004. Improving air quality in megacities: Mexico City case study. Annals of the New York Academy of Sciences 1023: 142-158. https://doi.org/10.1196/annals.1319.006

Peralta O, Ortínez-Álvarez A, Torres-Jardón R, Suárez-Lastra M, Castro T, Ruiz-Suárez LG. 2021. Ozone over Mexico City during the COVID-19 pandemic. Science of The Total Environment 761: 143183. https://doi.org/10.1016/j.scitotenv.2020.143183

Sahu SK, Tyagi B, Pradhan C, Beig G. 2019. Evaluating the variability, transport and periodicity of particulate matter over smart city Bhubaneswar, a tropical coastal station of eastern India. SN Applied Sciences 1: 1-12. https://doi.org/10.1007/s42452-019-0427-2

Şen Z. 2012. Innovative trend analysis methodology. Journal of Hydrologic Engineering 17: 1042-1046. https://doi.org/10.1061/(ASCE)HE.1943-5584.0000556

Şen Z. 2017. Innovative trend significance test and applications. Theoretical and Applied Climatology 127: 939-994. https://doi.org/10.1007/s00704-015-1681-x

Shen J, Valagolam D, McCalla S. 2020. Prophet forecasting model: A machine learning approach to predict the concentration of air pollutants (PM_2.5, PM₁₀, O₃, NO₂, SO₂, CO) in Seoul, South Korea. PeerJ 8: e9961. https://doi.org/10.7717/peerj.9961

Sierra A, Vanoye AY, Mendoza A. 2013. Ozone sensitivity to its precursor emissions in northeastern Mexico for a summer air pollution episode. Journal of the Air & Waste Management Association 63: 1221-1233. https://doi.org/10.1080/10962247.2013.813875

SSA. 2014. Norma Oficial Mexicana NOM-025-SSA1-2014. Salud ambiental. Valores límite permisibles para la concentración de partículas suspendidas PM₁₀ y PM_2.5 en el aire ambiente y criterios para su evaluación. Secretaría de Salud, Mexico. Diario Oficial de la Federación, August 20.

Suárez G, Caballero-Jiménez GV, Novelo-Casanova DA. 2019. Active crustal deformation in the Trans-Mexican Volcanic Belt as evidenced by historical earthquakes during the last 450 years. Tectonics 38: 3544-3562. https://doi.org/10.1029/2019TC005601

Taylor SJ, Letham B. 2017. Forecasting at scale. American Statistician 72: 37-45. https://doi.org/10.1080/00031305.2017.1380080

Torrence C, Compo GP. 1998. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society 79: 61-78. https://doi.org/10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2