## 1. Introduction

Clean air is considered to be a basic requirement for human health and well being. Various chemicals are emitted into the air from both, natural and anthropogenic sources. In spite of the introduction of cleaner technologies in industry, energy production and transport, air pollution remains a major health risk and tighter emission controls are being enforced by many governments. Atmospheric particles – aerosols – are some of the key components of the atmosphere. They influence the energy balance of the Earth's surface, visibility, climate and environment as a whole [1-3]. According to World Health Organization (WHO), ozone, particulate matter(PM), heavy metals and some hydrocarbons present the priority pollutants in the troposphere [4].Public health can also be indirectly affected by deposition of air pollutants in environmental media and uptake by plants and animals, what results in entering of chemicals into the food chain or drinking water, and thereby constituting additional sources of human exposure. A number of epidemiological studies have demonstrated that acute and chronic health effects are related to the inhalable PM_{10} (aerodynamic diameter less than 10 μm) exposure in the urban environment, and some data also seem to indicate possible seasonal effects of the particulate matter on human health [5-10]. This is especially important for urban aerosols, whose variety of size and composition make complete characterization a difficult task. Particulate matter pollution is nowadays one of the problems of the most concern in great cities, not only because of the adverse health effects, but also of the reducing atmospheric visibility and affect to the state of conservation of various cultural heritages [11]. Therefore, the measurement of the levels of atmospheric particulate matter is a key parameter in air quality monitoring throughout the world.

As a result of health and environmental impacts of PM, more rigorous regulations are in force in the USA and European countries. PM standards, issued by European Commission (EC), have included PM_{10} monitoring and limit values in the Air Quality Directive in 1999 [12]. Directive established in the first stage, annual limit value of 40 μg m^{-3} and 24h limit value of 50 μg m^{-3} (not to be exceeded more than 35 times in a calendar year) to be met by 2005, and in the second stage annual limit value of 20 μg m^{-3} and 24h limit value of 50 μg m^{-3} (not to be exceeded more than 7 times in a calendar year) to be met by 2010.The detaileddiscussion of these limit values, regulations and relations of new EU standards to US EPA standards can be found elsewhere [13,14].

Most of the trace metals are emitted in particulate form [15-17] and present in almost all aerosol size fractions, but are mainly accumulated in the smaller particles [18]. This has a great effect on the toxicity of metals since the degree of respiratory penetration depends on particle size. Trace metals are persistent and widely dispersed in the environment and interacting with different natural components result in toxic effects on the biosphere [19,20]. The studies of the transport and mobilization of trace metals up to now have attracted attention of many researchers [21-24]. Trace metals are released into the atmosphere by human activities, such as combustion of fossil fuels and wood, high temperature industrial activities and waste incinerations. The combustion of fossil fuels constitutes the principal anthropogenic source for Be, V, Co, Ni, Se, Mo, Sn, Sb, and Hg. It also contributes to anthropogenic release of Cr, Mn, Cu, Zn, and As. High percentages of Ni, Cu, Zn, As, and Cd are emitted from industrial metallurgical processes. Exhaust emissions from gasoline may contain variable quantities of Ni, Cu, Zn, Cd, and Pb [25,26]. In addition to the PM mass concentration limit values, also based on health impact criteria, recent European Union standards set target (Ni, As, Cd) and limit (Pb) values for metals [12]. Environmental technologies may have to be adopted in specific industrial spots to reach the target values. For aimed reduction of PM_{10} levels and their trace metal contents detailed knowledge of sources and their respective contribution to the PM_{10} levels, is required.

One of the main difficulties in air pollution management is to determine the quantitative relationship between ambient air quality and pollutant sources. Various receptor models have been used to identify aerosol sources and estimate their contributions to PM_{10}and trace elements concentrations at receptor sites [27-32]. The combined application of different types of receptor models could possibly solve the limitations of the individual models, by constructing a more robust solution based on their strengths.

Processes in the atmosphere represent a complex problem due to the simultaneous influence of several independent factors. Similarly to other air pollutants, PM_{10} concentrations are random variables influenced by the emission level, meteorological conditions and topography. Each area is specific and the required emission reduction to meet Air Quality Standard (AQS) is different. Information about the frequency distribution of pollutants is useful for developing air pollution control strategies. When the specific probability function of an air pollutant is known, it is easy to predict the required emission reduction, the frequency of exceedance of the AQS, as well as the return period [33]. Many types of probability distributions have been used to fit air pollutant concentrations. These include lognormal distribution [34,35], Weibull distribution [36], gamma distribution [37] and type V Pearson distribution [38]. The Weibull probability model is suitable for presenting air quality data in a “pseudo-lognormal” distribution [39,40]. In the present study, daily average mass concentrations of PM_{10} were taken from 2003 to 2005 in order to estimate the parameters for three theoretical distributions. Lognormal, Weibull and type V Pearson distributions were chosen to fit the PM_{10} data in Belgrade. The distributional parameters were estimated by the maximum likelihood method. Based on the fitted distributions the minimum reduction required in order to achieve AQS was estimated using rollback equation [41-43]. Moreover, since the tail of theoretical distributions diverged in the high concentration region, a two-parameter exponential distribution was used to fit high PM_{10} concentrations and estimate exceedances and return period more precisely [44]. These methods can reasonably predict the return period and exceedances in the succeeding period and can also be used for predicting source emission reduction.

In many countries, the weekly variation in the concentration of atmospheric aerosols and meteorological variables in urban areas has been reported as evidence of anthropogenic influence on the weather system. Of the primary sources road transport is the dominant contributor to high particulate level in the air in most major cities. This study examines the temporal and weekly variations in meteorological variables and PM_{10} mass concentrations in Belgrade. The cumulative frequencies and quantile-quantile plot for Sundays against workdays were constructed. The results showed that Sunday concentrations were significantly reduced from those during the week.

The Mediterranean region and particularly the Balkan Peninsula have been under the influence of Saharan dust transport and deposition over millennia. Identification of the concentration, composition, origin, transport and geographical distribution of PM_{10} in Mediterranean atmosphere has been the subject of research activities since the last two decades as it is heavily affected by two contrasting sources; mineral dust (mainly from Sahara Desert) [45] and various anthropogenic (from industrialized/semi-industrialized countries) emissions.Deviation from lognormal distribution of atmospheric aerosols is a positive indicator of possible transport processes [33]. Since we are particularly interested in transportation of trace metals content in PM_{10} samples, the quantiles forsome trace elements were analyzed by the quantile-quantile P-P slope test. This was used for fitting the quantiles of the expected theoretical cumulative probability for normal, lognormal and Weibull distribution with the quantiles of the cumulative probabilities calculated for the experimental dataset [46].

In addition, we investigated thetransport pathways and potential sources of PM_{10}concentrations based on backward trajectories and PM_{10} concentration records from 2004 to 2006. Cluster analysis was used to reveal the major pathways for different seasons as well as corresponding statistical analysis related to different clusters. Hybrid receptor models Potential Source Contribution Function (PSCF) and Concentration Weighted Trajectory (CWT) were used for identification of source regions. The PSCF values can be interpreted as a conditional probability describing the spatial distribution of probable geographical source locations by analysis of trajectories arriving at the sampling site [47]. Since the PSCF method is known to have difficulties in distinguishing strong from moderate sources, the CWT model that determines the relative significance of potential sources has been additionally performed[29,48].

In this review, we report some of the results of the integral monitoring of air quality in Belgrade urban area in order to evaluate the impact of airborne trace metals on the pollution load for the period from 2003 to 2006. Some of the results concerning suspended particle mass and trace metal concentrations in ambient air of the Belgrade will be presented, with the aim to examine elemental associations and to indicate the main sources of trace and other metals in the city. The results of this long-term project of the pollution monitoring could be used as the baseline data for analysis of health risks due to inhalation of suspended aerosols, and to provide scientific evidence for setting up an air pollution control strategy. This information is crucial in environmental quality assessment, and can lead to the determination of a possible exceedance of the critical loads.

## 2. Experimental methods and procedures

Sampling of particulate matter PM_{10} started in the very urban area of Belgrade (Hs = 117 m,

The second sampling point was on the roof of the Rector’s Office building of Belgrade University, at a height of about 20 m. As this sampling point is in the very city centre, on the rooftop where the airflow is not blocked by any direction, it can be considered as representative for urban-background concentrations. Suspended particles were collected on pre-conditioned and pre-weighed Pure Teflon filters (Whatman, 47 mm diameter, 2 µm pore size) and Teflon-coated Quartz filters using MiniVol air sampler (Airmetrics Co. Inc., 5 l min^{-1} flow rate) provided with PM_{10} cutoff inlet. Particulate matter mass concentration was determined by weighting of the filters using a semi-micro balance (Sartorius, R 160P), with a minimum resolution of 0.01 mg. Loaded and unloaded filters (stored in Petri dishes) were weighed after 48 hours conditioning in a desiccator, in the clean room at a relative humidity of (45-55)% and a temperature of (20 ± 2) ^{0}C. Quality assurance was provided by simultaneous measurements of a set of three ‘‘weigh blank’’ filters that were interspersed within the pre- and post- weighing sessions of each set of sample filters and the mean change in “weigh blank” filter mass between weighing sessions was used to correct the sample filter mass changes. After completion of gravimetric analysis, PM samples were digested in 0.1 N HNO_{3} on an ultrasonic bath. An extraction procedure with dilute acid was used for the evaluation of elements which can become labile depending on the acidity of the environment [49]. This procedure gives valid information on the extractability of elements, since the soluble components in an aerosol are normally dissolved by contact with water or acidic solution in the actual environment.

The elemental composition (Al, V, Cr, Mn, Fe, Ni, Cu, Zn, Cd, and Pb) of the aerosol samples was measured by the atomic absorption spectroscopy method (AAS). Depending on concentration levels, samples were analyzed for a set of elements by flame (FAAS) (Perkin Elmer AA 200) and graphite furnace atomic absorption spectrometry (GFAAS) using the transversely-heated graphite atomizer (THGA; Perkin Elmer AA 600) with Zeeman-effect background correction. The THGA provided a uniform temperature distribution over the entire tube length, rapid heating and an integrated L’vov platform, which gave an improved signal/interference ratio and high analytical sensitivity. Analyte injection (20 µl) and the atomization were done in five steps controlled by the appropriate software and auto-sampler. For calibration, standard solutions containing all metals of interest were prepared using Merck certified atomic absorption stock standard solutions containing 1000 mg l^{-1} metal in 0.5 N HNO_{3} and Milli-Q quality deionised water, with no matrix modifier addition. Details on sampling procedures and PM analysis are given in detail elsewhere [50-55].

## 3. Estimation of emission source reduction

Assuming unchanged spatial distribution of emission sources, meteorological conditions and noreactive species, according to the rollback equation [36], the emission source reduction R (%) required to meet AQS can be estimated as:

In this equation E{C}_{s} is the mean (expected) concentration of the distribution when the extreme value equals C_{s} (the concentration of the AQS), E{C_{p}} is the mean concentration of the actual distribution and C_{b} the background concentration. If C_{s} = 125 µg m^{-3}, the PM_{10} daily average concentration is not exceeded more than once per year (P[PM_{10}>C_{s}] = 1/365 = 0.00274), then E{C}_{s} is the expected daily PM_{10} average concentration of a distribution where the probability of a concentration exceeding 125 µg m^{-3}is equal to 0.00274.

Frequency distribution of PM_{10}is a very usefultool to estimate how frequentlya critical concentrationlevel is exceeded. In order to obtain the best frequency distribution of PM_{10}, three theoretical functions were used: lognormal, Weibull and type V Pearson distributions. The distribution parameters were estimated by the maximum likelihood method. Probability density functions as well as corresponding equations for estimating distribution parameters are presented in Table 1.

The appropriateness of each distribution was assessed by the Kolmogorov-Smirnov (K-S) test. The K–S statistic is defined as the maximum difference between the sample cumulative probability and the expected cumulative probability i.e.

where _{α/2} acceptable for the K-S test at a certain significance level α. If D is less than D_{α/2}, one can conclude that the hypothesis that there is no difference between observed and expected distributions has been accepted at the α level of significance. Generally speaking, if D value declines the goodness of fit increases.

The fitted results of three theoretic distributions and the measured data for PM_{10} are presented in Figure 1.Once the distribution parameters were determined, K-S test could be used to assess which type of distribution is more appropriate for representing the PM_{10} distribution. From Figure 1 it is clear that the Weibull distribution is inappropriate for representing PM_{10} distribution (D = 0.073), while type V Pearson distribution is the most suitable one (D = 0.038).

Figure 2 shows the relation between exceeding probability and PM_{10}concentrations for different distribution functions. It was found that the probabilities of exceeding the AQS (125 μg m^{-3}) were 0.12, 0.096 and 0.088 for lognormal, Weibull and type V Pearson distribution, respectively. This means that the number of days exceeding the AQS in the following year would be 43, 35 and 32, respectively. The actual probability for this period was 0.075 (27 days) and it is clear that all three distributions overestimated the exceeding probability although type V Pearson distribution most closely represented the true PM_{10} data.

After determining the most appropriate distribution function for PM_{10}, the emission source reduction required to meet AQS can be predicted from a rollback equation.The complementary distribution function for calculating the probability of variable*x*_{i}exceeding the critical value *x* for lognormal and type V Pearson distributions can be found in a statistical textbooks.It was found that the value of E{C}_{s} (the average concentration of lognormal distribution where the probability of a concentration exceeding 125 μg m^{-3}
equals 0.00274) was 32.1 μg m^{-3}. Therefore, the mean PM_{10} concentration should be reduced from the observed value of 68.4 μg m^{-3} to 32.1 μg m^{-3} and the estimated emission reduction can be calculated from equation 1 being equal to 53.1%. For type V Pearson distribution the estimated source reduction required to meet AQS is 58.6%. Therefore, the source emissions should be controlled much more to reduce PM_{10}concentrations and meet AQS in the future period [41].

Previous theoretical distributions can give good result for estimating the mean concentration and required reduction of source emission. However, the fitted results of the parent distributions are not accurate enough in the high concentrations region. Therefore, a two-parameter exponential distribution was applied to predict the return period and exceedances of the critical PM_{10} concentration.

A two-parameter exponential distribution derived from extreme value theory [44] represents the cumulative frequency distribution of high concentrations over a specific percentile.

where _{10} concentration exceeding the specific percentile. From the complete data set of PM_{10} (2003 - 2005), concentrations exceeding 80-th percentile were chosen to fit the two-parameter exponential distribution. The estimated cumulative probability can be calculated from the chosen high PM_{10} concentration,

where N_{1} is the size of the chosen high PM_{10} concentration and _{1} values. The relation between variate

and the parameters*b*_{m} and *x*_{c}, can be calculated as:

where *f* is the chosen specific percentile. After determining return period the days of exceedances within one year can be calculated in order to develop air control strategy.

Figure 3 shows the fitted theoretical line of the variate *y*_{n}and PM_{10} concentration over 80-th percentile, *x*_{n.} The fitted equation for the theoretical two-parameter exponential distribution is *x*_{c} = 125 μg m^{-3}, F_{L} and the return period can be calculated as:_{10} is 125 μg m^{-3} and the allowed exceeding probability is once per year, then the expected return period is 365 days and F_{L}(125 μg m^{-3}) = 0.986. In addition, the PM_{10} concentration expected to be equaled or exceeded once per year can be calculated from equation (3). The expected concentration is 337 μg m^{-3} and it can be seen that a reduction of 212 μg m^{-3}(about 63%) needs to be achieved in Belgrade in order to meet AQS (125 μg m^{-3}).

Although this technique can predict the probability of exceedance of the AQS and return period for the highest concentrations level it does not allow reproducing the time pattern of the estimated concentration levels. Based on the observed time series data set an empirical modelling approach can be applied in order to derive mathematical formulation suitable for the description of the feature of the concentration event [56-58].

## 4. Estimating the frequency distribution of PM_{10} by the statistics of wind speed

Since the distributions of air pollutant concentrationsare influenced by meteorological conditions, such as wind speed, Simpson et al. [59] linked airquality data and wind speed to the AtmosphericTurbulence and Diffusion Laboratory (ATDL)model.The concentration of air pollutant data at cumulativeprobability, *p*, is inversely proportional to thewind speeds *v*at probability (100-p) when thefrequency distributions of both data are log-normallydistributed and the geometric standard deviation(shape factor) is the same e.i.

where K depends on the average area source strength of pollutant and atmospheric stability factor (which varies with atmospheric conditions). This model was used to predict the maximum level of air pollutant. Finally, the frequencydistribution of air pollutants can be estimated from thefrequency distribution of wind speed [60]. Here we demonstrate the relationship between the frequencydistributions of wind speed and daily PM_{10} measured in Belgrade during 2003 year. Threetheoretical frequency distributions, (lognormal,Weibull and gamma) were used tofit the measured PM_{10} data.The distributional parameterswere estimated by the method of maximum-likelihood and the K–S test was used to determine whichtype of distribution is appropriate to represent thefrequency distribution of PM_{10}and windspeed. The fitted results of PM_{10} are shown in Figure 4 and it can be seen that the lognormal distribution is the most appropriate one to represent the statistical character of PM_{10} (D_{max} = 0.06). The fitted results of lognormal distributions for PM_{10} and wind speed, with corresponding parameters, are presented in Figure 5. If is found that the shape factors of lognormal distributions σ_{g} are similar, therefore equation 8 can be used to estimate the distribution of PM_{10} from the distribution of wind speed. The average K value can be calculated from the regression line using method of least-squares for the wind speed data and PM_{10} concentrations from the 30-th to 80-th percentiles (the K value is almost constant in that range). The obtain regression is

with the correlation coefficient R = 0.99 and estimated K = 141.5 µg m^{-2} s^{-1}.

When the atmospheric condition is unchanged and the K value is estimated, the concentration of PM_{10} at different percentiles can be obtained. Observed and estimated frequency distribution of PM_{10} is shown in Figure 6. The difference between the observed and estimated concentration is significant in the higher cumulative probability range, because the inconsistency of the both, wind speed and PM_{10} data with lognormal distribution in this region. The result of K-S test shows that the difference between observed and estimated distribution is not significant at the 5% significance level, so the obtain results are acceptable.

It is clear that there is no universal distribution that can be used to represent the frequency distribution of PM_{10}concentrations. This analysis shows that if the distribution type and shape factors are the same for the wind speed data and air pollutant concentration, there exists simple relationship between two data sets and the frequency distribution of the air pollutant can be reasonably estimated. It can be useful for understanding statistical characteristics of air quality in different region.

## 5. Weekend effect and PM_{10} temporal variations

Weekly cycles of the mass concentrations of anthropogenic aerosols have been observed in many regions and the intensity of these cycles vary significantly depending on region and season [61,62]. In addition, ozone concentrations tend to be higher on weekends than on weekdays, despite the fact that lower emissions of ozone precursors are expected on weekends.This phenomenon is known as the “weekend effect” [63]. The mechanisms for weekend effect on ozone formation are still not well understood. In a lot of papers that are dealing with weekend effect in detail, the diminution of NOx, CO, PM_{10} and volatile organic compound (VOCs) concentrations can be noticed during the weekend [64,65]. Although the relationship between the variation in aerosol concentration and weather variables has not been clearly understood, many studies have revealed significant weekday–weekend differences in weather variables such as the diurnal temperature range [66], and clouds and precipitation [67,68].

Generally speaking, several source types make significant contribution to particle level in Belgrade: combustion processes, mixes road dust, traffic, metal processing and transport [30,50]. Like in many other major cities with a few millions inhabitants, of the primary sources, fossil fuel combustion and traffic contributions play the most significant role [51].

Since traffic and commercial pattern in urban area show a marked change at weekends, it is expected that particulate concentrations will decrease. Markovic et al. [69] reported the mean daily concentrations of SO_{2}, NO_{2}, PM_{10}, CO and O_{3} during the weekdays and weekend in November/December 2005 case study in Belgrade urban area (Figure 7). The results explicitly show that daily average concentrations of SO_{2}, NO_{2}, PM_{10} and CO show significant weekly cycles with the largest values around midweek and smallest values in weekend. With the exception of ozone, concentrations of airborne pollutants were higher during weekdays. It has been suggested that the weekend condition of lower particle concentrations resulted in less absorption of sunlight, which in turn leaded to enhanced photochemical reaction and enhanced ozone formation [70]. Besides, the higher concentrations of PM_{10}, NO_{x}, CO and SO_{2} during the weekdays influence on the increased deposition of O_{3} on the particles as well as on the increased consummation of O_{3} in the reactions with NOx, CO and SO_{2} what also can have an impact.

In order to investigate differences in PM_{10} mass concentrations between Sundays and weekdays quantile-quantile (Q-Q) plots were constructed [71]. Q-Q plots are useful in detection the distribution differences in data set and allow comparison of concentration distribution of the entire range of values, from the smallest to the largest. If two cumulative frequency diagrams (*x* and *y*) are plotted, then, corresponding to any frequency*p*, there are two quantile values *p*. Obviously, if *x* and *y* are identically distributed, the plot of

The cumulative frequency diagram and quantile-quantile plot for Sundays against weekdays (weekdays were defined as Monday to Friday) were constructed and are shown in Figures 8 and 9. It can be seen that Sunday concentrations were reduced from those during the weekdays. These figures are based on complete PM_{10} data set for the period 2003-2005 in Belgrade urban area continuously recorded by the Public Health Institute of Belgrade.

We have analysed mean diurnal variations of PM_{10} mass concentrations and the characteristic pattern based on data set for the year 2005 is shown in Figure 10. It can be seen that the maximum concentrations occur in the morning and late evening which can be significantly related to traffic pattern and anthropogenic influence. Significant seasonal variations of some metals (Ni, V, Al and Fe) have already been reported [30] indicating the impact of fossil fuel combustion and local heating units. Monthly variation of PM_{10} mass concentrations for the year 2005 is presented in Figure 11 showing significant differences during the winter and summer period.

Analysis of air quality data in the frequency domain contributes to the understanding of periodic behaviours and yields information about temporal and spatial scales of the underlying mechanisms [72]. The results obtained in many studies show that the frequency analysis of air quality time series is a powerful technique that can provide important information about the nature of the processes behind the measurements. This methodology could be applied to determine representative background concentrations of air pollutants by removing short-term fluctuations associated with influence of local emission sources [73].

## 6. Transport analysis

Processes in the atmosphere represent a complex problem due to the simultaneous influence of several independent factors (meteorological conditions, pollutant emission level and topography) and thus the environmental data are random variables that follow natural lognormal distribution [33]. Deviation from lognormal distribution of atmospheric aerosols is a positive indicator of possible transport processes. In this study the quantiles for some trace metals in PM_{10} samples were analyzed by the quantile-quantile P-P slope test. It is used for fitting the quantiles of the expected theoretical cumulative probability for normal, lognormal and Weibull distribution with the quantiles of the cumulative probabilities calculated for the experimental dataset.

Fittings of quantiles of the experimental databasewith the expected quantiles for theoretical normal,lognormal and Weibull distributions for Al, Mn and Vare shown in Figure 12.

The P-P slope test results show that the Weibull distribution is the most appropriate for V and Ni suggesting their possible regional transport. The other metals followed lognormal distribution (except Cu and Zn data which are not successfully described by any of theoretical distribution function) what is a positive signal that an active sources exist in the vicinity.For unknown metal sources, hybrid models that incorporate wind trajectories Potential Source Contribution Function (PSCF) and Concentration Weighted Trajectory (CWT) were used to resolve source locations [29].

### 6.1.Potential. Source Contribution Function (PSCF)

Air parcel back trajectories, ending at the receptor site, are represented by segment endpoints. Each endpoint has two coordinates (latitude, longitude) representing the central location of an air parcel at a particulate time. To calculate PSCF, the whole geographic region of interest is divided into an array of grid cells whose size is dependent on the geographical scale of the problem so that PSCF will be a function of locations as defined by the cell indices *i* and *j*. If a trajectory end point lies at a cell of address (*i, j*), the trajectory is assumed to collect material emitted in the cell. Once aerosol is incorporated into the air parcel, it can be transported along the trajectory to the receptor site. The PSCF value can be understood as a conditional probability that describes the spatial distribution of probable source locations. The staying time of all trajectories in a single grid cell is *n*_{ij}, and *m*_{ij} is the staying time in the same cell that corresponds to the trajectories that arrived at a receptor site with pollutant concentrations higher than a pre-specified criterion value (average PM_{10} and metal concentrations were used in this study). The PSCF value for the *ij-*th cell is then defined as

Cells related to the high values of potential source contribution function are the potential source areas. However, the potential source contribution function maps do not provide an emission inventory of a pollutant but rather show those source areas which emissions can be transported to the measurement site.To remove the large uncertainty caused when a grid cell has small staying time and large PSCF values, usually a weight function W(*n*_{ij}) multiplies the PSCF value to better reflect the uncertainty in the values for these cells [27]. In this study the grid covers area of interest defined by (35^{0}- 50^{0})N and (10^{0}- 30^{0})Ewith cells 0.5^{0}x0.5^{0} latitude and longitude.

Air masses back trajectories were computed by the HYSPLIT (HYbrid Single Particle Lagrangian Integrated Trajectory) model [74] through interactive READY system [75]. Daily 48h back trajectories, started from the center of Belgrade (44.804^{0}N 20.478^{0}E) at 12:00 UTC each day, were evaluated for six different heights above the starting point at ground level (350, 500, 750, 1000, 1500 and 2000 m). The total number of days was 280, therefore, the total number of end points was 6x48x280 or 80640. In average, there is 80640/1200 or about 70 end points per cell.

In the present study the weight function W(*n*_{ij}) was defined as:

The weighting function reduced the PSCF values when the total number of the endpoints in a particular cell was less than about two times the average value of the end points per each cell.

### 6.2. Concentration Weighted Trajectory (CWT)

In the current PSCF method, grid cells having the same PSCF values can result from samples of slightly higher concentrations than defined by the criterion or extremely high concentrations. As a result, larger sources can not be distinguished from moderate sources. Therefore, a method of weighting trajectories with associated concentrations (CWT - concentration weighted trajectory) was developed. In this procedure, each grid cell gets a weighted concentration obtained by averaging sample concentrations that have associated trajectories that crossed that grid cell as follows:

*i,j*),

*i,j*) associated with the

For air mass trajectory visualization and statistical analysis, a software application called TrajStat was used in which clustering, PSCF and CWT methods were included [76].

### 6.3. PSCF and CWT Results

Some studies have already reported a long–range transport of PM from western countries over Balkan Peninsula which is sporadically (mostly in spring and summer) associated with African dust outbreaks [77,78]. Based on the trace metals data set contents in PM_{10} for the period 2003-2005 the possible transport process over Belgrade was investigated. The PSCF and CWT plots for PM_{10}, V, Al and Mn are presented in Figure 13. It can be seen that PSCF plot for PM_{10} shows high probability area in the northeast and west region. CWT plot indicated that the high concentrations of PM_{10} are related to the sources located in the west, northwest and northeast direction distinguishing major sources from moderate ones by calculating concentration gradients. The highest PSCF and CWT values for V were similarly distributed in the northeast region. Since V followed Weibull distribution, transportation process is expected and this region the most likely influences receptor site. Aluminum and Mn are metals related to dominant local sources. As can be seen from Figure 13, there are no significant emission region that Al can be transported from, while Mn transportation process can be expected from southeast direction.

### 6.4. Clustering of modelled backward trajectories

During 2004-2008, 48-hour backward trajectories of air masses over Belgrade were daily calculated at the previously defined six different altitudes. Transport paths of air masseswith similar history and origin were obtained by cluster analysis. Since we were concerned on the directions of the trajectories, the angle distance between back trajectories [76] has been used as the cluster model. The angle distance between two backward trajectories was defined as

where

The variables x_{0} and y_{0} define the position of the studied site and*i*-th segment for trajectories 1 and 2, respectively. In order to investigate possible seasonal variation in transportation process of PM_{10} over Belgrade, we analysed summer and winter season separately. Cluster analysis resolved the four main directions of air mass flows over Belgrade during both periods. A cluster number can be deduced through visual inspection and comparison of the mean-trajectory maps. The representative trajectories for each cluster group are presented in Figure 14. As can be seen, the four classes were related to west southwest, east, west northwest and north directions. Trajectories within different classes of airflow had distinct effects on the PM_{10} concentrations. The mean PM_{10} concentrations and frequency of air back trajectories occurrences within each cluster group are presented in Figure 15. During the summer period the most frequently arriving direction of air masses is from west northwest direction (40%) and evenly distributed among the others cluster groups. Despite the high frequency of occurrences, the contribution from this direction to PM_{10} concentration measured at the receptor site is the least. High PM_{10} concentrations are associated with the trajectories coming from west southwest and east directions. Also, the highest number of air back trajectories from these directions is associated with the observed PM_{10} concentrations above critical level of 50 µg m^{-3}. During the winter period the most frequently arriving direction of air masses is from west southwest direction (37%) and north direction (26%). Although the frequency of air back trajectories coming from east is the lowest, their contributions to the high PM_{10} concentrations are the most significant. The highest average contribution to the observed PM_{10} is from the cluster group representing the arrival direction from west southwest.

In these directions there are several possible emission sources: coal-fired thermal power plants – “Nikola Tesla A and B”, Obrenovac (west southwest); “Kostolac A and B”, Kostolac (east southeast); “Kolubara”, Veliki Crljeni (south southwest); as well the steel industry complex in Smederevo (eastsoutheast).

Obviously this kind of analysis can reveal the different transport pattern of particulate matter during summer and winter season. Cluster analysis of air back trajectories can be useful tool for investigation of the major transport pathways and their contributions to the high PM_{10} values. Additional analysis can be used to estimate background and long distance transportation of particulate matter [79,80].

## 7. Conclusion

The atmospheric processes are very complex, and pollutant concentrations are inherently random variables because of their dependence on the fluctuation of a variety of meteorological conditions and emission intensity. When sets of air pollution data are available and we aim to predict future trends, we have to know some characteristics of the measured variables, therefore, various statistical characteristics can be determined (maximum and minimum concentration, standard deviation, frequencies of exceedances of some critical level, etc) and assigned to the pollutant concentration. The distribution function is one of the most important characteristic of the random variables. Although there is a priori no reason to expect that “atmospheric” distribution should adhere to a specific probability distribution in some cases, there may be a natural probability distribution which relates emission levels to pollution levels. A great number of papers dealing to this problem have been published and in this study three theoretical distributions (lognormal, Weibull and type V Pearson) were used to fit PM_{10} concentrations in the Belgrade urban area while two-parameter exponential distribution was used to fit the high percentile PM_{10} daily concentration region. Based on the fitted distributions the minimum reduction required to attain the daily AQS was estimated using the rollback equation. These methods can reasonably predict the return period and exceedances in the succeeding period. The results suggested that the required PM source reduction ranged from 53% to 63% in Belgrade area. Recently measurements indicated that some progress in this direction has already been done and the PM_{10} mass concentrations have been reduced during the last few years.

In addition, it is demonstrated that the distribution of PM_{10} can be successfully estimated from the wind speed data if the corresponding two distributions followed the same type of distribution and have similar shape factors.

Daily and monthly variations of PM_{10} were analyzed and confirmed the existing seasonal variation. Weekend effect has also been observed and the differences between weekdays and weekend were analyzed by Q-Q plot and cumulative distribution function.

The analysis of quantiles of real and theoretical distributions for trace metals in PM_{10} indicated dominant local emission sources for the mostof metals. To identify the source locations associated with thehigh fluxes of trace metal contents in PM_{10} (Al, Mn and V), two trajectory-based models, PSCF and CWT, were used. Both PSCF and CWT resolved similar source locations for those elements. The results indicated that the source areas contributing to receptor site were located mostly in north east and west directions.

In addition, trajectory clustering was further used to investigate the transport pathways during summer and winter season. Four clusters were generated for both periods and corresponding average contributions to PM_{10} mass concentrations. Although the frequency of air back trajectories coming from east was the lowest during the summer period, their contributions to the average PM_{10} concentration were the most significant (41.7 μgm^{-3}) followed by the trajectories coming from west south west (40.4 μgm^{-3}). The highest average contribution to the observed PM_{10}concentration during the winter season was from the cluster group representing the arrival direction from west southwest (74.4 μgm^{-3}). Obviously this kind of analysis can reveal the different transport pattern of particulate matter during summer and winter season.

Nowadays commercial instruments with online measurement are helpful and can provide a large base of measurement. Still, researchers have to analyses these measurements and mathematical models and statistical analysis are some of the useful methods for developing air pollution control strategy related to specific geographic region.