Winter and Summer PM_2.5 Land Use Regression Models for the City of Novi Sad, Serbia

Abstract

In this study, we describe the development of seasonal winter and summer (heating and non-heating season) land use regression (LUR) models for PM_2.5 mass concentration for the city of Novi Sad, Serbia. The PM_2.5 data were obtained through an extensive seasonal measurement campaign conducted at 21 locations in urban, urban/industrial, industrial and background areas in the period from February 2020–July 2021. At each location, PM_2.5 samples were collected on quartz fibre filters for 10 days per season using a reference gravimetric pump. The developed heating season model had two predictors, the first can be associated with domestic heating over a larger area and the second with local traffic. These predictors contributed to the adjusted R² of 0.33 and 0.55, respectively. The developed non-heating season model had one predictor which can be associated with local traffic, which contributed to the adjusted R² of 0.40. Leave-one-out cross-validation determined RMSE/mean absolute error for the heating and non-heating season model were 4.04/4.80 μg/m³ and 2.80/3.17 μg/m³, respectively. For purposes of completeness, developed LUR models were also compared to a simple linear model which utilizes satellite aerosol optical depth data for PM_2.5 estimation, and showed superior performance. The developed LUR models can help with quantification of differences between seasonal levels of air pollution, and, consequently, air pollution exposure and association between seasonal long-term exposure and possible health risk implications.

1. Introduction

Particulate matter pollution is known to be the single largest environmental health risk in Europe [1] and various adverse health effects are linked to the exposure to particulate matter (PM) [2,3]. With the continuous rapid growth of urban areas, dense air quality monitoring within these environments becomes increasingly important since significant contrasts in ambient air quality and level of PM_2.5 can exist within different city microenvironments and consequently manifest different degree of health risks [4,5].

However, sparse referent monitoring networks and limited geographic coverage, which makes capturing intra-city variabilities difficult, are the reality of most countries in the world [6]. Nowadays there is a wide palette of different techniques used for processing air quality data with the goal of estimating spatial variations in measured pollutants, and assessing exposure, each having different requirements for data availability, going from more simple spatial interpolation methods to the more advanced dispersion models that integrate chemical transport modelling.

Data-driven statistical modelling of air pollution generally has several advantages in the specific context of land use regression (LUR) modelling. LUR modelling has a solid foundation since it relies on directly measuring air pollution using typical reference instruments, and it allows for approximate modelling of pollution sources/sinks or their proxies on both smaller and larger scales of spatial variation but does so without introducing complexities of the physical–chemical modelling and dependencies of extensive datasets. This makes, in a way, LUR models robust and capable of quantifying the level of air pollution and air pollution exposure even at the personal level [6,7,8,9,10,11,12,13,14,15,16,17,18,19]. Many LUR models were developed for densely populated cities where high levels of air pollution were detected. Those studies focused on the development of either local or regional LUR models [20,21,22]. A well-known LUR model targeting multiple European cities was developed in the framework of the European Study of Cohorts for Air Pollution Effects—ESCAPE project [14,23,24]. The most commonly used target variables for LUR modelling in the context of air pollution are the PM_2.5 fraction of particulate matter pollution, and the concentration of gaseous pollutant, NO₂ [25]. These models can often provide a high percentage of explained variability. Recent publications that utilized LUR models for predicting PM_2.5 have implemented different advanced machine learning algorithms, satellite data or other approaches in combination with basic LUR modelling principles in order to enhance model prediction efficiency and ability to be used for a larger study area [15,17,18,26,27,28,29,30,31,32].

Monitoring air pollutants for LUR development is often time-limited to several sampling campaigns during the year, lasting one or two weeks [33,34,35]. Seasonal PM_2.5 LUR models have been developed in Asia [36], Canada [37], Iran [38] and Colombia [34]. Previously, in Serbia, adaptations of the European LUR model developed by Wang et al. in 2014 within the ESCAPE study [24] and its validation with the development of supporting average daily traffic intensity models, and data fusion simulation were performed in Belgrade for PM_2.5 by Davidović et al. [39] and for NO₂ by Davidović et al. [40]. Some aspects of the application of the European LUR model developed by Wang et al. [24] in the Novi Sad municipality were shown by Dmitrašinović et al. [41]. LUR models are usually developed based on 20–80 sites [24]; however, in the ESCAPE project, 40–80 and 20–40 sites in 36 study areas in Europe were used for NO_X and PM_2.5 model development [42]. The smaller the number of measurement sites used for assessing the level of pollutant concentration, the more the accuracy of the developed LUR model will depend on the position of the individual measuring points [43].

LUR model development typically has the following components: the target variable specifications (specific air quality measurements (e.g., NO₂, PM_2.5) we want to predict or estimate), specific air pollution campaign data or existent air pollution monitoring network data, land use and land cover maps (e.g., forests, residential areas, roads within the study area), distance to feature/buffer zones (used to model the correlations those features have with pollutant concentrations at a specific location), and other predictor variables [44]. Once a statistically significant relationship between predictor variables and a target pollutant is established, this relationship can be used to estimate pollutant concentrations at various unmonitored sites. This is possible as long as all the predictor variables can be calculated for those sites [18,26,31,44,45,46].

In this study, we had a goal to develop seasonal winter and summer (heating and non-heating season) PM_2.5 land use regression (LUR) models for the city of Novi Sad, Serbia. One of the main motivations behind such a modelling effort is filling in the knowledge gap on spatial and seasonal variations in PM_2.5 concentrations since the network of reference monitoring stations in the city is rather unevenly distributed. In addition, in most cities in Serbia, PM concentrations show drastic differences in heating and non-heating seasons, thus the origin and quantification of this potential difference in the city of Novi Sad needs confirmation. Once in place, these models can help with the quantification of the level of air pollution, air pollution exposure, and the association between long-term exposure and various health risks, even in microenvironments not covered by regulatory monitoring. Furthermore, these models can also be used to obtain high-resolution seasonal base maps for PM_2.5 concentrations in the city of Novi Sad.

2. Materials and Methods

2.1. Study Area

Novi Sad, Serbia’s second largest city after the capital, Belgrade, is an urban-industrial agglomeration divided by the Danube River into two distinct regions. The city encompasses an area of 702.7 km² and experiences a regional climate that transitions from moderately continental to continental. Novi Sad sits at an elevation of approximately 80 m above sea level [47]. According to the latest publicly available population census data from 2022 [48], the city has 367,121 inhabitants.

Urban air quality in Novi Sad is influenced by several factors, including residential heating using wood, coal or natural gas, intensive traffic in urban areas, low energy efficiency of plants in the energy sector and industry, diffuse pollution from agriculture, and others [49]. Fine particulate matter pollution in Novi Sad frequently exceeds the proposed daily and annual WHO limits of 15 µg/m³ and 5 µg/m³, respectively [50]. In Novi Sad, PM_2.5 is measured by two networks, the national network running under the Serbian Environmental Protection Agency (SEPA) (one measuring site, urban) and a local network established by the City Administration for Environmental Protection (two measuring sites, urban and industrial) [51]. The study area and sampling sites were categorized into urban zone (URB), industrial zone (IND), border of urban and industrial zone (URB/IND) and background zone (BCG). Figure 1 shows a photo of each of the sampling location types.

2.2. Sampling Campaign

The sampling campaign was conducted during winter (15 February–18 March 2020, 8 December 2020–19 February 2021) and summer (24 July–18 September 2020, 6–16 July 2021) periods and PM_2.5 was measured at 21 locations. At each location, PM_2.5 was collected for 10 days using a reference gravimetric pump Sven Leckel Model LVS3, with standard PM_2.5 inlet and 2.3 lpm flow rate. Gravimetric pumps were set to sample air for 48 h and collect particles on quartz fiber filters (Whatman QM-A, 47 mm). We had a total of 4 gravimetric pumps, thus particle concentrations during the winter and summer seasons were collected in batches of 4 sites simultaneously, for 10 days at each site per season. Each sample was collected for 48 h, resulting in 5 samples at each site per season. This cycle was performed until all measurements were completed at all 21 sites. Based on these results, seasonal averages were calculated according to the ESCAPE exposure manual [52].

Figure 2a shows 21 measuring locations, divided by colour according to type (URB, IND, BCG, URB/IND). Of these locations, 18 can be classified as URB-type locations with different levels of traffic intensity. During the heating season, air quality strongly depends on heating modalities present in urban microenvironments (individual heating (IH) using wood, coal, gas and electricity, or district heating by thermal power plants (DH), see

Table 1. Locations used for PM_2.5 sampling in two seasons: location types, coordinates and sampling periods.

Location Number	Measuring Location	Coordinates		Location Type/Traffic Type/Heating Type Options	Sampling Period
		Latitude	Longitude
1	Veselina Maslese	45.26928	19.83304	URB/vicinity of roundabout, traffic site/DH area	8 September–18 September 2020 8 December–18 December 2020
2	Temerinski put	45.29313	19.83101	URB/traffic site/IH area	26 January–5 February 2021 6 July–16 July 2021
3	Svetozara Savkovica	45.28565	19.82087	URB, IND/resident area, low traffic/IH area	15 February–25 February 2020 24 July–3 August 2020
4	Sopocanska	45.26471	19.83378	URB/resident area, low traffic/IH area	26 February–7 March 2020 4 August–14 August 2020
5	Sirine I	45.22417	19.86177	URB/resident area, low traffic/IH area	8 September–18 September 2020 26 December 2020–5 January 2021
6	Save Ljubojeva	45.26074	19.7998	URB/strict resident area, no traffic/DH area	26 February–7 March 2020 4 August–14 August 2020
7	Rumenacki put	45.24968	19.82467	IND/high traffic site, vicinity of roundabout/IH area	8 March–18 March 2020 15 August–25 August 2020
8	Rackog	45.23406	19.88453	URB/traffic site/IH area	15 February–25 February 2020 24 July–3 August 2020
9	Pavla Ivica	45.25226	19.83725	URB/resident area, low traffic/DH area	4 August–14 August 2020 11 January–21 January 2021
10	Partizanska	45.24932	19.87728	URB/traffic site/IH area	8 September–18 September 2020 8 December–18 December 2020
11	Mise Dimitrijevica	45.24814	19.83448	URB/traffic site/DH area	15 February–25 February 2020 24 July–3 August 2020
12	Juraja Krizanica	45.27573	19.80066	URB/traffic site/IH area	8 March–18 March 2020 15 August–25 August 2020
13	Jevrejska	45.22818	19.89271	URB/high traffic site, vicinity of intersection/DH area	26 August–5 September 2020 11 January–21 January 2021
14	Janka Cmelika	45.26472	19.80861	URB/traffic site, low traffic/DH area	15 February–25 February 2020 24 July–3 August 2020
15	Hajduk Veljkova 4	45.24167	19.80375	URB/traffic site, vicinity of intersection/border of DH and IH area	8 March–18 March 2020 15 August–25 August 2020
16	Hajduk Veljkova 11	45.25749	19.8238	URB/traffic site/DH area	8 September–18 September 2020 9 February–19 February 2021
17	Jirecekova	45.26078	19.80913	BCG/vicinity of the river, no traffic/DH area	26 August–5 September 2020 8 December–18 December 2020
18	Darinke Grujic	45.24124	19.85262	URB/residential area, low traffic/IH area	26 August–5 September 2020 8 December–18 December 2020
19	Bulevar PP	45.23936	19.78643	URB/traffic site/border of DH area and IH area	8 March–18 March 2020 15 August–25 August 2020
20	Beogradski kej	45.25608	19.85465	URB/traffic site/DH area	26 January–5 February 2021 6 July–16 July 2021
21	Balzakova	45.23787	19.83439	URB/residential area, low traffic/DH area	26 February–7 March 2020 4 August–14 August 2020

). Location no. 7 is in the IND zone, location no. 3 is on the very border of the URB/IND zone, and location no. 17 is in the vicinity of the local background station, thus it is treated as a BCG site. The background location is situated near the Danube River, a residential area heated by a DH and near a pedestrian zone, with no traffic or domestic heating nearby. In Figure 2b, DH area borders and coverage areas (orange, pink, blue, and yellow colour) are shown and every study site is numbered (from 1–21), following the site sequence given in Table 1. Sampling sites are coloured green and white, that is, sites marked with the green colour are within areas that are within the city DH area, and sites marked with the white colour, are within areas with mixed means of individual heating (denoted as IH area in subsequent text).

2.3. Model Development and Evaluation

The database containing all PM_2.5 measurements from the 21 sites and a set of potential predictor variables was used for LUR model development. ESCAPE methodology [24,38,53,54] was applied, including enforcement of potential predictor direction of effect, as listed in

Table A1. List of potential predictor variables, units, buffer size and direction of effect of target variable PM_2.5.

Land Use/Land Cover Database	Database Name/Tag	Predictor Variable	Unit	Buffer Size	Direction of Effect
Corine land cover	Class1		m2	300, 1000, 3000 m	+
	Class2				+
	Class3				-
	Class4				-
	Class5				-
Open street maps	landuse = building	building_area	m2	300, 1000, 3000 m	+
	landuse = park	park_area			-
	landuse = industrial	industrial_area			+
	landuse = forest	forest_area			-
	landuse = farmland	suburban_area			+
	landuse = residential	resident_area			+
Local traffic density	A_class		number of vehicles per day−1 m	25, 50, 100, 200, 300, 500, 1000, 1500 m	+
	B_class				+
	C_class				+
Open street map (Road length and classes)	ROADLENGTH	ROADLENGTH	m	25, 50, 100, 200, 300, 500, 100, 1500 m	+
	PRIMARYROADLENGTH	PRIMARYROADLENGTH			+
	SECONDARYROADLENGTH	SECONDARYROADLENGTH			+
	TERTIARYROADLENGTH	TERTIARYROADLENGTH			+
	RESIDENTROADLENGTH	RESIDENTROADLENGTH			+
	DISTINVMAJOR1 and DISTINVMAJOR2	DISTINVMAJOR1 and DISTINVMAJOR2	m−1, m−2	Inverse distance and squared inverse distance	+
	DISTINVSECOND1 and DISTINVSECOND2	DISTINVSECOND1 and DISTINVSECOND2			+
	DISTINVTERT1 and DISTINVTERT2	DISTINVTERT1 and DISTINVTERT2			+
	DISTINVNEAR1 and DISTINVNEAR2	DISTINVNEAR1 and DISTINVNEAR2			+
Distance to the river	River Danube		m	NA	+
Population	Total population around measuring location, within a buffer		Number of inhabitants	300, 1000, 3000 m	+

Note. Class1—artificial area, Class2—agricultural areas, Class3—forest and semi natural areas, Class4—wetlands, Class5—water bodies, suburban_area—area of households, outside strictly urban areas, with household heating on solid fuels and private central heating, WIN_A—mean daily winter traffic densities for class A on all roads in a buffer (sum (traffic density × length of all road segments)), WIN_B—mean daily winter traffic densities for class B on all roads in a buffer (sum (traffic density × length of all road segments)), WIN_C—mean daily winter traffic densities for class C on all roads in a buffer (sum (traffic density × length of all road segments)), ROADLENGTH—length of all roads in a buffer, PRIMARYROADLENGTH—length of all primary roads in a buffer, SECONDARYROADLENGTH—length of all secondary roads in a buffer, TERTIARYROADLENGTH—length of all tertiary roads in a buffer, RESIDENTROADLENGTH—length of all resident roads in a buffer, DISTINVMAJOR1 and DISTINVMAJOR2—distance to the nearest primary road, DISTINVSECOND1 and DISTINVSECOND2—distance to the nearest secondary road, DISTINVTERT1 and DISTINVTERT2—distance to the nearest tertiary road, DISTINVNEAR1 and DISTINVNEAR2—distance to the nearest road. NA—Not applicable [53].

. The initial step in model development was calculating univariate regression between each potential predictor variable (Table A1) and measured PM_2.5 target variable (

Table 2. Mean seasonal PM_2.5 concentrations (µg/m³) for all sites during heating and non-heating seasons.

Site No.	Location Type	Winter PM2.5	Summer PM2.5
1	URB	25.57	19.56
2	URB	50.21	22.31
3	URB/IND	43.77	14.22
4	URB	38.19	10.85
5	URB	34.17	10.49
6	URB	32.89	10.71
7	IND	39.65	17.14
8	URB	33.55	13.72
9	URB	25.25	12.19
10	URB	35.64	20.36
11	URB	25.45	22.48
12	URB	36.81	13.73
13	URB	34.16	20.62
14	URB	27.92	16.15
15	URB	35.63	14.41
16	URB	28.37	17.82
17	BCG	27.39	12.22
18	URB	32.5	16.15
19	URB	36.12	14.81
20	URB	30	20.32
21	URB	26.2	10.47

Note. URB—urban location, IND—industrial location, URB/IND—location at the border of the urban and industrial zone, BCG—background location.

) in order to determine the strength of their correlation by observing the Pearson correlation coefficient (R) and p-value. The predictor variable with the highest R value was used in the model, and the supervised stepwise forward procedure was continued in an iterative manner. Predictor variables that were not correlated with PM_2.5 measurements, had p > 0.05, did not have all observations for the study sample (N = 21) or had a predictor direction of effect opposite to the initially defined direction of effect were rejected. Additionally, for the buffer type predictor variables, if in the set of variables (e.g., corresponding to road length variable ROADLENGTH_50m) one had the highest correlation with PM_2.5, then the variable with the smaller buffer radius was not considered since values of this variable were already included within the larger buffer (e.g., variable ROADLENGTH_25m was not taken into account since the larger buffer, variable ROADLENGTH_50m, is included). One other criterion was if, e.g., variable ROADLENGTH_100m also had a significant correlation with PM_2.5, then in the regression model, an additional ring buffer was defined as the difference between ROADLENGTH_100m and ROADLENGTH_50m. After all remaining suitable predictor variables were added one by one into the set of possible predictor variables, the change in the percentage of explained variance of PM_2.5 (adjusted R², adj. R²) was observed. During the stepwise iterative procedure for the LUR model development, the following steps were followed in order to avoid overfitting and to ensure the choice of a stable model. In each step, the best candidate predictor is included in the model only if the adjusted R² increased by more than 1%, otherwise, if there is no such predictor, the iterative procedure stops and the model is finalized. The direction of effect (i.e., sign of the coefficient corresponding to the particular predictor) for the newly included predictor must be in accordance with the predetermined sign, and furthermore, all coefficients in the new regression equation must remain the same, i.e., in accordance with their predetermined direction of effect [53]. Adding predictors to the model stops when there are no remaining predictors that can increase adj. R² by at least 1% compared to the previous regression model. As a final step, variables that are not statistically significant (p > 0.10), were eliminated from the model, while those with p ≤ 0.10 were kept in the model.

Evaluation of the final winter and summer PM_2.5 LUR models was conducted by examining Cook’s distance (Cook’s D), heteroscedasticity (White test) and normality of residual distribution.

For model validation, we used the Leave-One-Out Cross-validation (LOOCV), in which data from one site is left out, the LUR model is rebuilt based on the remaining points, the error for the left-out point is calculated, and then the total error is averaged for all observation sites. The benefit of such a leave-one-out analysis is that, in a way, all measuring points are used for both model development and model validation. The performance of the developed model was assessed by observing adj. R², RMSE and MAE values [38,54,55].

3. Results

3.1. PM_2.5 Measurements

Mean seasonal PM_2.5 concentrations during heating and non-heating seasons and their descriptive statistics are given in Table 2 and

Table 3. PM_2.5 (µg/m³) descriptive statistics for heating and non-heating seasons.

Statistics	Winter PM2.5	Summer PM2.5
Median	33.55	14.81
Mean	33.31	15.75
Std. Deviation	6.44	3.99
Minimum	25.25	10.47
Maximum	50.21	22.48
Variance	41.43	15.92

. Note that, as mentioned in the Materials and Methods section, each location had a total of five 48 h measurements per season that were used for inferring seasonal concentrations. Some details and comments about these individual measurements are stated in Appendix D.

3.2. Development and Evaluation of Seasonal LUR Models

3.2.1. Development and Evaluation of Winter LUR Model

Among all the potential predictor variables offered for the development of the winter PM_2.5 LUR model (PM_2.5_LUR_WINTER) from Table A1, variables ROADLENGTH_50m, park_area_3000m, suburban_area_3000m, class3_area_1000m and resident_area_300m had the highest R; however, class3_area_1000m and resident_area_300m were not used in further calculation due to insufficient number of observations (different from zero). Correlations between PM_2.5 and potential predictors expressed through the R values and the number of observations for these potential predictors are given in

Table 4. Correlation coefficient, R, and the number of available observations (different from zero), N, for the initial set of possible predictors for the winter LUR model.

Predictors	R	N
ROADLENGTH_50m	0.458 *	21
Class3_1000m	−0.778 *	7
park_area_3000m	−0.548 *	21
suburban_area_3000m	0.602 **	21
resident_area_300m	0.497 *	18

The distinction between moderate and strong correlations was made based on p-value using the following criteria: p < 0.05 (moderate correlation, marked by *), and p < 0.01 (strong correlation, marked by **).

.

Tested models considered within a forward stepwise regression process along with values of R², adj. R², F-test, DF, p-value, VIF and RMSE are given in

Table 5. Models tested for the winter season.

Potential Predictors	R	adj. R2	F-Test	DF	p Value	VIF	RMSE
suburban_area_3000m	0.602	0.328	10.77	1. 19	0.004	1.00	5.02
suburban_area_3000m+park_area_3000m	0.624	0.322	5.75	2. 18	0.012	2.05	4.91
suburban_area_3000m+ROADLENGTH_50m	0.770	0.548	13.12	2. 18	0.000	1.00	4.01

Note. R—coefficient of multiple correlation, adj. R²—corrected coefficient of multiple determination—percentage of explained variance, F-Test—value of the F quotient, DF—number of degrees of freedom, VIF—variance inflation factor (multicollinearity index), RMSE—root mean square error.

. Considering the values of adj. R² and RMSE, the third model in which the predictor variables are suburban_area_3000m and ROADLENGTH_50m is noticeably better than the first and second proposed models, and as such it is optimal for predicting winter PM_2.5. The values of B (unstandardized) coefficient and the p-value for optimal model final predictors are B = 0.0000019, p = 0.001 for suburban_area_3000m and B = 0.024, p < 0.01 for ROADLENGTH_50m. Unstandardized coefficients in the final regression equations are used due to different measuring scales of predictors. The VIF criteria for multicollinearity and Cook’s D criteria used for model development were <3 and <1 [56].

The regression equation for this final PM_2.5_LUR_WINTER model is given by Equation (1). By applying the final PM_2.5_LUR_WINTER regression formula, it is possible to explain about 55% of the variance in the PM_2.5 concentration.

PM_2.5_LUR_WINTER = 24.89 + suburban_area_3000m × 1.9 × 10⁻⁶ + ROADLENGTH_50m × 2.4 × 10⁻²

(1)

The evaluation process for the model gave Cook’s D observation values below 1 and White’s test was not statistically significant (χ²(5) = 1.21, p = 0.944), with a p-value higher than the level of significance (0.05) [57], showing the absence of heteroscedasticity from the tested data (Figure 3b). Values of skewness (SK = −0.077) and kurtosis coefficients (KU = −1.09) are within the acceptable range ±1 [58,59], thus the distribution of residuals for the winter model can be treated as near normal as presented on residual plots in Figure 3a–d.

Validation of the model using LOOCV is presented in Appendix B.

3.2.2. Summer LUR Model Development and Evaluation

For the development of a summer LUR model (PM_2.5_LUR_SUMMER), all potential predictor variables, their correlations with target PM_2.5 expressed through their R values and the number of observations are shown in

Table 6. Correlation coefficient, R, and the number of available observations (different from zero), N, for the initial set of possible predictors for the summer LUR model.

Predictors	R	N
ROADLENGTH_25m	0.497 *	21
ROADLENGTH_50m	0.658 **	21
ROADLENGTH_100m	0.494 *	21
ROADLENGTH_200m	0.507 *	21
ROADLENGTH_500m	0.510 *	21
PRIMARYROADLENGTH_100m	0.951 *	4
PRIMARYROADLENGTH_200m	0.981 *	4
PRIMARYROADLENGTH_300m	0.994 **	4
PRIMARYROADLENGTH_500m	0.813 *	7
SECONDARYROADLENGTH_500m	0.574 *	16
SECONDARYROADLENGTH_1000m	0.680 **	21
SECONDARYROADLENGTH_1500m	0.479 *	21
DISTINVSECOND1	0.491 *	21
DISTINVSECOND2	0.442 *	21
Class2_300m	0.956 *	4
industrial_area_300m	0.970 *	4

The distinction between moderate and strong correlations was made based on p-value using the following criteria: p < 0.05 (moderate correlation, marked by *) and p < 0.01 (strong correlation, marked by **).

. The same elimination criteria used in the winter model were applied to these potential predictors. Buffer-type predictor variables that significantly correlate with the PM_2.5 but also have a corresponding predictor with a smaller buffer radius that has a lower correlation were not included in the model. For example, the variable SECONDARYROADLENGTH_500m was not included in the analysis as a predictor because the variable SECONDARYROADLENGTH_1000m is more strongly correlated to PM_2.5 and has a larger buffer radius. When a variable of the same class has a weaker correlation with the dependent PM_2.5 variable compared to the same class with a smaller buffer radius, then a new ring buffer variable is derived. For example, from the difference of variables ROADLENGTH_300m and ROADLENGTH_50m, the variable ROADLENGTH_50_300 is considered for inclusion in the model.

Tested models considered within a forward stepwise regression process are given in

Table 7. Tested models for the summer season.

Potential Predictors	R	adj. R2	F-Test	DF	p Value	VIF	RMSE
SECONDARYROADLENGTH_1000m	0.680	0.435	16.38	1.19	0.001	1.00	2.85
SECONDARYROADLENGTH_1000m + ROADLENGTH_50m	0.767	0.542	12.83	2.18	0.000	1.38	2.50
SECONDARYROADLENGTH_1000m + ROADLENGTH_50m + ROADLENGTH_50_300	0.797	0.570	9.85	3.17	0.001	1.81	2.35
SECONDARYROADLENGTH_1000m + ROADLENGTH_50m + ROADLENGTH_50_300 + ROADLENGTH_300_500	0.812	0.574	7.73	4.16	0.001	3.89	2.27
SECONDARYROADLENGTH_1000m + ROADLENGTH_50m + ROADLENGTH50_300 + DISTINVSECOND1	0.797	0.544	6.95	4.16	0.002	1.84	2.35
SECONDARYROADLENGTH_1000m + ROADLENGTH_50m + SECONDARYROADLENGTH_1000_1500m	0.836	0.623	9.27	4.16	0.000	2.78	2.14
SECONDARYROADLENGTH_1000m + ROADLENGTH_50m + ROADLENGTH_50_300 + DISTINVSECOND_2_1	0.797	0.544	6.96	4.16	0.002	1.84	2.35

Note. R²—coefficient of multiple correlation, adj. R²—corrected coefficient of multiple determination—percentage of explained variance, F-test—value of the F quotient, DF—number of degrees of freedom, VIF—variance inflation factor (multicollinearity index), RMSE—root mean square error.

. According to the values of adj. R², adj. ΔR², RMSE and VIF coefficients, the third model, which included the predictor variables SECONDARYROADLENGTH_1000m, ROADLENGTH_50m and ROADLENGTH_50_300, was optimal. This model had lower RMSE and higher adj. R² values compared with the first two models, and multicollinearity was not present, as in the fourth model. The fifth model has a negative adj. ΔR² coefficient value (compared to the third model), while in the sixth and seventh models, there are B coefficient values that are not in accordance with the direction of effect. By applying the final step in model development, that is, statistical significance criteria, variables ROADLENGTH_50_300 and SECONDARYROADLENGTH_1000m were removed from the PM_2.5_LUR_SUMMER model because their p > 0.05. Thus, in the third optimal model, ROADLENGTH_50m is the only significant predictor of the positive direction, with a B coefficient of 0.020 and p < 0.01.

This final PM_2.5_LUR_SUMMER model was shown to be statistically significant (F-test = 14.49, p = 0.001), with the regression equation represented by Equation (2). By applying this regression equation, the predictor explained 40.3% of the variation in the summer PM_2.5 concentration.

PM_2.5_LUR_SUMMER = 12.43 + ROADLENGTH_50m × 2.0 × 10⁻²

(2)

Evaluation of the PM_2.5_LUR_SUMMER model also showed that there were no observations with Cook’s distances greater than 1.00, that White’s test was not statistically significant (χ² (2) = 3.28, p = 0.198) (Figure 4b) and that values of skewness (SK = 0.853) and kurtosis (KU = −0.803) coefficients were within the acceptable range of ±1 [58,59], implying near normal residual distribution as can be seen on the residual plots (Figure 4a–d).

Validation of the model using LOOCV is described in Appendix B.

3.2.3. Simple Seasonal PM_2.5 Model Based on Aerosol Optical Depth Data

In this section, we show a simple and direct linear regression approach for connecting aerosol optical depth (AOD) data for the city of Novi Sad and PM_2.5 ground-level data.

In recent years, successful efforts have been made to utilize satellite-based aerosol information like observations of Aerosol Optical Depth (AOD) to improve the spatial mapping and prediction of ground-level PM_2.5. The relationship between AOD and PM is complex and depends on aerosol type and size, relative humidity, planetary boundary layer height and the vertical structure of aerosol distribution. Statistical approaches including single and multiple linear regression and more complex machine learning approaches for estimating PM based on satellite data have demonstrated some advantages compared to more traditional physical retrieval methods [60]. Various satellite products including Copernicus Atmosphere Monitoring Service (CAMS), Global Near-Real-Time provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) [61] and MCD19A2.061: Terra & Aqua MAIAC Land Aerosol Optical Depth provided by NASA LP DAAC at the USGS EROS Center [62], can be used to infer ground level PM_2.5 concentrations using aerosol optical depth measurements originating from satellites. Such estimates of PM_2.5 can be conducted using statistical approaches, which can range in complexity from linear regression models (e.g., generalized linear model [63] multiple linear regression [64]) to more complex models based on machine learning approaches [65].

Google Earth Engine was used to access and process the ECMWF/CAMS and MODIS datasets. The spatial resolution of the CAMS data is approximately 0.1 degrees (~10 km at the equator), while the MODIS AOD data from the MCD19A2 product has a native spatial resolution of 1 km. The process involved filtering the data by date, selecting relevant bands for AOD, averaging the data over the specified periods and resampling the data using bilinear interpolation. Figure 5 shows a scatterplot and regression line for seasonal models that use simple univariate linear regression with no intercept to quantify the relation between AOD data and seasonal PM_2.5 ground-level data (which was also used for the seasonal LUR models presented in this paper). Validation of the model using the LOOCV approach is described in Appendix C.

4. Discussion

In most cities in Serbia, PM concentrations show drastic differences during heating and non-heating seasons [66], and the origin and simple quantification of this difference in the city of Novi Sad can be seen using results from reference monitoring and emission sources [67]; however, these lack spatial resolution. This difference was also analysed using a mobile monitoring approach (this increases spatial resolution) for PM₁₀, PM_2.5 and ultrafine aerosol fraction, and was performed using research-grade instruments, TSI Nanoscan 3910 and TSI OPS 3330 [68,69]. These underpinning reports and research studies gave motivation for the development of the seasonal LUR model, which used reference instruments and enables quantification of differences between seasonal levels of air pollution, and consequently air pollution exposure and association between seasonal long-term exposure and possible health risk implications over the whole city of Novi Sad.

The initial set of predictors for the development of the Novi Sad winter PM_2.5 LUR model (Table 4) had a mixture of source and sink predictors. The initial set contained the following predictors that model PM_2.5 sources: ROADLENGTH_50m, which is a proxy for traffic intensity on nearby roads; suburban_area_3000m, which can be used as a proxy for individual heating sources in a site’s surrounding area, and similarly, resident_area_300m. The initial set contained the following PM_2.5 sinks: Class3_1000m, which is the number of forests and seminatural areas in a 1000 m buffer that models the amount of aerosol deposition, and similarly, park_area_3000m. Target PM_2.5 seasonal concentration was inferred based on the measuring campaign that was conducted during the heating season months from December to March. These raw measurements were corrected to obtain seasonal averages using the ESCAPE procedure [52].

The iterative development procedure resulted in the three models shown in Table 5. It is interesting to note that the best model in terms of explained variance contains only source predictors (suburban_area_3000m, ROADLENGTH_50m), while the second-best model contains a source (suburban_area_3000m) and a sink (park_area_3000m) predictor. While it is evident that the seasonal concentration depends on a complex interplay of source and sink processes, the extent of available data and applied statistical procedure for the model development favoured the source-only model.

The predictor group suburban area can be used as a proxy of air pollution sources related predominantly to household areas that are situated next to the urban city centre and DH area and quantifies the contribution of IH during the winter. The predictor that explained the most variance in PM_2.5 concentration during the winter season in Novi Sad was suburban area in a buffer of 3000m, with adj. R² = 0.328. Buffers with smaller radii did not show statistical significance when explaining PM_2.5 pollution. The second significant predictor variable, total road length within a buffer of 50 m, is a very rough proxy for variations in traffic density in the close vicinity, adding to the total explained variance of adj. R² = 0.548.

Generally, each developed LUR model contains unique definitions of potential independent variables, depending on the study area specifications, available data, and pollutants of interest [70]. In cases where data on traffic intensity are not available or do not exist and thus cannot be linked to the local and central road bases, road classification and their lengths within the research area can be used as proxies for the potential predictor variables [23,36,71]. Similarly, traffic emissions can be used as a proxy for major arterials, railroads or street lengths in cases where these data do not exist [72]. Seasonal LUR models that model higher winter concentrations of PM_2.5 and the potential influence of the heating combustion process were developed for several cities in China [35,73,74]. A winter model developed by Shi et al. [35] explained 61% of winter PM_2.5 variability, and predictors in the final LUR model were average building floor areas within 3000 m that explained the most variation, followed by water area within 3000 m and average building heights within 50 m (predictors with negative coefficients). Other types of predictors such as meteorological factors, roads, land use types, various points of interest for better temporal resolution and elevation were used for LUR model development in Beijing [73] and explained between 54–86% of PM_2.5 winter variability. The highest percentage of explained PM_2.5 variability among the mentioned heating season models was 94%, with the explanatory variables precipitation, density of roads, elevation, average air pressure and vegetation [74]. The LUR model created by Saucy et al. [75] explained 29% of the variance with the use of population/building density within 25 and 300 m buffers, construction sites within a 100 m buffer, length of bus routes in a 300 m buffer and distance to the nearest waste burning sites as explanatory variables. Road lengths, as significant predictors within smaller buffer zones, were also reported in winter LUR models developed for Calgary, Canada [37], where the total length of major roads within a buffer of 200 m was significant, noting that in that study buffer zones started from 100 m due to local city fabric. The value of adj. R² for this winter model was 45%, similar to [36]. The best-explained variability of PM_2.5 in the winter season was 90%, reported in [38], which used six independent predictor variables.

The initial step of the model development produced a list of predictors for the Novi Sad summer PM_2.5 LUR model (Table 6) that had only source predictors serving as proxies for traffic intensity. Target PM_2.5 seasonal concentration was inferred based on a measuring campaign that was conducted during the non-heating season months of July–September. These raw measurements were corrected to obtain seasonal averages using the ESCAPE procedure [52].

The iterative development procedure resulted in the seven models shown in Table 7. The procedure started with a univariate model with SECONDARYROADLENGTH_1000m as a predictor. This predictor was calculated over a large buffer area and varied with the number of secondary roads in the part of the city encircled by the buffer. The iterative procedure then continued and arrived at the bivariate model with SECONDARYROADLENGTH_1000m and ROADLENGTH_50m (number of roads in a very small buffer around the site) as predictors. Note that so far, all models derived from the stepwise procedure have statistically significant predictors; however, also note that VIF is increasing. Finally, following the stepwise ESCAPE procedure, we arrived at the trivariate model with SECONDARYROADLENGTH_1000m, ROADLENGTH_50m and ROADLENGTH_50_300 (number of roads in a circular ring with inner and outer diameters of 50 and 300 m, respectively) and explains a total of adj. R² = 0.570. Here, we stop the ESCAPE procedure since no new predictors can be added that increase adj. R². Note also that the VIF parameter is further increased, indicating a mild collinearity between the predictors. It is known that such collinearity between the predictors does not reduce the explained variance; however, it does reduce the statistical significance of the predictors. The ESCAPE procedure mitigates this by omitting, in the final step of the model development, predictors that are not statistically significant. For purposes of completeness, we have also, in Appendix B, added validation of the first two models in the stepwise procedure, since these were fully statistically significant but were not reached by direct application of the ESCAPE procedure.

In the summer seasonal model developed for the city of Novi Sad, ROADLENGTH_50m was the only significant predictor variable included in the model, with an explained variance of PM_2.5 concentration of 40%. The road length predictor within a buffer of 50 m had a higher value of R = 0.658 than in the winter. This can be explained by the fact that small-scale traffic variations are more easily observed due to fewer variations in emission sources of PM_2.5 in urban environments in the summer season. Few existing studies have efficiently used road length without linked traffic, but good road classification is needed for their application as potential predictor variables [70,71].

In comparison to our PM_2.5_LUR_SUMMER model and achieved percentage of explained variance, models in [32,37,38] explained 72%, 50% and 59% of PM_2.5 variance, respectively. In all three studies, predictor variable total lengths of local roads were included in the final LUR models. The summer model developed by Bertazzon et al. [37] included the sum of several road segment lengths in a buffer of 750 m as a final predictor. Contrary to previously mentioned models and our LUR_SUMMER model, the model developed by Jeffrey [72] used explanatory variables based on road classes and explained 97% of PM summer variability using additional predictor variables, including temperature, average air pressure, altitude, humidity and precipitation. The winter model’s final predictors developed by [35] were the same as for their summer season model, which explained 52% of summer variability. The lowest explained PM_2.5 variability for summer among reviewed articles was 36% and the explanatory variables included in the final model were the length of railways in a 1000 m buffer, distance to the nearest open grills, population/building density, waste burning sites and distance to nearest refuse transfer station [75].

Both the winter and summer LUR models outperform the AOD satellite data-driven approach, as evidenced by the consistently lower RMSE values for the winter and summer datasets (see Appendix B and Appendix C). This suggests that incorporating LUR predictors enables better capturing of local-scale variations in PM_2.5 concentrations compared to AOD satellite-derived data. While satellite observations can offer broader coverage, their limited ability to resolve fine-scale variations results in higher RMSE values, highlighting the importance of integrating LUR-specific predictors for more accurate air quality assessments.

5. Conclusions

The developed PM_2.5_LUR_WINTER and PM_2.5_LUR_SUMMER models are seasonal land use regression models. These models represent an advancement in the field of land use regression (LUR) for PM_2.5 in the region.

The PM_2.5_LUR_WINTER model includes two predictors; the first one can be associated with individual heating sources in the city. Along with the suburb_3000m predictor as a proxy for combustion processes in individual heating sources, predictor ROADLENGTH_50m stood out as the predictor of traffic originating PM_2.5 in the winter model, quantifying the influence on PM_2.5 particulate matter emission during the winter due to close vicinity roads. The PM_2.5_LUR_SUMMER model emphasizes the influence of traffic as the most important emission source in the summer, though only one predictor ROADLENGTH_50m was statistically significant following the final step of the ESCAPE procedure. The lack of street-level traffic data probably affected the achieved percentage of variability explained by both models and including currently unavailable street-level traffic model would be the straightforward path for improvement of the model.

The results obtained in this study could be of direct use in the assessment of air pollution-related health risks in the municipality and contribute scientifically and practically to local authorities dealing with air quality control. Since models showed good predictive accuracy, developed seasonal PM_2.5 LUR models for Novi Sad can be used in information products as base maps of fine particulate matter pollution.

As part of further research opportunities, we plan to refine developed models through additional measurement campaigns, utilizing low-cost sensor data and adding more exact traffic and heating predictors for emissions instead of proxies, thus increasing the interpretability of the model. Low-cost sensor data can be used to increase the spatial and temporal resolution of the seasonal models by fusing measurement data with modelled pollution data. Another possible direction of future research that may bring improvements to heating and non-heating season model predictive efficiency is the inclusion of machine learning algorithms and more complex models utilizing satellite data since a simple aerosol optical depth model gave promising results in this study.