Long-term trends (2005–2016) of source apportioned PM2.5 across New York State

A R T I C L E I N F O


Introduction
Starting in the 1970s, many developed countries including the United States implemented air quality management policies to protect human health and the environment. Several mitigation strategies were adopted across the United States to improve air quality (AQ), including implementation of legislation and regulations and the requirement of increasingly stringent emissions abatement technologies and standards. Since the 1970s, many air pollutant concentrations have declined across the United States (e.g., Hogrefe et al., 2004;Frost et al., 2006;Tagaris et al., 2007;Pye et al., 2009;Dallmann and Harley, 2010;Parrish et al., 2011;Brown-Steiner et al., 2016;Duncan et al., 2016;Nopmongcol et al., 2016;Rattigan et al., 2016;Emami et al., 2018;Masiol et al., 2018a;Squizzato et al., 2018a). Additionally, economic drivers have also contributed to decreases in air pollutant concentrations over the past decades (e.g., Tong et al., 2016;Squizzato et al., 2018a).
New York State (NYS) has had generally downward trends for several gaseous criteria pollutants (nitrogen oxides, SO 2 , CO) and PM 2.5 mass concentration across 54 monitoring sites. However, ozone concentrations showed slight increases . Downward trends for the particle number concentrations measured in Rochester have also been reported .
Currently, the PM 2.5 National Ambient Air Quality Standard (NAAQS) for PM 2.5 have been attained across most of the U.S. However, there were non-attainment areas in California, Idaho, Ohio, and Pennsylvania (US EPA, 2018b). Despite the overall improvement in air quality and the attainment of the NAAQS, ambient PM 2.5 pollution was still estimated to cause 88,400 (95% confidence interval = 67,000-115,000) deaths in 2015 across the United States (Cohen et al., 2017). Recent population-based studies have reported that even exposures to low PM 2.5 mass concentrations, both acute and chronic, may contribute to increased morbidity and mortality (Shi et al., 2015). Recent work by Zhang et al. (2018) found that as PM 2.5 has decreased across NYS, rates of hospitalizations and emergency department visits for cardiovascular diseases have decreased. However, in recent years, the rates of hospitalizations per unit mass of PM 2.5 have increased suggesting increased toxicity of the PM 2.5 with respect to certain diseases such as ischemic heart disease and congestive heart failure.
A challenge currently faced by policy-makers is the identification of additional steps to further improve AQ and reduce risks to public health. However, before new policies can be designed and implemented, there is the need to determine if: (i) past and current regulations and resulting actions have produced detectable and successful AQ improvements, and (ii) the past trends in air pollution represent feedbacks of specific mitigation measures or the synergy of implementation of multiple policies at different scales rather than being a reflection of recent economic events and drivers, and/or the direct/ indirect effects of climate change . These tasks are challenging for the following reasons: • Economic changes. Over time, different policies may have differently affected many sectors (e.g., manufacturing, power generation, road transport, maritime shipping, nonroad activities, fuel sulfur content, etc.), so that the effects of changes in individual sectors cannot be resolved; • Time duration. Many changes were phased-in with different timings and/or multiple programs (by introducing increasingly stringent standards). In addition, some mitigation measures experienced delays or were planned to roll out over long periods (e.g., Clean Air Interstate Rule, Cross-States Air Pollution Rule; nonroad diesel fuel sulfur) due to required changes in production processes, subsequent extensions to complete the changes, credit flexibilities, and time for fleet turnover. Consequently, effects on AQ evolved over the time.
• Target pollutants. Policies have mostly affected the emissions of PM precursors, i.e. SO 2 , nitrogen oxides (NO + NO 2 =NO x ), and volatile organic compounds (VOCs), rather than primary emissions that had been largely controlled in the past. Thus, the effects on AQ at a given site reflect the interplay of changes in local emissions and regional/ transboundary transport; • Source signatures. The effects of sources that emit PM 2.5 with similar chemical fingerprints or the same aerosol precursors may make it more difficult to resolve the trends for specific sources; • Changing chemistry. Changes in emission scenarios at one time point may produce changes in the atmospheric chemistry that may be taken into account in future AQ actions. For example, the sharp decrease in SO 2 emissions in the past decade has resulted in a significant decline in acidic sulfate available to react with ammonia, and producing a shift of the NH 4 -NO 3 equilibrium toward increased particulate nitrate as observed in the U.S. and Europe, and forecasted for China (e.g., West et al., 1999;Erisman and Schaap, 2004;Pye et al., 2009;Wang et al., 2013).
Thus, trends in PM sources need to be investigated using long time series of data. To do so, the contributions of individual PM 2.5 source types need to be resolved and quantified by chemometric methods that solve the mixture resolution problem, i.e. able to identify primary and secondary sources. Source apportionment (SA) using receptor models (RMs) can identify the major sources of air pollution and determine their contributions (Hopke, 2016). The U.S. EPA first included RMs as part of State Implementation Plans (SIPs) approximately 3 decades ago (US EPA, 2018a). Positive matrix factorization (PMF) provides a rigorous approach to apportionment through its weighted, non-negative least squares fitting to the data and error analyses to assess measurement error and rotational ambiguity . It has become the most widely used RM since the release of the freely-available EPA versions in the mid-2000s (Hopke, 2016).
The PMF identification and quantification of the major sources of PM 2.5 mass concentrations between 2005 and 2016 was performed for the 8 sites across NYS for which chemical speciation data are available. Those results were presented by Squizzato et al. (2018b). In the current work, the source apportioned PM 2.5 mass concentrations were used to investigate the trends in the PM 2.5 sources. The identified trends have been interpreted based on changes in AQ regulations across the U.S. and in New York as well as economic drivers of changing emissions. Trends in source apportioned-PM 2.5 mass concentrations were quantified using different approaches and then correlated with the available emission inventories (EIs) for the U.S. and Canada. Individual or aggregated EIs were considered. The changes in potential source areas of precursors of secondary sources were then evaluated using a new trajectory ensemble method: differential concentration-weighted trajectory (DCWT) analysis. These results provide information on which source emissions were substantially reduced over the past decade, identify those sources that did not change during this period, and identify those sources whose emissions increased.

Materials and methods
2.1. Overview of sites, sampling, analytical methods, and data consistency Details of the data repositories, data handling, blank and OC artifact corrections, PMF strategies, factor interpretation and source identification are extensively reported in Squizzato et al. (2018b) and in the Supplementary Information Sections S1-S4. Briefly, chemical speciation data for PM 2.5 during 2005-2016 were retrieved from the EPA Chemical Speciation Network (CSN; https://www.epa.gov/outdoor-airquality-data). Long-term data are available at 8 sites across NYS (Fig.  S1). Six sites are representative of city-wide air pollution in Albany (ALB), Buffalo (BUF), Rochester (ROC), and New York City (NYC with 3 sites: The Bronx (BRO), Manhattan (MAN), and Queens (QUE)). Two sites (Pinnacle State Park, PIN, and Whiteface base, WHI) are categorized as rural, i.e. in areas not directly influenced by large anthropogenic sources. Details of the sampling sites are provided in the supplementary materials (Section S1); maps are provided in Figs. S2-S9. Fig. S10 reports the monthly averages of air temperature recorded at weather stations close to the sampling sites.
Daily (24-h) samples were collected every third or sixth day and analyzed for species that would provide mass closure (Section S2): (i) elemental carbon (EC) and organic carbon (OC) determined by thermooptical analysis, (ii) major inorganic ions by ion chromatography, and (iii) elements by energy-dispersive X-ray fluorescence. Details of sampling methods, analytical protocols, and quality assurance/quality control are summarized in Solomon et al. (2014) and Section S2. Methods and analytical protocols for EC and OC changed between 2007 and 2009 (see Section S3): the sampling methods transitioned from MetOne SASS to URG 3000N samplers and the analytical protocol from NIOSH4050-like thermal optical transmission (TOT) to OC/EC IM-PROVE_A and thermal optical reflectance (TOR). In order to help the identification of gasoline and diesel traffic, single EC and OC fractions were used as input for the PMF analysis, i.e. EC1, EC2, EC3 for EC for the IMPROVE_A protocols and OC1, OC2, OC3, OC4 for both protocols, plus OP (pyrolized organic carbon) under the IMPROVE_A protocol. Blank correction strategies, including the handling of evaporative/absorbing/shrinking OC artifacts are summarized in Section S3.

PMF analysis and source interpretation
PM 2.5 chemical species define chemical fingerprints (profiles) that may be attributed to specific sources. PMF apportions the measured mass of an atmospheric pollutant at a given site to potential emission sources by solving a mass balance equation (Hopke, 2016). A quick overview of PMF is provided elsewhere (Paatero and Tapper, 1994;Paatero, 1997;Hopke, 2015Hopke, , 2016. U.S. EPA PMF version 5 was used. Beyond the determination of the trends in the source contributions presented in this study, the ultimate goal of the source apportionment analysis will be the detection of the impacts of PM 2.5 sources on hospital admissions for a variety of health outlines. These latter results will be presented elsewhere. Under this view, having several period-specific profiles would add additional variation into both the trend and exposure assessment analyses. In addition, without detailed organic tracers included in the PMF analysis, it would be difficult to be certain that the variations were due to changes in the source profiles relative to changes in the sampling and analytical methods (e.g., changes in sampling systems, OC/EC analytical methods and a change in laboratories that changed the protocols for the elemental analyses from filterbased systems to secondary target systems). With all those changes being made in the study period, the only practical way to obtain good source contribution estimates for the trends and exposure assessments was to analyze the data over long periods. Hence, the best compromise was to run PMF analyses over two periods, i.e. pre-and post-OC/EC change switch. Details are reported in Section S4. The PMF protocols adopted in this study followed the guidelines and best practices reported in the literature (e.g., Reff et al., 2007;Belis et al., 2014;Paatero et al., 2014;Brown et al., 2015;Masiol et al., 2017b). The detailed interpretations of the sources were provided by Squizzato et al. (2018b) and summarized in Section S4.

Trend analyses
Data analysis was performed in R 3.4.2 (R Core Team, 2018). Trends in PMF contributions were analyzed similarly to Masiol et al. (2018a), i.e. adopting three different approaches applied on the monthly-averaged data: • The shape of trends and their seasonal variations were investigated through the seasonal-trend decomposition time series procedure based on 'Loess' (STL, Cleveland et al., 1990). STL was previously used in several AQ studies (e.g., Carslaw, 2005;Bigi and Harrison, 2010;Masiol et al., 2014). Missing monthly-averaged data were linearly interpolated. STL was performed in the robust mode, as data were generally not normally distributed (Shapiro-Wilk test at p < 0.05). Confidence intervals (95%) for STL were then assessed by bootstrap (n = 5000).
• The linear trends over the entire period were assessed by the Mann-Kendall trend test (p < 0.05) and the Theil-Sen nonparametric estimator of slope (MK-TS, Mann, 1945;Theil, 1950;Sen, 1968;Kendall, 1975) implemented in the "openair" package (Carslaw and Ropkins, 2012). This technique assumes monotonic linear trends and is therefore useful to estimate the interannual trends. Missing data such as a break in the time series at the Bronx site were ignored.
• Since the MK-TS approach does not take into account the shape of trends (well depicted by STL), the presence of possible breakpoints splitting the trends in segments with different slopes was investigated with piecewise regression. The most plausible locations of breakpoints were estimated using the "segmented" package (Muggeo, 2003(Muggeo, , 2008, which uses an interactive procedure requiring approximate starting dates for the breakpoints and implementing bootstrap restarting (Wood, 2001) to estimate slopes and breakpoint positions with standard errors. The breakpoints were detected by visual analysis of the data and STL shapes. Piecewise regression was used only when these criteria were met: (i) the statistical significance of piecewise slopes was p < 0.05; (ii) the incremental F-test between the simple linear and the piecewise regressions was statistically significant at p < 0.05; (iii) the detected segments were computed over a relevant number of monthly data (at least 1 year); (iv) the errors in breakpoint locations were reasonably low (usually lower than 1 year, but a few exceptions were accepted); and (v) the resulting adjusted r 2 was over 0.3. When the piecewise regressions did not meet these criteria, the results of simple linear regressions were provided. STL analyses were only performed for the entire time series. However, MK-TS and piecewise regression were performed for the whole period and each season separately (winter: Dec-Jan-Feb; summer: Jun-Jul-Aug; transition: remaining 6 months) (Emami et al., 2018;Masiol et al., 2018a) to account for the seasonal patterns of some sources.

Differential concentration-weighted trajectory models
One-hour time-resolved back-trajectories were computed between 2005 and 2016 with the HYSPLIT version 4 model (Stein et al., 2015;Rolph et al., 2017) using NCEP/NCAR Reanalysis 1 × 1 deg. gridded meteorological data (Kalnay et al., 1996). HYSPLIT set-up: 120 h with a starting height of 500 m agl., and the vertical velocity model. Starting points were set at each city coordinates, except for NYC sites where a single starting point was set at a centroid location with respect to BRO, MAN and QUE. Source contributions were available for integrated 24-h samples; thus, they were matched with the 24 trajectories calculated for that day (Kim and Hopke, 2004a).
Trends of secondary sources may reflect changes in emissions over the source areas for the PM precursors. Long-range transported contributions are usually evaluated by trajectory ensemble methods based on back-trajectories (Lupu and Maenhaut, 2002;Hsu et al., 2003;Zhou et al., 2004;Pekney et al., 2006;Kabashnikov et al., 2011;Fleming et al., 2012;Brereton and Johnson, 2012;Squizzato and Masiol, 2015;Hopke, 2016). In this study, the concentration weighted trajectory (CWT) model (Hsu et al., 2003;Zhou et al., 2004) was used. Briefly, each grid cell ij in a grid domain was used to compute the weighted concentration obtained by averaging sample concentrations that have associated trajectories passing the grid cell according to: (1) where i and j are the coordinates of the grid cell, k the trajectory index, N the number of trajectories, C k the PMF source contribution measured at the site upon arrival of the trajectory k, and τ ijk is the residence time of trajectory k in the ij cell. Changes in CWT values over time can be then evaluated by using a differential CWT (DCWT) analysis defined as the difference between the CWT concentrations estimated in the two time different periods p 1 and p 2 , = − Δp p p 2 2 . A similar approach was presented in Masiol et al. (2018b), but a different trajectory ensemble method (probability source concentration function) was used. The grid resolution was 1°latitude x 1°longitude. The DCWT values were only computed for grid cells with at least 10 endpoints for both periods. The CWT and DCWT calculations were performed in R 3.4.2.

Relationship with emission inventories
The source contributions of major secondary sources were then used to perform correlation analyses with EIs available for North America. U.S. emission data for SO 2 and NO x were retrieved from various repositories (details provided in Section S3) including: (i) annual average US EPA air pollutant emissions trends (APETs, US EPA, 2018c) following Tier I sectors and aggregated by State; (ii) annual average national emission inventories (NEIs, released every third year: 2005(NEIs, released every third year: , 2008(NEIs, released every third year: , 2011(NEIs, released every third year: , 2014; US EPA, 2018d) following the Tier I sectors and their sub-categorization Tier II aggregated by State or County; (iii) clean air markets (CAM) program (US EPA, 2018e) supplying monthly emission data only for SO 2 and NO x from large stationary sources across the U.S., such as power plants and aggregated by State. Yearly-resolved Canadian EIs were provided by the Canada air pollutant emission inventory (APEI) historical trends (Environment Canada, 2018). These data were aggregated by province.
Correlation analyses were computed using different approaches based on PM 2.5 sources identified by PMF (Section S3): • Since they were of regional origin, source areas for secondary sulfate and nitrate were identified by DCWT analyses and from previous studies (Masiol et al., 2017a;Masiol et al., 2018b). Annuallyaveraged PMF contributions were then correlated with the sum of the Tier I APETs for the PM precursors (SO 2 and NO x for secondary sulfate and nitrate, respectively) among the U.S States identified by DCWT (Table S1).
• Since source areas for secondary aerosols also encompass Canadian regions, some selected NEI Tier II sectors were matched with their equivalent APEI sectors to estimate the combined U.S. and Canadian emissions for a limited set of aggregated sectors (Table S2).
• Secondary sources are also matched with CAM data. In this case, monthly averaged source contributions were matched with the CAM actual emissions data for the same set of U.S. States/Provinces selected for APETs (Table S1).

Results and discussion
The timeline of the main emissions reductions actions implemented at federal and/or national levels are summarized in Fig. S12. From 2005 to 2016, regional and local emissions have changed significantly. Vehicle emissions were significantly reduced (Maricq, 2007;Bishop and Stedman, 2008;Bishop et al., 2012;McDonald et al., 2013;May et al., 2014) with improved emissions controls and cleaner fuels. In mid-2007, all new heavy-duty diesel trucks were required to have particulate filters. NO x controls were added in January 2010. However, credit flexibilities allowed the sale of engines with NO x emissions greater than the 2010 limit (0.2 g/bhp-hr) through model year 2014. The light-duty vehicle (< 6000 lbs gross weight) Tier II emissions program (extended to vehicles up to 10,000 lbs) required an average sulfur in fuel standard of 30 ppm (from 120 ppm) with a sulfur cap of 80 ppm (reduced from 300 ppm), which was phased-in from 2004 to 2009. On-road diesel fuel sulfur content was reduced from < 500 ppm to ultralow sulfur diesel fuel (ULSD, < 15 ppm) in 2006 such that by October 1, 2006, 75% of the on-road diesel fuel was ULSD. ULSD fuels were subsequently required for nonroad vehicles by 2010 and for locomotives and marine vessels by 2014. Additionally, all distilled oil sold in New York State for any purpose (including building heating) were required to be ULSD after July 1, 2012. In New York City, there was substantial use of No. 6 (residual oil) or No. 4 (mixture of No. 2 and No. 6) oils for large building space heating. After 2011, the use of No. 6 oil was disallowed to ensure that by 2015, only No. 4 or cleaner oils could be burned. There still could be up to 1500 ppm S in No. 4 oil (Kheirbek et al., 2014) so that this source of SO 2 has not been completely eliminated within the period of 2005-2016. However, it is likely that there was a significant reduction in PM 2.5 from space heating systems by 2015 (Cleanheat, 2018).
Among the 8 sites, PMF identified and apportioned 12 sources. Seven sources were detected at all the sites: secondary sulfate (SS), secondary nitrate (SN), spark-ignition vehicles (GAS), diesel (DIE), road dust (RD), biomass burning (BB), and OP-rich aerosol (OP). The OP factor was detected after the changes in the EC/OC sampling and analysis protocols. The other 5 sources were road salt (RS, only detected in upstate cities, ALB, BUF, ROC), residual oil (RO, only detected at NYC sites), fresh sea salt aerosol (FSS, only detected at NYC sites), aged sea salt (AGSS, detected over NYC sites and rural sites), and industrial (IND, only detected in Buffalo). Site-specific profiles and temporal variability of each source were extensively discussed in Squizzato et al. (2018b).
These trend analyses generated many plots and tables. Thus  Table S3 and summarized in Fig. S25. Piecewise slopes and breakpoints are reported in Table S4 along with the adjusted r 2 values, root mean squared errors, and mean absolute errors expressed in μg/m 3 . Analogous tables calculated for each season are provided for winter (Tables S5 and S6), transition (Tables S7 and S8) and summer (Tables S9 and S10). Additional figures are also provided by sources for each season in the supplemental information (SI) file. The summary of the seasonal MK-TS results is shown in Fig. 2.

Secondary sulfate (SS)
The median linear MK-TS slopes computed over the whole study period (Fig. S25) show statistically significant (p < 0.001) downward trends for SS at all sites, ranging from −0.19 μg/m 3 /y (ROC) to −0.36 μg/m 3 /y (BRO, QUE) representing reductions of −6 to −9%/y. The shapes of the trends were relatively constant through the years, as evidenced by STL and the absence of breakpoints ( Fig. 1 and Fig. S13; Table S4). The SS profiles include organic carbon in addition to ammonium and sulfate. There are known associations between secondary organic aerosol (SOA) and ammonium sulfate (e.g., Kleindienst et al., 1999;Jang et al., 2003;Murphy et al., 2006;Zhang et al., 2007;Hallquist et al., 2009;Seinfeld and Pandis, 2016) that appear as the presence of organic carbon in the SS source profiles . However, there was a notable change in the distributions of the SS contributions between the pre-EC/OC change period and following the OC protocol change with wider distributions before the change. There was more OC associated with the SS profiles prior to the protocol changes, compared to after with the difference represented by the new OP factor that was resolved only after the protocol changes.
Seasonally, the sharpest decreases in the MK-TS slopes were observed during summer (Fig. 2, Figs. S26-S28). This result reflects the higher sulfate concentrations expected in warmer seasons due to the increased concentrations of atmospheric oxidants (mainly hydroxyl radical) and the consequent higher conversion rate of SO 2 to sulfate. For example, Masiol et al. (2017a) reported that the sulfate/SO 2 ratios in NYC during 2009/2010 ranged from < 0.3 in winter to 0.8-1 during warm periods.
The CWT analyses calculated for the two rural sites (Fig. S29) show that the upper Ohio River Valley and the Tennessee Valley were major source areas of SS similar to previous reports in NYC (Masiol et al., 2017b) and Rochester (Emami et al., 2018;Masiol et al., 2018a,b). The contributions clearly declined over time. The DCWT calculated between 2005/2007 and 2014/2016 (Fig. 3) shows a substantial decrease over the eastern United States. The sequential DCWTs between consecutive periods (Fig. S30) also show small areas with apparent increases in SO 2 emissions between 2005/2007 and 2018/2011 in Florida, Texas and northcentral States where oil fracking operations were increasing in intensity. The differences detected in the seasonal MK-TS (Fig. 2) are more evident for the more southerly sites (NYC and PIN) comparted to BUF, ROC, ALB, and WHI. These sites may reflect their relative closeness to SO 2 sources such as coal-fired power plants in Pennsylvania. Figure S31 shows the coefficients of determination (r 2 ) computed for the correlations between the SS contributions at each site and the various emissions calculated from the EIs across the potential source area identified by CWT results. This region included 31 eastern U.S. states and 3 Canadian provinces (Table S1). Most of the APETs and aggregated NEIs/APEIs were highly correlated (r 2 > 0.6) with almost all the economic sectors at all sites. Alternatively, the correlations with the CAM data were high only in summer at the NYC sites and PIN. The APETs and NEIs correlations were computed over a limited number of points (12 years). However, these results suggest that the drop in SS can be attributed to the reductions in emissions that resulted from a combination of imposed controls and economic drivers including the 2008 recession and the change in the relative costs of coal and natural gas that led to substantial fuel switching (Squizzato et al., 2018a).

Secondary nitrate (SN)
Negative annual trends were found for SN ( Fig. S25; Table S3), ranging from −0.02 μg/m 3 /y (rural sites) to approx. −0.2 μg/m 3 /y (BRO, MAN). On a percentage basis, the decrements ranged from −5 to −11%/y (−3%/y at BUF, but these trends were not statistically significant: p > 0.05). The annual STL trends were generally linear over all the study period. However, one breakpoint was detected between 2008 and 2009 in NYC (BRO and MAN).
Seasonally, the downward trends were steeper in winter (range 0 to −0.33 μg/m 3 /y, Fig. 2, Figs. S32 to 34). This result is strongly affected by the seasonality of particulate nitrate given high ammonium nitrate volatilization (e.g., Hering et al., 1988;Hering and Cass, 1999) at higher ambient temperatures. The CWT analyses revealed potential source areas north central states and the central Canadian provinces (Fig. S35) (Table S1). Most of the APETs and aggregated NEIs/APEIs were highly (r 2 > 0.6) or moderately (0.4 < r 2 < 0.6) correlated with almost all economic sectors for all sites, except Buffalo (r ≤ 0.16). The correlations with CAM data were higher in winter due to the seasonality of SN concentrations.
Analogous to the SS trends, the decline in SN can be attributed to the previously discussed mitigation measures and economic changes. However, the correlations with EIs were lower than for SS, reflecting the nature of SN that is more affected by local sources than SS (Zhao et al., 2007).
In Buffalo, the absence of a clear downward trend and no correlation with any of the EIs suggest that source changes had less effect on the concentrations than at the other sites across NYS. This result could have been due to: (i) fewer samples (1 sample every sixth day), and (ii) the presence of a local industrial source accounting for a substantial fraction of the nitrate . The correlation of the combined SN + industrial source contributions were still not well correlated with the estimated emissions. The absence of trends for SN in BUF remains unclear and will require more detailed study.

Traffic-related sources (GAS, DIE, RD)
GAS was the only source that showed upward annual trends at all urban sites on both annual and seasonal bases (Figs. 1, 2, Figs. S38 to 40; Tables S3, S5, S7, and S9). Annual urban trends ranged from 0.02 μg/m 3 /y (ROC, not statistically significant) to ∼0.2 μg/m 3 /y (ALB, BRO, MAN), while essentially null trends were estimated at the rural sites. The overall increase in spark-ignition contributions is consistent with the increase in registered vehicles (Fig. S41) in the past decade. Although there was a further decrease in fuel sulfur content over the period (Fig. S12), spark-ignition emissions still represent a major PM 2.5 source in the urban areas. Slightly higher slopes were observed during summer at almost all sites (Fig. S40). This result is consistent with the observed increase in secondary organic aerosol . Spark-ignition vehicle emissions represent a potential source of secondary organic aerosol precursors in the urban areas (Bahreini et al., 2012;Gordon et al., 2014Gordon et al., , 2013Hayes et al., 2013), particularly with the recent focus on intermediate volatility organic compounds (IVOCs) (Zhao et al., 2014(Zhao et al., , 2105(Zhao et al., , 2016. Annual MK-TS trends for DIE were all statistically significant (Table  S3) and similar patterns were observed in all seasons (Figs. S42-S44). Generally, DIE showed decreasing trends, except in ROC and at WHI.  Tables SI1-3-5. M. Masiol et al. Atmospheric Environment 201 (2019) 110-120 Downward slopes varied between −0.03 μg/m 3 /y in MAN and −0.13 μg/m 3 /y in BUF. The decreasing diesel contribution may have results from the implementation of the Clean Heavy-Duty Bus and Truck Rule and the gradual replacement of old trucks with new, cleaner vehicles as well as retrofit programs such as on public buses in NYS. Cleaner non-road fuels were also phased in for railroad diesel engines and ships docking in Albany and NYC. The upward trend in WHI should be viewed with caution as it is likely related to the concurrent effects of changes in sampling collection and analytical protocols for EC/OC and the near-detection limit concentrations for EC (Rattigan et al., 2016). However, the piecewise regression analysis in ROC and WHI (Table S4) show that the slopes can be split in two segments with increasing concentrations (0.12 and 0.05 μg/m 3 /y, respectively) through 2008 or 2009 followed by decreasing slopes from mid-2012 (−0.13 and −0.02 μg/m 3 /y, respectively). This period represents a time when the 2008 regression reduced the need for transport, so an additional incremental reduction may have occurred because of economic forces. RD likely represents a combination of different sources, i.e. local soil resuspension, road surface abrasion, and local construction activities . Although RD was detected at all sites, the contributions may differ site-to-site resulting in different trends across NYS on both on annual and seasonal bases (Fig. 2, Figs. S45 to S47). Significant downward annual trends were observed at ROC, MAN, QUE, and PIN (Table S1), with similar slopes detected throughout the seasons (Tables S5, S7, and S9). Different seasonal patterns of BB were observed among the sites. The deseasonalized STL trends, MK-TS slopes, and piecewise regressions performed by season are shown in Figs. S48-S50. The three upstate sites (ALB, BUF, and ROC) showed the highest concentrations during winter and a slight increase during summer. The rural and NYC sites showed a clear increase during summer. Recreational activities (e.g., barbeques) or long-range transport from wildfires were the likely summer sources of BB across the state, while BB for domestic heating was the major winter source for the upstate sites. BB was estimated to contribute up to 30% of average PM 2.5 mass concentration in winter in Rochester (Wang et al., 2012). In addition, Monroe and Eire County, where ROC and BUF are located, respectively, were the main consumer of wood for burning in NYS (91151 and 98515 wood mass burned tons, respectively) (NYSERDA, 2016). Under this view, the decreasing trends observed during winters at BUF and ROC (Table S5, Fig. S48) can be ascribed to a decrease in the use of wood burning for domestic heating. Similarly, the decreasing trend observed at BRO and QUE (Table S5, Fig. S48) can be associated with the variation in wood burning for heating. The sampling sites are located close to residential areas, and they can be affected by residential emissions of the New Jersey shore and the Suffolk County (79700 wood mass burned tons in 2011) (NYSERDA, 2016).
Conversely, PIN showed a significant positive MK-TS slope (0.02 μg/ m 3 /y). The largest significant slope was observed at PIN during transition period (Table S7, Fig. S49). This result was likely the result of the increased use of wood to supplement domestic heating in areas where other fuels (natural gas, propane) were difficult to retrieve (many rural areas do not have access to gas) or when the fuel prices were high (Loughlin and Dodder, 2014;NYSERDA, 2016).

OP-rich (OP)
OP was only observed when the IMPROVE protocol was used for the OC/EC analyses, therefore trend analyses start between 2007 and 2009. Overall, OP contributions showed a significant decrease at the urban sites ranging between −0.07 μg/m 3 /y at QUE and −0.26 μg/m 3 /y at MAN. The rural sites did not exhibit statistically significant trends. Seasonally, no differences were found for MK-TS trends, making the decrements constant throughout the year (Tables S5, S7, and S9; Figs. 51-53).
The OP factor represented a fraction of secondary aerosol. A similar factor was resolved in a number of prior source apportionment studies using IMPROVE carbon fractions (Kim and Hopke, 2004b;Kim et al., 2004). The source profile also contains potassium (a marker of biomass combustions) and sulfate (due to the heterogeneous acidic catalyzed reaction between the sulfuric acid and VOCs; Jang et al., 2003). Prior studies have suggested that it included emissions from large wild fires (Kim and Hopke, 2004c). However, the CWT maps (Fig. S54) show that the higher contributions occurred when air masses arrived from a wide area across the continental U.S. and do not reveal any specific source area. The DCWT (Fig. S55) shows a strong increase in concentrations between 2009/2010 and 2011/2013, but a decrease in the most recent period. The causes of this behavior are unclear since no statistically significant differences in air temperature were measured among the three periods that would have substantially changed the biogenic VOC emission rates. However, 2011 was characterized by having the largest area of burned woodlands in the last several decades (1,791,469) (https://www.ncdc.noaa.gov/societal-impacts/wildfires/). Correlation analysis was not possible given the short period when OP was detected and the resulting small number of available data points.
3.6. Road salt, fresh and aged sea salt (RS, FSS, AGSS) The salt-related sources have similar fingerprints, but FSS accounts for more chlorine than any other source. AGSS accounts for a substantial portion of the sodium, but little chlorine since chlorine was almost totally depleted by reactions with gaseous acids. RS had a strong seasonal pattern, appearing primarily in the winter since it was used on icy roads. RS was only detected in the upstate sites (ALB, BUF, ROC), where cold and snowy winters require de-icing salts on the roads. It did not appear at the rural sites since they are not close enough to a road network. FSS was only found in NYC sites, and AGSS was only resolved at the NYC and rural sites. Although the trend analyses revealed some significant trends on both annual and seasonal bases (Tables S3, S5, S7, S9, S46 to S54), the quantitative changes were small for all three sources. A trend analysis for the combined contributions of the three salt-related sources did not show any observable trend. Thus, total saltrelated contributions were essentially constant across NYS.

Residual oil (RO)
RO was only detected in NYC sites and was attributed to the use of No.6/4 residual oil for domestic/commercial heating in large building heating systems. However, other sources with similar fingerprint are present in the area, including a large oil refinery just south of Elizabeth, NJ (∼35 km WSW of QUE), ship emissions in the NY-NJ harbors, and the container port at the Port of Elizabeth, NJ. On an annual basis, the contributions showed a downward trend across the 3 NYC sites ranging from −0.01 μg/m 3 /y (QUE) to −0.05 μg/m 3 /y (BRO). Slightly higher slopes were observed during the winter and transition months (Tables S5, S7, and S9; Figs. S65-S67).

Industrial source in BUF
An "industrial" source was detected only in BUF. Considering the number of former and current industrial facilities in the Buffalo area, this factor represented a mixed source of industrial emissions (e.g., steel processing plant and coke production) and resuspension from various locations including the former secondary lead smelters located close to the monitoring site. The trend analysis did not identify any significant annual trends. However, the piecewise regression detected three breakpoints. Industrial contributions increased until June 2006 (0.23 μg/m 3 /y) followed by sharp declines at the end of 2007 when the 2007-2008 recession began (−0.39 μg/m 3 /y). A slight but significant increase was observed in the post-period (0.07 μg/m 3 /y) (Table S4), but the IND contributions were still lower relative to the pre-recession period. Similar patterns were observed during the transition and summer months but not in winter (Figs. S68-S70).

Conclusions
The United States experienced significant emissions reductions during the last two decades due the implementation of multiple mitigation strategies designed to control emissions from power plants, vehicles, and space heating. In addition, economic drivers such as a major recession and the changing relative costs of fuels for electricity generation also affected pollutant concentrations. Negative trends were detected across NYS for secondary sulfate (from −0.19 μg/m 3 /y at ROC to −0.36 μg/m 3 /y at BRO and QUE) and secondary nitrate (from −0.02 μg/m 3 /y at the rural sites to approx. −0.2 μg/m 3 /y at BRO and MAN) reflecting the effectiveness of the policies aimed to control gaseous precursors emissions (i.e. control of power plant emissions, NO x control systems for vehicles, use of ultralow [ < 15 ppm] sulfur fuels). However, spark-ignition vehicles, a major PM 2.5 source, experienced upward annual trends at all the urban sites with slopes ranging from 0.02 μg/m 3 /y (ROC, not statistically significant) to ∼0.2 μg/m 3 /y (ALB, BRO, MAN). The increase was consistent with the increase of the secondary organic aerosol and of the number of registered vehicles in NYS. Other sources exhibited different trends among the various sites.
The decreasing trends found across NYS were particularly important for the secondary sources that were the real target of the past mitigation policies adopted in the U.S. However, relating the trends to specific mitigation policies was complicated by (i) the synergy of multiple mitigation measures at different time-scales and involving different economic sectors, (ii) the phase-in of multiple programs and timings for roll out of plans or turnover in production processes, and (iii) the economic factors including the 2008 recession and the historically low price of fracked natural gas relative to the costs of other fossil fuels.
Although it is not possible to disaggregated the causal factors, this study shows that trends in single PM 2.5 sources are generally consistent across NYS. Policy-makers must identify those steps needed to further improve AQ and protect public health. Although it is not possible to fully separate the relative effects of policy and economics, the trend analysis of this long time series of PM 2.5 apportioned sources can be a useful tool to examine changes in source contributions and help direct the development of future AQ improvement strategies.