Primary Sources of Polycyclic Aromatic Hydrocarbons to Streambed Sediment in Great Lakes Tributaries Using Multiple Lines of Evidence

Abstract Polycyclic aromatic hydrocarbons (PAHs) are among the most widespread and potentially toxic contaminants in Great Lakes (USA/Canada) tributaries. The sources of PAHs are numerous and diverse, and identifying the primary source(s) can be difficult. The present study used multiple lines of evidence to determine the likely sources of PAHs to surficial streambed sediments at 71 locations across 26 Great Lakes Basin watersheds. Profile correlations, principal component analysis, positive matrix factorization source‐receptor modeling, and mass fractions analysis were used to identify potential PAH sources, and land‐use analysis was used to relate streambed sediment PAH concentrations to different land uses. Based on the common conclusion of these analyses, coal‐tar–sealed pavement was the most likely source of PAHs to the majority of the locations sampled. The potential PAH‐related toxicity of streambed sediments to aquatic organisms was assessed by comparison of concentrations with sediment quality guidelines. The sum concentration of 16 US Environmental Protection Agency priority pollutant PAHs was 7.4–196 000 µg/kg, and the median was 2600 µg/kg. The threshold effect concentration was exceeded at 62% of sampling locations, and the probable effect concentration or the equilibrium partitioning sediment benchmark was exceeded at 41% of sampling locations. These results have important implications for watershed managers tasked with protecting and remediating aquatic habitats in the Great Lakes Basin. Environ Toxicol Chem 2020;39:1392–1408. © 2020 The Authors. Environmental Toxicology and Chemistry published by Wiley Periodicals LLC on behalf of SETAC.


INTRODUCTION
Representing 84% of the fresh surface water in North America (US Environmental Protection Agency 2015), the Great Lakes are an invaluable natural resource to the United States and Canada. However, a history of industrial, agricultural, and household pollution has left a legacy of contaminated sediments in many areas of the Great Lakes and their tributaries. Some of the contaminants, such as polychlorinated biphenyls (PCBs) and DDT, are primarily historical, because regulations in recent decades have resulted in major reductions in their source contributions. Other contaminants, such as polycyclic aromatic hydrocarbons (PAHs), have historical and modern sources, and therefore continue to enter and accumulate in the Great Lakes and their tributaries today (Baldwin et al. 2016(Baldwin et al. , 2017. A recent study of organic compounds in water samples from Great Lakes tributaries found that among 15 classes of organic contaminants, including herbicides and insecticides, PAHs posed the greatest risk to aquatic organisms (Baldwin et al. 2016).
The PAHs are a class of >100 organic compounds composed of 2 or more fused aromatic rings. They are widespread contaminants with sources both historical and modern, and both natural and anthropogenic. Petrogenic PAHs, which form at low temperatures over geologic time scales, come from refined petroleum products (asphalt, diesel, gasoline, home heating oil, motor oil, and lubricants) and unprocessed coal and crude oil, among others (Pietara et al. 2010). Pyrogenic sources, in contrast, form at high temperatures during incomplete combustion of carbon-based material (Pietara et al. 2010). Natural sources of pyrogenic PAHs include volcanic eruptions and wildfires. Anthropogenic sources of pyrogenic PAHs include residential wood burning, exhausts from diesel and gasoline engines, and emissions from coal-fired power plants and coke-ovens, creosote, and coal tar from pavement sealants and former manufactured gas plants (Mahler et al. 2005;Neff et al. 2005;Pietara et al. 2010).
A number of PAHs are carcinogenic, mutagenic, teratogenic, and/or toxic to aquatic organisms (Eisler 1987). As a result, the International Joint Commission has identified carcinogenic PAHs as a "critical pollutant" in the Great Lakes (Agency for Toxic Substances and Disease Registry 2008), making PAHs a priority for reduction and elimination efforts. The PAHs are listed as contaminants of concern at 61% of Great Lakes Areas of Concern (US Environmental Protection Agency 2013). To that end, the US Environmental Protection Agency (USEPA) and other groups have spent more than $550 million cleaning up contaminated sediment in the Great Lakes region since 2002, with a primary focus on PAHs, PCBs, and metals (US Environmental Protection Agency 2002). With such resources being devoted to sediment clean-up, it is important to understand the distribution, magnitude, and sources of PAHs in the Great Lakes Basin.
The objectives of the present study were to 1) assess the occurrence and potential adverse biological effects of PAHs in recently deposited sediments in Great Lakes tributaries across 6 US states, and 2) identify the most important sources of PAHs to Great Lakes tributaries using multiple lines of evidence. The lines of evidence were 1) land-use analysis, 2) parent/alkylated ratios, 3) high-molecular-weight/low-molecular-weight (HMW/ LMW) ratios, 4) PAH profiles, 5) principal component analysis (PCA), 6) positive matrix factorization (PMF) receptor modeling, and 7) mass fraction analysis.

Site selection
Streambed sediment samples were collected from June to July 2017 from Great Lakes tributaries in Minnesota, Wisconsin, Michigan, Indiana, Ohio, and New York (USA). Between 1 and 7 locations within 26 tributary watersheds were sampled for a total of 71 sampling locations (Table 1, Figure 1, and Supplemental Data, Table S1). Locations were selected to represent watersheds with a range of land uses, and were from 0.7 to 100% urban. Watershed drainage areas ranged from 3.5 to 16 300 km 2 , and population densities ranged from 2.8 to 2260 persons/km 2 .

Streambed sediment sample collection and analysis
Sample collection and analysis methods are detailed in the Supplemental Data. To summarize, streambed sediment sample collection was performed either by boat or while wading in the stream. Fine-grained sediments (silts) were targeted at all locations. Sediment was collected to a depth of 15 cm using a push core sampler (WaterMark ® Universal Core Head Sediment Sampler, Forestry Suppliers) with a 70-mm (2 3/4 inch) outer diameter, a 66.7-mm (2 5/8 inch) inner diameter, and 1.6-mm (1/16 inch) wall polycarbonate tubing (Forestry Suppliers; or United States Plastic). The depth of 15 cm was selected to focus on recently deposited sediments, and thus modern versus historical PAH sources. The sediment was emptied into a stainless-steel pan and split vertically; then one-half was transferred to a baked amber-glass bottle for use in the present study. The other half was used for a separate study that is not described in the present study. Samples were stored in the dark on ice. All sediment processing equipment was field-cleaned between sampling locations by scrubbing with detergent (Alconox ® ) water followed by 3 rinses each of tap and deionized water. A new core tube was used at each location. Laboratory detection limits for PAHs ranged from 0.17 to 70.5 µg/kg (median 0.48 µg/kg), varying by compound and by analytical batch. The detection limit for TOC was 0.05%. Only 5.4% of PAH results were below the detection limit, and no TOC results were below the detection limit. Zeros were used as conservative substitutes for PAH concentrations below the detection limit in summations of total sample concentrations.
Duplicate samples were collected at 8 locations, resulting in 288 duplicate pairs (8 duplicate pairs for each of the 36 compounds). The PAH concentrations in 22 of the duplicate pairs were below the detection limit in both samples; in 4 duplicate pairs, concentrations were below the detection limit in 1 of the 2 samples. Relative percentage differences (RPDs) were calculated for all remaining duplicate pairs (i.e., those with concentrations above the detection limit in both samples). The median RPD was <20% for 24 of the 36 PAH compounds. The compounds with the highest RPDs (medians of 20.7-42.4%) were generally those with the lowest molecular weights and occurring at the lowest concentrations (parent and alkylated naphthalene, parent and alkylated fluorene, acenaphthylene, and acenaphthene). The TOC duplicates (n = 5) had RPDs < 4%.
Field blanks were collected at 5 locations by pouring organic-free water (OmniSolv ® ) through a core tube into the stainless-steel pan, and then into a baked amber-glass jar. The majority of field blank PAH concentrations (71%) were below the detection limit; the 29% of blank concentrations above the detection limit ranged from 0.30 to 5.20 ng/L. There were no instances of blank concentrations above the detection limit when accompanying environmental sample concentrations were below the detection limit. All TOC blanks (n = 5) were

Predicted toxicity using sediment quality guidelines
Potential PAH-related toxicity of streambed sediments to aquatic organisms was assessed using the following sediment quality guidelines: the probable effect concentration (PEC; 22 800 µg/kg for ΣPAH 16 ), the threshold effect concentration (TEC; 1610 µg/kg for ΣPAH 16 ), and the sum equilibriumpartitioning sediment benchmark toxicity unit (ΣESBTU; Ingersoll et al. 2001;US Environmental Protection Agency 2003;Kemble et al. 2013). The PEC quotients (PECQs) and the TEC quotients (TECQs) were computed for each sample by dividing the ΣPAH 16 concentration in the sample by the PEC and TEC, with adverse effects to benthic organisms predicted at PECQs > 1.0, and unlikely at TECQs < 1.0 (Ingersoll et al. 2001).
The ΣESBTU approach accounts for the biological availability of individual PAH compounds in a mixture, and is applicable across sediment types (US Environmental Protection Agency 2003). To compute the ΣESBTU, TOC-normalized concentrations of 35 PAHs (listed in the Supplemental Data, Table S2) were divided by compound-specific final chronic values and summed. Streambed sediments with ΣESBTUs < 1.0 are expected to be nontoxic to benthic organisms from PAHs, whereas sediments with ΣESBTUs > 1.0 are expected to have adverse effects from PAHs.
Toxicity related to contaminants other than PAHs was not considered in the present study.   Sources of PAHs to Great Lakes tributaries-Environmental Toxicology and Chemistry, 2020;39:1392-1408

Geographic information system methods
Area-Characterization Toolbox (Price 2017), which allows for standard GIS summaries in nested basins (Price 2017). Specifically, mean 2011 percentage imperviousness and mean 2012 parking lot abundance were calculated using the Feature Statistics to Table tool, and 2012 land-use category percentages were determined using the Tabulate Features to Percent tool (Falcone 2015;Homer et al. 2015;Falcone and Nott 2019). Land-use data were further summarized into major areas: transportation (code 21 in the NAWQA Anthropogenic Land Use Trends Dataset; Falcone 2015), commercial (code 22), industrial/military (code 23), residential (codes 25-26), and total urban (codes 21-27). Notably, the calculated parking lot abundance in each basin was derived from (and therefore overlaps with) the land-use categories. For example, an area categorized as commercial likely includes (and does not exclude) adjacent parking areas.
Watershed population summaries were determined by creating a mosaic of 2010 state census block polygons, calculating the population density/census block polygon, intersecting census block and watershed boundaries, calculating the areas of the intersected polygons, multiplying these areas by the population density, and summing the result by watershed (US Census Bureau Geography Division 2010a, 2010b, 2010c, 2010d.

Identification of PAH sources
Multiple lines of evidence were used to identify the most likely source of PAHs to sediment samples. This approach mitigates uncertainties of individual methods and strengthens the overlapping conclusions (Larsen and Baker 2003;O'Reilly et al. 2014). The individual methods used were described previously (Baldwin et al. 2017), and are briefly described in the present study. Many of the methods used a subset of the 36 PAHs analyzed. The number of PAHs used and the reason for subsetting are described for each method.

Land-use analysis
Relations between different urban land uses and streambed sediment ΣPAH 16 concentrations were assessed by using Spearman correlation with a significance level (p value) of 0.05. Spearman correlation results were unaffected by treatment of results below the detection limit: substitutions with 0, ½ × detection limit, and 1 × detection limit were each tested, and all yielded the same correlation values.

Parent/alkylated and HMW/LWM compounds
Ratios of the mass concentrations of parent and alkylated PAHs were used to differentiate between petrogenic and pyrogenic sources. Petrogenic sources are generally dominated by alkylated compounds, whereas pyrogenic sources are generally dominated by parent compounds (Neff et al. 2005). Parent/alkylated ratios were computed using only the compounds for which both the parent and alkylated forms were measured. The 8 parent compounds included for this analysis were benz[a]anthracene, chrysene, pyrene, fluoranthene, fluorene, naphthalene, anthracene, and phenanthrene. The 18 alkylated compounds were the C1 to C3 or C1 to C4 alkylated forms of the 8 parents (Supplemental Data, Table S2).
Ratios of the mass concentrations of LMW (2-3 rings) and HMW (>3 rings) PAHs (Supplemental Data, Table S2) were also used to differentiate between petrogenic and pyrogenic sources. Petrogenic sources are generally dominated by LMW compounds, whereas pyrogenic sources are generally dominated by HMW compounds (Crane 2014 perylene. The proportional concentration of each compound was calculated as the fraction of the ΣPAH 12 concentration; each 12-compound profile was summed to 1.0. Samples with concentrations below the detection limit were excluded from analysis of PAH profiles (and PCA, described in the PCA section following) because concentrations of all compounds are necessary to properly characterize the PAH profiles for these analysis methods.
The sample profiles were compared quantitatively with 12-compound proportional concentration profiles of different sources from the literature (Table 2 and Supplemental Data,  Table S3). Benzo[e]pyrene was not reported in the profiles of creosote-treated railway ties (representing weathered creosote) and creosote product (unweathered). For those profiles, following Li et al. (2003), the benzo[e]pyrene concentration was assumed to be the same as that of benzo[a]pyrene. The similarity between source and sample profiles was evaluated using the chi-square (χ 2 ) statistic, calculated as the square of the difference in proportional concentrations of individual compounds, divided by the mean of the 2 values, summed for the 12 PAHs (Van Metre and Mahler 2010). A lower χ 2 indicates greater similarity between source and sample profiles. A profile was not computed for one site, Kalamazoo River at New Richmond, MI (MI-KAL), because of concentrations below the detection limit.

PCA
The PCA was performed using the same 12-compound PAH profiles just discussed, with data standardized to have a mean of 0 and unit variance. Euclidean distances in n-dimensional space were computed between sources and samples in the space defined by the principal components that accounted for ≥10% of the variability. Sources with the shortest Euclidean distance to the samples were considered to be most similar to the samples. The PCA computation was done using the prcomp function from the stats package in R (R Core Development Team 2015).

PMF receptor model
The PMF is a multivariate receptor modeling tool that decomposes a matrix of speciated sample data into 2 matricessample contributions and factor profiles (Norris et al. 2014;Brown et al. 2015). The theory and detailed methods of PMF have been described previously (Paatero and Tapper 1994;Paatero 1997;Norris et al. 2014). The PMF was run using USEPA PMF Ver 5.0.14, a graphical user interface for the Multilinear Engine Program ME-2. Unlike receptor models based on weighted linear regression (e.g., the USEPA Chemical Mass Balance Model), which require selection of the sources contributing to a sample prior to the regression analysis, the USEPA PMF determines profiles based on sample concentrations and the number of sources or factors specified by the user. The PMF does not identify the factor names; it is up to the user to match the factor profiles to known sources using measured source profiles employing approaches like the chi-square analysis we describe.
A strength of the PMF model is that it considers the uncertainty of individual compound concentrations and the total sample concentration. The uncertainty of individual compound concentrations was computed using Equation 1 (Qi et al. 2016): where DL is the detection limit. The error term (i.e., measurement error) was calculated as the median relative percentage difference between duplicate samples. The error term was therefore specific to each compound, varying from 0.12 to 0.20. The uncertainty on the total sample concentration (the combined standard uncertainty, U c (x)) was computed as the root sum of the squares of the individual compound uncertainties (Equation 2): No extra modeling uncertainty was included in the PMF analysis.
The model was run in robust mode with 29 compounds (species), all of which were considered "strong" based on signalto-noise ratios > 1.0 (concentrations were ≥2 × uncertainties). Two compounds, C4-chrysenes and C3-fluorenes, were omitted from the PMF analysis because of low detection frequencies. In addition, naphthalene and its alkylated forms (C1-4) were omitted because they were poorly predicted by PMF (observed vs predicted r 2 values of 0.004-0.23), likely because of their low molecular weight and high volatility. The omission of naphthalenes reduced (i.e., improved) the model's Q value (a goodness-of-fit parameter) by 1344. A total sample concentration variable was also included as a species but was downweighted to "weak" (Norris et al. 2014).
Three samples (Indiana Harbor Canal at East Chicago, Indiana [IN-IHC], Geddes Brook at Fairmount, NY [NY-GBF], and Cuyahoga River at Independence, OH [OH-CRI]) were omitted from the PMF analysis because their unique PAH sources were not similar to the other samples. These sites had high scaled residuals for one or more PAHs and/or abnormally high concentrations. Samples with one or more concentrations below the detection limit were excluded from the PMF analysis for consistency with other methods used in the present study (i.e., profile correlations and PCA), which resulted in 14 samples being omitted. A total of 54 samples were included in the final PMF runs. The final PMF input files are provided in the Supplemental Data, Tables S4 and S5.
Multiple PMF runs were performed to determine the best solution, varying the number of factors (2-4) and the convergence criteria (default vs relaxed; see the Supplemental Data for discussion of relaxed convergence criteria) between runs. Each run consisted of 50 base runs with a random seed value. The best solution within each run was the one that minimized Q robust , a goodness-of-fit statistic calculated by dynamically downweighting points for which the uncertaintyscaled residual was >4.0 (Paatero 1997). The solution was then tested for rotational ambiguity (the existence of multiple similar solutions, which may invalidate the chosen solution) using the USEPA PMF's displacement analysis, and for disproportionate effects of a small set of observations using the USEPA PMF's bootstrapping analysis. Bootstrapping was performed with 100 runs, with a block size of 2 and a minimum r value of 0.6. The final solution was the one that maximized the number of factors to account for the most important sources while still Sources of PAHs to Great Lakes tributaries-Environmental Toxicology and Chemistry, 2020;39:1392-1408 minimizing rotational ambiguity and disproportionate effects of small sets of observations. The final solution was further tested for rotational ambiguity by adding constraints using ratios derived from the 12-compound source profiles from the literature (Supplemental Data, Table S3). Constraints are further discussed in the Supplemental Data. Finally, to assess potential bias in the goodness-of-fit toward low or high concentration samples, scaled residuals were plotted against ΣPAH 29 .
To match each of the PMF factors to a specific PAH source, the factors were compared with source profiles from the literature (Table 2 and Supplemental Data, Table S3) using the χ 2 statistic. Although PMF was run using 29 compounds, the χ 2 was computed using only the 12 compounds available in the literature source profiles.

Mass fractions
Mass fractions analysis has been used as a tool to eliminate potential PAH sources as the primary source to environmental samples because of unlikely or even impossible mass fractions of source material required to achieve environmental PAH concentrations (Ahrens and Depree 2010;Baldwin et al. 2017).
The mass fraction for each sample/source combination was calculated by dividing the ΣPAH 16 concentrations in samples by mean PAH concentrations of potential sources gathered from the literature ( Table 3). The ΣPAH 16 concentrations were used (rather than ΣPAH 36 ) because ΣPAH 16 is commonly reported in the literature, and thus the number of source concentrations included in the analysis could be maximized. The mass fraction therefore represented the hypothetical percentage mass of source material in a given sediment sample, assuming negligible contributions from other sources. We also assumed that PAHs are confined to the organic carbon fraction of the sediment sample, thereby providing an upper limit on the amount of source material (in percentage mass) possible in a given sample. If the source mass fraction was greater than the percentage of TOC in the sample, it was determined that the source could not possibly be the primary contributor of PAHs to the sediment (Table 4). If the mass fraction was <TOC but >½TOC, it was determined that the source was possible but unlikely to be the primary contributor of PAHs to the sediment, based on the assumption that >50% of the mass of organic carbon in a sediment sample is likely comprised of materials other than PAHs. If the mass fraction was <½TOC, it was determined that the source could be a primary contributor of PAHs to the sediment.

Observed concentrations relative to sediment quality guidelines
Concentrations of ΣPAH 16 in streambed sediment ranged from 7.4 to 196 000 µg/kg (mean 13 300 µg/kg; median 2600 µg/kg; Figure 2 and Table 5). Concentrations were especially high at  2 sites, NY-GBF and IN-IHC, where ΣPAH 16 was >2 × higher than at other sites. The percentage of TOC ranged from 0.1 to 8.8, with a mean and median of 1.7 and 1.2, respectively. All PAH and TOC results are provided in the Supplemental Data ,  Table S6. The TEC was exceeded in 62% of samples, with a median TECQ of 1.6 ( Table 5). The PEC was exceeded in 18% of samples, with a median PECQ of 0.1. The highest PECQs were at NY-GBF (8.6) and IN-IHC (5.9). The ΣESBTU benchmark value of 1.0 was exceeded in 38% of samples, with a median ΣESBTU of 0.4. The highest ΣESBTUs were at NY-GBF (10.5), Tributary to Buck Creek at Wyoming, MI (MI-TBC; 5.2), and Lower River Rouge at Wayne, MI (MI-LRW; 5.1).

Identification of PAH sources
Land-use analysis. The 6 urban land-use categories (percentage area of parking lot, major transportation, commercial, industrial/military, residential, and urban total), percentage impervious, and population density were all significantly related to sediment ΣPAH 16 concentrations (p < 0.05). Population density was the most strongly correlated with ΣPAH 16 concentrations, with a Spearman's rank correlation coefficient (r) of 0.66 (Supplemental Data, Figure S1), followed by parking lot land use (r = 0.64), percentage impervious (r = 0.63), commercial land use (r = 0.62), and total urban land use (r = 0.61). The categories least correlated with ΣPAH 16 concentrations were residential land use (r = 0.55), industrial/ military land use (r = 0.48), and major transportation land use (r = 0.48).
Parent/alkylated and HMW/LWM compounds. The ratio of parent to alkylated compounds was >1.0 in 78% of the 71 streambed sediment samples, indicating a pyrogenic PAH source (Neff et al. 2005). The median ratio of parent/alkylated compounds was 2.4 across all sites (Table 5).
The HMW compounds were dominant over the LMW compounds in 90% of the streambed sediment samples, another indicator of a dominantly pyrogenic source of PAHs to most samples (Crane 2014). The median ratio of HMW/LMW compounds was 5.4 across all sites.
Seven samples had HMW/LMW and parent/alkylated ratios <1.0. Unlike the majority of the samples in the present study, the primary source of PAHs to these 7 samples was likely petrogenic. Four of the 7 samples were from the Rocky River basin (OH-WBR, OH-RRB, OH-RRS, and OH-EBR) where especially thin zones of depositional sediment may have resulted in inadvertent sampling of bank material; the remaining 3 samples were from the Cuyahoga (OH-TCD), Maumee (OH-SCE), and Rouge (MI-LRH) River basins. The ΣPAH 16 concentrations of these samples were generally low, ranging from 81.0 to 2010 µg/kg (median 315 µg/kg). PAH profiles. The 12-compound PAH profiles of the individual streambed sediment samples were generally similar, especially among those samples with a pyrogenic signature (i.e., those with ratios of parent/alkyl compounds ≥1.0 and/or  Table 1. Watersheds are ordered by maximum PAH concentration; sites are ordered upstream to downstream within each watershed. PEC = consensus-based probable effect concentration; TEC = consensus-based threshold effect concentration. HMW/LMW compounds ≥1.0; Figure 3 and Comparisons with 12-compound PAH source profiles from the literature showed that sediment samples were most similar to the profile of coal-tar-sealed pavement dust (CTD7), with a median χ 2 statistic of 0.07 (Figure 4). Other sources with profiles similar to sediment samples included vehicle/traffic average (VAVG; median χ 2 0.11), traffic tunnel air (TUN1; median χ 2 0.15), coal combustion average (CCB1; median χ 2 0.15), and pine combustion #2 (PIN2; median χ 2 0.17). Sources with profiles least similar to sediment samples included oak combustion (OAKS; median χ 2 0.98), creosote product (CRE4; median χ 2 0.88), used motor oil #1 and 2 (UMO1 and 2, median χ 2 0.74 and 0.56, respectively), creosote-treated railway ties (CRE2; median χ 2 0.68), tire particles (TIRE; median χ 2 0.65), and residential heating (RESI; median χ 2 0.58).
PCA. Principal components 1 to 4 each explained >10% of the total variance in the dataset and together explained 80% of the total variance. Consistently, Euclidean distances between streambed sediment samples and PAH sources (computed for all paired combinations of principal components 1-4; plotted in the Supplemental Data, Figure S2), identified coal-tar-sealant pavement dust (CTD7) as the most similar source to streambed sediment samples ( Figure 5). Other sources showing similarity to streambed sediment samples in some principal component combination graphs were vehicle/traffic average (VAVG), pine combustion #2 (PIN2), and coal combustion average (CCB1). Sources most distant from streambed sediment samples (i.e., largest Euclidean distances) were creosote-treated railway ties  Table 1. TOC = total organic carbon; ΣPAH 16 = sum concentration of US Environmental Protection Agency 16 priority pollutant PAH compounds; HMW = high molecular weight; LMW = low molecular weight; PECQ = consensus-based probable effect concentration quotient; TECQ = consensus-based threshold effect concentration quotient; ΣESBTU = sum equilibrium partitioning sediment benchmark toxicity units; NA = not measured or computed.   Table S7). The 3-factor solutions (with default and relaxed convergence criteria) maximized the number of factors without violating model diagnostics recommendations described in the USEPA PMF User Guide (Norris et al. 2014) and listed in the Supplemental Data, Table S7. Relaxing the convergence criteria had little impact on the diagnostics of the 3-factor solution (Supplemental Data, Table S7), so the default convergence criteria were used in the final solution. Scaled residuals plotted against ΣPAH 29 indicated no bias in the goodness-of-fit of the final solution (Supplemental Data, Figure S5). The proportional concentration profiles of factors 1 to 3 in the final 3-factor PMF solution were generally similar to one another (Supplemental Data, Figure S6), with some exceptions noted below. Factor 1 contributed 48% of the total PAHs. The 12-compound profile of factor 1 was similar to 2 sources, coaltar-sealed pavement dust (CTD7; χ 2 0.064) and vehicle/traffic average (VAVG; χ 2 0.069; Table 6). Because the χ 2 value of coaltar-sealed pavement dust was only slightly lower than that of vehicle/traffic average (VAVG; χ 2 difference of 0.005), factor 1 is interpreted to be a mixture of the 2 sources (and possibly others). Factor 2 contributed 31% of the total PAHs and, compared with factor 1, had slightly higher proportions of HMW compounds (Supplemental Data, Figure S6). The 12-compound profile of factor 2 was most similar to that of coal-tar-sealed pavement dust (CTD7; χ 2 0.133), followed by traffic tunnel air (TUN1; χ 2 0.163; Table 6). Compared with factor 1, the wider gap (0.03) in χ 2 values between the top 2 sources lends more confidence in attributing factor 2 to a single source: coaltar-sealed pavement dust. Factor 3 contributed 21% of the total PAHs and, compared with factors 1 and 2, had higher proportional concentrations of phenanthrene and alkylated phenanthrenes/anthracenes. The 12-compound profile of factor 3 was most similar to that of coal-tar-sealed pavement dust (CTD7; χ 2 0.042) followed by vehicle/traffic average (VAVG; χ 2 0.082; Table 6). As with factor 2, the relatively large difference (0.04) in χ 2 values between the top 2 sources provides some confidence in attributing factor 3 to coal-tar-sealed pavement dust. Because the χ 2 analysis of factors 2 and 3 attributed both factors to coal-tar-sealed pavement dust, their contributions were combined. Thus, based on the PMF results, coal-tar-sealed pavement dust was the dominant PAH source to 70% of samples (excluding samples omitted from the PMF analysis), contributing an average of 68% of the total PAHs to each sample (minimum 16%, median 68%, maximum 100%). The rest of the PAHs were from a mixture of sources, including coal-tar-sealed pavement dust and vehicle/traffic average.
The χ 2 statistics between PAH sources and some of the other model runs (i.e., the 3-factor model with relaxed convergence FIGURE 5: Euclidean distances between sources and samples for principal component analysis components 1 through 4 using 12compound polycyclic aromatic hydrocarbon (PAH) profiles. Boxes = 25th to 75th percentiles; dark line = median; whiskers = 1.5 × the interquartile range (IQR); circles = value outside the 1.5 × the IQR. criteria, and the 2-factor models with default and relaxed convergence criteria) are provided in the Supplemental Data, Table S8. Contributions to individual sediment samples were estimated for each PMF factor in the final 3-factor solution ( Figure 6).

Mass fractions. For most sediment samples, mass fractions
analysis considerably narrowed the list of potential primary PAH sources. Many potential primary PAH sources to streambed sediment samples had PAH concentrations in the lowest approximately 10th percentile (Figure 7), but the number of potential primary sources diminished rapidly with increasing sample concentrations. The PAH sources such as asphalt-sealed pavement dust, wood combustion, and traffic tunnel and road dust could not be the primary sources to the upper approximately  (Table 2 and Supplemental Data, Table S3). For site abbreviations, see Table 1. CT dust = coal-tar-sealed pavement dust.

DISCUSSION
Multiple lines of evidence were used to determine the most likely source of PAHs to sediment samples from Great Lakes tributaries. The results of the individual source identification methods for each site are summarized in Figure 8. The likely dominant sources to each site determined by using PCA, profiles analysis, and the PMF model are listed. For PMF, the dominant source was coal-tar-sealed pavement dust (CTD7) where the sum contribution of factors 2 and 3 exceeded 50%, and the dominant source was a "mix" where the contribution of factor 1 exceeded 50%. Sources determined to be impossible based on mass fractions analysis are struck-through. Gray bars indicate parent/alkyl ratios and HMW/LMW ratios, with values <1.0 indicative of a petrogenic source, and values >1.0 indicative of a pyrogenic source. For samples with majority agreement (i.e., >50% of the identification methods agree) between the multiple lines of evidence, the most likely primary PAH source is identified in the rightmost column ("weight-ofevidence top source"). The weight-of-evidence top source has a check mark for samples with unanimous agreement across the multiple lines of evidence, indicating greater confidence. Some of these lines of evidence are more powerful than others; although assigning weights to them would be subjective, the weakest lines of evidence are likely the parent/alkyl and HMW/ LMW ratios, and the most powerful is likely mass fractions.
For 35 of the sampled sites (49%), unanimous agreement across all lines of evidence indicated that coal-tar-sealed pavement dust was the most likely primary source of PAHs. At an additional 22 sites (31%), coal-tar-sealed pavement dust was identified as the most likely primary source by the majority of (but not all) methods. At some sites, some portion of the coal-tar signature may have been from former manufactured gas plant contamination. Sampling surficial sediment from tributary streambeds was meant to minimize historical contributions, but it is possible that historical sediments were reworked and deposited at the surface.
Vehicle emission-related sources (VAVG and DVEM) were identified as the most likely sources at 5 sites (10%). At 8 sites there was no majority agreement on the most likely primary source. Creosote, despite having a very high PAH concentration, was not identified as the likely primary source at any site because its unique PAH profile differed considerably from the profiles of sediment samples. Likewise, coal combustion, a common PAH source in urban areas, was not identified as a primary source at any site.
The land-use analysis lends some support to the conclusion that coal-tar-sealed pavement dust was the most likely primary source of PAHs. Parking land use was found to better correlate with sediment PAH concentrations than major transportation land use (i.e., major roads), despite the similarities between these 2 categories: they are impervious surfaces made from materials that accumulate tire and brake particles, motor oil, exhaust from diesel and gasoline engines, and atmospheric deposition of PAHs. An important difference between parking areas and major roads is that pavement sealants are commonly used on parking areas (if made of asphalt) but are not typically used on major roads.
The methods used in the present study provide an evidence-based approach for identifying the most likely sources contributing to each sample, but each method has limitations and uncertainties including PAH source and sediment concentrations, inability of source profiles from the literature to capture variability in PAH sources in the study region (some of the source profiles are decades old or from other countries), variability in data quality in literature source profiles, the potential misinterpretation of results if the analysis lacks an important PAH source, and the potential for weathering to affect PAH profiles in sediment samples (Baldwin et al. 2017). The lack of consensus among the methods for 51% of the sites highlights these uncertainties and the value of using multiple lines of evidence, which mitigates uncertainties of individual methods and strengthens the overlapping conclusions (Larsen and Baker 2003;O'Reilly et al. 2014).
There is a need for an updated, comprehensive set of source profiles to improve future source-identification studies. The present study used 12-compound source profiles gathered from several different studies, with varying and often unknown measurement uncertainties. An updated list of source profiles, collected and analyzed using consistent methods, would provide a better understanding of uncertainties in the profiles and on the analyses reliant on the profiles. Including a greater number of compounds (≥16) may help differentiate between sources with similar profiles.
A limitation of our study was the assumption that a single sample was used to represent the PAH concentration and profile at each location. The PAH concentrations in streambed sediment are not spatially homogenous, as illustrated by the duplicate samples in the present study, which had median RPDs of up to 42.4% for individual compounds. However, despite the variability in concentrations between duplicate field samples, the PAH profiles of duplicate samples were quite similar (Supplemental Data, Figure S7). In fact, the PAH profiles were similar not only within a site, but at most of the sampled locations across the Great Lakes Basin (Figure 3). This finding suggests that, although multiple samples at each location would have shown a potentially wide range in concentrations, the PAH profiles may not have differed substantially.
Coal-tar-sealed pavement dust has been identified as the likely primary source of PAHs to streambed sediments elsewhere in the central and eastern United States, including in urban and suburban lakes, streams, and stormwater ponds in Austin (TX; Mahler et al. 2005  For PMF, the dominant source was coal-tar-sealed pavement dust (CTD7) where the sum contribution of factors 2 and 3 exceeded 50%, and the dominant source was "mix" where the contribution of factor 1 exceeded 50%. Sources determined to be impossible based on mass fractions analysis are struck-through. Gray bars indicate parent/alkyl ratios and high-molecular-weight to low-molecular-weight (HMW/LMW) ratios, with values <1.0 indicative of a petrogenic (PETRO) source, and values >1.0 indicative of a pyrogenic (PYRO) source. The "weight-of-evidence top source" (right column) is identified where >50% of the identification methods agree, with a check indicating unanimous agreement. NA = ≤50% of identification methods agree. Site MI-KAL is omitted because of concentrations below the detection limit. Site abbreviations are defined in Table 1 The ubiquity of the coal-tar-sealed pavement dust profile in sediments across the Great Lakes Basin and elsewhere raises the question of whether that profile may actually represent "urban background" (i.e., the mixture of common urban PAH sources such as coal and wood combustion and vehicle emissions). However, traditional urban background sources cannot account for the high PAH concentrations measured in many of the samples in our study. This point was demonstrated using mass fractions analysis and is further supported by simply comparing PAH concentrations from the present study with concentrations in areas where coal-tar sealants are not used (i.e., outside of the eastern and central United States and Canada). Compared with a ΣPAH 16 mean of 13 300 µg/kg and maximum of 196 000 µg/kg for the present study, studies of urban streams, canals, drains, and lakes in Portland  Boonyatumanond et al. 2006) reported mean ΣPAH 16+ concentrations of 663 to 5570 µg/kg (Figure 9). The maximum ΣPAH 16+ concentrations in these studies were all <9000 µg/kg, with the exception of Delhi storm drain sediments, which had a maximum of 19 300 µg/kg. Although not a comprehensive examination of PAH concentrations in urban sediments worldwide, these studies indicate that ΣPAH 16 concentrations are typically <9000 µg/kg where the primary PAH sources are urban background related. Even in the end-member case of Delhi storm drains, which likely approach the urban background maximum, ΣPAH 16 concentrations do not exceed 20 000 µg/kg. Thirty-two percent of the sites in the present study had ΣPAH 16 concentrations >9000 µg/kg, and 18% had ΣPAH 16 concentrations >20 000 µg/kg. Based on these comparisons, either 1) urban background PAH concentrations are substantially higher in the Great Lakes Basin than in these other locations in China, India, Thailand, and elsewhere; or 2) another PAH source, unrelated to traditional urban background sources and not present in these other locations, is elevating PAH concentrations above typical urban background concentrations in Great Lakes tributaries. A study of global atmospheric PAH emissions from coal, petroleum, and biofuel consumption and transformation (i.e., traditional urban background sources) reported an annual atmospheric PAH emission density (annual emission rate/land area) of 4.3 kg/km 2 /y -1 for the United States (excluding Alaska and Hawaii; Zhang and Tao 2009). In comparison, the reported atmospheric PAH emission density for China was 12.2 kg/km 2 /y -1 , and for India 30.2 kg/km 2 /y -1 . Given the much lower PAH emission density in the United States compared with China and India, it seems logical that sediments in the United States would tend to have lower rather than higher PAH concentrations. The fact that PAH concentrations in our study commonly exceeded the maximum concentrations measured in urban sediments in China and India suggests that a different PAH source, unrelated to traditional urban background sources and not widely present in these other locations, is elevating PAH concentrations in Great Lakes tributaries. Based on the multiple lines of evidence in the present and previous studies (Mahler et al. 2005;Van Metre and Mahler 2010;Watts et al. 2010;Yang et al. 2010;Pavlowsky 2013;Crane 2014;Baldwin et al. 2017;Valentyne et al. 2018), that source is most likely coal-tarsealed pavement dust.
Supplemental Data-The Supplemental Data are available on the Wiley Online Library at https://doi.org/10.1002/etc.4727. did not collect, generate, evaluate, or use the environmental data described in the present study.
Data Availability Statement-Data are provided in the Supplemental Data and are also available online at https:// doi.org/10.5066/F7P55KJN and at https://waterdata.usgs. gov/nwis.