We collected and analysed 58 samples of groundwater from wells in the Barcellona-Milazzo Plain, one of the most important coastal aquifers of Sicily (Italy), to determine major, minor, and trace element concentrations. In this area, geogenic and anthropogenic sources of heavy metals and other pollutants co-act, making the individuation of the main pollution sources difficult. Our work was aimed at the application of geostatistical criteria for discriminating between these pollution sources. We used probability plots for separating anomalous values from background concentrations, which were plotted on maps and related to possible sources of pollutants. Our results indicate that hydrothermal fluid circulation and the water–rock interaction of country rocks that host mineralized ore deposits generate a significant flux of heavy metals to groundwater, as well as anthropogenic sources like intense agriculture and industrial activities. In particular, NO₃, F, and Ni exceed the Maximum Admitted Concentrations (MACs) established by the WHO and Italian legislation for drinking-water. The spatial distributions of geogenic and anthropogenic sources were so deeply interlocked that their separation was not easy, also employing geostatistical tools. This complex scenario makes the implementation of human health risk mitigation actions difficult, since the flow of pollutants is in many cases controlled by simple water–rock interaction processes.

Vulcano, the southernmost of the Aeolian island arc (Italy), is characterized by a shallow coastal aquifer resulting from the mixing of seawater, meteoric recharge and volcanogenic fluids. The aquifer has been intensively studied during the last decades, but a comprehensive hydrogeological model has never been developed due to the lack of direct information about the litho-stratigraphic columns of the wells and the depth of water bearing levels. We present and discuss here the time and spatial analysis of water table elevation, temperature and electric conductivity data, acquired during the last 20 years in 33 wells located at Vulcano Island, with the aim of developing a groundwater circulation scheme able to fit the field observations. We retrieved a circulation scheme characterized by an intricate geometry of flow paths driven by horizontal and vertical permeability variations, accounting for the strong variability of geochemical data evidenced in this area by the related scientific literature. Extending these results to a general context, particular care must be taken in approaching the study of aquifers in volcanic islands, because a strong, small spatial scale variability of the hydrogeochemical parameters is expected, and a reliable knowledge of the local conditions is required for developing successful groundwater circulation schemes.

The characterization of a groundwater body involves the construction of a conceptual model that constitutes the base knowledge for monitoring programs, hydrogeological risk assessment, and correct management of water resources. In particular, a detailed geological and geophysical approach was applied to define the alluvial Caronia Groundwater Body (CGWB) and to reconstruct a hydrogeological flow model. The analysis of the CGWB, located in north-eastern Sicily, was initially approached through a reanalysis of previous stratigraphic (boreholes) and geophysical (vertical electrical soundings and seismic refraction profiles) data, subsequently integrated by new seismic acquisitions, such as Multichannel Analysis of Surface Waves (MASW) and horizontal-to-vertical seismic ratio (HVSR). The analysis and reinterpretation of geoelectrical data allowed the construction of a preliminary 3D resistivity model. This initial modeling was subsequently integrated by a geophysical data campaign in order to define the depth of the bottom of the shallow CGWB and the thickness of alluvial deposits. Finally, a preliminary mathematical model flow was generated in order to reconstruct the dynamics of underground water. The results show that integration of multidisciplinary data represent an indispensable tool for the characterization of complex physical systems.

Oxygen isotopic composition is useful for individuating recharge areas of groundwater bodies by the comparison with those of local rainfalls. While on a global scale general relationships, such as the isotopic vertical gradient or continentality effects, efficiently describe spatial variations of the isotopic signature, hydrogeological applications need spatial models that are more focused on the effects of local topographic structures and/or subsoil geology. This work presents a case study in northeastern Sicily (Italy) characterized by complex geological and orographic structures, in which isotopic composition of rainfalls is governed by orographic effects and the varying initial composition of humid air masses. We used a black box approach, comparing the average isotopic composition of rain collected from a network of eight samplers with their spatial descriptors (elevation, latitude and longitude). We obtained the best correlation with the simultaneous use of all these variables, applying their multiple linear correlation equation to transform the 1 × 1 km digital elevation model (DEM) of the study area into a digital isotopic model (DIM). The reliability of the DIM was confirmed by its good agreement with the oxygen isotopic composition contour map of the local groundwater.

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.

The integration of various geophysical methodologies is considered a fundamental tool for accurately reconstructing the extent and shape of a groundwater body and for estimating the physical parameters that characterize it. This is often essential for the management of water resources in areas affected by geological and environmental hazards. This work aims to reconstruct the pattern and extent of two groundwater bodies, located in the coastal sectors of the North-Eastern Sicily, through the integrated analysis and interpretation of several geoelectrical, seismic and geological data. These are the Sant’Agata-Capo D’Orlando (SCGWB) and the Barcelona-Milazzo (BMGWB) Groundwater Bodies, located at the two ends of the northern sector of the Peloritani geological complex. These two studied coastal plains represent densely populated and industrialized areas, in which the quantity and quality of the groundwater bodies are under constant threat. At first, the resistivity models of the two groundwater bodies were realized through the inversion of a dataset of Vertical Electrical Soundings (VES), constrained by stratigraphic well logs data and other geophysical data. The 3D resistivity models obtained by spatially interpolating 1D inverse VES models have allowed for an initial recognition of the distribution of groundwater, as well as a rough geological framework of the subsoil. Subsequently, these models were implemented by integrating results from active and passive seismic data to determine the seismic P and S wave velocities of the main lithotypes. Simultaneous acquisition and interpretation of seismic and electrical tomographies along identical profiles allowed to determine the specific values of seismic velocity, electrical resistivity and chargeability of the alluvial sediments, and to use these values to constrain the HVSR inversion. All this allowed us to recognize the areal extension and thickness of the various lithotypes in the two investigated areas and, finally, to define the depth and the morphology of the base of the groundwater bodies and the thickness of the filling deposits.