Developing a Black Carbon-Substituted Multimedia Model for Simulating the PAH Distributions in Urban Environments

A multimedia fugacity model with spatially resolved environmental phases at an urban scale was developed. In this model, the key parameter, organic matter, was replaced with black carbon (BC) and applied to simulate the distributions of phenanthrene (Phe), pyrene (Pyr) and benzo[α]pyrene (BaP) in Nanjing, China. Based on the estimated emissions and measured inflows of air and water, the Phe, Pyr and BaP concentrations in different environment media were calculated under steady-state assumptions. The original model (OC-Model), BC-inclusive model (dual C-Model) and improved model (BC-Model) were validated by comparing observed and predicted Phe, Pyr and BaP concentrations. Our results suggested that lighter polycyclic aromatic hydrocarbons (PAHs) were more affected by BC substitution than their heavier counterparts. We advocate the utilization of sorption with BC in future multimedia fate models for lighter PAHs based on the comparison of the calculated and observed values from measured and published sources. The spatial distributions of the Phe, Pyr and BaP concentrations in all phases were rationally mapped based on the calculated concentrations from the BC-Model, indicating that soil was the dominant sink of PAHs in terrestrial systems, while sediment was the dominant sink of PAHs in aquatic systems.

leads to enhanced sorption of planar and aromatic organic compounds [13][14][15] . The formation process of BC is similar to PAHs, as the two compounds are both formed by fossil fuel combustion (traffic, industry, coal, and oil); part of their formation results from the combustion of biomass (forest fires and residential wood burning) 1,16 . BC consists of the combustion-derived carbon fraction, including residues, of initial fuel (char) as well as highly condensed carbonaceous products (soot) 17 , which are ubiquitous in the environment. The fraction of BC to total organic carbon (TOC) ranges from 5-18% in sediment samples around the world 18 , but most PAHs are bound to the BC fraction 19 . Generally, BC contents are approximately 1-15% of the TOC, and therefore, in several cases, BC can be expected to more strongly contribute to overall sorption than all the other organic matter constituents 16 . Prevedouros et al. 12 developed a BC-inclusive multimedia model for predicting the fate of organics and found that lighter PAHs were more affected by BC inclusion than their heavier counterparts.
Based on the above analysis, this study attempted to improve the multimedia urban model (MUM) structure 20 , which included developing and evaluating a spatially resolved solution. We compared the simulated results among the original model (OC-Model), BC-inclusive model (dual C-Model) 12 , and BC-substituted model (BC-Model), which assumed that BC is entirely responsible for the capacity of solids to sorb organic chemicals, to determine the best model to simulate the environmental distribution of PAHs in Nanjing, China. Air, vegetation, water, soil, sediment, and organic film ("pure" film plus particles are covered on impervious surfaces) were included in the model 20 . The model was based on the steady-state, level III fugacity model of Mackay 21 , and the spatially distributed emissions of PAHs in the Nanjing urban areas from January 2012 -December 2013 were estimated and then used as the representative data for contemporary emissions. This study should provide a useful tool for modeling the distribution of PAHs in urban environments.

Methods
Study area descriptions. Nanjing was selected as a prototype city. The city is an ancient capital of six dynasties with a history of more than 2500 years and a profound cultural background. As a main port along the Yangtze River, Nanjing is a complex industrial base dominated by electronics, automobile, and chemical industries. Each environmental medium in Nanjing has been heavily contaminated by many types of pollutants, including PAHs 1,22,23 . The study region was divided into 36 grid cells (4 km × 4 km). Each grid cell was constructed with six connected phases representing air, water, soil, sediment, vegetation, and organic film. Figure 1 shows the location of the study area and model segment.
Model structure. Like other fugacity models, three key variables, including fugacity (f, Pa), fugacity capacity [Z, mol/(m 3 Pa)], and transfer rate coefficients [D, mol/(Pa h)], were employed in this model to describe the environmental distribution and fate processes. Z is specific to the properties of the contaminant and the environmental media. D represents the inter-media transport and transformation processes, including diffusion, advection, and degradation, and is related to the chemical transfer flux. The concentration of the contaminant within a phase (c, mol/m 3 ) is a product of f and Z. The original and developed Z values for each phase are listed in the Supporting Information (Tables S1-S3). In addition, air, water, soil, sediment, vegetation, and organic film were defined as the six bulk phases for modeling the distribution of PAHs in the area. The Supporting Information (Tables S4 and  S5) shows the D values and steady-state mass-balance equations used in the model. The spatial variations of air, water, soil, sediment, vegetation, and organic film were taken into consideration, and these environmental media were divided into 36 individual cells, each 4 × 4 km 2 in area. The major chemical transfer, transformation, and fate Parameterization. The parameters used in the model calculations include various environmental parameters, physical-chemical properties of Phe, Pyr and BaP, and flux matrixes among various media. Most of the data were obtained from the literature or from the model default value 4,24 . Information on the dimensions and areas of the various environmental media for the segments were derived from remote sensing discrimination. The physical-chemical properties of the three chemicals used as input data into this model were derived from Prevedouros et al. 12 , and the K BC (BC-water partition coefficient) and K BC-A (BC-air partition coefficient) of the three chemicals were derived from Prevedouros et al. 2 . The physical-chemical properties of Phe, Pyr and BaP are shown in Table S6. The rain rate and runoff rate were derived from local meteorological and hydrological data. In June 2014, 69 composite surface soil samples were collected from the Nanjing urban areas. The distribution of soil sampling sites is illustrated in Fig. 3. For each 10 m × 10 m sampling site, five soil subsamples (four corners and one center) were taken and bulked together to form one composite sample. All the samples were air-dried at room temperature for one week, sieved to 20-, 60-, and 100-mesh size particles after removing stones, residual roots, and other materials, and then stored in amber glass containers at −4 °C until analysis. The soil organic matter (SOM) contents were determined by the combustion oxidation-titration method 25 . The soil BC contents were determined by the method described by He and Zhang 26 using a Heraeus CHN-O-Rapid elemental analyzer (GmbH, Hanau, Germany). The SOM and soil BC content results for each grid in the study area are shown in Table S7. In this research, the concentration of BC aerosol was derived from Tian et al. 27 . The content of BC in water suspended particulates was obtained from Huang and Zhang 28 . The concentration of BC in sediment was derived from Huang and Zhang 29 .
The air and water flows were described by matrixes of flow rates for their inter-regional movements. Airflow rates were calculated from the annual mean wind speed 30 . The external source concentrations of PAHs in the air and validation data were measured using polyurethane foam passive air samplers deployed throughout the city for a 3 month period in 2014. Twenty water samples were collected from rivers and lakes in the Nanjing urban areas to explore the external source concentrations of PAHs and to construct validation data. The advections in water and air for each grid in the study area are shown in Table S7. The distributions of the air and water sampling sites are illustrated in Fig. 3. The data for the water flow rates and water flux matrixes were collected from the Nanjing City Water Resources Bulletin 31 . Parameter values for the study area summarized in Tables S8, Table S9 and S10 contain the parameter values for chemical transfer.
Estimated emissions. In this study, seven major sources were considered: coking production, primary aluminum production, transport petroleum combustion, non-transport petroleum combustion, industrial coal combustion, domestic coal combustion, and biomass burning. However, biomass burning, such as straw and fire wood burning, has rarely occurred in recent years in our research area. Primary aluminum is also not produced in the Nanjing urban areas, and domestic coal combustion has been replaced by coal gas and natural gas in these urban areas. Therefore, only four sources were considered in this research: transport petroleum combustion, non-transport petroleum combustion, industrial coal combustion, and coking production. The emissions from these sources were then distributed across the grid cells in the model for the gridded region.
The emission estimates in our study were based at the urban street level (known as the township level division in China), but no energy consumption data were available at this level. Zhang et al. 32 found that the consumption of transport petroleum, non-transport petroleum and industrial coal linearly correlated with secondary plus tertiary GDP s (GDP 23 ). Energy consumption data on a provincial scale as well as data of the county and street level during January 2012 -December 2013 were collected from yearbooks [33][34][35][36][37][38] . The main coking production data were collected from official published sources 38 . Thus, the energy consumption for each street was estimated, the data were integrated into each grid using a Geographic Information System, and the emissions were calculated using the emission factors 39 . Figure 4 illustrates the spatial distribution of the Phe, Pyr and BaP emissions in the Nanjing urban areas. The total emissions of Phe, Pyr and BaP in the Nanjing urban areas were estimated as 7.05t, 1.16 and 1.22t during January 2012 -December 2013, respectively.

Sensitivity analysis and uncertainty analysis. Sensitivity analysis is an important methodology used to
test the sensitivity of the model to parameter changes and identify the most influential parameters for the model output. In this research, sensitivity analysis was performed for all the inputs parameters. Each parameter was individually adjusted by ± 10%, and the sensitivity coefficient (S) was calculated by the following formula: where X i (i is the index of the input variable) denotes the input parameter of the model, Y denotes the output of the model when the test parameter is adjusted ± 10%, and Δ refers to perturbations. Monte Carlo simulations were performed to evaluate the overall uncertainty of predictions based on the probability distributions of the input parameters. This method computes results based on repeated random sampling and statistical analysis. The model can be formulated as follows: where g(r) is the random variable, f(r) is the distribution density function, and <g> is the mathematical expectation of g(r). The approximate formula can be converted into: However, most of the inputs had negligible influence on the output of the model. For this reason, all input parameters and their uncertainty were estimated according to their distribution, for which coefficients of variation were calculated to describe the distribution of the parameters based on the literature data. In this study, key input parameters were selected according to the sensitivity coefficient (S) > 0.2 37 . The values of the key input parameters were randomly selected within their respective probability distributions. The Monte Carlo simulations were repeatedly run 10000 times using Oracle Crystal Ball to obtain the distribution of the output.

Results and Discussion
Model Validation. Two procedures were used to test the validity of the model. One procedure assessed the calculated predictions by comparing the observations and predictions of the PAH concentrations in the bulk phase (air, water, soil, vegetation, and sediment). The other compared the prediction results among the BC-Model, dual C-Model and OC-Model.
As shown in Fig. 5, almost no difference in the Phe concentrations in the air phase was observed among the OC-Model, dual C-Model and BC-Model. The maximum difference between the calculated and the measured concentrations was 0.14 log-units for Phe in the air phase. However, large differences in specific cells, such as 2, 17 and 21, were found for Phe in the water phase. The results could be attributed to ignore the PAHs from external water sources when simulating the fate of PAHs because of the water areas are very small in these cells.  (Fig. S1); the available observed data for the Phe concentrations were collected from published data 40,41 . The independently measured mean Phe concentrations were 591.4 ng/g and 5.05 ng/g for vegetation and sediment, respectively, compared to 627.00 ng/g and 1.18 ng/g for the BC-Model, 1.63 ng/g and 0.5 ng/g for the dual C-Model, and 1.62 ng/g and 0.1 ng/g for the OC-Model, respectively. Obviously, the pre- Additionally, no difference in the Pyr concentrations in the air phase among the three models was found, and the predicted values versus the observed values were similar. In the water phase, the distribution of the calculated Pyr concentration was similar to that of the Phe concentration. In the soil phase, the calculated Pyr concentrations in the OC-Model and dual C-Model were the same but lower than the calculated value in the BC-Model, which was closer to the measured values in the most grids. However, their concentration differences fell within 1 log-unit. The calculated Pyr concentration results in the other phases are illustrated in Fig. S2. The differences of the calculated Pyr concentrations in vegetation, sediment and organic film among the three models also fell within 1 log-unit. Compared to the observed values from published data 40,41 , the absolute differences between the calculated and measured mean concentrations in vegetation fell within 1 log-unit, but the difference was greater than 1 log-unit in the sediment phase. In all cases, the difference between the modeled concentrations of the BC-Model, dual C-Model and OC-Model fell within 0.5 log-unit for BaP. Compared to the observed values, almost no difference in BaP modeled concentration was observed in the air phase. The average difference between the calculated and measured concentrations fell within 0.5 and 0.1 log-units for BaP in water and soil, respectively, which are acceptable for this type of study. The calculated BaP concentration results in the other phases are illustrated in the Supporting Information (Fig. S3). Compared to the observed values from published data 40,41 , the absolute differences between the calculated and measured mean concentrations in vegetation and sediment are both approximately 0.5 log-unit. The scatter plots of predicted versus observed values are illustrated in Fig. S4.
Based on the results of the abovementioned analysis, the BC-Model did not cause any considerable variations to the air levels and water levels of the grids containing PAHs from external water sources when simulating the fate of PAHs. However, the BC-Model caused increases in other environmental phase concentrations related to increase solid partitioning, except for BaP. In other words, the difference of the calculated concentration decreased with an increasing molecular weight of PAH in the three models. Our results are similar to the findings of Prevedouros et al. 42 , which indicated that lighter PAHs were more affected by BC inclusion than their heavier counterparts in all cases. In this study, we simulated the distributions and inventories of Phe, Pyr and BaP using the BC-Model, which provided more accurate simulations compared to the observed values from measured and published sources. The predicated results are as follows.
Spatial Distribution. The spatial distributions of the Phe, Pyr and BaP concentrations in all phases were mapped based on the calculated results. As shown in Figs S5-S7, the spatial distributions of the Phe, Pyr and BaP concentrations in air, soil, vegetation, and film are very similar, and the patterns in water and sediment are alike. The higher Phe concentrations in the former group were mainly distributed in the central north-south axis, while the lowest Phe concentrations were distributed around the north of Purple Mountain (segment 23). This phenomenon could be due to traffic and population distribution in the study areas and effected by the chemical properties of Phe. The spatial distribution of Pyr was similar to that of Phe, which is attributed to their similar physical-chemistry properties. As shown in Fig. S7, the BaP concentration patterns in air, soil, vegetation, and film are very similar to the emission data, indicating that the areas at risk of BaP exposure are narrowly distributed around the emission sources. A recognizable difference in segment 4 occurred between the emission data, as shown in Fig. S7_a and S7_b, which is probably due to the influence of advective flow. Compared to the emission map, the spatial distributions of Phe, Pyr and BaP in water and sediment showed unique characteristics. The higher Phe, Pyr and BaP concentrations were mainly distributed in segments containing the Yangtze River flow, indicating that rivers are also an important risk diffusion pathway for PAHs. Inventory. The model predictions for the inventory in each grid are shown in Tables S11-S13. In many cases, most of the total Phe, Pyr and BaP resided in the soil, except for some segmentations mainly located in the river cells, where the total amounts of Phe, Pyr and BaP are predominantly in the sediment. Thus, soil and sediment serve as the predominant PAH sinks, as soil and sediment possessed the highest PAH retention times in the multimedia phases. In addition, the sediment and soil volumes were also important factors for their high storage capacity. The highest Phe, Pyr and BaP concentrations among all phases were in organic film, but the total amounts were low because of the small volume of organic film.
Sensitivity and uncertainty of the model. As no differences in the source, distribution and contribution of the input parameters were observed between different segments, the sensitivity and uncertainty analysis results in one cell could be representative on the entire region 43 . Segment 27 was selected as an example for the sensitivity and uncertainty analyses. The sensitivities and uncertainties of the modeled concentrations in all the phases to the input parameters of Phe, Pyr and BaP in the OC-Model, dual C-Model and BC-Model were analyzed. The Supporting Information (Tables S14-S16) shows the absolute values of the sensitivities to the key parameters controlling the Phe, Pyr and BaP concentrations in all phases. The key parameters for the various phases were apparently different. The soil phase was affected by the highest number of key parameters in all phases, while the water was least affected. The key parameters were selected for Monte Carlo simulations according to the results of the sensitivity analysis. All the distributions of the key parameters were well fitted to log-normal distributions according to comparisons between the statistics of the fitted distribution. The results of the uncertainty analysis for the Phe, Pyr and BaP concentrations in all phases are illustrated in the Supporting Information (Tables S17-S19). The model outputs of the Phe, Pyr and BaP concentrations in soil and sediment had apparently higher uncertainties than those in the other phases due to the increased numbers of key parameters in the soil and sediment phases. Little difference of the key parameters was observed for the three chemicals among the BC-Model, dual C-Model and OC-Model. In addition, large differences of the value ranges in soil, vegetation, film and sediment phases were found for lighter PAHs in the BC-Model, dual C-Model and OC-Model, but the value ranges were similar for heavy PAHs. However, the model validation results are the same even under extreme conditions when the two extreme uncertainty values fell within 1 log-unit.

Conclusion
In this study, a spatially resolved MUM model was developed. The key parameter, organic carbon, was displaced with BC, and model predictions were validated by comparing the observations and predictions of PAH concentrations. Our results suggested that lighter PAHs were more affected by BC substitution than their heavier counterparts. Based on the comparison of the calculated and observed values from measured and published sources, we advocate the utilization of sorption with BC in future multimedia fate models for lighter PAHs. The spatial distributions of the Phe, Pyr and BaP concentrations in all phases were rationally mapped based on the calculated concentrations of the BC-Mode, indicating that soil was the dominant sink of PAHs in terrestrial systems, while sediment was the dominant sink of PAHs in aquatic systems.