This study presents the benefits of the application of multivariate techniques for the hazard assessment of a heavily polluted marine ecosystem. The study area, named Mar Piccolo, near Taranto city (Southern Italy), is a shallow marine basin located nearby an industrial compromised area, declared by the national government as Contaminated Site of Environmental Interest (SIN) due to the presence of long-lasting large industrial settlements that have severely impacted the marine environment. Besides the anthropogenic pressures, the marine basin is characterized by high productivity of several species at different trophic levels of the food chain, that confers to the bottom sediments an unusually high organic matter content. The latter is even enhanced by the presence of freshwater springs in the marine basin. The dynamism of the ecosystem demands for advanced evaluation tools for its correct characterisation. Multivariate ANOVA allied by Hierarchical Clustering are applied in this research to provide a readable picture of the quality status of the sediments, aiming at identifying different loading factors and getting insights into their simultaneous effects. The innovative approach adopted circumvents inefficient time-consuming procedures, usually required by the conventional univariate analyses of each parameter selected for sediment characterisation. A comprehensive hazard assessment was possible thanks to a clear graphical representation of the hazard distribution that supported the identification of the hazard controlling factors, confirming the efficiency of the adopted approach. The tools proposed herein can thus be recommended for the decision makers in investigating and interpretation of the quality status of complex polluted marine eco-systems.
Mar Piccolo is a shallow marine basin constituted by two inlets (Figure 1), near the city of Taranto (Ionian Sea, Southern Italy). It represents a very complex marine ecosystem, both from ecological and social–economic points of view. Water, sediments, and biota of this basin have been affected by a wide spectrum of anthropogenic pressures (among others, the most important steel factory in Europe, the largest Italian Navy shipyard, an oil refinery, and shipbuilding activities) (Cardellicchio et al., 2016). These stress factors have been responsible for a severe environmental contamination, mainly due to heavy metals and persistent organic pollutants with their consequent possible transfer to the aquatic trophic chain, also considering the widespread fishing and mussel farming activities currently operating in this area (Petronio et al., 2012, Giandomenico et al., 2016).
Besides the anthropogenic factors, the peculiarity of the area is enhanced by the presence of different submarine springs, discharging fresh water into the basin, that influence the equilibrium of the marine eco-system, conferring to it the typical features of dynamic transition environment (Zuffianò et al., 2016). In addition, in the first inlet of the Mar Piccolo, which is the most subjected to anthropogenic pressures, are present the most important submarine fresh water that introduce terrigenous material that makes this system very complex. The first inlet is also connected to the open sea through two channels, a natural one named Porta Napoli Channel, and an artificial one named the Navigable Channel, created 130 year ago to allow the transit of warships heading to the Italian Navy naval base, located in the southern part of the inlet. All these fluxes and connections promote a limited water circulation within the basin that, nevertheless, endorses a transport flux of contaminants associated to surface fine-sediment as well as in the very deep-water column strata. Such fluxes govern the contaminant plume migration, affecting the pathways for receptor exposure and control the hazard distribution and the hot spot areas. Bellucci et al., 2016 have identified an advective transport from the First to the Second Basin, considered as the leading transfer mechanism of fine particles and adsorbed contaminants. Considering that organic matter degradation is a dynamic engine behind benthic dynamics, deciphering the influence and the impact of these factors becomes essential in determining the hazard evolution in the whole marine ecosystem (Arndt et al., 2013).
The strong dynamism of the system that seems to be subjected to factors loading simultaneously, requires advanced characterisation tools that go beyond the conventional hazard degree distribution. Although the analysis of the contaminant concentrations provides important information about the hazard degree, it results inefficient to describe the evolution process of the pollution degree, hindering the identification of the factors controlling the distribution of contaminants, which is the most important information for undertaking correct solutions for environmental remediation.
The application of multivariate statistical methods in environmental data exploration is considered a great opportunity, since it supports a quick assessment to obtain even quicker responses (Mali et al., 2017). The chemometric models offer a useful way to understand many features of the analysed ecosystem, providing: i) similarities and dissimilarities between sampling sites; ii) detection of site outliers in the monitoring net and their explanation; iii) detection of hindered factors, which could be interpreted as anthropogenic or natural sources responsible for the chemical content of the environmental samples, iv) trends in the behaviour of pollutants; v) distribution of the contribution of each identified source to the total species mass. The aim of this study is to identify and distinguish anthropogenic factors influencing the contaminant transport and the relative hazard degree by using multivariate techniques. To this purpose, two multivariate tools were exploited: Hierarchical Cluster Analyses (CA) and Multivariate ANOVA (Analyses of Variance). The long-term aim of the research is to provide reliable protocols/tools for decision makers and economic, health and political solutions.
Materials and methods
Sediment sampling and pre-treatment
The data utilized for this analysis are referred to the investigation campaign conducted by ISPRA (Institute for Environmental Protection and Research) (ISPRA 2010) on 209 sites located in both inlets for a total of 1021 sediment samples. The sediment samples were collected in the period May 2010 – September 2011, according to the sampling strategy reported in Table S1. Core-samples were collected through vibro-corer PF1, equipped with a liner and a support vessel and provided of the differential GPS system for positioning of sampling cores. The sample cores were up to 1.50 - 3.50 m deep. Each core was divided into several 50 cm length sub-cores. Aliquots of wet-sediment samples were transferred from the liner into cleaned plastic bags and stored at 4 °C during transportation to laboratory, where they were stored at –20 °C before analysis. After that, samples were dried at 40 °C for 48 h. Each sediment was classified according to Shepard (1954) into four classes: gravel: > 2 mm; sand: 2-0.063 mm; silt: 0.063-0.002 mm; clay < 0.002 mm. A set of ASTM sieves was used for the granulometric separation. For the organic compound determination, wet samples were used. Foreign bodies (shells, wood) and any coarse fraction were removed prior analysis. Samples without sieving process have been used for determination of hydrocarbons with C < 12 extracted through mini-corer (syringes), according methods in Table S2.
Chemical and physical-chemical analyses
In this study, the following specific variables were considered: i) Heavy Metals and metalloids (As, Cd, Cu, Cr, Hg, Ni, Pb, V, Zn, Al and Fe); ii) Persistent organic pollutants in terms of sum of Polychlorinated Biphenyls (PCBs), sum of Polycyclic Aromatic Hydrocarbons (PAHs) and sum of Total Hydrocarbons (THC), the latter divided into two sub-groups: sum of Light molecular Total Hydrocarbons LTHC (VOCs and gasoline range organics, the latter including hydrocarbons from C5H12 to C12H26) and sum of High molecular Total Hydrocarbons, HTHC (hydrocarbon compounds between C12H26 and C40H82).
In addition, some physical-chemical parameters are also included. Grain size content (in terms of coarse%, sand%, mud% and clay%) and Total Organic Carbon (TOC). For the 209 sampling points, sediment cores were collected at six different depths (0.10 cm; 10-30 cm; 30-50 cm; 100-120 cm; 180-200 cm; 280-300 cm), and grouped in three classes: surficial sediment for those collected up to 50 cm depth; intermedium sediment layer for those collected in the range of 50-120 cm; and depth layer strata for samples collected at more than 120 cm depth. The whole analytical procedures for the selected variables are reported in Table S2.
Multivariate analyses have been widely exploited for generating spatial pollution maps and for identifying anthropogenic stressors in contaminated area (Mali et al., 2019; 2017 a, b, c). Two kinds of software were utilized for statistical analysis: Soft Independent Modelling of Class Analogy (SIMCA) 10.2 for Hierarchical Cluster Analyses (CA) and STATISTICA 10.0 for two-way factorial Analyses of Variance (ANOVA) and for normality-test of the raw and log-transformed data. Two-way ANOVA was preferred with respect to the three ANOVA in order to reduce the risk of getting false positive effects with increasing the independent factors (Mc Killup, 2012). The Type II ANOVA model was chosen because the levels assigned to main independent factors were randomly selected. In addition, the model was considered suitable for the purpose of the work.
CA is a well-known and widely classification approach used for environmental purposes with its hierarchical and non-hierarchical algorithms (Massart and Kauffman, 1983; Einax J. et al., 1997) It is widely recognized that CA allows to organize the whole set of variables into two or more mutually exclusive unknown groups/clusters based on combination of internal variables. A preliminary step of data scaling was performed (logarithmic transformation), where normalized dimensionless numbers replaced the real data values. The similarity between the objects in the space variable are determined on the base of Euclidean distance.
In this paper, the selected dependent variables are the concentration of each contaminant considered (As, Cd, Cu, Ni, Pb, Zn, Ni, V, Fe and Al, HTHC, TOC within sediments), while the categorical independent factors (Table S1) considered are the following:
organic matter content, expressed in terms of TOC concentration, considering four classes/levels: “low”, “medium”, “high” and “very high” for TOC concentration ranging from 0-1%w, 1-3%w, 3-5%w and > 5%w, respectively:
the grain-size features, considering three classes: sediments with predominance of mud fractions (ϕ < 63 micron), sediment with predominance of coarse fractions (ϕ > 63 micron) and mixed granulometric features (ϕ mixed grain size fraction).
Results and discussion
The data obtained by chemical analyses showed different levels of contamination within inner area of the first basins due to different pressures. Table S3 reports statistical data of contaminant concentration. In order to graphically present the spatial contamination pattern within the study area, we exploited the two-way ANOVA, Analyses of Variance of the contaminant concentrations considering two categorical spatial variables, the “depth” and “basin” type. The “depth”, as previously reported is referred to the depth core sampling and include three classes: “surficial” (S) referred to samples collected at 0-10cm, 20-30cm and 30-50 cm depth, “intermediate” (I) including samples depth of 100-120 cm depth, and “deep” (P) including samples depth of 180-200 cm and 280-300 cm. The “Basin” is referred to the areas mostly influenced by the main industrial activities and human pressure loading in the study area. We define six “basins” categories, as following:
“channel” including samples located at the entrance of the channels; “citro” indicates samples located around two main fresh-springs within the first inlet, as well as the samples in the Galeso river flow; “Ars” refers to samples located nearby the main impacted zone where the Navy is situated (Bellucci et al., 2016); “Mit” belongs to all core-samples located in the mussel-culture activity. Being some samples collected on the second inlet, we created also a group namely “2nd Sen” including all samples located within this inlet. Lastly, the samples countersigned with “In” indicates all remaining samples not included in previous groups.
Graphs of two-way ANOVA are considered meaningful for representing the spatial contamination (Figure 2). The left ordinate reports the average concentrations of contaminants (log-transformed values) considering different basin areas (in abscissa) and different depth classes (in right ordinate).
The analysis reveals that As, Hg, Pb, Cu, Cd and Zn display higher concentrations in surficial layers (blu lines) with respect to the deepest layers in almost all basins (Fig. 2a displays the trends for Cr, taken as an example. The graphs for the other elements are reported in Figure S1).
Hg, Pb and As concentrations have high values within Navy area, where high levels of Ptot are also registered. High contamination levels are also observed in some deepest layers (I and P) as revealed by large vertical bars (at 0.95 confidence intervals). Zn and Cu show the highest concentrations in the area of Mussels-culture area (Mit), while area around fresh-springs, “Citri”, display high contamination by Cu, Zn, Pb, Cd and organic pollutants (PAHs and HTHC). A peculiar trend is displayed by samples collected at “Channels” basin, that resulted to be strongly contaminated by HTHC, Pb and Cd. These samples are also characterized by very high concentration of TOC and nutrients (in terms of Ntot). It seems that fuel losses or oil spilling come from channels within basin. The contaminants of most concern resulted to be As, Pb, Hg and Cu that display high concentration not only in superficial layers but also in almost all intermedium levels (I) (red lines) of the basins considered.
As to the lithogenic metals (Al, Cr, V, Fe and Ni) distributions, similar average concentrations within all basins are observed, irrespective to depth. In almost all cases, the distribution of these metals resulted associated to granulometric distribution. The higher the grain size, the lower the metal concentration (except for Fe, which resulted to be irregularly high in the Navy irrespective to grain size). This circumstance suggests that grain size controls lithogenic metals variability: thus, normalisation of their concentration to grain size is needed in order to evaluate their real enrichment (Covelli and Fontolan, 1997). Indeed, it is supposed that the presence of channels nourishes the inner basins with sediments from the open sea and favours the introduction of terrigenous allochthone organic matter. It is well known that both phyllosilicates and organic matter have a strong chemical affinity for trace elements and organic pollutants (Mali et al., 2019), that are usually concentrated in fine (2–60 μm) fractions.
In order to identify similarities among contaminants, Hiearchical Cluser Analysis (CA) was also performed and the results were represented in the dendrogram shown in Figure 3. The data before CA are transformed and scaled properly. By comparing results reported in clusters, it is apparent that CA grouped contaminants in two clusters, the first one being constituted by Al, Fe, Cr, V and Ni, that are all lithogenic elements, and the second one including toxic metals, such as Hg, As, Cd, Pb, Cu and Zn. The first group resulted associated, as mentioned, with finest sediment (mud), while the second one, constituted by sulphophilic elements (those elements displaying high affinity for sulphur), resulted to be associated with TOC, Ntot and Ptot. The tight relationships between Pb, Hg, and organic markers in the cores indicate that these metals are derived from industrial combustion emissions. This correlation suggest the potential occurrence of methylation reactions promoted by organic matter and nutrients, and consequently potential leaching of metal contaminants from sediment to water column, especially for Hg and As that are known to be more subjected to methylations processes with respect to other sulphophilic elements (Mali et al., 2017c; Lofu et al., 2016).
The role of organic matter in chemistry of marine sediments is worldwide accepted. Thus, we decided to analyse the trend of Organic Matter (OM) content in the sediments. The OM trend is usually expressed in terms of TOC - Total Organic Carbon, although the OM in sediments includes also significant amounts of hydrogen, oxygen, nitrogen and sulphur. High concentrations of TOC are found within surficial and intermedium depth of “Citro”, “Channels” and “Mit” basins. The TOC concentration reached values of 8% in the first centimetres, resulting relatively higher than those found for the Adriatic and Ionic sea-basin, measured in sites away from the coast (organic carbon less than 2%) (Bartolini et al., 2016). This circumstance suggests the presence of additional allochthone source of organic carbon in these sites. In shallow marine basin, near the coast, terrigenous contribution of organic matter (allochthone) might occur. The allochthone OM, added to the organic matter naturally present in the marine sediments (autochthonous), makes more complex the processes that control the distribution of contaminants (Fostner, 1994; Mayer, 1994). In addition, it is known that Ntot plays an important role as a source of nutrients decreasing with depth due to the remineralization of organic matter and non-biological oxidation (Gao et al., 2012).
The ANOVA elaborations using spatial categorical variables (basin type and depth core) and the Cluster Analyses allowed to understand the influence of several factors on the contamination pattern and to distinguish groups of contaminants that characterize each basin and each depth layer.
However, when the hazard assessment is required, it is necessary to understand if simultaneous effects occur and at what extent these synergistic or antagonist effects control the hazard distribution pattern. For this reason, we used two-way ANOVA aiming to determine the factors that most significantly control the chemical compositions of marine sediments. The independent factors selected for ANOVA analysis and discussed in the present paper are: organic matter content, expressed in terms of TOC range concentration (TOC level) and grain size distribution, whose classes are reported in Table S1. We performed the two-way ANOVA considering pair of independent factors (A, B, AxB), i.e. TOC Level and Grain Size. The homogeneity of variance was tested using Leven test prior to analyses and data not respecting the required assumption were excluded. The results of two-way ANOVA (4x3) considering TOC as factor A and grain size as factor B, are summarized in Table S4, reporting the values of variance ratios (F), probabilities (p), and effect size (η2). The obtained values F(144, 5770.3)=3.6603, p = 0.0000 and η2 = 0.45 indicated that the contamination pattern resulted significantly controlled by the combined effects of TOC and Grain Size deducted either by the partial η2 determined, and graphically (presence of numerous intersections). The two way ANOVA confirms the presence of two groups of metals, the first one constituted by Pb, Hg, Cd, Cu, Zn, As (Figure 4 and Figure 5) and the second one constituted by Cr, Al, V, Fe and Al (Figure 6). Some slight differences are noticed also for distribution of contaminants of the first metal-group: indeed the mean concentrations of Pb and Hg in the samples with mixed granulometric distribution (red lines) and with high and medium TOC levels resulted higher with respect to the mean concentrations determined within the other two grain size fractions (θ<63 μm and θ>63 μm). On the other side, the concentration of As, Cd, Zn and Cu resulted always higher in the finest sediments than in the other two fractions (mix and coarse fractions), except for samples having medium TOC levels, for which it resulted higher than the one found in high and low TOC content sediments.
In Figure 6 the variability of concentration of two representative elements (Cr and Al) for lithogenic cluster metals is reported. As displayed graphically, all data confirm the multiple effects of multiple factors in the spatial distribution of contaminants within Mar Piccolo (and the complexity of the whole ecosystem). The organic matter content and the sediment grain size features affect the hazard distribution within the basin and control the location of the hot spot areas.
In this paper, the spatial distribution of multi-elements within a contaminated sea basin is discussed using two complementary chemometric techniques, the Hierarchical Cluster Analysis and the factorial ANOVA. The obtained results reflected a complex pattern that arises from the superposition of natural and anthropogenic processes and by multiple factors acting simultaneously on a local scale. The main processes and factors that control the contamination pattern in the Mar Piccolo resulted to be: i) the original composition/mineralogy of the sediment parent material; ii)the nature and origin of organic matter within different depth layers and the large scale of anthropogenic enrichment of elements by multiple sources; iii) The relative role of these factors and processes on the sediment composition is assessed through coupling of multifactorial ANOVA.
Looking at the spatial pattern of the accumulation of heavy metals, As, Hg and Pb concentrations appear most clearly related to hydrocarbons and fuel losses, while Cd, Cu, Zn seem to derive mostly from atmospheric inputs and contaminated material uncontrolled dumped in the basin. Moreover, for evaluation of future scenarios and associated environmental risks in the study area, it would, be important to assess the bioavailability of these elements. These different patterns show the importance of a spatially representative overview to infer the historical anthropogenic sources of these elements.
The present paper proves that chemometric methods are useful tools for assessing and modeling contamination patters on a large scale of polluted areas, contributing to efficiently and economically friendly monitoring their quality status.
ISPRA is gratefully acknowledged for chemical analyses and the availability of dataset.
Matilda Mali http://orcid.org/0000-0002-2523-2807
Leonardo Damiani http://orcid.org/0000-0002-3125-3284