Random Forests Assessment of the Role of Atmospheric Circulation in PM₁₀ in an Urban Area with Complex Topography

Abstract

This study presents the assessment of the quantitative influence of atmospheric circulation on the pollutant concentration in the area of Kraków, Southern Poland, for the period 2000–2020. The research has been realized with the application of different statistical parameters, synoptic meteorology tools, the Random Forests machine learning method, and multilinear regression analyses. Another aim of the research was to evaluate the types of atmospheric circulation classification methods used in studies on air pollution dispersion and to assess the possibility of their application in air quality management, including short-term PM₁₀ daily forecasts. During the period analyzed, a significant decreasing trend of pollutants’ concentrations and varying atmospheric circulation conditions was observed. To understand the relation between PM₁₀ concentration and meteorological conditions and their significance, the Random Forests algorithm was applied. Observations from meteorological stations, air quality measurements and ERA-5 reanalysis were used. The meteorological database was used as an input to models that were trained to predict daily PM₁₀ concentration and its day-to-day changes. This study made it possible to distinguish the dominant circulation types with the highest probability of occurrence of poor air quality or a significant improvement in air quality conditions. Apart from the parameters whose significant influence on air quality is well established (air temperature and wind speed at the ground and air temperature gradient), the key factor was also the gradient of relative air humidity and wind shear in the lowest troposphere. Partial dependence calculated with the use of the Random Forests model made it possible to better analyze the impact of individual meteorological parameters on the PM₁₀ daily concentration. The analysis has shown that, for areas with a diversified topography, it is crucial to use the variability of the atmospheric circulation during the day to better forecast air quality.

1. Introduction

The abundant air pollution with particulate matter (PM) is a serious environmental and social problem in many regions all over the world [1,2,3,4]. Exposure to ambient PM concentration with a diameter below 10 μm (PM₁₀) increases the possibility of preterm birth [5], deaths from respiratory disease [6] and also causes lung irritation, cellular damage, coughing asthma, and cardiovascular diseases [7]. High PM concentrations in urbanized areas are the consequence of the interaction of many factors, including anthropogenic and natural sources of air pollution, chemical and physical reactions between primary and secondary pollutants, and dispersion conditions determined by atmospheric circulation types, meteorological conditions, and meso- and microclimatic features of the analyzed area [1,8,9]. Numerous studies confirm that atmospheric circulation is an important factor determining the level of air pollution in the lower troposphere, especially in urbanized and industrial areas, which are characterized by elevated pollution emissions [3,10,11,12,13]. The atmospheric circulation processes contribute not only to the dispersion of pollution but also to its transport over great distances from emission sources [14]. Previous research indicated that the duration of air pollution episodes is mainly influenced by the processes of atmospheric blocking and atmospheric stagnation, which contribute to the accumulation of pollutants near the ground, especially during wintertime periods [3,15,16]. Analyses of future climate change indicate that the occurrence of an increase in air stagnation cases is expected [17,18]. Studies performed for Thessaloniki showed that smog episodes can also occur often under weak flow conditions, with warm air advection, as a consequence of stabilization of the lower troposphere and limited vertical mixing [19]. A similar effect of reducing the available mixing volume caused by warm air advection occurs often in mountain valleys and is linked to foehn occurrence [20,21]. It is worth mentioning that similar local weather conditions can occur under very different rearrangements of the large-scale flow; therefore, there is a need to study the relations between unfavorable local pollutant dispersion conditions and large-scale atmospheric circulation, with the impact of particular environmental features.

The problem of the increased PM₁₀ concentration level in the European Union is common for all the nations which are members; in 2018, the daily PM₁₀ concentration limit (50 μg⋅m⁻³) was exceeded in numerous cities in Poland, Bulgaria, the Czech Republic, Croatia, Hungary, Italy and Slovenia [22]. In Poland, poor air quality is a problem, especially in southern regions, both in urban agglomerations and many small localities, where the highest number of days with exceedance of the daily limit of PM₁₀ concentration on the national scale is noted [23]. The problem of air quality in Southern Poland also concerns the area of Kraków, the second largest city in Poland in terms of the number of inhabitants, where air pollution has been a serious and unsolved environmental and social problem for several decades [24,25,26]. In the Małopolska region, where Kraków is located, the main source of PM₁₀ is the emission from the municipal and housing sector (78.9% of the annual emission), from transportation (5%), and from industry (7.8%). During recent decades, there have been many actions aimed to reduce local emissions of PM₁₀ and SO₂ from different sectors. Those actions include liquidation of solid fuel boilers, thermal modernization of buildings, installation of renewable energy sources, modernization of public transport and heating networks or the expansion of bicycle routes. As a result, the air quality in the city has gradually improved, although the PM₁₀ daily limits are still exceeded during the cold seasons [23]. In addition, on 1 September 2019, the prohibition of solid fuels usage in individual heating devices in Kraków was introduced, which could partially contribute to the reduction of PM concentration level during the cold season. The atmospheric circulation conditions play a crucial role in determining the air quality in the city, as it is located in the Wisła River valley, in an area of very diversified relief. The properties of planetary boundary layer (PBL) are strongly modified both by the relief and the synoptic situation, and so are the air pollution’s dispersion conditions which in turn affect the concentration of pollutants. Studies of fog occurrence for Kraków city for the period from 1965 to 2015 indicated that fog occurred usually on days with non-advective anticyclonic types Ca and Ka or cyclonic and anticyclonic advection from sector S-SW (types SWa and Sc) according to Niedźwiedź classification [27]. This indicated that the majority of winter fogs at Kraków might be related to air pollution from heating during frosty anticyclonic winter weather. Research of long-term variability of the cloud for Kraków has shown that the greatest cloudiness and one of the smallest variabilities are associated with cyclonic situations involving northerly and northeasterly advection. The relationships between the cloudiness and atmospheric circulations were stronger during the cold half of the year than during the warm half, when the radiation factor plays a major role [28]. One of the situations when the influence of atmospheric circulation on air pollution dispersion is well visible is the occurrence of foehn winds, which can worsen or improve the air quality in the city [21].

High weather variability in Poland is associated with frequent movement of low and high-pressure systems [29]. Studies of circulation types for Kraków in the 20th century were summarized by Z. Ustrnul [30]. Significant variation in the annual incidence of individual circulation types according to the Niedźwiedź classification was found; the most frequent were anticyclonic non-directional types (high-pressure center and anticyclonic wedge or ridge); they constituted 15% of all cases during the year. The second most common types were those with the advection of air masses from the western sector (SW-W-NW) (both during anticyclonic and cyclonic situations), with the frequency reaching a total of 40% during the year. The least frequent were the cyclonic types, with the advection of air from the North and the eastern sector. In the individual seasons (from spring to winter), there was quite large variability in dominant circulation types. A slight positive trend was observed for circulation types with air advection from the West. In the 20th century, circulation types were characterized by high inter-annual frequency variability and the absence of distinct, characteristic periods, with the prevalence of certain types.

Assessment of the role of atmospheric circulation on PM₁₀ concentration level during a particular period is highly challenging because many factors including emission level, changeable weather conditions, microclimatic features and chemical and physical processes affect the air quality levels.

Recently, there has been a growing interest in the application of machine learning techniques in statistical analysis [31,32] and forecasting air quality over the wide temporal and spatial scale [33,34]. The most commonly used machine learning tools include Artificial Neural Network [35], Deep Neural Network, Extreme-Gradient Boosting [36] and Random Forests [33]. Machine learning techniques have been successfully applied to assess population exposure to poor air quality in metropolitan areas [33], downscaling of air pollutants at a higher resolution [37], and also meteorological normalization used in air quality trend analysis [31]. Random Forests, which is a machine learning method based on constructing decision trees, is widely used for regression and classification. One of the main advantages of this method, besides its being accurate and straightforward in implementation, is the simple and intuitive way of accessing variables that are important in the process of training the model in complex and nonlinear problems. The research was undertaken in order to assess the quantitative influence of atmospheric circulation on the pollutant concentration in the area of Kraków, and to compare the results for the study period mentioned with the research from earlier decades. The research was also aimed to evaluate two main groups of the classification methods of atmospheric circulation types’ used in studies on air pollution dispersion (described in Section 2.5), and to assess the possibility of their application in air quality management. In our study, we were focused more on the interpretation of the importance of variables by the Random Forests technique than predictions of the model itself. Kraków is an adequate study area for such considerations as, on one hand, it is located in diversified environmental conditions, and on the other hand, relatively long series of air quality measurements are available. Two different classifications of atmospheric circulation were used in the present study: a manual classification by Niedźwiedź [29] and an automatic classification by Lityński [38,39]; those two different classifications were used with the aim of minimizing the risk of misinterpretation of the results, and both classification methods were widely used by different groups of researchers in studies of atmospheric types over Central Europe [14,40,41]. Previous studies concerning the influence of atmospheric circulation on air quality in Krakow, with the application of Niedźwiedź classification, have been summarized in the monograph by J. Godłowska [42]. The research indicated that during air masses advection from S-SW sector and non-directional circulation types (high-pressure center and anticyclonic wedge or ridge), wind speed near the ground was reduced, which in consequence could lead to the occurrence of a high-level PM₁₀ concentration during the cold season. Studies of atmospheric stability in Kraków using SODAR for the period 1994–1999 showed that, for types with air masses advection from a direction between 135° and 225°, the anticyclonic wedge and the high-pressure center, the duration of elevated thermal inversions in the cold season was the longest.

Understanding the relationship between the pollutants’ concentration and the atmospheric circulation, at the synoptic and local scale, is crucial for forecasting air pollution episodes and minimizing the negative impact of air pollution on the health of city residents and on the condition of the natural environment.

2. Material and Methods

Data used in the present study come from different sources and cover the period October 2000–September 2020 (additionally, the period October 2021–December 2021 was selected for operational tests of the predictive model). Data analyses were realized for two sub-periods: cold half-year (October–March) and warm half-year (April–September), owing to significant differences in air pollution emissions and concentrations, which have been observed in Kraków in those sub-periods. First, data sets are described, and their basic statistical features are shown, then the methods combining data from different data sets are presented. Air quality in Kraków was characterized with data on PM₁₀ as the allowed concentrations of that pollutant are exceeded much more frequently than the concentrations of other pollutants. The authors are aware that the division of the year proposed above is only one of the options available, as, in particular years, the frequency of circulation types and meteorological conditions may change significantly; however, such a division was also used in other climatological studies [43,44].

2.1. Study Area

Kraków is the second largest city in Poland, located in the Małopolska (Lesser Poland) region, with an area of 326.8 km² and the number of inhabitants reaching almost 800,000 [45]. The Kraków agglomeration consists of the city itself and the highly populated towns and villages which surround it, and the total number of inhabitants is estimated to exceed 1 million. The city’s area belongs to three different geographical regions and geological structures presented at Figure 1, i.e., the Polish Uplands, the Western Carpathians, and the basins of the Carpathian Foredeep in between. The central part of the city is located in the Wisła River valley, at an altitude of about 200 m a.s.l. In the western part of Kraków, the valley is as narrow as 1 km; however, in the eastern part of the city, the valley widens to about 10 km and there is a system of river terraces (Figure 1b). The hilltops bordering the city to the north and the south reach about 100 m above the river valley floor, similar to the hilltops in the western part of the valley which means that the city is located in a semi-concave landform (open only to the east) and sheltered from the prevailing western winds (Figure 1b). The local scale processes linked to the impact of relief include, for example, katabatic flows, cold air pool formation, frequent air temperature inversions, much lower wind speed in the valley floor than at the hilltops [46]. All the factors mentioned contribute to the poor natural ventilation of the city, and one of its consequences is the occurrence of high PM₁₀ concentration levels, especially during heating seasons.

2.2. Instrumental Meteorological Data

Weather data for Kraków were obtained from the meteorological station located in the Wisła River valley (Balice). The station is administered by the Institute of Meteorology and Water Management—National Research Institute (IMWM-NRI). Measurements of air temperature in the vertical profile of the valley were performed by the Jagiellonian University (JU) at the television mast of EMITEL company, located in the western part of the city (Bokwa, 2010). Measurements of meteorological parameters at the point administered by JU and IMWM-NRI were realized in accordance with WMO guidelines [47]. The location of measurement points and details on weather data used in the study are included in Figure 1 and

Table 1. Location of meteorological and air quality stations in Kraków and its vicinities, and elements used in the study.

No.	Station	Lat N	Lon E	Altitude (m a.s.l.)	Manager of the Station	Landform	Parameters	Data Availability Period	Data Resolution
1	Balice	50.08	19.80	237	IMWM-NRI	Valley bottom	V, D, T, RH, C, PP	1960-currently	1 h, 3 h and 1 day
2	TV mast: 2 m a.g.l. 100 m a.g.l.	50.05	19.90	222 272 322	JU	Valley bottom	T	1.01.2010-currently	3 h
3	Krasińskiego St	50.06	19.93	207	NIEP	Valley bottom	PM10	1.01.2000-currently	1 day

Explanations: V—wind speed (m⋅s⁻¹), D—wind direction, T—air temperature (°C), RH—relative humidity (%), C—cloudiness (oktas), PP—atmospheric precipitation (mm), PM₁₀—mean daily PM₁₀ concentration (µg·m⁻³).

. The measurements from the TV mast were crucial in the analysis of ground thermal stratification in the Wisła River valley at the local scale. These measurements were not available for the whole study period.

2.3. Atmospheric Reanalysis

In order to analyze the stratification of the lower troposphere for the period October 2000–September 2020, ERA5 reanalysis provided by European Center for Medium Range Weather Forecasts [48] was used in this study. Air temperature, relative humidity and wind components data from pressure levels 975, 925 and 850 hPa were applied, at 00:00, 6:00, 12:00 and 18:00 UTC, respectively, for grid point representing Kraków (geographical coordinates 50° N and 20° E). The pressure level 1000 hPa was not used in the analysis owing to the fact that, in some cases, this level could be below the ground level.

2.4. Air Quality Measurements

Data on PM₁₀ concentrations in Kraków come from the databases of the National Inspectorate of Environmental Protection (NIEP) [49]. The methodology for measuring PM₁₀ concentration was realized in accordance with the guidelines of the European Parliament and of the Council included in Directive 2008/50/EC [50]. Daily data from the measurement point located in Krasińskiego St. for the period October 2000–September 2020 were used (Table 1). The measurement point is located in a street canyon, in the city center, at the bottom of the Wisła River valley, with a very busy municipal transportation route and intensive traffic. A comparison of mean daily PM₁₀ concentration from Krasińskiego St. and two other air quality stations: Kurdwanów district and Bulwarowa St., located in the eastern and northern part of the city, for the common period 2010–2020 (3556 days), confirmed a high correlation between measurements from all those points. For the analysis of daily PM₁₀ concentration the Pearson correlation coefficient was used, for three pairs of stations: Krasińskiego–Bulwarowa, Krasińskiego–Kurdwanów and Bulwarowa–Kurdwanów, where correlation coefficients were close to 0.93. The station in Krasińskiego St. is characterized by an increased level of daily PM₁₀ concentration in comparison with other measurement points in Kraków during the year.

The period October 2000–September 2020 is suitable for showing the seasonal, long-term variability of the pollutants’ concentration resulting from changes in the level of air pollutants emissions as well as fluctuations in circulation conditions.

Appendix A summarized air quality measurements used in this study, with special focus on the variability of PM₁₀ daily concentration in cold and warm half-years, number of days with exceedance of selected concentration limits and deseasonalized trend observed in the multiyear period.

2.5. Atmospheric Circulation Classification

According to the suggestions of many authors [29,51,52,53], more than one classification of circulation types was used. Owing to the methodological approach, two classifications of circulation types have been chosen. Each of them represents a different group of atmospheric circulation classifications. Therefore, they differ essentially in many features and, above all, in the method of distinguishing individual types. The first classification included is the traditional, manual approach often used in Poland and developed by T. Niedźwiedź [29]. The second one is an objective classification according to Lityński’s original concept [38]. Taking into account these 2 different classifications allowed for a more objective look at the impact of circulation and its changes in the analyzed 20-year period on the state of the atmosphere, including the concentration of pollution in the study area.

Classification by Lityński is based on an automatic approach, which may be considered, in simplified terms, as the objective one. In the literal sense, it is not like that, because it is based on arbitrarily imposed criteria; however, this approach is different from the manual and obviously subjective one proposed by Niedźwiedź. Both classifications are based on different input data sources. The division of Lityński uses numerical data (currently grid data), while the division of Niedźwiedź is based on the assessment of synoptic daily maps (charts). The spatial scale is also different in both classifications. The Niedźwiedź classification is a typical mesoscale one, while the Lityński classification characterizes the circulation on a larger scale. To sum up, both classifications have a different synoptic approach, and their application seems advisable.

Detailed information circulation classification characteristics for both approaches with analysis of multiyear trends for individual atmospheric patterns are summarized in Appendix B.

2.6. Atmospheric Stratification Determination

Data on air temperature and relative humidity provided by the European Center for Medium Range Weather Forecasts for atmospheric pressure levels 975, 925 and 850 hPa representing the point with geographical coordinates 50° N and 20° E were used to determine the presence of low (layer 975–925 hPa) or upper (925–850 hPa) inversion layers. The atmospheric stratification gradient was determined as the difference between the two nearest levels (layers 975–925 hPa and 925–850 hPa). Lower limits for the occurrence of air temperature inversion and air relative humidity inversion in the lower troposphere were set to 0 °C and 10%, respectively.

Additionally, with the aim of analyzing the near-ground thermal inversion layer, the measurements in the vertical profile (2 m to 100 m a.g.l.) obtained from the TV mast are used. The period of the day has been divided into two sub-periods of equal length:

-

daytime period: from 6 to 17 UTC;

-

nighttime period: from 18 to 5 UTC the next day.

The near-ground thermal gradient was calculated as the difference between lower and upper measurement points. The lower limit of the occurrence of thermal inversion between two levels of the TV mast was set equal to +1 °C. The condition was checked for each time period separately, and then summed up for day and night periods for individual days. Data for the TV mast station were available for the period from January 2010 to September 2020.

2.7. Data Analysis

The data set created to assess the influence of meteorological conditions on air quality includes:

-

Meteorological observations from Balice synoptic station with 6-h resolution: air temperature, relative air humidity, wind speed and direction, cloudiness, the 6-h sum of atmospheric precipitation;

-

air temperature, relative air humidity and wind speed and direction at three pressure levels obtained from ERA5 reanalysis (975, 925 and 850 hPa); differences between neighboring pressure levels of air temperature, relative air humidity, wind speed and wind direction (layers 975–925 hPa and 925–850 hPa) with 6-h resolution;

-

mean daily PM₁₀ concentration from previous day;

-

difference of mean daily PM₁₀ concentration between current day and previous day (used for determining PM₁₀ decrease);

-

day of week;

-

atmospheric circulation types on a certain day according to Niedźwiedź and Lityński classification.

With the aim of investigating the relation between PM₁₀ concentration and meteorological conditions, the Random Forests algorithm was used, which is an ensemble machine learning method based on constructing many decision trees. This method combines a large number of small decision trees into new predictors, and therefore is able to make a better prediction. By using this method, it is possible to assess which variables have the highest importance in machine learning. In our study, we compared results from multilinear regression with stepwise selection and the Random Forests method. Studies of variable selection for Random Forests models were conducted with use of the Boruta method available in package Pomona on GitHub repository [54,55]. In order to provide the best of hyperparametric values, repeated leave-group-out cross-validation (LGOCV) was used. The resampling method LGOCV was available in the function trainControl in the caret R package. For the multilinear regression model, the stepwise Akaike Information Criterion (AIC) algorithm was used [56], which is available in the function stepAIC in the MASS R package. The meteorological database from Balice station and ERA-5 reanalysis and PM₁₀ daily concentration at the previous day were used as an input to models that were trained to predict daily PM₁₀ concentration and its day-to-day changes on a randomly sampled 75% of data. The remaining 25% was used to validate models. Two sets of input data, which differ in the time resolution of meteorological parameters listed above (6-h resolution data and daily averages obtained from 6-h resolution data), were used in the studies. This analysis aimed to answer the question of whether the increase of the temporal resolution of parameters describing weather conditions during the day would improve model accuracy. The hyperparameters tuning and selection of crucial variables was done separately for each Random Forests and multilinear regression model. The plots with variable importance are presented, for clarity only, for the most important parameters. With the aim of determining the partial relationship between daily averages of individual meteorological parameters and the level of the daily PM₁₀ concentration, the optimized Random Forests model was used. Partial dependencies were obtained with use of the function partial_dependence available in the open-source package edarf in the R environment.

In both half-years, days with the worst air quality and days with a significant improvement of air quality in relation to the previous day were selected. That choice of this selection of analyses is due to the fact that such situations are important in terms of the inhabitants’ health protection, but also for various environmental effects. In both groups of the cases selected, very high concentrations of PM₁₀ occur, so the analyses should support the assessment of the atmospheric circulation and weather conditions which contribute to such situations. In the second group of cases, the analyses should additionally support the assessment of the conditions favorable for a sudden decrease of the PM₁₀ concentrations, due to the change of the dispersion conditions. Owing to the fact that the distribution of PM₁₀ daily concentration differs significantly between both half-years (see Figure 2), some criteria of the cases delimitation in both sub-periods differ, too.

The criteria used to distinguish the two groups of days are the following:

-

Group 1: days with high PM₁₀ concentration against the background of a particular half-year, which meet two conditions: daily PM₁₀ concentration is greater than the upper quartile in the selected half-year (see

Table A3. Height of PM₁₀ daily concentration upper quartile in warm and cold half-years.

Year	Cold Half-Year (μg⋅m−3)	Warm Half-Year (μg⋅m−3)
2001	47	41
2002	106	94
2003	137	60
2004	116	60
2005	162	73
2006	145	71
2007	134	77
2008	155	69
2009	123	67
2010	135	57
2011	137	53
2012	131	49
2013	130	47
2014	103	43
2015	115	54
2016	99	50
2017	99	38
2018	87	50
2019	82	42
2020	71	32

) and greater than 50 or 40 μg⋅m⁻³ during cold or warm half-year, respectively.

The number of days meeting the above conditions is 842 and 837 for cold and warm half-years, respectively.

-

Group 2: days characterized by significant PM₁₀ concentration decrease in relation to the previous day, which meet three conditions: the decrease is greater than 25% of the concentration on the previous day, the decrease of PM₁₀ daily concentration is equal at least to 20 or 10 μg⋅m⁻³, in cold or warm half-years, respectively, and days assigned to Group 1 are omitted.

The number of days meeting the above conditions is 634 and 461 for cold and warm half-years, respectively.

With the aim of better understanding the role of atmospheric circulation and local topography on the weather conditions over the analyzed region the distribution of the daily meteorological parameters from Balice station and vertical profiles from ERA5 reanalysis have been analyzed for individual circulation types for cold and half-years separately.

In the next step, the types of atmospheric circulation were assigned to the days from both groups and half-years), and then weather conditions for particular atmospheric circulation types were analyzed. The meteorological parameters used in this part of the study were selected based on the parameter importance obtained from the machine learning analyses. Data used in the study include wind speed and air temperature near the ground, atmospheric precipitation, vertical gradient of air temperature and relative air humidity, air temperature, relative humidity and wind speed at pressure level 925 hPa and wind speed difference in layer between 925 and 975 hPa. The aim of that step was to see whether there are significant differences in weather conditions during a certain atmospheric circulation type occurring in the two groups of days described above and the remaining days.

Selection of half-year for further analyses.

Comprehensive research for both half-years indicated that the problem of air quality in the warm half-year is insignificant compared to the cold half-year. In the analyzed multi-year period, an average share of days with exceedance of the PM₁₀ daily limit in warm half-year was twice as low as in the cold half-year (34% and 67%, respectively). Furthermore, during the period 2016–2020, the average share of days with exceedance of the PM₁₀ daily limit in the warm half-year was equal to 14%, with the lowest values in 2017 and 2020 (8% and 2%, respectively). Days with exceedance of daily PM₁₀ limit occurred mostly in April and September (early spring and early autumn periods), while in June-July such cases almost did not occur. Therefore, only the cold half-year has been described in detail in this paper.

Selection of atmospheric circulation classification.

The analysis of the influence of atmospheric circulation on the dispersion conditions and the level of PM₁₀ concentration was performed for both circulation type classifications (Lityński and Niedźwiedź classifications) for two half-years with particular attention to selected days with the worst air quality and a significant improvement in air quality.

The aim of the research was to determine which circulation type classification better separates the circulation patterns that negatively affect dispersion conditions from the patterns, from those which positively influence the air quality in an urbanized valley. The analysis with the use of both circulation classifications for both half-years and for both groups showed similar dependencies. In order to determine which type of circulation classification is more appropriate for the analysis of air quality in the cold half-year, the Gini coefficient [57] was determined for Niedźwiedź classification (11 and 21 types) and for Lityński classification (27 types). The Gini coefficient has been widely used to measure the inequality among values of a frequency distribution [58,59], the value ranges from 0 to 1. The zero value of the Gini coefficient indicates full uniformity of the distribution. The zero value of the Gini coefficient indicates the perfect equality of the distribution, while the greater the Gini coefficient refers to the greater spread of the distribution. The value of the Gini coefficient was similar for the two Niedźwiedź classifications (0.346 and 0.347), a slightly lower value was obtained for the Lityński classification, equal to 0.338. Owing to this fact and for the better clarity of the article, the paper presents only the analysis for 11 types of Niedźwiedź classification. Similar studies concerning the relation between circulation type classifications and smog days, using Gini coefficient, were conducted in COST Action 733 for air pollution in winter in Polish urban areas [12]. The analysis of the Gini coefficient for individual cold half-years also showed that the variability variation in air quality for individual types of circulation was greater for the types of Niedźwiedź classification (11 and 21 types) than for the Lityński classification (maximum difference was equal 0.037 and average difference equal to 0.010). An additional argument in favor of the selected circulation classification by Niedźwiedź is that it was designed to be most suitable for Southern Poland, while the Lityński circulation classification describes the atmospheric circulation in Central Europe [38].

Schematic representation of the scientific analysis has been presented in Figure 2.

3. Results

The analysis was focused on two groups of days determined for cold half-year:

-

Group 1—days with the highest daily concentration; of PM₁₀;

-

Group 2—days with the greatest decrease day by day in the concentration of PM₁₀.

For both groups of days, the most frequent circulation types according to the Niedźwiedź classification were selected. Weather conditions during days from both groups were compared with remaining days for selected atmospheric circulation patterns.

With the aim of estimating the impact of individual meteorological parameters on air quality, the results of ensemble machine learning methods were used.

All the analyses of the influence of atmospheric circulation on air quality in the light of PM₁₀ confirmed the significant role of circulation types. In the last 20 years, despite a significant reduction in emissions, which is the result of administrative pro-ecological activity, there are still serious smog episodes, when the average daily concentration of PM₁₀ in the cold half-year period may exceed 100 µg⋅m⁻³. This, of course, applies to non-advective circulation types, although surprisingly high dust concentrations may also occur during the types from southerly advection. The results of the performed analyses, i.e., circulation type vs. PM₁₀ concentration is included in the Appendix D.

Random Forests Analyses

At the first step the Random Forests and multilinear regression models were built to predict daily PM₁₀ concentration with the use of two different meteorological data resolution sets (6-h resolution and daily averages from 6-h resolution data). The Boruta variable selection method was applied for Random Forests models, and it showed that for the model which uses daily averages, the following parameters: daily cloudiness, wind direction changes in the layer between 925 hPa and 850 hPa and wind direction at 850 hPa were unnecessary. For the Random Forests model which uses meteorological parameters with 6-h resolution 28 of 111 selected variables were rejected by the Boruta method, including both circulation types, wind direction at 850 hPa, wind direction and wind speed change in layer between 925 and 850 hPa, relative air humidity at 2 m a.g.l. at 0, 12 and 18 UTC, wind direction at 925 hPa at 0, 6 and 18 UTC and day of week. The results of both Random Forests models were similar, the average value of mean absolute error (MAE) and root-mean square error (RMSE) for both models were equal 19.6 μg⋅m⁻³ and 26.9 μg⋅m⁻³, respectively. The Random Forests models analysis for specific measured PM₁₀ concentration ranges 0–50 (25% testing data), 50–100 (40% testing data) and 100–200 (27% testing data) indicated that the RMSE error was equal respectively to 18, 20, and 32 μg⋅m⁻³. The group of observations with PM₁₀ concentration exceeding 200 μg⋅m⁻³ was relatively small (5% testing data), and the RMSE for this group was the largest equal to 42 μg⋅m⁻³. An example plot presenting comparison observations with the Random Forests model forecast is included in Appendix C, Figure A6a). In Figure 3 the most important parameters affecting daily PM₁₀ concentration are presented. Air quality on the previous day (PM₁₀ daily concentration) was the most important parameter for both models. For the sake of clarity of Figure 3, this parameter was not presented in the chart owing to the large differences between the GINI Index for this parameter and the next one. The analysis of variable importance for both Random Forests models confirmed the similarity of the results. Apart from the parameters whose significant influence on air quality is well established (air temperature and wind speed at the ground and air temperature gradient), the key factor was also the gradient of relative air humidity and wind shear in the lowest troposphere (layer between 975 and 925 hPa; Figure 3a).

With the aim of analyzing the significance of individual parameters concerning the model accuracy, tests by removing a single variable from the model were performed (Figure 4). These studies also confirmed that the most important parameter was the PM₁₀ concentration level on the previous day. The lack of this variable in the model affected on MAE increased by almost 25% compared to the forecast results where all parameters were included.

Studies of multilinear regression for both data groups were done for the same teaching and testing sets. Variable criterion with use of Akaike algorithm showed that for data with daily averages of meteorological parameters, six variables were excluded in the analysis: relative air humidity at 2 m a.g.l., wind speed at three pressure levels (975 hPa, 925 hPa and 850 hPa), relative humidity gradient in layer between 925 and 850 hPa and relative humidity at 850 hPa. The results obtained for multilinear models were slightly worse for the Random Forests model which used the same data set, RMSE and MAE for multilinear models were equal to 29 μg⋅m⁻³ and 21.4 μg⋅m⁻³, respectively. A comparison of results obtained from four Random Forests and multilinear regression models is presented in the Taylor diagram in Figure 5. The analysis of results presented at Figure 5 indicates that there are slight differences between the two model groups in Pearson correlation coefficient and standard deviation of predicted values. In the case of multilinear regression and meteorological parameters with 6-h resolution, the results of the multilinear model were close to the previous multilinear model. The application of the Akaike method to select the best parameters for the multilinear model caused a reduction of the parameters from 111 to 44 variables. The selection of crucial variables did not improve model performance, RMSE and MAE were equal to 29 μg⋅m⁻³ and 21.4 μg⋅m⁻³, respectively. In this case the application of the Random Forests model for numerous variables showed better results than multilinear regression. The comparison of observations with the multilinear model forecast with the use of daily averages of meteorological parameters presented at Figure A6b) in Appendix C indicates that the model underestimates PM₁₀ daily concentration for values below 25 μg⋅m⁻³. In contrast, Random Forests models more often overestimate PM₁₀ concentration than multilinear regression in the range between 0 and 50 μg⋅m⁻³ (Figure A7 in Appendix C).

In the second part, Random Forests and multilinear regression models were used to predict day-to-day changes of PM₁₀ daily concentration. The same database as presented above was used in these studies, including measurements from the Balice station, ERA-5 reanalysis, two circulation types, day of the week, month, day of the year and PM₁₀ daily concentration at the previous day. The analysis of Random Forests and multilinear regression models showed similar results, the values of RMSE and MAE were equal on average to 30 μg⋅m⁻³ and 20 μg⋅m⁻³, respectively. For this case hyperparameter tuning and parameter selection for Random Forests models did not significantly improve model accuracy (change of RMSE and MAE did not exceed 3%). It is worth to mentioning that variable selection with the use of the Boruta method for the Random Forests model which uses data with 6-h resolution were reduced from 110 to 27 variables. An example plot presenting comparison observations of day-to-day changes with the Random Forests model forecast is included in Appendix C, Figure A8b). In the case of a multilinear linear regression model with the same data set, stepwise the Akaike method selected 46 variables from 110 available. Verification of four models presented with the use of a Taylor diagram in Figure 6b indicates that the differences between them are negligible; however, comparison of density curves for Random Forests and multilinear regression with observations indicates that also the Random Forests model often predicts day-to-day PM₁₀ concentration changes in a range between −10 and 25 μg⋅m⁻³ (Figure A9 in Appendix C). For both models groups the most important parameter was the PM₁₀ concentration level on the previous day. Figure 6 presents parameters importance for Random Forests models affecting day-to-day PM₁₀ concentration changes.

Analyses of the importance of variables presented at Figure 6 showed that the most important parameters were air temperature and wind speed near the ground and at the closest pressure level 975 hPa. It is worth mentioning that relative air humidity and relative air humidity gradient were more crucial parameters than air temperature gradient in the layer between 975 and 925 hPa concerning the prediction of day-to-day PM₁₀ changes (Figure 6a). Two additional sensitivity tests were performed. Firstly, one for the whole period, without dividing the data set into two half-years and another one with training the model with data from 2000 to 2015 and testing it with data from 2016 to 2020 (also without dividing data set into cold and warm half-years). In both cases we achieved a similar order of importance of parameters for both Random Forests models as in Figure 4, while scores were slightly improved, e.g., with a decrease of RMSE around by 7 μg⋅m⁻³ for both models. It can be explained by adding a warm half-year to the data set that is characterized by lower values of PM₁₀ concentration level. Results obtained by multilinear linear regression models were slightly worse for both tests, mean differences between Random Forests models were equal to 3 μg⋅m⁻³ for RMSE. As mentioned before, the main motivation for using Random Forests was to determine which meteorological parameters should be considered for further analysis, but a comparison of the accuracy of our forecasts with similar models (both physical and based on machine learning) shows also the good predictive potential of such an approach [60,61,62].

The optimized Random Forests model built to predict daily PM₁₀ concentration based on the daily averages presented above was used to analyze the partial relationship between individual meteorological parameters and the level of the PM₁₀ daily concentration. Figure 7 presents the partial dependence of predicted PM₁₀ daily concentration for selected meteorological parameters. The plots show the value of the lower and upper quartiles in the analyzed cold half-years (dashed vertical lines). Detailed analysis of the daily thermal gradient between 925 and 975 hPa, has shown that predicted PM₁₀ concentration did not differ significantly for the positive values of vertical gradients (the number of such days in all analyzed cold half-years did not exceed 25%—Figure 7b). For the range of the daily gradient values between −3 °C and 0 °C in the layer between 925 and 975 hPa, the influence of this parameter on the daily value of the PM₁₀ concentration increases significantly. During the days without absence of the elevated inversion in the layer between 850 and 925 hPa, the predicted value of PM₁₀ concentration was close to the minimum value (low statistical importance). When the daily temperature gradient in the layer between 850 and 925 hPa decreases below −3 °C, a significant increase in the predicted pollutant concentration can be observed. The plots of the dependence of the air humidity at the ground (Figure 7d) and at the pressure level of 925 hPa (Figure 7e) on predicted PM₁₀ concentration has shown a different relationship. For the days when daily relative humidity of 2 m a.g.l. exceeds 80%, there is a linear increase of the predicted daily PM₁₀ concentration. On the other hand, with a decrease in relative humidity at the height of 925 hPa the predicted PM₁₀ concentration increases. The plots of predicted PM₁₀ concentration from the relative humidity gradient in the layer between 925 and 975 hPa have shown gradual deterioration of air quality with the decrease of humidity gradient in the range of from 0 to −25%. The relationship between the average daily wind speed at 10 m a.g.l and predicted air pollution level presented at Figure 7g indicates a strong decrease of PM₁₀ concentration for the wind speed in the range from 0 up to 5 m⋅s⁻¹; above this value the increase of wind speed did not significantly improve the air quality in the city. The relationship between the vertical wind shear in the layer between 925 and 975 hPa and the pollution concentration is presented in Figure 7h,i. When wind shear increases wind speed in the vertical profile, the crucial point is the exceedance of 5 m⋅s⁻¹. For such situations the increase in the speed difference between layer 925 and 975 hPa does not significantly improve the air quality. When the wind shear is associated with a significant change in the wind direction between the level of 925 and 975 hPa, an increase in the difference in wind directions negatively affects the predicted air quality in the valley.

4. Discussion

The statistical analysis of the impact if meteorological conditions and atmospheric circulation types on air quality in Kraków with the use of machine learning methods made possible an objective selection of crucial parameters influencing the pollutant concentration level. In the studies presented, we have compared results from multilinear regression and the Random Forests method to predict daily PM₁₀ concentration and its day-to-day changes for two sets of input data which differed in the temporal resolution of meteorological parameters. The application of 6-h resolution meteorological data in comparison with daily averages to predict daily PM₁₀ concentration and its day-to-day changes showed similar results for both methods. This confirms the statement that the use of daily averages of meteorological parameters is sufficient to predict PM₁₀ daily concentration a day ahead. Studies of the importance of variables’ in predicting PM₁₀ day-to-day changes the with the use of 6-h resolution data indicated that the number of crucial parameters was significantly lower than for predicting PM₁₀ daily concentration for the Random Forests model (equal to 27 from 111 possible variables). It is also worth mentioning that in the case of predicting PM₁₀ concentration with the use of numerous variables (6-h data resolution), the multilinear regression model significantly reduced the number of variables while the Random Forests model selected more variables as important. Analyses of model performance showed that for this case Random Forests results were better than for multilinear regression models.

Additional tests with the use of measurements from October 2000 to September 2020 were carried out to estimate the impact of changes in meteorology and emissions during the recent cold season (October 2021 to December 2021). Eight different sets of training data were prepared for tests for Random Forests models and multilinear regression models (sets of two temporal resolutions: 6-h and daily averages; data for cold half-years only and for both half-years; training data for period from October 2000 to September 2020 and from January 2015 to September 2020). Analysis of the shorter period for model training was done to answer the question how models’ performance is affected. An analysis of Taylor diagram plots (Figure A10 in Appendix C) indicated that the multilinear regression models had a higher standard deviation than the observations, indicating an excessively high variability of predicted PM₁₀ concentration. The results obtained from Random Forests models were closer to the observations than the results from linear regression models (lower value of RMSE and standard deviation closer to the observation’s). Secondly, using shorter periods for model training showed better results (lower RMSE and lower overestimation of PM₁₀ concentration), however in individual smog episodes PM₁₀ daily concentration was underestimated in comparison with forecasts using a longer training data set (Figure A11 in Appendix C). Analysis of the time course of predicted PM₁₀ concentration showed that the multilinear model in some cases overestimated the PM₁₀ decrease, while the Random Forests model performed better.

The conducted experimental studies (based on data from 2021) have shown that such analyses should take into account the circular data with greater resolution than the daily one. Considering only one type of circulation for the whole day does not make it possible to take into account the dynamics of circulation changes, which is particularly high over Central Europe in winter. On the basis of the existing classifications, especially the local one, a method of assessing the atmospheric circulation should be developed, taking into account the daily course (at least with 3 h resolution).

Furthermore, more detailed analysis of the importance of individual parameters on PM₁₀ daily concentration level was available with the use of the Random Forests model. The results obtained were used in further analysis to investigate the dispersion of selected parameters for individual types of circulation.

The study on the influence of atmospheric circulation patterns on air quality made it possible to distinguish the dominant circulation types during which the probability of occurrence of poor air quality (Group 1) and a significant improvement in air quality conditions (Group 2) was the highest. Days with the high PM₁₀ concentration at cold half-year, occurred mostly during the advection of air masses from the S-SW sector, non-directional anticyclonic situations (Ca + Ka type) and also anticyclonic situations with air advection from the W-NW sector. Such days were characterized by lower wind speed and air temperature at ground level and greater stability of the atmosphere during the day and night periods in comparison with days not assigned to both special groups (remaining days). According to the Mann–Whitney U test, the distribution of daily sums of precipitation was similar for dominant circulation types for days in Group 1 and remaining days, but the frequency of precipitation was lower for days in Group 1. Furthermore, during the daytime for days with high PM₁₀ concentration, a local minimum of relative air humidity at level 925 hPa occurred frequently. The partial dependence of meteorological parameters obtained from the Random Forests model has also confirmed the negative effect of strong negative relative humidity gradient in layers between 925 and 975 hPa on air quality. In this case, advection of dry air masses at a height of 925 hPa, especially frequent for the S + SWa and W + NWa types, contributed to the increase in the stability of the atmosphere in the valley and resulted in a longer persistence of humid cold air pool. During the winter, when foehn wind occurs, there is often the advection of warm and dry air masses above the analyzed region [63]. Additionally, previous studies pointed out that circulation types S + SWa, Ca + Ka enhance the occurrence of fog in Kraków which confirms the poor air pollution dispersion conditions linked to those circulation types [64]; the occurrence of haze and fog episodes is studied widely in different parts of the world in the context of air pollution control [65].

The second case studied consisted of days characterized by a significant improvement of air quality (Group 2). A significant reduction in the PM₁₀ daily concentration occurred mostly for three circulation patterns (air advection from W-NW sector and cyclonic type with differentiated air advection—Cc + Bc type). These days were characterized by increased wind speed and a greater share of days with precipitation in comparison with remaining days. Atmospheric stratification (relative air humidity and air temperature gradient) was similar for dominant circulation patterns for days assigned to Group 2 and to remaining days. Near-ground temperature inversion for days in Group 2 during daytime almost did not occur, which was confirmed by local measurements and ERA5 reanalysis. It is also worth mentioning that during these days, the local maximum of relative air humidity in the layer 925 hPa during daytime occurred frequently.

Days with anticyclonic conditions with air advection from the W-NW sector are characterized by changeable weather conditions, which contributed to improvement or deterioration of air quality conditions during the cold half-year. The improvement is observed when air masses could penetrate into the valley and remove the cold air pool, while deterioration can be seen when air masses pass over the valley. Studies conducted over the region of the Dead Sea Valley using a high resolution WRF model [66] indicated that foehn wind intrusion into the valley depends on synoptic and mesoscale conditions which affect the vertical structure of the lower troposphere. For the cases with a high stable layer over the Dead Sea Valley, the foehn reached the valley floor, while during a low stable layer, it did not.

Studies of air quality in cold and warm half-years have shown that weak wind speed is one of the most important factors which deteriorates air quality. Owing to this fact circulation patterns which are characterized by weak wind speed, caused by the interaction of local orography (air advection from S-SW and W-NW sectors) and also atmospheric stagnation were the most important (non-advective types with anticyclonic situation, according to Niedźwiedź: type Ca + Ka). The high importance of wind speed on air quality was confirmed by numerous previous studies [3,67,68]; however, studies of the partial dependence of meteorological individual parameters for daily PM₁₀ level showed nonlinear dependency [31].

Furthermore, studies of individual meteorological parameters have shown that vertical wind shear can worsen but also improve air quality in the valley. An increase of wind speed difference between the layers 925 and 975 hPa had a positive impact on air quality. On the other hand, strong wind shear associated with a change of the wind direction in vertical profile affects the deterioration of air quality by reducing the height of the mixing layer during the daytime. The study of the PM₁₀ concentration vertical profiles in Kraków presented in the work of Sekula et al. 2021 [69] indicate that this phenomenon often occurs at the valley bottom height (approx. 100 m a.g.l.). During the cold half-year, poor dispersive conditions are more frequent than in the warm half-year, which in combination with high rates of emission from the residential sector led to accumulation of pollutants inside PBL. Analysis of the deseasonalized multiyear PM₁₀ trend has shown that in the decade 2011–2020 a negative trend was observed which may be linked to the positive trend of air temperature.

According to the Random Forests model, adding a vertical gradient between neighboring pressure fields improved the quality of the PM₁₀ level forecast. Other studies concerning application of machine learning methods in air quality forecasting confirmed that meteorological parameters like, wind speed, air temperature, relative air humidity and atmospheric precipitation were important factors affecting air quality [70]; however, studies of the effect of atmospheric precipitation on the concentration of particulate matter showed that it mainly washes out coarse particles while having little effect on fine particles [71]. Attention should also be paid to the representation of the atmospheric stability in the machine learning models; in this research we can distinguish two approaches: the application of planetary boundary layer height [31,72] or a more complex one with the application of meteorological conditions at different atmospheric pressure levels [73]. Owing to the fact that the estimation of planetary boundary layer height in numerical atmospheric models still requires validation and further development [74,75] we would like to suggest applying vertical profiles of atmosphere rather than PBL height in studies.

5. Conclusions

The analysis of air quality conditions in the multiyear period has proved that wind speed, air temperature, atmospheric stability connected to relative humidity gradient and air temperature gradient at lower troposphere and the occurrence of precipitation significantly influences pollutant concentrations. Apart from the non-directional anticyclonic conditions which affect air stagnation, air advection from the S-SW sector, strongly modified by local topography, has usually caused an increase of PM₁₀ levels in the study area. Studies have shown that for the region analyzed the direction of air advection and its intensity is of greater importance than the type of pressure system concerning the impact on PM₁₀ levels. Certain types of circulation can be indicated as significant both in terms of improving the dispersion of pollution and its deterioration; this is the result of the modification of large-scale processes by orography and near-ground atmospheric conditions. For example, air masses advection from the W-NW sector may strengthen near-ground thermal inversion and reduce wind speed in the valley, but it can also break thermal inversion in the valley and topographically channel the air flow. Research has indicated that particular types of circulation may affect the deterioration of air quality conditions in the cold half-year. During these circulation types lasting for a few or several days, a continuous increase of air pollution can be observed. Sometimes it leads to extremely high values of PM₁₀ concentrations (e.g., types S-SW, W + NWa).

The analysis of the number of days with PM₁₀ levels exceeding the daily limit in the study period showed that the emission reduction contributed to a significant improvement in the air quality in the city; however, the occurrence of days with poor air quality in the future is very likely due to the strong influence of meteorological conditions on that element. The number of days with low thermal inversion in the 975–925 hPa layer in a cold half-year turned out to be particularly important. Significant factors influencing the improvement of air quality in the cold half-year were the occurrence of longer rainfall (rainfall during the day and night), high daily wind speed in the valley and negative air temperature gradient.

One of the limitations in the studies presented above is the assignment of a single circulation type to the whole day; on days when an atmospheric front or the pressure center passes over a certain area, the meteorological conditions may change significantly during the day. Therefore, for the detailed analysis of atmospheric circulation, the daily fluctuations of circulation conditions should be taken into account. Currently, studies on application of the Convolutional Neural Network in automatic classification of atmospheric circulation according to the Niedźwiedź classification with the use of ERA5 reanalysis are conducted in our research group. The first results obtained are promising, however some model optimizations are still necessary.

Analysis of hourly PM₁₀ concentration data and meteorological parameters with the use of cross correlation function have shown the occurrence of delayed time response of PM₁₀ concentration level in the city to the change in meteorological conditions. For instance, at Krasińskiego station, the delayed time response obtained for the PM₁₀ level for the wind speed, wind gusts, air temperature, as well as the ground thermal gradient was equal to 2 h; however, it should be mentioned here that air quality in the city may vary significantly on spatial and temporal scale as it was presented in other studies [8]. Further research on the intra-city spatial dependency from meteorological parameters and circulation patterns is necessary from the point of view of habitability and health risk.

Previous studies indicated that each technique of circulation patterns classification has some limitations e.g., there is a problem of equally sized classes, separation of different types, seasonal or inter-annual variability of a class frequency. In the case of the subjective classifications, there are high inter-annual variability and larger long-term trends of the frequencies of the types’ in comparison to the automated circulation classification methods; however, subjective classification includes important expert knowledge concerning analyzed geographical regions, which is difficult to formulate in precise rules for automated classification methods [52]. In conclusion the authors would suggest using local circulation classification methods in studies for different regions, owing to the effect of topography on modifications of atmosphere dynamics; however, because of the obvious limitations of the use of manual approaches, connected to their subjective nature, it seems that the best solution would be to use a local subjective classification of circulation types, which could be automated. Such approaches are known in the literature, although, they were applied in larger spatial scales [51].

Owing to the fact that machine learning methods create great opportunities in air quality studies, further development works are planned using the Random Forests method to analyze and forecast air quality on a larger spatial scale (e.g., cities in Central Europe) by supplementing the model with additional data such as land cover, topography, and turbulence parameters, as well as the results from operational forecasts of numerical air quality models to improve model accuracy. The further step in air quality studies will be an application of multi-step time series forecasting to model daily cycle of air pollution but also to predict daily pollution levels for three days ahead by using weather forecast and air pollution levels on the current day [75,76]. The next direction of development of the current research focuses on the analysis of spatial and temporal variability of pollution for large cities using data from air quality stations as well as non-governmental air quality systems. The first tests of using convolutional neural networks to determine air pollution level with respect to circulation patterns over larger domains are very promising, and so it is also planned to further investigate those methods.