Development of a versatile source apportionment analysis based on positive matrix factorization: a case study of the seasonal variation of organic aerosol sources in Estonia

. Bootstrap analysis is commonly used to capture the uncertainties of a bilinear receptor model such as the positive matrix factorization (PMF) model. This approach can estimate the factor-related uncertainties and partially assess the rotational ambiguity of the model. The selection of the environmentally plausible solutions, though, can be chal-lenging, and a systematic approach to identify and sort the factors is needed. For this, comparison of the factors between each bootstrap

Abstract. Bootstrap analysis is commonly used to capture the uncertainties of a bilinear receptor model such as the positive matrix factorization (PMF) model. This approach can estimate the factor-related uncertainties and partially assess the rotational ambiguity of the model. The selection of the environmentally plausible solutions, though, can be challenging, and a systematic approach to identify and sort the factors is needed. For this, comparison of the factors between each bootstrap run and the initial PMF output, as well as with externally determined markers, is crucial. As a result, certain solutions that exhibit suboptimal factor separation should be discarded. The retained solutions would then be used to test the robustness of the PMF output. Meanwhile, analysis of filter samples with the Aerodyne aerosol mass spectrometer and the application of PMF and bootstrap analysis on the bulk water-soluble organic aerosol mass spectra have provided insight into the source identification and their uncertainties. Here, we investigated a full yearly cycle of the sources of organic aerosol (OA) at three sites in Estonia: Tallinn (urban), Tartu (suburban) and Kohtla-Järve (KJ; industrial). We identified six OA sources and an inorganic dust factor. The primary OA types included biomass burning, dominant in winter in Tartu and accounting for 73 % ± 21 % of the total OA, primary biological OA which was abundant in Tartu and Tallinn in spring (21 % ± 8 % and 11 % ± 5 %, respectively), and two other primary OA types lower in mass. A sulfur-containing OA was related to road dust and tire abrasion which exhibited a rather stable yearly cycle, and an oil OA was connected to the oil shale industries in KJ prevailing at this site that comprises 36 % ± 14 % of the total OA in spring. The secondary OA sources were separated based on their seasonal behavior: a winter oxygenated OA dominated in winter (36 % ± 14 % for KJ, 25 % ± 9 % for Tallinn and 13 % ± 5 % for Tartu) and was correlated with benzoic and phthalic acid, implying an anthropogenic origin. A summer oxygenated OA was the main source of OA in summer at all sites (26 % ± 5 % in KJ, 41 % ± 7 % in Tallinn and 35 % ± 7 % in Tartu) and exhibited high correlations with oxidation products of a-pinene-like pinic acid and 3-methyl-1, 2, 3-butanetricarboxylic acid (MBTCA), suggesting a biogenic origin.

Introduction
Particulate matter of an aerodynamic diameter smaller than 10 µm (PM 10 ) has been extensively explored at many sites around the globe due to its various adverse effects on human health and climate. In Europe, several monitoring net-works have been measuring PM 10 for long time periods, and an increasing trend in concentrations from north to south was noticed (Fuzzi et al., 2015;Putaud et al., 2010). Despite this, some northern European countries are still suffering from PM 10 daily limit exceedances (European Environment Agency report No 13/2017No 13/ , 2017, and according to modeling studies following the "current legislation" scenarios, some of these sites will remain exposed to high PM 10 standards until 2030 (Kiesewetter et al., 2015). Therefore, understanding the origins of the pollutants can play a crucial role in future abatement policies.
While large efforts have been devoted to the investigation of the sources and the chemical composition of PM 10 , and in particular the organic fraction in western and central Europe, measurements in eastern Europe and the Baltic region are scarce. More specifically, the organic aerosol (OA) composition in Estonia has received little attention so far. PM 2.5 organic aerosol source apportionment was extensively studied by Elser et al. (2016), who performed mobile lab measurements during March 2014 at two different sites, Tallinn and Tartu. They found similar sources of OA at both sites where residential biomass burning, traffic and long-range transported OA were the major sources of OA. They also found a localized residentially influenced OA factor, which was connected to cooking activities and possibly to coal and waste burning. While the long-range transported OA dominated during nighttime and during several events when polluted air masses were transported from northern Germany, the remaining factors were important during daytime. These results provided insights into the spatial resolution of OA. Nevertheless, they were limited to short time periods; hence, the seasonal variation of the pollutants remains unknown. Residential wood combustion and traffic were also presented as important sources of PM in previous long-term air pollution studies in Estonia (Urb et al., 2005;Orru et al., 2010). However, they did not provide any quantitative source apportionment for OA.
The offline aerosol mass spectrometer (AMS) technique was recently developed by Daellenbach et al. (2016), where aqueous filter extracts are measured after nebulization with an Aerodyne high-resolution time-of-flight aerosol mass spectrometer (HR-ToF-AMS; Canagaratna et al., 2007), and the resulting organic mass spectra are analyzed with positive matrix factorization (PMF; Paatero, 1997). This technique has significantly increased our capability to investigate and identify the seasonal behavior of OA sources at several sites around the globe Bozzetti et al., 2017a). In addition, this technique allows for OA measurements of different size fractions overcoming the limitation given by the transmission window of the AMS, resulting in quantifying sources from the coarse mode, such as primary (i.e., OA directly emitted in the atmosphere), biological  or sulfur-containing primary OA sources .
PMF is widely used to analyze ambient aerosol measurement data by decomposing the input aerosol mass spectra into factor concentration time series and factor profiles. To do so, PMF iteratively solves the bilinear Eq. (1), where X i,j represents the measured input data matrix in which i is the number of samples and j is the chemical species measured, G i,k represents the concentration time series matrix in which k is the number of factors, F k,j represents the factor profile matrix, and E represents the residual matrix. The goal of PMF is to solve Eq. (1) such that the object function Q (Eq. 2) is minimized. In Eq. (2), U represents the corresponding error matrix: Bilinear models suffer from rotational ambiguity, that is, a mathematically similar goodness of fit (Henry, 1987), leading to uncertainties in extracting the contributions of different OA sources. Additional modeling errors may occur due to the user subjectivity in analyzing natural phenomena, when, for example, selecting the number of interpretable factors or estimating the error matrix. The bootstrap analysis (Davison and Hinkley, 1997), a resampling technique of the original data and error matrices, has been widely adopted to assess, to a certain extent, the rotational ambiguity related to PMF analysis (Brown et al., 2015). For each bootstrap iteration a random number of samples are selected with repeats from the original input matrices (base case) to recreate new input matrices with the same dimensions (bootstrap iteration) that will be analyzed with PMF. As a result, the bootstrap analysis results in a great number of solutions whose environmental representativeness has to be assessed, which requires a systematic approach in relating the separated factors to specific sources. This factor classification has been typically based on the contributions of certain markers (e.g., C 2 H 4 O + 2 to identify biomass burning OA or CO + 2 to identify secondary OA; Daellenbach et al., 2017Daellenbach et al., , 2018Bozzetti et al., 2017a;Vlachou et al., 2018) in the case of the application of PMF to AMS data. However, in datasets with similar factor profiles, for example two oxygenated OA factors or more, this type of sorting can become challenging, especially when a bootstrap iteration does not provide the expected factors. To assess the quality of the different bootstrap solutions, users typically compare the factor time series to external marker measurements, when available, and discard suboptimal solutions using a set of acceptance criteria (Ulbrich et al., 2009;Q. Zhang et al., 2011;Norris et al., 2014;Bozzetti et al., 2017a, b;. In this study, we propose a novel technique for evaluating the selection of the PMF solutions generated through a large number of bootstrap iterations. The method is based on the examination of the correlation matrix between the base case and bootstrap iterations for both time series and profiles, without setting a priori a defined threshold in the correlation coefficient. We assess the performance of the technique (accuracy and probability of false positive and false negative results) by comparing the factors' time series to the available specific markers. We applied the technique to an unprecedented dataset of 150 PM 10 filter samples from three sites in Estonia covering a full year (September 2013-September 2014), where anthropogenic and natural emissions of primary and secondary organic aerosols could be extracted.

Sampling sites
The samples were collected at three different sites in Estonia: Tallinn, Tartu and Kohtla-Järve (KJ). Tallinn is the capital and the largest city of Estonia located on the northern coast facing the Gulf of Finland. The measurement station is located about 9 km from the city center, in the sub-district Õismäe (59 • 24 50.927 N, 24 • 38 57.207 E; 8.5 m a.s.l.). Tartu is the second largest city of Estonia, located in the southeastern part of the country, in the valley of the Emajõgi River, a location that favors temperature inversions and the trapping of air pollutants. The measurement station (58 • 22 14.183 N, 26 • 44 5.517 E; 39.5 m a.s.l.) is positioned in the city center. According to Orru et al. (2010) traffic and local heating are important sources of air pollution at these sites. In both cities the fleet of cars increases in contrast to the limited street network capacity. Moreover, the local heating is more pronounced in Tartu compared to Tallinn. KJ is a coastal industrial city located in the northeastern part of Estonia (59 • 24 34.513 N, 27 • 16 43.166 E; 55.5 m a.s.l.). The main industries are related to large production of petroleum products, oil shale processing and electricity generation. As it is not a highly populated area, residential heating or traffic is not as important as in the other two cities.
The measurements were performed with 24 h integrated PM 10 quartz fiber filter samples from KJ (31 August 2013 to 25 August 2014), Tallinn (5 September 2013 to 1 September 2014) and Tartu (5 September 2013 to 31 September 2014; see Tables S1 and S2 and Fig. S1 in the Supplement for details). PM was collected onto 15 cm diameter quartz filters, using a high volume sampler (500 L min −1 ). After exposure, the filter samples were wrapped in lint-free paper, sealed in polyethylene bags and stored at −20 • C.

Major ionic species, sugar and acid analyses
The additional filter measurements, performed to corroborate and support the source apportionment, are listed in Table 1.

Offline AMS technique
The offline AMS technique was established by Daellenbach et al. (2016) and is briefly described in the following. From each filter sample, four punches of 16 mm in diameter were extracted in 15 mL of ultrapure water (18.2 M cm at 25 • C with total organic carbon < 3 ppb). The liquid extracts were inserted into an ultrasonic bath for 20 min at 30 • C. The ultrasonicated samples were then filtered through a nylon membrane syringe of 0.45 µm. Out of the resulting solutions, aerosols were generated in Ar (≥ 99.998 % volume, Carbagas, Gümligen, Switzerland) via an apex Q nebulizer (Elemental Scientific, Inc., Omaha, NE, USA) operating at 60 • C and subsequently directed into the AMS after being dried by a Nafion dryer.
The technique was performed on 150 filter samples in total: 39 from KJ, 69 from Tallinn and 42 from Tartu (Tables S1  and S2). The resulting organic spectra were analyzed by PMF with the use of the multilinear engine 2 (ME-2; Paatero, 1999). The interface for the data processing was provided by the source finder toolkit (SoFi version 4.9; Canonaco et al., 2013) for Igor Pro (WaveMetrics, Inc., Portland, Oregon, USA).
The NH 4 NO 3 artifact on the CO + 2 signal (Pieber et al., 2016) was also accounted for. For a thorough description of the artifact and the correction procedure, the reader is referred to Pieber et al. (2016) and to Daellenbach et al. (2017).

PMF input and uncertainties
As stated in the introduction, the input data matrix for the PMF is the organic mass spectrum data matrix X i,j which consists of a combination of factor profiles, and time series and the input error matrix U include the blank variability and the measurement uncertainties. Before using the PMF algorithm, all the fragments with a signal-to-noise ratio (SNR) below 0.2 were removed from the input matrix, and the ones with SNR below 2 were down weighted according to the recommendations of Paatero and Hopke (2009). Note that the PMF input matrix X i,j included the data from all three sites.
To obtain quantitative results, both data and error matrices were multiplied by the externally measured water-soluble OC (WSOC) times the OM : OC ratio retrieved from the factor profiles of the matrix X i,j .
Even though traffic is expected to be one of the primary sources of air pollutants, especially in Tallinn and Tartu, a clear hydrocarbon-like OA (HOA) which mainly contains non-water-soluble compounds could not be identified due the extraction procedure used. To assess a potential effect of the water-soluble HOA (WSHOA) on the PMF results, we estimate the WSHOA contribution to the different fragments in the data matrix. The calculation was based on the time series of the concentration of elemental carbon (EC) and the averaged high-resolution HOA reference factor profiles from Crippa et al. (2013) and Mohr et al. (2012)  Liquid chromatographyelectrospray ionization mass spectrometry (Jacob et al., 2019) Organic acids (e.g., benzoic acid and pinic acid)

(Tallinn)
Sunset EC-OC analyzer (Birch and Carry, 1996) with the EU-SAAR 2 protocol (Cavalli et al., 2010) Organic carbon (OC) and elemental carbon ( 4 for more information on the recoveries). The calculated WSHOA data matrix was then subtracted from the original data matrix. The PMF output did not change with the subtraction of WSHOA, even though the calculated concentration of the latter was a high estimate due to the assumption that EC only originates from traffic. A thorough apportionment of EC and the calculation of HOC will be discussed in Sect. 4.1. The variability in the AMS input dataset was best described by a seven-factor solution, which will be thoroughly described in Sect. 3.1 below. To assess the stability of the PMF solution and the sensitivity of the model for the WSHOA subtraction, we performed 250 bootstrap runs by perturbing the HOA : EC and R HOA parameters within their errors (1 standard deviation, σ ), assuming a normal distribution. Note that this number of runs was restricted by computational limitations.
The sorting of the factors and the concomitant selection of the retainable solutions out of the 250 runs was based on the correlation (R Spearman -R s ) of the time series between the base case (which is the PMF result before bootstrapping) and each bootstrap iteration n as well as the correlation (R s ) of each factor profile between the base case and each iteration n. In the following, the sorting based on the time series is denoted as "ts", and the sorting based on profiles is denoted as "pr". The criteria, followed for the selection or rejection of each solution, are described in Sect. 3.2.

Interpretation of PMF factors
As already mentioned, the variability in the water-soluble organic mass spectra was best explained by a seven-factor solution, which we refer to as the base case. The factors found were as follows: 1. An oil-related OA was found that was rich in hydrocarbons (Fig. 1a) and showed elevated concentrations, mainly in KJ (Fig. 1b).
2. A sulfur-containing OA (SCOA) factor was found with a pronounced peak at m/z 79 (CH 3 SO + 2 ; Fig. 1a) and rather stable contributions at all sites (Fig. 1b). As this factor was significantly enriched in the coarse mode , mainly found in urban areas influenced by traffic at other European sites  and clearly associated with fossil carbon , we have previously related it to asphalt abrasion or tire wear.
3. A summer oxygenated OA (SOOA) was found with enhanced m/z 43 (C 3 H 2 O + ) and 44 (CO + 2 ) peaks (Fig. 1a), which was highest during summer at all sites (Fig. 1b). 4. A winter oxygenated OA (WOOA) was found with enhanced peaks at m/z 28 (CO + ) and 44 (CO + 2 ; Fig. 1a), dominating in winter at all sites (Fig. 1b). Both of these oxygenated factors (SOOA and WOOA) were also found in different European sites and were connected to non-fossil sources, biogenic and anthropogenic, respectively Daellenbach et al., 2018). Such a distinction was also found in Canonaco et al. (2015), where an online ACSM was used.

5.
A factor with a significantly pronounced peak at m/z 44 (CO + 2 ; Fig. 1a) and elevated concentrations in summer ( Fig. 1b) was found, which was identified as dust. This factor will be more thoroughly examined in Sect. 4. 6. A primary biological OA (PBOA) was found which exhibited high contributions of the fragment C 2 H 5 O + 2 (at m/z 61; Bozzetti et al., 2016;Fig. 1a) and increased concentrations during late spring and summer at all sites (Fig. 1b). 7. A biomass burning OA (BBOA) was found with a characteristic peak at m/z 60 (C 2 H 4 O + 2 ; Alfarra et al., 2007;Fig. 1a) and elevated concentrations during late fall and winter in Tallinn and especially in Tartu (Fig. 1b).

PMF uncertainty analysis: factor sorting and solution selection
The framework for factor sorting and solution selection proceeded as follows: 1. A correlation matrix was composed, including all the correlations between base-case factor time series (profiles) represented in rows and bootstrap iteration n time series (profiles) represented in columns, demonstrating the R s per correlation (Fig. 2).
2. Factors were sorted according to the highest correlation of their time series (profiles) with the base-case factor time series (profile).
3. Solutions were discarded if any of the correlation coefficients occurring in the matrix diagonal were not statistically significantly higher than at least one of the coefficients in the respective column or row (significance level α = 0.05). These selection criteria have two implications: (1) every factor separated in a bootstrap run should correspond to a unique factor of the seven factors separated in the base case, and (2) all factors that could be identified in the base case have one unambiguous corresponding factor in the bootstrap run. We have statistically evaluated the comparison between the Spearman coefficients by treating them as if they were Pearson coefficients (Myers et al., 2006) and by applying the standard Fisher's z transformation and subsequent comparison using a t test. This approach was reported to be more robust with respect to the type I error (false positive) than ignoring the non-normality and using Pearson instead of Spearman coefficients.
An example of the correlation matrix of a retained solution is shown in Fig. 2a for time series and Fig. 2b for profiles. Meanwhile, Fig. 2c and Fig. 2d represent examples of a bootstrap iteration (n = 140) where the solution was rejected because SCOA was not resolved, based on both ts and pr analysis, respectively. In Fig. 2c, the highest correlation between the time series of factor 2 and the factors of the base case was found with dust instead of SCOA, yet it was much weaker (R s = 0.34) than the correlation between factor 7 and dust (R s = 0.79). Therefore, factor 7 could be identified as dust, and factor 2 could not be identified as an interpretable factor. In Fig. 2d, factor profile 2 correlated most with the base-case SCOA profile; however this correlation was not significantly higher than the correlation between factor 2 and dust. Moreover, base-case SCOA correlated better with factor 6, which was related to SOOA, indicating that SCOA could not be unambiguously related to factor 2.
To validate the selection of the solutions, we compared the factors of each bootstrap iteration with an externally measured marker, more specifically BBOA with levoglucosan, PBOA with cellulose, WOOA with phthalic acid and SOOA with 3-methyl-1,2,3-butanetricarboxylic acid (MBTCA). The retained solutions exhibited the highest correlations between the external markers and the respective factors (red markers in Fig. 3). To seal the validity of the retained solutions, we also compared the R s between the base-case factors and their respective external marker with the R s between the bootstrap iteration factors and their respective external markers. We performed the bootstrap technique for a second time on the time series of the base-case factors and the respective external markers 1000 times. The resulting R s coefficients are represented in probability density functions (PDFs) indicated as red curves in Fig. 3, centered at 0.8 for BBOA (Fig. 3a), 0.7 for WOOA (Fig. 3c) and 0.9 for SOOA (Fig. 3d) and a broader one centered at 0.45 for PBOA (Fig. 3b). In all four cases, the retained solutions, either coming from the ts or the pr approach, spanned around the center of each PDF (Fig. 3 and Fig. S2 for the pr), and most of the solutions where at least one factor was not resolved (black markers in Fig. 3) were not included within the PDF  boundaries. The general agreement between the PDFs and the retained solutions ratified the solution selection approach; however, there were still some cases of possible misallocation of retained or rejected solutions (for example, a few black markers appearing at the center of the PDF; Fig. 3c).
To assess whether the retained bootstrap solutions share the same quality with the base-case solution with regard to correlations between factors and markers, we calculated the probability of type I (false positive) and type II (false negative) errors associated with the solution selection approach (Fig. 4). The analysis entailed a quantitative comparison of the Spearman coefficients obtained between markers and factor time series from the bootstrap iterations, R sboots , with the respective Spearman coefficients obtained between markers and factor time series from the base case R sbase , considering the same samples as in the corresponding bootstrap iterations. The comparison between R sboots and R sbase was performed by applying a Fisher transformation followed by a t test. We defined true positive and false negative as the red (retained solutions) and black (non-resolved factors) points, respectively, lying within the PDF boundaries with regard to the total number of red and black points within the PDF boundaries. True negative and false positive were defined as the black and red points lying outside the PDF boundaries with respect to the total number of red and black points outside the PDF boundaries. The limits between false positive and false negative were set by 2 standard deviations from the one-to-one line. The percentages of the accuracy and the probability of false positive or false negative cases are compiled in Table S3. Sorting based on profiles seemed less reliable and is more prone to false negative solutions (Tables S3  and S4), as the profiles often look similar, and therefore the R s exhibits high values for all factors (Fig. 2b and d). On the contrary, sorting based on time series showed clearer results, as the R s spanned over a greater range of values ( Fig. 2a  and c). Still the ts method produced false negative solutions, for example 53 % for PBOA due to the combination of (i) a limited number of points available for cellulose and (ii) the representativeness of the marker time series after the resampling (bootstrap analysis). Note that PBOA was important only during a few days in spring, and therefore it is possible that these days were not always selected in the resampling process. The SOOA on the other hand exhibited 0 % false negative and 16 % false positive cases, always demonstrating high R sboots and R sbase values.
However, in either of the two methods ts and pr, the Fishertransformed correlation coefficient rendered the selection of the solution evident, and eventually the two sorting methods yielded a very similar retained solution space (Figs. 4 and 5). Figure 5 depicts the correlation between the averaged common retained solutions and the averaged retained solutions coming from either the ts or the pr sorting method for the example of BBOA (correlations for the other factors are shown in Fig. S3). There is a minor deviation from the oneto-one line for the standard deviation scatterplot (Fig. 5b) for the ts sorting method. However, as soon as the solutions were weighted according to the correlation between external marker and bootstrap run time series, then the deviation decreased (blue markers in Fig. 5). The weighting factor w i was calculated as where SE is the standard error resulting from the regression between the external marker and bootstrap iteration.

Estimation of traffic contribution to EC and OC
We estimated above that the traffic contribution to WSOA (< 1 %) can be neglected and has little effect on the PMF outputs. However, traffic might be an important source of EC and OA, which is assessed in the following.
To estimate the percentage of EC coming from traffic (EC tr ), we used the ME-2 model (with SoFi standard version 6.399; Canonaco et al., 2013), assuming that the sources of EC are biomass burning, resuspension of road dust, industrial emissions from the oil shale factories and traffic. Here the input data matrix included the time series of water-soluble biomass burning OA (WSBBOA), watersoluble sulfur-containing OA (WSSCOA), water-soluble oil OA (WSOilOA) and EC. PMF was conducted 1000 times, varying the initial seed to solve Eq. (4): EC =EC tr + EC bb + EC oil + EC rrd = EC tr + a · WSOA bb Here, EC bb , EC oil and EC rrd represent the EC concentration time series related to the primary sources biomass burning, oil shale industry, and resuspension of road dust and tire wear, respectively, while a, b and c are the EC : WSOA ratios characteristic of the emissions from the same sources and were obtained as outputs of the model. This new methodology, based on PMF, is especially pertinent, as unlike other inversion techniques, it sets positive constraints on a, b and c and offers the possibility of resolving the contributions of factors for which no constraints are available, here EC tr . We found that EC tr contributed 57 % ± 5 % to the total EC (on a yearly average), while 36 % ± 5 % of EC was attributed to biomass burning, 4 % ± 1 % to road dust resuspension and 3 % ± 1 % to the oil shale emissions (Fig. S4). The contribution of EC from fossil fuel combustion (EC ff ) measured at a site similar to Tartu, i.e., an Alpine valley in Magadino, southern Switzerland, in 2014, was in agreement with our EC tr contribution, with a yearly average of 55 % ± 7 %. Also, in Zurich, an urban site, EC ff ranged from 40 % to 55 % during the winter of 2012 (Zotter  , 2014). From the EC tr contribution, we estimated that the HOC (obtained by multiplication of the EC tr time series with the HOC : EC ratio) contributed 4 % to the total OC on a yearly average.

Scaling to organic carbon
All the retained solutions (in total 62 %) were averaged per factor, and their seasonal behavior as well as their correlations, with the external markers, are presented in Sect. 4.4. Note that all the water-soluble factors were recovered following the method described in Daellenbach et al. (2016) and Vlachou et al. (2018). The recoveries were calculated by fitting Eq. (5): where OC i,n is the OC concentration per bootstrap run n per sample i, R k is the recovery R per factor k, and WSOC i,n,k is the water-soluble OC concentration calculated based on the measured WSOC and the OM : OC per factor. From the OC i,n the part of hydrocarbon-like OC and the inorganic carbon related to dust were removed. The carbonate-carbon investigation and calculation is discussed in Sect. 4.3. To define the recoveries, we used a non-negative multilinear fit. The starting points of the fitting for each R k with the exception of R oil were obtained from the literature Daellenbach et al. 2016;Vlachou et al., 2018) and were randomly varied within their literature range with an increment of 10 −4 . The final distributions of the recoveries are shown in the Supplement (Fig. S5). The recoveries in this study were all shifted to the lower end of the recoveries reported in the literature. While the reason remains unclear, the water solubility of OA is specific to the dataset; therefore we can expect differences to other datasets. Moreover, we re-measured, in a different laboratory, the OC concentrations from a subset of 21 samples covering all sites and all seasons. The agreement between the two differently determined OC concentrations was excellent ( Fig. S6; slope = 0.93, R 2 = 0.99).

Exploration of the dust factor
Mineral dust can contain a significant amount of inorganic carbon in the form of CO 2− 3 . The water extracts used in the offline AMS technique have a pH that is always < 8. Therefore, the CO 2− 3 in our samples is all solubilized into HCO − 3 (pKa (HCO − 3 /CO 2− 3 ) = 10.33; Bruice, 2010). This is shown to release CO 2 when it undergoes thermal decomposition on  the AMS vaporizer (at 600 • C; Bozzetti et al., 2017b). Thus, the contribution of CO + 2 to organic aerosols is overestimated and the fraction coming from the inorganic carbon should be removed from the OA spectra.
To remove the influence of the inorganic dust, we estimated the carbonate-carbon concentrations [C_CO 3 ] corrected for the relative ionization efficiency, as discussed in the Supplement. This estimated C_CO 3 concentration was compared to measured carbonate on a subset of 19 filter samples. While the agreement between measured and estimated concentrations is poor for Tartu, a decent agreement was found for KJ and Tallinn, especially given the large uncertainties in both variables (Fig. S7). On an absolute basis, PMF seems to overestimate the C_CO 3 concentrations by ∼ 20 % compared to the measured concentrations.
Ca 2+ is one of the most common constituents of mineral dust and can therefore be used to trace this source. The time series of C_CO 3 and Ca 2+ , available for the entire set of samples, displayed in Fig. 6, show that the two variables exhibit similar trends (except for Tallinn). Despite the large errors in the [C_CO 3 ] estimates, an uncertainty weighted correlation analysis (Supplement and Fig. S8) shows that [C_CO 3 ] and [Ca 2+ ] correlate statistically significantly (R = 0.4, p < 10 −5 ) with a slope of 0.2, consistent with C_CO 3 and Ca 2+ being in the form of calcium bicarbonate.
We have validated the identity of the dust factor even further by measuring the same subset of 19 filters with the offline AMS technique before and after fumigation with HCl (as described in Zhang et al., 2016). The comparison of the mass spectra of fumigated and non-fumigated samples is illustrated in Fig. 7 for two samples: with high and low contribution of dust. In the example of KJ (5 June 2014), where the dust factor exhibited the highest contribution, f 44 was substantially decreased after fumigation (Fig. 7a). In the case of Tallinn (19 January 2014), where the dust factor concentration was negligible, f 44 remained stable after fumigation (Fig. 7b). Overall, the comparison of the f 44 modeled (= f 44 total -f 44 org ) from the initial dataset of 150 filters and the f 44 measured (= f 44 non_fum -f 44 fum ) from the subset of 19 filters showed consistent results (Fig. 8). Taken together these results provide strong confidence on the nature of the dust factor extracted by PMF.

Seasonal variation of organic aerosol sources
The sources were quantified after removing the contribution of the dust factor from the total OA. All the factor concentrations with their uncertainties averaged per season are presented in Table S3. In general, the relative uncertainties decreased with increasing concentrations per factor (Fig. S9). For concentrations above ∼ 1 µg m −3 the percentage error became more important than the error related to noise and was thus more stable for all factors. Consistent with the factor separation and uncertainty analysis above, the factors that were well separated, such as SOOA, exhibited low relative uncertainty (0.15), while the factors that were more difficult to extract, such as BBOA, exhibited higher relative uncertainty (0.45).
BBOA exhibited high concentrations in Tallinn during winter (on average 3.7 ± 2.7 µg m −3 ) and fall (1.2 ± 0.9 µg m −3 ) and in Tartu (8.4 ± 3.9 and  3.8 ± 1.9 µg m −3 , winter and fall, respectively; Fig. 9a; Table S5). In both cities BBOA and the marker levoglucosan correlated very well (Fig. 9b) confirming the identity of the factor. For KJ the concentrations of BBOA were lower in winter (1.3±0.8 µg m −3 ) as expected, due to the low number of residents and low biomass burning activities in the region. At all sites the levoglucosan-to-BBOA ratio (0.08 for KJ, 0.05 for Tallinn and 0.05 for Tartu) assessed in this study was consistent with the one reported in the neighboring country, Lithuania (Bozzetti et al., 2017a). Potassium (K + ), which is often used as a wood-burning marker, especially for ash, also correlated well with BBOA for Tallinn (R 2 = 0.80) and Tartu (R 2 = 0.58). Different BBOA : K + ratios were used in the past to describe burning conditions Daellenbach et al., 2018), and higher values were linked to inefficient burning conditions. Here, the BBOA : K + ratio (14.3 for Tallinn and 18.1 for Tartu) was in agreement with the one found at southern Alpine valley sites (Magadino and San Vittore, Switzerland; Daellenbach et al., 2018) where older stoves are still used. In Estonia more than 80 % of households use old types of stoves for heating purposes (Maasikmets et al., 2015); therefore, BBOA could be linked to inefficient residential wood burning.
The recognition of PBOA as described in Sect. 3.2 was also supported by the high correlations of this factor with cellulose and erythritol (Fig. 9d). Cellulose is related to plant debris and is typically used as a marker for primary biological aerosols , while erythritol, among other sugar alcohols, reflects airborne biogenic detritus (Graham et al., 2002). The seasonal behavior of PBOA was very similar to the respective behavior of both markers (Fig. 9c), with average spring concentrations of 0.2 ± 0.2 µg m −3 for KJ, 1.2 ± 0.8 µg m −3 for Tallinn and 0.7 ± 0.4 µg m −3 for Tartu.
SOOA exhibited a clear yearly cycle at all sites, with the highest concentrations witnessed in summer and early fall (in summer on average 1.8 ± 0.7 µg m −3 for KJ, 2.8 ± 0.6 µg m −3 for Tallinn and 2.1 ± 0.4 µg m −3 for Tartu; Fig. 9e). In previous studies this factor showed an exponential increase (Fig. S10) with temperature and was linked to terpene oxidation products Leaitch et al., 2011). In another study in an Alpine valley site, this factor was also found to be 79 % non-fossil, supporting the connection to biogenic secondary OA . Here, SOOA not only correlated with temperature (R s = 0.77; Fig. S10) but also with two oxidation products of a-pinene, i.e., with pinic acid and even better with MBTCA (Fig. 9f). The latter was shown to be produced by reaction of pinonic acid with the OH radical (Müller et al., 2012).
WOOA was more pronounced during fall and winter at all sites, with average concentrations in winter of 1.4 ± 0.5 µg m −3 for KJ, 2.2 ± 0.8 µg m −3 for Tallinn, and 1.5 ± 0.5 µg m −3 for Tartu (Fig. 9g). In earlier studies WOOA was characterized as non-fossil  and was identified based on its correlation with NH + 4 coming mainly from NH 4 NO 3 in the wintertime Daellenbach et al., 2017). It was also related to long-range transported oxygenated OA stemming from anthropogenic emissions during winter, such as biomass burning. Also, WOOA demonstrated high correlations with two anthropogenic organic acids, benzoic and phthalic acid (Fig. 9h), formed via the photo-oxidation of aromatic hydrocarbons such as toluene and naphthalene and therefore suggested to be tracers for anthropogenic sources (Kawamura and Yasui, 2005;Deshmukh et al., 2016). Recently, Bruns  (2016) found that aromatic compounds, such as benzene and naphthalene emitted by wood combustion, can indeed produce highly oxygenated SOA. Taking all the above into account, it was concluded that WOOA might be linked mostly to aged wood-burning OA. Figure 10 illustrates the relative contributions per factor per site to the total OA (all averaged contributions per season per factor with their uncertainties shown in Table S6). Out of the primary sources, the major contributor was BBOA during winter and fall (on average and 1 standard deviation: 39 % ± 16 % and 27 % ± 13 % in Tallinn and 73 % ± 21 % and 53 % ± 14 % in Tartu). However, in KJ during winter and fall, WOOA was the dominant source (36 % ± 14 % and 39 % ± 13 %, respectively), indicating that for this site, regional transport of OA is important. In spring, PBOA was the major source in Tartu (21 % ± 8 %), while in Tallinn BBOA and SOOA were the dominating sources during that season (30 % ± 14 % and 18 % ± 5 %, respectively). This could be due to the fact that temperatures in early spring are still low (2 • C on average in Tallinn in March) and the wood burning for residential heating is still widely used. Towards the end of spring (15 • C on average in Tallinn in May) the rising temperature favors the biogenic emissions. In KJ, the most important source was oil OA (36 % ± 14 % in spring), most possibly coming from the oil shale industries in the region. The presence of the oil factor at the other two sites could be an indication that this factor is mixed with coal or waste burning, as also found by Elser et al. (2016). Besides this, the oil OA profile resembled the coal profile identified in the city of Cork, Ireland (Dall'Osto et al., 2013). During the summer months and early fall, SOOA prevailed over all sources at all sites, with 26 % ± 5 % in KJ, 41 % ± 7 % in Tallinn and 35 % ± 7 % in Tartu. Even though KJ is highly industrialized, SOOA can still be related to the production of secondary OA from biogenic volatile organic compounds. The least-significant source, especially in Tartu, with rather stable seasonal behavior, was SCOA. The yearly average contribution of SCOA was 12 % ± 4 % in KJ, 14 % ± 5 % in Tallinn and 4 % ± 2 % in Tartu. Although it is generally found that the secondary sources prevail over the primary ones, here the primary sources seem to dominate the secondary ones. This is also observed at other European sites such as Payerne , for coarse particles) or Magadino, especially in winter (Vlachou et al., 2018, for PM 10 ), as well as in Beijing, China , for PM 1 ).

Conclusions
The offline AMS technique was applied to a set of 150 filter samples covering a yearly cycle at three sites in Estonia. The uncertainties of the PMF solution were assessed by bootstrap analysis. In order to identify the factors, the Spearman R (R s ) coefficients between base-case time series and bootstrap run time series (ts) as well as base-case profiles and bootstrap run profiles (pr) were monitored. The results showed that the retained solution space if one follows the ts or the pr sorting method is very similar. Weighting with the R s between external markers and bootstrap run increased our confidence towards the solution space. The source apportionment results revealed four primary OA sources, two secondary OA and a dust factor. The dust factor was identified by measurements of calcium carbonate as well as by acidification with HCl of a selected batch of filters. Out of the primary sources, three had an anthropogenic influence. BBOA was mainly present in winter and autumn in Tallinn and Tartu, the two largest cities of Estonia, where residential heating activities are common. SCOA was mostly important in winter in Tallinn and KJ, in contrast to Tartu. The third anthropogenic primary factor was oil OA, which exhibited the highest concentrations in KJ, as expected. The reason why this factor was evident in Tallinn and Tartu could be that it may include coal combustion for residential heating purposes. PBOA was the only primary OA not related to anthropogenic emissions and prevailed in spring at all sites. The two oxygenated OA factors were separated according to their seasonal behavior: WOOA was linked to anthropogenic wood-burning activities, as it dominated in winter and autumn at all sites and also correlated with phthalic and benzoic acid. SOOA was significant during summer at all sites and was related to biogenic emissions and strong aging, as it was highly correlated with a second-generation oxidation product of a-pinene.
Author contributions. AV performed the offline AMS measurements, data curation and formal analysis; AV and IEH wrote the paper; AT assisted with the formal analysis; HL provided expertise on AMS measurements; FC and KD provided expertise on software; JLJ performed additional measurements including major ions, sugars, organic acids, cellulose, OC-EC and WSOC; MCM performed the carbonate measurements; MM and ET organized the filter sampling; and UB, ASHP and IEH were involved with the supervision and conceptualization. All authors commented on the paper and assisted in the interpretation of the results.
Competing interests. The authors declare that they have no conflict of interest.