PAH contamination in Beijing’s topsoil: A unique indicator of the megacity’s evolving energy consumption and overall environmental quality

To improve its air quality, Beijing, the capital of China, has implemented high-cost pollution control measures mainly focused on shifting its energy mix. However, the effectiveness of these measures has long been questioned, especially given the recent problem of severe haze. The main study objectives are to achieve independent, although indirect, information on Beijing’s air pollution by measuring the level of polycyclic aromatic hydrocarbon (PAH) contamination in topsoil and to examine how soil contamination reflects energy consumption. Soil sampling data from two years, 2004 and 2013, were used. The key findings are as follows: 1) although the total PAH content in the topsoil did not significantly decrease from 2004 to 2013, the composition changed considerably; 2) as of 2013, vehicle emissions replaced coal combustion as the leading source of soil PAHs, which validates the existing policy measures regarding vehicle purchasing and traffic volume; 3) the regional transport of atmospheric pollutants, as indicated by the contribution of coking sources in 2013, is not negligible; and 4) appropriate policy measures are needed to control the growing practice of burning biomass. Overall, this study demonstrates that the PAH contamination in topsoil represents an informative indicator of Beijing’s energy consumption and overall environmental quality.


PAH concentrations.
Our study focused on the downtown and suburban districts of Beijing. All 16 PAH compounds were detected in the soil samples (45 samples in each campaign) collected within the study area. In the downtown districts, the median concentrations of ∑ 16 PAHs in 2004 and 2013 were 699.44 ng/g and 500.27 ng/g, respectively. In the suburban districts, the median concentrations were 183.20 ng/g and 177.11 ng/g, respectively, indicating a much lower level of PAH contamination than in the downtown districts. Figure 1 illustrates the concentrations of the 16 compounds in the two periods. As the concentrations have highly skewed distributions, the median concentrations rather than mean values are presented in the figure. To examine whether the concentrations had a statistically significant change from 2004 to 2013, Mann-Whitney U tests were performed, and the results are presented in Table 1. In the table, "All" indicates that both downtown soil samples and suburban soil samples were considered, while "Downtown" and "Suburban" refer to the respective subsets of samples. Although the content of ∑ 16 PAHs did not change significantly (see Table 1), the 16 PAH compounds showed different trends over time (see Fig. 1 and Table 1). Most of the 3-ring and 4-ring compounds increased over time, whereas most of the 5-ring and 6-ring compounds decreased over time. However, the 2-ring compounds showed no general tendencies. Because 2-ring PAHs are the least persistent compounds, soil samples collected in the field campaigns may not reflect the long-term average conditions of such compounds. The implications of the changes are two-fold: first, the source characteristics of soil PAHs varied, and second, the environmental risk associated with the soil contamination changed. The following subsections further discuss these issues.
Source apportionment results. The source profiles derived for both 2004 and 2013 are given in Fig. 2.
Each subplot in the figures illustrates the loading of a source factor on different PAH compounds. The PAHs with relatively high loadings are represented by the darker bars. The interpretation of the source factors was based on the prior knowledge that coal combustion, vehicle emissions, coking and biomass burning are the major origins of PAHs in Beijing. In 2004 (Fig. 2a), Factor 1 has high loadings on ANT, PYR, BaA, CHR and BaP. These compounds, as well as FLA, BbF and BkF, are markers for coal combustion [32][33][34][35][36] . Therefore, Factor 1 in Fig. 2a can be interpreted as the coal combustion source. Factor 2 is dominated by IcdP, DahA and BghiP. Both BghiP and IcdP have been identified as typical tracers of a vehicular source of PAHs 35,36 , and DahA is also associated with traffic emissions 37 . Factor 2 in Fig. 2a can therefore be attributed to vehicle emissions. Factor 3 is mostly associated with FLO, which is a typical marker for coke oven sources 34,36 . Therefore, Factor 3 in Fig. 2a represents coking sources. Factor 4 has high loadings of ACY, which is the dominant PAH compound released during biomass burning 38,39 . Factor 4 in Fig. 2a can therefore be interpreted as the biomass burning source. Similarly, Factors 1, 2, 3 and 4 in Fig. 2b can be interpreted as coal combustion, vehicle emissions, coking and biomass burning, respectively. Figure 3 shows the contributions of the different sources as quantified by the PMF analyses. It is clear that the source composition experienced a significant change from 2004 to 2013. In 2004, the leading PAH source was coal combustion (42.8%), followed by vehicle emissions (28.2%), coking (16.8%) and biomass burning (12.2%). However, in 2013, the contributions of vehicle emissions (38.8%) and biomass burning (18.3%) increased, while those of coal combustion (30.2%) and coking (12.7%) decreased. The order of the four sources changed to vehicle emissions > coal combustion > biomass burning > coking. Spatial distributions of risk level. We used the carcinogenic equivalent concentration of BaP (denoted as BaP eq ) as a risk indicator. The BaP eq concentration mainly reflects the magnitude of coal combustion and vehicle emission sources. The average BaP eq values of all samples were 139.48 ng/g and 310.42 ng/g for 2013 and 2004, respectively, indicating a significant reduction of the potential risk associated with soil contamination. Figure 4 illustrates the spatial distributions of risk level in downtown Beijing based on spatially interpolated BaP eq values. The interpolation was finished in ArcGIS for Desktop (version 10.2.1) using the inverse distance weighting method. As shown in the figure, the risk level substantially decreased from 2004 to 2013 for most of the downtown area, and its spatial pattern changed significantly. In 2004, two contamination hotspots were in the southeastern and northern regions (Fig. 4a). The south-eastern region was once an industrial center, and PAH Ring number All Downtown Suburban Table 1. Results of the Mann-Whitney U tests. "↓ ↓ " and "↑ ↑ " indicate highly significant (p-value equals to 0.01) decreases and increases between 2004 and 2013, respectively; "↓ " and "↑ " indicate significant (p-value equals to 0.05) decreases and increases, respectively; and "-" indicates no significant changes.
the northern region is a major residential area. Intensive coal consumption occurred in both regions during the early 2000s. The two hotspots disappeared in 2013 (Fig. 4b), which should be attributed to Beijing's pollution prevention efforts, as explained below. Many high-pollution plants in the southeastern region, such as Beijing Coking and Chemical Works, Beijing Chemical Works, and Corrosion-proof Railway Sleeper Plant, were relocated outside of the downtown area before the 2008 Olympic Games. On the other hand, during the past decade, the primary energy sources for residential heating in downtown Beijing have shifted away from being dominantly coal to what now being mainly natural gas and electricity. The new hotspot in the southwestern region may be related to increased vehicle emissions in this region. A possible cause of this increase could be the establishment of a second-hand vehicle market in 2002. This market has developed into the largest market of used vehicles in Beijing. The traffic volume in this region has substantially increased due to the expansion of this market, which may explain the new hotspot in 2013.

Discussion
The PAH contamination in topsoil can be viewed as an indicator of local and regional energy consumption. According to the Beijing Statistical Yearbooks (     Many studies have been conducted on the atmospheric PAHs in Beijing from 2004 to 2013. Source apportionment has also been attempted. A study 11 found that fossil fuel uses (mainly coal combustion and vehicle emissions) contributed approximately 65% of the atmospheric PAHs in 2004. A study in 2008 44 estimated that coal combustion and vehicle emissions accounted for 69% of atmospheric PAHs. A study in 2011 45 reported that the combined contribution of coal combustion and vehicle emissions was 83% during non-heating seasons. Furthermore, as discussed in the introduction, atmospheric PAHs are mainly distributed in fine particles smaller than 2.5 μ m. The UNEP and BEPB's joint report, as mentioned in the introduction, estimated that motor vehicle and coal consumption contribute 31.1% and 22.4% of PM 2.5 in Beijing, respectively. All the above results are consistent with our findings on soil PAHs (Fig. 3). In addition, we analyzed the correlation between topsoil and atmospheric PAH content in China's large cities based on data from the literature [46][47][48][49][50][51] . Figure 6 demonstrates a strong correlation and a clear logarithmic relationship between PAH content in the soil and air. However, because the sample size is small, the regression function presented in Fig. 6 may not be considered a general relationship; more data are needed to derive a function that is widely applicable. The above evidence reinforces the theory that the PAH contamination in topsoil is closely related to air pollution and could provide valuable, independent information about air pollution.
In summary, this study, based on soil sampling campaigns in 2004 and 2013, investigated the PAH contamination in Beijing topsoil over the past decade. The study results show that although the total PAH content in the topsoil did not significantly decrease from 2004 to 2013, the PAH composition changed considerably, decreasing environmental risk. The source apportionment analyses revealed a notable decrease in the contribution of coal combustion to PAH contamination as well as a significant increase in the contribution of vehicle emissions. Vehicle emissions have become the leading cause of PAH contamination in soil. The source apportionment results clearly reflect the change in Beijing's energy mix and are consistent with previous studies on Beijing's atmospheric PAHs. Therefore, we argue that the evolution of PAH contamination in topsoil is indicative of the long-term variation of Beijing's air quality.
The study results have many important policy implications. First, from the perspective of soil contamination, Beijing's tremendous efforts to control pollution over the past decade, mainly aimed at shifting its energy mix, had a positive effect. Second, although our investigation into soil contamination did not directly address the issue of PM 2.5 , it offered a unique perspective that could be used to determine whether vehicle emissions or coal is a more important source of Beijing's PM 2.5 . Our results support the viewpoint that vehicle emissions are currently the most important combustion source. We further argue that the control measures on vehicle purchasing and traffic volume in Beijing, which have been highly controversial, are relevant although not necessarily cost effective. Third, straw burning is a very traditional waste disposal practice in rural China and is a major form of biomass burning in Beijing and its vicinities. Whether straw burning significantly contributes to Beijing's haze events has been of great public concern. Our study offers indirect evidence of the increasing effects of straw burning (refer to Fig. 3). Appropriate policy measures are therefore desired to control these burning practices. Fourth, although coking has been eliminated in Beijing, its influence on PAH contamination is still considerable due to the rapid growth of coke production in Beijing's surrounding areas. An important implication of this finding is that local efforts may not be sufficient to tackle Beijing's pollution problems, and regional collaboration is important.
Overall, our study demonstrated that PAH contamination in topsoil serves as a key indicator of a megacity's energy consumption and overall environmental quality and should therefore receive continued investigation.

Study area and sampling locations.
Beijing is situated in the north of the North China Plain (Fig. 7) with a sub-humid continental monsoon climate. Its mean annual temperature is approximately 11.6 °C, and its mean annual precipitation is approximately 546 mm. Among the sixteen administrative districts of Beijing, the central six districts (the blue polygons in Fig. 7) represent the downtown area, and the other ten districts are suburban and rural areas. The total population of Beijing was 20.69 million in 2012, of which 59% was in the downtown districts. As the population, traffic, and industrial activities are mainly distributed in Beijing's plain area, this study considered the six downtown districts and the four suburban districts (the red polygons in Fig. 7) adjacent to them.
The soil samples were collected in two field campaigns. The later campaign was organized in June 2013. A total of 45 topsoil samples were obtained, 28 of which were within the downtown districts and 17 of which were from the four suburban districts (see Fig. 7). Most of the samples were collected in the plain area, and only two were collected in the low-mountain area near the edge of the plain. The sampling locations cover representative land uses in Beijing, including residential areas (10 samples), transportation areas (i.e., large roads, 7 samples), commercial areas (5 samples), recreational areas (i.e., urban parks, 5 samples), industrial sites (8 samples), agriculture land (8 samples) and forests (2 samples). The earlier campaign was conducted from April to June 2004, and its details were reported in a previous study 47 . The 2004 campaign investigated a much larger area than that of our study. From the 2004 study, we considered only the 45 samples from the areas that were also evaluated in the 2013 study: 20 in the downtown districts and 25 in the suburban districts. Compared with the 2013 campaign, the sampling locations in 2004 were more uniformly distributed spatially, and the major land uses were well represented. Both sample sets adequately covered the entire study domain, and they comprise a unique two-year dataset.
Collection and processing of soil samples. The soil sample collection methods and processing used in 2013 were consistent with those of the 2004 study 47 . The sampling was conducted on dry days with no significant antecedent rainfall events. In the mostly impervious downtown area, samples were collected from green spaces. Our pre-sampling investigation ensured that all of the sampling plots were not notably disturbed for at least one year. At each sampling plot, plant materials were first removed from the soil, and soil samples of the top 10 cm were taken from multiple spots (four corners plus the center) and mixed. The collected samples were stored in stainless steel containers and transferred back to the laboratory in a timely manner. In the laboratory, the samples were air-dried at room temperature. Any large pieces of debris (e.g., gravel or plant residues) were removed from the dried samples. The samples were then ground sufficiently fine to pass through a 1-mm stainless steel sieve. After pretreatment, the soil was air-dried at room temperature and then mixed together completely. The soil was stored in sealed stainless steel containers for further processing. Ten grams of soil sample mixed with 10 g of baked anhydrous sodium sulfate was extracted using accelerated solvent extraction (Dionex ASE 300). Extractions were performed using 34 mL of solvent (1:1 mixture of hexane and dichloromethane). The extracts were concentrated, solvent-exchanged to hexane, and passed through a silica column chromatography. The column was eluted with 25 mL hexane and then 50 mL of dichloromethane:hexane (3:2) mixture. The combined collected solvent was evaporated using a rotary evaporator. The final extracts were concentrated to 1 mL, and known quantities of internal standards 2-fluorobiphenyl and p-terphenyl-d 14 were added and then transferred into vials before analysis.
Laboratory analysis and quality control. All the samples were analyzed for the following 16 USEPA pri- The laboratory sample analysis and quality control used in 2004 were described in a previous paper 47 and, the same methods were used for the samples in 2013. The PAH concentrations in the extracts were determined by an Agilent 6890 gas chromatograph (Agilent, USA) equipped with a 5973N mass selective detector. The identity of each PAH was confirmed using a standard PAH mixture (610/525/550 in methanol from Chem Service, USA) containing the 16 PAHs. The standard PAH mixture was analyzed in the GC/MS at full scan mode. An HP-5 silica fused capillary column (30 m × 250 μ m inner diameter × 250 μ m film thickness) was used with a carrier gas of helium at a constant flow rate of 1 mL min −1 . The oven temperature was programmed to increase from 60 °C to 280 °C at a rate of 5 °C min −1 and maintained at 280 °C for 20 min. All solvents, hexane and dichloromethane, were HPLC-grade (CNW Technologies GmbH, Germany). Silica gel (60-80 mesh, Beijing Chemical Reagent Co., China) was heated at 450 °C for 4 hours, kept in a sealed desiccator, and reactivated at 130 °C for 16 hours prior to use. Sodium sulfate was pre-baked at 650 °C for 4 hours and then kept in a sealed desiccator. The working standard solution was prepared by diluting a stock standard from a commercial mixed standard (J&K Chemical Ltd., USA) with hexane.
All data were subject to strict quality control procedures. Quantification was performed using an internal calibration method. The samples were spiked with five surrogate standards (naphthalene-d 8 , acenaphthene-d 10 , phenanthrene-d 10 , chrysene-d 12 , and perylene-d 12 ), and the surrogate recoveries were 65% to 127%. For each batch of samples, two procedural blanks were processed, and all data were blank-corrected. The relative percent difference between individual PAH compounds identified in method duplicate samples was < 15%. The detection limits were 0.30 ng/g to 0.60 ng/g based on a 10 g soil sample.
Positive Matrix Factorization. Positive matrix factorization (PMF) is a principal component analysis (PCA)-based receptor model with non-negative constraints: where X n×m is a matrix of sample concentrations that can be decomposed into G n×p , a matrix of source contributions, and F p×m , a matrix of source profiles, plus a residual matrix, E n×m . n and m denote the number of samples and the number of pollutant species, respectively, and p is the number of influential sources. The physical nature of the influential pollution sources can be interpreted from F p×m based on prior information such as measured profiles of specific sources and emission inventories. Each row in F p×m represents an influential source (i.e., a factor), and the sum of the row elements (i.e., factor loadings) quantifies the overall contribution of the source. PMF minimizes the following objective function 52 : where S ij represents the "uncertainty" in the jth species for the ith sample and E ij (i.e., elements of the residual matrix E n×m ) represents the error corresponding to S ij .
In this study, the EPA PMF 5.0 program 53 was used to conduct PMF analyses. The program calculates s ij based on the method detection limit (MDL) of targeted pollutant species as: Because NAP is a relatively volatile and unstable compound, it was excluded from the source apportionment, which was aimed at determining the long-term characteristics of PAH sources. Two PMF models were established for the sample sets of 2004 and 2013. p was varied between 3 and 6 and eventually set to 4, for which the Q values calculated by Eq. (2) (429.9 for 2013 and 438.3 for 2004) best approximate the theoretical Q value (equal to 435) calculated by Eq. (4).

Risk levels.
We used the carcinogenic equivalent concentration of BaP (denoted as BaP eq ) as a risk indicator 56 .
BaP eq can be calculated as follows: where TEF i is the toxic equivalency factor for the ith PAH compound (i.e., PAH i ). The TEF values are 0.001 for NAP, ACY, ACE, FLO, PHE, FLA and PYR; 0.01 for ANT, CHR and BghiP; 0.1 for BaA, BbF, BkF and IcdP; and 1 for BaP and DahA 57 . In China, the existing Environmental Quality Standards for Soils (GB15618-1995) does not have content limits for soil PAHs. A new version (GB15618-2008) of this national standard is currently under discussion but has not been officially released yet. This unreleased standard defines four content limits for BaP in soils: I) BaP ≤ 10 ng/g (background condition); II) 10 ng/g < BaP ≤ 100 ng/g (good for agricultural uses); III) 100 ng/g < BaP ≤ 500 ng/g (good for residential land uses); and IV) 500 ng/g < BaP ≤ 1000 ng/g (good for industrial and commercial uses). In this study, we compared the derived BaP eq values to these limits and defined five levels of risk: I) BaP eq ≤ 10 ng/g (no risk); II) 10 ng/g < BaP eq ≤ 100 ng/g (low risk); III) 100 ng/g < BaP eq ≤ 500 ng/g (medium risk); IV) 500 ng/g < BaP eq ≤ 1000 ng/g (high risk); and V) BaP eq > 1000 ng/g (extremely high risk).