Application of Geostatistical Tools to the Geochemical Characterization of the Peloritani Mts (Sicily, Italy) Aquifers

Abstract

Sources of groundwater contaminants in inhabited areas, located in complex geo-tectonic contexts, are often deeply interlocked, thus, making the discrimination between anthropic and natural origins difficult. In this study, we investigate the Peloritani Mountain aquifers (Sicily, Italy), using the combination of probability plots with concentration contour maps to retrieve an overall view of the groundwater geo-chemistry with a special focus on the flux of heavy metals. In particular, we present a methodology for integrating spatial data with very different levels of precision, acquired before and during the “geomatic era”. Our results depict a complex geochemical layout driven by a geo-puzzle of rocks with very different lithological natures, hydraulically connected by a dense tectonic network that is also responsible for the mixing of deep hydrothermal fluids with the meteoric recharge. Moreover, a double source, geogenic or anthropogenic, was individuated for many chemicals delivered to groundwater bodies. The concentration contour maps, based on the different data groups identified by the probability plots, fit the coherency and congruency criteria with the distribution of both rock matrices and anthropogenic sources for chemicals, indicating the success of our geostatistical approach.

1. Introduction

Access to safe drinking-water, sanitation and hygiene is fundamental to human health and well-being as well as to livelihoods, school attendance and dignity, helping to create resilient communities living in healthy environments [1]. The global water demand has increased as the population has grown six-fold during the past 100 years [2] as a consequence of the economic progress and of the change of the main water usages [3].

It is estimated that water demand will continue to increase particularly in economically emerging countries [3], both for industrial (20%) and domestic (10%) practises, even if agriculture remains the largest consumer (70%) [3,4]. Groundwater is one of the most important sources for drinking water: its consumption covers over a third of the global demand [5]; however, in arid and semi-arid regions, characterized by limited precipitation and high evapotranspiration, it is often the most important water resource [6].

Water quality is worsening due to the impact of chemicals and pollutants, especially in lower-middle income countries not only for the growth of population and economy but also for the lack of wastewater management systems [3]. The great majority of evident water-related health problems are the result of microbial contamination. However, a significant number of severe health concerns may occur as a consequence of chemical contamination [1] whose sources are anthropogenic, especially in countries experiencing a rapid economic development, and geogenic, due to the dissolution of minerals composing the Earth’s crust ([7] and the references therein). From the topological point of view, sources of groundwater contaminants can be distinguished between punctual, as sewage discharges or industrial plants, or sparse, as the leaching of fertilizers and pesticides used in intensive agriculture; punctual sources are easy to identify, whereas sparse sources are much more difficult to monitor [8].

In densely inhabited areas, located in complex geo-tectonic contexts, these sources are deeply interlocked, which makes their discrimination difficult, even when employing efficient geostatistical tools ([9], and the references therein; as a consequence of this, remediation actions are often difficult [10,11].

Peloritani Mts, roughly corresponding to North-Easternmost Sicily (Italy) (Figure 1), are an ideal place for these studies due to their very complex orographic and geo-tectonic setting, where a densely populated area is deeply exploited for industrial and agricultural purposes.

Dongarrà et al. [13] investigated groundwater from the south-eastern Peloritani area to recognise the processes influencing their composition. They found a general geochemical fingerprint of these waters dominated by the carbonate dissolution and the hydrolysis of aluminosilicates. High As, Pb, Sb and Zn concentrations, often exceeding the recommended limits for human consumption, were found and attributed to the weathering of ore deposits associated to metamorphic lithotypes.

Cangemi et al. [11] applied geostatistical tools for discriminating between geogenic and anthropogenic sources of heavy metals in the groundwater of the Barcellona–Milazzo plain, located in the central-northern sector of the Peloritani Mts. Their results confirmed fluxes of heavy metals from mineralized ore deposits to groundwater and also highlighted the concurrence of anthropogenic heavy metal sources from intense agriculture and industrial activities.

Apollaro et al. [14] applied the concept of the Natural Background Level (NBL) for discriminating between natural and anthropic sources of As in groundwater bodies of the Calabria region (southern Italy), whose southern sector has the same geological nature of the Peloritani Mts. They combined aquifer-based preselection criteria with multivariate, non-parametric geostatistical tools and the Probability Kriging method, for mapping the probability of exceedance of As with respect to NBL, with the aim of supporting contamination risk mitigation strategies.

The present study has the aim of combining the information from these studies with previously unreleased data, applying homogeneous geostatistical criteria for their analysis and interpretation and retrieving an overall view on the groundwater geochemistry of the Peloritani Mts, with a special focus on the flux of heavy metals. We present and discuss a methodology for integrating data characterized by very different spatial precisions, which is of general interest for applications requiring the merging of new information with “pre-geomatic” datasets. Another aspect of general interest of this study is the application of geostatistical criteria, based on probability plots [15,16], for fixing “objective” thresholds separating geochemical anomalies from the background environmental noise.

2. Study Area Setting

The Peloritani Chain (PC) is a continuous mountain range that entirely occupies the north-easternmost termination of Sicily, reaching its maximum elevation at Montagna Grande (1374 m a.s.l.). On the Ionian side, it reaches directly the coastline, whereas, on the Tyrrhenian side, it decreases in the Barcellona–Milazzo coastal plane. This orographic asymmetry results in the distribution of population and related productive activities. The main town is Messina, the capital city of the homonymous province, with about 235,000 inhabitants, located on the Ionian side immediately south-eastward of the northern termination of Sicily.

The second most populated area is the coastal Tyrrhenian area, (about 191,000 inhabitants), which also hosts industries, petroleum refineries and it is exploited by intensive agriculture making use of pesticides and fertilizers. A different situation characterizes the coast of the Ionian side, with smaller cities populated by a total of 56,000 inhabitants. The inner, mountainous territory is scarcely populated with small cities and villages accounting for an overall population of about 88,000 inhabitants.

The PC represents the geometrically highest element of the Sicilian Fold and Thrust Belt [17,18]; it crosses the Ionian and the Southern Tyrrhenian Seas, limited by the so-called Drepano Thrust Front [19,20], continuing in North Africa as the Kabilyan chain [21]. It is the south-westward termination of the Calabro–Peloritan Arc (CPA), linking the NW-SE trending Apennines with the E-W trending Siculo–Maghrebid chain [22,23].

It consists of a succession of S-SE verging tectonic units (nappes) originating from the delamination of part of the original south-eastern European margin and/or of the Austroalpine–Alpine domain [19,23,24]. Deep offshore and onshore seismic profiles show the tectonic superimposition of the PC on the sedimentary units belonging to the original Sicilide basin, currently outcropping in the Nebrodi Mountains, which delimit westward along the PC [12,25].

The PC constitutes a crystalline Paleozoic basement, covered by Meso-Cenozoic sedimentary units, piled up since the Upper Oligocene–Early Miocene [12,26].

The terrains outcropping in the studied area, illustrated in the map of Figure 1, from the top to the bottom of the tectono-stratigraphic column are [12,26,27]:

• A Neogene-Quaternary Cover, topped by alluvial deposits with thickness up to several tens of metres (Middle Pleistocene–Holocene), marine terrigenous deposits, calcarenites and clays (Late Pliocene–Middle Pleistocene), clayey-sandy-calcarenitic deposits (Pliocene–Pleistocene), sands, conglomerates and sandy clays (Serravallian–Early Messinian);
• The Calabride Chain Units, belonging to the European margin and thrusted onto the Alpine Tethys complex during the Balearic stage and the counter-clockwise rotation of the Sardinia Block (Late Eocene–Early Miocene). Locally, The Calabride Chain Units consist of the Floresta calcarenites and M. Pilo marls (Late Burdigalian–Early Serravallian), Antisicilide Unit composed of varicoloured clays with olistolithes of calcarenites and sandstones (Late Cretaceous), Capo D’Orlando flysch (Late Oligocene–Early Burdigalian), Pre-Variscan and Variscan metamorphic basement overthrusted with its sedimentary coverage;
• The Apenninic–Maghrebian Chain Units, originated by the deformation of sedimentary successions originally deposited upon the oceanic crust and carbonate platforms located upon the continental crust. Locally, the Apenninic–Maghrebian Chain Units consist of Flysch-type sequences of the Alpine Tethydes and the Oligocene–Miocene terrigenous cover of the Panormide succession.

In the study area, several polymetallic ores occur in the Ercinian metamorphic rocks. These are located within the Mandanici Unit, consisting of marble layers and graphitic schists with intercalation of ankerite, associated with fluorite, barite, galena and pyrite. Above and laterally to marbles, Cu-Sb-Ag-As, Ni and Bi, with traces of Pb, Zn, W and Au sulphasalt-sulphide mineralizations are present [28], in mineral phases, such as galena, pyrite, chalcopyrite, sphalerite, pyrrhotite, arsenopyrite, the tetrahedrite group of minerals [(Cu,Fe,Ag,Zn)₁₂Sb₄S₁₃] stibnite [29], and scheelite associated with tourmaline [30,31]. Here, a hydrothermal system is recognised along the Tyrrhenian (Terme Vigliatore area) and the Ionian (Ali Terme area) coasts.

Two main kinds of aquifers can be recognized in the area: huge groundwater bodies hosted in the coastal alluvial planes, mainly concentrated in the northern sector of the PC, and minor aquifers at higher elevations, linked to the presence of less permeable (or impermeable) layers interposed between permeable successions. The coastal alluvial aquifers are hosted in porous media, receive a significant lateral recharge from the inner mountain range, as evidenced by isotopic data [32] and are intensively exploited through a huge number of wells distributed along the whole area.

The aquifers of the inner mountain range have minor dimensions and are distributed all along the altitudinal gradient, due to the complex thrust and belt structure articulated in numerous and small tectonic blocks by an impressively developed fault network (Figure 1), which also ensures their mutual hydraulic interconnection. These aquifers are mainly hosted in fractured media and feed springs with flow rates from few tenths to several tens of litres per second.

3. Materials and Methods

Data used in this work refer to major, minor and trace element concentrations and partial pressure of dissolved CO₂ of groundwater from two different sources, characterized by a very different spatial precision. The first dataset is composed of 142 among springs and wells, sampled in the years 2004–2005 within the framework of a work commissioned by the Water and Wastes Department (DAR) of the Sicilian Regional Government to INGV, Sezione di Palermo [33]. Part of these data were discussed in previous papers [11,34], where readers interested to details about sampling criteria and analytical methods can find all the related information. The position of these sampling points was determined using single frequency (L1) GPS receivers, with an average precision of 10 m.

The second dataset is from the paper by Dongarrà et al. [13], and consists of 89 out of the 103 sampling points originally reported; 14 points were not used in the present work due to the impossibility of retrieving their correct positions. The reason is that coordinates of the sampling points were not reported in the original paper, which only showed their locations on a simplified, hand-drafted map at a small nominal scale (1:200,000); the only relevant geographic information tabled as a number was the elevation associated to the identifier of the points.

We adopted the following procedure for transforming the alphanumeric dataset in a georeferenced dataset:

(a)

The original map was georeferenced with the software Qgis, version 3.10, associating to the symbols indicating the position of the main cities and villages their UTM-WGS84 coordinates, as resulting from Google Earth^tm, and using these as Ground Control Points.

(b)

Using, as a base, the previously georeferenced map, we created a shapefile, topologically built as a point vector, containing the location of sampling points associated to their identifiers.

(c)

The previous shapefile was exported as a Keyhole Markup Language (.KML) file, visualized in Google Earth^tm and used for associating elevations to the points, which were inserted as a new field in the attribute table of the shapefile.

(d)

The alphanumeric dataset and the attribute table of the shapefile were ordered for increasing elevations, grouping the data in 3 different subsets following the classification reported in [13] for reducing interpretation ambiguities.

(e)

Elevations as in the original paper were compared to elevations from Google Earth^tm, creating a bijective correspondence between the two datasets, and allowing to attribute the correct identifier to the X,Y,Z positions of the shapefile. The 14 missing points were generated during this step, and are due to differences between the two diverse elevations too large for not generating ambiguities.

The georeferenced dataset, populated with a total of 241 points, 142 from the DAR-INGV project and 89 from [13], was used for plotting concentration maps of selected parameters using the kriging algorithm implemented in Golden Software Surfer, release 20.2.218.

4. Results

Coordinates and field physical-chemical parameters (pH and Eh), concentrations of major and minor and trace elements discussed in this study are reported as online supplementary material in the Table S1, S2 and S3, respectively.

4.1. Physical-Chemical Parameters and Concentration of Major Elements

The groundwater showed pH values between 2.7 and 8.6 and Eh between −29 and 549 mV, consistent with the interaction processes between groundwater and country rocks.

Following the Piper diagram of Figure 2, the relative abundance of the main dissolved ions indicates that the groundwater composition falls along the Ca²⁺–Mg²⁺ side of the diagram, with variable amount of Na⁺+K⁺, and along the HCO₃⁻–SO₄²⁻ side, with a trend toward the Cl⁻ corner.

4.2. Minor and Trace Element Concentrations

Minor and trace element concentrations range from few µg L⁻¹ for elements like U, Mn, Pb, Cu, As, Cr and Sb, up to more than one hundred thousand for NO₃. Their average, minimum and maximum values are synthesized in Figure 3 and compared to the Maximum Admitted Concentrations (MAC) established by the WHO [35] and the Italian National Legislative decree (D. Lgs.) 31/2001 [36].

The average concentration values are always below both the MACs, whereas the maxima of some chemical species, such as NO₃, Ni, U, As and Sb, exceed them.

4.3. Partial Pressures of Dissolved CO₂

Data of the partial pressure of dissolved CO₂ vary in the range from 0.002 to 0.0864 atm. Higher values are recorded in the alluvial plains of the northern Tyrrhenian sector and, secondarily, immediately northward of Montagna Grande, the highest elevation of the Peloritani Mts.

5. Discussion

The complex and highly fragmented geotectonic setting of the Peloritani Mts, accompanied by the widespread presence of punctual and diffused sources of anthropogenic pollutants, requires a rigorous geostatistical approach for discriminating among the different origins for chemicals dissolved in groundwater. The shape of the point cloud representing the major ion composition of groundwater, in the Piper diagram of Figure 2, reflects the variable chemical nature of the lithologic matrix hosting the aquifers.

The most populated geochemical facies is represented by carbonates with calcium as the dominant cation, associated to variable amounts of magnesium and to a secondary enrichment trend in sodium and potassium. Similarly, the carbonate–bicarbonate pair dominates among anions, whose relative proportions are modulated by two distinct enrichment trends, oriented toward the chlorine and sulfate corners, respectively.

These trends, for both cations and anions, reflect the concurrence, at very variable extents, of different water–rock interaction processes, fostered by the huge tectonic fragmentation: dissolution of carbonates and evaporites, ionic exchanges with clays and leaching of metamorphic minerals. Moreover, and mainly in coastal areas, the intrusion of the saline wedge and the leaching of the sea spray deposited on the ground represent an additional source for halides.

Regarding minor and trace element concentrations, we considered those for which MAC are reported for the WHO [1,35] and Italian legislation [36]. The average values were always below the different MACs, indicating that no hazardous contaminations of groundwater were found, neither geogenic nor anthropogenic. For some species, like NO₃⁻, Ni, U, As and Sb, the maximum recorded values fell above the MACs.

Normal probability plots (PP) were used to identify different population of data for the different analyzed species, pointing to their possible different sources. Data are normally distributed plot along a straight line; therefore, a distribution articulated in segments showing different slopes indicates the presence of various populations with different origin. Figure 4 shows the PP for all the considered species, whereas

Table 1. Main identified populations with relative percentages and threshold values for selected trace elements (as µg L⁻¹).

	I %	Threshold Value	II %	Threshold Value	III %	Threshold Value	IV %
As	50.93	1	45.33	11	3.74
Ba
Cr	91.95	0.5	8.05
Cu	67.87	1	25.80	3	6.33
Ni	90	4.5	10
NO3−	80.1	16,000	12.57	36,000	3.67	44,000	3.66
Pb	40.94	0.1	50.33	0.3	8.73
Sb	84.87	1	15.13
Se	85.1	1.5	9.47	2.5	5.43
U	91.87	5	8.13

summarizes the threshold values (inflection points) separating the different populations and their relative abundances.

Each species shows a polynomial distribution (Figure 4): Ni, Cr, U and Sb are divided into two populations; As, Cu, Se and Pb into three populations; and NO₃⁻ into four populations. Barium does not show significant inflection points, indicating that a single population is present and consequently a unique source to groundwater. Generally, the populations with highest concentrations represent a percentage comprised between 3.66% and 8.73%.

We used the concentrations of the inflection points as “objective” thresholds for grouping, into different clusters, the contour lines of the concentration distribution maps presented hereafter, as a draft for checking the “spatial coherency” among points pertaining to a given data population: in other words, the coherency between the spatial distribution of the concentration values forming a population and one of the possible sources for that chemical species.

5.1. Spatial Distributions of Nitrates, Heavy Metals and PCO₂

The concentration contour maps of the chemical species described in Figure 3 are presented in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14; the last figure illustrates the spatial distribution of PCO₂, which is useful for clarifying the relationship among tectonics, contributions of deep, hydrothermal fluids and meteoric circulation. We used the combination of the thresholds listed in Table 1 with the MACs (both by WHO and Italian legislation) indicated in Figure 3 for clustering the contour lines (see colour scales in the top left corner of each figure); the tectonic lineaments (white lines), location of main mineral deposits and sampling points exceeding the MACs are also reported in the maps.

The geogenic nature for the As source is highlighted in the map of Figure 5, showing that the higher concentrations are found in points associated to metamorphic rocks (see Figure 1 for comparison) and, conversely, the lower concentrations are present where the sedimentary deposits outcrop diffusely.

Of the three points showing concentrations higher than the MAC, the westernmost two are located along major faults oriented SE-NW, emphasizing the possible contribution of thermalised waters drained by these discontinuities. The third one, located on the Ionian side of the Peloritani Mts., falls in an area where several ore bodies are present: the high As concentration here detected, about 44 μg L⁻¹, is due to a particularly low pH, about 2.7, a condition under which this metal is highly soluble.

In contrast, and as expected, the anthropogenic origin of nitrates is evidenced in Figure 6, showing that the highest concentrations were recorded in the coastal areas, and, in particular, at the north in the Barcellona–Milazzo Plain, as previously described by Cangemi et al. [11], who attributed their origin to the wide use of fertilisers in agriculture, eventually associated to urban wastewaters.

A geogenic nature and a specific association to mineral bodies is seen in the Sb (Figure 7) and Cr (Figure 8) spatial distributions. The highest concentrations were found in points located on the Ionian side, hosting the highest number of mineralized ore bodies. The other elements, Pb (Figure 9), Ni (Figure 10) and Cu (Figure 11) are related to more articulated sources, both concentrated and dispersed, as well as both anthropogenic and geogenic. The lower Pb concentrations were mainly recorded in the inner part of the mountain chain; however, outside this area, higher values were found where both sedimentary and metamorphic rocks outcrop (Figure 1).

Punctual geogenic sources are associated to the occurrences of metalliferous ores in the Ionian side of the Peloritani Mts, consisting of Fe-Pb-Zn-Cu-Sb-As-Ag mineralization, accompanied by quartz, fluorite, barite and siderite [31]. Anthropogenic sources are linked to the heavy urbanized coastal areas, both on the Tyrrhenian and Ionian sides, and to the industrial and electric power plants of the Barcellona–Milazzo Plain [11].

A similar behavior is shown by Ni distribution (Figure 10), with lower concentrations characterizing the inner sectors of the mountain chain, and maxima with points exceeding the related MAC mainly located in the easternmost part of the studied area and linked to mineral deposits. The same is for Cu (Figure 11) distribution, with higher concentrations in the coastal areas and lower concentrations in the inner zones. The maximum concentrations are found again close to the mineralization of the Ionian side.

Among the anthropogenic sources, road traffic and industrial activities (refineries or electrical plants) emit atmospheric pollutants (Ni, Cr, Pb, Co and Cu), which reach groundwater bodies after the wet deposition of particulate matter [37]. Relative maxima of Se concentrations (Figure 12) are very scattered with the exception of the mineralized area of the Ionian side previously described.

A different interpretation can be invoked for the Ba distribution (Figure 13). The highest concentrations were found on the Tyrrhenian coastal area, close to the city of Messina and in a SE-NW oriented belt north-eastwards of Montagna Grande. This distribution is similar to that of the dissolved CO₂ partial pressure (Figure 14). These anomalies are due to the uprising of deep CO₂, driven by active tectonic structures, as described by Cangemi et al., [34]. The anomalous areas are delimited northwards by the Aeolian–Tindari–Letojanni Fault System (ATLFS, [38]), southwards by its onshore segment and eastwards by the seismogenic source credited for the 1908 Stretto di Messina earthquake.

Figure 14. Spatial distribution of PCO₂, reporting locations of sampling points and tectonic lineaments (white lines). The thick red lines, with associated red dotted areas, represent the main faults (lines) and seismogenic sources (areas) by [39].

Figure 14. Spatial distribution of PCO₂, reporting locations of sampling points and tectonic lineaments (white lines). The thick red lines, with associated red dotted areas, represent the main faults (lines) and seismogenic sources (areas) by [39].

5.2. Concentration Contour Maps as a Basic Tool for Groundwater Resource Management

The ultimate aim of the geochemical maps presented and discussed in this work is for their use as a basic support for decision makers in establishing governance rules for groundwater resource management. A fundamental tool, considered in all the EU directives devoted to groundwater resource management, is the monitoring of the changes over time of their quantitative and qualitative characteristics [33]. This aspect is crucial now, due to the impressive acceleration of the tangible effects of climatic changes, for maintaining the resilience of the water resource supply systems to these perturbations.

A well-performing monitoring procedure of groundwater resources should be based on (i) the correct individuation of the functional relationships (energy and matter exchange processes) between a particular system (or sub-system) and the external universe (here including confining hydro-structures, the hydrological cycle and the deep hydrothermal circulation) and (ii) fixing a “zero time condition” for following any change over time. These two conditions are difficult to achieve with extremely complex, mixed natural–anthropic contexts, such as the Peloritani Mts Area, as described in this work.

The correct definition of the functional relationships among the different hydro-structures was pursued applying geostatistical tools for the discrimination between geogenic and anthropogenic sources of contaminants. In addition, we extended our analysis to parameters not usually considered in groundwater quality studies, such as dissolved pCO₂, which must be considered in this specific study case as a fundamental proxy of the hydrothermal, non-meteoric component of the groundwater circulation system.

The other relevant point is the procedure adopted for fixing the “zero time condition”. It is likely that a new sampling and analytic campaign, carried out in a time span as short as possible, on all the springs and wells relevant for ensuring the right comprehension of the studied system is the best choice for accomplishing this goal. However, at the same time, this is also a costly solution that may not be affordable, especially in low-income countries.

As an alternative, we proposed a different approach, aimed to recover the best heritage of knowledge acquired in studies conducted during the past decades, highly valuable concerning the quality of the analytical information content, but affected by serious geo-localization issues regarding these data.

6. Conclusions

The Peloritani Mts area is a paradigmatic example of how important applying “objective” geostatistical tools is for discerning among different sources of chemicals for groundwater in complex geo-tectonic settings hosting civil and industrial settlements variably distributed in these territories. The recognition of the proper sources for a given dissolved element/molecule is here made difficult by a geo-puzzle of rocks with very different lithological natures, from carbonate alluvia to metamorphic rocks hosting metalliferous deposits, hydraulically connected by a dense network of faults and fractures.

The fault network is also responsible for the upraising of hydrothermal, deep fluids, which mix with the main meteoric recharge of the aquifers, introducing an additional degree of complexity in this geosystem. Moreover, many chemicals delivered to groundwater bodies, such as the most of heavy metals, have a double potential source, geogenic or anthropogenic, adding further uncertainty to their discrimination.

Geostatistical tools, specifically pairing probability plots to distribution maps of concentration of the studied elements, can be very useful in the correct identification of the processes driving groundwater geochemistry. Probability plots allow the identification of different populations of data, avoiding subjective classifications that may lead to a potential misinterpretation of the available information. The spatial representation of these families of data, using concentration contour maps, is resolutive for validating their physical coherency.

A potential source for a given chemical can be considered correctly identified if the spatial distribution of the different families of its concentration values is (i) represented by coherent sets of contour lines, e.g., assemblages of relative minima and maxima encompassing different points and (ii) is congruent with that of the potential sources for that chemical. Conversely, if the families identified by the probability plots generate random assemblages of contour lines, with peaks and troughs centered on single points (random, high-frequency spatial noise), these data families are a mere statistical accident deprived of any physical meaning.

The concentration contour line maps presented and described in this work, and based on the different data families identified by the probability plots, are spatially coherent and congruent with the spatial assemblages of the rock matrices and of the anthropogenic sources for the chemicals delivered to the groundwater bodies, indicating the successful use of our geostatistical approach. The spatial coherency of the concentration contour line maps, not affected by high frequency spatial noise, implies that we used a valid method for attributing the position information to the old data. This last consideration finds a general application to procedures aimed to recover good quality geochemical data acquired in the pre-geomatic era, which cannot be properly used if their spatial properties (positions in the X,Y,Z topographic space) remain unknown.