Multivariate Analysis for Assessing Sources , and Potential Risks of Polycyclic Aromatic Hydrocarbons in Lisbon Urban Soils

Urban soils quality may be severely affected by polycyclic aromatic hydrocarbons (PAHs) contamination, as is the case of Lisbon (Portugal). However, to conduct a risk assessment analysis in an urban area can be a very difficult task due to the patchy nature and heterogeneity of these soils. Thus, the present study aims to provide an example on how to perform the first tier of a risk assessment plan in the case of urban soils using a simpler, cost effective, and reliable framework. Thus, a study was conducted in Lisbon to assess the levels of PAH, their potential risks to the environment and human health, and to identify their major sources. Source apportionment was performed by studying PAHs profiles, their relationship with potentially toxic elements, and general characteristics of soil using multivariate statistical methods. Results showed that geostatistical tools are useful for evaluating the spatial distribution and major inputs of PAHs in urban soils, as well as to identify areas of potential concern, showing their usefulness in risk assessment analysis and urban planning. Particularly, the prediction maps obtained allowed for a clear identification of areas with the highest levels of PAHs (close to the airport and in the city center). The high concentrations found in soils from the city center should be a result of long-term accumulation due to diffuse pollution mostly from traffic (through atmospheric emissions, tire debris and fuel exhaust, as well as pavement debris). Indeed, most of the sites sampled in the city center were historical gardens and parks. The calculation of potential risks based on different models showed that there is a high discrepancy among guidelines, and that risks will be extremely associated with the endpoint or parameters used in the different models. Nevertheless, this initial approach based on total levels was useful for identifying areas where a more detailed risk assessment is needed (close to the airport and in the city center). Therefore, the use of prediction maps can be very useful for urban planning, for example, by crossing information obtained with land uses, it is possible to define the most problematic areas (e.g., playgrounds and schools). Minerals 2019, 9, 139; doi:10.3390/min9030139 www.mdpi.com/journal/minerals Minerals 2019, 9, 139 2 of 17


Introduction
Polycyclic Aromatic Hydrocarbons (PAHs) are a large group of semivolatile, chemically stable, and hydrophobic compounds which may be highly persistent in the environment, accumulating in soils and through the food chain.Therefore, the contamination of soils by PAHs is mainly a result of long-term pollution, and in areas with a continuous input of these contaminants, even though at low levels, it may result in high concentrations in soils as a result of accumulation over the years [1][2][3][4].Due to the carcinogenic and/or mutagenic potential of PAHs, the contamination of soils may cause deleterious effects on ecosystems, and may also jeopardize human health [2].
Typically, the first step for assessing soil quality is based on the comparison of the total content of contaminants with threshold concentrations, in order to calculate their potential hazard [2,5,6].At this initial stage of risk assessment, multivariate and geostatistical tools can be very useful for site characterization, for example, to evaluate the extent of contamination and identify major sources.Indeed, source apportioning is an essential step for controlling and/or reducing the inputs of contaminants to soils.
Soils, and especially urban soils, can be very heterogeneous, and consequently their properties and concentrations of pollutants may vary remarkably over short distances, making it difficult to have a homogeneous data set.In addition, due to their patchy nature it is not always possible to sample at the desirable location, and consequently it can be difficult to obtain a "real picture" of the area.However, unknown concentrations can be estimated from data taken at specific locations by applying geostatistical procedures such as kriging interpolation [7].This has been successfully used to produce illustrative contour plots of contaminant distribution on the basis of scattered observed concentration data [8][9][10][11][12].Therefore, kriging interpolation of spatial data, together with multivariate statistical methods can be useful tools for extracting additional information from the data set, thus reducing uncertainties and minimizing costs.In addition, geostatistical methods can be very useful in risk management, for example, identifying areas potentially affected considering a specific risk model.
A previous study showed that Lisbon urban area can represent a potential risk to human health due to high levels of PAHs existing in some areas [3].Therefore, in our present study, we aimed to better understand the spatial distribution and major inputs of these compounds in Lisbon urban soils, by using a data set with a higher number of samples, which allowed the development of prediction maps.In addition, we provided a practical example on how to perform the first tier of a risk assessment plan in the case of urban soils using the framework proposed by Cachada et al. [2].In addition to a simple comparison with soil screening guidelines, other specific models for assessing risks to the environment or to human health were also used.

Materials and Methods
The study was undertaken in Lisbon urban area, which is the biggest Portuguese town, and has been well described [3,[13][14][15].A total of 97 topsoil samples (0-10 cm) were collected considering different land uses (roadsides, airport, ornamental gardens, parks and open spaces, schoolyards, and playgrounds), and a broad geographical coverage of the urban area (Figure 1).The sampling location map was produced with ArcGis ® Software (version 9.3, ESRI, Redlands, CA, USA).
Univariate statistical methods and the 95% upper confidence level (UCL) of the mean (used for human health risk assessment) were calculated as described by Cachada et al. [3].Factor analysis (FA), was used to obtain underlying factors describing the covariance among variables.The factors extracted resulted in groups with common characteristics, and thus interpreting the data in accordance with their hypothetical origin (natural, anthropogenic or mixed).Due to differences in units of measurement, FA was conducted in a correlation matrix.The "Kaiser-Meyer-Olkin Measure of Sampling Adequacy test" and the "Bartlett test of Sphericity" were performed in order to measure if the distribution of results was adequate to conduct FA.The extraction method was the principal component analysis (PCA) and the rotation method was the Varimax with Kaiser Normalization.Missing values were considered as equal to half of the detection limit.Statistical tests were performed in SPSS© (V.24, IBM, New York, NY, USA), with exception of the UCL (95%) that was obtained using the statistical package ProUCL ® available by United States Environmental Protection Agency (USEPA).
The variogram (γ(h)) technique, which is the basic instrument of geostatistics, was used to measure the spatial variability of variables [7].The ordinary kriging interpolation analysis with the spherical model was used to construct the filled contour maps, using ArcGis ® Software (version 9.3).
The calculation of risks was performed as described by Cachada et al. [2].Briefly, as a first approach of risk assessment, simple hazard quotients (HQ) were determined by comparing the total concentrations in soils (Csoil) with the acceptable soil screening level (SSL).In order to calculate ecological risks, the Toxic Units (TU) were obtained by calculating the ratio between the concentration in soils and the predicted no effect concentration (PNEC).The toxic units of the mixture (TUm), were also calculated assuming the concentration addition model (CAM).We also applied the Dutch methodology in which cumulative risks of mixtures in one sample are expressed as multisubstance potentially affected fraction (msPAF).This model allows quantifying the fraction of species that is expected to be locally affected beyond the selected effect level.Finally, the Toxic Pressure of the mixture (TPm) was also calculated using the approach suggested by Jensen and Mesman [5], and by Cachada et al. [2], which assumes a dissimilar toxic mode of action for all compounds, and therefore the response addition model (RAM) was applied.
Regarding the human health risk assessment, first, the PAH toxic potency of soil samples, expressed as BAP-equivalents (BAPeq), was calculated as suggested by the Canadian Council of Ministers of the Environment (CCME) [17].Cancer and mutagenic risks associated with PAHs exposure were calculated based on the deterministic approach from USEPA Risk Assessment Guidance for Superfund (RAGS) methodology [18,19], and previously described by Cachada et al. [2,3].The same routes of exposure were considered (ingestion, dermal absorption, and inhalation of particulates), as well as the land use scenarios (residential, recreational, and occupational), and the exposure parameters.Results are expressed as total lifetime cancer risk (TLCR), which corresponds to the sum of risk for the three exposure routes.The target risk value considered was one-in-one-million (10 −6 ).

General Characterization and Levels of PAHs in Lisbon Urban Soils
Results for general parameters can be found in Table S1, whereas PTE results can be found in Table S2.These results remain very similar to those documented in a previous study [13], where it was concluded that soil properties are in line with the type of soil.In the same study it was observed that Lisbon soils had very high concentrations of elements that normally occur in basaltic or calcareous rocks such as Ca, Co, Cr, and Ni [20], whereas the relatively high concentrations of Cu, Pb, Zn, and Hg were mostly due to anthropogenic sources (mainly traffic and industry).Concentrations of the Sum16PAHs obtained in surface soils ranged between 6.3 and 73,395 µg kg −1 , with a mean value of 2645 µg kg −1 and a median of 559 µg kg −1 (Figure 1 and Table 1).Compared to results from our previous study conducted in the same area (Table 1), the mean value increased dramatically, mostly due to the presence of two hotspots located at the airport (sample 96 and 97; Figure 1 and Figure S1).Regarding other Portuguese urban areas (Table 1), the higher levels found in Lisbon can be explained by their long-term accumulation in urban soils and the presence of a high number and more intense sources (e.g., traffic and industry).Due to the above-mentioned hotspots found in this study, Lisbon area showed the highest maximum values compared to other cities around the world (Table 1).Consequently, the mean value was also one of the highest, being exceeded only by Glasgow.Even so, the median concentration was still slightly lower than those observed in other European cities, reflecting the presence of several hotspots, in addition to the ones mentioned above (Figure S1).On the other hand, the median concentration observed in Lisbon was higher than those found in most Asian cities and the climate has been addressed as the major reason for the differences found [3].Indeed, several authors [21,22] argue that Asian cities have a warm and more humid climate, which enhances biodegradation, volatilization, and photo-oxidation, leaching of semivolatile organic compounds into the groundwater.The use of the ordinary kriging interpolation method (Figure 1) allowed a better understanding of the spatial distribution of PAHs in Lisbon urban soils.It is clearly shown that areas with the highest PAHs concentrations were the airport and the area close to the river (city center), which are the oldest parts of the city, with higher population density.Concentrations of the Sum16PAHs obtained in surface soils ranged between 6.3 and 73,395 µg kg −1 , with a mean value of 2645 µg kg −1 and a median of 559 µg kg −1 (Figure 1 and Table 1).Compared to results from our previous study conducted in the same area (Table 1), the mean value increased dramatically, mostly due to the presence of two hotspots located at the airport (sample 96 and 97; Figure 1 and Figure S1).Regarding other Portuguese urban areas (Table 1), the higher levels found in Lisbon can be explained by their long-term accumulation in urban soils and the presence of a high number and more intense sources (e.g., traffic and industry).Due to the above-mentioned hotspots found in this study, Lisbon area showed the highest maximum values compared to other cities around the world (Table 1).Consequently, the mean value was also one of the highest, being exceeded only by Glasgow.Even so, the median concentration was still slightly lower than those observed in other European cities, reflecting the presence of several hotspots, in addition to the ones mentioned above (Figure S1).On the other hand, the median concentration observed in Lisbon was higher than those found in most Asian cities and the climate has been addressed as the major reason for the differences found [3].Indeed, several authors [21,22] argue that Asian cities have a warm and Minerals 2019, 9, 139 5 of 17 more humid climate, which enhances biodegradation, volatilization, and photo-oxidation, leaching of semivolatile organic compounds into the groundwater.The use of the ordinary kriging interpolation method (Figure 1) allowed a better understanding of the spatial distribution of PAHs in Lisbon urban soils.It is clearly shown that areas with the highest PAHs concentrations were the airport and the area close to the river (city center), which are the oldest parts of the city, with higher population density.
The location of some of the specific sources of PAHs, identified previously by Cachada et al. [21] (railway, crematoriums, hospital waste incinerator, the docks, and shipyard activity), are shown in Figure 1.Nevertheless, it is known that traffic is one of the most important diffusion sources of PAHs into urban soils, through atmospheric emissions, tire debris and fuel exhaust, as well as pavement debris.Due to the tendency of PAHs to accumulate in soils and to their persistence, the high concentrations found in soils from the city center are essentially a result of long-term accumulation in historical parks and gardens due to point and diffuse pollution.Most sites sampled in the city center were historical gardens and parks such as for example samples 9-18, 22-23, and 54-56.a Sum15PAHs; b Sum23PAHs; ND = Not Detected.

Source Apportionment of PAHs in Lisbon Urban Soils
The observed PAHs profiles (Figure S2) agreed with those from our previous study [3], and also with those found in literature for urban soils (e.g., Chung et al. [26]; Ma et al. [27]).The high molecular weight (HMW) PAHs (≥4 rings) dominated the profiles, with the 5 + 6-ring PAHs showing a median value of 47%, and the 4-ring compounds a median of 43%; whereas the 2 + 3-ring group only showed a median percentage of 9%.The most abundant compounds were FLA, PYR, BBF, and BGHI (Figure S2).FLA and PYR are usually associated with fossil fuel (especially diesel), coal and biomass combustion, but also with tire debris [28,29].BBF, in addition to fossil fuel, can be related with asphalt debris, whereas BGHI is normally related with light-duty vehicle emissions [28].
Nevertheless, it should be kept in mind that these composition profiles can be weathered due to volatilization or degradation of low molecular weight (LMW) PAHs (<4 rings).Moreover, HMW tend to accumulate close to emission sources, whereas LMW PAHs can be transported over longer distances.Therefore, soil profiles may not reflect emission profiles [28,[30][31][32].For this reason, and due to the presence of mixed sources in urban areas, the use of ratios for selected pairs of isomers is commonly applied to identify sources of PAHs.However, despite the several advantages that have been attributed to this methodology (confounding factors such as differences in volatility, water solubility or adsorption are minimized), some authors refer that they have limited usefulness when addressing specific PAHs sources in soils [3,31,33].The most commonly used isomer ratios, due to their stability, are FLA/(FLA + PYR) and IND/(IND + BGHI) [33,34].In the present study, the ratio FLA/(FLA + PYR) pointed to coal and biomass combustion in 96% of samples, but according to the ratio IND/(IND + BGHI) only 28% of samples had this origin.Therefore, it is not possible to identify specific sources based on these ratios.Nevertheless, considering all the ratios typically used (those previously referred together with ANT/(ANT + PHE) and BAA/(BAA + CRY)), they were coincident pointing to combustion sources.
Significant correlations (ρ < 0.01) were observed between the SumPAHs and As, Ca, Cu, Hg, Pb, and Zn, with the highest value observed for Pb (0.80) and the lowest for Ca (0.45).Factor analysis (FA), including the Sum16PAHs and PTEs, previously identified as having a telluric (Al, Fe, Co, Cr, and Ni) or anthropogenic (Ca, Cu, Hg, Pb, and Zn) origin [13], resulted in two factors (Table 2).Arsenic was not considered in the present FA due to its low communality value.
Factor 1 can be described as the telluric origin of PTEs, since it has high loadings of Al, Fe, Cr, and Ni; whereas factor 2 is related to the anthropogenic origin of both PTEs and PAHs.The presence of Ca in factor 2 can be related to its origin in coal combustion, cement production, and incineration.Copper and Zn, which are normally associated with traffic contamination in urban areas, also showed some influence on factor 1 and this behavior is reflected on the loading plot (Figure S3).In this loading plot of both factors, it is possible to clearly observe how contaminants are grouped, and to note the presence of two groups for anthropogenic elements: one formed by Hg, Pb, and PAHs, and the other formed by Ca, Cu, and Zn.Despite the common anthropogenic sources of these PTEs and PAHs, the separation of these two groups could be related to other specific sources or to the influence of soil properties on the retention of these contaminants.
Figure 2 shows the spatial distribution of both factor scores, being clear which are the areas most affected by each source.Southwestern was the area with high factor scores for factor 1 (Figure 2a), and high levels of Co, Cr, and Ni (high loadings on factor 1) have been previously related to the basalts from the Lisbon Vulcanic Complex [13].The scores of factor 2 (Figure 2b) had a similar spatial distribution to the Sum16PAHs (Figure 1), with two major hotspot areas clearly identified: the airport and the city center.The reasons for the distribution observed for factor score 2 are the same as referred for PAHs, since the major sources of PAHs and anthropogenic PTEs are common: traffic, railway, shipyard industry, former industry, incineration, and crematoriums [3,13].
In order to understand the role of soil properties on the distribution of PTEs and PAHs, the relationship between them was studied.The SumPAHs were significantly correlated (ρ < 0.01) with TC, OC, OM, and silt, yet values were low (the strongest correlation was observed for TC, with a value of 0.483).A new FA was performed (Table 2), this time including two parameters considered representative of soil properties (OC and clay), and results remained similar to the first FA.However, it is interesting to note that clay showed a high loading on factor 1 (telluric inputs), whereas OC had similar loadings in both factors, strengthening its role on retention of PTEs (both with telluric or anthropogenic origin) and PAHs.Therefore, both soil properties and specific sources should influence the spatial distribution of PTEs and PAHs in Lisbon urban soils.

Risk Assessment of PAHs in Lisbon Soils
The first approach to assess potential problems concerning contamination of Lisbon soils with PAHs, consisted in comparing total levels with generic soil quality guidelines, i.e., guidelines not specifically related with human or ecosystems health.The Dutch guidelines refer three values for the sum of 10 PAHs: a target or background value (1500 µg kg −1 ), a maximum value for residential land use (6800 µg kg −1 ), and an intervention value (40,000 µg kg −1 ) [6].Considering these three levels, a prediction map was obtained which showed areas that are likely to exceed the Dutch guidelines (Figure 3).The former target value normally considered for PAHs (1000 µg kg −1 ; VROM [35]) is also shown in the prediction map, since it is still commonly used.
Since the Dutch guidelines were developed for a standard soil with 10% OM, the concentrations in Lisbon soils were corrected to this content.Observing this figure, it is possible to conclude that, with exception of the airport area, levels of PAHs can be considered safe for residential land use (<6800 µg kg −1 ).Yet, according to these guidelines, a considerable area is considered "slightly contaminated" (1500-6800 µg kg −1 ), and in this case, the possibilities for sustainable soil management should be determined [6].
Considering other generic guidelines, hazard quotients (HQs) were calculated.Regarding precaution levels from Germany [2,36], 11% of sites (samples 4, 5, 12, 23, 31, 53, 66, 69, 71, 96, and 97) exceeded this value for BAP and for the Sum16PAHs, indicating a certain chance of future soil problems.Compared to the Finish thresholds [2,36], the values were exceeded in 18% for BAP (samples 9, 11, 16, 17, 22, and 49, in addition to the ones referred for German guidelines), and less than 5% for the other compounds.According to the Finish guidelines, a site-specific risk assessment should be carried out in these areas.

Risk Assessment of PAHs in Lisbon Soils
The first approach to assess potential problems concerning contamination of Lisbon soils with PAHs, consisted in comparing total levels with generic soil quality guidelines, i.e., guidelines not specifically related with human or ecosystems health.The Dutch guidelines refer three values for the sum of 10 PAHs: a target or background value (1500 µg kg −1 ), a maximum value for residential land use (6800 µg kg −1 ), and an intervention value (40,000 µg kg −1 ) [6].Considering these three levels, a prediction map was obtained which showed areas that are likely to exceed the Dutch guidelines (Figure 3).The former target value normally considered for PAHs (1000 µg kg −1 ; VROM [35]) is also shown in the prediction map, since it is still commonly used.
Since the Dutch guidelines were developed for a standard soil with 10% OM, the concentrations in Lisbon soils were corrected to this content.Observing this figure, it is possible to conclude that, with exception of the airport area, levels of PAHs can be considered safe for residential land use (<6800 µg kg −1 ).Yet, according to these guidelines, a considerable area is considered "slightly contaminated" (1500-6800 µg kg −1 ), and in this case, the possibilities for sustainable soil management should be determined [6].
Considering other generic guidelines, hazard quotients (HQs) were calculated.Regarding precaution levels from Germany [2,36] value for BAP and for the Sum16PAHs, indicating a certain chance of future soil problems.Compared to the Finish thresholds [2,36], the values were exceeded in 18% for BAP (samples 9, 11, 16, 17, 22, and 49, in addition to the ones referred for German guidelines), and less than 5% for the other compounds.According to the Finish guidelines, a site-specific risk assessment should be carried out in these areas.
This example gives a clear demonstration how differences in the soil quality guidelines between countries result in very different efforts and costs when assessing potential risks.Very protective levels may result in higher costs, considering that a site-specific RA should be conducted in sites that exceed those levels.This further demonstrates the utmost importance of how each country defines its own guideline values, however, a common methodology to obtain these values should be implemented in order to prevent high discrepancies among countries.The prediction maps, as shown in Figure 3, can be very useful at this point in order to identify areas in which a site-specific RA should be conducted.Before starting a site-specific RA, two major issues must be addressed: land use and area dimensions.Land use will define receptors potentially at risk and exposure pathways, and in the case of urban areas either an ecological risk assessment (ERA) or a human health risk assessment (HHRA) may have to be performed.In addition, land uses can be subdivided into levels of sensitivity (e.g., land uses where there is an evident likelihood of children ingesting soil or areas with ecological importance can be considered more sensitive).The dimension of affected areas is an important issue in urban areas since, due to the fragmentation of soils, in some cases the areas affected may be very small and therefore the effects may be reduced.With this approach is it is possible to prioritize areas that may need special attention during urban planning.In addition, the number of sites to be further studied can be substantially reduced.

Potential Risks to the Environmental Health (ERA)
Potential risks to the environment were first calculated using HQs, considering different guidelines for protection of environmental health.Comparing levels of PAHs found in Lisbon soils with the quality criteria based on ecotoxicological effects recommended by Danish EPA guidelines [2,37], it was observed that 18% of samples exceeded the guideline values for Sum5PAHs and 28% for BAP.Considering the Spanish guidelines for protection of soil organisms [2,38], 14% of samples had a HQ > 1 for BAP and 8% for FLA.Regarding the Finish guidelines [2,36], values were exceeded for FLA, BAA, and SumPAHs, but in all cases less than 5% of samples can be considered a potential The prediction maps, as shown in Figure 3, can be very useful at this point in order to identify areas in which a site-specific RA should be conducted.Before starting a site-specific RA, two major issues must be addressed: land use and area dimensions.Land use will define receptors potentially at risk and exposure pathways, and in the case of urban areas either an ecological risk assessment (ERA) or a human health risk assessment (HHRA) may have to be performed.In addition, land uses can be subdivided into levels of sensitivity (e.g., land uses where there is an evident likelihood of children ingesting soil or areas with ecological importance can be considered more sensitive).The dimension of affected areas is an important issue in urban areas since, due to the fragmentation of soils, in some cases the areas affected may be very small and therefore the effects may be reduced.With this approach is it is possible to prioritize areas that may need special attention during urban planning.In addition, the number of sites to be further studied can be substantially reduced.

Potential Risks to the Environmental Health (ERA)
Potential risks to the environment were first calculated using HQs, considering different guidelines for protection of environmental health.Comparing levels of PAHs found in Lisbon soils with the quality criteria based on ecotoxicological effects recommended by Danish EPA guidelines [2,37], it was observed that 18% of samples exceeded the guideline values for Sum5PAHs and 28% for BAP.Considering the Spanish guidelines for protection of soil organisms [2,38], 14% of samples had a HQ > 1 for BAP and 8% for FLA.Regarding the Finish guidelines [2,36], values were exceeded for FLA, BAA, and SumPAHs, but in all cases less than 5% of samples can be considered a potential hazard.Similarly, considering the Canadian guidelines for environmental health protection [2,17], values were exceeded in less than 5% of samples for PYR, BAA, BBF, BKF, and IND.In the case of Finish and Spanish guidelines, values were corrected for the properties of the soils tested.
The second approach to identify potential risks to the environment consisted in the calculation of toxic units (TU).PNEC values used in these calculations were taken from the EU risk assessment report for coal-tar pitch [2,39], which were derived by applying an assessment factor to the lowest available non-observed effect concentration (NOEC), and the TU calculated for each individual compound can be found in Table 3.Since values were derived for a standard soil with 2% OC and 3.4% OM, they were first corrected for the properties of the soils tested.Risks were always below 1 for the LMW compounds (NP, ACY, ACE, FLU, and PHE), with the exception of ANT in one sample.For all other compounds, at least two samples showed a ratio above 1, and for BAA, BAP, IND, and BGHI, unacceptable effects on organisms were likely to occur in more than 10% of samples.BAP is the compound with higher percentage of samples exceeding the target value of 1, followed by BAA, IND, and BGHI.Considering the toxic units of the mixture (TUm), in 60% of samples, unacceptable effects on organisms were likely to occur as a result of the exposure to the mixture of PAHs (Table 3).Indeed, the median value of the TUm was 1.8 and it can reach values as high as 93.
Since using PNEC values did not provide a quantification of the expected impact, TU were also calculated using a hazard concentration (HC) value based on the species sensitivity distributions (SSD) modeling (Verbruggen, 2012).In this case, values corresponding to the HC5, i.e., to the 5th percentile of a SSD based on chronic no observed effect concentrations.The ratio between the concentration in sample and the HC5 value gives a quantification of the likelihood, and a characteristic of the extent of effects.Considering this, the results presented in Table 3 correspond to the percentage of samples for which individual PAHs are at concentrations that will affect a given percentage of species (more than 5%), and whose harms will be similar to the ecotoxicological data used to derive the SSD.
When the ratio is equal to 1, it means that the concentration of the compound in the sample is equal to the concentration that has been predicted to affect 5% of the species of the ecosystem.At this level, the number of species that is affected by the toxicant is equal to 5% and the risk on adverse effects is equal to or lower than 5%.For some samples this ratio was exceeded for individual compounds, and considering the TUm, the two highest values were 42 and 31 (observed in the airport samples), showing that both locals had a PAHs load that could clearly affect a percentage greater than 5% of the species of the soil ecosystem.All other values were below 20, but it should be noted that for almost 40% of samples the probability of a random species being affected by the ∑16PAHs will be greater than 5%.
Comparing results in Table 3 it is clear that the PNEC values are much more protective than the HC5 values, which is mostly based on the fact that the former are based on less information, thus they need to be more conservative.Moreover, the TU calculated based on the PNEC values tended to increase with the increase on ring number, whereas the TU calculated using the HC5 were higher for 4-ring compounds (FLA and PYR), which is probably related to the way these later values were derived, i.e., levels in soils were predicted based on the equilibrium partition theory and therefore only the pore water uptake was considered [40].
Table 3. Minimum, median, and maximum values of toxic units (TU) for individual compounds and toxic units of the mixture (TUm), using the predicted no effect concentration (PNEC) values and the hazard concentration (HC5) levels.Percentage of samples with TU and TUm above 1 is also shown.The multisubstance potentially affected fraction (mSPAF) of species was calculated and the toxicity threshold considered was the HC5 derived by the Dutch authorities, as referred above [40].The minimum value for mSPAF was 2.6%, the median 38%, and the maximum 98%. Figure 4 shows the spatial distribution of mSPAF values calculated using this most protective value, being possible to identify areas with higher risk.In the two darker areas of Figure 4, more than 50% of species are affected by the mixture of toxicants beyond the toxicity threshold considered.

Compound
Finally, the toxic pressures (TP) were also calculated.In this case, the effect level selected was the different soil screening levels (SSL), and TPs were calculated for individual compounds.Thus, the environmental quality guidelines from The Netherlands [40], Spain [38], and Canada [17], which were all obtained based on assessment factors, were used in these calculations.Considering the maximum permissible concentration (MPC) values from The Netherlands, which corresponds to concentrations that should protect all species in the ecosystem from adverse effects, the compounds with the highest median TPs were BAP and BAA (Table 4).With exception of NP, ACE, FLU, and PHE, at least one sample presented a TP above 50%.Considering the serious risk concentration (SRC), concentrations at which functions or species will be seriously affected, the TP of individual compounds were very low and the maximum value was 22% for BGHI in one of the airport samples.Using the Spanish guidelines for protection of soil organisms, the higher median value was observed for BAP, with 15 samples showing TPs higher than 50%, followed by FLA.Considering the Canadian guidelines for environmental health protection, TPs were very low, and only 3 samples showed TP values higher than 50% for PYR, BAA, BBF, BKF, and IND.
The Toxic Pressures of the mixture (TPm) were then calculated.Considering the Dutch MCP values, a very high median value (77%) was observed, with 16 of the 97 samples showing a TP of 100%.Using the SRC values, the median TPm was only 1.4%, with just one sample showing a value higher than 50% (Table 4 and Figure S4).Considering the other two guidelines, the median value was lower than 20% in both cases, where 17 samples presented a TPm higher than 50% for the Spanish guidelines and 16 for the Canadian (Table 4 and Figure S4).
only the pore water uptake was considered [40].
The multisubstance potentially affected fraction (mSPAF) of species was calculated and the toxicity threshold considered was the HC5 derived by the Dutch authorities, as referred above [40].The minimum value for mSPAF was 2.6%, the median 38%, and the maximum 98%. Figure 4 shows the spatial distribution of mSPAF values calculated using this most protective value, being possible to identify areas with higher risk.In the two darker areas of Figure 4, more than 50% of species are affected by the mixture of toxicants beyond the toxicity threshold considered.Finally, the toxic pressures (TP) were also calculated.In this case, the effect level selected was the different soil screening levels (SSL), and TPs were calculated for individual compounds.Thus, the environmental quality guidelines from The Netherlands [40], Spain [38], and Canada [17], which were all obtained based on assessment factors, were used in these calculations.Considering the maximum permissible concentration (MPC) values from The Netherlands, which corresponds to concentrations that should protect all species in the ecosystem from adverse effects, the compounds with the highest median TPs were BAP and BAA (Table 4).With exception of NP, ACE, FLU, and PHE, at least one sample presented a TP above 50%.Considering the serious risk concentration (SRC), concentrations at which functions or species will be seriously affected, the TP of individual compounds were very low and the maximum value was 22% for BGHI in one of the airport samples.Using the Spanish guidelines for protection of soil organisms, the higher median value was observed    All the calculations presented above correspond to a first tier of the risk assessment and a site-specific ecological RA should then be performed in areas identified as being of potential concern.The areas defined as potentially affected will be strongly dependent on the ecotoxicological endpoint and/or guidelines selected.There is a high variability of protection values (Figure S4), which may depend on the ecotoxicological endpoint considered (e.g., MPC vs. SRC values), but there is also a high variability among countries, which may be due to minor differences in the scientific methodology, or due to political and social assumptions considered during the derivation process.In order to obtain similar protection values within European countries, the possibility of implementation of a common ERA approach has been suggested [41], but not yet implemented.

Potential Risks to Human Health (HHRA)
Regarding the human health protection guidelines from several countries (Canada, Denmark, Italy, Spain, and The Netherlands), only less than 5% of samples showed levels above the recommended [2,17,[36][37][38].The exceptions were for the Sum5PAHs (15% considering the Danish soil quality criteria) and BAP (31% for Danish guidelines, 10% for Canadian, and 15% for Dutch).
In Figure 5a it is possible to observe the distribution of BAPeq concentrations, based on its statistical distribution: minimum (0.59 µg BAPeq kg −1 ), quartiles (23,74 and 203 µg BAPeq kg −1 ), the upper outlier limit (468 µg BAPeq kg −1 ), and the maximum value (7,653 µg BAPeq kg −1 ).It is interesting to note that the distribution map is very similar to the one of the distributions of the Sum16PAHs (Figure 1).The mean value for the BAPeq in Lisbon was 347 µg BAPeq kg −1 , and compared to other cities around the world, these levels were lower than the ones found in Shanghai (428 µg BAPeq kg −1 ), but higher than levels found in Beijing (181 µg BAPeq kg −1 ) or in Tarragona (124 µg BAPeq kg −1 ) [21,42,43].
Minerals 2019, 9 FOR PEER REVIEW 12 recommended [2,17,[36][37][38].The exceptions were for the Sum5PAHs (15% considering the Danish soil quality criteria) and BAP (31% for Danish guidelines, 10% for Canadian, and 15% for Dutch).In Figure 5a it is possible to observe the distribution of BAPeq concentrations, based on its statistical distribution: minimum (0.59 µg BAPeq kg -1 ), quartiles (23,74 and 203 µg BAPeq kg -1 ), the upper outlier limit (468 µg BAPeq kg -1 ), and the maximum value (7,653 µg BAPeq kg -1 ).It is interesting to note that the distribution map is very similar to the one of the distributions of the Sum16PAHs (Figure 1).The mean value for the BAPeq in Lisbon was 347 µg BAPeq kg -1 , and compared to other cities around the world, these levels were lower than the ones found in Shanghai (428 µg BAPeq kg -1 ), but higher than levels found in Beijing (181 µg BAPeq kg -1 ) or in Tarragona (124 µg BAPeq kg -1 ) [21,42,43].Although only the Canadian guidelines refer to the use of BAPeq, we decided to compare the distribution of the SumBAPeq with all the guidelines values of BAP.Therefore, Figure 5b shows the distribution of SumBAPeq in which the intervals chosen were based on the several guidelines available: 15 µg kg −1 was the USEPA soil screening level for residential areas; 100 µg kg −1 corresponded to the Italian value for residential areas and to the Danish soil quality criteria; 200 µg kg −1 corresponded to the Spanish value for residential areas; 600 µg kg −1 was the Canadian limit for protection of human health considering a target total lifetime cancer risk (TLCR) of 10 −6 ; 1000 µg kg −1 corresponded to the Danish cut-off criteria; 2000 µg kg −1 corresponded to the German value for playgrounds and to the Finish lower guideline value considering a target TLCR of 10 −5 ; 4000 µg kg −1 corresponded to the Although only the Canadian guidelines refer to the use of BAPeq, we decided to compare the distribution of the SumBAPeq with all the guidelines values of BAP.Therefore, Figure 5b shows the distribution of SumBAPeq in which the intervals chosen were based on the several guidelines available: 15 µg kg −1 was the USEPA soil screening level for residential areas; 100 µg kg −1 corresponded to the Italian value for residential areas and to the Danish soil quality criteria; 200 µg kg −1 corresponded to the Spanish value for residential areas; 600 µg kg −1 was the Canadian limit for protection of human health considering a target total lifetime cancer risk (TLCR) of 10 −6 ; 1000 µg kg −1 corresponded to the Danish cut-off criteria; 2000 µg kg −1 corresponded to the German value for playgrounds and to the Finish lower guideline value considering a target TLCR of 10 −5 ; 4000 µg kg −1 corresponded to the German value for residential areas; and 5300 µg kg −1 was the Canadian limit for protection of human health considering a target TLCR of 10 −5 [17,[36][37][38]44].No sample reached the value of 10,000 µg kg −1 , which is the maximum value for commercial areas, according to the Italian guidelines, and the German limit for parks and recreational areas.
In the case of human health guidelines for carcinogenic compounds, the major differences between countries are likely to be due to the target TLCR defined, normally ranging from 1 in 10,000 (10 −4 ) to 1 in 1,000,000 (10 −6 ).In addition, differences are also related to the land use scenarios considered, correction or not for soil type, pathways and conditions of exposure, inclusion or not of non-soil sources of exposure, or due to the toxicological data considered [2,36].Therefore, when choosing the guidelines to be used in this first step of RA, it is advised to first understand how they were derived.
Only in the white areas of the map (Figure 5b), the levels were below any of the guidelines and in the light-yellow area, values were above the most conservative value (from USEPA).In Figure 5b some of the playgrounds and schools sampled were also identified.By crossing this information, it is possible to identify areas in which a more detailed assessment should be made, since these are sensitive land uses regarding human exposure.
Regarding cancer and mutagenic risk assessment, the results for residential, worker, and recreational exposure to PAHs are shown in Table 5. Risks for residential exposure are high, being the reasonable maximum exposure, given by the upper confidence level of the mean (UCL 95%), higher than 10 −6 for cancer risk and higher than 10 −5 for mutagenic risks.In both cases, the values were higher than the target excess individual lifetime risk of 10 −6 advised by USEPA [18].Figure 6 shows the prediction map of TLCR for residential land use, being possible to identify areas which may represent a potential risk.The most contaminated areas in Lisbon represent a risk of 1.2 × 10 −4 , and are located nearby the airport and at lower extent in the southwest.In Figure 6 it is also possible to cross this information about probable risks with the population density of the Lisbon districts and conclude that in some of the most populated areas the risks are between 10 −6 and 10 −5 .If samples from the airport are excluded, the cancer risk, given by the UCL (95%), drops to 6.5 × 10 −6 and for mutagenic risks to 2.5 × 10 −5 , but it is still higher than the target value.
Figure 6 shows the prediction map of TLCR for residential land use, being possible to identify areas which may represent a potential risk.The most contaminated areas in Lisbon represent a risk of 1.2 × 10 −4 , and are located nearby the airport and at lower extent in the southwest.In Figure 6 it is also possible to cross this information about probable risks with the population density of the Lisbon districts and conclude that in some of the most populated areas the risks are between 10 −6 and 10 −5 .If samples from the airport are excluded, the cancer risk, given by the UCL (95%), drops to 6.5 × 10 −6 and for mutagenic risks to 2.5 × 10 −5 , but it is still higher than the target value.Cancer risks were also calculated for occupational land use (outside workers, e.g., gardeners, basic sanitation technicians, and farmers) and it was observed that PAHs levels in Lisbon may represent some concern (Table 5), being the UCL (95%) higher than 10 −6 .In recreational areas, cancer risks were also calculated, with the most contaminated site showing a value higher than 10 −6 , but the UCL (95%) was lower than the target value.However, mutagenic risks were higher, with the most contaminated site representing a risk higher than 10 −5 , whereas the UCL (95%) was higher than 10 −6 .

Conclusions
In urban areas, the assessment of soil quality and the characterization of potential risks is not an easy task, mainly due its heterogeneity and complexity, and the enhanced patchy nature of these areas.The collection of a high number of samples with a broad geographic coverage and analysis of the total contents of PAHs can be a good initial approach.
As shown in this study, geostatistical tools can be very helpful to obtain a characterization of the entire urban area or to identify sources of contamination.The prediction maps obtained allowed for a clear identification of areas with the highest levels of PAHs.It was also possible to build the maps of the factor scores obtained from the factor analysis, which show the areas most affected by telluric and anthropogenic inputs.It was clearly shown that the areas with the highest PAHs concentrations were the ones close to the airport and the river (city center).The latter area is the oldest part of the city and has a high population density.The high concentrations found in soils from the city center should be a result of long-term accumulation due to diffuse pollution mostly from traffic (through atmospheric emissions, tire debris and fuel exhaust, as well as pavement debris).Indeed, most of the sites sampled in the city center were historical gardens and parks.
The identification of contaminant sources, a critical step in risk assessment and management, in the case of PAHs, is normally performed using multivariate statistical methods or simple profiles of compounds and isomer ratios.Nevertheless, in the specific case of soils, the source apportionment based on individual compounds and isomer ratios are not free of some serious drawbacks since profiles are normally weathered and due to the presence of mixed sources.

Minerals 2019, 9 FOR PEER REVIEW 4 Figure 1 .
Figure 1.Locations of sampling sites and contour map interpolated by ordinary kriging of the distribution of the Sum16PAHs.The class limits correspond to the minimum, the quartiles (25, 50, and 75), the upper outlier limit, and the maximum value.Some of the most important sources are identified.

Figure 1 .
Figure 1.Locations of sampling sites and contour map interpolated by ordinary kriging of the distribution of the Sum16PAHs.The class limits correspond to the minimum, the quartiles (25, 50, and 75), the upper outlier limit, and the maximum value.Some of the most important sources are identified.

Figure 2 .
Figure 2. Contour maps (interpolated by ordinary kriging) for the first (a) and second (b) factor scores obtained in the factor analysis.

Figure 2 .
Figure 2. Contour maps (interpolated by ordinary kriging) for the first (a) and second (b) factor scores obtained in the factor analysis.

Figure 4 .
Figure 4. Contour map (interpolated by ordinary kriging) for the multisubstance potentially affected fraction (mSPAF) calculations using the concentration addition model and the maximum permissible concentration levels (HC5); the class limits correspond to the 25, 50, and 75% of potentially affected species.

Figure 4 .
Figure 4. Contour map (interpolated by ordinary kriging) for the multisubstance potentially affected fraction (mSPAF) calculations using the concentration addition model and the maximum permissible concentration levels (HC5); the class limits correspond to the 25, 50, and 75% of potentially affected species.

Table 4 .
Minimum, median, and maximum value of toxic pressures (TP) for individual compounds and the toxic pressures of the mixture (TPm) according to different guidelines: Dutch maximum permissible concentrations (MPC) and serious risk concentration (SRC), Spanish, and Canadian.

Figure 5 .
Figure 5. Contour map for SumBAPeq levels in Lisbon area, which was interpolated by ordinary kriging; the class limits correspond to the minimum, the quartiles (25, 50, and 75), the upper outlier limit, and the maximum value (a), and considering the different guideline values for human health protection (b).

Figure 5 .
Figure 5. Contour map for SumBAPeq levels in Lisbon area, which was interpolated by ordinary kriging; the class limits correspond to the minimum, the quartiles (25, 50, and 75), the upper outlier limit, and the maximum value (a), and considering the different guideline values for human health protection (b).

Figure 6 .
Figure 6.Contour map interpolated by ordinary kriging for the cancer risks considering a residential land use; the class limits correspond to the risk limits suggested by different guidelines.The map also shows the 53 districts of the city and their population density.

Figure 6 .
Figure 6.Contour map interpolated by ordinary kriging for the cancer risks considering a residential land use; the class limits correspond to the risk limits suggested by different guidelines.The map also shows the 53 districts of the city and their population density.

Table 1 .
Median, mean concentration, and range of polycyclic aromatic hydrocarbons (PAHs) (µg kg −1 ) in urban soils around the world, as well as the number of inhabitants (inhab) in each city.

Table 2 .
Factor loadings of PAHs and potentially toxic elements (PTEs), after Principal Component extraction and Varimax rotation.Factor analysis (FA), including organic carbon (OC) and clay, is also shown.Coefficients higher than 0.5 are marked in bold.Communalities, eigenvalues, and the cumulative explained variance (CEV in %) by the factors are also shown.

Table 5 .
Statistical data, including the upper confidence level (UCL) of the mean, of total cancer and mutagenic risks for different land uses of Lisbon soils.