Hydrogeological Behaviour and Geochemical Features of Waters in Evaporite-Bearing Low-Permeability Successions: A Case Study in Southern Sicily, Italy

Featured Application: The research suggests an approach for refining the guidelines to be used in studying heterogeneous media and planning optimal monitoring networks and protocols for several anthropogenic purposes (e.g., environmental monitoring of landfills or contaminated sites managing). Abstract: Knowledge about the hydrogeological behaviour of heterogeneous low-permeability media is an important tool when designing anthropogenic works (e.g., landﬁlls) that could potentially have negative impacts on the environment and on people’s health. The knowledge about the biogeochemical processes in these media could prevent “false positives” when studying groundwater quality and possible contamination caused by anthropogenic activities. In this research, we ﬁrstly reﬁned knowledge about the groundwater ﬂow ﬁeld at a representative site where the groundwater ﬂows within an evaporite-bearing low-permeability succession. Hydraulic measurements and tritium analyses demonstrated the coexistence of relatively brief to very prolonged groundwater pathways. The groundwater is recharged by local precipitation, as demonstrated by stable isotopes investigations. However, relatively deep groundwater is clearly linked to very high tritium content rainwater precipitated during the 1950s and 1960s. The deuterium content of some groundwater samples showed unusual values, explained by the interactions between the groundwater and certain gases (H 2 S and CH 4 ), the presences of which are linked to sulfate-reducing bacteria and methanogenic archaea detected within the saturated medium through biomolecular investigations in the shallow organic reach clayey deposits. In a wider, methodological context, the present study demonstrates that interdisciplinary approaches provide better knowledge about the behaviour of heterogeneous low-permeability media and the meaning of each data type.


Introduction
The Mediterranean region was affected by a pervasive "salinity crisis" during the Messinian, when it was progressively restricted and partially isolated from the Atlantic Ocean by a combination Appl. Sci. 2020, 10 of tectonic and glacio-eustatic processes [1][2][3][4] that resulted in the deposition of large volumes of evaporitic sediments [5,6]. The Messinian Sicilian Basins are very important geologic systems for analysing these evaporite successions, in view of their lateral variations and subsequent deformation [7]. Their syn-tectonic evolution is in fact fundamental to reconstructing the timing and geometry of the propagating thrust belt [8,9]. In particular, progressive filling of sub-basins led to the formation of aquifer systems characterised by the coexistence of very low permeability clay successions and evaporitic lenses/horizons [10]. From a hydrochemical point of view, both low-and high-salinity groundwaters may be found associated with the mineralogical features described above. Several hydrogeological and geochemical studies have been carried out in such settings. Some of these studies provided preliminary characterisations through the analysis of single sampling campaigns [11][12][13], whereas others investigated in greater detail the hydrogeological behaviour and the hydrochemical evolution over time based on more prolonged and multidisciplinary approaches [14][15][16]. Worldwide, a variety of conceptual models have been proposed for groundwater circulation in evaporite deposits. In some systems, the hydrogeology and hydrochemistry of the evaporitic aquifers are significantly influenced by deep ascending regional fluids [17][18][19][20]; conversely, at other sites, the hydrogeology and hydrochemistry of the studied system are influenced by the mixing of saline groundwater (flowing through evaporitic rocks) with more diluted waters flowing through nearby porous aquifers [21], with the possible influence of surface waters [22,23] or peat layers [24].
The main objective of this study was to refine knowledge about the hydrogeological behaviour of such systems. Because of the expected complexity, the study was carried out by merging hydraulic head measurements, isotopic analyses, and microbial community investigations, as to acquire a broad spectrum of complementary hydraulic and biogeochemical information.
The hydraulic head measurements were carried out to analyse the groundwater flow field from a three-dimensional perspective and investigate the existence of possible vertical flow components influencing the groundwater pathways and residence times.
The isotopic analyses were conducted to refine knowledge about the groundwater origins (stable isotopes of oxygen and deuterium) and residence times (tritium). Isotopes, together with chemical features, are among the major groundwater tracers traditionally and widely employed in hydrogeological studies (e.g., [25][26][27][28][29]).
The microbial community investigations were performed to analyse the isotopic signatures of the groundwater from a biogeochemical perspective, as to avoid incorrect interpretations of geochemical data. Microbial communities are influenced by the physicochemical features of the environment in which they live and are excellent investigative tools in several hydrogeological scenarios. The efficacy of using microorganisms for specific hydrogeological purposes has been verified in several settings that are partially comparable with that addressed in this study, such as karstified media (e.g., [30,31]), low-permeability media (e.g., [32,33]), and high-salinity groundwaters (e.g., [34]).
The usefulness and the efficacy of a coupled isotopic-microbiological approach has been verified in other complex hydrogeological settings (e.g., [35]), taking into consideration how microorganisms migrate in the subsurface (e.g., [36]).
Taking into consideration the key items highlighted in previous papers, this study was devoted mainly to examining in depth the role of certain hydraulic heterogeneities in influencing the groundwater flow field, as well as the role of microbial communities in influencing the isotopic signature of groundwater. Based on this goal, and according to the successful results of former studies performed at the same scale in similar setting [14][15][16], this first step of the research was carried out within a relatively narrow experimental site (in the order of 1 km 2 ). In fact, minimising the extension of the study site at this stage, allowed to minimise the distance between investigation boreholes and wells, and maximise the opportunity to understand the heterogeneity degree of the studied system from the geological, hydrogeological and biogeochemical points of view.

Study Area
The study area is located in south-western Sicily, between the small towns of Siculiana and Montallegro ( Figure 1). It geologically corresponds (Figure 2a; [37]) to the south to southeast-vergent Sicilian Fold and Thrust Belt (FTB, [38]) bordered by the Gela thrust front and Kabilian-Calabrian thrust front. The FTB is an element of the collisional complex of Sicily, which also includes a late Pliocene-Quaternary foredeep (Gela foredeep), onlapping the frontal sector of the thrust belt in the southern part of the island, its offshore area, and the Pelagian-Iblean foreland with its African crust. In particular, the southern and central parts of the FTB (Figure 2b), are characterised by the outcropping of the most complete evaporitic succession of the "Gessoso-Solfifera" formation of the Messinian age [39].
The evaporitic deposits ( Figure 3) mainly consist of limestone, gypsum, salt, and numerous intercalations of clays, marls, and carbonates. Gypsum strata are composed of various kinds of selenite (branching, banded, and massive) and detrital gypsum layers alternating with marl and carbonate.

Study Area
The study area is located in south-western Sicily, between the small towns of Siculiana and Montallegro ( Figure 1). It geologically corresponds (Figure 2a; [37]) to the south to southeast-vergent Sicilian Fold and Thrust Belt (FTB, [38]) bordered by the Gela thrust front and Kabilian-Calabrian thrust front. The FTB is an element of the collisional complex of Sicily, which also includes a late Pliocene-Quaternary foredeep (Gela foredeep), onlapping the frontal sector of the thrust belt in the southern part of the island, its offshore area, and the Pelagian-Iblean foreland with its African crust. In particular, the southern and central parts of the FTB (Figure 2b), are characterised by the outcropping of the most complete evaporitic succession of the "Gessoso-Solfifera" formation of the Messinian age [39].  The evaporitic deposits ( Figure 3) mainly consist of limestone, gypsum, salt, and numerous intercalations of clays, marls, and carbonates. Gypsum strata are composed of various kinds of selenite (branching, banded, and massive) and detrital gypsum layers alternating with marl and carbonate. The Messinian units overlie diatomitic laminites (Tripoli formation [Fm.]: TRP) and clay deposits of late Tortonian to Early Messinian age, and the succession is overlain by Pliocene marly calcilutites (Trubi Fm.: TRB). The whole sequence [41] is locally overlain by middle-upper Pliocene marly clays (Monte Narbone Fm.: NRB). Sandy clays and arenites (Montallegro Fm.: MNT), clays breccias (BRC), and turbiditic calcarenites (Agrigento Fm.: GRG) are the more recent Pleistocene units.
In greater detail, a first and a second sedimentary cycle [42] are distinguished in the Messinian evaporites, which are separated by an angular unconformity. The first cycle comprises [43]: carbonate deposits (Calcare di Base Fm.: BEC), massive selenite (Cattolica Fm.: CTL), and a salt unit (Clastic Casulfates, Mg and K salts, and halite unit: SLT). The younger second cycle is mainly characterised by In greater detail, a first and a second sedimentary cycle [42] are distinguished in the Messinian evaporites, which are separated by an angular unconformity. The first cycle comprises [43]: carbonate deposits (Calcare di Base Fm.: BEC), massive selenite (Cattolica Fm.: CTL), and a salt unit (Clastic Ca-sulfates, Mg and K salts, and halite unit: SLT). The younger second cycle is mainly characterised by thin gypsum layers (balatino and selenite: GPQ) and marls layers (upper evaporites, Pasquasia Fm.: GPQ) levels interbedded with detrital mud, silt, sandstones, and conglomerates overlain by siliciclastic sediments (mainly silty clays; Arenazzolo Fm.: RNZ), characterised by marked organic matter levels.
The complexity of the geological setting of the area is increased to be a consequence of the Plio-Pleistocene tectonic phases that generated the fold and thrust belts and high-angle faults with subsequent lateral contacts between the Messinian gypsum units and the older clayey deposits. In addition, the wide distribution of soluble rocks in this sector of Sicily is responsible for intense karstic processes, producing a great variety of either epigeous and hypogeous landforms [41].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 24 The complexity of the geological setting of the area is increased to be a consequence of the Plio-Pleistocene tectonic phases that generated the fold and thrust belts and high-angle faults with subsequent lateral contacts between the Messinian gypsum units and the older clayey deposits. In addition, the wide distribution of soluble rocks in this sector of Sicily is responsible for intense karstic processes, producing a great variety of either epigeous and hypogeous landforms [41].
From a morpho-structural standpoint, the landscape is characterised by a system of gentle anticlines and synclines extending parallel to the NW-SE trend of the FTB and large blocks of evaporitic rocks encircled by clayey deposits, with frequent relief inversion produced by differential erosion phenomena [39]. The folded relief is also intersected by perpendicular faults connecting different parts of the evaporites bodies and/or clayey deposits. The area is characterised by a series of slightly sloped fluvial terraces, frequently cut in their upper parts by long gullies and V-shaped valleys. Karst processes have affected gypsum areas, resulting in a great variety of surface landforms, such as poljes, karren, and gypsum bubbles. On the clayey slopes, morphogenetic processes produced landforms such as shallow landslides, rills, and gullies [45].
Following the above-mentioned general geological setting, in the study area, the Mio-Pliocene From a morpho-structural standpoint, the landscape is characterised by a system of gentle anticlines and synclines extending parallel to the NW-SE trend of the FTB and large blocks of evaporitic rocks encircled by clayey deposits, with frequent relief inversion produced by differential erosion phenomena [39]. The folded relief is also intersected by perpendicular faults connecting different parts of the evaporites bodies and/or clayey deposits. The area is characterised by a series of slightly sloped fluvial terraces, frequently cut in their upper parts by long gullies and V-shaped valleys. Karst processes have affected gypsum areas, resulting in a great variety of surface landforms, such as poljes, karren, and gypsum bubbles. On the clayey slopes, morphogenetic processes produced landforms such as shallow landslides, rills, and gullies [45].
Following the above-mentioned general geological setting, in the study area, the Mio-Pliocene deposits consisting of an evaporite series and pre-evaporite strata are involved in a wide range of fold structures on a NNW-SSE axis and are capped by weakly deformed Pleistocene deposits and Quaternary terrains. The shapes of the structures are irregular fashion because of the different stratigraphic thicknesses of the evaporitic rocks, which were also controlled by active thrust and fold structures at the time of deposition [7] and by the significant competence contrasts between various formations described in Figure 3, such as coarse-grained gypsum, muds, silty clays, and silts. In general, the evaporitic deposits of the second cycle (gypsum and marls; GPQ1; GPQ2 and GPQ3) locally overlie the massive selenites of the CTL (first cycle) and the SCD with an erosional unconformity, and are in turn overlain by the silty clays of the RNZ.
To reconstruct the local geology, the stratigraphic logs of 18 geognostic boreholes and 22 piezometers were taken into consideration, and a geological map ( Figure 4) was constructed, which is marked by the outcropping of several components of the whole succession. Based on the available stratigraphic data, two interpretative profiles extending E-W and N-S were also drawn ( Figure 5).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 24 general, the evaporitic deposits of the second cycle (gypsum and marls; GPQ1; GPQ2 and GPQ3) locally overlie the massive selenites of the CTL (first cycle) and the SCD with an erosional unconformity, and are in turn overlain by the silty clays of the RNZ.
To reconstruct the local geology, the stratigraphic logs of 18 geognostic boreholes and 22 piezometers were taken into consideration, and a geological map ( Figure 4) was constructed, which is marked by the outcropping of several components of the whole succession. Based on the available stratigraphic data, two interpretative profiles extending E-W and N-S were also drawn ( Figure 5). The two profiles depict the general expected setting of the area, with a folded regular sequence marked by first cycle (CTL) gypsum layers overlaying the diatomites (TRP) levels or directly in contact on the Tortonian clayey deposits (SCD). The second cycle evaporitic layers (GPQ1, GPQ2, GPQ3, and RNZ) cover in sequence the CTL only along the north-south profile, and elsewhere directly overlie the SCD. At the base, the whole area is characterised by the continuous thick Tortonian clayey layer (SCD). In the southern flank of the anticline, a thick sandy/sandy-clayey The two profiles depict the general expected setting of the area, with a folded regular sequence marked by first cycle (CTL) gypsum layers overlaying the diatomites (TRP) levels or directly in contact on the Tortonian clayey deposits (SCD). The second cycle evaporitic layers (GPQ1, GPQ2, GPQ3, and RNZ) cover in sequence the CTL only along the north-south profile, and elsewhere directly overlie the SCD. At the base, the whole area is characterised by the continuous thick Tortonian clayey layer (SCD). In the southern flank of the anticline, a thick sandy/sandy-clayey (MNT) layer outcrops, covering the RNZ. In the central sector of the study area, along the main drainage surficial line, the evaporitic terrains are covered by a large patch of mainly alluvial quaternary deposits and characterised by high organic matter content.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 24 (MNT) layer outcrops, covering the RNZ. In the central sector of the study area, along the main drainage surficial line, the evaporitic terrains are covered by a large patch of mainly alluvial quaternary deposits and characterised by high organic matter content.

Materials and Methods
An interdisciplinary approach was employed in order to analyse the groundwater flow field and the hydrogeological behavior of the studied system from a three-dimensional perspective. The possible influence of vertical flow on groundwater pathways and residence times was explored by merging hydrogeological and isotopic investigations. Microbial community analyses were conducted to investigate the isotopic signature of the groundwater from a biogeochemical perspective, to ensure that the results of the isotopic analyses were correctly interpreted, with emphasis on the deuterium content.
Graphic images were obtained using QGIS, Adobe Illustrator and Microsoft Excel.

Hydrogeological Investigations
From May 2011 to November 2019, the hydraulic head at the piezometers was measured on a weekly to monthly basis using a water level meter to reconstruct the groundwater flow net and its evolution over time. To specifically investigate the vertical zonation of the hydraulic head and the isotopic/microbiological features within the heterogeneous medium, one multilevel groundwater monitoring system (cluster type) was drilled in August 2019: the C5 cluster ( Figure 6), which included three piezometers screened at different depths (C5a, C5b, and C5c). The C5 cluster was drilled close to a former and abandoned piezometer (30 m deep and fully screened), in which the sampled groundwater was characterised by high tritium content (36.5 tritium units [T.U.]), to also deeply investigate the origin of this apparent anomaly. The other piezometers ranged in depth between 15 and 40 m, depending on the local depth of the groundwater head.

Materials and Methods
An interdisciplinary approach was employed in order to analyse the groundwater flow field and the hydrogeological behavior of the studied system from a three-dimensional perspective. The possible influence of vertical flow on groundwater pathways and residence times was explored by merging hydrogeological and isotopic investigations. Microbial community analyses were conducted to investigate the isotopic signature of the groundwater from a biogeochemical perspective, to ensure that the results of the isotopic analyses were correctly interpreted, with emphasis on the deuterium content.
Graphic images were obtained using QGIS, Adobe Illustrator and Microsoft Excel.

Hydrogeological Investigations
From May 2011 to November 2019, the hydraulic head at the piezometers was measured on a weekly to monthly basis using a water level meter to reconstruct the groundwater flow net and its evolution over time. To specifically investigate the vertical zonation of the hydraulic head and the isotopic/microbiological features within the heterogeneous medium, one multilevel groundwater monitoring system (cluster type) was drilled in August 2019: the C5 cluster ( Figure 6), which included three piezometers screened at different depths (C5a, C5b, and C5c). The C5 cluster was drilled close to a former and abandoned piezometer (30 m deep and fully screened), in which the sampled groundwater was characterised by high tritium content (36.5 tritium units [T.U.]), to also deeply investigate the origin of this apparent anomaly. The other piezometers ranged in depth between 15 and 40 m, depending on the local depth of the groundwater head. Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 24

Isotopic Investigations
Groundwater samples for stable isotope (δ 18/16 O and δ 2/1 H) and tritium ( 3 H) analyses were collected four times (January, June, September and December 2019), together with the hydraulic head measurements. During the same observation period, rainwater samples for δ( 18/16 O), δ( 2/1 H) and tritium analyses were collected monthly (depending on the amount of precipitation) from one rain sampler installed at the test site ( Figure 4). The rainfall was collected using ten-liter polyethylene bottles containing about 300 mL of Vaseline oil to prevent evaporation. Oil contamination was carefully avoided by syringing the water samples out of the bottles. All the samples were transported to the laboratory in a refrigerated box.
Stable isotope analyses were carried out at the Isotope Geochemistry Laboratory of the University of Parma, Italy, using a water equilibrator at 18 °C equipped with a Finnigan Delta XP spectrometer. For oxygen isotope determination, 5 cm 3 of water was equilibrated with pure CO2, whereas for hydrogen isotopes, 5 cm 3 of water was equilibrated with pure H2 (platinum wire was used as a catalyser of gas-liquid water equilibration). The isotope ratio is expressed as: where A is 18 O or 2 H, B is 16 O or 1 H, R is the ratio of the isotopic abundances, i the sample of interest, ‰ = 10 −3 , and RF is the primary international standard of reference (in our case VSMOW-SLAP). The analyses of 3 H were carried out at the Isotope Geochemistry Laboratory of the University of Trieste, Italy. Following other authors (e.g., [46]), the samples followed the procedure of the preventive electrolytic enrichment of tritium was applied to the samples to decrease measurement errors. In this process, the 250 g of the water sample was expected to be reduced to 20 g by electrolysis.
The analyses for the determination of the tritium activity were carried out according to the procedures provided by the International Atomic Energy Agency [47].

Isotopic Investigations
Groundwater samples for stable isotope (δ 18/16 O and δ 2/1 H) and tritium ( 3 H) analyses were collected four times (January, June, September and December 2019), together with the hydraulic head measurements. During the same observation period, rainwater samples for δ( 18/16 O), δ( 2/1 H) and tritium analyses were collected monthly (depending on the amount of precipitation) from one rain sampler installed at the test site ( Figure 4). The rainfall was collected using ten-liter polyethylene bottles containing about 300 mL of Vaseline oil to prevent evaporation. Oil contamination was carefully avoided by syringing the water samples out of the bottles. All the samples were transported to the laboratory in a refrigerated box.
Stable isotope analyses were carried out at the Isotope Geochemistry Laboratory of the University of Parma, Italy, using a water equilibrator at 18 • C equipped with a Finnigan Delta XP spectrometer. For oxygen isotope determination, 5 cm 3 of water was equilibrated with pure CO 2 , whereas for hydrogen isotopes, 5 cm 3 of water was equilibrated with pure H 2 (platinum wire was used as a catalyser of gas-liquid water equilibration). The isotope ratio is expressed as: where A is 18 O or 2 H, B is 16 O or 1 H, R is the ratio of the isotopic abundances, i the sample of interest, % = 10 −3 , and RF is the primary international standard of reference (in our case VSMOW-SLAP). The analyses of 3 H were carried out at the Isotope Geochemistry Laboratory of the University of Trieste, Italy. Following other authors (e.g., [46]), the samples followed the procedure of the preventive electrolytic enrichment of tritium was applied to the samples to decrease measurement errors. In this process, the 250 g of the water sample was expected to be reduced to 20 g by electrolysis.
The analyses for the determination of the tritium activity were carried out according to the procedures provided by the International Atomic Energy Agency [47]. The analytical prediction uncertainty was ±0.1% for δ 18 O, ±1% for δ 2 H, and ±0.5 T.U. for 3 H.

Microbiological Analyses: 16S Ribosomal RNA Gene Next Generation Sequencing (NGS)
The microbiological survey was carried out in December 2019 on samples from the multilevel piezometers C5a-c to refine knowledge about the bacterial community in the groundwater. This survey was limited to cluster C5 because it was carried out to investigate the possible influence of microbial activity on some isotopic features, with emphasis on deuterium content in the C5 groundwater.
Water samples (1 L) were filtered through sterile mixed esters of cellulose filters (S-Pak TM membrane filters, 47 mm diameter, 0.22 µm pore size, Millipore Corporation, Billerica, MA, USA) within 24 h after collection. Bacterial DNA extraction from the sample filters was performed using the commercial kit FastDNA SPIN Kit for soil and the FastPrep ® Instrument. Once the DNA extraction was complete, the quantity and integrity of the DNA were evaluated by electrophoresis in 0.8% agarose gel, containing 1 µg/mL of Gel-Red TM , in the running buffer TAE 1X (40 mM tris base, 20 mM acetic acid, and 1 mM EDTA pH 8), and amplification reactions were performed by polymerase chain reaction (PCR). The 16S rDNA profiles of the bacterial communities in the samples were obtained by NGS technologies at the Genprobio Srl Laboratory (Parma, Italy). Partial 16S rRNA gene sequences were obtained from the extracted DNA by PCR, using the primer pair Probio_Uni and/Probio_Rev, which targets the V3 region of the bacterial 16S rRNA gene sequence [48], and the primer pair ArchV46 for archaeal 16S rRNA [49]. Amplifications were carried out using a Verity Thermocycler (Applied Biosystems Foster City, CA, USA), and PCR products were purified by the magnetic purification using Agencourt AMPure XP DNA Purification Beads (Beckman Coulter Genomics GmbH, Bernried, Germany) to remove primer dimers. Amplicon checks were carried out as previously described [48]. Sequencing was performed using an Illumina MiSeq sequencer with MiSeq Reagent Kit v3 chemicals. The fastq files were processed using a custom script based on the QIIME software suite [50]. Quality control retained sequences with lengths between 140 and 400 bp and mean sequence quality scores >20, whereas sequences with homopolymers >7 bp and mismatched primers were omitted. To calculate downstream diversity measures, operational taxonomic units (OTUs) were defined at 100% sequence homology using DADA2 [51]; OTUs not encompassing at least two sequences of the same sample were removed. All reads were classified to the lowest possible taxonomic rank using QIIME2 [50,52] and a reference dataset from the SILVA database v132 [53]

Hydrogeological Features
The hydraulic head fluctuated over the year from tens of centimeters to several meters, depending on the area. A recession period was typically observed from late spring to early autumn, whereas recharge was typically observed from early autumn to spring (Figure 7). Overall, the head measurements suggested relatively slow recharge, according to the low bulk permeability of the studied system, which does not favor rapid percolation of fresh-infiltration waters within the unsaturated zone and corresponding short-term fluctuations of the groundwater surface. The low bulk permeability of the studied medium is in agreement with the findings of other authors who performed pumping tests in similar Messinian successions in Southern Italy (transmissivity on the order of 10 −5 m 2 /s; [14,16]).
Based on the hydrogeological monitoring and survey results, the general groundwater flow scheme in the area was reconstructed (Figure 8a). The groundwater flows from the highlands towards the topographic depressions, with the hydraulic gradient ranging from~30% to~3%, respectively. The lowest head was consistently measured within the eastern depression, in the area of piezometers S3-S5, and no significant modifications of the groundwater flow net were observed over time, during the observation period. Therefore, all the piezometers used to analyse the isotopic signature from a hydrogeological perspective maintained the same hydraulic relationships throughout the year.
This finding implies that all the results discussed in the isotopic section were not influenced at all by modifications of the groundwater flow field. Based on the hydrogeological monitoring and survey results, the general groundwater flow scheme in the area was reconstructed (Figure 8a). The groundwater flows from the highlands towards the topographic depressions, with the hydraulic gradient ranging from ∼30% to ∼3%, respectively. The lowest head was consistently measured within the eastern depression, in the area of piezometers S3-S5, and no significant modifications of the groundwater flow net were observed over time, during the observation period. Therefore, all the piezometers used to analyse the isotopic signature from a hydrogeological perspective maintained the same hydraulic relationships throughout the year. This finding implies that all the results discussed in the isotopic section were not influenced at all by modifications of the groundwater flow field.
In light of the sedimentological features and common lateral inhomogeneity of the evaporitic layers (Figure 3), simple expected hydrogeological and hydrochemical behaviors are difficult to define for the sequence that outcrops at the key site of the C5 cluster. In fact, the RNZ layer, which overlies the second cycle gypsum evaporitic bodies, is composed of a very chaotic sandy-clayey deposits, including dispersed gypsum blocks; consequently, a moderate primary permeability is expected, despite the high clay content. In addition, the first cycle evaporitic layers and the basal Tortonian clays are expected to be characterised by high and low permeability, respectively, whereas the role of the second cycle gypsum bodies, which are to be classified as permeable on a small scale, could be strongly controlled, overall, by the lack of lateral continuity, resulting in very slow hydraulic circuits. As observed at the C5 cluster, the hydraulic head significantly varied with depth ( Figure 8b). In greater detail, it progressively decreased with increasing depth, showing the highest head in the shallow system and the lowest one into the piezometer C5c, screened within a gypsum horizon. The head zonation was confirmed based on the seasonal groundwater oscillation, with coherent hydraulic In light of the sedimentological features and common lateral inhomogeneity of the evaporitic layers (Figure 3), simple expected hydrogeological and hydrochemical behaviors are difficult to define for the sequence that outcrops at the key site of the C5 cluster. In fact, the RNZ layer, which overlies the second cycle gypsum evaporitic bodies, is composed of a very chaotic sandy-clayey deposits, including dispersed gypsum blocks; consequently, a moderate primary permeability is expected, despite the high clay content. In addition, the first cycle evaporitic layers and the basal Tortonian clays are expected to be characterised by high and low permeability, respectively, whereas the role of the second cycle gypsum bodies, which are to be classified as permeable on a small scale, could be strongly controlled, overall, by the lack of lateral continuity, resulting in very slow hydraulic circuits.
As observed at the C5 cluster, the hydraulic head significantly varied with depth ( Figure 8b). In greater detail, it progressively decreased with increasing depth, showing the highest head in the shallow system and the lowest one into the piezometer C5c, screened within a gypsum horizon. The head zonation was confirmed based on the seasonal groundwater oscillation, with coherent hydraulic heads shifts, consistent with the fact that the analysed medium behaves as a continuum at the observation scale of the present study. This finding is also in agreement with the expectation for saturated heterogeneous media comparable with the studied media. In these media, if a lower geological formation has a higher hydraulic conductivity than the overlying layer, it acts as the major conduit of flow [54]. The highest head in C5a was compatible with the phreatic surface reconstructed within the whole studied area, the morphology of which reflects the local topography (Figure 8a), as expected in such low-permeability systems.

Hydrogen and Oxygen Isotopes
Waters from all the piezometers P, S, and C shown in Figure 8a  Regarding the affinity of the sampled groundwaters, all groups plotted along the local meteoric water lines available in the scientific literature [55][56][57][58][59][60][61][62], being consistent with those of rainwater collected in the study area during the observation period; this finding suggests a meteoric origin for the analysed waters ( Figure 9), with no influence of ascending fluids coming from deep reservoirs, the isotopic signature of which are notably different from those of actual precipitations [35].
The δ 2/1 H and δ 18/16 O values are correlated with the depth of the phreatic surface ( Figure 10). The piezometers showing the deepest phreatic surface (>25 m below ground, m b.g.) belong to Group A; Group B includes phreatic surface at 23.5 to 17 m b.g.; Group C may be divided into two different sub-groups: one with a depth lower than 6 m b.g., and another one with a level ranging from 14.5 to 17 m b.g. Wells C5b and C5c are excluded because their water levels are representative of the hydraulic head measured at different depths below ground and is not referable to the local phreatic surface (coinciding with the water level measured in C5a).
Regarding the affinity of the sampled groundwaters, all groups plotted along the local meteoric water lines available in the scientific literature [55][56][57][58][59][60][61][62], being consistent with those of rainwater collected in the study area during the observation period; this finding suggests a meteoric origin for the analysed waters (Figure 9), with no influence of ascending fluids coming from deep reservoirs, the isotopic signature of which are notably different from those of actual precipitations [35]. ; Group C may be divided into two different sub-groups: one with a depth lower than 6 m b.g., and another one with a level ranging from 14.5 to 17 m b.g. Wells C5b and C5c are excluded because their water levels are representative of the hydraulic head measured at different depths below ground and is not referable to the local phreatic surface (coinciding with the water level measured in C5a). Regarding tritium, the rainwater samples had high 3 H contents (6.2 to 10.8 T.U.), whereas the groundwater samples showed a wide range of tritium contents (Figure 11). Taking into consideration the tritium content detected in local rainwater during the observation period, as well as data already available for Southern Italy (e.g., 4.5 T.U. in [63]; 5.0 T.U. in [16]), three groundwater types have been detected at the study site: (i) groundwater related to rapid and/or short pathways within the aquifer system (4 to 10 T.U.), (ii) groundwater related to longer and/or prolonged pathways, with higher mean-residence times (<4 T.U.), (iii) groundwater related to very prolonged pathways (>10 T.U., corresponding to at least 300 T.U. at the beginning of 1960s), the recharge of which, according to findings in other sites in Italy [64], is linked mainly to rainwater precipitated during the 1950s and 1960s, when tritium was spiked in the atmosphere by nuclear weapons testing [65], reaching up to several thousand T.U. in rainwater (e.g., [66]). Figure 11 reports the relation between 10 3 δ 2/1 H and T.U. The evident and strongly significant correlation (R 2 = 0.998, and null hypothesis Ho(slope=0) < 0.001) for the data of the C5 cluster is an intriguing problem. Different explanations for this interesting correlation are proposed below. Regarding tritium, the rainwater samples had high 3 H contents (6.2 to 10.8 T.U.), whereas the groundwater samples showed a wide range of tritium contents (Figure 11). Taking into consideration the tritium content detected in local rainwater during the observation period, as well as data already available for Southern Italy (e.g., 4.5 T.U. in [63]; 5.0 T.U. in [16]), three groundwater types have been detected at the study site: (i) groundwater related to rapid and/or short pathways within the aquifer system (4 to 10 T.U.), (ii) groundwater related to longer and/or prolonged pathways, with higher mean-residence times (<4 T.U.), (iii) groundwater related to very prolonged pathways (>10 T.U., corresponding to at least 300 T.U. at the beginning of 1960s), the recharge of which, according to findings in other sites in Italy [64], is linked mainly to rainwater precipitated during the 1950s and 1960s, when tritium was spiked in the atmosphere by nuclear weapons testing [65], reaching up to several thousand T.U. in rainwater (e.g., [66]).
1960s, when tritium was spiked in the atmosphere by nuclear weapons testing [65], reaching up to several thousand T.U. in rainwater (e.g., [66]). Figure 11 reports the relation between 10 3 δ 2/1 H and T.U. The evident and strongly significant correlation (R 2 = 0.998, and null hypothesis Ho(slope=0) < 0.001) for the data of the C5 cluster is an intriguing problem. Different explanations for this interesting correlation are proposed below.  Figure 11 reports the relation between 10 3 δ 2/1 H and T.U. The evident and strongly significant correlation (R 2 = 0.998, and null hypothesis H o(slope=0) < 0.001) for the data of the C5 cluster is an intriguing problem. Different explanations for this interesting correlation are proposed below.

Next-Generation Sequencing Results
The microbial community of groundwater sampled in the C5 cluster was analysed to explore the possible influence of microbial activity (H 2 S and CH 4 production) on the unusual deuterium content detected in C5-waters. The 16S rRNA gene sequences of both bacteria and archaea identified in this study have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive under the accession number PRJNA631131. Based on analysis of the sequences obtained from the C5 cluster with NGS technologies, the broad diversity of the microbial community in terms of relative abundance of bacterial genera has been highlighted (Table 1).
Regarding the main aim of the present study, the analyses showed a bacterial community characterised by aerobic and anaerobic/facultative anaerobic genera. Many of these genera are common in bacterial communities found in salty groundwater and marine environments. Halophilic bacteria can be considered moderate halophiles or extreme halophiles based on the concentration of salts present in water and soils (e.g., Alcanivorax [67] and Marinomonas [68] in Table 2); in contrast, halotolerant bacteria, despite growing better without NaCl, can tolerate high quantities of salt [69] (e.g., Arcobacter [70] and Pseudomonas [71] in Table 2). In addition, some of the bacterial genera can exploit sulfate reduction (Desulfatiglans [72], Desulfovibrio [73], and Shewanella [74] in Table 2) as a metabolic pathway. Sulfate reduction is a form of anaerobic respiration, typical of sulfate-reducing bacteria (SRB) or sulfur-reducing bacteria, and as such typically occurs in anoxic environments and leads to the production of hydrogen sulfide (H 2 S); however, some studies have found that this process can also occur in the presence of oxygen [71,75]. Regarding the archaeal component, the most abundant genus found in the analysed groundwater samples (C5a-c) was Candidatus Nitrosopumilus [76] (Table 3). This genus is typically found in marine environments where ammonia can be oxidised to nitrite. Moreover, halophilic genera were identified, such as Halococcus [77] (Table 3), a genus capable of living in the presence of very high concentrations of NaCl (up to 4.5 M NaCl). Methanogenic genera, such as Methanobacterium [78], Methanimicrococcus [79], and Methanolobus [80] were also found ( Table 3). Methanogenesis is a strictly anaerobic metabolic pathway that leads to the formation of methane, and s performed only by methanogenic archaea.
Thus, the biomolecular investigations demonstrated that C5-groundwaters are characterised by bacterial (Desulfatiglans, Desulfovibrio, and Shewanella) and archaeal (Methanobacterium, Methanimicrococcus, and Methanolobus) genera able to produce H 2 S and CH 4 , respectively, therefore suggesting that microbial activity is a factor of utmost importance to be taken into consideration when analysing the isotopic signature of local groundwater.

Hydrogeological Model
Geological, hydrogeological, isotopic and biomolecular investigations allowed characterisation of the hydrogeological functioning of the studied heterogeneous low-permeability system, as well as refinement of knowledge about the residence time of groundwater within such systems.
As expected, taking into consideration previous studies in similar geologic settings [16], the studied system is characterised by vertical and discontinuous heterogeneity [81], caused by various factors such as (i) the coexistence of different geological formations, (ii) their folded regular sequence, and (iii) the fractured (karstified) evaporitic intrastrata/lenses that locally increase the bulk permeability of the low-permeability system. In contrast with other heterogeneous systems, no fault zones influence fluid flow, neither enhancing [82][83][84] nor partially or totally impeding it [85][86][87].
When the hydrogeological system is analysed at the basin scale, the whole system can be depicted as a continuous medium, and the groundwater flow net is a smoothed replica of the local topography. Groundwater flow clearly occurs from the highlands towards the topographic depressions.
The vertical heterogeneity of the system causes the hydraulic head to significantly vary with depth. Close to the C5a-c cluster, the hydraulic head decreases with increasing depth, with the lowest head into the karstified gypsum layer (C5c). This head distribution agrees with that observed by Petrella et al. [16] within the peninsular Southern Italy, where the more transmissive and discontinuous localised media drain the surrounding lower permeability rocks up-gradient (lower hydraulic head), whereas they are characterised by higher head down-gradient and feed the surrounding lower permeability media ( Figure 12). Because the vertical heterogeneity is discontinuous, it is unable to provide a continuous solution within the stratigraphic sequence over an extended area. Therefore, in contrast with other systems characterised by vertical heterogeneity, from the hydraulic point of view [64,[88][89][90], no perched aquifers overlie deeper ones, and a unique saturated zone can be assumed.
Overall, the heterogeneous distribution of lower-and higher-permeability rocks within the studied system, as well as the groundwater flow field described above (Figure 12), cause the mean residence time of groundwater to not necessarily be in agreement with the only groundwater flow net depicted in Figure 8. The coexistence of slower and faster flow velocities and the existence of vertical flow, cause a sort of patchwork distribution of tritium content, with higher levels detected at depths extending into evaporitic interstrata/lenses, linked mainly to rainwater precipitated during the 1950s and 1960s, when tritium was spiked in the atmosphere by nuclear weapons testing.
hydraulic head), whereas they are characterised by higher head down-gradient and feed the surrounding lower permeability media ( Figure 12). Because the vertical heterogeneity is discontinuous, it is unable to provide a continuous solution within the stratigraphic sequence over an extended area. Therefore, in contrast with other systems characterised by vertical heterogeneity, from the hydraulic point of view [64,[88][89][90], no perched aquifers overlie deeper ones, and a unique saturated zone can be assumed.  Conceptual hydrogeological sketch of heterogeneous media characterised by evaporite-bearing successions. The blue lens is a discontinuous evaporitic layer immersed in a clay sequence (grey). Blue dashed lines are the equipotential lines. Blue arrows are the groundwater flow lines. The C5a-c cluster is also shown. The symbols K1 and K2 are the hydraulic conductivities of the evaporitic lens and clay sequence, respectively.
The stable isotopes confirmed the local origin of the groundwater. However, as indicated above, the isotopic features of the samples coming from C5a-c cluster are unusual for groundwaters, and may be due to the phenomena discussed hereafter. Based on data reported by [91] for partitioning between liquid and gaseous water and by Galley et al. [92] for gaseous water and H 2 S gas, the following approximate value for the fractionation factor of H 2 O liquid vs. H 2 S gas at 21 • C is obtained: However, based on the data reported by Horibe and Craig [93], for H 2 O liquid vs. CH 4 gas, the value: Thus, liquid water is strongly enriched in heavy isotopes with respect to H 2 S and CH 4 . Based on the presence of sulfate-reducing bacteria and methanogenic archaea in the analysed C5-waters, we suggest two possible explanations for the high hydrogen isotope values: (1) interaction of water with gaseous H 2 S and/or CH 4 coming from the surrounding environment, and (2) generation of H 2 S and CH 4 by reduction of oxidised sulfur and carbon species involving the H component of the water (Equations (4) and (5)) [94,95]: Here, we consider the case of continuous seepage of H 2 S gas and/or CH 4 gas through groundwater with the exchange of hydrogen isotopes between the gases and water. At the highest δ 2/1 H = −8.3% (C5c), using Equations (2) and (3), we obtain: At least for CH 4 , the 2/1 δ values are, for example, in the range of the values characterising the gases related to oil fields worldwide [96].
Next, we consider a closed system under strongly reducing conditions in which H 2 S and/or CH 4 are produced by reduction at the expense of S of sulfate and/or of C of carbonate species in solution and of H of the water molecules (d n H 2 S = −d n H 2 O ). How much H 2 S and/or CH 4 are generated to produce the hydrogen isotope variation observed in the investigated groundwater (i.e., 2/1 δ H tot from −23.5% to −8.3% )? We consider 2/1 δ in the water at the beginning of reduction (in our case, 2 In this case, for n H2O,o = 100, we obtain n CH4 = 0.5 * 100 * 0.0636 = 3.2 (moles). Thus, a few moles per cent of CH 4 formed is sufficient to greatly increase the hydrogen delta values of the groundwater.
It is noteworthy that whereas CH 4 gas is only slightly soluble, H 2 S gas is strongly soluble in water. The weight ratio total S(II-) as H 2 S to total C(IV-) as CH 4 is about 122 ± 3 in pure water at the temperature of 21 • C for partial pressure of these gases ranging from 10 −4 to 3 atm, according to the PHREEQC speciation program [97]; moreover, salinity does not greatly change this result. Already at about 3-4 mg/L of total H 2 S, the waters typically smell; thus, the absence of a bad smell would indicate that H 2 S is absent or present in only a small amount that supports a prevalent long-term role of methane seepage or the relevant role of methane generation.
As described in a previous paragraph, there is an evident correlation between δ 2/1 H and T.U. Two different explanations can be proposed for this behavior: (1) The chemical component carrying 3 H may also be the cause of increasing δ 2/1 H. This carrier could be methane formed by degradation of organic matter in an environment where tritium may reach high values, such as sediments rich in organic matter that acquired a high tritium concentration at the time of atomic weapon experiments. For instance, Eyrolle et al. [98] reported high tritium concentrations in the organic portion of the sediments from the Loire River in Western France.
These findings "demonstrate that tritium from global atmospheric fallout stored in a sedimentary reservoir for decades as organically bound forms". However, the presence of a modern organic substance in deep sediments is impossible in undisturbed settings. (2) The chemical component carrying 3 H may be different from that of increasing δ 2/1 H. Thus, the two processes are separate but correlated. The enrichment in 2 H is caused by exchange with methane, whereas the enrichment in 3 H is attributed to more or less prolonged groundwater pathways. In this scenario, the deepest investigated groundwater (higher 3 H content in C5c) is characterised by a longer residence time and a major link to rainwater precipitated during the 1950s and 1960s. At the same time, a heterogeneous distribution of organic matter and methanogenic archaea within the saturated medium could lead to (i) different amounts of methane and (ii) a more or less prolonged interaction between the gas and groundwater at different depths. Moreover, considering the decrease of hydraulic head with depth and the downward groundwater flow component, the C5a groundwater type (characterised by relatively low δ 2/1 H and 3 H values) may progressively dilute the deepest C5c groundwater type, thus causing a mixing between two waters with extreme δ 2/1 H and 3 H values, namely C5a and C5c.
With the currently available data, we cannot choose between these hypotheses, but the second one seems the most reasonable. Further research is planned to be carried out to refine knowledge about this issue, including the possibility of finding a similar phenomenon in other portions of the heterogeneous medium.

Concluding Remarks
The in-depth analysis carried out at the study site made it possible to refine knowledge about the hydrogeological behavior of evaporite-bearing low-permeability successions at a representative site. In more details:

•
Hydrogeological investigations and isotopic analyses demonstrated the coexistence of relatively brief to very prolonged groundwater pathways, but everywhere recharged by local precipitation; • The deuterium content of some samples showed unusual values for groundwater; however, the coupled isotopic-biomolecular approach allowed to demonstrate that these unusual values are due to the interactions between the groundwater and certain gases (H 2 S and CH 4 ), the presences of which are linked to sulfate-reducing bacteria and methanogenic archaea; • At this stage of the research, these unusual isotopic content was detected in only one sampling sub-area (C5-cluster); nevertheless, this discrepancy between C5-groundwaters and those sampled in the other available wells is not surprising when taking into consideration the (hydro)geological heterogeneity of the investigated system, as well as the results obtained in a similar setting from the biogeochemical points of view [14][15][16]; in this previous case, a "patchwork" of geochemical and microbiological sub-environments was detected within the whole heterogeneous system, therefore demonstrating the coexistence of different microbial communities (with different microbial activity) also at a metric or decametric scale.
Taking into consideration the interesting results obtained at this stage concerning the heterogeneous distribution of biogeochemical processes, and the importance of moving from speculations to statements concerning the role of H 2 S/CH 4 in influencing the local anomaly in deuterium content, a second step of the research is in progress. During this new phase, the study area will be significantly expanded, and the actual data set will be integrated also with measures of gas concentration and their isotopes, organic carbon percentage in both sediments and groundwater.
In a wider, and methodological context, the present study further demonstrates the importance of interdisciplinary approaches when studying such complex hydrogeological systems. An interdisciplinary approach combining well-established investigation methods (in this case, hydrogeological, isotopic and biomolecular methods), allows refinement of knowledge about the hydrogeological features and behavior of heterogeneous media characterised by evaporite-bearing successions, and can help ensure correct interpretation of the meaning of each data type, avoiding misunderstandings and misinterpretations. Thus, this work can be used to refine the guidelines used to study heterogeneous media through purpose-designed interdisciplinary approaches and plan optimal monitoring networks and protocols for several anthropogenic purposes (e.g., environmental monitoring of landfills or contaminated sites).