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.

The circulation of groundwater is influenced by several factors, including permeability changes due to the dynamics of the stress field acting along tectonic discontinuities. Open faults can act as preferential pathways for the escape of crustal and subcrustal gases, such as CO₂, towards the surface, allowing their chemical interaction with meteoric fluids. Here, we present the first map of dissolved CO₂ partial pressure in the groundwater of Sicily, Italy. Based on the chemical analysis of 557 samples of groundwater, which were collected in wells and springs, we calculated the related CO₂ partial pressure (PCO₂) using PHREEQC software. The spatial distribution of the calculated PCO₂ values highlights a general homogeneity at the regional scale, interrupted by positive anomalies linked to the main seismogenic or volcanic areas. Faults work as preferential escape pathways for deep CO₂. The subsequent CO₂–water–rock interaction determines the chemical composition of circulating water. As a consequence, groundwater composition can be successfully used to identify anomalous degassing areas.