## 2 Materials and Methods
### PMF Modelling  
One of the methods that allows for the identification of particulate matter sources and the quantification of their contribution is the positive matrix factorization method (PMF) introduced by [PERSON] et al. [36] Due to the need to have many samples (at least several dozen), it is often classified as a statistical method. The basic assumption of the method is a constant relative share of components that characterize a given source. This share is called the profile of a given source. The PMF method, based on the matrix of the sample of chemical species of particulate matter, calculates the matrix of participation of a given number of factors and profiles of these factors [37]. The PMF receptor model solves the set of equations:  
\[x_{ij}=\sum_{k=1}^{P}g_{ik}f_{kj}+e_{ij}, \tag{2}\]  
where \(x_{ij}\) is an element of the concentration matrix X (i sample index and j species index), gik is an element of the contribution matrix G with p sources (k is the index of the sources), fij is an element of F the source and, finally, eij is an element of the residual matrix E [38]. The PMF multivariate statistical method decomposes the concentration matrix (X) into the source contribution (G) and source profile (F) matrixes in such a way that G and F obtain non-negative values only, ensuring the physical meaning of the model. The profile determines the share of individual components in each factor and is the basis for the physical assignment of a given factor to identify the sources of particulate matter. In this method, the number of factors is set arbitrarily. In practice, modeling is carried out for a different number of factors, and the number for which the determination of sources is unambiguous is finally selected [26]. In this work, the EPA PMF 5.0 software, developed by the US Environmental Protection Agency (US EPA), was used. As input, not only is the concentration matrix (X) required but also the appropriate uncertainty matrix (U). For a given number of factors (p), the matrixes G and F are adjusted by minimizing the objective function Q, which is defined as follows:  
\[Q=\sum_{j=1}^{m}\sum_{i=1}^{n}\frac{e_{ij}^{2}}{u_{ij}^{2}}, \tag{3}\]  
where uij is an element of the uncertainty matrix (U), m is the number of species and n is the number of samples. A general optimization method where the measured value is 'weighted' by its uncertainty requires an accurate estimate of the uncertainty. The following 20 chemical species identified in elemental and ions analyses of the PM10 samples were used: S, K, Ti, Cr, Mn, Fe, Co, Cu, Zn, Br, Rb, Sr, Pb, Cl\({}^{-}\), NO\({}_{3}\)\({}^{-}\), \(\text{WO}_{4}\)\({}^{3-}\), \(\text{NH}_{4}\)\({}^{+}\), \(\text{Na}\)\({}^{+}\), \(\text{Mg}_{2}\) and \(\text{Ca}^{2+}\). All of these data were classified as'strong', which resulted in a high signal-to-noise ratio. In the present work, if the concentration is less than or equal to the detection limit (LOD) for a given element, the uncertainty is set at 5/6 LOD and the concentration replaced by 1/2 LOD [26; 39]. The missing data were substituted by median values, and the corresponding uncertainties were set at four times the LOD median value [26; 39]. After the factorization run, the PMF software provides the possibility of analyzing the factorization stability by the \"Fpeak Bootstrap Method\" [40]. Modeling was carried out together for the data from 2018 and 2020 years. The following values were obtained for Q robust 18,274; Q true 64,574.