Air pollution in the shadow of global crises: lessons from a small city in the Western Balkans

Abstract

This study investigates how energy disruption, stemming from COVID-19 pandemic impacts, geopolitical instability, and rapid energy transitions, shaped air quality outcomes in Kičevo, a small city in the understudied pollution landscape of the Western Balkans. The research analyses a five-year dataset from 2019 to 2023, encompassing the concentration of atmospheric pollutants (CO, NO₂, O₃, SO₂, PM₁₀, and PM_2.5) and meteorological parameters, to elucidate the complex interactions between emission sources, meteorological conditions, and anthropogenic activities. The air quality data were deweathered using a machine learning method to isolate the effects of emission sources from meteorological influences. The findings reveal significant variations in pollutant concentrations, with notable anomalies observed in 2020 and 2022. These anomalies were primarily driven by elevated activities in both near and distant power plants in the area, corresponding to shifts in energy production and consumption linked to the COVID-19 crisis. Another contributing factor was the reopening of lignite power plants in neighbouring Greece, undertaken due to instability in renewable energy supply. Greece's strategy of rapidly transitioning to renewables and natural gas proved premature given the lack of adequate storage capacity, leading to a fallback reliance on fossil fuels as a backup plan. The study highlights the impact of socio-economic factors, particularly the substantial demographic declines due to emigration, on reduced emissions from local sources such as residential heating and transportation.

---

Environmental significance

Global disasters such as pandemics, wars and energy insecurity can trigger sudden shifts in emissions with major air quality impacts. Using a five-year dataset from North Macedonia, we show how COVID-19 lockdowns and regional energy shortages drove local and transboundary pollution anomalies, including cross-border pollutant transport. By combining deweathering and trajectory analyses, we highlight the links between energy instability and atmospheric emissions. These findings emphasise that resilient, sustainable energy transitions are vital not only for climate goals but also to protect populations from heightened exposure to harmful pollutants during times of crisis.

1. Introduction

Air pollution represents a significant global threat with far-reaching implications for human health and the environment.¹ It is linked to a myriad of health concerns, including respiratory diseases, cardiovascular disorders, and elevated mortality rates.² Prolonged exposure to pollutants such as particulate matter (PM), nitrogen oxides (NO_x), and sulphur dioxide (SO₂) has been demonstrated to exacerbate chronic illnesses, reduce life expectancy, and impair lung function, particularly in vulnerable populations such as children and the elderly.^3–5 Approximately seven million premature deaths are attributed to poor air quality by the World Health Organization.⁶ Beyond human health, air pollution causes climate change, degrades ecosystems, and diminishes agricultural productivity and biodiversity.^7,8

Given its multidimensional impact, a robust understanding of the sources, dynamics, and mitigation of air pollution is essential to inform effective policy responses. While extensive research has emerged globally, ranging from the evaluation of Clean Air Zones in the UK⁹ to large-scale reviews of policy impacts in rapidly developing urban environments such as the UAE,⁵⁴ Iran,¹⁰ and Chile⁵³ these efforts tend to concentrate on high-income or geopolitically significant regions. Comparatively little attention has been directed toward southeastern Europe and the Western Balkans, where industrial legacy, urbanisation, and weak enforcement of environmental regulations have created a persistent and complex air quality crisis.

The Western Balkans, a subregion in southeastern Europe, comprises countries that are not currently part of the European Union (EU), yet it is of strategic importance due to its location and growth potential. The region in question comprises Albania, Bosnia and Herzegovina, Kosovo, Montenegro, Serbia, and North Macedonia. Since the dissolution of Yugoslavia in the 1990s, the Western Balkans has experienced a series of notable political and economic transitions, conflicts, developments, and setbacks.

The Western Balkans has experienced a notable decline in population, which is substantially influenced by economic failures.¹¹ The implementation of neoliberal economic policies during the 1990s resulted in deindustrialisation and a stagnant gross domestic product (GDP), thereby exacerbating economic instability. These conditions have coincided with high levels of emigration, particularly among younger residents moving abroad in search of better opportunities.¹² The demographic crisis is further intensified by low fertility rates and an ageing population, with projections indicating a significant reduction in the size of the population by 2050.¹³ In addition, the region's economic growth has been constrained by structural weaknesses and external economic shocks, which have dampened trade and investment prospects. As a result, the Western Balkans are confronted with a difficult demographic outlook, grappling with the challenge of retaining their working-age population. However, it is possible that this decline in population may have a beneficial impact on the air quality of the region if air pollution emissions per capita remain similar.

Research focused on air pollution in the Western Balkans remains notably sparse when compared to studies conducted in the rest of Europe, the Middle East, East Asia, and North America. Significant contributions to this field include the work of Belis et al.,¹⁴ who used a combination of chemistry transport models (CTM), local emission inventories, meteorological fields and boundary conditions to simulate the dispersion process and chemical reactions of urban air pollution in the Danube and Western Balkans regions. Results revealed that over 60% of particulate matter with diameters smaller than 2.5 microns (PM_2.5) originated from anthropogenic sources, including energy production (22%), agriculture (19%), residential combustion (16%), and road transport (7%). Belis et al.¹⁵ further assessed the health impacts and costs attributable to urban air pollution in the Western Balkans, studying levels of ozone (O₃), nitrogen dioxide (NO₂), particulate matter with diameters smaller than 10 microns (PM₁₀), and PM_2.5 in 28 urban areas across Albania, Bosnia and Herzegovina, Kosovo, Montenegro, and Serbia. This review of available literature underscores a critical gap in air quality research within the region, particularly regarding the fundamental drivers of air pollution and their consequent impacts.

The city of Kičevo in the Western Balkans (North Macedonia) has experienced several socio-economic, geographical and local governance conditions that have the potential to influence the impact of and response to air pollution.¹⁶ Despite the potential impact of these factors on air quality, there is a notable lack of comprehensive air pollution research studies focused on Kičevo itself. This scarcity of local research reflects a broader trend in the Western Balkans, where air pollution studies are generally limited. While not specific to Kičevo, Samara et al.¹⁷ employed robotic chemical mass balance (RCMB) receptor modelling to estimate the impact of fugitive dust emissions from opencast lignite mines on ambient PM₁₀ levels. Their findings indicate that these mines contribute between 9% and 42% to the PM₁₀ levels in the Western North Macedonia region. Rizos et al.¹⁸ employed Hidden Markov Models to identify the principal sources of PM₁₀ in Western Macedonia (Greece), thereby confirming the existence of significant anthropogenic emissions, particularly from the energy and industrial sectors. Srbinovska et al.¹⁹ investigated the influence of meteorological variables and the spatial distribution of green infrastructure on PM pollution in North Macedonia. To this end, they employed PM data from low-cost sensors installed in Skopje. The results indicated that there was no significant correlation between PM levels and ambient temperature, pressure, and humidity. It is noteworthy that the concentration levels of other criteria air pollutants, including NO₂, SO₂, O₃, and CO, have not been the subject of investigation in either of the studies reviewed above.

This research paper presents a comprehensive analysis of air pollution in the designated study area, building upon the aforementioned rationale. The investigation examines the temporal concentration patterns of criteria air pollutants and interprets the influence of both proximal and distant emission sources. To enhance the precision of our impact assessment, we applied a deweathering technique to the pollutant concentration data to reduce the influence of confounding meteorological variables. This approach isolates the effects of emission source variations and provides a more robust understanding of pollution dynamics, with the subsequent section detailing its theoretical underpinnings and practical application in this study.

Although geographically limited, Kičevo offers a representative setting to explore local air pollution dynamics within a broader regional context marked by comparable environmental and socio-economic challenges. By systematically analysing temporal trends and their driving factors, this research yields nuanced insights into the sources and behaviours of pollution that may resonate across similar urban environments in the Western Balkans, thereby contributing to the limited environmental literature on this under-researched part of Europe.

2. Material and methods

2.1. The study area

Kičevo (41°31′04″N 20°57′56″E) is a city in the western part of North Macedonia, located in a valley with a population of around 40 000 according to the 2021 census.²⁰ The location of the city in Europe and North Macedonia is presented in Fig. 1. The city is situated 112 km from the North Macedonian capital Skopje, 213 km from the Adriatic Sea, and 271 km from the Aegean Sea. Kičevo's economy is primarily driven by industry and mining. The coal-fired power station TEC Oslomej, a substantial source of electricity in the region, is located northeast of the city, as indicated by the blue arrow in Fig. 1.

Fig. 1 The city of Kičevo at North Macedonia within the Western Balkans, Southeastern Europe region. The blue arrows indicate the location of TEC Oslomej power plant.

As is the case with most cities in the Western Balkans, the population of Kičevo has experienced a decline in recent years. The municipality's population decreased from 56 734 in 2002²¹ to 39 669 in 2021, a decline that is mainly attributed to economic factors. The broader economic context of North Macedonia, including Kičevo, is characterised by a GDP of €12.5 billion in 2023, with real GDP per capita at €6618.²²

Kičevo is distinguished by a continental climate, characterised by pronounced seasonal fluctuations. The city typically experiences warm summers and cold winters. Wood burning represents a significant source of heating for the Kičevo residential sector, with this practice occurring from late October through March. The mean temperature in Kičevo ranges from approximately −1 °C in January, the coldest month, to around 24 °C in July, the warmest month. The seasonal changes are pronounced, with significant temperature variations between summer and winter. The Boundary Layer Height (BLH) is closely related to these seasonal temperature variations, and the BLH plays a crucial role in air pollution dispersion. Fig. 2(a) and (b) provide the temporal variation of BLH in Kičevo from 2020 to 2023. The source and methodology for obtaining the BLH data are detailed in a subsequent section.

Fig. 2 The variation of boundary layer height over the studied (a) years and (b) seasons defined as winter: December–February, spring: March–May, summer: June–August, and autumn: September–November).

Fig. 2(a) shows the average daily BLH values over the four-year period. The plot highlights substantial day-to-day variability, with BLH ranging from approximately 100 m to over 1000 m. This indicates considerable fluctuations in the effective atmospheric mixing volume on short timescales. Despite this variability, a recurring annual cycle is observable, with recurring peaks and troughs throughout each year.

Fig. 2(b) presents the seasonal averages of BLH for each year from 2020 to 2023. This plot confirms the seasonal patterns suggested in Fig. 2(a). Summer consistently exhibits the highest average BLH, typically exceeding 600 m, facilitating greater vertical mixing and potential dispersion of pollutants. Spring follows with slightly lower BLH values, generally around 550–600 m. Winter and autumn show the lowest BLH, often below 400 m, which will contribute to the accumulation of pollutants near the surface due to reduced vertical mixing.

2.2. Air quality, meteorology and boundary layer analysis

This study utilises data from the Kičevo air quality monitoring station, situated near the city centre, measuring PM_2.5, PM₁₀, NO₂, SO₂, CO, and O₃, as well as meteorological parameters. These pollutants are monitored using automated instruments consistent with the State Air Quality Monitoring System of North Macedonia: a pulsed fluorescence analyser for SO₂ (Thermo Electron Model 43C), a chemiluminescence analyser for NO/NO₂/NOx (Model 42C), a non-dispersive infrared analyser for CO (Model 48C), a UV photometric analyser for O₃ (Model 49C), and a β-ray absorption monitor for particulates (FH 62 I–R). Meteorological parameters are recorded via integrated sensor modules monitoring conditions such as temperature, humidity, wind speed, and wind direction.⁵⁵ The data for the years 2019 to 2023 were analysed in this study, while the PM_2.5 data were only available from 2021 onwards. These years were selected because consistent data for the analysed pollutants became available starting from 2019, with earlier data either unavailable or incomplete. The station offers a representative sample of the air quality experienced by the local population, situated 75 m from the closest freeway and 7 m from the closest residential road. The monitoring station is situated 8 km to the north of the thermoelectric plant TEC Oslomej. The data is publicly accessible online via the air quality portal, published by the Ministry of Environment and Physical Planning of North Macedonia. No collaboration was undertaken in the collection or verification of the data, and the data were used as presented online. The temporal variations in air quality and meteorological data were analysed using the sophisticated R package OpenAir.²³

This study also used ERA5 reanalysis data to analyse the relationship between the Boundary Layer Height (BLH) and temperature. ERA5, produced by the European Centre for Medium-Range Weather Forecasts (ECMWF), provides high-resolution atmospheric data on a global scale.^24,25 The dataset encompasses a four-year period, with hourly temporal resolution and a spatial resolution of 0.25 ° × 0.25 °. BLH and 2 metre temperature data were extracted from the ERA5 dataset for the study area for the years 2019–2023. The analysis focused on examining the diurnal and seasonal variations of BLH, as well as its relationship with temperature, following methodologies similar to those employed in previous studies.^26,27

2.3. Deweathering and weather normalization

A significant limitation of numerous air quality studies is the failure to consider the intrinsic correlations between air quality and meteorological phenomena. Meteorological conditions have the potential to influence both the relevance of emission sources and pollutant concentrations.^28,29 Fluctuations in weather may serve to either mitigate or amplify the impact of emission alterations on air quality. De-weathering refers to the statistical removal of meteorological influences (such as wind speed, temperature, humidity, and boundary layer height) from air quality datasets. This process isolates the contribution of emissions and long-term trends from short-term variability caused by weather, allowing for a clearer understanding of the effectiveness of pollution control measures and underlying emission patterns. The application of machine learning (ML) methodologies has been identified as a promising approach for the deweathering of air quality datasets.³⁰

To remove the influence of meteorological variability from observed pollutant concentrations, a machine learning–based deweathering approach was applied following the framework proposed by.^29,31 For each pollutant, a gradient boosting machine (GBM) model was trained using meteorological predictors (air temperature, wind speed, wind direction, and month) together with temporal features capturing long-term trend and weekly cycles. Missing values in both predictors and response variables were imputed using an iterative model-based procedure, and all continuous variables were normalised prior to model training. The modelling was implemented using the rmweather package, which constructs an ensemble of decision trees sequentially to minimise prediction error, with hyperparameters tuned to balance predictive accuracy and generalisation.

Deweathered concentrations were derived via a meteorological resampling strategy, whereby the trained model produced 1000 predictions for each time point by randomly sampling observed meteorological conditions while holding temporal variables constant. The median of these predictions represents the normalised pollutant concentration, effectively removing meteorological influences. Model performance was assessed using the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE), while the distribution of resampled predictions provided uncertainty estimates.

A randomly selected 70% subset of the dataset was used for model training, while the remaining 30% served as a test set for model evaluation. This analytical approach has been successfully implemented in previous studies, including those conducted by Vu, T.V., et al.³² and Ghaffarpasand, et al.³³ All computational analyses were conducted using the ‘rmweather’ package within the R statistical environment.

2.4. Analysis of short- and long-distance emission sources

An analysis of short-distance emission sources was conducted utilising bivariate polar plots (BPPs), a methodology widely used for evaluating the influence of localised emissions on ambient air quality.^34–36 This analytical technique has gained prominence due to its efficacy in source characterisation when high-resolution temporal data for both air quality parameters and meteorological conditions are available. BPPs utilise a radial-axis representation to visualise statistical metrics of atmospheric pollutant concentrations as a function of wind speed and direction. This approach generates a polar coordinate system, whereby the radial component corresponds to wind speed and the angular component represents wind direction. The resulting graphical output displays the calculated statistical values, such as mean pollutant concentrations, as colour-coded cells within this polar framework.

The ‘percentileRose’ function was also employed to supplement BPPs, facilitating a more comprehensive examination of short-distance emission sources and an investigation into the influence of seasonality. It plots Conditional Probability Functions (CPF), which suggest the probability of elevated pollutant concentrations occurring from specific wind directions.³⁷ The CPF is calculated as the ratio of samples exceeding a specified high concentration threshold to the total number of samples within a given wind sector. The function enables straightforward comparison of multiple pollutants on a single plot, utilising a 0–1 probability scale.

The application of BPPs and ‘percentileRose’ offers significant advantages in the assessment of the spatial contribution of various emission sources to fluctuations in pollutant levels. The integration of meteorological data with air quality measurements allows us to identify potential source regions and their relative impacts on local air quality, thereby enhancing our understanding of pollutant transport and dispersion mechanisms in urban environments.

The investigation of long-range emission sources was conducted through analysis of air mass trajectories traversing the study area. This approach utilised the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model to compute 72 hours backward trajectories at hourly resolution.^38–40 The trajectory calculations were initiated at an elevation of 100 metres above ground level, a carefully selected altitude designed to ensure the receptor's position within the atmospheric boundary layer throughout the study period, while mitigating immediate surface influences. Meteorological parameters essential for the analysis, including ambient temperature, relative humidity, atmospheric pressure, wind speed, and wind direction, were acquired at hourly intervals from the nearest observational site within the National Oceanic and Atmospheric Administration (NOAA) Integrated Surface Database (ISD). These data were accessed and processed using the ‘WorldMet’ R package, thereby ensuring the generation of a robust and comprehensive meteorological dataset for the study.

To characterise the predominant air mass pathways, a cluster analysis was implemented. The process began with the determination of the optimal number of clusters using the gap statistic method.⁴¹ This technique compares the observed within-cluster dispersion to that of randomly generated data, effectively identifying the number of clusters that maximises data separation. Following this optimisation step, the K-means algorithm, a robust unsupervised machine learning approach, was applied to categorise air mass back-trajectories into distinct groups based on their spatial characteristics.⁴² The K-means method iteratively partitions n observations into k clusters, assigning each data point to the cluster with the nearest centroid. This process continues until convergence, minimising the within-cluster sum of squares and optimising cluster cohesion.⁴³

To further refine the analysis, the elbow method was employed to examine the explained variance as a function of cluster numbers, complemented by silhouette analysis, which assesses clustering quality based on the similarity of observations within their assigned clusters relative to other clusters.⁴⁴ This multi-faceted approach to cluster determination and validation ensures a rigorous and comprehensive categorisation of air mass trajectories, providing a solid foundation for subsequent analyses of pollutant transport and source attribution.

3. Results and discussions

3.1 Air quality analysis

Fig. 3 shows the temporal evolution of the monthly-averaged concentrations of monitored air pollutants in Kičevo over the period 2019–2023. The solid lines represent smoothed trends derived using the LOESS (locally estimated scatterplot smoothing) method to illustrate underlying seasonal and long-term variations. LOESS provides a visual representation of underlying seasonal and long-term variations in the data but does not directly test statistical significance. The data reveal notable fluctuations in CO and PM_2.5 levels, while a declining trend is observed for NO₂, O₃, and, to some extent, PM₁₀. In contrast, SO₂ concentrations exhibit a marked increase, particularly during the 2021–2022 interval, suggesting a shift in emission patterns or contributing sources during that period.

Fig. 3 Temporal variation of monthly-averaged concentrations of (a) CO, (b) NO₂, (c) O₃, (d) SO₂, (e) PM₁₀, and (f) PM_2.5 in Kičevo, North Macedonia, from 2019 to 2023. The solid lines represent smoothed trends derived using locally estimated scatterplot smoothing (LOESS), which highlights underlying seasonal and long-term patterns by attenuating short-term fluctuations. The smoothing serves as a visual aid and does not indicate statistical significance.

The observed increase in SO₂ concentrations in early 2021 in Kičevo likely reflects changes in regional energy and industrial activities. Coal-fired power plants, including TE Oslomej and TE Bitola, are major contributors to SO₂ emissions in North Macedonia, and periods of increased coal combustion or changes in fuel quality can cause temporary spikes in ambient SO₂ levels. Such variations may also be influenced by broader regional emission patterns, consistent with trends reported for North Macedonia during this period.^45,46 This suggests that the notable rise in SO₂ AQI during early 2021 reflects a combination of local industrial emissions and regional energy practices, although further detailed local studies would be required to precisely identify the specific causes.

Fig. 4 illustrates the temporal evolution of the Air Quality Index (AQI) for each studied pollutant in Kičevo between 2019 and 2023. Unlike the conventional aggregated AQI, which combines several pollutants into a single composite value, we calculated individual pollutant AQIs to enable a clearer examination of pollutant-specific trends and potential health implications. Each pollutant's measured concentration was converted to its corresponding AQI value using the standard methodology defined by the United States Environmental Protection Agency.⁴⁷ This pollutant-by-pollutant approach provides a more transparent representation of the underlying variability and avoids potential bias arising from incomplete datasets that could distort a unified AQI. Consequently, it allows a more robust interpretation of temporal dynamics and exceedance patterns for each pollutant.

Fig. 4 Temporal variation and trend of the Air Quality Index (AQI) of the studied pollutants in Kičevo, North Macedonia, from 2019 to 2023. Data on PM_2.5 was only available from 2021 onward. Red dashed lines represent linear regression trends fitted to the monthly mean AQI values, and grey shaded areas indicate the 95% confidence intervals.

A detailed examination of the AQI time series reveals contrasting temporal dynamics across monitored pollutants. The red dashed lines represent linear regression trends fitted to the monthly mean AQI values, with grey shaded areas indicating the corresponding 95% confidence intervals. The trends reveal statistically significant decreases in NO₂ (slope = −0.92, p < 0.001, R² = 0.35) and O₃ (slope = −3.18, p < 0.001, R² = 0.32), indicating a statistically significant downward trend over the study period. In contrast, CO shows a relatively stable trend with a slight decrease, though the trend is not statistically significant (slope = −0.39, p = 0.043, R² = 0.07). The data show a notable and statistically significant increase in SO₂ AQI values (slope = 0.34, p < 0.001, R² = 0.76), particularly from February through September 2022, a period that coincided with shifts in regional energy practices and industrial activities that may have influenced local air quality. PM₁₀ (slope = −0.80, p = 0.11, R² = 0.04) exhibits similar variability but a less pronounced, non-significant reduction. PM_2.5 shows considerable short-term variability superimposed on an overall decreasing trend, though this trend is not statistically significant (slope = −5.29, p = 0.14, R² = 0.06). Despite the lack of statistical significance in its temporal trend, the AQI analysis clearly indicates that PM_2.5 remains the dominant pollutant with respect to health impacts.

Fig. 5 presents the temporal variations of the studied pollutants, analysed on both monthly and diurnal scales. To facilitate the interpretation of fluctuations in air pollutant levels, the study period was stratified into three phases relative to the COVID-19 pandemic: the pre-pandemic period (January 2019 -January 2020), the pandemic period (February 2020–April 2022), and the post-pandemic period (May 2022–December 2023). Further details regarding this stratification are provided in Table S1 within the SI. In North Macedonia, restrictive measures such as curfews, travel bans, and business closures were first implemented in March 2020 and gradually lifted by early 2022. The post-pandemic period, however, coincides with the onset of the Russia–Ukraine conflict, which led to significant changes in energy consumption and mobility across the region. This overlap is acknowledged as a confounding factor that may influence air quality recovery trends during that time. It should be noted that PM_2.5 data were available only from 2021, covering the latter part of the pandemic and the post-pandemic period.

Fig. 5 Time trend analysis of (a) CO, (b) NO₂, (c) O₃, (d) SO₂, (e) PM₁₀, and (f) PM_2.5 pollutants for Kičevo 2019–2023. The left- and right-hand side figures represent the monthly and daily variations, respectively. Data on PM_2.5 was only available from 2021 onward.

A notable decline in NO₂ and O₃ concentrations is evident in the post-pandemic period compared to pre-pandemic baseline levels. During the pandemic, NO₂ concentrations decreased substantially, particularly during lockdown phases, while CO and PM₁₀ levels remained relatively stable. In contrast, SO₂ concentrations exhibited a marked increase during both the pandemic and post-pandemic periods, with post-pandemic levels exceeding pre-pandemic concentrations by more than a factor of two. NO₂ concentrations during the post-pandemic period remained below pre-pandemic levels, despite a partial recovery from pandemic minima. Changes in the seasonal cycle of NO₂ during this period likely reflect a combination of factors, including altered commuting patterns, continued teleworking, changes in residential heating behaviour, and broader demographic trends such as population decline due to emigration. These combined influences suggest that post-pandemic NO₂ variability was driven by more than meteorological factors alone, consistent with recent European studies indicating persistent mobility-related emission changes following the pandemic, e.g. 48.

Pollutants primarily associated with wood-burning and coal-based heating systems, including CO, PM₁₀, and PM_2.5, demonstrate elevated concentration levels during the winter months. The temporal distribution of PM_2.5 concentrations in Kičevo exhibits a distinct U-shaped pattern when analysed monthly, indicating pronounced seasonal fluctuations in particulate matter pollution. This phenomenon can be attributed to a confluence of meteorological, anthropogenic, and biogenic factors. During the winter months, PM_2.5 levels reach their zenith due to the intensification of residential heating activities, the prevalence of temperature inversions that impede vertical mixing and trap pollutants in the lower atmosphere, and elevated vehicle emissions associated with cold-start conditions. Conversely, the summer months are characterised by a marked reduction in PM_2.5 concentrations, which can be ascribed to the cessation of heating-related emissions and improved atmospheric dispersion conditions (see Fig. 2(b)). Ozone (O₃) concentrations exhibit a distinctive seasonal pattern, with peak levels occurring during the summer months. This phenomenon can be attributed to an increase in photochemical reactions, which are driven by elevated levels of solar radiation. Increases in photochemical reactions under higher solar radiation have been documented in regional studies. For example, Oikonomakis⁴⁹ found that solar “brightening” in Europe between 1990 and 2010 led to increased photolysis rates and consequently higher surface ozone formation. Alexandri et al.⁵⁰ report that over Greece, tropospheric NO₂ and aerosol concentrations strongly modulate the amount of solar radiation reaching the surface, which in turn influences photochemistry. These findings suggest that in the Western Balkans, elevated solar radiation during clear-sky periods likely enhanced NO₂ photolysis and other photochemical processes, contributing to the seasonal shifts in NO₂ and ozone cycles observed post-pandemic.

3.1.1 Deweathered data analysis. The exclusion of meteorological influences from the air pollution datasets enables a more robust isolation of non-meteorological drivers of air quality variability. Fig. 6 illustrates the temporal evolution of deweathered pollutant concentrations over the study period, obtained after accounting for key meteorological parameters (wind speed, wind direction, and temperature) and temporal factors (weekday/weekend effects and seasonal cycles). The effectiveness of the normalisation is quantitatively supported by the predictive performance of the deweathering models, reported in Table S2 in the SI, with moderate explanatory power achieved for gaseous pollutants, indicating that a substantial fraction of observed variability is attributable to systematic meteorological and activity-related influences that have been successfully removed. In contrast, lower performance for particulate matter, reflects the more complex interplay of local resuspension, secondary formation, and regional transport processes that are not fully captured by the available predictors.

Fig. 6 Deweathering time series plots of (a) CO, (b) NO₂, (c) O₃, (d) SO₂, (e) PM₁₀, and (f) PM_2.5 pollutant for the city of Kičevo (2019–2023). PM_2.5 data was only available for after 2021.

Results reveal divergent trends across various atmospheric constituents during the pandemic (February 2020–April 2022). Looking at the general trends, CO and SO₂ exhibit negligible fluctuations over the pandemic, while NO₂, O₃, and PM₁₀ show substantial declines in concentrations. These reductions are consistent with decreased anthropogenic activity during lockdown periods, particularly reduced traffic-related emissions, as also indicated by deweathered NO₂ trends. However, changes in ozone concentrations likely reflect a complex interaction between precursor availability, atmospheric chemistry, and meteorological conditions, rather than a direct response to emission reductions alone.⁵¹

Fig. 6 illustrates notable anomalies (spikes) in pollutant concentrations, particularly evident around 2020 and 2022, with pronounced elevations in CO, NO₂, SO₂, PM₁₀ and PM_2.5 levels and a significant reduction in O₃. To investigate the source of these anomalies, Fig. 7 presents a temporal analysis of electricity generation patterns alongside coal and fuel oil (Mazut) consumption at the TEC Oslomej power plant, which is situated in proximity to the city (see Fig. 1), spanning from 2018 to 2023. The analysis draws upon primary data obtained from official power plant reports, provided directly by the facility's director. These reports, originally documented in Macedonian and available in printed format, contain comprehensive records of electricity production and coal consumption metrics. The raw data underwent systematic review and was subsequently transformed into a structured dataset aligned with the research objectives.

Fig. 7 Monthly trends in (a) electricity generation, (b) coal consumption, and (c) fuel oil (Mazut) consumption at the TE Oslomej power plant from March 2018 through March 2024. Data points indicate observed values, while the solid line represents the trend with a 70% confidence interval shaded around it. The data used to make Fig. 7 are tabulated in Table S3 in the SI.

For the operation of TPP Oslomej, coal is primarily supplied from the active surface mine “Oslomej–West”, which is situated directly next to the power plant. In contrast, the so-called “Old Mine,” which is located within the same mining basin but lies slightly farther from the power plant and which also forms part of the former “Oslomej–East” deposit, is no longer an active extraction site. Historically, this mine was one of the original lignite sources for the plant, but in recent years, including the period analysed here, no coal has been excavated from it. Its depletion and closure underscore the reliance of the power plant on the remaining Oslomej–West reserves and highlight the transition challenges faced by the region's energy sector.

The elevated pollutant concentrations observed during 2020 coincided with unusually cold winter conditions, as indicated by exceptional monthly mean low boundary layer heights, below ∼300–400 m, see Fig. 2(b), compared to typical summer values exceeding ∼600 m, which favoured pollutant accumulation through reduced vertical mixing. COVID-19 containment measures increased time spent indoors, likely amplifying residential energy demand. Fig. 7 shows concurrent increases in electricity generation and coal and fuel oil consumption during winter 2020 relative to the preceding year, temporally aligned with the strongest deweathered concentration anomalies observed across primary pollutants. For example, accomplished electricity generation in January 2020 (38 894 MWh) was nearly an order of magnitude higher than in January 2019 (4000 MWh), accompanied by elevated coal and fuel oil consumption (Fig. 7). Taken together, these observations suggest that the combination of adverse meteorological conditions and increased energy demand likely contributed to higher pollutant concentrations, although the relative contribution of each factor cannot be fully disentangled with the available data.

A second anomaly in 2022 coincides with the reactivation of the Ptolemaida coal-fired power plant in northern Greece, implemented to secure domestic energy supply amid regional fuel shortages and elevated natural gas prices. This period was characterised by increased coal-based electricity generation and associated SO₂ emissions. Given the well-established long-range transport potential of SO₂, elevated concentrations observed in Kičevo are consistent with regional transport influences. Supporting this interpretation, HYSPLIT back-trajectory analysis (Fig. 10) indicates a dominant southern air-mass origin (68.5%), encompassing northern Greece.

Pollution spikes in 2020 and 2022 differ significantly in their causes and duration. The 2020 winter crisis, marked by pandemic-induced supply chain disruptions, increased residential heating demand, and reliance on fossil fuels to secure energy supply, was treated as a temporary emergency, leading to short-term increases in power generation, as shown in Fig. 7(a). Plant operators expected the high energy demand to be brief and planned accordingly. A key distinction between the 2020 and 2022 anomalies lies in their underlying drivers. The 2020 spike was a short-term adjustment in energy production to address shifts in demand associated with the public health emergency and was treated as a temporary deviation. In contrast, the 2022 event was rooted in long-term energy planning and crisis response, with implications for sustained emission levels detectable across borders. The latter reflects the shifting energy landscape in Southeastern Europe, where coal remains a fallback option during supply crises, with measurable consequences for regional air quality.

Both anomalous periods were associated with increased fossil-fuel combustion for energy production, which is known to emit a complex mixture of gaseous pollutants and particulate matter that can influence regional atmospheric chemistry. However, the specific role of VOC emissions and their effects on ozone formation could not be quantified in this study due to the lack of direct measurements. During the 2022 period, the regional energy crisis was accompanied by a documented increase in fuel prices,⁵⁶ which may have contributed to reduced mobility and lower NO₂ concentrations, although the magnitude of this effect cannot be directly assessed with the available data.

The post-pandemic recovery period revealed distinct temporal patterns in pollutant concentrations. CO concentrations exhibited an immediate upward trend, consistent with the rapid resumption of vehicular traffic and transportation activities. In contrast, NO₂ and SO₂ demonstrated a more gradual, synchronous increase, suggesting a delayed recovery in industrial operations relative to transportation sectors. Throughout this recovery phase, ozone concentrations maintained relative stability, indicating that the proportional relationship between NO₂ emissions and other ozone precursor pollutants remained largely unchanged during this period.

3.2. Analyses of short- and long-distance emission sources

3.2.1. Analyses of short-distance emission sources. Fig. 8 depicts the temporal and directional variations of six major air pollutants (CO, NO₂, O₃, SO₂, PM₁₀, and PM_2.5) in Kičevo across three distinct time periods: pre-pandemic (January 2019–January 2020) pandemic (February 2020–April 2022), and post-pandemic (May 2022–December 2023).

Fig. 8 Bivariate polar plots of (a) CO, (b) NO₂, (c) O₃, (d) SO₂, and (e) PM₁₀ pollutant for the city of Kičevo (2019–2023). Left, middle, and right polar plots correspond with pre-pandemic, pandemic and post-pandemic period, respectively.

For the gaseous pollutants, CO and NO₂ exhibit broadly similar directional behaviours. During the pandemic, both pollutants show higher concentrations under easterly and southeasterly winds, suggesting the persistence of local emission influences even during restricted mobility periods. In the post-pandemic period, their concentrations slightly declined relative to the pandemic phase, though values remained above pre-pandemic levels under some wind sectors, reflecting gradual recovery of traffic and heating activity.

O₃ exhibits relatively higher concentrations under southwestern winds during the pre-pandemic period, with less pronounced directional contrasts during the pandemic and post-pandemic phases, suggesting a stronger influence of regional background air masses rather than local emission changes. SO₂ pollution roses indicate elevated concentrations under northerly winds before the pandemic, consistent with emissions from nearby industrial and residential sources. During the pandemic, SO₂ levels under northern winds remained moderate, showing only a limited reduction. In the post-pandemic period, concentrations increased again, with new peaks emerging under southern and southwestern winds, suggesting both renewed operation of local sources and episodic influence from long-range transport.

Prior to the pandemic, PM₁₀ concentrations were elevated under southeastern winds, with a reduction during the pandemic. After the pandemic, southerly winds became more influential, associated with moderately higher PM₁₀ levels, potentially exceeding pre-pandemic values. For PM_2.5, the highest concentrations are consistently observed under southwestern winds, a pattern that persists across all periods and indicates stable influence from urban and transport-related emissions. Overall, particulate matter shows greater directional stability compared to gaseous pollutants, though with modest post-pandemic increases in magnitude.

The ‘percentileRose’ plots, presented in Fig. S2 in the SI, further reinforce the spatial and seasonal patterns observed in the main analysis. CO (Fig. S2(a)) shows a pronounced winter peak, particularly from easterly and south-easterly directions, reflecting intensified combustion activities and stagnant meteorological conditions. A similar seasonal peak is observed for NO₂, consistently elevated from the eastern sectors, indicating persistent emission sources. O₃, as expected, displays a contrasting pattern, with summer maxima from the north and northeast, consistent with enhanced photochemical activity. For SO₂, elevated wintertime concentrations are notable from the northwestern and northeastern directions, again suggesting high energy demand and local activity of power plants near the city. PM₁₀ and PM_2.5 levels follow the general winter enhancement trend, with directionality broadly aligned to that seen in the pollution rose plots. These findings support the existence of both fixed and seasonally variable short-range emission sources around Kičevo.

The hourly variation in NO₂ concentration levels, as depicted in Fig. 9, provides a comprehensive visual representation for analysing the interplay between short-distance emission sources, meteorological phenomena, and socio-economic variables. The heat map incorporates data from 2020 onwards due to substantial data gaps in earlier periods. Three notable observations emerge.

Fig. 9 Heatmap of the hourly variation of NO₂ concentration levels (µg m⁻³) over the years 2020–2023, highlighting diurnal and seasonal trends.

Firstly, elevated NO₂ concentrations are observed in February, March, and April compared to other months. This pattern can be attributed to the transition period between winter and spring, when boundary-layer heights (discussed in Fig. 2) remain relatively low and temperatures are still cool, limiting vertical dispersion and slowing photochemical removal of NO₂ as stated from Okionomakis et al.⁴⁹ Although boundary-layer heights in late spring and summer are sometimes lower on average, higher solar radiation and stronger photochemical activity during those months enhance NO₂ oxidation and lead to lower ambient concentrations, which explains why May–July values are reduced. Increased residential heating and traffic activity in late winter and early spring may further amplify emissions during this period.

Secondly, distinct nocturnal NO₂ peaks are observed in October, November, and December, though in 2021 some peaks occurred during summer months. These nighttime increases may be associated with periods of reduced atmospheric mixing under low BLH or stagnant wind conditions rather than specific economic activities. A general temporal trend is evident, with the intensity of these peaks decreasing from 2020 to 2023, possibly reflecting broader socio-economic or behavioural changes affecting emissions.

Overall, the observed temporal patterns of NO₂ likely arise from a combination of meteorological conditions (boundary-layer height, temperature, wind speed, and possible inversions) and anthropogenic factors (heating, traffic, and local industrial emissions). The relative contributions of these factors remain uncertain but collectively explain the persistence of elevated NO₂ in late winter and early spring.

Finally, while elevated concentrations of NO₂ have been a recurring phenomenon, the summer of 2023 stands out with particularly pronounced peaks. This can be attributed to the surge in tourism and the influx of returning emigrants during the holiday season. Although similar seasonal patterns were observed in previous years, the post-pandemic recovery appears to have amplified this effect, making 2023 a period of heightened environmental concern. Statistical data support this interpretation, showing a 30.1% overall increase in the number of tourists from January to June 2023, including a striking 46.8% rise in foreign tourists.⁵² This robust resurgence of international travel likely contributed disproportionately to NO₂ emissions, as foreign visitors tend to rely heavily on rental vehicles and other transport services. A comparative analysis of NO₂ levels between 2023 and earlier years confirms that the summer of 2023 represents an anomalous period, underscoring the environmental implications of intensified tourism activity.

While temporal co-variation between meteorological indicators, energy production, and pollutant concentrations is evident, the lack of sector-specific emission inventories and activity data prevents formal attribution or quantification of individual contributions.

3.2.2. Analysis of long-distance emission sources

To assess the influence of long-distance pollutant transport on Kičevo's air quality, a clustering analysis of air mass back-trajectories was performed (see Section 2.4 for methodological details). The optimal number of clusters was determined using the gap statistic method,⁴¹ which compares the observed within-cluster dispersion to that of randomly generated datasets. As shown in Fig. S3 in the SI, the gap values plateau after reaching a maximum, indicating that three clusters best capture the structure of the data.

The spatial distribution and directional pathways of these three dominant air mass clusters are illustrated in Fig. 10, while their relative yearly contributions across the study period are summarised in Table 1. Importantly, the proportions of these clusters remained relatively stable throughout the years, suggesting persistent atmospheric circulation patterns affecting the city.

Fig. 10 Cluster analysis of the averaged 72 h air mass back-trajectories arriving at the city of Kičevo, North Macedonia, for the years 2019–2023. The numbered red arrows show the locations of 1. Oslomej, 2. TEC Bitola and 3. Ptolemais power plants.

Table 1 Yearly relative contributions (%) of the three dominant air mass clusters identified over the city of Kičevo, North Macedonia, for the years 2019–2022

Year	Cluster 1	Cluster 2	Cluster 3
2019	9.4	69.3	21.3
2020	14.9	68.1	17.0
2021	12.5	68.3	19.2
2022	11.5	67.8	20.7

---

The most dominant cluster, representing 68.5% of all trajectories, originates from the south, passing over the North Pindus and Tomorr National Parks and the broader region around Thessaloniki, Greece. Along this pathway lie key emission sources, including the Ptolemaida lignite complex in northern Greece and TEC Bitola coal-fired power plants in southern North Macedonia (locations marked in Fig. 10). Emissions from these facilities have been shown to contribute substantially to SO₂ and particulate matter pollution, with long-range transport influencing downwind urban areas⁵⁷ and.⁴⁶ The elevated SO₂ concentrations observed in Kičevo (Fig. 8(d)) can therefore be, at least in part, attributed to transboundary transport from these coal-fired power plants).

Another significant cluster arrives from the northeast, carrying air masses from the Skopje metropolitan area. Skopje is a major urban and industrial hub with dense traffic, metal processing industries, chemical plants, and residential heating. These sources emit a complex mixture of pollutants, notably SO₂, NO₂, and PM₁₀, which are likely transported towards Kičevo along this path. The remaining 11.2% of trajectories originate from the west, crossing Albania and particularly the region around Tirana, the capital. This cluster, although less frequent, may still carry urban and industrial emissions from Tirana, potentially impacting air quality in Kičevo.

4. Conclusion

Urban air quality is influenced by a wide range of complex factors, many of which extend beyond local emissions. Socio-economic conditions, long-range pollutant transport, and unexpected events all play significant roles in shaping air quality patterns. Against this broader context, this study presents a comprehensive analysis of air quality trends in Kičevo, North Macedonia, from 2019 to 2023. The results reveal that air pollution in this small Western Balkans city is shaped by a complex interplay of meteorological conditions, energy production patterns, residential heating practices, transportation behaviour, and demographic shifts, particularly ongoing emigration. Significant temporal and spatial variability in pollutant concentrations was observed. Winter months consistently exhibited elevated levels of PM_2.5, PM₁₀, NO₂, SO₂, and CO, driven by residential heating, combustion of solid fuels, and unfavourable atmospheric conditions such as low boundary layer height and temperature inversions. In contrast, ozone concentrations peaked in the summer, influenced by photochemical activity and precursor availability.

A key finding of this study is the apparent impact of demographic change on local emissions. The population decline in Kičevo, driven by emigration and ageing, has likely contributed to a reduction in heating and some transport-related emissions, particularly on a per capita basis, with the most notable changes occurring in the post-pandemic period. However, the number of registered diesel cars increased significantly between 2020 and 2022 (Fig. S1), suggesting a structural change in vehicle ownership. This may reflect a growing preference for second-hand diesel vehicles, longer use of older cars, or shifts in commuting patterns, which could offset some reductions in transport emissions. Consequently, while population decline influences total emissions, changes in the vehicle fleet composition and usage patterns add complexity to the net effect on local air quality. Directional and seasonal analyses further identified the south and southeast quadrants, including emissions from nearby thermal power plants and long-range transport corridors, as dominant sources of pollutant inflows.

The period also featured two notable air pollution anomalies linked to domestic energy production patterns. In 2020, short-term increases in power generation during the COVID-19 lockdowns led to temporary pollution spikes. In 2022, longer-term shifts toward domestic coal and fuel oil usage, driven by market instability, resulted in sustained increases in emissions, particularly from TE Oslomej. Despite this, rising fuel prices during the same year appear to have reduced traffic volumes, potentially offsetting some NO₂ emissions. Weekly and seasonal patterns highlight the dual importance of traffic and household energy use as determinants of local air quality. These findings emphasise the tight coupling between energy systems and air quality. Periods of instability in energy supply, whether driven by public health crises, economic shocks, or geopolitical tensions, were consistently mirrored in local and regional pollution levels. Reliance on coal and fuel oil during times of energy insecurity not only slowed decarbonisation but also increased exposure to harmful pollutants with measurable cross-border effects. From a scientific perspective, this suggest that energy system vulnerability is a key driver of air quality degradation. Building resilient, diversified, and cleaner energy infrastructures is therefore critical not only for climate mitigation but also as a strategy to safeguard public health.

Future research should investigate the health impacts of long-term exposure, particularly in relation to mortality and respiratory illness rates in the Polog and southwestern regions (Kičevo, Gostivar, Tetovo). More broadly, these results underscore that scientific understanding of the interactions between energy systems, socio-economic dynamics, and meteorology is essential for guiding both urban planning and energy transitions in the Western Balkans and similar regions.

Author contributions

Albinota Nuredini: conceptualization, data curation, methodology, original draft preparation, visualization, writing–reviewing and editing, Omid Ghaffarpasand: conceptualization, methodology, supervision, Writing–reviewing and editing, Dimitrios Bousiotis: conceptualization, methodology, supervision, writing–reviewing and editing Francis Pope: conceptualization, supervision, methodology, writing–reviewing and editing, funding acquisition.

Conflicts of interest

No conflicts of interest.

Data availability

The air quality data, upon which the analysis is based upon, is publicly accessible online via the air quality portal, published by the Ministry of Environment and Physical Planning of North Macedonia.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5ea00116a.

Acknowledgements

Albinota Nuredini thanks the McCall McBain Foundation and the Clean Air Fund for the McCall McBain Clean Air Fellowship which funded them to study for their Masters in ‘Air Pollution Management and Control’. This research has been supported by the UKRI ISPF institutional support grant (ODA) funding.

References

1. K. Liu, et al., Air pollution and individual risk preference: Evidence from China, Energy Econ., 2024, 136, 107738.
2. C. Tao, et al., Estimating neighborhood-based mortality risk associated with air pollution: A prospective study, J. Hazard. Mater., 2024, 475, 134861.
3. K. Liu, et al., Ambient air pollution and Children's health: An umbrella review, Atmos. Pollut. Res., 2024, 15(6), 102108.
4. F. A. Rahnemaei, et al., Prenatal exposure to ambient air pollution and risk of fetal overgrowth: Systematic review of cohort studies, Ecotoxicol. Environ. Saf., 2024, 280, 116526.
5. R. Liang, et al., Air pollution exposure, accelerated biological aging, and increased thyroid dysfunction risk: Evidence from a nationwide prospective study, Environ. Int., 2024, 188, 108773.
6. WHO, World Health Organization report on the global status for road safety 2023. World Health Organization, 2023.
7. K. Arshad, N. Hussain, M. H. Ashraf and M. Z. Saleem, Air pollution and climate change as grand challenges to sustainability, Sci. Total Environ., 2024, 928, 172370.
8. A. Talaiekhozani, et al., Evaluation of emission inventory of air pollutants from railroad and air transportation in Isfahan metropolitan in 2016, J. Air Pollut. Health, 2017, 2(1), 1–18.
9. B. Liu, et al., Assessing the Impacts of Birmingham's Clean Air Zone on Air Quality: Estimates from a Machine Learning and Synthetic Control Approach, Environ. Resour. Econ., 2023, 86(1), 203–231.
10. O. Ghaffarpasand, et al., Real-world evaluation of driving behaviour and emission performance of motorcycle transportation in developing countries: A case study of Isfahan, Iran, Urban Clim., 2021, 39, 100923.
11. J. Petrović and J. Ateljević, Neo-Liberalism, Depopulation and Economic Stagnation in the Balkans, J. Balkan Near East. Stud., 2024, 26(4), 411–431.
12. A. Heil, Depopulation Disaster: The Balkans And Its Creeping Demographic Crisis. RadioLiberty, 2020, Available from: https://www.rferl.org/a/depopulation-disaster-the-balkans-and-its-creeping-demographic-crisis/30477952.html.
13. T. Judah, B. Bye-bye: A region in critical demographic decline. Balkan insight, 2019. 10: p. 14.
14. C. A. Belis, et al., Urban pollution in the Danube and Western Balkans regions: The impact of major PM2.5 sources, Environ. Int., 2019, 133, 105158.
15. C. A. Belis, et al., Assessment of health impacts and costs attributable to air pollution in urban areas using two different approaches. A case study in the Western Balkans, Environ. Int., 2023, 182, 108347.
16. European Bank for Reconstruction and Development, Regional Gasification Project. 2022.
17. C. Samara, et al., Chemical characterization and receptor modeling of PM10 in the surroundings of the opencast lignite mines of Western Macedonia, Greece, Environ. Sci. Pollut. Res., 2018, 25(13), 12206–12221.
18. K. Rizos, et al., A machine learning modelling approach to characterize the background pollution in the Western Macedonia region in northwest Greece, Atmos. Pollut. Res., 2023, 14(10), 101877.
19. M. Srbinovska, et al., Quantifying the impact of meteorological factors and green infrastructure location on particulate matter (PM) mitigation in Republic of North Macedonia using sensor collected data, Meas. Sens., 2023, 27, 100819.
20. State Statistical Office of the Republic of North Macedonia, Census of Population, Households and Dwellings in the Republic of North Macedonia, 2021 - First Dataset, State Statistical Office of the Republic of North Macedonia: Skopje, North Macedonia, 2022.
21. State Statistical Office of the Republic of Macedonia, Census of Population, Households and Dwellings in the Republic of Macedonia, 2002 - Final Data. State Statistical Office, Dame Gruev - 4, Skopje: Skopje, Macedonia, 2005.
22. European Commission, Communication from the Commission to the European Parliament, the Council, the European Economic and Social Committee and the Committee of the Regions in 2021 Communication on EU Enlargement Policy, C.S.W. DOCUMENT and N.M. Report, Editors. 2021: Strasbourg.
23. D. C. Carslaw and K. Ropkins, Openair—an R package for air quality data analysis, Environ. Model. Softw., 2012, 27, 52–61.
24. Copernicus, ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus climate change Service climate data store (CDS), 2017. 15(2): p. 2020.
25. H. Hersbach, et al., The ERA5 global reanalysis, Q. J. R. Meteorol. Soc., 2020, 146(730), 1999–2049.
26. J. Guo, et al., A merged continental planetary boundary layer height dataset based on high-resolution radiosonde measurements, ERA5 reanalysis, and GLDAS, Earth Syst. Sci. Data, 2024, 16(1), 1–14.
27. D. J. Seidel, et al., Climatology of the planetary boundary layer over the continental United States and Europe, J. Geophys. Res. Atmospheres, 2012, 117, D17106.
28. S. Grange, et al., Random forest meteorological normalisation models for Swiss PM10 trend analysis, Atmos. Chem. Phys., 2018, 18, 6223–6239.
29. S. K. Grange and D. C. Carslaw, Using meteorological normalisation to detect interventions in air quality time series, Sci. Total Environ., 2019, 653, 578–588.
30. J. Friedman, R. Tibshirani and T. Hastie, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Springer-Verlag, New York New York, 2009.
31. S. K. Grange, et al., Random forest meteorological normalisation models for Swiss PM10 trend analysis, Atmos. Chem. Phys., 2018, 18(9), 6223–6239.
32. T. V. Vu, et al., Assessing the impact of clean air action on air quality trends in Beijing using a machine learning technique, Atmos. Chem. Phys., 2019, 19(17), 11303–11314.
33. O. Ghaffarpasand, et al., The impact of urban mobility on air pollution in Kampala, an exemplar sub-Saharan African city, Atmos. Pollut. Res., 2024, 15(4), 102057.
34. S. K. Grange, A. C. Lewis and D. C. Carslaw, Source apportionment advances using polar plots of bivariate correlation and regression statistics, Atmos. Environ., 2016, 145, 128–134.
35. E. Jang, et al., Spatial and temporal variation of urban air pollutants and their concentrations in relation to meteorological conditions at four sites in Busan, South Korea, Atmos. Pollut. Res., 2017, 8(1), 89–100.
36. S. Sooktawee, et al., Characterising particulate matter source contributions in the pollution control zone of mining and related industries using bivariate statistical techniques, Sci. Rep., 2020, 10(1), 21372.
37. L. L. Ashbaugh, W. C. Malm and W. Z. Sadeh, A residence time probability analysis of sulfur concentrations at grand Canyon National Park, Atmos. Environ., 1985, 19(8), 1263–1270.
38. R. Draxler, HYSPLIT4 user's guide. NOAA Tech. Memo. ERL ARL-230: p. 1–35.
39. R. Draxler, HYSPLIT4 user's guide NOAA Tech. Memo. ERL ARL-230 (Silver Spring, MD: NOAA Air Resources Laboratory), 1999.
40. A. F. Stein, et al., NOAA's HYSPLIT atmospheric transport and dispersion modeling system, Bull. Am. Meteorol. Soc., 2015, 96(12), 2059–2077.
41. R. Tibshirani, G. Walther and T. Hastie, Estimating the number of clusters in a data set via the gap statistic, J. R. Stat. Soc. Ser. B Stat, 2001, 63(2), 411–423.
42. J. A. Hartigan and M. A. Wong, Algorithm AS 136: A k-means clustering algorithm, J R Stat Soc Ser C Appl Stat, 1979, 28(1), 100–108.
43. S. Dorling, T. Davies and C. Pierce, Cluster analysis: A technique for estimating the synoptic meteorological controls on air and precipitation chemistry—Method and applications. Atmospheric Environment. Part A, General Topics, 1992, 26(14), 2575–2581.
44. P. J. Rousseeuw, Silhouettes: a graphical aid to the interpretation and validation of cluster analysis, J. Comput. Appl. Math., 1987, 20, 53–65.
45. EEA, North Macedonia– air pollution country fact sheet 2024. European Environment Agency, 2024.
46. Bankwatch, New report – deadly legal breaches by Western Balkan coal plants increased in 2022. CEE Bankwatch Network, 2023.
47. EPA, U., Technical Assistance Document for the Reporting of Daily Air Quality - the Air Quality Index (AQI). United States Environmental Protection Agency, 2016. EPA-454/B-16-002.
48. G. Varga, et al., Non-uniform tropospheric NO2 level changes in European Union caused by governmental COVID-19 restrictions and geography, City Environ. Interac., 2024, 22, 100145.
49. E. Oikonomakis, et al., Solar “brightening” impact on summer surface ozone between 1990 and 2010 in Europe – a model sensitivity study of the influence of the aerosol–radiation interactions, Atmos. Chem. Phys., 2018, 18(13), 9741–9765.
50. G. Alexandri, A. K. Georgoulias and D. Balis, Effect of Aerosols, Tropospheric NO2 and Clouds on Surface Solar Radiation over the Eastern Mediterranean (Greece), Remote Sens., 2021, 13(13), 2587.
51. X. B. Li, et al., Vertical changes in volatile organic compounds (VOCs) and impacts on photochemical ozone formation, Atmos. Chem. Phys., 2025, 25(4), 2459–2472.
52. A. Dimoski, Country sees 20.6% more tourists in 2023 compared to 2022: statistics. MNA, Skopje: Media Information Agency (MIA), 2024, Available from: https://mia.mk/index.php/en/story/country-sees-20.6-more-tourists-in-2023-compared-to-2022-statistics.
53. L. Villacura, L. F. Sánchez, F. Catalán, A. R. Toro and M. A. Leiva, An overview of air pollution research in Chile: bibliometric analysis and scoping review, challenger and future directions, Heliyon, 2024, 10(3), e25431,  DOI:10.1016/j.heliyon.2024.e25431.
54. H. Akasha, O. Ghaffarpasand and F. D. Pope, Air pollution and economic growth in Dubai: a fast-growing Middle Eastern city, Atmos. Environ.: X, 2024, 21, 100246,  DOI:10.1016/j.aeaoa.2024.100246.
55. Ministry of Environment and Physical Planning (Republic of North Macedonia). Air quality in the Republic of North Macedonia. Skopje: Ministry of Environment and Physical Planning. https://air.moepp.gov.mk/?page_id=221.
56. Fuelo.net. Fuel prices in North Macedonia on 5 October 2025 [Internet]. https://mk.fuelo.net/prices/date/2025-10-05?lang=de.
57. V. Evagelopoulos, P. Begou and S. Zoras, In-Depth Study of PM2.5 and PM10 Concentrations over a 12-Year Period and their Elemental Composition in the Lignite Center of Western Macedonia, Greece, Atmosphere, 2022, 13(11), 1900,  DOI:10.3390/atmos13111900.

---