Source attribution of insoluble light-absorbing particles in seasonal snow across northern China

Seasonal snow samples obtained at 46 sites in 6 provinces of China in January and February 2010 were analyzed for a suite of chemical species and these data are combined with previously determined concentrations of insoluble light-absorbing particles (ILAP), including all particles that absorb light in the 650–700 nm wavelength interval. The ILAP, together with 14 other analytes, are used as input to a positive matrix factorization (PMF) receptor model to explore the sources of ILAP in the snow. The PMF analysis for ILAP sources is augmented with backward trajectory cluster analysis and the geographic locations of major source areas for the three source types. The two analyses are consistent and indicate that three factors/sources were responsible for the measured light absorption of snow: a soil dust source, an industrial pollution source, and a biomass and / or biofuel burning source. Soil dust was the main source of the ILAP, accounting for∼ 53 % of ILAP on average.


Introduction
Snow is the most reflective natural surface on earth and thus a major component of the global albedo.Small amounts of insoluble light-absorbing particles (ILAP) deposited in snow can reduce the snow albedo significantly (Warren and Wiscombe, 1980).One such ILAP, perhaps the single most important, is black carbon (BC).Black carbon particles result from incomplete combustion of organic materials such as fossil fuels or biomass (Bond et al., 2004;Bond and Bergstrom, 2006).BC aerosols have been assessed as a significant anthropogenic forcing for twentieth century climate change and as one of the largest sources of uncertainty in various assessments of radiative forcing for climate change (e.g., Jacobson, 2001;Hansen et al., 2005;IPCC, 2007;Bond et al., 2013).While much of the effort to assess the impact of ILAP on snow albedo has centered on the Arctic snow fields (e.g., Clarke and Noone, 1985;Grenfell et al., 2002;Forsström et al., 2009;Perovich et al., 2009;Doherty et al., 2010), it has been suggested that the effect of ILAP on snow may be more important at midlatitudes because snow at these lower latitudes is exposed to more sunlight and may be closer to ILAP sources (Qian et al., 2009;Ye et al., 2012).While currently under-studied, radiative forcing due to ILAP in snow in middle latitudes may also impact the timing and rate of snow melt (e.g., Qian et al., 2009;2011).
Given the importance of the impact of ILAP on snow albedo, it is urgent to understand the sources of the ILAP in the snow.This is by no means a straightforward issue.While most studies have focused only, or primarily, on snow albedo reduction due to BC, recent studies in the Arctic have suggested that in a substantial portion of the ILAP in snow is not BC (Doherty et al., 2010).Further, there are in fact multiple sources for the ILAP, with biomass burning as the predominant source in the Arctic (e.g., Hegg et al., 2010).
For these reasons, an assessment of the sources of ILAP in snow over a broad midlatitude venue would be timely and useful.A field campaign recently conducted across northern China in January and February of 2010 provides an excellent opportunity for such an assessment (Huang et al., 2011).This campaign is particularly valuable because of the potential diversity of the sources of ILAP in Chinese snow, including rapidly increasing industrial emissions.the composition and potential sources of airborne aerosols in China suggest both temporal and spatial variations in the sources of snow ILAP (e.g., Cao et al., 2006;Zhang et al., 2008;Yang et al., 2009), but quantification of the contributions of various sources to the ILAP is as yet sparse and source attribution of ILAP in snow has not yet appeared in the literature -to our knowledge.In this paper, we present source attribution of ILAP in seasonal snow across northern China.The analysis presented here compliments the earlier study of Wang et al. (2013).In this earlier study, measurements are presented of ILAP in snow and an assessment made of the contribution of three light absorbing species to the total absorption by ILAP.In this study, we determine the sources of the ILAP found in the snow (e.g., industrial emissions) using additional analyses of ILAP chemical composition coupled with positive matrix factorization (PMF) modeling.

Snow sample collection
Seasonal snow was collected in January and early February 2010 on a road trip encompassing 40 sites in the provinces of Inner Mongolia, Heilongjiang, Jilin, and Liaoning.A shorter trip in mid-February 2010 obtained snow samples at six sites in Qinghai and Gansu Provinces.The sampling locations are shown in Fig. 1, and are numbered following the temporal order of the sampling.Also shown in Fig. 1 are the locations of desert regions within China that are prime potential sources of soil dust (e.g., Zhao et al., 2003;Zhang et al., 2003).Areas of heavy industrial emissions, based on the study of Zhang et al. (2009), as well as open fires detected by the MODIS Fire Mapper during the deposition period for the snow sampled (see discussion in Sect.2.4), are also shown.In order to discuss and contrast different regions, we divided the 46 sites into 7 groups by geographical region.The sampling sites were chosen to be far from local sources of pollution, so that the data could represent large areas.An exception is site 4, which is only 30 km from the large city of Baotou and was selected to investigate the influence of the city.In the semiarid regions, the snow was thin and patchy, so that sampling was sometimes possible only in drifts that were much deeper than the average snow depth.At each site, snow pits were dug and snow samples were gathered from two vertical profiles that were noted by "left" and "right" to partially encompass potential small-scale variability and to assess the precision of our measurements.For the analysis presented here, only samples from the top 6 cm of the snow pack are utilized (see Sect. 2.4).The snow was stored in plastic bags and kept frozen during the trip until the time of filtration in the temporary laboratories used en route.On occasion, some melting occurred during transport, necessitating correction for losses to the plastic bags (Wang et al., 2013).The filtering technique used was that pioneered by Clarke and Noone (1985) and also employed by Doherty et al. (2010), Huang et al. (2011), andYe et al. (2012).Each sample was transferred into a clean glass beaker, melted rapidly in a microwave oven and then immediately filtered through a 0.4 µm Nuclepore polycarbonate filter to separate out the black carbon and other water-insoluble species from the soluble analytes.Aliquots of 40-50 ml meltwater, taken both before and after filtration, were treated with biocide to prevent bacterial-induced decay of potential organic analytes, saved in small bottles, and refrozen for later chemical analysis for the source attribution.The filters were also returned to the laboratory for measurements with the Integrating Sphere/Integrating Sandwich (ISSW) spectrophotometer as reported by Wang et al. (2013), and for chemical composition analyses (Hegg et al., 2009(Hegg et al., , 2010)).

Chemical speciation
The filter samples were analyzed for BC and other lightabsorbing aerosols using the ISSW spectrophotometer (Grenfell et al., 2011, andDoherty et al., 2010) previously employed for analysis of Arctic snow samples.The values used here are taken from Wang et al. (2013).In this analysis, we only use the parameter C max BC (ng g −1 ), which is defined as the mass of BC per mass of snow assuming that all of the aerosol light absorption in the spectral interval 650-700 nm is due to BC (Grenfell et al., 2011), i.e., it is the mass of the BC calibration aerosol (with a mass absorption cross section of 6.3 m 2 g −1 at 550 nm and 5.0 m 2 g −1 at 650 nm -the value at 650 nm was used in the conversion) which would produce the equivalent of the observed light absorption.We used this parameter in preference to others reported in Wang et al. (2013) because it is the most closely coupled to the absorption measurements and involves the fewest assumptions.It can be considered as a scaled, total absorption by all ILAP since it varies linearly with the aerosol light absorption on the filters.The partitioning between various ILAP sources is left to the PMF analysis.Various correction factors, for example for filter undercatch and losses to the walls of sample containers, were determined and applied to all ILAP (and other species) concentrations (Wang et al., 2013).
The laboratory analysis of the water samples generally followed the procedures described by Hegg et al. (2009Hegg et al. ( , 2010)).In the laboratory, the snow aliquots, both filtered and unfiltered, were melted and analyzed for standard anions via ion chromatography with conductivity detection (IC-COND), for carbohydrates via both liquid chromatography-mass spectroscopy (LC-MS) and ion chromatography with pulsed amperometric detection (IC-PAD), and for various elements via inductively coupled plasma-optical emission spectroscopy (ICP-OES).Details on these procedures are given in Gao et al. (2003).The unfiltered snow samples corresponding to the filtered samples were also analyzed to assess the possible impact of filtering on the soluble analytes.No significant difference between filtered and unfiltered samples was observed for the soluble analytes (e.g., sulfate, nitrate, oxalate, etc.), but for some insoluble analytes normally associated with soil dust (e.g., Ca, Mg), differences were observed.Hence, values for these analytes were derived from unfiltered snow aliquots acidified to 0.1 N HNO 3 and leached for 48 h prior to analysis by ICP-OES. .
Using the techniques described above, we obtained the concentrations (ppbm) of two carbohydrates (levoglucosan and xylitol), four standard anions (chloride, nitrate, sulfate and oxalate) and six elements (Ba, Ca, K, Mg, Si and Cu), for which the sampling statistics are shown in Fig. 2. With the exception of Cu, the statistics suggest no anomalous samples.For Cu, on the other hand, there are clearly outliers skewing the distribution function.Examination of the individual samples reveals that the skewness is associated with a single sample, from site 14, for which an exceptionally high Cu concentration was obtained.Several other trace metal concentrations (e.g., Fe, Zn) were also quite high at this site but the ILAP value was similar to that from adjacent sites.We have included the sample in the analysis since there is no compelling reason to reject it.Test runs of the receptor model with and without the sample show that it has no significant impact on the results.

Source attribution using positive matrix factorization
Receptor models, developed for source characterization, provide source attribution by identifying and quantifying source profiles and contributions based on the internal variance structure of a data set.Like most receptor models, the PMF receptor model is a multivariate factor analysis tool that decomposes a matrix of speciated sample data into factor contributions and factor profiles.The factor profiles then need to be interpreted by an analyst as to what source types they represent using available ancillary information such as wind direction analysis, emission inventories and measured source composition profiles.The PMF method solves the factor analysis problem using a least squares analysis with data point weighting and non-negativity constraints (Paatero andTapper, 1993, 1994;Paatero, 1997).The particular version of PMF employed here is the US EPA PMF 4.1 model.
While not yet widely utilized, the model differs from previous EPA versions (e.g., 3.0) mainly in its expanded treatment of error analysis which facilitates the choice of the optimum PMF solution (cf.Norris et al., 2010), and has been successfully employed in receptor modeling of ambient aerosol sources (Wang et al., 2012).Two data sets (i.e. a concentration/physical property data set and an associated uncertainty data set) are required as input for the PMF model.The equation-based approach to construction of the uncertainty files was used (Norris et al., 2008).Uncertainty estimates and detection limits for each of the snow water species were based on analysis of replicate standards with uncertainties calculated as twice the standard error of the mean for each analyte species.Uncertainties for ILAP were estimated as per Doherty et al. (2010).Normally PMF is applied to a time series at an observation site, i.e., temporal variance is analyzed.
In this paper, we treat northern China as a receptor site whose various geographic subregions are differentially impacted by the potential sources, i.e., we analyze spatial variance.Although atypical, such an application of receptor modeling is perfectly consistent with the mathematical structure of the technique and has in fact been successfully utilized by others (e.g., Paatero et al., 2003;Chen et al., 2007;Hegg et al., 2010).

Back trajectory cluster analysis and MODIS fire locations map
An important question is whether the observed aerosols arise entirely from local sources or if there are contributions from distant sources through long-range transport.Back trajectory modeling has been widely used as a standard tool for determining the source regions and transport patterns of air parcels observed at receptor sites, and has been used specifically in northern China (e.g., Huang et al., 2007Huang et al., , 2008;;Wang et al., 2010).In this study, 10-day back trajectories arriving at 500 m above ground level, initialized every 6 h (00:00, 06:00, 12:00, 18:00 UTC) for a representative site of each of the seven regions studied, were generated using the NOAA HYbrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model, version 4 (Draxler and Rolph, 2012; http://ready.arl.noaa.gov/HYSPLIT.php),and then used for further analysis.Each back trajectory contains endpoints describing the hourly location of an air mass in latitude and longitude coordinates.In order to compact the otherwise unwieldy amount of information in the individual back trajectories, and to facilitate its use for delineating source regions, cluster analysis has been applied to the back trajectory data.
Cluster analysis is a multivariate statistical technique in which a set of trajectories is subdivided into groups such that the trajectories in the same cluster are more similar than in other clusters (Kalkstein et al., 1987).It thus allows differentiation of the main routes by which airborne pollutants reach particular receptor sites.In this paper, cluster analysis is used to divide trajectories for each receptor region into seven cluster groups.Gridded meteorological data combined with the geographic information system (GIS)-based software TrajStat (Wang et al., 2009) is used and is based on Ward's hierarchical clustering algorithm as described and applied in Ward (1963) and Fu et al. (2010).
For the potential fire source distribution, we utilize data from the MODIS radiometer.MODIS is carried on both Terra and Aqua satellites and MODIS fire observations are made four times a day.The MODIS daily 1 km level 3 hotspot/fire products over East Asia and Asian Russia used in this study were from the Fire Information for Resource Management System (FIRMS; Davies et al., 2009), which was developed by the University of Maryland and is available online: http://earthdata.nasa.gov/firms.Industrial emissions were mapped from data available at http://www.cgrer.uiowa.edu/EMISSIONSDATA new/.The desert area information is available at the TrajStat website (Wang et al., 2009) as is the TrajStat software to plot all the geographic data (http://www.meteothinker.com/TrajStatProducts.aspx).
An issue common to both the back trajectory analysis and utilization of the FIRMS data is the time period over which to assimilate the data.To estimate this, the sampling depths of each snow sample utilized are first noted (depths varied but were always within the top 6 cm of the snow pack for this study).Then, daily precipitation and temperature (to aid in assessing the precipitation type) data from the China Meteorological Service were used to calculate at what calendar date the snow corresponding to the sample depth for each site began to accumulate (data available at: http://cdc.cma.gov.cn).Due to the non-urban nature of most of the sites, the reporting station closest to the site had to be used in lieu of data actually at the site.The time period over which to accumulate the FIRMS data and run the back trajectories was then the interval from 10 days prior to the starting accumulation date (to take into account the trajectory period) to the sampling date.However, a cautionary note is in order here.The snow in Regions 2 and 3 was commonly sparse and non-uniform (Wang et al., 2013).Because of this, sampling was sometimes done from snow drifts and the chronology of the deposition is unreliable.Nevertheless, the surface nature of the samples coupled with the meteorological station data suggest that the snow was in fact deposited approximately over the estimated time intervals.These intervals ranged from 2.5 months for Region 1 to 4.3 months for Region 7.

PMF model results
The PMF 4.1 model was employed with the Chinese snow data set in the model's robust mode (statistical outliers downweighted) with factor numbers ranging from two to six and always with six or more random seeds.A standard evaluation of the Q robust parameter as a function of the number of factors suggested that a three-factor solution is optimal, i.e., Q robust was a minimum for the three factors and closest to Q theoretical (Reff et al., 2007).The diagnostic regression for ILAP with the three-factor solution had a quite satisfactory R 2 value of 0.77.Rotational ambiguity was evaluated by varying the FPEAK parameter (Norris et al., 2008), and the three-factor solution was also found to be rotationally stable.
Factor, or source, profiles for the three factors are shown in Fig. 3. Interpretation of the profiles in terms of possible sources is never entirely unambiguous but seems relatively straightforward here.The first factor contains virtually all of the levoglucosan and non-soil K, both well-known markers for biomass and biofuel aerosol emissions (e.g., Simoneit, 2002;Hegg et al., 2009;Kehrwald et al., 2012).Indeed, levoglucosan has been found in Asian ice cores in regions of frequent biomass burning (Kawamura et al., 2012).We therefore consider it associated with a biomass and / or biofuel burning source.The second factor includes ∼ 70 % of the sulfate (essentially identical to non-sea-salt sulfate at these locations).Sulfate is, of course, a well-known industrial emission, particularly in China with its relatively large number of coal-fired power plants.These plants supply ∼ 63 % of China's power requirements (Cao et al., 2006), and account for most of China's SO 2 emissions.Industrial emissions are in fact the dominant source of atmospheric sulfate in northern China (Wang et al., 2011;Zhao et al., 2012).The only other species with over 50 % of its mass loaded onto the second factor is oxalate, normally considered a biomass burning emission (cf.Hegg et al., 2010).However, it does have other sources such as motor vehicle emissions (e.g., Kawamura et al., 1987) and is a major secondary organic aerosol constituent, as sulfate is a major secondary inorganic aerosol constituent (e.g., Kawamura et al., 1999).Indeed, it is now widely accepted that oxalate as well as sulfate has a major source in clouds and precipitation (cf.Sorooshian et al., 2007;Kaku et al., 2006).We speculate that the relatively high loading of oxalate on the sulfate factor is due in part to a similarity in the atmospheric pathways to deposition as well as some common sources.Based largely on the sulfate loading, but consistent with the oxalate, we interpret the second factor as an industrial/urban aerosol source.The third factor encompasses the major mass fractions of a number of well-known soil dust constituents including Si, Mg, K and Ca.Furthermore, the ratios of Ca/K and K/Mg, 5 and 0.63, respectively, are reasonably close to those reported for dust aerosol from the northern and western regions of China, namely 4 and 1 (Zhang et al., 2003).Hence, association of factor 3 with a soil dust source is quite plausible.The mean source attribution of the ILAP to the three sources just discussed is given, for each of the seven regions specified in Fig. 1, in Fig. 4.Not surprisingly, for regions 7, 1 and 2, situated either in or adjacent to major dust source ar- eas, the dominant source of ILAP is soil dust.Regions 3 and 4 have a somewhat less marked dust presence though dust is still the main ILAP source.For regions 5 and 6, situated more distantly from dust sources and relatively close to both major industrial regions and in areas where fires were prevalent during the snow deposition period (see Fig. 1), the dust source no longer dominates though it is still always important.For the data set as a whole, i.e., for the entire northern China region, soil dust contributed 53 % of the ILAP, industrial emissions 27 %, and biomass and/ or biofuel burning 20 %.This finding is somewhat different from that implicit in the results presented by Wang et al. (2013), particularly for the northeastern sampling sites where BC was estimated to contribute ∼ 80 % of the ILAP absorption.While our actual sources are not directly comparable to the species partitioning of ILAP absorption, our results would imply as much as ∼ 60 % of the light absorption would be due to BC in much of the northeast.Given the substantial uncertainty of the Wang et al. (2013) partitioning, this discrepancy is mild.

Atmos
Our results are substantially different from previous data sets we have examined with the same methodology, encompassing northern Siberia, Greenland, northern Canada and the North Pole region (Hegg et al., 2009(Hegg et al., , 2010)).This difference arises not only due to different geographic distributions of source regions (and their strengths) compared to the respective receptor sites, but also, of course, due to differing transport trajectories for the source emissions.To explore this latter issue, back trajectory analysis is employed.

Back trajectory cluster analysis
The dominance of the soil dust source for ILAP illustrated in of sample collection (Huang et al., 2011) and more objective estimates of the non-BC fraction of the ILAP based on the absorption Ångstrom exponent analysis, as described in Wang et al. (2013), suggest a strong non-BC source of ILAP in the snowpack for the northern China region as a whole.On the other hand, the finding that soil dust is a dominant constituent of the snow ILAP is quite different from our findings in other venues using the same methodology (Hegg et al., 2009;2010).Hence, we employ back trajectory cluster analysis as discussed in Sect.2.4 to independently test the plausibility of this finding.
An example of the results of the cluster analysis, in this instance for Region 7, is given in Fig. 5a.The trajectories shown are the mean trajectories of each of the clusters that encompassed 10 % or more of the total trajectories.Together, they account for ∼ 90 % of the total trajectories run for this region (536 total trajectories for Region 7, 336 for Region 5 and 400 for Region 4).Virtually all of the cluster trajectories have encountered arid source regions for soil dust within a few days travel time of the receptor region.Conversely, very few fire regions and essentially no built-up industrial areas are encountered.Hence, the trajectories are quite consistent with the dominance of the dust source as shown in Fig. 4. Another example of the cluster analysis, this time for the more complex source mix of Region 5, is shown in Fig. 5b.In this case, the back trajectories encounter fire areas and industrial emissions areas as well as dust source regions very near the receptor sites, in accord with the roughly equal partitioning of the ILAP between the three sources.This consistency is characteristic of most of the results for the other source regions as well (results of the back trajectory analysis for regions 1, 2, 3, and 6 are given in the Supplement).However, a few anomalies are also apparent.This is illustrated in Fig. 5c, which shows results for the Region 4 analysis.While sources for the substantial fire component of the ILAP shown in Fig. 4 are evident, rationalization of the other sources is less clear, in particular for the largest component of the ILAP, soil dust.In this instance, it is likely that very local soil dust, essentially dirt, has been deposited on the snow.And in fact the nature of some of the sampling sites in this region renders this hypothesis quite feasible (see Fig. 2b of Huang et al., 2011).Furthermore, given the necessity of this rationalization for this region, it is also quite likely that very local sources of soil dust impacted numerous sites in the other regions as well, as suggested by Huang et al. (2011) and Wang et al. (2013).

Conclusions
Snow samples obtained at 46 sites across northern China in early 2010 were analyzed for a suite of chemical species and these data collated with the ILAP concentrations measured by Wang et al. (2013).This data set was successfully employed in a receptor model (EPA PMF4.1) to identify the main sources of the ILAP.It was found that ILAP in the snowpack originated from three sources: soil dust (53 %), industrial/urban aerosol (27 %), and biomass and biofuel burning (20 %).While our results suggest a somewhat more predominant role for soil dust in the northeastern area of China than is suggested by the analysis of Wang et al. (2013), this outcome is essentially consistent with the assessments made by Huang et al. (2011) and Wang et al. (2013), which suggested that soil dust might well prove to be the highest contributor to snow ILAP across northern China.This situation is appreciably different from that which we have previously found for the Arctic utilizing the same methodology (Hegg et al., 2009(Hegg et al., , 2010)), but can be easily rationalized in terms of the relative proximity of the two geographic areas to soil dust source regions.Back trajectory cluster analysis coupled with the identification of source areas for the identified ILAP sources was largely consistent with the PMF analysis but suggests that local sources of soil dust are likely very important across northern China.

Fig. 1 .
Fig. 1.Locations of the 46 sampling sites in the 2010 campaign.The different symbols show different regions.Also shown are areas of deserts, of relatively high industrial emissions and of fires.The fires were detected by the MODIS FIRMS system over the snow accumulation period.Note that the industrial emissions are gridded at 30 × 30 resolution.

Fig. 2 .
Fig. 2.Sampling statistics for ILAP and other species used in the PMF analysis.Mean concentrations are given by the dots, the quartile range is defined by the boxes, the medians by the horizontal lines within the boxes and the two standard deviation sample range (95 % percentile) by the error bars.

Fig. 3 .
Fig. 3. Factor profiles for the three-factor solution to the PMF inversion of the China snow data set.

Fig. 4 .
Fig. 4. Average source apportionment for each of the seven regions of northern China discussed in the text.Note that the apportionment is of the light absorption by insoluble impurities in the snow.

Fig. 5 .
Fig. 5. Back trajectory analysis for three of the regions discussed in the text.A seven cluster analysis back trajectory was done on each of the regions with mean trajectories shown for each cluster that encompassed 10 % or more of the back trajectories.(a) Region 7 (Qinghai-Tibet Plateau), 536 total trajectories; (b) Region 5 (industrial northeast), 336 total trajectories; and (c) Region 4 (clean northeast), 400 total trajectories.Note that, as in Fig. 1, the fires shown are those detected over the entire snow accumulation period.