Source Apportionment and Toxic Potency of Polycyclic Aromatic Hydrocarbons (PAHs) in the Air of Harbin, a Cold City in Northern China

: A total of 68 PUF samples were collected seasonally from 17 sampling sites in Harbin, China from May 2016 to April 2017 for analyzing 15 congeners of gaseous polycyclic aromatic hydrocarbons ( Σ 15 PAHs). An improved non-negative matrix (NMF) model and a positive matrix factorization (PMF) model were used to apportion the sources of PAHs. The carcinogenic risk due to exposure to PAHs was estimated by the toxicity equivalent of BaP (BaPeq). The results showed that the average concentration of Σ 15 PAHs was 68.3 ± 22.3 ng/m 3 , and the proportions of 3-ring, 4-ring, 5-ring, and 6-ring PAHs were 64.4%, 32.6%, 2.10%, and 0.89%, respectively. Among the six typical functional areas in Harbin, the Σ 15 PAHs concentrations were 98.1 ± 76.7 ng/m 3 , 91.2 ± 76.2 ng/m 3 , 71.4 ± 75.6 ng/m 3 , 67.9 ± 65.6 ng/m 3 , 42.6 ± 34.7 ng/m 3 , and 38.5 ± 38.0 ng/m 3 in the wastewater treatment plant, industrial zone, business district, residential area, school, and suburb, respectively. During the sampling period, the highest concentration of Σ 15 PAHs was in winter. The improved NMF model and PMF model apportioned the PAHs into three sources including coal combustion, biomass burning, and vehicle exhaust. The contributions of coal combustion, biomass burning, and vehicle exhausts were 34.6 ± 3.22%, 48.6 ± 4.03%, and 16.8 ± 5.06%, respectively. Biomass burning was the largest contributor of Σ 15 PAHs concentrations in winter and coal combustion contributed signiﬁcantly to the concentrations in summer. The average Σ BaPeq concentration was 0.54 ± 0.23 ng/m 3 during the sampling period, high concentrations occurred in the cold season and low levels presented in the warm period. Vehicle exhaust was the largest contributor to the Σ BaPeq concentration of PAHs in Harbin.


Introduction
Polycyclic aromatic hydrocarbons (PAHs) have received considerable attention in recent decades because of their relatively high burden on the environment and strong toxic potency [1,2]. PAHs are a diverse class of aromatic compounds containing at least two benzene rings, which have carcinogenic, teratogenic, and mutagenic toxicity [1,3]. These compounds are mainly derived from the incomplete combustion of carbonaceous materials, associated with anthropogenic activities in particular [4,5]. The anthropogenic emission sources primarily include the combustion of fossil fuels and biomass for domestic and

Sample Collection
The air samples were collected at 17 sampling sites of 6 types of areas in Harbin by PUF disks. The locations of the sampling sites were showed in Figure 1. The 6 types of areas include industrial zone, business district, residential area, school, sewage treatment plant, and suburb. There were no obvious pollution sources around each sampling site, and the installation height of the PUF sampler was more than 1.5 m. The structure and size of the PUF sampler used in this study were the same as those widely used for passive sampling of atmospheric persistent organic pollutants (POPs) [17]. The diameter, thickness, and density of the PUF disks were 14 cm, 1.35 cm, and 0.0213 g/cm 3 , respectively. Passive air samples were collected seasonally from May 2016 to April 2017 according to summer, autumn, winter, and spring, and a total of 68 PUF disk samples were collected. The PUF disks were Soxhlet-extracted with acetone for 24 h and dichloromethane (DCM) for 24 h, respectively. The extracted PUF disks were wrapped with aluminum foil after being incandesced at 450 °C for 4 h, placed in disk-shaped aluminum cans, sealed in polyethylene bags, and placed in a −20 °C freezer before the sampling campaign. The PUF disks were used according to the standard procedures [18]. After a 90-day sampling period, the PUF disks were wrapped, sealed, and stored in a −20 °C freezer in the lab until further analyses. Six extracted PUF disks were exposed to the air for a few seconds during the sampling period and then stored at −20 °C in a freezer for use as field blank samples.

Analytical Procedure
Determinations of PAHs followed the methods introduced in previous studies [4,16,19]. Briefly, the exposed PUF disks were spiked with a mixture of surrogate standards (naphthalene-d8, acenaphthene-d10, phenanthrene-d10, chrysene-d12, and perylene-d12) and Soxhlet-extracted with dichloromethane for 24 h. Activated copper granules were added into the collection flask to remove sulfur. The extracts were concentrated with a rotary evaporator and solvent-exchanged into n-hexane with a volume of ~0.5 mL. Purification was carried out by a multilayer column containing anhydrous sodium sulfate, neutral silica gel, and neutral alumina from top to bottom. PAHs were eluted with 20 mL of a mixture solvent of dichloromethane and n-hexane (1:1, v/v). The eluent solvent was blown down to approximately 25 μL, and hexamethylbenzene was added to each concentrated sample as an internal standard before the instrumental analysis.
The PAHs were analyzed by a gas chromatography (Agilent 7890a, Agilent Technologies Inc., CA, USA) equipped with a capillary column (DB-5MS, 30 m × 0.25 mm × 0.25 μm) and a mass spectrometer (MSD, Agilent 5975c, Agilent Technologies Inc., CA, USA). The parameter setting was also obtained from previous studies [4,16,19]. Briefly, each extracted sample (1 μL) was injected in a splitless mode with a solvent delay time of 10 min. High-purity helium was applied as the carrier gas with a flow rate of 1.83 mL/min. The temperatures of the injector and transfer lines were 280 °C and 300 °C, respectively. The initial oven temperature was set at 60 °C for 5 min, elevated to 280 °C at a heating rate of 3 °C/min, The PUF disks were Soxhlet-extracted with acetone for 24 h and dichloromethane (DCM) for 24 h, respectively. The extracted PUF disks were wrapped with aluminum foil after being incandesced at 450 • C for 4 h, placed in disk-shaped aluminum cans, sealed in polyethylene bags, and placed in a −20 • C freezer before the sampling campaign. The PUF disks were used according to the standard procedures [18]. After a 90-day sampling period, the PUF disks were wrapped, sealed, and stored in a −20 • C freezer in the lab until further analyses. Six extracted PUF disks were exposed to the air for a few seconds during the sampling period and then stored at −20 • C in a freezer for use as field blank samples.

Analytical Procedure
Determinations of PAHs followed the methods introduced in previous studies [4,16,19]. Briefly, the exposed PUF disks were spiked with a mixture of surrogate standards (naphthalene-d 8 , acenaphthene-d 10 , phenanthrene-d 10 , chrysene-d 12 , and perylene-d 12 ) and Soxhlet-extracted with dichloromethane for 24 h. Activated copper granules were added into the collection flask to remove sulfur. The extracts were concentrated with a rotary evaporator and solvent-exchanged into n-hexane with a volume of~0.5 mL. Purification was carried out by a multilayer column containing anhydrous sodium sulfate, neutral silica gel, and neutral alumina from top to bottom. PAHs were eluted with 20 mL of a mixture solvent of dichloromethane and n-hexane (1:1, v/v). The eluent solvent was blown down to approximately 25 µL, and hexamethylbenzene was added to each concentrated sample as an internal standard before the instrumental analysis.

Quality Assurance and Quality Control
A procedural blank, a field blank, and a real sample selected randomly were run with each batch of 10 real samples to assess potential contamination during the analysis. The results showed that the change of PAHs concentrations in duplicates was less than 15%. The target compounds in the procedural and field blanks were not detected except for Nap. The deuterated PAHs were added to samples and blanks to monitor the analytical and sampling procedure. The deuterated PAHs recoveries were 38 ± 36%, 83 ± 15%, 108 ± 26%, 91 ± 12%, and 92 ± 17% for naphthalene-d 8 , acenaphthene-d 10 , phenanthrened 10 , chrysene-d 12 , and perylene-d 12 , respectively. The method detection limits (MDLs) were defined as the average concentrations of target PAH compounds in field blanks plus three times the standard deviation of the blanks. The MDLs of target PAH compounds were in a range from 0.01 to 3.7 pg/m 3 . All the PAHs concentrations reported in this study were not corrected with surrogate recoveries. Nap was not discussed throughout the text because of its low recovery, as reported in previous studies [21,22].

Methods of Source Apportionment
PMF 5.0-the newest version of the PMF model released by the U.S. EPA-and an improved NMF model were used to apportion sources of gaseous PAHs in Harbin in this study. The PMF model was developed based on ME-2 algorithm and detailed information can be found in previous literature [23]. The improved NMF model has a similar function with the PMF model [23], namely, an original matrix (V) is decomposed into two matrices (W and H) with each coefficient no less than zero, which can be expressed as follows: where p is the number of factors and e is the simulated residual of the components in row i and column j. This model decomposes the matrix by minimizing the objective function (Q NMF ), which can be expressed as follows: where u is the matrix of the uncertainty and the other symbols have the same meaning as the symbols in Equation (1) [23,24]. The improved NMF model consists of four algorithms, which are a multiplicative update method (MU) [25], an optimal gradient method (OG) [26], a highly efficient monotonic fixed-point algorithm (FP) [27], and a conjugate gradient algorithm (CG) [28]. The four algorithms were modified by adding an uncertainty matrix to decompose a data matrix to adapt to the source apportionment [13]. The four algorithms in the improved NMF model were compiled using Matlab R2016b software [13].
The improved NMF model and the PMF model used the same input data and iterations in this study. The input data of PAHs concentration was a matrix of 68 × 15. The uncertainty matrix was made as follows: when the concentration of the jth PAH component of the ith sample was lower than the detection limit of the component, the corresponding element value was 5/6 times the detection limit in row i and column j of the uncertainty matrix; otherwise, the corresponding element value in the uncertainty matrix was calculated by the following formula: where Uc is the relative uncertainty of each PAH component (%), C is the concentration of PAHs component (ng/m 3 ), and MDL is the detection limit of PAH component (ng/m 3 ).
In order to maximize the chance of reaching the globally optimal solution and compare the results, the improved NMF model and the PMF model were run 100 times for a final solution by using a new random seed or starting point for each iteration [23]. Principal component analysis (PCA) was used to pre-estimate the number of sources before the source apportionments by using the improved NMF model and the PMF model. The assessment was because PCA is often used to analyze the interrelationships among a large number of variables and to explain these variables in terms of a smaller number of variables with a minimum loss of information [16]. SPSS 19.0 for Windows was used for the PCA analysis.

Toxic Potency Assessment of PAHs
The potential health risks of exposure to PAHs were assessed by the BaPeq concentrations. The total BaPeq concentration contributed by 15 PAH components were calculated by the equation as follows: where C i and TEF i are the measured concentration and the toxicity equivalent factor of ith PAH component. The TEF values of PAHs components used in this study were given by previous studies [16,29].
Among the six typical functional areas in Harbin, the Σ 15 PAHs concentrations in the wastewater treatment plant, industrial zone, business district, residential area, school and suburb were 98.1 ± 76.7 ng/m 3 , 91.2 ± 76.2 ng/m 3 , 71.4 ± 75.6 ng/m 3 , 67.9 ± 65.6 ng/m 3 , 42.6 ± 34.7 ng/m 3 , and 38.5 ± 38.0 ng/m 3 , respectively. The Kruskal-Wallis test showed that PAH concentrations had no significant differences among the six typical areas of Harbin. The high concentration of PAHs suggests that their sources were widespread; as reported, the wastewater treatment plant was a heavy emission source of PAHs [31,32] and the relatively high emission intensity of PAHs in the industrial zone, business district, and residential areas was caused by human activities in these areas [33][34][35][36]. Generally, the components of PAHs in the six areas had a similar pattern. The percentage of 3-ring PAHs to the Σ 15 PAHs (61.1%) in the wastewater treatment plant was lower than the average proportion, while the percentage of 4-ring PAHs (36.6%) was higher than the average level. Meanwhile, the proportion of 3-ring PAHs (75.8%) in residential area was higher than the average, while the percent of 4-ring PAHs (21.8%) was lower than the average. The slight difference may be due to differences in their sources [31][32][33][34][35]. The seasonal variation of gaseous concentrations of PAHs in the atmosphere of Harbin from summer 2016 to spring 2017 is listed in Table 2. The average concentrations of Σ 15 PAHs in Harbin in summer, autumn, winter of 2016, and spring of 2017 were 28.2 ± 26.2 ng/m 3 , 40.6 ± 23.1 ng/m 3 , 153.5 ± 38.1 ng/m 3 and 57.7 ± 10.3 ng/m 3 , respectively. Obviously, the concentration of Σ 15 PAHs was about 3-5 times higher in winter than in other seasons. The most significant difference of PAH concentrations between winter and summer was found in residential and suburban areas. The ratios of PAH concentrations in winter to summer were about 13 in the two areas. The Residential area is an area where the heating demand is relatively concentrated in winter, and there are many household heating cases in the suburbs. It is concluded that the combustion emission of heating fuels is the main reason for the significant increase of PAHs concentration in winter [3,36,37]. In China, about 30% of the population lives in areas that need heating in winter [38], indicating the significance of reducing PAH emissions from the heating sector.

Source Apportionment of PAHs in the Atmosphere
The PCA model developed for the atmospheric PAHs explained 89.7% of the total variance when 2 principal components (PC) were considered. PC1 (78.4% explained variance, Figures 2 and 3) gave high positive scores and loadings for all PAH components and PC2 (11.3% explained variance, Figures 2 and 3) clearly associated concentrations of Fla, Pyr, and Chr. The overwhelming interpretability of PC1 indicates that the concentration changes of PAHs in the city were controlled primarily by fewer significant emission sources [39], and the primary sources were not well apportioned by the PCA analysis [40]. Based on the pre-estimation results, the improved NMF model and PMF model were used to apportion sources of PAHs in the Harbin atmosphere. Three factors were identified after iterative testing from 3 to 5 factors in model exercises. The Q and Q(robust) values calculated by the four algorithms in the improved NMF model were the same (922.96 for Q and 901.55 for Q(robust), see Table 3). The Q and Q(robust) values estimated by the PMF model were 933.83 and 930.88, which were slightly higher than that calculated by the improved NMF model. The BS-DISP estimation by using the PMF model scenario with three factors showed that the largest decrease of Q(robust) were less than 1% of the Q value of corresponding base model. The minor difference indicated that the results of the base model could be treated as the global optimum solutions.   The three factors identified by the algorithms of MU, OG, FP, and CG in the improved NMF model and the PMF model showed similar patterns. Factor 1 was biomass burning, typically characterized by high loadings of low molecular weight PAHs such as Acy, Ace, Flu, Phe, and Ant [16,41]. Factor 2 was considered as vehicle exhaust emissions, characterized by high loads on high molecular weight PAHs such as BghiP, DahA, IDP, BaP, BkF, and BbF [16,41,42]. Factor 3 showed high proportions to the moderate molecular   The three factors identified by the algorithms of MU, OG, FP, and CG in the improved NMF model and the PMF model showed similar patterns. Factor 1 was biomass burning, typically characterized by high loadings of low molecular weight PAHs such as Acy, Ace, Flu, Phe, and Ant [16,41]. Factor 2 was considered as vehicle exhaust emissions, characterized by high loads on high molecular weight PAHs such as BghiP, DahA, IDP, BaP, BkF, and BbF [16,41,42]. Factor 3 showed high proportions to the moderate molecular  The three factors identified by the algorithms of MU, OG, FP, and CG in the improved NMF model and the PMF model showed similar patterns. Factor 1 was biomass burning, typically characterized by high loadings of low molecular weight PAHs such as Acy, Ace, Flu, Phe, and Ant [16,41]. Factor 2 was considered as vehicle exhaust emissions, characterized by high loads on high molecular weight PAHs such as BghiP, DahA, IDP, BaP, BkF, and BbF [16,41,42]. Factor 3 showed high proportions to the moderate molecular weight PAHs of Fla, Pyr, BaA, and Chr, which were treated as the contribution related to coal combustion [16,41,43].
The contributions of three factors to PAHs in the atmosphere apportioned by the improved model and the PMF model are shown in Figure 4. These results indicate that the models used in this study identified similar resolutions of source apportionments. The average contributions of biomass burning, coal combustion, and vehicle emission to Σ 15 PAHs concentration obtained by the two models were 48.6 ± 4.03%, 34.6 ± 3.22%, and 16.8 ± 5.06%, respectively. The molecular diagnostic ratios of PAHs, such as Ant/(Ant + Phe), Fla/(Fla + Pyr), BaA/(BaA + Chr), and IDP/(IDP + BghiP) ( Table 4), were often used to qualitatively differentiate sources of PAHs in past studies [44][45][46]. The ratios of Ant/(Ant + Phe), Fla/(Fla + Pyr), BaA/(BaA + Chr), and IDP/(IDP + BghiP) of PAHs in this study were 0.04 ± 0.03, 0.62 ± 0.03, 0.25 ± 0.08, and 0.82 ± 0.04, respectively. Among them, the ratios of Fla/(Fla + Pyr) and IDP/(IDP + BghiP) indicated that biomass burning and coal combustion were the main contributors to PAHs in the atmosphere in Harbin [44,45], the ratios of Ant/(Ant + Phe) and BaA/(BaA + Chr) suggested that vehicle exhaust was the main contributor of PAHs [44,47]. The similar sources identified by the improved NMF model and PMF model and the method of diagnostic ratios indicate that the models used in this study can comprehensively apportion PAHs sources.
Atmosphere 2021, 12, x FOR PEER REVIEW 9 of 13 weight PAHs of Fla, Pyr, BaA, and Chr, which were treated as the contribution related to coal combustion [16,41,43]. The contributions of three factors to PAHs in the atmosphere apportioned by the improved model and the PMF model are shown in Figure 4. These results indicate that the models used in this study identified similar resolutions of source apportionments. The average contributions of biomass burning, coal combustion, and vehicle emission to Σ15PAHs concentration obtained by the two models were 48.6% ± 4.03%, 34.6% ± 3.22%, and 16.8% ± 5.06%, respectively. The molecular diagnostic ratios of PAHs, such as Ant/(Ant + Phe), Fla/(Fla + Pyr), BaA/(BaA + Chr), and IDP/(IDP + BghiP) (Table 4), were often used to qualitatively differentiate sources of PAHs in past studies [44][45][46]. The ratios of Ant/(Ant + Phe), Fla/(Fla + Pyr), BaA/(BaA + Chr), and IDP/(IDP + BghiP) of PAHs in this study were 0.04 ± 0.03, 0.62 ± 0.03, 0.25 ± 0.08, and 0.82 ± 0.04, respectively. Among them, the ratios of Fla/(Fla + Pyr) and IDP/(IDP + BghiP) indicated that biomass burning and coal combustion were the main contributors to PAHs in the atmosphere in Harbin [44,45], the ratios of Ant/(Ant + Phe) and BaA/(BaA + Chr) suggested that vehicle exhaust was the main contributor of PAHs [44,47]. The similar sources identified by the improved NMF model and PMF model and the method of diagnostic ratios indicate that the models used in this study can comprehensively apportion PAHs sources.   The improved NMF model and the PMF model identified that contributions in summer, autumn, winter, and spring were 10.0 ± 3.57%, 32.5 ± 3.92%, 70.3 ± 4.43%, 12.5 ± 3.44% for biomass burning; 67.6 ± 5.33%, 57.9 ± 2.82%, 15.1 ± 3.46%, 55.2 ± 7.80% for coal combustion; and 22.4 ± 7.51%, 9.59 ± 3.54%, 14.7 ± 4.34%, 32.3 ± 10.1% for vehicle exhaust, respectively. It is already clear that the contribution of biomass burning increased significantly in winter and the source became the preponderant contributor of gaseous PAHs in the city during winter due to heating needs. The emission source of PAHs should be paid more attention in China because about 30% of the total population lives in areas that need heating in winter [38]. Coal combustion emitted more moderate molecular weight PAHs, such as Fla, Pyr, BaA, and Chr [16,41,43]. These moderate molecular weight PAHs are more favorable to partition into a gaseous form in the atmosphere under the condition of a high-temperature environment in summer [48]. The factors lead to coal combustion as the main contribution source of gaseous PAHs in Harbin in summer. Table 1 lists the ΣBaPeq concentrations calculated using Equation (4). The average ΣBaPeq concentration was 0.54 ± 0.23 ng/m 3 . The average ΣBaPeq concentrations in the industrial area, sewage treatment plant, business district, residential area, school, and suburb were 0.89 ± 0.77 ng/m 3 , 0.65 ± 0.54 ng/m 3 , 0.63 ± 0.78 ng/m 3 , 0.53 ± 0.51 ng/m 3 , 0.28 ± 0.23 ng/m 3 , and 0.23 ± 0.20 ng/m 3 , respectively. In terms of season, the average ΣBaPeq concentration reached its maximum value in winter, which was 1.20 ± 0.52 ng/m 3 , followed by spring (0.47 ± 0.13 ng/m 3 ), autumn (0.27 ± 0.13 ng/m 3 ), and summer (0.23 ± 0.12 ng/m 3 ), as listed in Table 2. Figure 4 displays the contribution of three sources to the ΣBaPeq concentration in Harbin identified by the improved NMF model and the PMF model. The contribution proportions of biomass burning, coal combustion, and vehicle exhaust were 28.7 ± 4.69%, 14.8 ± 2.80%, and 56.5 ± 6.17%, respectively. The contributions indicate that vehicle exhaust was the largest contributor to potential health risks of exposure to PAHs although the source contributed weakly to PAH concentration in the atmosphere in Harbin.

Conclusions
This study comprehensively investigated ∑ 15 PAHs in passive air samples collected from 17 sampling sites in six functional areas in Harbin, China from May 2016 to April 2017. The results showed that the average concentration of Σ 15 PAHs in 68 collected samples was 68.3 ± 22.3 ng/m 3 . It was dominated by 3-ring PAHs (64.4%), followed by 4-ring (32.6%), 5-ring (2.10%), and 6-ring (0.89%). The components of Phe, Fla, and Flu had the highest concentrations among the Σ 15 PAHs. The Σ 15 PAHs concentrations were 98.1 ± 76.7 ng/m 3 , 91.2 ± 76.2 ng/m 3 , 71.4 ± 75.6 ng/m 3 , 67.9 ± 65.6 ng/m 3 , 42.6 ± 34.7 ng/m 3 , and 38.5 ± 38.0 ng/m 3 in the wastewater treatment plant, industrial zone, business district, residential area, school, and suburb, respectively. The improved NMF model and the PMF model analyzed three source of PAHs. Among the three sources, the average contribution of biomass burning, coal combustion, and vehicle exhaust were 48.6 ± 4.03%, 34.6 ± 3.22%, and 16.8 ± 5.06%, respectively. The contribution of biomass burning was significantly increased in winter, while that of coal combustion was significantly higher in summer. Vehicle exhaust was the largest contributor to the potential health risks of exposure to PAHs in Harbin. As an important city with a population of over a million in northeast China, Harbin should not only pay attention to the atmospheric concentration level of PAHs and their sources, but also pay enough attention to other pollutants with significant impacts on health risks and their sources.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no competing financial interest.