Source Apportionment Assessment of Marine Sediment Contamination in a Post-Industrial Area (Bagnoli, Naples)

: The area of Bagnoli (Gulf of Naples, central Tyrrhenian Sea) has been heavily exposed to pollution for over a century due to the presence of industrial sites along its coastline. The aim of this study is to analyze contaminant concentrations (i.e., heavy metals and hydrocarbons) in seabed sediments through a statistical multivariate approach. Multivariate methods permit us to describe the pollution dynamics a ﬀ ecting the area and distinguish between anthropogenic and natural pollution sources. Additionally, the association between contamination patterns and the wave climate characteristics of the gulf (i.e., wave period, direction, height, power, and energy) is investigated. The study conﬁrms that the main contamination source in the Bagnoli bay is anthropogenic activities (i.e., former steel plant and sewage discharges) for the majority of investigated pollutants. It also provides evidence, however, for the potential co-existence of multiple anthropogenic and geogenic sources of arsenic and other metals that may be originating also from the water-rock interaction and submarine volcanic emissions in the Phlegraean area.


Introduction
The contamination of marine environments represents an increasing global concern because of the potential risks to both human health and along the coast heavily affecting the marine ecosystems. The Mediterranean Sea, due to reduced circulation and the presence of multiple industrial inputs along the coastline, is particularly vulnerable to environmental impacts and risks. Moreover, there is evidence that contamination may persist long after the end of industrial activities [1]. The historical industrial district and metallurgical production at the Bagnoli steel factory (ILVA), active for roughly a century, has exposed the marine sediments of the Gulf of Pozzuoli (GoP) to pollution by heavy metals and polycyclic aromatic hydrocarbons (PAHs). This area became a key site for twentieth-century Italian economic in sediments and to assess the related response on the biological assemblages that inhabit the seabed. Some correlation analyses were carried out between pollutant concentrations and sediment granulometry [31]. Other studies examined the distribution patterns of the meiofauna and the diversity and abundance of microorganisms inhabiting the sediments of GoP in relation to environmental variation and chemical pollution [8,33,34]. Most PAH inputs in the environment are linked to anthropogenic activity (e.g., wastes from industrialized and urbanized areas, off-shore petroleum hydrocarbons production or petroleum transportation) [35]. One of the main issues is the connection between pollutant concentrations and their possible source (natural or anthropic) combined with the influence of wave climate on their concentration patterns. Within this project, the present research aims to apply a robust statistical approach (PCA/FA) to demonstrate its practical application for assessing the main contamination sources, distinguishing among natural or anthropic/industrial contributions, and finding correlation between the contaminant concentrations and wave hydrodynamics in the area. The workflow framework of the study is shown in Figure 1.
Water 2020, 12, x FOR PEER REVIEW 3 of 20 that inhabit the seabed. Some correlation analyses were carried out between pollutant concentrations and sediment granulometry [31]. Other studies examined the distribution patterns of the meiofauna and the diversity and abundance of microorganisms inhabiting the sediments of GoP in relation to environmental variation and chemical pollution [8,33,34]. Most PAH inputs in the environment are linked to anthropogenic activity (e.g., wastes from industrialized and urbanized areas, off-shore petroleum hydrocarbons production or petroleum transportation) [35]. One of the main issues is the connection between pollutant concentrations and their possible source (natural or anthropic) combined with the influence of wave climate on their concentration patterns. Within this project, the present research aims to apply a robust statistical approach (PCA/FA) to demonstrate its practical application for assessing the main contamination sources, distinguishing among natural or anthropic/industrial contributions, and finding correlation between the contaminant concentrations and wave hydrodynamics in the area. The workflow framework of the study is shown in Figure 1.

Study Area and Geological Settings
The study area, Bagnoli, is located within the Bay of Pozzuoli, in the western sector of the urban territory of Naples. Close to the Bagnoli site, the already mentioned Phlegraean fields is a large "super volcano" situated to the west of Naples, Italy. Declared a regional park in 2003, the area of the caldera consists of 24 craters and volcanic edifices, most of them submerged under the bay with very high hydrothermal activity and effusive gaseous manifestations [19] (Figure 2). From a geological point of view, this large caldera is in a state of quiescence. The Phlegraean fields experienced extreme volcanic activity in the last 39,000 years. The two main events occurred 35,000 years ago (the 'Campanian ignimbrite eruption'), followed by caldera collapse, and 15,000 years ago (the 'Neapolitan Yellow Tuff', NYT eruption). The NYT eruption formed the central and eastern part of the Bagnoli-Fuorigrotta plain, which is constituted by a sequence of pyroclastic and volcanic material that thickens along the Agnano crater. Agnano volcanic products dominates the western part of the plain. The last eruptive event occurred in A.D. 1538 and gave rise to the Mt. Nuovo cone. From then on, only bradyseismic and hydrothermal activity (in the Solfatara and Pisciarelli area) are present in the Figure 1. Workflow scheme of the study. This work basically generates three main results boxes: the first and the second box include bivariate correlations analysis, respectively between PC/Fs and distances from sewage discharge and PC/Fs and distances from thermal spring, while the third box includes the hierarchical cluster analysis (HCA) and Kruskal Wallis test.

Study Area and Geological Settings
The study area, Bagnoli, is located within the Bay of Pozzuoli, in the western sector of the urban territory of Naples. Close to the Bagnoli site, the already mentioned Phlegraean fields is a large "super volcano" situated to the west of Naples, Italy. Declared a regional park in 2003, the area of the caldera consists of 24 craters and volcanic edifices, most of them submerged under the bay with very high hydrothermal activity and effusive gaseous manifestations [19] (Figure 2). From a geological point of view, this large caldera is in a state of quiescence. The Phlegraean fields experienced extreme volcanic activity in the last 39,000 years. The two main events occurred 35,000 years ago (the 'Campanian ignimbrite eruption'), followed by caldera collapse, and 15,000 years ago (the 'Neapolitan Yellow Tuff', NYT eruption). The NYT eruption formed the central and eastern part of the Bagnoli-Fuorigrotta plain, which is constituted by a sequence of pyroclastic and volcanic material that thickens along the Agnano crater. Agnano volcanic products dominates the western part of the plain. The last eruptive event occurred in A.D. 1538 and gave rise to the Mt. Nuovo cone. From then on, only bradyseismic and hydrothermal activity (in the Solfatara and Pisciarelli area) are present in the area. The ground deformation in this area was recognized as caused mainly by the actions of long periods of uplift, eruptive activity, and possible deep magma movements [36][37][38]. According to some studies [39], the Phlegraean fields activity has been predominantly explosive and characterized by an interaction between magma and water. There is evidence of both interaction of the magma with sea water as well as intra-caldera lake water and deep-seated aquifers in some Plinian volcanic events. High concentrations of metals and metalloids, such as As, Mn and Fe are found in groundwater and in the sea sediments near the main tectonic and hydrothermal activity [3,18,31,40,41]. area. The ground deformation in this area was recognized as caused mainly by the actions of long periods of uplift, eruptive activity, and possible deep magma movements [36][37][38]. According to some studies [39], the Phlegraean fields activity has been predominantly explosive and characterized by an interaction between magma and water. There is evidence of both interaction of the magma with sea water as well as intra-caldera lake water and deep-seated aquifers in some Plinian volcanic events. High concentrations of metals and metalloids, such as As, Mn and Fe are found in groundwater and in the sea sediments near the main tectonic and hydrothermal activity [3,18,31,40,41].

Figure 2.
Location of the study area, Bagnoli bay, Naples, Italy.
Since the early 1900's Bagnoli has been affected by the presence of various industries. The most important industrial activity was through ILVA steelworks, formerly Italsider. The ILVA steel plant was located along the Bagnoli coastline and was characterized by two long piers (still present now) which served as a berth for large ships that carried raw materials such as coal, iron ore and limestone to the steelworks. In the early 1960's, in order to enlarge the plant, the area between the two piers was filled with industrial waste coming from the steelworks. A new coastline was designed, giving more space for the industrial activities, but causing consistent damage to the marine ecosystems. Furthermore, the presence of nine discharge points along the coastline, combining both industrial and civil sewage, may have contributed to the deterioration of the quality of the environment. Some of these collectors are no longer in use since the steelworks stopped its production in 1990.
The Bagnoli Bay is also characterized by the presence of some submarine thermal springs, mainly located close to the former industrial area. Some offshore thermal springs are also present in the GoP.

Data Availability: Sampling, Sewage Discharge, and Wave Information
In this study, a total of 126 sampling points were considered as shown in Figure 3 (black dot points). These samples were part of an explorative campaign carried out in 2017 in the GoP by the Stazione Zoologica Anton Dohrn [42,43]. On the sea bottom, 94 sediment samples were collected using drilling (0-50 cm) and 32 sediment samples were collected by bucket (surface layer). In each point, the concentrations of heavy metals and PAHs were measured. A total of 11 heavy metals (i.e., Since the early 1900's Bagnoli has been affected by the presence of various industries. The most important industrial activity was through ILVA steelworks, formerly Italsider. The ILVA steel plant was located along the Bagnoli coastline and was characterized by two long piers (still present now) which served as a berth for large ships that carried raw materials such as coal, iron ore and limestone to the steelworks. In the early 1960's, in order to enlarge the plant, the area between the two piers was filled with industrial waste coming from the steelworks. A new coastline was designed, giving more space for the industrial activities, but causing consistent damage to the marine ecosystems. Furthermore, the presence of nine discharge points along the coastline, combining both industrial and civil sewage, may have contributed to the deterioration of the quality of the environment. Some of these collectors are no longer in use since the steelworks stopped its production in 1990.
The Bagnoli Bay is also characterized by the presence of some submarine thermal springs, mainly located close to the former industrial area. Some offshore thermal springs are also present in the GoP.

Data Availability: Sampling, Sewage Discharge, and Wave Information
In this study, a total of 126 sampling points were considered as shown in Figure 3 (black dot points). These samples were part of an explorative campaign carried out in 2017 in the GoP by the Stazione Zoologica Anton Dohrn [42,43]. On the sea bottom, 94 sediment samples were collected using drilling (0-50 cm) and 32 sediment samples were collected by bucket (surface layer). In each point, the concentrations of heavy metals and PAHs were measured. A total of 11 heavy metals (i.e., Al, As, Cd, Cr, Cu, Fe, Hg, Ni, Pb, V, and Zn) and 18 PAHs (i.e., naphthalene, anthracene, phenanthrene, acenaphthylene, acenaphthene, fluorene, fluoranthene, pyrene, benzo(a)anthracene, chrysene, benzo(b)fluoranthene, benzo(a)pyrene, benzo(k)fluoranthene, indeno(1,2,3-cd)pyrene, benzo(ghi)perylene, dibenz(ah)anthracene, benz(j)fluoranthene, benz(e)pyrene) indicated by Environmental Protection Agency (EPA) as important toxicological contaminants were considered. All the samples and the analytical determination derive from the series of studies conducted under the framework of the ABBaCo project [4,33]. Contaminants general statistics are reported in Table 1. A total of 14 sewage discharge points were mapped and selected in the study area (from 1 to 9C, in Figure 3). Some of these discharge points were supposed to be inactive (e.g., Points 4, 7 and 8). However, it is well recognized that due to poor maintenance and some illegal discharges in sewage disposal these points, especially when heavy rains occur, release water carrying waste and sewage to the coastline. For this reason, all the discharge points shown in the figure have been considered as active. It should also be noted that the discharge points labeled as "C" in Figure 3 are control points chosen as the nearest to discharge point (i.e., 9C is the control point of discharge point 9). Representing natural contaminant sources spread along the Bagnoli bay seabed, 13 thermal springs were identified [18] (Figure 4). In addition to the presence of sewage discharges and thermal springs, concentration patterns can be influenced by marine dynamics that remobilize sediments and accumulate them in places where there are steady-state conditions. In order to include the analysis, the main effects of the wave- Representing natural contaminant sources spread along the Bagnoli bay seabed, 13 thermal springs were identified [18] (Figure 4). Representing natural contaminant sources spread along the Bagnoli bay seabed, 13 thermal springs were identified [18] (Figure 4). In addition to the presence of sewage discharges and thermal springs, concentration patterns can be influenced by marine dynamics that remobilize sediments and accumulate them in places where there are steady-state conditions. In order to include the analysis, the main effects of the wave- In addition to the presence of sewage discharges and thermal springs, concentration patterns can be influenced by marine dynamics that remobilize sediments and accumulate them in places where there are steady-state conditions. In order to include the analysis, the main effects of the wave-induced hydrodynamics inside the GoP area in the analysis, we incorporated wave properties generated by a numerical model of the gulf.

Numerical Wave Modelling
The wave climate analysis was carried out using the hindcast data from the European Centre for Medium-Range Forecasts (ECMWF). Significant wave height (H s ), mean period (T m ), and mean direction (Θ m ) for the time period ranging from January 1979 to December 2018 were extracted from the wave model (WAM) of the ERA-interim archive (ERA: ECMWF Atmospheric Reanalysis), available for download online [44]. Both the offshore energetic patterns and the nearshore water conditions have been studied by means of the MIKE 21 SW coastal propagation model [45]. The model has been previously calibrated and validated applying a multi-collocation-based estimation approach as described in [46]. Such studies perform an optimization procedure for: 1.
The WAM offshore hindcast data, consisting in an amplification of each value of WAM dataset by means of an "enhancement factor". The enhancement factor has been obtained by comparison of WAM hindcast dataset at a point located offshore the study area and the time series obtained by transposition of the available data from an offshore wave buoy record. Then, the fictitious WAM time series has been used as input for the numerical model.

2.
The nearshores wave propagation model output, applying a two-step calibration strategy by comparison with measurements from acoustic Doppler current profilers and a set of innovative low-cost drifter-derived GPS-based wave buoys [47] located both inside and outside the GoP.
It is worth noting that spreading and dispersion studies and studies of water quality or ecological systems in marine areas are generally carried out by means of coastal hydrodynamic models coupled with sediment transport and particle tracking models. Generally, such modelling is computationally demanding. Therefore, in the perspective to undertake multi-year wave hindcast studies more quickly, a spectral modeling of wave propagation has been used. In order to assess the hypothesis that the bay's hydrodynamics is a factor of influence on pollution patterns, a single wave scenario able to represent the yearly average condition in the bay in terms of wave energy flux, had to be selected. Such a sea state, here defined as "energy equivalent", was obtained considering the wave power content of each wave from the whole dataset of 40-year wave records. Operatively, from the 6-h data of H s , T m , θ m provided by the fictious WAM model, a 6-h wave power and a wave energy dataset were obtained. For a specific sea state, the average wave energy per unit area is proportional to the square of the significant wave height, H s , according to the known relationship [48]: where ρ is the sea water density, assumed equal to 1025 kg/m 3 , g is the gravity acceleration (equal to 9.81 m/s 2 ). For a given spectrum, the significant wave height computation is based on zero-order moment of the spectral function and readily estimated as follows: and the wave characteristic/statistic periods can be defined as: where T e , T 01 are the energy period and the spectral mean period, respectively, with the spectral moment being defined as: where n = −2, −1, 0, 1, 2, . . . The wave power in irregular waves can be computed as: Following a conservative approach, according to [49,50], the energy period in the present study has been assumed as 1.14 T m . Equation (6) can be written in the approximate deep water expression: For each j-th year, an average wave power (P J ) has been computed. Considering that the energy period can be computed directly using the approximate formula: then, an "energetic" yearly significant wave height, H e , can be calculated as: for j = 1, . . . , 40, and a correspondent energy period can be estimated as follows: The 40-year average for these parameters gives H e = 0.93 m and T ee = 4.3 s. Regarding the representative direction, D e , a vectoral analysis about the energy flows (carried out comparing records of ADCP installed in the bay and correspondent offshore waves) indicate as the offshore wave direction can be considered 217 • N. Finally, the "energy equivalent scenario" with H e , T ee and D e has been propagated applying the MIKE 21 SW model. The seabed was performed by interpolating at the grid nodes the information provided by the General Bathymetric Chart of the Oceans (GEBCO) database [51]. A gross verification of the 10-m isobath, as representative of the nearshore region in a wave propagation model, has been applied by comparison with a recent bathymetric campaign [19].

Statistical Analysis: Multivariate Analysis and Bivariate Correlation
All the statistical computations were implemented using the statistical software IBM SPSS Statistics ® (version 26.0), while spatial queries, distance calculation and data mapping were managed through the open source environmental software QGIS (version 3.8.3, https://qgis.org/en/site/). PCA [52] is a powerful technique to investigate patterns of correlations among a set of variables. PCA in fact allows to extract a set of new variables, uncorrelated, called principal components (herein after PCs), each aggregating the variables more strongly correlated with the PC and among themselves. PCA seeks in fact a linear combination of input variables, which extracts the maximum variance of the bivariate correlation matrix. Following this, with a second linear combination, orthogonal from the first one, PCA extracts the maximum variance of the remaining variance, and so on. PCs represent all the linear combinations of the original variables weighted by their contribution in explaining the variance, in a specific orthogonal dimension. PCA was turned into a factor analysis (FA) to reduce the contribution of the less significant parameters within each component, by extracting a new set of vari-factors through rotating the axes defined by the PCA extraction. The Varimax rotation criterion was used to rotate the PCA axes allowing us to maintain the axes, orthogonality. The number of factors to be retained was chosen based on the "eigenvalue higher than 1" criterion (i.e., all the factors that explained less than the variance of one of the original variables were discarded).
To assess the variability within the whole set of contaminants, PCA/FA was performed considering all 11 metals, and all 18 PAHs (Table 1). For each of the N sampling points, PCA enabled to calculate principal component scores aggregating the information of the different parameters' concentration values. Moreover, in order to assess correlations between the sampled contaminants concentrations and the position of sewage discharges and thermal springs, a distance matrix was created through the QGIS vector toolbox, containing the pairwise distance between each sampling point and each discharge point. The measurements-discharges distance matrix was then associated to the component loadings and a bivariate Pearson's linear correlation analysis was performed. Finally, in order to assess the possible effects of the waves' dynamics on the concentration patterns found with PCA, a Kruskal Wallis test was performed to test the significance of hydrodynamics' effect on the contaminant concentrations, likewise similar studies [53][54][55][56]. Five wave characteristics were considered, such as the direction, the height, the mean peak period, the energy and the resulting power of the wave motion inside the bay. These characteristics were subdivided into classes, identified based on a hierarchical cluster analysis (HCA) [52] and tested as factors on the vari-factors assumed as dependent variable in the Kruskal Wallis test. Therefore, the KW null hypothesis was the equality of the vari-factor medians across wave characteristics classes. The level of significance assumed was P-value lower than 0.05.

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. F3 accounts for 6.8% of the variance and it is loaded by Cr, Ni and V. 4.
F4 accounts for 5.5% of the variance and it is loaded by As.
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).

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. 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.

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). 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].

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.
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].

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.  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.
Water 2020, 12, x FOR PEER REVIEW 13 of 20 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. 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. 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). 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. 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.).
Water 2020, 12, x FOR PEER REVIEW 14 of 20 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.). 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.

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 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.

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.
Water 2020, 12, x FOR PEER REVIEW 15 of 20 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. 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., Asenriched 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 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.

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.