*3.1. PCA and PCs Extraction*

PCA/FA allowed us to extract four rotated PCs accounting for about 86% of cumulative variance from the original 29 variables. Table 2 shows the matrix of the Factor loadings, which represent the correlation of each variable with each rotated PCs.


**Table 2.** Factor loadings matrix: the factor loading values higher than |0.5| are highlighted in bold and indicated as the significant loadings for the corresponding factor. Total variance explained by the rotated components (Fs) are shown.


**Table 2.** *Cont*.

Based on factor loadings, the Factors can be interpreted as follows:


PCA/FA highlights that PAHs concentrations are split into two different pollution components and this is consistent with their mobility characteristics. The molecular weight is in fact considered as an index of the dispersion of a pollutant: the heavier a substance is, in terms of molecular weight, the lower their environmental mobility is (Supplementary Materials: Table S1 shows the PAH molecular weights).

#### *3.2. Factors vs. Sewage Discharges*

After the factors identification, a bivariate correlation analysis with the distances from the existent discharges was performed. In Table 3, the Pearson's linear correlation coefficients are shown. It can be observed that the significant correlations are all negative, outlining that higher is the distance of the sampling point from the discharge point, the lower is the pollutant concentration. It is also worthwhile to observe that the points labeled as C, control points of the specific discharge points, show the same pattern of correlation of actual discharge points.

**Table 3.** Correlation coefficients between the vari-factors obtained by principal component analysis (PCA)/factor analysis (FA) and their distances from the sewage discharges. The correlations highlighted in bold and with a superscript are significant. The number of data are 125.


\* Correlation is significant at level 0.05 (two-tailed). \*\* Correlation is significant at level 0.01 (two-tailed).

The results show that the pollutants included in F1 are inversely correlated with the distance from the discharges 5 and 5C (i.e., the discharge and the control point of the former industrial area). No significant correlations were instead found between the pollutants loaded on F2 and the distances from the discharges, while F3, including Cr, Ni and V, correlates inversely with 8C, 9, 9C discharge points. These discharges are located on the southern side of the Nisida isthmus and are the discharges of the Coroglio Plant 8 and 9 (Figure 3). Finally, F4 showed significant correlations with distances from all discharges except 9 (and 9C). Being F4 loaded only by As, it may be speculated that As concentrations did not depend on specific discharge points, but it may rather derive from the whole coastal area and it is distributed along the whole gulf. The geology of the land surrounding the gulf (i.e., pyroclastic rocks) had in fact characteristics that were consistent with a geogenic origin of the As. Moreover, this analysis suggested that the presence of As in sediments is not necessarily due to the industrial spillage, but it may have a possible natural origin related to the geology of the area.

Based on these results, F2 was the only pollution component, which was not found correlated with any discharge point. This can be a consequence of the higher mobility of these compounds and their tendency to disperse.

#### *3.3. Factors vs. Thermal Springs*

As mentioned before, F4 (representative of As contamination pattern) resulted to be inversely correlated to the majority of the discharge points located along the Bagnoli coastline, suggesting that the As contamination in the gulf might not be due to a single source, but it may be caused by a presence of multiple sources (i.e., discharged into the gulf through groundwater). In addition, As might be released through submarine thermal springs that were clustered on the Pozzuoli Bay seabed. In order to investigate this possibility, another correlation analysis was performed. Being As loaded only on F4, the other three components were excluded from this analysis (Table 4).


**Table 4.** Correlations of As with respect to the position of natural thermal springs present on Bagnoli bay seabed. Significant correlations are highlighted in bold. The number of data are 125.

\* Correlation is significant at the 0.05 level (2-tailed). \*\* Correlation is significant at the 0.01 level (2-tailed).

Some positive correlations between As and the distance from Thermal Springs 10 and 11 were found, whereas inverse correlations were found with T5, T6, T12 and T13. Such an ambiguous pattern might be attributed to the fact that only some of those hydrothermal springs were fully active as was shown by recent studies [16,18,19]. As some of those springs are located in the open sea, the positive correlations might be the effect of the bathymetry of the bay [19], where the seabed sinks rapidly from the coastal area, determining the coastal accumulation of the sediments enriched with As.

It is also worth noting that these springs were in the proximity of discharge points (e.g., Thermal Springs 5 and 6 which are close to the "Conca d'Agnano", Discharge 2, or Thermal springs 12 and 13, which are close to Discharges 8 and 9, shown in Figure 3). All these drains are likely to be affected by groundwater and thermal water infiltrations which may possibly contain As [57].

#### *3.4. Pollution Patterns and Wave Hydrodynamics*

In order to assess whether the wave hydrodynamics have an influence on pollutant concentrations in sediments, the wave characteristics estimated at each sampling point were interpolated and pollutant vari-factors were tested against five classes of wave characteristic (i.e., wave height, peak period and mean direction) through the Kruskal Wallis test. Figure 5 shows the bathymetry implemented in the wave numerical model, the significant wave height and the energy period as results of the wave propagation of the 40-year averaged energy equivalent scenario.

**Figure 5.** (**a**) Bathymetry implemented in the wave propagation model; (**b**) resulting significant wave height in the bay after propagation of the 40-year averaged energy equivalent sea state; (**c**) the same of (**b**) but in terms of energy period.

Besides wave height, peak period and mean direction, the energy content and the energy flux were also calculated. Figure 6 showed the distribution of the wave energy and wave power at each sampling point.

**Figure 6.** Distribution of the waves' energy (**a**) and the waves' power (**b**) at each sampling point.

HCA results suggested to subdivide the wave hydrodynamics profiles of the Gulf into four classes of increasing hydrodynamics/energy (see Table 5). These classes have been tested as fixed factors in the Kruskal Wallis test design.


**Table 5.** Characteristics of the four wave hydrodynamics classes.

The different pollution components (i.e., F1, F2, F3 and F4) showed different results:

1. F1, loaded by the heavier PAH compounds and by some heavy metals and F4 loaded by arsenic, were both found significantly influenced by wave hydrodynamics (test results were respectively H: 12.9; df: 3, *P* < 0.01; and H: 51; df: 3, *P* < 0.01).

2. F2, loaded by the lighter PAHs, and F3, loaded by chromium, nickel and vanadium were not found influenced by the wave hydrodynamics (*P* > 0.05).

Moreover, pairwise comparisons where significance values have been adjusted by the Bonferroni correction for multiple tests allowed to assess that the F1 pollution component had the highest concentration values associated with classes of intermediate hydrodynamics (*P* < 0.05, see Figure 7a.) whereas F4 (and therefore arsenic) was found mostly associated with a low hydrodynamics (*P* < 0.05, see Figure 7b.).

**Figure 7.** Box plots showing (**a**) F1 concentration values (heavy PAHs) vs. wave hydrodynamics classes; (**b**) F4 concentration values (As) vs. wave hydrodynamics classes.

Therefore, the PAHs' F1 component appears to be affected by the wave-induced currents, being the concentration pattern of these pollutants in sediments dependent also on the shallow water hydrodynamic parameters in the bay. On the other hand, the wave hydrodynamics was not found to influence the PAHs' F2 component whose higher mobility had probably contributed to homogenize their contamination level in the bay. Likewise, F3 pollutants (i.e., chromium, nickel and vanadium) do not any the wave-induced pattern in the bay. Arsenic (F4) instead was found strongly associated to the areas closest to the coast and characterized by a low wave hydrodynamics.

#### **4. Discussion**

This study confirms that a multivariate statistical analysis (PCA/FA) approach can be extremely effective in assessing the apportionment of contaminant sources, even in a site with a complex geological characteristic and a long historical industrial development. With respect to the previous studies (e.g., [3,4,6,7,18]) which used similar methods to investigate the sediments or biota ([31,33,34]), this study considers a larger portion of the Gulf of Bagnoli, and it investigates the relationship of metals/PAHs with both the hydrodynamics of the Gulf and potential discharges/sources. Additionally, while studies in literature considered mostly the sum of total PAHs, in this study 18 PAH compounds are considered individually as their contamination pattern in the gulf. This aspect of the analysis allows us to determine that PAHs of higher molecular weight follow the pattern of metals such as Cd, Cu, Hg, Pb, Zn and Fe and are mainly located near the former discharges of the ILVA steelwork plant. PAHs of lower molecular weight are more dispersed in the gulf. Moreover, while heavier compounds seemed to be influenced by the wave climate of the bay, lighter compounds were apparently much less influenced by it.

Heavy metals such as chromium, nickel and vanadium were found more concentrated near the Nisida northwestern coast and in offshore waters (Figure 8) and did not appear to be strongly influenced by the wave hydrodynamics of the outer gulf. These findings were coherent with the findings of other studies ([4,6,58,59]) even though the measurements in this study derived mostly from drilling and did not include sediment cores.

**Figure 8.** Representation of factor scores inside the bay. The interpolation of factor scores was made with the IDW (inverse distance weighting) method.

Armiento et al. [4] within the same ABBaCo framework had analyzed both drilling (i.e., top layer sediments or TL) and sediment cores and found that that TL samples were characterized by a lower contamination than sediment cores. It is important to note that the correlation patterns found in TL samples reflected what was found in studies that analyzed historical contamination patterns through sediment cores.

The association of Cd, Cu, Hg, Pb, Zn and PAHs was described by Romano et al. through the analysis of sediment cores in a study [6] that determined the historical contamination pattern from industrial activity. The same study found that As, Cr, Ni and V were associated with deeper levels (209–299 cm) of the cores having with the highest percentages of clay fraction, suggesting a prevalent natural contribution for these elements. The only contaminant that is almost completely uncorrelated with the other pollutants is arsenic. This metalloid is significantly present in the bay showing a gradient from offshore to the inshore zones. Arsenic variability was found significantly correlated with both the distance from discharge sites located along the coastline and with some thermal springs present on the seabed. This pattern suggests that it is impossible to attribute the origin to a single source, even when the anthropogenic or the geogenic origin is slightly dominant.

Finally, our study confirms that arsenic is not correlated with other elements and highlights that the spatial pattern of its contamination might be dependent on both a land-driven origin (e.g., As-enriched groundwater, see [58]) and the presence of subaerial/submarine geothermal springs [59]).

The correlation of arsenic with the sewage discharge sites suggests a potential effect due to multiple anthropic activities: the former ILVA steelworks used arsenic in their production cycle but glass factories, located along the coastline in the past, may have contributed to the contamination due to their use of arsenic in glass production. Furthermore, the correlation with thermal springs might account for a "natural" arsenic source and conditions are not infrequent where discharges are characterized by a mixture of anthropogenic and geogenic sources. One example of this is the Conca d'Agnano discharge point which releases a mix of water from Agnano lake (i.e., thermal) and municipal wastewaters. Sewage discharge might also contain As-enriched groundwaters due to gas-water-rock interactions inside the aquifer [57,60]. The diffusion of arsenic in the bay was found to be influenced by the wave hydrodynamics, suggesting that arsenic dispersion in sediments might be attributed to various sources (anthropic, natural, or "mixed") and further diffused due to the particular marine dynamics in the Bagnoli Bay.

### **5. Conclusions**

The results of this study contribute to the reconstruction of the contamination history of the Bagnoli coastal zone, determined by the presence and past activity of industrial sites. Due to the high urbanization of the Neapolitan coastal zone it is impossible to find an unpolluted area-characterized by the same natural environmental framework of the impacted industrial site—to use as a control area. In the present study, a PCA/FA assessment was carried out in the Bagnoli Bay in order to demonstrate the utility of the robust statistical tool in a complex environmental/industrial scenario to investigate the concentrations' variability and their relationship with the locations of sewers, industrial sites, thermal springs, and the wave action inside the GoP. This study confirms the importance of performing statistically robust multidimensional analysis to support the source apportionment assessment of marine sediment contamination. PCA/FA, considering several parameters, allowed better discrimination among the many contamination components affecting the Bagnoli Gulf area and proved to be a very powerful tool in a complex environment with a mixture of effects due to anthropogenic and natural sources. The results of the analyses confirm that the main contamination source is anthropogenic activities (i.e., former steel plant and sewage discharges) but it also suggests the existence of multiple anthropogenic and geogenic sources of arsenic and other metals that might be originating from the volcanic rocks present in the Phlegraean area.

These findings suggest the need to define a "natural background level" (NBL) for the area for arsenic and other heavy metals to distinguish the natural from the anthropogenic component of the contamination. The source apportionment assessment presented here may help define such NBLs, and also facilitate decisions on the contamination control and the remediation management of the area, permitting public authorities to apply knowledge-based management actions.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/8/2181/s1, Table S1: The molecular weights of PAHs, split into to their principal components.

**Author Contributions:** Conceptualization and methodology: all authors contributed equally. Data interpretation, A.A., R.S., S.G.; resources, D.V., L.M., P.C., G.A.; writing—original draft preparation, S.G. and L.C., A.A.; writing—review and editing, all authors contributed equally. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Italian Ministry for Education, University and Research (MIUR) through the ABBaCo project, grant number C62F16000170001.

**Acknowledgments:** We thank "Anton Dohrn" Zoological Station, ISPRA (Institute for the Environmental Protection and Research), CNR (National Research Center), INGV (National Institute of Geophysics and Volcanology) and INVITALIA for providing the information about Bagnoli Bay and about available datasets PAHs and heavy metals concentrations. We thank also Università degli Studi Luigi Vanvitelli for providing the characterization of the Bagnoli Bay wave climate. Further, we gratefully acknowledge Karen Holmberg for having read the paper, giving useful suggestions and correcting the English version. We also thank two anonymous reviewers for useful suggestions, which helped to clarify and improve the paper.

**Conflicts of Interest:** All authors declare no conflicts of interest.
